Plane_notch.py

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

# Script for plane notched specimen with Gurson


# Material law creation
# ===================================================================================
BulkLawNum1 = 1 
# material parameters - bulk law
rho = 7600. 		# kg/m3:	density
young = 210.e3	# Pa: 	module de young
nu    = 0.3		# -: 	poisson coefficient
sy0   = 354.41 	#  	initial yield stress
h_hard  = 500. 
h_exp = 0.13		#	exponential coefficient

HardenLawNum1 = 2
Harden1 = SwiftJ2IsotropicHardening(HardenLawNum1, sy0, h_hard, h_exp)
# - gurson
q1    = 1.5
q2    = 1.
q3    = 1.5 
fVinitial = 5e-3
# - coalesence
fC     = 0.0#0.1#0.005
ff     = 1./q1#0.3#0.225
ffstar = 1./q1

NucleationLawNum1 = 3
Gdn1 = LinearNucleationLaw(NucleationLawNum1, 0.2, 0.005, 1.)

# - non-local law
l = 1.
cl = l*l
ClLawNum1 = 4	
Cl_Law1 = IsotropicCLengthLaw(ClLawNum1, cl)

# material law creation
BulkLaw1 = NonLocalPorousCoupledDG3DMaterialLaw(BulkLawNum1, young, nu, rho, q1,q2,q3, fVinitial, 1., Harden1, Cl_Law1,1e-8,False,1e-8)
BulkLaw1.setOrderForLogExp(9)
BulkLaw1.setYieldSurfaceExponent(20)
BulkLaw1.setNucleationLaw(Gdn1)

BulkLaw1.setSubStepping(True,3)

CoalesNum1 = 5
Coales1 = ThomasonCoalescenceLaw(CoalesNum1)
BulkLaw1.setCoalescenceLaw(Coales1)

voidEvoLawNum1 = 6
W0 = 1
lambda0 = 1
kappa = 1.5
voidEvoLaw1 = SpheroidalVoidStateEvolutionLawWithAspectRatioEvolution(voidEvoLawNum1,lambda0,W0,kappa)
BulkLaw1.setVoidEvolutionLaw(voidEvoLaw1)

BulkLaw1.setNonLocalMethod(2)


# Solver parameters - shift
# ===================================================================================
ftime =1.0 		# Final time
r_impl = 1.0

# Solver parameters - implicit
# ===================================================================================
soltype = 1 		# StaticLinear=0 (default, StaticNonLinear=1, Explicit=2, # Multi=3, Implicit=4, Eigen=5)
nstep = 260 		# Number of step
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 = 15		# Maximum number of iterations
StepIncrease = 2	# Number of successfull timestep before reincreasing it
StepReducFactor = 2.0 	# Timestep reduction factor
NumberReduction = 100	# Maximum number of timespep reduction: max reduction = pow(StepReducFactor,this)

fullDg = bool(0)        # O = CG, 1 = DG
dgnl = bool(0)		# DG for non-local variables inside a domain (only if fullDg)
eqRatio = 1.0e3		# Ratio between "elastic" and non-local equations 
space1 = 0 		# Function space (Lagrange=0)
beta1  = 30.0		# Penality parameter for DG

# Domain creation
## ===================================================================================
numPhysVolume1 = 700 		# Number of a physical volume of the model in .geo
numDomain1 = 1001 		# Number of the domain
field1 = dG3DDomain(numDomain1,numPhysVolume1,space1,BulkLawNum1,fullDg,2,3)
#field1.stabilityParameters(beta1) 			# Adding stability parameters (for DG)
#field1.setNonLocalStabilityParameters(beta1,dgnl) 		# Adding stability parameters (for DG)
#field1.gaussIntegration(0,-1,-1)
#field1.matrixByPerturbation(1,1,1,1e-8) 		# Tangent computation analytically or by pertubation
#field1.strainSubstep(2,10)
field1.setNonLocalEqRatio(eqRatio)


