add benchmarks for thomason

dG3D/benchmarks/CMakeLists.txt

@@ -11,6 +11,8 @@ add_subdirectory(CohesiveBand_cube)
 add_subdirectory(CohesiveBand_Plate)
 add_subdirectory(Gurson_Cube)
 add_subdirectory(Gurson_TwoHole)
+add_subdirectory(Thomason_Cube)
+add_subdirectory(Thomason_planeStrain)
 add_subdirectory(nonLocalDamageParallel)
 add_subdirectory(nonLocalMFH)
 add_subdirectory(seebeck)

dG3D/benchmarks/Thomason_cube/CMakeLists.txt

+# test file
+
+set(PYFILE cube.py)
+
+set(FILES2DELETE 
+  *.csv
+  disp*
+  stress*
+)
+
+add_cm3python_test(${PYFILE} "${FILES2DELETE}")

dG3D/benchmarks/Thomason_cube/cube.geo

+// Cube.geo
+mm=1.0e-3;	// Units
+n=1;		// Number of layers (in thickness / along z)
+m = 1;		// Number of element + 1 along y
+p = 2; 		// Number of element + 1 along x (min 3 to have a crack)
+L=2.*mm;	// Cube lenght
+sl1=L/n;
+
+// Face z = 0 (point in anti-clockwise order)
+Point(1)={0,0,0,sl1};
+Point(2)={L,0,0,sl1};
+Point(3)={L,L,0,sl1};
+Point(4)={0,L,0,sl1};
+Line(1)={1,2};
+Line(2)={2,3};
+Line(3)={3,4};
+Line(4)={4,1};
+Line Loop(5) = {3, 4, 1, 2};
+Plane Surface(6) = {5};
+
+// Cube extrusion
+Extrude {0, 0, L} {
+  Surface{6}; Layers{n}; Recombine;
+}
+
+// Physical entities
+	// Volume
+Physical Volume(29) = {1};
+	// Surface
+Physical Surface(30) = {19}; // Left face x = 0
+Physical Surface(31) = {27}; // Right face x = L
+Physical Surface(32) = {6};  // Back face z = 0
+Physical Surface(33) = {28}; // Front face z = L
+Physical Surface(34) = {15}; // Down face y = 0
+Physical Surface(35) = {23}; // Up face y = L
+
+// Mesh
+Transfinite Line {2, 4} = m Using Progression 1;
+Transfinite Line {1, 3} = p Using Progression 1;
+Transfinite Surface {6};
+Recombine Surface{6};
+
+
+// To save cube elongation
+Physical Point(41) = {1}; // Point on left face (0,0,0)
+Physical Point(42) = {4}; // Point on left face (0,L,0)
+Physical Point(43) = {2}; // Point on right face (L,0,0)
+Physical Point(44) = {3}; // Point on right face (L,L,0)
+
+
+

dG3D/benchmarks/Thomason_cube/cube.py

