[GenericResidualLaw] - using on geometry C10

dG3D/benchmarks/GenericResi/MaterialProperties.csv

+0.00000101;1374.26432899812;0.282603662323953;7.89098459151026;10.8430874308007;16603.3622833876;0.580533151556069;902.792150484622;129.853782892108;0.835450700887024;7903.61824701311;51.0754184856761;1399.632483218;0.858564039808116;11.6372096193557;7.66712095541471;2.45515873181254;21.5787270812983;6.52802930750955;1520.01598877666;8.9652462603245;152249.139847035;0.179387657001735;3.70474938213917;0.217536394479992;29.6963765349203;0.00293378675867645;14.3658751637066;299.09400060762;0.77710464112875;23.0771920299935;139.528343366004;0.000654827439341387;9648.94772906659;0.7;0.49;0;0;0.1;3;0.00725;0.08;7.5;0.07; 2.0

dG3D/benchmarks/GenericResi/genericRM.py

@@ -20,8 +20,8 @@ resfile = "USF_Results.csv"
 
 # geometry
 # geofile="cube.geo"
-# meshfile="150001_10_128.msh" # name of mesh file
-meshfile="cube.msh" # name of mesh file
+meshfile="150001_10_128.msh" # name of mesh file
+# meshfile="cube.msh" # name of mesh file
 # solver
 sol = 2  # Gmm=0 (default) Taucs=1 PETsc=2
 soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
@@ -57,8 +57,8 @@ law1.setLawForFRes(_cp);
 
 pertLaw1 = dG3DMaterialLawWithTangentByPerturbation(law1,1.e-8)
 
-nfield = 10 # number of the field (physical number of surface)
-# nfield = 51
+# nfield = 10 # number of the field (physical number of surface)
+nfield = 51
 # myfield1 = dG3DDomain(1000,nfield,space1,lawnum+1,fullDg,3,0,1)
 
 myfield1 = dG3DDomain(1000,nfield,space1,lawnum+1,fullDg,3,0,0)
@@ -92,22 +92,22 @@ mysolver.options("-pc_type lu")
 
 # # mysolver.displacementBC("Face",5678,0,0.0)
 # # mysolver.displacementBC("Face",5678,1,0.0)
-mysolver.displacementBC("Face",5678,2,0.)
-mysolver.displacementBC("Face",2376,0,0.)
-mysolver.displacementBC("Face",1265,1,0.)
+# mysolver.displacementBC("Face",5678,2,0.)
+# mysolver.displacementBC("Face",2376,0,0.)
+# mysolver.displacementBC("Face",1265,1,0.)
 """
 mysolver.displacementBC("Volume",nfield,0,0.0)
 mysolver.displacementBC("Volume",nfield,1,0.0)
 mysolver.displacementBC("Volume",nfield,2,0.0)
 """
 
-mysolver.displacementBC("Face",1234,2,-1.e-10) #@Mohib
+# mysolver.displacementBC("Face",1234,2,-1.e-10) #@Mohib
 
-# mysolver.displacementBC("Face",100,0,0.)
-# mysolver.displacementBC("Face",110,1,0.)
-# mysolver.displacementBC("Face",120,2,0.)
+mysolver.displacementBC("Face",100,0,0.)
+mysolver.displacementBC("Face",110,1,0.)
+mysolver.displacementBC("Face",120,2,0.)
 
-# mysolver.displacementBC("Face",121,2,-0.0001)
+mysolver.displacementBC("Face",121,2,-0.0001)
 
 
 #thermal BC

dG3D/benchmarks/GenericResi/genericRM_VEVP.py

