model.py

#coding-Utf-8-*-
from gmshpy import *
from dG3Dpy import*

#script to launch PBC problem with a python script

#DEFINE MICRO PROBLEM

# micro-material law
lawnum1 = 11 # unique number of law


rho   = 7850e-9
young = 28.9e3
nu    = 0.3 
sy0   = 99.
h     = young/20.

harden = LinearExponentialJ2IsotropicHardening(1, sy0, h, 0., 10.)
cl     = IsotropicCLengthLaw(1, 1e-2)
damlaw = SimpleSaturateDamageLaw(1, 50.,1.,0.,1.)
law1   = NonLocalDamageJ2HyperDG3DMaterialLaw(lawnum1,rho,young,nu,harden,cl,damlaw)


micromeshfile="micro.msh" # name of mesh file
#micromeshfile="square.msh" # name of mesh file

# creation of part Domain
nfield = 11 # number of the field (physical number of entity)
myfield1 = dG3DDomain(1000,nfield,0,lawnum1,0,2,1)
#myfield1.setNonLocalStabilityParameters(beta1,True)
myfield1.setNonLocalEqRatio(1.e3)

microBC = nonLinearPeriodicBC(1000,2)
microBC.setOrder(1)
microBC.setBCPhysical(1,4,3,2)
# periodiodic BC
method = 0 # Periodic mesh = 0 Langrange interpolation = 1 Cubic spline interpolation =2
degree = 5 
addvertex = False
microBC.setPeriodicBCOptions(method, degree,addvertex)

# DEFINE MACROPROBLEM
matnum1 = 1;
macromat1 = dG3DMultiscaleMaterialLaw(matnum1, 1000)
macromat1.loadModel(micromeshfile);
macromat1.addDomain(myfield1)
macromat1.addMaterialLaw(law1);
macromat1.addMicroBC(microBC)

macromat1.setNumStep(3)
macromat1.setTolerance(1e-6,1e-10)
macromat1.setSystemType(1)
macromat1.Solver(2)
#for i in range(0,1):
macromat1.setViewAllMicroProblems(True,0)

macromat1.setBlockDamageAfterFailureOnset(True)
macromat1.setBlockDamageAfterFailureOnsetTolerance(0.001)

macromat1.Scheme(1)
macromat1.setSameStateCriterion(1e-16)
macromat1.stressAveragingFlag(True) # set stress averaging ON- 0 , OFF-1
macromat1.setStressAveragingMethod(0)
macromat1.tangentAveragingFlag(True) # set tangent averaging ON -0, OFF -1
macromat1.setTangentAveragingMethod(2,1e-6);

macromat1.internalPointBuildView("Green-Lagrange_xx",IPField.STRAIN_XX);
macromat1.internalPointBuildView("Green-Lagrange_yy",IPField.STRAIN_YY);
macromat1.internalPointBuildView("Green-Lagrange_xy",IPField.STRAIN_XY);
macromat1.internalPointBuildView("Green-Lagrange_zz",IPField.STRAIN_ZZ);
macromat1.internalPointBuildView("Green-Lagrange_yz",IPField.STRAIN_YZ);
macromat1.internalPointBuildView("Green-Lagrange_xz",IPField.STRAIN_XZ);
macromat1.internalPointBuildView("sig_xx",IPField.SIG_XX);
macromat1.internalPointBuildView("sig_yy",IPField.SIG_YY);
macromat1.internalPointBuildView("sig_xy",IPField.SIG_XY);
macromat1.internalPointBuildView("sig_xz",IPField.SIG_XZ);
macromat1.internalPointBuildView("sig_zz",IPField.SIG_ZZ);
macromat1.internalPointBuildView("sig_yz",IPField.SIG_YZ);
macromat1.internalPointBuildView("sig_VM",IPField.SVM);

macromat1.internalPointBuildView("F_xx",IPField.F_XX);
macromat1.internalPointBuildView("F_yy",IPField.F_YY);
macromat1.internalPointBuildView("F_xy",IPField.F_XY);
macromat1.internalPointBuildView("F_yx",IPField.F_YX);

macromat1.internalPointBuildView("P_xx",IPField.P_XX);
macromat1.internalPointBuildView("P_yy",IPField.P_YY);
macromat1.internalPointBuildView("P_xy",IPField.P_XY);
macromat1.internalPointBuildView("P_yx",IPField.P_YX);

macromat1.internalPointBuildView("Equivalent plastic strain",IPField.PLASTICSTRAIN);
macromat1.internalPointBuildView("Damage",IPField.DAMAGE);