+#coding-Utf-8-*-
+from gmshpy import*
+#from dG3DpyDebug import*
+from dG3Dpy import*
+from math import*
+
+# Script for testing cube problem with Gurson
+
+# Material law creation
+# ===================================================================================
+BulkLawNum1 = 1 
+
+# material parameters - bulk law
+rho = 7850. 		# kg/m3:	density
+young = 210.0*1.e9	# Pa: 	module de young
+nu    = 0.3		# -: 	poisson coefficient
+sy0   = 400.0e6 	# Pa: 	initial yield stress
+h_hard  = 100.0e6	# Pa:	hardenning coefficient
+h_exp = 1.0		# -: 	exponential coefficient
+HardenLawNum1 = 2
+Harden1 = PowerLawJ2IsotropicHardening(HardenLawNum1, sy0, h_hard, h_exp)
+# - gurson
+q1    = 1.5
+q2    = 1.
+q3    = 1.
+fVinitial = 0.01
+# - nucleation
+fn     = 0.1;
+sn     = 0.1;
+epsilonn = 0.3;
+NucleationLawNum1 = 3
+Gdnexp1 = ExponentialNucleationLaw(NucleationLawNum1, fn, sn, epsilonn);
+# - non-local law
+cl = 2.0e-4 		# m2: 		non-local parameter (== NL-length^2) 
+ClLawNum1 = 4	
+Cl_Law1 = IsotropicCLengthLaw(ClLawNum1, cl)
+# material law creation
+BulkLaw1 = NonLocalPorousThomasonDG3DMaterialLaw(BulkLawNum1, young, nu, rho, fVinitial, 5., Harden1, Cl_Law1)
+BulkLaw1.setOrderForLogExp(9)
+BulkLaw1.setNucleationLaw(Gdnexp1)
+BulkLaw1.setPredictorCorrectorParameters(2,50,1.)
+BulkLaw1.setYieldSurfaceExponent(20)
+
+
+
+
+
+# Solver parameters
+# ===================================================================================
+soltype = 1 		# StaticLinear=0 (default, StaticNonLinear=1, Explicit=2, 
+			# Multi=3, Implicit=4, Eigen=5)
+nstep = 250   		# Number of step
+ftime =1.0   		# Final time
+tol=1.e-6  		# Relative tolerance for NR scheme
+tolAbs = 1.e-10		# Absolute tolerance for NR scheme
+nstepArch=nstep/250		# Number of step between 2 archiving
+nstepArchEnergy = 1	# Number of step between 2 energy computation and archiving
+nstepArchForce = 1	# Number of step between 2 force archiving
+MaxIter = 25		# Maximum number of iterations
+StepIncrease = 3	# Number of successfull timestep before reincreasing it
+StepReducFactor = 5.0 	# Timestep reduction factor
+NumberReduction = 4	# Maximum number of timespep reduction: max reduction = pow(StepReducFactor,this)
+fullDg = bool(1)        # O = CG, 1 = DG
+dgnl = bool(1)		# DG for non-local variables inside a domain (only if fullDg)
+eqRatio = 1.0e6		# Ratio between "elastic" and non-local equations 
+space1 = 0 		# Function space (Lagrange=0)
+beta1  = 30.0		# Penality parameter for DG
+
+# Domain creation
+## ===================================================================================
+numPhysVolume1 = 29 		# Number of a physical volume of the model in .geo
+numDomain1 = 1000 		# Number of the domain
+field1 = dG3DDomain(numDomain1,numPhysVolume1,space1,BulkLawNum1,fullDg,3,1)
+field1.stabilityParameters(beta1) 			# Adding stability parameters (for DG)
+field1.setNonLocalStabilityParameters(beta1,dgnl) 		# Adding stability parameters (for DG)
+field1.setNonLocalEqRatio(eqRatio)
+field1.gaussIntegration(0,-1,-1)
+#field1.matrixByPerturbation(1,1,1,1e-8) 		# Tangent computation analytically or by pertubation
+
+
+# Solver creation
+# ===================================================================================
+mysolver = nonLinearMechSolver(numDomain1) 		# Solver associated with numSolver1
+geofile="cube.geo"
+meshfile= "cube.msh" 			# name of mesh file
+mysolver.createModel(geofile,meshfile,3,1)
+#mysolver.loadModel(meshfile)		# add mesh
+mysolver.addDomain(field1) 		# add domain
+mysolver.addMaterialLaw(BulkLaw1) 	# add material law
+mysolver.Solver(2) 			# Library solving: Gmm=0 (default) Taucs=1 PETsc=2
+# solver parametrisation
+mysolver.Scheme(soltype) 			# solver scheme	 				
+mysolver.snlData(nstep,ftime,tol,tolAbs) 	# solver parameters
+mysolver.snlManageTimeStep(MaxIter,StepIncrease,StepReducFactor,NumberReduction) #timestep
+# solver archiving
+mysolver.stepBetweenArchiving(nstepArch) 	# archiving frequency
+mysolver.energyComputation(nstepArchEnergy)	# archiving frequency for energy
+mysolver.lineSearch(bool(0))		# lineSearch activation
+mysolver.options("-ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps")
+
+# Boundary conditions
+# ===============================
+tot_disp = 1.e-3 # Max disp == elongation 50 pourcent
+mysolver.displacementBC("Face",30,0,0.)		# face x = 0
+mysolver.displacementBC("Face",31,0,tot_disp) # face x = L
+mysolver.displacementBC("Face",34,1,0.)		# face y = 0
+mysolver.displacementBC("Face",32,2,0.)		# face z = 0
+
+mysolver.initialBC("Volume","Position",29,3,fVinitial)
+
+# 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("epl",IPField.PLASTICSTRAIN, 1, 1)
+mysolver.internalPointBuildView("coalescence",IPField.COALESCENCE, 1, 1)
+mysolver.internalPointBuildView("chi",IPField.LIGAMENT_RATIO, 1, 1)
+
+mysolver.archivingForceOnPhysicalGroup("Face", 31, 0, nstepArchForce)
+mysolver.archivingNodeDisplacement(43,0, nstepArchForce)
+mysolver.archivingIPOnPhysicalGroup("Volume",29, IPField.LOCAL_POROSITY,IPField.MEAN_VALUE);
+mysolver.archivingIPOnPhysicalGroup("Volume",29, IPField.CORRECTED_POROSITY,IPField.MEAN_VALUE);
+mysolver.archivingIPOnPhysicalGroup("Volume",29, IPField.DAMAGE,IPField.MEAN_VALUE);
+mysolver.archivingIPOnPhysicalGroup("Volume",29, IPField.LIGAMENT_RATIO,IPField.MEAN_VALUE);
+
+# Solving
+# ===========
+mysolver.solve()
+
+
+
+
+
+
+# Test
+# ===========
+check = TestCheck()
+check.equal(1.966555e+00,mysolver.getArchivedForceOnPhysicalGroup("Face", 31, 0),1.e-5)
+
+import csv
+
+data = csv.reader(open('IPVolume29val_LOCAL_POROSITYMean.csv'), delimiter=';')
+porosity = list(data)
+check.equal(3.412025e-01,float(porosity[-1][1]),1e-5)
+
+data = csv.reader(open('IPVolume29val_LIGAMENT_RATIOMean.csv'), delimiter=';')
+porosity = list(data)
+check.equal(9.995000e-01,float(porosity[-1][1]),1e-5)
+
+
+
+
+