+#coding-Utf-8-*-
+from gmshpy import *
+#from dG3Dpy import*
+from dG3DpyDebug import*
+
+
+# CONFIG: GP level NR Solver Tollerance
+GP_TOLL = 1.e-4
+# CONFIG: Log approximation order for HENKY Strains
+log_order = 5
+# CONFIG: GP level sub stepping  
+n_GP_subSteps = 5
+# CONFIG: FE level NR Solver Tollerance
+tol=1.e-4 # 1.e-5
+# CONFIG: Point wise time manager tollerances
+PWT_toll = 1.e-12
+PWT_save_toll = 1.e-9
+
+system = 2
+#script to launch beam problem with a python script
+
+resfile = "USF_Results.csv"
+
+mat = open("MaterialProperties.csv", 'r')
+params = mat.readline().split(';')
+
+# VEVP Law parameters
+rho = float(params[0])
+E = float(params[1]) 
+nu = float(params[2]) 
+g1 = float(params[3]) 
+k1 = float(params[4]) 
+kk1 = float(params[5])
+kk2 = float(params[6])
+kk3 = float(params[7])
+kk4 = float(params[8])
+kk5 = float(params[9])
+kk6 = float(params[10])
+kk7 = float(params[11])
+kk8 = float(params[12])
+gg1 = float(params[13])
+gg2 = float(params[14])
+gg3 = float(params[15])
+gg4 = float(params[16])
+gg5 = float(params[17])
+gg6 = float(params[18])
+gg7 = float(params[19])
+gg8 = float(params[20])
+eta = float(params[21])
+s = float(params[22])
+alpha = float(params[23])
+beta = 1.5*float(params[24]) 
+sy0c = float(params[25])
+c0 = float(params[26])
+hhc0 = float(params[27])
+hc = float(params[28])
+m = float(params[29])
+sy0t = float(params[30])
+t0 = float(params[31])
+hht0 = float(params[32])
+ht = float(params[33])
+
+# Saturation law parameters
+l = float(params[34])
+cl = float(params[35])
+H = float(params[36])
+Dinf = float(params[37])
+pi = float(params[38])
+alpha2 = float(params[39])
+
+# Failure law parameters
+pf = float(params[40]) # 0.117#0.01
+af = float(params[41])#0.55
+bf = float(params[42])#0.55
+cf = float(params[43])
+
+# control flag 
+d_flag = float(params[44])
+
+mat.close()
+
+# material law
+lawnum1 = 11 # unique number of law
+
+# Isotropic Hardening law - Compression
+hardenc = LinearExponentialJ2IsotropicHardening(1, sy0c, c0, hhc0, hc)
+# Isotropic Hardening law - Tension
+hardent = LinearExponentialJ2IsotropicHardening(2, sy0t, t0, hht0, ht)
+
+law1   = HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw(lawnum1,rho,E,nu, GP_TOLL)
+law1.setCompressionHardening(hardenc)
+law1.setTractionHardening(hardent)
+
+# CONFIG: Log approximation order for HENKY Strains
+law1.setStrainOrder(log_order)
+
+# Yield Criterion and flow rule for VP
+etac = constantViscosityLaw(1,eta)
+law1.setViscosityEffect(etac,s)
+
+law1.setYieldPowerFactor(alpha)
+law1.setNonAssociatedFlow(True)
+law1.nonAssociatedFlowRuleFactor(beta)
+
+# VE related parameters
+law1.setViscoelasticMethod(0)
+law1.setViscoElasticNumberOfElement(8)
+
+# VE Bulk props
+law1.setViscoElasticData_Bulk(1, kk1, k1)
+law1.setViscoElasticData_Bulk(2, kk2, k1*0.1)
+law1.setViscoElasticData_Bulk(3, kk3, k1*0.01)
+law1.setViscoElasticData_Bulk(4, kk4, k1*0.001)
+law1.setViscoElasticData_Bulk(5, kk5, k1*0.0001)
+law1.setViscoElasticData_Bulk(6, kk6, k1*0.00001)
+law1.setViscoElasticData_Bulk(7, kk7, k1*0.000001)
+law1.setViscoElasticData_Bulk(8, kk8, k1*0.0000001)
+
+# VE Shear Props
+law1.setViscoElasticData_Shear(1, gg1, g1)
+law1.setViscoElasticData_Shear(2, gg2, g1*0.1)
+law1.setViscoElasticData_Shear(3, gg3, g1*0.01)
+law1.setViscoElasticData_Shear(4, gg4, g1*0.001)
+law1.setViscoElasticData_Shear(5, gg5, g1*0.0001)
+law1.setViscoElasticData_Shear(6, gg6, g1*0.00001)
+law1.setViscoElasticData_Shear(7, gg7, g1*0.000001)
+law1.setViscoElasticData_Shear(8, gg8, g1*0.0000001)
+
+
+law2   = GenericResidualMechanicsDG3DMaterialLaw(lawnum1+1, rho, resfile, True)
+law2.setMechanicalMaterialLaw(law1);
+# law1.setApplyReferenceF(False);
+_cp = constantScalarFunction(1.0);
+law2.setLawForFRes(_cp);
+pertLaw1 = dG3DMaterialLawWithTangentByPerturbation(law2,1.e-8)
+
+
+# geometry
+# geofile="cube.geo"
+meshfile="150001_10_128.msh" # name of mesh file
+# meshfile="cube.msh" # name of mesh file
+# solver
+sol = 2  # Gmm=0 (default) Taucs=1 PETsc=2
+soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
+nstep = 4000   # number of step (used only if soltype=1)
+ftime =4600.   # Final time (used only if soltype=1)
+nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)
+fullDg = 0 #O = CG, 1 = DG
+space1 = 0 # function space (Lagrange=0)
+beta1  = 10.
+eqRatio =1.e8
+
+
+# nfield = 10 # number of the field (physical number of surface)
+nfield = 51
+# myfield1 = dG3DDomain(1000,nfield,space1,lawnum+1,fullDg,3,0,1)
+
+myfield1 = dG3DDomain(1000,nfield,space1,lawnum1+1,fullDg,3,0,0)
+
+myfield1.matrixByPerturbation(1,1,1,1.e-8)
+myfield1.stabilityParameters(beta1)
+# CONFIG: GP level sub stepping
+myfield1.strainSubstep(2, n_GP_subSteps)
+# myfield1.setConstitutiveExtraDofDiffusionEqRatio(eqRatio) 
+# myfield1.setConstitutiveExtraDofDiffusionAccountSource(fieldSource, mecaSource)
+
+#myfield1.setConstitutiveExtraDofDiffusionStabilityParameters(beta1,bool(fullDg))
+
+# creation of Solver
+mysolver = nonLinearMechSolver(1000)
+# mysolver.createModel(geofile,meshfile,3,2)
+mysolver.loadModel(meshfile)
+mysolver.addDomain(myfield1)
+mysolver.addMaterialLaw(law2)
+mysolver.addMaterialLaw(pertLaw1)
+mysolver.Scheme(soltype)
+mysolver.Solver(sol)
+mysolver.snlData(nstep,ftime,tol,tol/100.)
+mysolver.stepBetweenArchiving(nstepArch)
+mysolver.options("-ksp_type preonly")
+mysolver.options("-pc_type lu")
+
+
+# CONFIG: Create Point wise time manager and assign tollerances
+timePlan = PointwiseTimeManager(PWT_toll, PWT_save_toll)
+
+ttotal = ftime/nstep
+for i in range(nstep):
+  ttotal = ttotal + i*ftime/nstep
+  timePlan.setNumStepTimeInterval(ttotal,1)
+  timePlan.put(ttotal)
+
+# Assign solver level time stepping plan
+mysolver.setTimeManager(timePlan)
+
+system = system 
+control = 0  
+
+mysolver.setSystemType(system)
+mysolver.setControlType(control)
+
+
+mysolver.options("-ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps -mat_mumps_icntl_24 1 -mat_mumps_icntl_13 1 -mat_mumps_icntl_23 5000")
+
+mysolver.setMessageView(True)
+
+microBC = nonLinearPeriodicBC(1000,3)
+
+microBC.setOrder(1)
+microBC.setBCPhysical(100,110,120,101,111,121)
+method = 0	#[0] Periodic mesh = 0
+degree = 3 # [3]Order used for polynomial interpolation
+addvertex = 0  # Polynomial interpolation by mesh vertex = 0, Polynomial interpolation by virtual vertex
+microBC.setPeriodicBCOptions(method, degree, bool(addvertex))
+microBC.setDeformationGradient(2.,1.,1.,1.,2.,1.,1.,1.,2.)
+
+fctxx = piecewiseScalarFunction()
+fctyy = piecewiseScalarFunction()
+fctzz = piecewiseScalarFunction()
+fctxy = piecewiseScalarFunction()
+fctyx = piecewiseScalarFunction()
+fctxz = piecewiseScalarFunction()
+fctzx = piecewiseScalarFunction()
+fctyz = piecewiseScalarFunction()
+fctzy = piecewiseScalarFunction()
+
+
+ttotal = ftime/nstep
+ffx = 0.0174481263816847/nstep
+ffz = 0.300288204913361/nstep
+for i in range(nstep+1):
+  ttotal = ttotal + i*ftime/nstep
+  ffz = ffz + i * 0.300288204913361/nstep
+  fctxx.put(ttotal ,0.)
+  fctyy.put(ttotal ,0.)
+  fctzz.put(ttotal ,ffz)
+  fctxy.put(ttotal ,0.)
+  fctyx.put(ttotal ,0.)
+  fctxz.put(ttotal ,0.)
+  fctzx.put(ttotal ,0.)
+  fctyz.put(ttotal ,0.)
+  fctzy.put(ttotal ,0.)
+
+
+microBC.setPathFunctionDeformationGradient(0,0,fctxx)
+microBC.setPathFunctionDeformationGradient(1,1,fctyy)
+microBC.setPathFunctionDeformationGradient(2,2,fctzz)
+microBC.setPathFunctionDeformationGradient(0,1,fctxy)
+microBC.setPathFunctionDeformationGradient(1,0,fctyx)
+microBC.setPathFunctionDeformationGradient(0,2,fctxz)
+microBC.setPathFunctionDeformationGradient(2,0,fctzx)
+microBC.setPathFunctionDeformationGradient(1,2,fctyz)
+microBC.setPathFunctionDeformationGradient(2,1,fctzy)
+
+mysolver.addMicroBC(microBC)
+
+
+mysolver.stressAveragingFlag(True) # set stress averaging ON- 0 , OFF-1
+mysolver.setStressAveragingMethod(0) # 0 -volume 1- surface
+
+mysolver.tangentAveragingFlag(False) 
+mysolver.setTangentAveragingMethod(2,1e-6) # 0- perturbation 1- condensation
+
+
+
+
+
+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("FRes_XX",IPField.FfAM_XX, 1, 1)
+mysolver.internalPointBuildView("FRes_YY",IPField.FfAM_YY, 1, 1)
+mysolver.internalPointBuildView("FRes_ZZ",IPField.FfAM_ZZ, 1, 1)
+mysolver.internalPointBuildView("FRes_XY",IPField.FfAM_XY, 1, 1)
+mysolver.internalPointBuildView("FRes_ZX",IPField.FfAM_ZX, 1, 1)
+mysolver.internalPointBuildView("FRes_YZ",IPField.FfAM_YZ, 1, 1)
+
+
+# mysolver.internalPointBuildView("temperature",IPField.TEMPERATURE, 1, 1)
+# mysolver.internalPointBuildView("qx",IPField.THERMALFLUX_X, 1, 1)
+# mysolver.internalPointBuildView("qy",IPField.THERMALFLUX_Y, 1, 1)
+# mysolver.internalPointBuildView("qz",IPField.THERMALFLUX_Z, 1, 1)
+# mysolver.archivingForceOnPhysicalGroup("Face", 1234, 2)
+# mysolver.archivingForceOnPhysicalGroup("Face", 5678, 2)
+# mysolver.archivingForceOnPhysicalGroup("Face", 1234, 3)
+# mysolver.archivingForceOnPhysicalGroup("Face", 5678, 3)
+# mysolver.archivingNodeDisplacement(1,3,1)
+# mysolver.archivingNodeDisplacement(5,2,1)
+
+mysolver.solve()
+
+# check = TestCheck()
+# check.equal(4.431763e+08,mysolver.getArchivedForceOnPhysicalGroup("Face", 5678, 3),1.e-6)
+# check.equal(3.010568e+02,mysolver.getArchivedNodalValue(1,3,mysolver.displacement),1.e-6)
+# check.equal(0.000000e+00,mysolver.getArchivedNodalValue(5,2,mysolver.displacement),1.e-6)
+