Select Git revision
polynomialBasis.h
-
Christophe Geuzaine authoredChristophe Geuzaine authored
cube.py 7.93 KiB
#coding-Utf-8-*-
from gmshpy import*
from dG3Dpy import*
from math import*
# Script for plane notched specimen with Gurson
# Material law creation
# ===================================================================================
BulkLawNum1 = 11
rho = 7850
young = 205.4e3 #MPa
nu = 0.3
sy0 = 856.4 #MPa
h = 3399.
p1 = 0.012
n1 = 0.084
p2 = 0.056
n2 = 0.058
harden_Law1 = LinearFollowedByMultiplePowerLawJ2IsotropicHardening(BulkLawNum1,sy0,h,p1,n1,p2,n2)
# - non-local law
l = 1e3
cl = l*l
Cl_Law1 = IsotropicCLengthLaw(BulkLawNum1, cl)
# - gurson
q1 = 1.5
q2 = 1.
q3 = 1.5
fVinitial = 5e-3 # initial porosity
ffstar = (q1- sqrt(q1*q1-q3*q3))/(q3*q3) # critical porosity 0.66667
print("critical porosity = ", ffstar)
# material law creation
BulkLaw1 = NonLocalPorousCoupledDG3DMaterialLaw(BulkLawNum1, young, nu, rho, q1,q2,q3, fVinitial, 3., harden_Law1, Cl_Law1)
BulkLaw1.setOrderForLogExp(9)
BulkLaw1.setShearPorosityGrowthFactor(10.)
# - nucleation
function1 = LinearNucleationFunctionLaw(0, 0.5)
criterion1 = PorosityBasedNucleationCriterionLaw(0, 0.007, 1.)
BulkLaw1.setNucleationLaw(function1, criterion1)
BulkLaw1.setCrackTransition(True)
BulkLaw1.setPostBlockedDissipationBehavior(1)
interfaceLawNum1 = 2
lBand = 1. # m
Kp = 1e12 # N/m
cohLaw1 = PorosityCohesiveBand3DLaw(interfaceLawNum1,lBand,Kp,True)
fracLawNum1 = 6
fracLaw1 = FractureByCohesive3DLaw(fracLawNum1,BulkLawNum1,interfaceLawNum1)
# Solver parameters - implicit
# ===================================================================================
soltype = 1 # StaticLinear=0 (default, StaticNonLinear=1, Explicit=2, # Multi=3, Implicit=4, Eigen=5)
nstep = 1000 # Number of step
ftime =1.0 # Final time
tol=1.e-6 # Relative tolerance for NR scheme
tolAbs = 1.e-8 # Absolute tolerance for NR scheme
nstepArchIP=1 # Number of step between 2 archiving
nstepArchForce = 1 # Number of step between 2 force archiving
nstepArchEnergy = nstepArchForce # Number of step between 2 energy computation and archiving
MaxIter = 12 # Maximum number of iterations
StepIncrease = 2 # Number of successfull timestep before reincreasing it
StepReducFactor = 2 # Timestep reduction factor
NumberReduction = 100 # Maximum number of timespep reduction: max reduction = pow(StepReducFactor,this)
fullDg = True
eqRatio = 1.e2 # Ratio between "elastic" and non-local equations
space1 = 0 # Function space (Lagrange=0)
beta1 = 100.0 # Penality parameter for DG
# Domain creation
meshfile= "cube.msh" # name of mesh file
## ===================================================================================
field1 = dG3DDomain(11,11,0,fracLawNum1,fullDg,3,1)
field1.stabilityParameters(beta1) # Adding stability parameters (for DG)
field1.setNonLocalStabilityParameters(beta1,fullDg) # Adding stability parameters (for DG)
field1.setNonLocalEqRatio(eqRatio)
field1.setBulkDamageBlockedMethod(0)
field1.averageStrainBased(False)
# Solver creation
# ===================================================================================
mysolver = nonLinearMechSolver(1000) # Solver associated with numSolver1
mysolver.loadModel(meshfile) # add mesh
mysolver.addDomain(field1) # add domain
mysolver.addMaterialLaw(BulkLaw1) # add material law
mysolver.addMaterialLaw(cohLaw1) # add material law
mysolver.addMaterialLaw(fracLaw1) # add material law
# solver parametrisation
mysolver.Scheme(soltype) # solver scheme
mysolver.Solver(2) # Library solver: Gmm=0 (default) Taucs=1 PETsc=2
mysolver.snlData(nstep,ftime,tol,tolAbs) # solver parameters for imp. solving (numb. of step, final time and tolerance)
mysolver.snlManageTimeStep(MaxIter,StepIncrease,StepReducFactor,NumberReduction)
# solver archiving
mysolver.stepBetweenArchiving(nstepArchIP) # archiving frequency
mysolver.energyComputation(nstepArchEnergy) # archiving frequency for energy
# Boundary conditions
# ===============================
mysolver.initialBC("Volume","Position",11,3,fVinitial)
mysolver.displacementBC("Face",1,0,0.)
mysolver.displacementBC("Face",3,1,0.)
mysolver.displacementBC("Face",5,2,0.)
mysolver.sameDisplacementBC("Face",2,7,0)
mysolver.sameDisplacementBC("Face",4,7,1)
mysolver.sameDisplacementBC("Face",6,7,2)
method=0
mysolver.pathFollowing(True,method)
if method==0:
mysolver.setPathFollowingIncrementAdaptation(True,6,3)
mysolver.setPathFollowingControlType(0)
mysolver.setPathFollowingCorrectionMethod(0)
mysolver.setPathFollowingArcLengthStep(2e-1)
mysolver.setBoundsOfPathFollowingArcLengthSteps(0,5e-1);
elif method==1:
mysolver.setPathFollowingIncrementAdaptation(False,5)
mysolver.setPathFollowingLocalSteps(1e-1,1e-2)
mysolver.setPathFollowingLocalIncrementType(1);
mysolver.setPathFollowingSwitchCriterion(0)
T = 2
fact1 = 0.3; # (T-1./3.)/(T+2./3)
fact2 = 0.3; #
mysolver.pressureOnPhysicalGroupBC("Face",2,1e3,0.)
mysolver.pressureOnPhysicalGroupBC("Face",4,fact1*1e3,0.)
mysolver.pressureOnPhysicalGroupBC("Face",6,fact2*1e3,0.)
# Variable storage
# ===============================
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("damage",IPField.DAMAGE,1,1)
mysolver.internalPointBuildView("local_fV",IPField.LOCAL_POROSITY,1,1)
mysolver.internalPointBuildView("local_fV_max",IPField.LOCAL_POROSITY,1,IPField.MAX_VALUE)
mysolver.internalPointBuildView("epl",IPField.PLASTICSTRAIN, 1, 1)
mysolver.internalPointBuildView("triaxiality",IPField.STRESS_TRIAXIALITY, 1, 1)
mysolver.internalPointBuildView("LODE_PARAMETER",IPField.LODE_PARAMETER, 1, 1)
mysolver.internalPointBuildView("FIRST_PRINCIPAL_STRESS",IPField.FIRST_PRINCIPAL_STRESS, 1, 1)
mysolver.internalPointBuildView("SECOND_PRINCIPAL_STRESS",IPField.SECOND_PRINCIPAL_STRESS, 1, 1)
mysolver.internalPointBuildView("THIRD_PRINCIPAL_STRESS",IPField.THIRD_PRINCIPAL_STRESS, 1, 1)
mysolver.internalPointBuildView("COALESCENCE",IPField.COALESCENCE, 1, 1)
mysolver.internalPointBuildView("NONLOCAL_POROSITY",IPField.NONLOCAL_POROSITY, 1, 1)
mysolver.internalPointBuildView("CORRECTED_POROSITY",IPField.CORRECTED_POROSITY, 1, 1)
mysolver.internalPointBuildView("CORRECTED_POROSITY_MAX",IPField.CORRECTED_POROSITY, 1, IPField.MAX_VALUE)
mysolver.internalPointBuildView("POROSITY_AT_COALESCENCE",IPField.POROSITY_COALESCENCE_ONSET, 1, 1)
mysolver.internalPointBuildView("FAILED",IPField.FAILED, 1, 1)
mysolver.OneUnknownBuildView("nonlocalVar",3,1)
mysolver.archivingForceOnPhysicalGroup("Face", 1, 0, nstepArchForce)
mysolver.archivingForceOnPhysicalGroup("Face", 3, 1, nstepArchForce)
mysolver.archivingForceOnPhysicalGroup("Face", 5, 2, nstepArchForce)
mysolver.archivingNodeDisplacement(7,0, nstepArchForce)
mysolver.archivingNodeDisplacement(7,1, nstepArchForce)
mysolver.archivingNodeDisplacement(7,2, nstepArchForce)
mysolver.archivingIPOnPhysicalGroup("Volume", 11, IPField.NONLOCAL_POROSITY,IPField.MAX_VALUE,1)
mysolver.archivingIPOnPhysicalGroup("Volume", 11, IPField.LOCAL_POROSITY,IPField.MAX_VALUE,1)
mysolver.archivingIPOnPhysicalGroup("Volume", 11, IPField.CORRECTED_POROSITY,IPField.MAX_VALUE,1)
mysolver.archivingIPOnPhysicalGroup("Volume", 11, IPField.YIELD_POROSITY,IPField.MAX_VALUE,1)
mysolver.archivingIPOnPhysicalGroup("Volume", 11, IPField.VOLUMETRIC_PLASTIC_STRAIN,IPField.MAX_VALUE,1)
mysolver.archivingIPOnPhysicalGroup("Volume", 11, IPField.VOLUMETRIC_PLASTIC_STRAIN,IPField.MAX_VALUE,1)
mysolver.archivingIPOnPhysicalGroup("Volume", 11, IPField.SVM,IPField.MAX_VALUE,1)
mysolver.archivingIPOnPhysicalGroup("Volume", 11, IPField.LIGAMENT_RATIO,IPField.MAX_VALUE,1)
mysolver.solve()
check = TestCheck()
import csv
data = csv.reader(open('IPVolume11val_LOCAL_POROSITYMax.csv'), delimiter=';')
porosity = list(data)
check.equal(2.660465e-02,float(porosity[-1][1]),1e-6)