added3DCubeTestTVP

dG3D/README.txt

@@ -44,7 +44,7 @@ git clone https://username@gitlab.onelab.info/cm3/cm3Libraries.git cm3Libraries
 cd $CM3_DIR/cm3Libraries/dG3D && mkdir release && cd release
 
 # configuring
-cmake -DCMAKE_BUILD_TYPE=Release \
+cmake -DCMAKE_BUILD_TYPE=Release \  (ccmake ..)
 -DENABLE_MPI=ON \
 -DENABLE_FLTK=OFF \
 -DENABLE_EIGEN=OFF \
@@ -88,4 +88,3 @@ export PYTHONPATH=$CM3_DIR/cm3Libraries/dG3D/release:$CM3_DIR/cm3Libraries/dG3D/
 -ksp_type preonly
 -pc_type lu
 -pc_factor_mat_solver_type petsc
-    

dG3D/benchmarks/nonLinearTVP_cube/CMakeLists.txt

+# test file
+
+set(PYFILE cubeTVP.py)
+
+set(FILES2DELETE
+  for*.csv
+  ene*.csv
+  Nod*.csv
+  disp*
+  stress*
+  IPVol*
+)
+
+add_cm3python_test(${PYFILE} "${FILES2DELETE}")

dG3D/benchmarks/nonLinearTVP_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 = 1; 		// Number of element + 1 along x
+L=1.*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/nonLinearTVP_cube/cubeTVP.py

