SMP.py

#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)