addedBenchmarksForTVE_TVP_uniaxial

dG3D/benchmarks/CMakeLists.txt

@@ -262,4 +262,5 @@ add_subdirectory(multiScaleBodyForce_TgPert)
 add_subdirectory(multiScaleBodyForce_PF)
 add_subdirectory(multiScaleBodyForce_AnaPrev)
 add_subdirectory(honeycomb_compression)
-
+add_subdirectory(nonLinearTVE_uniaxial)
+add_subdirectory(nonLinearTVP_uniaxial)

dG3D/benchmarks/nonLinearTVE_uniaxial/CMakeLists.txt

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

dG3D/benchmarks/nonLinearTVE_uniaxial/line.geo

+mm = 1.;
+L = 1*mm;
+lsca = L;
+Point(1) = {0,0,0,lsca};
+Point(2) = {L,0,0,lsca};
+
+//+
+Line(1) = {1, 2};
+//+
+Physical Line(11) = {1};
+//+
+Physical Point(1) = {1};
+//+
+Physical Point(2) = {2};
+//+
+Transfinite Line {1} = 2 Using Progression 1.03;

dG3D/benchmarks/nonLinearTVE_uniaxial/runTVE.py

+#coding-Utf-8-*-
+
+from gmshpy import *
+# from dG3Dpy import*
+from dG3DpyDebug import*
+from math import*
+import csv
+import numpy as np
+
+# 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.33
+K = E/3./(1.-2.*nu)	# Bulk mudulus
+mu = E/2./(1.+nu)	  # Shear mudulus
+rho = 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 + 23 # 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)
+
+
+law1 = NonLinearTVEDG3DMaterialLaw(lawnum,rho,E,nu,Tinitial,Alpha,KThCon,Cp) #,True,1e-8)
+
+
+law1.setStrainOrder(11)
+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)
+
+# geometry
+geofile = "line.geo"
+meshfile = "line.msh"
+
+# 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
+
+
+ThermoMechanicsEqRatio = 1.e1 # setConstitutiveExtraDofDiffusionEqRatio(ThermoMechanicsEqRatio); WHAT is this parameter?
+thermalSource = True
+mecaSource = True
+myfield1 = dG3D1DDomain(11,11,space1,lawnum,0,0,1) #for non local
+myfield1.setConstitutiveExtraDofDiffusionAccountSource(thermalSource,mecaSource)
+myfield1.stabilityParameters(beta1)
+myfield1.ThermoMechanicsStabilityParameters(beta1,bool(1))    # Added
+myfield1.setThermoMechanicsEqRatio(ThermoMechanicsEqRatio)
+
+myfield1.setState(1) #UniaxialStrain=0, UniaxialStress =1, ConstantTriaxiality=2
+L = 1.
+x0 = 0;
+A0 = 0.99
+xL = L
+AL = 1
+a = (AL - A0)/(xL-x0)
+
+crossSection = linearScalarFunction(x0,A0,a)
+myfield1.setCrossSectionDistribution(crossSection)
+
+# creation of Solver
+mysolver = nonLinearMechSolver(1000)
+mysolver.createModel(geofile,meshfile,1,1)
+mysolver.addDomain(myfield1)
+mysolver.addMaterialLaw(law1)
+mysolver.Scheme(soltype)
+mysolver.Solver(sol)
+mysolver.snlData(nstep,ftime,tol)
+mysolver.stepBetweenArchiving(nstepArch)
+
+mysolver.snlManageTimeStep(15,5,2.,20) # maximal 20 times
+
+mysolver.displacementBC("Edge",11,1,0.)
+mysolver.displacementBC("Edge",11,2,0.)
+mysolver.displacementBC("Node",1,0,0.)
+
+method = 0
+withPF = False #True
+mysolver.pathFollowing(withPF,method)
+if method == 0:
+	mysolver.setPathFollowingIncrementAdaptation(True,4)
+	mysolver.setPathFollowingControlType(0)
+	mysolver.setPathFollowingCorrectionMethod(0)
+	mysolver.setPathFollowingArcLengthStep(1e-1)
+	mysolver.setBoundsOfPathFollowingArcLengthSteps(0,0.1);
+else:
+	# 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.)
+
+if withPF:
+	mysolver.forceBC("Node",2,0,1e2)
+else:
+	# mysolver.displacementBC("Node",2,0,7e-2*L)
+    mysolver.displacementBC("Node",2,0, 0.0834) # strain rate = 8.333e-2
+
+
+# mysolver.displacementBC("Node",2,0,7e-2*L)
+mysolver.initialBC("Volume","Position",11,3,Tinitial)
+fT = LinearFunctionTime(0,Tinitial,ftime,Tinitial);    # Linear Displacement with time
+# mysolver.displacementBC("Volume",11,3,fT)
+
+"""
+mysolver.internalPointBuildView("P_xx",IPField.P_XX, 1, 1)
+mysolver.internalPointBuildView("P_yy",IPField.P_YY, 1, 1)
+mysolver.internalPointBuildView("P_zz",IPField.P_ZZ, 1, 1)
+mysolver.internalPointBuildView("P_xy",IPField.P_XY, 1, 1)
+mysolver.internalPointBuildView("P_yx",IPField.P_YX, 1, 1)
+mysolver.internalPointBuildView("P_yz",IPField.P_YZ, 1, 1)
+mysolver.internalPointBuildView("P_zy",IPField.P_ZY, 1, 1)
+mysolver.internalPointBuildView("P_xz",IPField.P_XZ, 1, 1)
+mysolver.internalPointBuildView("P_zx",IPField.P_ZX, 1, 1)
+
+mysolver.internalPointBuildView("F_xx",IPField.F_XX, 1, 1)
+mysolver.internalPointBuildView("F_yy",IPField.F_YY, 1, 1)
+mysolver.internalPointBuildView("F_zz",IPField.F_ZZ, 1, 1)
+mysolver.internalPointBuildView("F_xy",IPField.F_XY, 1, 1)
+mysolver.internalPointBuildView("F_yx",IPField.F_YX, 1, 1)
+mysolver.internalPointBuildView("F_yz",IPField.F_YZ, 1, 1)
+mysolver.internalPointBuildView("F_zy",IPField.F_ZY, 1, 1)
+mysolver.internalPointBuildView("F_xz",IPField.F_XZ, 1, 1)
+mysolver.internalPointBuildView("F_zx",IPField.F_ZX, 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.archivingForceOnPhysicalGroup("Node", 1, 0)
+mysolver.archivingNodeDisplacement(2,0)
+
+# mysolver.archivingAverageValue(IPField.P_YY)
+# mysolver.archivingAverageValue(IPField.F_YY)
+mysolver.archivingAverageValue(IPField.P_XX)
+mysolver.archivingAverageValue(IPField.F_XX)
+
+
+mysolver.solve()
+
+check = TestCheck()
+check.equal(-8.135452e+01,mysolver.getArchivedForceOnPhysicalGroup("Node", 1, 0),1.e-4)

dG3D/benchmarks/nonLinearTVP_uniaxial/CMakeLists.txt

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

dG3D/benchmarks/nonLinearTVP_uniaxial/line.geo

+mm = 1.;
+L = 1*mm;
+lsca = L;
+Point(1) = {0,0,0,lsca};
+Point(2) = {L,0,0,lsca};
+
+//+
+Line(1) = {1, 2};
+//+
+Physical Line(11) = {1};
+//+
+Physical Point(1) = {1};
+//+
+Physical Point(2) = {2};
+//+
+Transfinite Line {1} = 2 Using Progression 1.03;

dG3D/benchmarks/nonLinearTVP_uniaxial/runTVP.py

+#coding-Utf-8-*-
+
+from gmshpy import *
+# from dG3Dpy import*
+from dG3DpyDebug import*
+from math import*
+import csv
+import numpy as np
+
+# 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.33
+K = E/3./(1.-2.*nu)	# Bulk mudulus
+mu = E/2./(1.+nu)	  # Shear mudulus
+rho = 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 + 23 # 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./0.78 #MPa
+hc = 30.
+sy0t = sy0c*0.78
+ht = 30.
+
+hardenc = LinearExponentialJ2IsotropicHardening(1, sy0c, hc, 10., 100.)
+hardent = LinearExponentialJ2IsotropicHardening(2, sy0t, ht, 10., 100.)
+# hardenc.setTemperatureFunction_h1(ShiftFactor_tempfunc)
+# hardenc.setTemperatureFunction_h2(ShiftFactor_tempfunc)
+# hardenc.setTemperatureFunction_hexp(ShiftFactor_tempfunc)
+# hardent.setTemperatureFunction_h1(ShiftFactor_tempfunc)
+# hardent.setTemperatureFunction_h2(ShiftFactor_tempfunc)
+# hardent.setTemperatureFunction_hexp(ShiftFactor_tempfunc)
+
+hardenk = PolynomialKinematicHardening(3,1)
+hardenk.setCoefficients(1,30.)
+# hardenk.setTemperatureFunction_K(ShiftFactor_tempfunc)
+
+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)
+law1.setViscosityEffect(eta,0.21)
+# eta.setTemperatureFunction(ShiftFactor_tempfunc)
+law1.setIsotropicHardeningCoefficients(1.,1.,1.)
+
+# geometry
+geofile = "line.geo"
+meshfile = "line.msh"
+
+# 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
+
+
+# creation of ElasticField
+#matrix
+#myfield1 = dG3D1DDomain(11,11,space1,lawnum,0,1) #for non local
+# myfield1 = dG3D1DDomain(11,11,space1,lawnum, False) #for non local
+#myfield1.matrixByPerturbation(1,1,1,1e-8)
+
+ThermoMechanicsEqRatio = 1.e1 # setConstitutiveExtraDofDiffusionEqRatio(ThermoMechanicsEqRatio); WHAT is this parameter?
+thermalSource = True
+mecaSource = True
+# myfield1 = ThermoMechanicsDG3DDomain(1000,11,space1,lawnum,fullDg,ThermoMechanicsEqRatio,1)       # Added
+myfield1 = dG3D1DDomain(11,11,space1,lawnum,0,0,1) #for non local
+myfield1.setConstitutiveExtraDofDiffusionAccountSource(thermalSource,mecaSource)
+myfield1.stabilityParameters(beta1)
+myfield1.ThermoMechanicsStabilityParameters(beta1,bool(1))    # Added
+myfield1.setThermoMechanicsEqRatio(ThermoMechanicsEqRatio)
+
+myfield1.setState(1) #UniaxialStrain=0, UniaxialStress =1, ConstantTriaxiality=2
+L = 1.
+x0 = 0;
+A0 = 0.99
+xL = L
+AL = 1
+a = (AL - A0)/(xL-x0)
+
+crossSection = linearScalarFunction(x0,A0,a)
+myfield1.setCrossSectionDistribution(crossSection)
+
+# creation of Solver
+mysolver = nonLinearMechSolver(1000)
+mysolver.createModel(geofile,meshfile,1,1)
+mysolver.addDomain(myfield1)
+mysolver.addMaterialLaw(law1)
+mysolver.Scheme(soltype)
+mysolver.Solver(sol)
+mysolver.snlData(nstep,ftime,tol)
+mysolver.stepBetweenArchiving(nstepArch)
+
+mysolver.snlManageTimeStep(15,5,2.,20) # maximal 20 times
+
+mysolver.displacementBC("Edge",11,1,0.)
+mysolver.displacementBC("Edge",11,2,0.)
+mysolver.displacementBC("Node",1,0,0.)
+
+method = 0
+withPF = False #True
+mysolver.pathFollowing(withPF,method)
+if method == 0:
+	mysolver.setPathFollowingIncrementAdaptation(True,4)
+	mysolver.setPathFollowingControlType(0)
+	mysolver.setPathFollowingCorrectionMethod(0)
+	mysolver.setPathFollowingArcLengthStep(1e-1)
+	mysolver.setBoundsOfPathFollowingArcLengthSteps(0,0.1);
+else:
+	# 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.)
+
+if withPF:
+	mysolver.forceBC("Node",2,0,1e2)
+else:
+	# mysolver.displacementBC("Node",2,0,7e-2*L)
+    mysolver.displacementBC("Node",2,0, 0.0834) # strain rate = 8.333e-2
+
+
+# mysolver.displacementBC("Node",2,0,7e-2*L)
+mysolver.initialBC("Volume","Position",11,3,Tinitial)
+fT = LinearFunctionTime(0,Tinitial,ftime,Tinitial);    # Linear Displacement with time
+# mysolver.displacementBC("Volume",11,3,fT)
+
+"""
+mysolver.internalPointBuildView("P_xx",IPField.P_XX, 1, 1)
+mysolver.internalPointBuildView("P_yy",IPField.P_YY, 1, 1)
+mysolver.internalPointBuildView("P_zz",IPField.P_ZZ, 1, 1)
+mysolver.internalPointBuildView("P_xy",IPField.P_XY, 1, 1)
+mysolver.internalPointBuildView("P_yx",IPField.P_YX, 1, 1)
+mysolver.internalPointBuildView("P_yz",IPField.P_YZ, 1, 1)
+mysolver.internalPointBuildView("P_zy",IPField.P_ZY, 1, 1)
+mysolver.internalPointBuildView("P_xz",IPField.P_XZ, 1, 1)
+mysolver.internalPointBuildView("P_zx",IPField.P_ZX, 1, 1)
+
+mysolver.internalPointBuildView("F_xx",IPField.F_XX, 1, 1)
+mysolver.internalPointBuildView("F_yy",IPField.F_YY, 1, 1)
+mysolver.internalPointBuildView("F_zz",IPField.F_ZZ, 1, 1)
+mysolver.internalPointBuildView("F_xy",IPField.F_XY, 1, 1)
+mysolver.internalPointBuildView("F_yx",IPField.F_YX, 1, 1)
+mysolver.internalPointBuildView("F_yz",IPField.F_YZ, 1, 1)
+mysolver.internalPointBuildView("F_zy",IPField.F_ZY, 1, 1)
+mysolver.internalPointBuildView("F_xz",IPField.F_XZ, 1, 1)
+mysolver.internalPointBuildView("F_zx",IPField.F_ZX, 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.archivingForceOnPhysicalGroup("Node", 1, 0)
+mysolver.archivingNodeDisplacement(2,0)
+
+# mysolver.archivingAverageValue(IPField.P_YY)
+# mysolver.archivingAverageValue(IPField.F_YY)
+mysolver.archivingAverageValue(IPField.P_XX)
+mysolver.archivingAverageValue(IPField.F_XX)
+
+
+mysolver.solve()
+
+check = TestCheck()
+check.equal(-3.089706e+01,mysolver.getArchivedForceOnPhysicalGroup("Node", 1, 0),1.e-4)