# Solver creation
# ===================================================================================
mysolver = nonLinearMechSolver(1002) 		# Solver associated with numSolver1
meshfile= "plane_notch.msh" 	# name of mesh file
geofile= "plane_notch.geo" 	# name of geo file
#mysolver.loadModel(meshfile)		# add mesh
mysolver.createModel(geofile,meshfile,3,2)
mysolver.addDomain(field1) 		# add domain
mysolver.addMaterialLaw(BulkLaw1) 	# 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*r_impl,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
#mysolver.lineSearch(bool(0))			# lineSearch activation

# path following
method = 0
withPF = True
mysolver.pathFollowing(withPF,method)
if method == 0:
	mysolver.setPathFollowingIncrementAdaptation(True,5,0.5)
	mysolver.setPathFollowingControlType(0)
	mysolver.setPathFollowingCorrectionMethod(0)
	mysolver.setPathFollowingArcLengthStep(1e-1)
	mysolver.setBoundsOfPathFollowingArcLengthSteps(0,1.);
elif method ==1:
	# time-step adaptation by number of NR iterations
	mysolver.setPathFollowingIncrementAdaptation(True,4)
	mysolver.setPathFollowingLocalSteps(1e-2,1.)
	mysolver.setPathFollowingSwitchCriterion(0.)
	mysolver.setPathFollowingLocalIncrementType(1); 
	mysolver.setBoundsOfPathFollowingLocalSteps(1.,3.)
elif method == 2:
	mysolver.setPathFollowingIncrementAdaptation(True,5,0.5)
	mysolver.setPathFollowingControlType(0)
	mysolver.setPathFollowingCorrectionMethod(0)
	mysolver.setPathFollowingArcLengthStep(1e-1)
	mysolver.setBoundsOfPathFollowingArcLengthSteps(0,1.);

# Boundary conditions
# ===============================
mysolver.displacementBC("Face",700,2,0.)
mysolver.displacementBC("Node",401,1,0.)
mysolver.displacementBC("Edge",1400,0,0.) 

if withPF == True:
	mysolver.sameDisplacementBC("Edge",400,401,0)
	mysolver.forceBC("Node",401,0,1e4) 
else:
	mysolver.displacementBC("Edge",400,0,3.)

stopCri = EndSchemeMonitoringWithZeroForceOnPhysical(0.3,"Edge",1400,0)
mysolver.endSchemeMonitoring(stopCri)

# 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("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("chi",IPField.LIGAMENT_RATIO, 1, 1)
mysolver.internalPointBuildView("chi_max",IPField.LIGAMENT_RATIO, 1, IPField.MAX_VALUE)
mysolver.internalPointBuildView("ASPECT_RATIO",IPField.ASPECT_RATIO, 1, 1)
mysolver.internalPointBuildView("SHAPE_FACTOR",IPField.SHAPE_FACTOR, 1, 1)

mysolver.internalPointBuildView("LIGAMENT_RATIO_COALESCENCE_ONSET",IPField.LIGAMENT_RATIO_COALESCENCE_ONSET, 1, IPField.MAX_VALUE)
mysolver.internalPointBuildView("ASPECT_RATIO_COALESCENCE_ONSET",IPField.ASPECT_RATIO_COALESCENCE_ONSET, 1, IPField.MAX_VALUE)
mysolver.internalPointBuildView("SHAPE_FACTOR_COALESCENCE_ONSET",IPField.SHAPE_FACTOR_COALESCENCE_ONSET, 1, IPField.MAX_VALUE)

mysolver.archivingForceOnPhysicalGroup("Edge", 1400, 0, nstepArchForce)
mysolver.archivingNodeDisplacement(401,0, nstepArchForce)
mysolver.archivingNodeDisplacementOnPhysical(0,100,1,nstepArchForce)
mysolver.archivingIPOnPhysicalGroup("Face", 700, IPField.DAMAGE,IPField.MAX_VALUE,nstepArchForce)
mysolver.archivingIPOnPhysicalGroup("Face", 700, IPField.LOCAL_POROSITY,IPField.MAX_VALUE,nstepArchForce)

mysolver.solve()

check = TestCheck()
check.equal(-8.920918e+02,mysolver.getArchivedForceOnPhysicalGroup("Edge", 1400, 0),1.e-4)