+#coding-Utf-8-*-
+
+from gmshpy import *
+# from dG3Dpy import*
+from dG3DpyDebug import*
+from math import*
+import csv
+import numpy as np
+import os
+
+# Load the csv data for relaxation spectrum
+with open('relaxationSpectrumBetter_N30.csv', newline='') as csvfile:
+    reader = csv.reader(csvfile, delimiter=';')
+    relSpec = np.zeros((1,))
+    i = 0
+    for row in reader:
+        if i == 0:
+            relSpec[i] = float(' '.join(row))
+        else:
+            relSpec = np.append(relSpec, float(' '.join(row)))
+        i += 1
+# print(relSpec)
+N = int(0.5*(i-1))
+relSpec[0:N+1] *= 1e-6    # convert to MPa
+
+#script to launch beam problem with a python script
+# material law
+lawnum = 11 # unique number of law
+
+E = relSpec[0] #MPa
+nu = 0.3 #0.33
+K = E/3./(1.-2.*nu)	# Bulk mudulus
+mu = E/2./(1.+nu)	  # Shear mudulus
+rho = 7850e-9 #905e-12 # Bulk mass   905 kg/m3 -> 905e-12 tonne/mm3
+Cp = rho*1900e+6 # 1900 J/kg/K -> 1900e+6 Nmm/tonne/K
+KThCon = 0.14  # 0.14 W/m/K -> 0.14 Nmm/s/mm/K
+Alpha = 0.6e-6 # 1/K
+
+Tinitial = 273.15 + 70 # K
+C1 = 136.17201827
+C2 = 256.58268053 # K
+
+# Temp Settings
+Tref = 273.15+20.
+TemFuncOpt = 0                  # 0-> Shift Factor based derivatives of Ki; 1-> Gibson/Huang based derivatives
+
+# Default tempfuncs
+E_G = 2.80717267e+09 * 1e-6
+E_R = 5.52466562e+08 * 1e-6
+Cp_G = rho*1700e+6
+Cp_R = rho*2100e+6
+Alpha_G = 0.2e-6
+Alpha_R = 1.e-6
+KThCon_G = 0.10
+KThCon_R = 1.5*0.12
+k_gibson = 3.36435338e-02
+m_gibson = 0.0
+Tg = Tref
+E_tempfunc = GlassTransitionScalarFunction(1, E_R/E_G, Tg, k_gibson, m_gibson, 0)  # can't be used
+Alpha_tempfunc = GlassTransitionScalarFunction(1, Alpha_R/Alpha_G, Tg, k_gibson, m_gibson, 0)
+Cp_tempfunc = GlassTransitionScalarFunction(1, Cp_R/Cp_G, Tg, k_gibson, m_gibson, 0)
+KThCon_tempfunc = GlassTransitionScalarFunction(1, KThCon_R/KThCon_G, Tg, k_gibson, m_gibson, 0)
+ShiftFactor_tempfunc = negExpWLFshiftFactor(C1,C2,Tref)
+
+sy0c = 10. #MPa
+hc = 300.
+
+sy0t = sy0c*0.78
+ht = 300.
+
+hardenc = LinearExponentialJ2IsotropicHardening(1, sy0c, hc, 0., 10.)
+hardent = LinearExponentialJ2IsotropicHardening(2, sy0t, ht, 0., 10.)
+hardenk = PolynomialKinematicHardening(3,1)
+hardenk.setCoefficients(1,370)
+
+law1 = NonLinearTVPDG3DMaterialLaw(lawnum,rho,E,nu,1e-6,Tinitial,Alpha,KThCon,Cp) # ,True,1e-8)
+
+law1.setCompressionHardening(hardenc)
+law1.setTractionHardening(hardent)
+law1.setKinematicHardening(hardenk)
+
+
+law1.setStrainOrder(11)
+
+law1.setYieldPowerFactor(3.5)
+law1.setNonAssociatedFlow(True)
+law1.nonAssociatedFlowRuleFactor(0.5)
+
+law1.setViscoelasticMethod(0)
+law1.setTestOption(3)
+law1.setShiftFactorConstantsWLF(C1,C2)
+law1.setReferenceTemperature(Tref)
+law1.setTempFuncOption(0)
+
+law1.setViscoElasticNumberOfElement(N)
+if N > 0:
+    for i in range(1, N+1):
+        law1.setViscoElasticData(i, relSpec[i], relSpec[i+N])
+
+law1.setTemperatureFunction_ThermalDilationCoefficient(Alpha_tempfunc)
+law1.setTemperatureFunction_ThermalConductivity(KThCon_tempfunc)
+law1.setTemperatureFunction_SpecificHeat(Cp_tempfunc)
+
+eta = constantViscosityLaw(1,3e4) #(1,3e4)
+law1.setViscosityEffect(eta,0.15) # 0.21)
+law1.setIsotropicHardeningCoefficients(1.,1.,1.)
+
+
+# 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)
+fullDg = 0 #O = CG, 1 = DG
+space1 = 0 # function space (Lagrange=0)
+beta1  = 100
+
+
+# geometry
+geofile = "cube.geo" # "line.geo"
+meshfile = "cube.msh" # "line.msh"
+nfield = 29
+
+pertLaw1 = dG3DMaterialLawWithTangentByPerturbation(law1,1e-8*(Tinitial)) # use with _pert*T0 (b/c In tangent_by_perturbation_TVP u use _pert*T0)
+ThermoMechanicsEqRatio = 1.e1 # setConstitutiveExtraDofDiffusionEqRatio(ThermoMechanicsEqRatio); WHAT is this parameter?
+thermalSource = True
+mecaSource = True
+myfield1 = ThermoMechanicsDG3DDomain(1000,nfield,space1,lawnum,fullDg,ThermoMechanicsEqRatio)       # Added
+myfield1.setConstitutiveExtraDofDiffusionAccountSource(thermalSource,mecaSource)
+myfield1.stabilityParameters(beta1)
+myfield1.ThermoMechanicsStabilityParameters(beta1,bool(1))    # Added
+
+
+# creation of Solver
+mysolver = nonLinearMechSolver(1000)
+mysolver.createModel(geofile,meshfile,3,1)
+mysolver.loadModel(meshfile)
+mysolver.addDomain(myfield1)
+mysolver.addMaterialLaw(law1)
+# mysolver.addMaterialLaw(pertLaw1)
+mysolver.Scheme(soltype)
+mysolver.Solver(sol)
+mysolver.snlData(nstep,ftime,tol)
+mysolver.stepBetweenArchiving(nstepArch)
+
+# BCs
+umax = 0.001
+fu_L = LinearFunctionTime(0,0,ftime,umax);    # Linear Displacement with time
+fu_L_2 = LinearFunctionTime(0,0,ftime,2*umax);
+fu_L_3 = LinearFunctionTime(0,0,ftime,0.6*umax);
+mysolver.displacementBC("Face",30,0,0.)		  # face x = 0   - Left Face fixed
+mysolver.displacementBC("Face",31,0,fu_L)     # face x = L   - Right Face displaced - compression direction
+mysolver.displacementBC("Face",31,1,fu_L_2)
+mysolver.displacementBC("Face",31,1,fu_L_3)
+mysolver.displacementBC("Face",34,1,0.)		  # face y = 0
+mysolver.displacementBC("Face",32,2,0.)		  # face z = 0
+
+# thermal BC
+mysolver.initialBC("Volume","Position",nfield,3,Tinitial)
+qmax = 0 #5.19e5                                # Heat Flux in kN(ms)⁻¹(mm)⁻¹  [5.19e-1 W/m²]
+fq_L = LinearFunctionTime(0,0,ftime,qmax);    # Linear Flux with time
+# mysolver.scalarFluxBC("Face",31,3,fq_L)
+fT = LinearFunctionTime(0,Tinitial,ftime,70+Tinitial);    # Linear Displacement with time
+# fT = LinearFunctionTime(0,Tinitial,ftime,Tinitial);    # Linear Displacement with time
+mysolver.displacementBC("Volume",nfield,3,fT)
+
+# Field-Output
+
+
+mysolver.internalPointBuildView("strain_zz", IPField.STRAIN_ZZ, 1, 1);
+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("epl",IPField.PLASTICSTRAIN, 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)
+
+
+# Print/Archive History Outputs
+
+# mysolver.archivingForceOnPhysicalGroup("Face", 31, 0, nstepArch)  # Force on the right force
+mysolver.archivingForceOnPhysicalGroup("Face", 30, 0, nstepArch)  # Force on the left face  === Reaction Force
+# mysolver.archivingForceOnPhysicalGroup("Node", 41, 0, nstepArch)  # Force on a point on the left face === Reaction Force
+mysolver.archivingNodeDisplacement(43, 0, nstepArch)              # Displacement of a point on the right face
+mysolver.archivingIPOnPhysicalGroup("Volume",nfield, IPField.TEMPERATURE,IPField.MEAN_VALUE);
+# mysolver.archivingForceOnPhysicalGroup("Face", 31, 3, nstepArch)
+# mysolver.archivingIPOnPhysicalGroup("Face",31, IPField.THERMALFLUX_X,IPField.MEAN_VALUE);
+
+# Solve
+mysolver.solve()
+
+check = TestCheck()
+check.equal(-1.062888e-04,mysolver.getArchivedForceOnPhysicalGroup("Face", 30, 0),1.e-4)

