Skip to content
Snippets Groups Projects
Select Git revision
  • f552f5892f68f7b9dc1f984e0e96c53cf29b9e30
  • master default protected
  • bl
  • pluginMeshQuality
  • fixBugsAmaury
  • hierarchical-basis
  • alphashapes
  • relaying
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • 3115-issue-fix
  • 3023-Fillet2D-Update
  • convert_fdivs
  • tmp_jcjc24
  • gmsh_4_14_0
  • gmsh_4_13_1
  • gmsh_4_13_0
  • gmsh_4_12_2
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
41 results

polynomialBasis.h

Blame
  • SMP.py 6.20 KiB
    #coding-Utf-8-*-
    from gmshpy import *
    
    #from dG3DpyDebug import*
    from dG3Dpy import *
    
    
    #script to launch beam problem with a python script
    
    # material law
    lawnum= 1 # unique number of law
    rho   = 1020.
    #young =Ex=Ey=Ez =2.1e11
    #nu    =Vxy=Vxz=Vyz= 0.3
    #MUxy=MUxz=MUyz=Ex/(2.*(1.+Vxy)) 
    cv=1.
    #Kx=Ky=Kz=51.9
    Kx=Ky=Kz=0.5
    #alphax=alphay=alphaz=12.e-6
    
    mu_groundState3=0.75e6
    Im3=5.
    pi=3.14159
    
    #tinitial = 273+22
    tinitial = 273+58#58  50
    Tr=310
    Nmu_groundState2= 0.045
    mu_groundState2=1.38e6
    Im2=6.3
    Sgl2=58.e6 
    Sr2=3.e2
    Delta=2.6
    m2=0.19
    epsilonr=5.2e-4
    n=2.1
    epsilonp02=5.2e-4
    alphar1=25.e-5
    alphagl1=13.e-5
    Ggl1=156.e6
    Gr1=13.4e6
    Mgl1=7.4e6
    Mr1=0.168e6
    Mugl1=0.35
    Mur1=0.49
    epsilon01=1.73e13
    Qgl1=1.4e-19
    Qr1=0.2e-21
    epsygl1=0.14 
    #epsyr1=0.
    d1=0.015
    Kb=1.3806488e-23
    m1=0.17
    V1=2.16e-27
    alphap=0.058
    Sa0=0.
    ha1=230.
    b1=5850e6
    g1=5.8
    phai01=0.
    Z1=0.083
    r1=1.3
    s1=0.005
    Sb01=0.
    Hgl1=1.56e6
    Lgl1=0.44e6
    Hr1=0.76e6
    Lr1=0.006e6
    l1=0.5
    be1=0.5
    #kb=1.3806488e-23
    v=0.6
    c0=1710.*rho 
    c1=4.1*rho
    
    # geometry
    geofile="cube.geo" # name of mesh file
    meshfile="cube.msh" # name of mesh file
    
    # solver
    sol = 2  # Gmm=0 (default) Taucs=1 PETsc=2
    soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
    nstep =3000# number of step (used only if soltype=1)#640
    #ftime1 =630   # Final time (used only if soltype=1)22 = t= 630
    #ftime  =640   # Final time (used only if soltype=1)22  t = 680
    ftime =360#58 330.
    #ftime =308.   # Final time (used only if soltype=1)
    tol=1.e-4   # relative tolerance for NR scheme (used only if soltype=1)
    nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)
    fullDg = True #O = CG, 1 = DG
    space1 = 0 # function space (Lagrange=0)
    beta1  = 40
    
    alpha = beta=gamma=0.
    #Tangent2
    
    #  compute solution and BC (given directly to the solver
    # creation of law
    
    #law1   = SMPDG3DMaterialLaw( lawnum,rho,Ex, Ey, Ez,Vxy,Vxz,Vyz,MUxy,MUxz,MUyz,alpha,beta,gamma,cv,t0,Kx,Ky,Kz,alphax,alphay,alphaz,mu_groundState3,Im3,pi,
    #Tr,Nmu_groundState2, mu_groundState2,  Im2,Sgl2, Sr2, Delta, m2,epsilonr,n,epsilonp02)
    law1   = SMPDG3DMaterialLaw( lawnum,rho,alpha,beta,gamma,tinitial,Kx,Ky,Kz,mu_groundState3,Im3,pi,
    Tr,Nmu_groundState2, mu_groundState2,  Im2,Sgl2, Sr2, Delta, m2,epsilonr,n,epsilonp02, alphar1, alphagl1, Ggl1, Gr1, Mgl1, Mr1, Mugl1, Mur1, epsilon01, Qgl1,
    		           Qr1, epsygl1 , d1,  m1, V1,  alphap, Sa0, ha1, b1, g1, phai01, Z1,r1, s1, Sb01, Hgl1, Lgl1, Hr1, Lr1, l1, Kb, be1,c0,v,c1)
    
    
    # creation of ElasticField
    nfield = 10 # number of the field (physical number of surface)
    #myfield1 = dG3DDomain(1000,nfield,space1,lawnum,fullDg)
    
    myfield1 = ThermoMechanicsDG3DDomain(1000,nfield,space1,lawnum,fullDg,1.e10)
    
    myfield1.matrixByPerturbation(0,0,0,1e-8)
    #myfield1.matrixByPerturbation(1,1,1,1e-8)
    myfield1.stabilityParameters(beta1)
    myfield1.ThermoMechanicsStabilityParameters(beta1,fullDg)
    # creation of Solver
    mysolver = nonLinearMechSolver(1000)
    mysolver.createModel(geofile,meshfile,3,1)
    mysolver.addDomain(myfield1)
    mysolver.addMaterialLaw(law1)
    mysolver.Scheme(soltype)
    mysolver.Solver(sol)
    mysolver.snlData(nstep,ftime,tol)
    mysolver.snlManageTimeStep(25, 3, 2, 10)
    mysolver.stepBetweenArchiving(nstepArch)
    # BC
    
    #mechanical BC
    
    #cyclicFunctionDisp=cycleFunctionTime(0., 0, ftime1, -0.000628, ftime , -0.00056);#22 working
    #cyclicFunctionDisp=cycleFunctionTime(0., 0, ftime/2., -0.000314, ftime , 0.);
    #cyclicFunctionDisp=cycleFunctionTime(0., 0, ftime1,  -0.000315, ftime,  -0.00028);#22 
    #cyclicFunctionDisp=cycleFunctionTime(0., 0, ftime1,  -0.0004, ftime,  -0.00044);#22
    #cyclicFunctionDisp=cycleFunctionTime(0., 0, ftime/2., -0.000613, ftime , 0.);#58   -0.000613
    cyclicFunctionDisp=cycleFunctionTime(0., 0, 2.*ftime/3., -0.000613, ftime , -0.000613/2.);#58   -0.000613
    
    mysolver.displacementBC("Face",1234,2,0.)
    mysolver.displacementBC("Face",2376,0,0.)
    mysolver.displacementBC("Face",1265,1,0.)
    mysolver.displacementBC("Face",5678,2,cyclicFunctionDisp)
    #mysolver.displacementBC("Volume",nfield,0,0)#when the michenisms terms dosen t have pressure
    #mysolver.displacementBC("Volume",nfield,1,0)#when the mechnisms terms dosen t have pressure
    
    
    #thermal BC
    cyclicFunctionTemp=cycleFunctionTime(0., tinitial, ftime,  tinitial);
    mysolver.displacementBC("Face",1234,3,cyclicFunctionTemp)
    mysolver.displacementBC("Face",5678,3,cyclicFunctionTemp)
    mysolver.displacementBC("Face",2376,3,cyclicFunctionTemp)
    mysolver.displacementBC("Face",1265,3,cyclicFunctionTemp)
    #mysolver.initialBC("Volume","Position",nfield,3,tinitial)
    mysolver.displacementBC("Volume",nfield,3,cyclicFunctionTemp)
    
    #electrical BC
    cyclicFunctionvolt=cycleFunctionTime(0., 0.0,  ftime,0.0);
    mysolver.initialBC("Volume","Position",nfield,4,cyclicFunctionvolt)
    mysolver.displacementBC("Volume",nfield,4,cyclicFunctionvolt)
    
    mysolver.internalPointBuildView("svm",IPField.SVM, 1, 1)
    mysolver.internalPointBuildView("sig_xx",IPField.SIG_XX, 1, 1)
    mysolver.internalPointBuildView("sig_yy",IPField.SIG_YY, 1, 1)
    mysolver.internalPointBuildView("sig_zz",IPField.SIG_ZZ, 1, 1)
    mysolver.internalPointBuildView("sig_xy",IPField.SIG_XY, 1, 1)
    mysolver.internalPointBuildView("sig_yz",IPField.SIG_YZ, 1, 1)
    mysolver.internalPointBuildView("sig_xz",IPField.SIG_XZ, 1, 1)
    mysolver.internalPointBuildView("temperature",IPField.TEMPERATURE, 1, 1)
    mysolver.internalPointBuildView("qx",IPField.THERMALFLUX_X, 1, 1)
    mysolver.internalPointBuildView("qy",IPField.THERMALFLUX_Y, 1, 1)
    mysolver.internalPointBuildView("qz",IPField.THERMALFLUX_Z, 1, 1)
    
    mysolver.internalPointBuildView("voltage",IPField.VOLTAGE, 1, 1)
    mysolver.internalPointBuildView("qx",IPField.ELECTRICALFLUX_X, 1, 1)
    mysolver.internalPointBuildView("qy",IPField.ELECTRICALFLUX_Y, 1, 1)
    mysolver.internalPointBuildView("qz",IPField.ELECTRICALFLUX_Z, 1, 1)
    
    mysolver.archivingForceOnPhysicalGroup("Face", 1234, 2)
    mysolver.archivingForceOnPhysicalGroup("Face", 5678, 2)
    mysolver.archivingForceOnPhysicalGroup("Face", 1234, 3)
    mysolver.archivingForceOnPhysicalGroup("Face", 5678, 3)
    mysolver.archivingIPOnPhysicalGroup("Volume",nfield, IPField.SIG_ZZ,IPField.MEAN_VALUE);
    mysolver.archivingIPOnPhysicalGroup("Volume",nfield, IPField.STRAIN_ZZ,IPField.MEAN_VALUE);
    mysolver.archivingIPOnPhysicalGroup("Volume",nfield, IPField.TEMPERATURE,IPField.MEAN_VALUE);
    
    
    mysolver.solve()
    
    check = TestCheck()
    check.equal(6.768912e-01,mysolver.getArchivedForceOnPhysicalGroup("Face", 1234, 2),1.e-5)