Select Git revision
composite.py

Ludovic Noels authored
composite.py 3.57 KiB
#coding-Utf-8-*-
from gmshpy import *
from dG3Dpy import*
#script to launch beam problem with a python script
# material law
lawnumEP = 2
rhoEP = 1000.
youngEP = 3.2e9
nuEP = 0.3
sy0EP = 25.e6
hEP = 7.10e9
lawnumF = 3
rhoF = 1000.
youngF = 320e9
nuF = 0.3
sy0F = 25.e16
hF = 7.10e9
lawcnumEP = 4
GcEP = 78.
sigmacEP = 83.e6
deltacEP = 2*GcEP/sigmacEP
KcEP = sigmacEP/(0.005*deltacEP)
lawcnumInt = 5
GcInt = 60.
sigmaInt = 30.e6
deltacInt= 2*GcInt/sigmaInt
KcInt = sigmaInt/(0.005*deltacInt)
beta =0.87 # ratio KII/KI
mu = 0. # friction coefficient ??
fsmin = 1.
fsmax = 1.
# geometry
meshfile="composite.msh" # name of mesh file
# solver
sol = 2 # Gmm=0 (default) Taucs=1 PETsc=2
soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
nstep = 500 # number of step (used only if soltype=1)
ftime =1 # Final time (used only if soltype=1)
tol=1.e-4 # relative tolerance for NR scheme (used only if soltype=1)
nstepArch=10 # Number of step between 2 archiving (used only if soltype=1)
fullDg =1 #O = CG, 1 = DG
space1 = 0 # function space (Lagrange=0)
beta1 = 100
# compute solution and BC (given directly to the solver
# creation of law
law1EP = J2LinearDG3DMaterialLaw(lawnumEP,rhoEP,youngEP,nuEP,sy0EP,hEP)
law2EP = LinearCohesive3DLaw(lawcnumEP,GcEP,sigmacEP,beta,mu,fsmin,fsmax,KcEP)
law3EP = FractureByCohesive3DLaw(6,lawnumEP,lawcnumEP)
lawF = J2LinearDG3DMaterialLaw(lawnumF,rhoF,youngF,nuF,sy0F,hF)
law2Int = LinearCohesive3DLaw(lawcnumInt,GcInt,sigmaInt,beta,mu,fsmin,fsmax,KcInt)
# creation of ElasticField
myfieldEP = dG3DDomain(10,78,space1,6,fullDg)
#myfieldEP.matrixByPerturbation(1,1,1,1e-8)
myfieldEP.stabilityParameters(beta1)
myfieldF = dG3DDomain(11,77,space1,lawnumF,fullDg)
#myfieldF.matrixByPerturbation(1,1,1,1e-8)
myfieldF.stabilityParameters(beta1)
#in // we need to create the ininterdomain after adding the domains
myinterfield = interDomainBetween3D(12,myfieldEP,myfieldF,lawcnumInt)
myinterfield.stabilityParameters(beta1)
#myinterfield.matrixByPerturbation(1,1e-10)
# creation of Solver
mysolver = nonLinearMechSolver(1000)
mysolver.loadModel(meshfile)
mysolver.addDomain(myfieldEP)
mysolver.addDomain(myfieldF)
mysolver.addMaterialLaw(law1EP)
mysolver.addMaterialLaw(law2EP)
mysolver.addMaterialLaw(law3EP)
mysolver.addMaterialLaw(lawF)
mysolver.addDomain(myinterfield)
mysolver.addMaterialLaw(law2Int)
mysolver.Scheme(soltype)
mysolver.Solver(sol)
mysolver.snlData(nstep,ftime,tol,1e-12)
mysolver.stepBetweenArchiving(nstepArch)
#mysolver.explicitSpectralRadius(ftime,0.5,0.)
#mysolver.dynamicRelaxation(0.1, ftime, 1.e-3,2)
#mysolver.explicitTimeStepEvaluation(nstep)
mysolver.options("-ksp_type preonly -pc_type lu -pc_factor_mat_solver_package petsc")
# BC
mysolver.displacementBC("Volume",77,2,0.)
mysolver.displacementBC("Volume",78,2,0.)
mysolver.displacementBC("Face",82,0,0.)
mysolver.displacementBC("Face",81,1,0.)
mysolver.displacementBC("Face",80,0,5e-6)
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("epl",IPField.PLASTICSTRAIN, 1, 1)
mysolver.archivingForceOnPhysicalGroup("Face", 82, 0, 1)
mysolver.archivingNodeDisplacement(86,0,1)
mysolver.solve()
check = TestCheck()
check.equal(1.056398e-08,mysolver.getEnergy(3),5.e-3)