dG3D/benchmarks/Thomason_planeStrain/CMakeLists.txt

+# test file
+
+set(PYFILE model.py)
+
+set(FILES2DELETE 
+  *.csv
+  disp*
+  stress*
+)
+
+add_cm3python_test(${PYFILE} "${FILES2DELETE}")

dG3D/benchmarks/Thomason_planeStrain/model.geo

+mm = 1.;
+scale = 1./3.;
+
+B = scale*7.5*mm;
+R = scale*12.5*mm;
+L = scale*32.5*mm;
+
+lsca1 = B*0.1;
+lsca2 = L*0.1;
+Point(1) = {0,0,0,lsca1};
+Point(2) = {L,0,0,lsca2};
+Point(3) = {L,B+R,0,lsca2};
+Point(4) = {B+R,B+R,0,0.7*lsca2};
+Point(5) = {B,B+R,0,0.5*lsca2};
+Point(6) = {B,B,0,0.4*lsca2};
+Point(7) = {0,B,0,lsca1};
+Point(8) = {B*0.2,0,0,lsca1};
+Point(9) = {B*0.2,B,0,lsca1};
+
+
+//+
+Symmetry {1, 0, 0, 0} {
+  Duplicata { Point{3}; Point{4}; Point{5}; Point{6}; Point{2}; Point{8}; Point{9}; }
+}
+//+
+Line(1) = {4, 3};
+//+
+Line(2) = {3, 2};
+//+
+Line(3) = {2, 8};
+//+
+Line(4) = {8, 15};
+//+
+Line(5) = {15, 14};
+//+
+Line(6) = {14, 10};
+//+
+Line(7) = {10, 11};
+//+
+Line(8) = {13, 16};
+//+
+Line(9) = {16,7};
+//+
+Line(10) = {7, 9};
+//+
+Line(11) = {9, 6};
+//+
+Circle(12) = {11, 12, 13};
+//+
+Circle(13) = {6, 5, 4};
+//+
+Line Loop(1) = {13,11, 8, 9, 10, 12, 1, 2, 3, 4, 5, 6, 7};
+//+
+Plane Surface(1) = {1};
+//+
+Physical Surface(11) = {1};
+//+
+Physical Line(1) = {5, 4, 3};
+//+
+Physical Line(2) = {2};
+//+
+Physical Line(4) = {6};
+//+
+Physical Point(1) = {14};
+//+
+Physical Point(2) = {2};
+//+
+Physical Point(3) = {3};
+//+
+Physical Point(4) = {10};
+//+
+Physical Point(5) = {7};

