first setup

dG3D/benchmarks/impactSphere/VEVP.py

+#coding-Utf-8-*-
+from gmshpy import *	
+from dG3Dpy import *			
+						
+import pickle
+import numpy as np
+
+
+def FillMatPara(x_E,x_P,VEVP):    #x = [E0,nu0,K1,G1,....,tk,tg, aP]
+  tk = np.power(10.0, x_E[int(len(x_E)-3)])
+  tg = np.power(10.0, x_E[int(len(x_E)-2)])
+  MatPara = []
+  MatPara.append([x_E[0], x_E[1]])
+  Ki = []
+  Gi = []
+  N = int((len(x_E)-5)/2)
+  for i in range(N):
+    a = 1.0/np.power(10.0, i)
+    Ki.append([np.power(10.0,x_E[2+2*i]), a*tk])
+    Gi.append([np.power(10.0,x_E[2+2*i+1]), a*tg])
+  MatPara.append(Ki)
+  MatPara.append(Gi)
+ 
+  
+  if(VEVP):
+    Compression = []
+    Tension = []
+    ViscoPlastic = []
+    
+    Compression.append(x_P[0])
+    Tension.append(x_P[4])
+    for i in range(1,3):
+      Compression.append(np.power(10.0, x_P[i]))
+      Tension.append(np.power(10.0, x_P[4+i]))
+    Compression.append(x_P[3])
+    Tension.append(x_P[7])  
+      
+    for i in range(4):  
+      ViscoPlastic.append(x_P[8+i])       
+      
+    MatPara.append(Compression)
+    MatPara.append(Tension) 
+    MatPara.append(ViscoPlastic) 
+  return MatPara
+############################################################################  
+
+
+
+
+def RunCase(LoadCase,MatPara, VEVP):
+  
+  nstep = 800 # number of step (used only if soltype=1)
+  # === Problem parameters ============================================================================ 
+  E0 = MatPara[0][0]
+  nu0 = MatPara[0][1]
+  Nt = len(MatPara[1])    #the number of branches in the generalized Maxwell model 
+  aP = 1.0 
+
+  if(VEVP):
+    indxplastic = 3 
+    sy0c = MatPara[indxplastic][0]  #MPa, compressive yield stress
+    hc = MatPara[indxplastic][1]
+    hc2 = MatPara[indxplastic][2]
+    kc = MatPara[indxplastic][3]
+  
+    sy0t = MatPara[indxplastic+1][0]*sy0c   #MPa, compressive yield stress
+    ht = MatPara[indxplastic+1][1]
+    ht2 = MatPara[indxplastic+1][2]
+    kt = MatPara[indxplastic+1][3]
+   
+    alpha = MatPara[indxplastic+2][0]
+    beta = MatPara[indxplastic+2][1]
+    eta = MatPara[indxplastic+2][2]
+    p = MatPara[indxplastic+2][3]
+    
+  #=============================================    
+  lawnum = 11
+  rho = 7850e-9
+  
+  if(VEVP):
+    law1 = HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw(lawnum,rho,E0,nu0)
+  else:
+    law1 = HyperViscoElasticDG3DMaterialLaw(lawnum,rho,E0,nu0)
+      
+  law1.setViscoelasticMethod(0)
+  law1.setViscoElasticNumberOfElement(Nt)
+  for i in range(Nt):
+    law1.setViscoElasticData_Bulk(i+1, MatPara[1][i][0], aP*MatPara[1][i][1])
+    law1.setViscoElasticData_Shear(i+1, MatPara[2][i][0], aP*MatPara[2][i][1])
+  law1.setStrainOrder(5) 
+  #-----------------------------------------------------------------------------
+  if(VEVP):  
+    hardenc = LinearExponentialJ2IsotropicHardening(1, sy0c, hc, hc2, kc)
+    hardent = LinearExponentialJ2IsotropicHardening(2, sy0t, ht, ht2, kt)
+  
+    law1.setCompressionHardening(hardenc)
+    law1.setTractionHardening(hardent)
+  
+    law1.setYieldPowerFactor(alpha)
+    law1.setNonAssociatedFlow(True)
+    law1.nonAssociatedFlowRuleFactor(beta)
+    etac = constantViscosityLaw(1,eta)
+    law1.setViscosityEffect(etac,p)
+
+  
+  if(LoadCase[0] == 0):
+    ft1 = 5.732	
+    eps = -0.0157 
+    if(VEVP):
+      ft1 = 106.5
+      eps = -0.296	
+    ftime = LoadCase[2]
+  else:
+    strainrate = LoadCase[1]
+    strain_end = LoadCase[2]
+    ftime = strain_end/strainrate	# Final time (used only if soltype=1)  
+ # print("######################",strain_end,strainrate,ftime)
+  # ===Solver parameters=============================================================================
+  sol = 2  # Gmm=0 (default) Taucs=1 PETsc=2
+  soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
+  tol=1.e-5  # relative tolerance for NR scheme (used only if soltype=1)
+  nstepArch=10 # 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						
+
+  # ===Domain creation===============================================================================
+  # creation of ElasticField
+  myfield1 = dG3D1DDomain(11,11,space1,lawnum, False) #for non local
+  myfield1.setState(1) #UniaxialStrain=0, UniaxialStress =1,ConstantTriaxiality=2
+  L = 5.
+  
+  # ===Solver creation===============================================================================
+  # geometry & mesh
+  geofile = "line.geo"
+  meshfile = "line.msh"
+  
+  # creation of Solver
+  mysolver = nonLinearMechSolver(1000)
+  #mysolver.createModel(geofile,meshfile,1,1)
+  mysolver.loadModel(meshfile)
+  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
+
+# ===Boundary conditions===========================================================================
+  mysolver.displacementBC("Edge",11,1,0.)
+  mysolver.displacementBC("Edge",11,2,0.)
+  mysolver.displacementBC("Node",1,0,0.)
+  
+  disp = PiecewiseLinearFunction()
+  disp.put(0.,0.) 
+  if(LoadCase[0] != 0):
+    disp.put(ftime, LoadCase[0]*strain_end*L)
+  else:  
+    disp.put(ft1,eps*L)  
+    disp.put(ftime,eps*L) 
+    
+  mysolver.displacementBC("Node",2,0,disp)	
+  mysolver.archivingAverageValue(IPField.F_XX)
+  mysolver.archivingAverageValue(IPField.P_XX)
+
+  # ===Solving=======================================================================================
+  mysolver.SetVerbosity(0)
+  mysolver.solve()   
+  return  
+
+
+#def test():
+#  print("It works!!")
+
+
+#MatPara=   
+#RunCase(L_Case,MatPara,VEVP)    	
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

