composite.py

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