dG3D/benchmarks/nonLinearTVP_cube/relaxationSpectrumBetter_N30.csv

+4.611542021776787639e+08
+2.238302906718514562e+08
+6.933795631409890950e+07
+5.612462715236927569e+07
+1.135300364704844803e+08
+8.198232567065705359e+07
+8.810335379910270870e+07
+1.004036562814406157e+08
+7.433283445541293919e+07
+9.264592353174345195e+07
+5.464365636767672747e+07
+3.155253068493874371e+07
+7.192803585835915804e+07
+6.030332288085080683e+07
+6.405368888230948895e+07
+5.431639856282526255e+07
+4.148370455185782164e+07
+4.886423797939539701e+07
+3.592593423976574838e+07
+6.220685617107152194e+07
+5.109722805920466781e+07
+1.927525339279997349e+07
+6.875568394042491913e+07
+6.375761663832499832e+07
+5.653656921271673590e+07
+5.447126032106360048e+07
+3.741191818174190074e+07
+6.998646951026493311e+07
+1.580391066094151512e+07
+7.113430347675971687e+07
+3.727696331510867923e+07
+2.214070727389127921e-12
+3.787529186956463470e-11
+5.705652310353853072e-11
+4.823563685188596367e-10
+2.868546847643579377e-09
+1.190414141551846151e-08
+9.698182287690671800e-08
+5.084635471550801825e-07
+2.608335587038697655e-06
+2.125642550198157774e-05
+2.712603521128600229e-05
+1.695535259534005365e-04
+1.120501111743022913e-03
+7.704679206158343736e-03
+6.780508131648937953e-02
+3.249279128365250568e-01
+1.575706002496040981e+00
+5.640665343564683631e+00
+3.250919857455923534e+01
+1.678403534861442665e+02
+3.637830877100003022e+02
+1.844418920184439685e+03
+1.326507094467443676e+04
+1.176262576267081022e+05
+3.862514253418788430e+05
+2.833541535623365082e+06
+9.998967646092619747e+06
+6.356490149844758213e+07
+2.536265878395002186e+08
+6.714104114903782606e+08