Select Git revision
import_perf.cpp
Forked from
gmsh / gmsh
Source project has a limited visibility.
model.py 7.01 KiB
#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 = nonLocalDamageDG3DDomain(1000,nfield,0,lawnum1,0,1e3,2,1)
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.setPeriodicity(1.,0,0,"x")
macromat1.setPeriodicity(0,1.,0,"y")
macromat1.setPeriodicity(0,0,1.,"z")
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,0.7)
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,0)
interdomain1.averageStrainBased(averageStrainBased)
interdomain1.setDistributedOtherRanks(True)
interdomain1.distributeOnRootRank(False)
interdomain1.fixInterfaceDofComponent(2)
# 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.displacementBC("Face",11,1,0.0)
#mysolver.displacementBC("Face",12,1,0.0)
mysolver.periodicBC("Edge",1,3,1,4,0)
mysolver.periodicBC("Edge",1,3,1,4,1)
#mysolver.displacementBC("Edge",1,1,0.0)
#mysolver.displacementBC("Edge",3,1,0.01)
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.552598e+01,mysolver.getArchivedForceOnPhysicalGroup("Edge", 4, 0),1.e-3)