dG3D/benchmarks/impactSphere/impact.py

+#coding-Utf-8-*-
+from gmshpy import *
+from dG3Dpy import*
+# import functions from VEVP.py
+from VEVP import *
+import math
+
+
+# geometry ===================================================================
+# Input .msh
+#meshfile = "square.msh" # name of mesh file
+
+meshfile = 'impactSphere.msh'
+geofile = 'impactSphere.geo' # name of mesh file
+order = 1
+
+#parameters
+nstep = 4000   # number of step
+ftime = 1e-3  # Final time v0 = 0.1mm/s, pour la demi structure v = v0/2
+nstepArch = 4 # Number of step between 2 archiving (used only if soltype=1)
+tol = 1e-4
+
+penalty= 1.e4
+
+# material law for impactor
+EImpactor = 2e5
+nuImpactor = 0.3
+mm=1.
+RextImpactor = 100*mm
+RintImpactor = 75.*mm
+lawnum3      = 3
+VolImpactor    = math.pi*(RextImpactor**3-RintImpactor**3)*4./3./8.
+MassImpactor = 5*(1e-3/mm)/4. #in tonnes if mm=1
+rhoImpactor  = MassImpactor/VolImpactor
+velocity     = -1.*(2.*9.81*1.8)**0.5*(mm/1.e-3)
+print("Density Impactor = ", rhoImpactor, "; Velocity Impactor = ", velocity)
+
+law3 = J2LinearDG3DMaterialLaw(lawnum3,rhoImpactor,EImpactor,nuImpactor,1000.*EImpactor,10.)
+
+#VEVP Law for plates as to generate the sve: line 70 000 of the MCMC
+x=[1374.2643289981168, 0.2826036623239532, 4.220196044260534, -0.06622730551411246, -0.2361729746734328, 1.0658488572731906, 2.955587774357835, 0.8846323146785388, 2.1134546059887156, 0.3900795755310005, -0.07807917229891632, 1.3340258222804844, 3.8978259552106973, 0.8147820955100668, 1.7082119337043935, 3.1818481562348735, 3.1460140132039993, 0.952562223405303, 1.0351529596236166, 0.8971311953490693, 0.04770622785629489, 29.696376534920287, -2.5325714558968597, 1.1573320881902116, 299.0940006076202, 0.7771046411287497, 2.1446624378399077, -3.183873130566032, 9648.947729066591, 3.7047493821391746, 0.2175363944799915, 5.182554847703984, 0.17938765700173523]
+
+	
+numVE = 21
+number_of_parametrs = len(x)
+x_E = list(x)[0:numVE]
+x_P = list(x)[numVE:number_of_parametrs] 
+VEVP = True # "Compression"
+MatPara=FillMatPara(x_E,x_P,VEVP)
+
+E0 = MatPara[0][0]
+nu0 = MatPara[0][1]
+Nt = len(MatPara[1])    #the number of branches in the generalized Maxwell model 
+aP = 1.0 
+
+indxplastic = 3 
+sy0c = MatPara[indxplastic][0]  #MPa, compressive yield stress
+hc = MatPara[indxplastic][1]
+hc2 = MatPara[indxplastic][2]
+kc = MatPara[indxplastic][3]
+  
+sy0t = MatPara[indxplastic+1][0]*sy0c   #MPa, compressive yield stress
+ht = MatPara[indxplastic+1][1]
+ht2 = MatPara[indxplastic+1][2]
+kt = MatPara[indxplastic+1][3]
+   
+alpha = MatPara[indxplastic+2][0]
+beta = 1.5*MatPara[indxplastic+2][1]
+eta = MatPara[indxplastic+2][2]
+p = MatPara[indxplastic+2][3]
+    
+lawnum1 = 11
+rho = 1200e-9
+  
+law1 = HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw(lawnum1,rho,E0,nu0)
+law1.setViscoelasticMethod(0)
+law1.setViscoElasticNumberOfElement(Nt)
+for i in range(Nt):
+  law1.setViscoElasticData_Bulk(i+1, MatPara[1][i][0], aP*MatPara[1][i][1])
+  law1.setViscoElasticData_Shear(i+1, MatPara[2][i][0], aP*MatPara[2][i][1])
+law1.setStrainOrder(5) 
+
+hardenc = LinearExponentialJ2IsotropicHardening(1, sy0c, hc, hc2, kc)
+hardent = LinearExponentialJ2IsotropicHardening(2, sy0t, ht, ht2, kt)
+  
+law1.setCompressionHardening(hardenc)
+law1.setTractionHardening(hardent)
+  
+law1.setYieldPowerFactor(alpha)
+law1.setNonAssociatedFlow(True)
+law1.nonAssociatedFlowRuleFactor(beta)
+etac = constantViscosityLaw(1,eta)
+law1.setViscosityEffect(etac,p)
+if (order==1):
+  law1.setUseBarF(True)
+
+# lattice  here to change
+lawnum2=2
+law2 = HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw(lawnum2,rho*0.02,E0*0.02,nu0)
+law2.setViscoelasticMethod(0)
+law2.setViscoElasticNumberOfElement(Nt)
+for i in range(Nt):
+  law2.setViscoElasticData_Bulk(i+1, MatPara[1][i][0]*0.02, aP*MatPara[1][i][1])
+  law2.setViscoElasticData_Shear(i+1, MatPara[2][i][0]*0.02, aP*MatPara[2][i][1])
+law2.setStrainOrder(5) 
+
+hardenc2 = LinearExponentialJ2IsotropicHardening(3, sy0c*0.02, hc*0.02, hc2*0.02, kc)
+hardent2 = LinearExponentialJ2IsotropicHardening(4, sy0t*0.02, ht*0.02, ht2*0.02, kt)
+  
+law2.setCompressionHardening(hardenc2)
+law2.setTractionHardening(hardent2)
+  
+law2.setYieldPowerFactor(alpha)
+law2.setNonAssociatedFlow(True)
+law2.nonAssociatedFlowRuleFactor(beta)
+eta2 = constantViscosityLaw(2,eta)
+law2.setViscosityEffect(eta2,p)
+
+if (order==1):
+  law2.setUseBarF(True)  #maybe no need if no locking
+
+
+# data law def ============================================================
+
+# solver info ============================================================
+sol = 2  # Gmm=0 (default) Taucs=1 PETsc=2
+soltype = 4 # StaticLinear=0 (default) StaticNonLinear=1
+fullDg = 0 #O = CG, 1 = DG
+space1 = 0 # function space (Lagrange=0)
+# solver info ============================================================
+
+# creation of Domain
+#Impactor
+nfield3 = 1001 # number of the field (physical number of volume in 3D)
+myfield3 = dG3DDomain(1000,nfield3,space1,lawnum3,fullDg)
+
+#Plates
+nfield1_1 = 51 # number of the field (physical number of volume in 3D)
+myfield1_1 = dG3DDomain(1001,nfield1_1,space1,lawnum1,fullDg)
+myfield1_1.strainSubstep(2,10)
+
+nfield1_2 = 52 # number of the field (physical number of volume in 3D)
+myfield1_2 = dG3DDomain(1001,nfield1_2,space1,lawnum1,fullDg)
+myfield1_2.strainSubstep(2,10)
+
+#lattice
+nfield2 = 53 # number of the field (physical number of volume in 3D)
+myfield2 = dG3DDomain(1001,nfield2,space1,lawnum2,fullDg)
+myfield2.strainSubstep(2,10)
+
+
+
+
+# creation of ElasticField ==============================================================
+
+# Solver ==============================================================
+mysolver = nonLinearMechSolver(1002)
+mysolver.createModel(geofile,meshfile,3,order)
+#mysolver.loadModel(meshfile)
+mysolver.addDomain(myfield1_1)
+mysolver.addDomain(myfield1_2)
+mysolver.addDomain(myfield2)
+mysolver.addDomain(myfield3)
+mysolver.addMaterialLaw(law1)
+mysolver.addMaterialLaw(law2)
+mysolver.addMaterialLaw(law3)
+mysolver.Scheme(soltype)
+mysolver.Solver(sol)
+mysolver.snlData(nstep,ftime,tol,tol/100.)
+mysolver.stepBetweenArchiving(nstepArch)
+mysolver.implicitSpectralRadius(0.05)
+#mysolver.options("-ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps")
+
+# Solver ==============================================================
+
+# BC ==============================================================
+# blocked in z bottom if no rigid contact with plane
+mysolver.displacementBC("Face",100,2,0.)
+
+# blocked in y side
+mysolver.displacementBC("Face",101,1,0.)
+mysolver.displacementBC("Face",1001,1,0.)
+# blocked in x side
+mysolver.displacementBC("Face",100,0,0.)
+mysolver.displacementBC("Face",1000,0,0.)
+
+mysolver.initialBC("Volume","Velocity",1001,2,velocity)
+
+# BC ==============================================================
+
+# Contact
+
+
+defoDefoContact1=dG3DNodeOnSurfaceContactDomain(1002, 2, 200, 2,2000, penalty, 3.*mm)
+mysolver.defoDefoContactInteraction(defoDefoContact1)
+
+
+#Contact
+
+
+# save
+mysolver.internalPointBuildView("svm",IPField.SVM, 1, nstepArch)
+# save stress tensor
+mysolver.internalPointBuildView("sig_xx",IPField.SIG_XX, 1, nstepArch)
+mysolver.internalPointBuildView("sig_yy",IPField.SIG_YY, 1, nstepArch)
+mysolver.internalPointBuildView("sig_zz",IPField.SIG_ZZ, 1, nstepArch)
+mysolver.internalPointBuildView("sig_xy",IPField.SIG_XY, 1, nstepArch)
+mysolver.internalPointBuildView("sig_yz",IPField.SIG_YZ, 1, nstepArch)
+mysolver.internalPointBuildView("sig_xz",IPField.SIG_XZ, 1, nstepArch)
+# save platic strain
+mysolver.internalPointBuildView("epl",IPField.PLASTICSTRAIN, 1, nstepArch)
+# save bottom force
+mysolver.archivingForceOnPhysicalGroup("Face", 100, 2,nstepArch)
+
+#save impactor displacement, velocity etc ...
+mysolver.archivingNodeDisplacement(1002, 2,nstepArch)
+mysolver.archivingNodeVelocity(1002, 2,nstepArch);
+mysolver.archivingNodeAcceleration(1002, 2, nstepArch);
+
+
+mysolver.solve()

dG3D/benchmarks/impactSphere/impactSphere.geo

@@ -57,6 +57,9 @@ impactor()+=newv; BooleanDifference(newv) = { Volume{impactortmp(0)}; Delete; }{
 surfimpactor[]    = Abs(Boundary{ Volume{impactor(0)}; });
 Printf("surf impactor ",surfimpactor());
 
+ptimpactor[] = PointsOf{ Surface{surfimpactor(1)}; };
+Printf("pts impactor ",ptimpactor());
+
 Characteristic Length { PointsOf{ Volume{ impactor[] };  } } = LImpactor;
 
 
@@ -94,6 +97,9 @@ Physical Surface(OXYIMPACTORUP) = {surfimpactor(1)};
 IMPACTOR=1001;
 Physical Volume(IMPACTOR) = {impactor(0)};
 
+PTIMPACTOR = 1002;
+Physical Point(PTIMPACTOR) = {ptimpactor(2)};
+
 
 //rigid contact
 X_contact=-0.1*mm;