dG3D/benchmarks/Thomason_planeStrain/model.py

+#coding-Utf-8-*-
+from gmshpy import*
+from dG3Dpy import*
+from math import*
+
+# Script for plane notched specimen with Gurson
+
+
+# Material law creation
+# ===================================================================================
+BulkLawNum1 = 11
+# 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
+
+# hardening
+Harden1 = SwiftJ2IsotropicHardening(BulkLawNum1, sy0, h_hard, h_exp)
+
+
+# - non-local law
+l = 0.4
+cl = l*l
+Cl_Law1 = IsotropicCLengthLaw(BulkLawNum1, cl)
+
+# - nucleation
+fn     = 0.1;
+sn     = 0.;
+epsilonn = 0.1;
+NucleationLawNum1 = 3
+Gdnexp1 = ExponentialNucleationLaw(NucleationLawNum1, fn, sn, epsilonn);
+
+# - gurson
+fVinitial = 0.05 # initial porosity
+
+# material law creation
+#BulkLaw1 = NonLocalDamageGursonDG3DMaterialLaw(BulkLawNum1, rho, young, nu, q1,q2,q3,fVinitial,Harden1, Cl_Law1, 1e-6,False,1e-8)
+
+BulkLaw1 = NonLocalPorousThomasonDG3DMaterialLaw(BulkLawNum1, young, nu, rho, fVinitial, 10., Harden1, Cl_Law1)
+BulkLaw1.setOrderForLogExp(7)
+BulkLaw1.setYieldSurfaceExponent(20)
+BulkLaw1.setNucleationLaw(Gdnexp1)
+BulkLaw1.setPredictorCorrectorParameters(0,20,1.)
+BulkLaw1.setSubStepping(bool(0),3)
+
+
+
+# Solver parameters - implicit
+# ===================================================================================
+soltype = 1 		# StaticLinear=0 (default, StaticNonLinear=1, Explicit=2, # Multi=3, Implicit=4, Eigen=5)
+nstep = 10 		# Number of step
+ftime =1.0 		# Final time
+tol=1.e-6 		# Relative tolerance for NR scheme
+tolAbs = 1.e-7		# 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. 	# Timestep reduction factor
+NumberReduction = 20	# Maximum number of timespep reduction: max reduction = pow(StepReducFactor,this)
+
+fullDg = False
+eqRatio = 1.0		# Ratio between "elastic" and non-local equations 
+space1 = 0 		# Function space (Lagrange=0)
+beta1  = 100.0		# Penality parameter for DG
+
+# Domain creation
+
+meshfile= "model.msh" 	# name of mesh file
+
+## ===================================================================================
+field1 = dG3DDomain(11,11,0,BulkLawNum1,fullDg,2,1)
+#field1.setBulkDamageBlockedMethod(-1)  # not block
+#field1.forceCohesiveInsertionAllIPs(True,0)
+field1.stabilityParameters(beta1) 			# Adding stability parameters (for DG)
+field1.setNonLocalStabilityParameters(beta1,fullDg) 		# Adding stability parameters (for DG)
+field1.setNonLocalEqRatio(eqRatio)
+#field1.matrixByPerturbation(1,1,1,1e-7)
+
+# 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
+
+# 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
+
+# path following
+method = 0
+withPF = bool(1)
+mysolver.pathFollowing(withPF,method)
+if method == 0:
+	mysolver.setPathFollowingIncrementAdaptation(True,3,0.3)
+	mysolver.setPathFollowingControlType(0)
+	mysolver.setPathFollowingCorrectionMethod(0)
+	mysolver.setPathFollowingArcLengthStep(0.01)
+	mysolver.setBoundsOfPathFollowingArcLengthSteps(1.e-5,0.1);
+else:
+	# time-step adaptation by number of NR iterations
+	mysolver.setPathFollowingIncrementAdaptation(True,3,0.3)
+	mysolver.setPathFollowingLocalSteps(1e-4,1.e-4)
+	mysolver.setPathFollowingSwitchCriterion(0.)
+	mysolver.setPathFollowingLocalIncrementType(1); 
+	mysolver.setBoundsOfPathFollowingLocalSteps(1.e-6,2.e-4)
+
+# Boundary conditions
+# ===============================
+mysolver.initialBC("Face","Position",11,3,fVinitial)
+mysolver.initialBC("Face","Position",12,3,fVinitial)
+mysolver.displacementBC("Face",11,2,0.)
+mysolver.displacementBC("Edge",1,1,0.)
+mysolver.displacementBC("Edge",4,0,0.) 
+
+if withPF == True:
+	mysolver.sameDisplacementBC("Edge",2,3,0)
+	mysolver.forceBC("Node",3,0,1e4) 
+else:
+	mysolver.displacementBC("Edge",2,0,1.e-2)
+
+
+# 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_max",IPField.PLASTICSTRAIN, 1, IPField.MAX_VALUE)
+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("POROSITY_AT_COALESCENCE",IPField.POROSITY_AT_COALESCENCE, 1, 1)
+mysolver.internalPointBuildView("FAILED",IPField.FAILED, 1, 1)
+mysolver.internalPointBuildView("chi",IPField.LIGAMENT_RATIO, 1, 1)
+mysolver.internalPointBuildView("chi_max",IPField.LIGAMENT_RATIO, 1, IPField.MAX_VALUE)
+mysolver.OneUnknownBuildView("nonlocalVar",3,1)
+
+mysolver.archivingForceOnPhysicalGroup("Edge", 4, 0, nstepArchForce)
+mysolver.archivingNodeDisplacement(3,0, nstepArchForce)
+mysolver.archivingNodeDisplacement(5,1, nstepArchForce)
+
+mysolver.archivingIPOnPhysicalGroup("Face", 11, IPField.NONLOCAL_POROSITY,IPField.MAX_VALUE,nstepArchForce)
+mysolver.archivingIPOnPhysicalGroup("Face", 11, IPField.LOCAL_POROSITY,IPField.MAX_VALUE,nstepArchForce)
+mysolver.archivingIPOnPhysicalGroup("Face", 11, IPField.CORRECTED_POROSITY,IPField.MAX_VALUE,nstepArchForce)
+
+mysolver.solve()
+
+
+
+
+
+check = TestCheck()
+check.equal(-2.002791e+02,mysolver.getArchivedForceOnPhysicalGroup("Edge", 4, 0),1.e-4)
+
+
+
+
+
+
+