macromat1.dirichletBC("Face",11,2,0.)


lcohNum = 13
lawCoh = TwoFieldMultiscaleCohesive3DLaw(lcohNum)
lawCoh.setCharacteristicLength(1.)
lawCoh.setLostSolutionUniquenssTolerance(0.05)

macromeshfile="model.msh" # name of mesh file

# solver
sol = 2  # Gmm=0 (default) Taucs=1 PETsc=2
soltype =1 # StaticLinear=0 (default) StaticNonLinear=1
nstep = 100  # number of step (used only if soltype=1)
ftime =1.   # Final time (used only if soltype=1)
tol=1.e-5  # relative tolerance for NR scheme (used only if soltype=1)
nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)


# creation of  macro part Domain

dim =2
beta1 = 50;
fullDG = False;

averageStrainBased = True

# non DG domain
nfield1 = 11
macrodomain1 = dG3DDomain(10,nfield1,0,matnum1,fullDG,dim)
macrodomain1.stabilityParameters(beta1)
macrodomain1.setDistributedOtherRanks(True)
macrodomain1.distributeOnRootRank(False)

nfield2 = 12
macrodomain2 = dG3DDomain(11,nfield2,0,matnum1,fullDG,dim)
macrodomain2.stabilityParameters(beta1)
macrodomain2.setDistributedOtherRanks(True)
macrodomain2.distributeOnRootRank(False)

# interface domain
beta1 = 1e-2;
interdomain1 =  interDomainBetween3D(12,macrodomain1,macrodomain2,lcohNum,matnum1)
interdomain1.stabilityParameters(beta1)
interdomain1.setMultipleFieldFormulation(True,1)
interdomain1.averageStrainBased(averageStrainBased)
interdomain1.setDistributedOtherRanks(True)
interdomain1.distributeOnRootRank(False)
interdomain1.fixInterfaceDofComponent(2)
interdomain1.setBulkDamageBlockedMethod(1)

# creation of Solver
mysolver = nonLinearMechSolver(1000)
mysolver.loadModel(macromeshfile)
mysolver.addDomain(macrodomain1)
mysolver.addDomain(macrodomain2)
mysolver.addDomain(interdomain1)
mysolver.addMaterialLaw(macromat1)

mysolver.addMaterialLaw(lawCoh)
mysolver.setMultiscaleFlag(True)
mysolver.Scheme(soltype)
mysolver.Solver(sol)
mysolver.snlData(nstep,ftime,tol)
mysolver.options("-ksp_type preonly -pc_type lu -pc_factor_shift_type NONZERO -pc_factor_mat_solver_package petsc")


# boundary condition
mysolver.displacementBC("Face",11,2,0.0)
mysolver.displacementBC("Face",12,2,0.0)

mysolver.periodicBC("Edge",1,3,1,4,0)
mysolver.periodicBC("Edge",1,3,1,4,1)

mysolver.displacementBC("Edge",2,0,0.02)
mysolver.displacementBC("Edge",2,1,0.02)
mysolver.displacementBC("Edge",4,0,0.0)
mysolver.displacementBC("Edge",4,1,0.0)


# archivage
mysolver.internalPointBuildView("Green-Lagrange_xx",IPField.STRAIN_XX, 1, 1);
mysolver.internalPointBuildView("Green-Lagrange_yy",IPField.STRAIN_YY, 1, 1);
mysolver.internalPointBuildView("Green-Lagrange_zz",IPField.STRAIN_ZZ, 1, 1);
mysolver.internalPointBuildView("Green-Lagrange_xy",IPField.STRAIN_XY, 1, 1);
mysolver.internalPointBuildView("Green-Lagrange_yz",IPField.STRAIN_YZ, 1, 1);
mysolver.internalPointBuildView("Green-Lagrange_xz",IPField.STRAIN_XZ, 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("sig_VM",IPField.SVM, 1, 1);
mysolver.internalPointBuildView("Green-Lagrange equivalent strain",IPField.GL_EQUIVALENT_STRAIN, 1, 1);
mysolver.internalPointBuildView("Equivalent plastic strain",IPField.PLASTICSTRAIN, 1, 1);

mysolver.archivingForceOnPhysicalGroup("Edge",4,0)


# solve
mysolver.solve()

check = TestCheck()
check.equal(-2.816619e+01,mysolver.getArchivedForceOnPhysicalGroup("Edge", 4, 0),1.e-3)