update testc ase

dG3D/benchmarks/DG_PRI6_resetBC/comparison.py

+import matplotlib.pyplot as plt
+import pandas as pd
+
+plt.figure()
+
+u = pd.read_csv("previousScheme_first/NodalDisplacementPhysical19Num4comp1_part2.csv",sep=";", header=None)
+f = pd.read_csv("previousScheme_first/force84comp1_part0.csv",sep=";", header=None)
+
+plt.plot(u.values[:,1],-f.values[:,1],"r.", label="first solve")
+
+u = pd.read_csv("previousScheme_second/NodalDisplacementPhysical19Num4comp1_part2.csv",sep=";", header=None)
+f = pd.read_csv("previousScheme_second/force84comp1_part0.csv",sep=";", header=None)
+
+plt.plot(u.values[:,1],-f.values[:,1],"g.", label="second solve")
+
+u = pd.read_csv("previousScheme/NodalDisplacementPhysical19Num4comp1_part2.csv",sep=";", header=None)
+f = pd.read_csv("previousScheme/force84comp1_part0.csv",sep=";", header=None)
+
+plt.plot(u.values[:,1],-f.values[:,1],"b.", label="third solve")
+
+u = pd.read_csv("NodalDisplacementPhysical19Num4comp1_part2.csv",sep=";", header=None)
+f = pd.read_csv("force84comp1_part0.csv",sep=";", header=None)
+
+plt.plot(u.values[:,1],-f.values[:,1],"k.", label="fourth solve")
+
+plt.xlabel("disp")
+plt.ylabel("force")
+
+plt.legend()
+plt.show()

dG3D/benchmarks/DG_PRI6_resetBC/twoHole.geo

@@ -70,3 +70,5 @@ Physical Volume(83) = {1};
 Physical Surface(84) = {37};
 //+
 Physical Surface(85) = {45};
+//+
+Physical Point(86) = {4};

dG3D/benchmarks/DG_PRI6_resetBC/twoHole.py

 #coding-Utf-8-*-
 from gmshpy import *
 from dG3Dpy import*
-
+import os
 #script to launch beam problem with a python script
 
 
@@ -20,7 +20,7 @@ meshfile="twoHole.msh" # name of mesh file
 # solver
 sol = 2  # Gmm=0 (default) Taucs=1 PETsc=2
 soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
-nstep = 25   # number of step (used only if soltype=1)
+nstep = 20   # number of step (used only if soltype=1)
 ftime =1.   # Final time (used only if soltype=1)
 tol=1.e-6 # relative tolerance for NR scheme (used only if soltype=1)
 nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)
@@ -34,7 +34,7 @@ law1.setUseBarF(True)
 
 # creation of ElasticField
 beta1 = 1e2
-fullDG = True;
+fullDG = False;
 averageStrainBased = False
 
 myfield1 = dG3DDomain(1000,83,0,lawnum1,fullDG,3)
@@ -61,7 +61,7 @@ mysolver.displacementBC("Face",84,0,0.)
 mysolver.displacementBC("Face",84,1,0.)
 mysolver.displacementBC("Face",84,2,0.)
 mysolver.displacementBC("Face",85,0,0.)
-mysolver.displacementBC("Face",85,1,5e-4)
+mysolver.displacementBC("Face",85,1,3e-4)
 mysolver.displacementBC("Face",85,2,0.)
 
 
@@ -76,30 +76,75 @@ mysolver.internalPointBuildView("epl",IPField.PLASTICSTRAIN, 1, 1)
 mysolver.internalPointBuildView("Damage",IPField.DAMAGE, 1, 1)
 
 mysolver.archivingForceOnPhysicalGroup("Face", 84, 1)
+mysolver.archivingForceOnPhysicalGroup("Face", 85, 1)
 mysolver.archivingNodeDisplacement(19,1,1)
+mysolver.archivingNodeDisplacement(86,1,1)
 
-mysolver.createRestartBySteps(15);
-
-t1 = mysolver.solve()
+mysolver.solve()
 
 check = TestCheck()
-check.equal(-2.720316e+03,mysolver.getArchivedForceOnPhysicalGroup("Face", 84, 1),1.e-4)
+check.equal(-2.612695e+03,mysolver.getArchivedForceOnPhysicalGroup("Face", 84, 1),1.e-4)
 
 nstep = 5   # number of step (used only if soltype=1)
 ftime =1.   # Final time (used only if soltype=1)
 tol=1.e-6 # relative tolerance for NR scheme (used only if soltype=1)
 mysolver.snlData(nstep,ftime,tol)
 
-fct = LinearFunctionTime(0,5e-4,ftime,0.)
+fval = mysolver.getArchivedForceOnPhysicalGroup("Face", 85, 1)
+print("force")
+A = 2e-5
+fct = LinearFunctionTime(0,fval/A,ftime,0.)
 mysolver.resetBoundaryConditions();
 
+mysolver.displacementBC("Face",84,0,0.)
+mysolver.displacementBC("Face",84,1,0.)
+mysolver.displacementBC("Face",84,2,0.)
+mysolver.displacementBC("Face",85,0,0.)
+mysolver.sameDisplacementBC("Face",85,86,1)
+mysolver.forceBC("Face",85,1,fct)
+mysolver.displacementBC("Face",85,2,0.)
+mysolver.solve()
+
+os.system("mv previousScheme previousScheme_first")
+
+
+nstep = 30   # number of step (used only if soltype=1)
+ftime =1.   # Final time (used only if soltype=1)
+tol=1.e-6 # relative tolerance for NR scheme (used only if soltype=1)
+mysolver.snlData(nstep,ftime,tol)
+mysolver.resetBoundaryConditions();
+fval = mysolver.getArchivedNodalValue(86,1, mysolver.displacement)
+fct = LinearFunctionTime(0,fval,ftime,6e-4)
 mysolver.displacementBC("Face",84,0,0.)
 mysolver.displacementBC("Face",84,1,0.)
 mysolver.displacementBC("Face",84,2,0.)
 mysolver.displacementBC("Face",85,0,0.)
 mysolver.displacementBC("Face",85,1,fct)
 mysolver.displacementBC("Face",85,2,0.)
+
 mysolver.solve()
 
-check = TestCheck()
-check.equal(2.631290e+03,mysolver.getArchivedForceOnPhysicalGroup("Face", 84, 1),1.e-4)
+os.system("mv previousScheme previousScheme_second")
+
+nstep = 5   # number of step (used only if soltype=1)
+ftime =1.   # Final time (used only if soltype=1)
+tol=1.e-6 # relative tolerance for NR scheme (used only if soltype=1)
+mysolver.snlData(nstep,ftime,tol)
+
+fval = mysolver.getArchivedForceOnPhysicalGroup("Face", 85, 1)
+print("force")
+A = 2e-5
+fct = LinearFunctionTime(0,fval/A,ftime,0.)
+mysolver.resetBoundaryConditions();
+
+mysolver.displacementBC("Face",84,0,0.)
+mysolver.displacementBC("Face",84,1,0.)
+mysolver.displacementBC("Face",84,2,0.)
+mysolver.displacementBC("Face",85,0,0.)
+mysolver.sameDisplacementBC("Face",85,86,1)
+mysolver.forceBC("Face",85,1,fct)
+mysolver.displacementBC("Face",85,2,0.)
+mysolver.solve()
+
+
+