first version

HyperElasticity/Integration_Lib.pro

+/*
+  Integration  method
+    GradGrad
+    CurlCurl
+*/
+
+
+Integration {
+  { Name GradGrad ;
+    Case { {Type Gauss ;
+        Case {
+          { GeoElement Point       ; NumberOfPoints  1 ; }
+          { GeoElement Line        ; NumberOfPoints  4 ; }
+          { GeoElement Triangle    ; NumberOfPoints  4 ; }
+          { GeoElement Quadrangle  ; NumberOfPoints  4 ; }
+          { GeoElement Tetrahedron ; NumberOfPoints  4 ; }
+          { GeoElement Hexahedron  ; NumberOfPoints  6 ; }
+          { GeoElement Prism       ; NumberOfPoints  9 ; } }
+      }
+    }
+  }
+  { Name CurlCurl ;
+    Case { {Type Gauss ;
+        Case {
+          { GeoElement Point       ; NumberOfPoints  1 ; }
+          { GeoElement Line        ; NumberOfPoints  4 ; }
+          { GeoElement Triangle    ; NumberOfPoints  4 ; }
+          { GeoElement Quadrangle  ; NumberOfPoints  4 ; }
+          { GeoElement Tetrahedron ; NumberOfPoints  4 ; }
+          { GeoElement Hexahedron  ; NumberOfPoints  6 ; }
+          { GeoElement Prism       ; NumberOfPoints  9 ; } }
+      }
+    }
+  }
+}
+
+/* --------------------------------------------------------------------------*/

HyperElasticity/Jacobian_Lib.pro

+/*
+  Jacobian methods
+    Vol
+*/
+
+/* I N P U T
+   ---------
+
+  GlobalGroup :
+  -----------
+    DomainInf                Regions with Spherical Shell Transformation
+
+  Parameters :
+  ----------
+    Val_Rint, Val_Rext       Inner and outer radius of the Spherical Shell
+                             of DomainInf
+
+*/
+
+/* --------------------------------------------------------------------------*/
+
+
+
+
+
+Group {
+  DefineGroup[ DomainInf ] ;
+  DefineVariable[ Val_Rint, Val_Rext ] ;
+}
+
+
+
+/* --------------------------------------------------------------------------*/
+
+
+
+
+Jacobian {
+  { Name Vol ;
+    Case { { Region DomainInf ;
+             Jacobian VolSphShell {Val_Rint, Val_Rext} ; }
+           { Region All ; Jacobian Vol ; }
+    }
+  }
+  { Name JSur ;
+    Case {
+      { Region All ; Jacobian Sur ; }
+    }
+  }
+  { Name JLin ;
+    Case {
+      { Region All ; Jacobian Lin ; }
+    }
+  }
+}
+
+
+
+/* --------------------------------------------------------------------------*/

HyperElasticity/Total_Lagrangian.pro

+//=====================================================================
+
+FunctionSpace {
+  For i In {1:3}
+    { Name H_u~{i} ; Type Form0 ;
+      BasisFunction {
+        { Name sn~{i} ; NameOfCoef un~{i} ; Function BF_Node ;
+          Support Region[ { Domain_Mecha} ] ; Entity NodesOf[ All ] ; }
+      }
+      Constraint {
+      { NameOfCoef un~{i} ;
+        EntityType NodesOf ; NameOfConstraint Displacement~{i} ; }
+      }
+    }
+  EndFor    
+}
+
+Formulation {
+  { Name Total_Lagrangian ; Type FemEquation ;
+    Quantity {
+      For i In {1:3}
+        { Name u~{i}  ; Type Local ; NameOfSpace H_u~{i} ; }
+      EndFor
+    }
+    Equation {
+      // The right-hand side
+      //===================
+      For i In {1:3}
+        Galerkin { [ P~{i}[ (TensorDiag[1., 1., 1.] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] , {d u~{i}} ] ;
+          In Domain_Mecha; Jacobian Vol; Integration GradGrad; }
+        For j In {1:3}          
+          Galerkin { [ C_Tot~{j}~{i}[ (TensorDiag[1.0, 1.0, 1.0] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] * Dof{d u~{i}}, {d u~{j}} ] ;
+            In Domain_Mecha ; Jacobian Vol ; Integration GradGrad ; }
+          Galerkin { [-C_Tot~{j}~{i}[ (TensorDiag[1.0, 1.0, 1.0] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] * {d u~{i}}, {d u~{j}} ] ;
+            In Domain_Mecha ; Jacobian Vol ; Integration GradGrad ; }
+        EndFor
+      EndFor
+          
+        Galerkin { [ - CompX[ -force_ext[] ] , {u~{1}} ] ;
+        In Domain_Mecha; Jacobian Vol; Integration GradGrad; }
+        Galerkin { [ - CompY[ -force_ext[] ] , {u~{2}} ] ;
+          In Domain_Mecha; Jacobian Vol; Integration GradGrad; }
+        Galerkin { [ - CompZ[ -force_ext[] ] , {u~{3}} ] ;
+          In Domain_Mecha; Jacobian Vol; Integration GradGrad; }
+
+    }
+  }
+}
+
+PostProcessing {
+  { Name Total_Lagrangian ; NameOfFormulation Total_Lagrangian ;
+    NameOfSystem Sys_Mec;
+    PostQuantity {
+      { Name u ; Value{ Term{ [ Vector[ {u~{1}}, {u~{2}}, {u~{3}} ]];
+	    In Domain_Mecha ; Jacobian Vol ; } } }
+      { Name um ; Value{ Term{ [ Norm[ Vector[ {u~{1}}, {u~{2}}, {u~{3}} ] ] ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+      { Name ux  ; Value{ Term{ [ {u~{1}} ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+      { Name uy  ; Value{ Term{ [ {u~{2}} ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+      { Name uz  ; Value{ Term{ [ {u~{3}} ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+
+      { Name PK2_xx; Value{ Term { [ CompXX[PK2[ (TensorDiag[1.0, 1.0, 1.0] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] ] ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+      { Name PK2_xy; Value{ Term { [ CompXY[PK2[ (TensorDiag[1.0, 1.0, 1.0] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] ] ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+      { Name PK2_xz; Value{ Term { [ CompXZ[PK2[ (TensorDiag[1.0, 1.0, 1.0] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] ] ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+      { Name PK2_yx; Value{ Term { [ CompYX[PK2[ (TensorDiag[1.0, 1.0, 1.0] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] ] ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+      { Name PK2_yy; Value{ Term { [ CompYY[PK2[ (TensorDiag[1.0, 1.0, 1.0] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] ] ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+      { Name PK2_yz; Value{ Term { [ CompYZ[PK2[ (TensorDiag[1.0, 1.0, 1.0] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] ] ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+      { Name PK2_zx; Value{ Term { [ CompZX[PK2[ (TensorDiag[1.0, 1.0, 1.0] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] ] ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+      { Name PK2_zy; Value{ Term { [ CompZY[PK2[ (TensorDiag[1.0, 1.0, 1.0] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] ] ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+      { Name PK2_zz; Value{ Term { [ CompZZ[PK2[ (TensorDiag[1.0, 1.0, 1.0] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] ] ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }        
+      
+      { Name F_tot_xx  ; Value { Term { [ CompXX[(TensorDiag[1.0, 1.0, 1.0] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+      { Name F_tot_xy  ; Value { Term { [ CompXY[(TensorDiag[1.0, 1.0, 1.0] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+      { Name F_tot_xz  ; Value { Term { [ CompXZ[(TensorDiag[1.0, 1.0, 1.0] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+      { Name F_tot_yx  ; Value { Term { [ CompYX[(TensorDiag[1.0, 1.0, 1.0] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+      { Name F_tot_yy  ; Value { Term { [ CompYY[(TensorDiag[1.0, 1.0, 1.0] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+      { Name F_tot_yz  ; Value { Term { [ CompYZ[(TensorDiag[1.0, 1.0, 1.0] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+      { Name F_tot_zx  ; Value { Term { [ CompZX[(TensorDiag[1.0, 1.0, 1.0] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+      { Name F_tot_zy  ; Value { Term { [ CompZY[(TensorDiag[1.0, 1.0, 1.0] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+      { Name F_tot_zz  ; Value { Term { [ CompZZ[(TensorDiag[1.0, 1.0, 1.0] + TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ) ] ] ;
+            In Domain_Mecha  ; Jacobian Vol; } } }
+    }
+  }
+}

HyperElasticity/beam.pro

+DefineConstant[
+  visu = {1, Choices{0, 1}, AutoCheck 0,
+    Name "Input/Solver/Visu", Label "Real-time visualization"}
+];
+
+Group {
+  Left = Region[{1}]; 
+  Top = Region[{2}]; 
+  Right = Region[{3}];
+  Domain_Mecha = Region[{100}]; // Mechanical domain
+}
+
+Flag_ExternalForce = 0;
+
+Function {
+  DefineFunction[ FT_F, PK2, C_Tot_xx, C_Tot_xy, C_Tot_xz, C_Tot_yx, C_Tot_yy, C_Tot_yz, C_Tot_zx, C_Tot_zy, C_Tot_zz,  E_g, nu_g];
+  
+  // 1. Parameters of the problem
+  //=============================
+  E  = 771.0e6; 
+  nu = 0.29; 
+
+  //Lame parameters
+  //===============
+  lambda = E/(2.0 * (1.0 + nu)); // First Lame parameter
+  mu = (nu * E)/((1.0 + nu) * (1.0 - 2 * nu)); // Second Lame parameter
+
+  //Other material properties
+  //=========================
+  rho_0 = 2700; // for Alu
+  cp_0 = 897; // for Alu
+
+  // Parameters used during resolution
+  //==================================
+  t_0 = 0.0;
+  t_end = 4.0e-0;
+  period = (t_end - t_0);
+  Freq = 1.0/period;
+  tol_abs[] = 1e-8;
+  tol_rel[] = 1e-8;
+  iter_max = 50;
+  dt_max[] = $var_Dt_Max;
+  dtime[] = $var_DTime;
+  v_TS[] = $var_TS;
+  num_steps = 32;
+  theta_value = 1.0;
+  
+  // Definition of parameters of the Dirichlet BCs, the dynamic problem and the external force
+  //==========================================================================================
+
+    force_amplitude = 1.8e03;
+    force_ext[] = Vector[0., force_amplitude, 0.0];    
+}
+
+//Include "smp_stent_KMat_Elasticity.pro";
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Function {
+
+  // 3. Kinetic information : second Piola Kirchhoff [PK2] stresses
+  //===============================================================
+  // Green Lagrange strain measure GL = E = 0.5 * (C - I) and the 2nd Piola-Kirchoff stress 
+  //=======================================================================================
+  FT_F[] = Transpose[$1] * $1; // C = (F^T * F)
+  GL[] = 0.5 * (FT_F[$1] - TensorDiag[1.0, 1.0, 1.0]);
+  GL_Trace[] = (CompXX[GL[$1]] + CompYY[GL[#1]] + CompZZ[GL[#1]]) ;
+  PK2[] = (lambda * GL_Trace[$1] * TensorDiag[1.0, 1.0, 1.0] + 2 * mu * GL[$1]); 
+
+  // 4. The nonlinear right hand side
+  //=================================
+  //RHS[] = PK2[$1] * Transpose[$1]; // The first Piola Kirchhoff : P = S * F^T
+  RHS[] = $1 * PK2[$1]; // The first Piola Kirchhoff : P = F S
+  
+  // x
+  P~{1}[] = Vector[CompXX[RHS[$1#1]], CompXY[RHS[#1]], CompXZ[RHS[#1]] ];
+  // y
+  P~{2}[] = Vector[CompYX[RHS[$1#1]], CompYY[RHS[#1]], CompYZ[RHS[#1]] ];
+  // z
+  P~{3}[] = Vector[CompZX[RHS[$1#1]], CompZY[RHS[#1]], CompZZ[RHS[#1]] ];
+
+  // 5. The stiffness matrix : C_Tot_ij = C_Mat_ij +  C_Geo_ij
+  //==========================================================
+
+  // 5.1. Geometric stiffness "C_Geo_ij" (see formula (4.69) of "Wriggers, P.,
+  // 2008. Nonlinear finite element methods. Springer Science & Business
+  // Media.")
+  //==========================================================
+  // xx
+  C_Geo~{1}~{1}[] = PK2[$1];
+  // yx
+  C_Geo~{2}~{1}[] = Tensor[ 0.0, 0.0, 0.0,   0.0, 0.0, 0.0,   0.0, 0.0, 0.0];
+  // zx
+  C_Geo~{3}~{1}[] = Tensor[ 0.0, 0.0, 0.0,   0.0, 0.0, 0.0,   0.0, 0.0, 0.0];
+  // xy
+  C_Geo~{1}~{2}[] = Tensor[ 0.0, 0.0, 0.0,   0.0, 0.0, 0.0,   0.0, 0.0, 0.0];
+  // yy
+  C_Geo~{2}~{2}[] = PK2[$1];
+  // zy
+  C_Geo~{3}~{2}[] = Tensor[ 0.0, 0.0, 0.0,   0.0, 0.0, 0.0,   0.0, 0.0, 0.0];
+  // xz
+  C_Geo~{1}~{3}[] = Tensor[ 0.0, 0.0, 0.0,   0.0, 0.0, 0.0,   0.0, 0.0, 0.0];
+  // yz
+  C_Geo~{2}~{3}[] = Tensor[ 0.0, 0.0, 0.0,   0.0, 0.0, 0.0,   0.0, 0.0, 0.0];
+  // zz
+  C_Geo~{3}~{3}[] = PK2[$1];
+
+
+
+    
+  // The material stiffness
+  //========================
+  // xx
+  C_Mat_mu~{1}~{1}[] = mu * Tensor[ 2 * CompXX[$1#1] * CompXX[#1] + CompXY[#1] * CompXY[#1]     + CompXZ[#1] * CompXZ[#1],
+                                    CompXX[#1] * CompXY[#1],
+                                    CompXX[#1] * CompXZ[#1],
+                                    CompXY[#1] * CompXX[#1],
+                                    CompXX[#1] * CompXX[#1]     + 2 * CompXY[#1] * CompXY[#1] + CompXZ[#1] * CompXZ[#1],
+                                    CompXY[#1] * CompXZ[#1],
+                                    CompXZ[#1] * CompXX[#1],
+                                    CompXZ[#1] * CompXY[#1],
+                                    CompXX[#1] * CompXX[#1]     + CompXY[#1] * CompXY[#1]     + 2 * CompXZ[#1] * CompXZ[#1] ];
+  // yx
+  C_Mat_mu~{2}~{1}[] = mu * Tensor[ 2 * CompYX[$1#1] * CompXX[#1] + CompYY[#1] * CompXY[#1]     + CompYZ[#1] * CompXZ[#1],
+                                    CompYX[#1] * CompXY[#1],
+                                    CompYX[#1] * CompXZ[#1],
+                                    CompYY[#1] * CompXX[#1],
+                                    CompYX[#1] * CompXX[#1]     + 2 * CompYY[#1] * CompXY[#1] + CompYZ[#1] * CompXZ[#1],
+                                    CompYY[#1] * CompXZ[#1],
+                                    CompYZ[#1] * CompXX[#1],
+                                    CompYZ[#1] * CompXY[#1],
+                                    CompYX[#1] * CompXX[#1]     + CompYY[#1] * CompXY[#1]     + 2 * CompYZ[#1] * CompXZ[#1] ];
+  // zx
+  C_Mat_mu~{3}~{1}[] = mu * Tensor[ 2 * CompZX[$1#1] * CompXX[#1] + CompZY[#1] * CompXY[#1]     + CompZZ[#1] * CompXZ[#1],
+                                    CompZX[#1] * CompXY[#1],
+                                    CompZX[#1] * CompXZ[#1],
+                                    CompZY[#1] * CompXX[#1],
+                                    CompZX[#1] * CompXX[#1]     + 2 * CompZY[#1] * CompXY[#1] + CompZZ[#1] * CompXZ[#1],
+                                    CompZY[#1] * CompXZ[#1],
+                                    CompZZ[#1] * CompXX[#1],
+                                    CompZZ[#1] * CompXY[#1],
+                                    CompZX[#1] * CompXX[#1]     + CompZY[#1] * CompXY[#1]     + 2 * CompZZ[#1] * CompXZ[#1] ];
+  // xy
+  C_Mat_mu~{1}~{2}[] = mu * Tensor[ 2 * CompXX[$1#1] * CompYX[#1] + CompXY[#1] * CompYY[#1]     + CompXZ[#1] * CompYZ[#1],
+                                    CompXX[#1] * CompYY[#1],
+                                    CompXX[#1] * CompYZ[#1],
+                                    CompXY[#1] * CompYX[#1],
+                                    CompXX[#1] * CompYX[#1]     + 2 * CompXY[#1] * CompYY[#1] + CompXZ[#1] * CompYZ[#1],
+                                    CompXY[#1] * CompYZ[#1],
+                                    CompXZ[#1] * CompYX[#1],
+                                    CompXZ[#1] * CompYY[#1],
+                                    CompXX[#1] * CompYX[#1]     + CompXY[#1] * CompYY[#1]     + 2 * CompXZ[#1] * CompYZ[#1] ];
+  // yy
+  C_Mat_mu~{2}~{2}[] = mu * Tensor[ 2 * CompYX[$1#1] * CompYX[#1] + CompYY[#1] * CompYY[#1]     + CompYZ[#1] * CompYZ[#1],
+                                    CompYX[#1] * CompYY[#1],
+                                    CompYX[#1] * CompYZ[#1],
+                                    CompYY[#1] * CompYX[#1],
+                                    CompYX[#1] * CompYX[#1]     + 2 * CompYY[#1] * CompYY[#1] + CompYZ[#1] * CompYZ[#1],
+                                    CompYY[#1] * CompYZ[#1],
+                                    CompYZ[#1] * CompYX[#1],
+                                    CompYZ[#1] * CompYY[#1],
+                                    CompYX[#1] * CompYX[#1]     + CompYY[#1] * CompYY[#1]     + 2 * CompYZ[#1] * CompYZ[#1] ];
+  // zy
+  C_Mat_mu~{3}~{2}[] = mu * Tensor[ 2 * CompZX[$1#1] * CompYX[#1] + CompZY[#1] * CompYY[#1]     + CompZZ[#1] * CompYZ[#1],
+                                    CompZX[#1] * CompYY[#1],
+                                    CompZX[#1] * CompYZ[#1],
+                                    CompZY[#1] * CompYX[#1],
+                                    CompZX[#1] * CompYX[#1]     + 2 * CompZY[#1] * CompYY[#1] + CompZZ[#1] * CompYZ[#1],
+                                    CompZY[#1] * CompYZ[#1],
+                                    CompZZ[#1] * CompYX[#1],
+                                    CompZZ[#1] * CompYY[#1],
+                                    CompZX[#1] * CompYX[#1]     + CompZY[#1] * CompYY[#1]     + 2 * CompZZ[#1] * CompYZ[#1] ];
+  // xz
+  C_Mat_mu~{1}~{3}[] = mu * Tensor[ 2 * CompXX[$1#1] * CompZX[#1] + CompXY[#1] * CompZY[#1]     + CompXZ[#1] * CompZZ[#1],
+                                    CompXX[#1] * CompZY[#1],
+                                    CompXX[#1] * CompZZ[#1],
+                                    CompXY[#1] * CompZX[#1],
+                                    CompXX[#1] * CompZX[#1]     + 2 * CompXY[#1] * CompZY[#1] + CompXZ[#1] * CompZZ[#1],
+                                    CompXY[#1] * CompZZ[#1],
+                                    CompXZ[#1] * CompZX[#1],
+                                    CompXZ[#1] * CompZY[#1],
+                                    CompXX[#1] * CompZX[#1]     + CompXY[#1] * CompZY[#1]     + 2 * CompXZ[#1] * CompZZ[#1] ];
+  // yz
+  C_Mat_mu~{2}~{3}[] = mu * Tensor[ 2 * CompYX[$1#1] * CompZX[#1] + CompYY[#1] * CompZY[#1]     + CompYZ[#1] * CompZZ[#1],
+                                    CompYX[#1] * CompZY[#1],
+                                    CompYX[#1] * CompZZ[#1],
+                                    CompYY[#1] * CompZX[#1],
+                                    CompYX[#1] * CompZX[#1]     + 2 * CompYY[#1] * CompZY[#1] + CompYZ[#1] * CompZZ[#1],
+                                    CompYY[#1] * CompZZ[#1],
+                                    CompYZ[#1] * CompZX[#1],
+                                    CompYZ[#1] * CompZY[#1],
+                                    CompYX[#1] * CompZX[#1]     + CompYY[#1] * CompZY[#1]     + 2 * CompYZ[#1] * CompZZ[#1] ];
+  // zz
+  C_Mat_mu~{3}~{3}[] = mu * Tensor[ 2 * CompZX[$1#1] * CompZX[#1] + CompZY[#1] * CompZY[#1]     + CompZZ[#1] * CompZZ[#1],
+                                    CompZX[#1] * CompZY[#1],
+                                    CompZX[#1] * CompZZ[#1],
+                                    CompZY[#1] * CompZX[#1],
+                                    CompZX[#1] * CompZX[#1]     + 2 * CompZY[#1] * CompZY[#1] + CompZZ[#1] * CompZZ[#1],
+                                    CompZY[#1] * CompZZ[#1],
+                                    CompZZ[#1] * CompZX[#1],
+                                    CompZZ[#1] * CompZY[#1],
+                                    CompZX[#1] * CompZX[#1]     + CompZY[#1] * CompZY[#1]     + 2 * CompZZ[#1] * CompZZ[#1] ];
+  
+  // 5.2.2. Material stiffness - the part resulting from "lambda" : C_Mat_lambda_ij
+  //===============================================================================
+  // xx
+  C_Mat_lambda~{1}~{1}[] = lambda * Tensor[ CompXX[$1#1] * CompXX[#1],   CompXX[#1] * CompXY[#1],   CompXX[#1] * CompXZ[#1],
+                                            CompXY[#1] * CompXX[#1],   CompXY[#1] * CompXY[#1],   CompXY[#1] * CompXZ[#1],
+                                            CompXZ[#1] * CompXX[#1],   CompXZ[#1] * CompXY[#1],   CompXZ[#1] * CompXZ[#1]];
+  // yx
+  C_Mat_lambda~{2}~{1}[] = lambda * Tensor[ CompYX[$1#1] * CompXX[#1],   CompYX[#1] * CompXY[#1],   CompYX[#1] * CompXZ[#1],
+                                            CompYY[#1] * CompXX[#1],   CompYY[#1] * CompXY[#1],   CompYY[#1] * CompXZ[#1],
+                                            CompYZ[#1] * CompXX[#1],   CompYZ[#1] * CompXY[#1],   CompYZ[#1] * CompXZ[#1]];
+  // zx
+  C_Mat_lambda~{3}~{1}[] = lambda * Tensor[ CompZX[$1#1] * CompXX[#1],   CompZX[#1] * CompXY[#1],   CompZX[#1] * CompXZ[#1],
+                                            CompZY[#1] * CompXX[#1],   CompZY[#1] * CompXY[#1],   CompZY[#1] * CompXZ[#1],
+                                            CompZZ[#1] * CompXX[#1],   CompZZ[#1] * CompXY[#1],   CompZZ[#1] * CompXZ[#1]];
+  // xy
+  C_Mat_lambda~{1}~{2}[] = lambda * Tensor[ CompXX[$1#1] * CompYX[#1],   CompXX[#1] * CompYY[#1],   CompXX[#1] * CompYZ[#1],
+                                            CompXY[#1] * CompYX[#1],   CompXY[#1] * CompYY[#1],   CompXY[#1] * CompYZ[#1],
+                                            CompXZ[#1] * CompYX[#1],   CompXZ[#1] * CompYY[#1],   CompXZ[#1] * CompYZ[#1]];
+  // yy
+  C_Mat_lambda~{2}~{2}[] = lambda * Tensor[ CompYX[$1#1] * CompYX[#1],   CompYX[#1] * CompYY[#1],   CompYX[#1] * CompYZ[#1],
+                                            CompYY[#1] * CompYX[#1],   CompYY[#1] * CompYY[#1],   CompYY[#1] * CompYZ[#1],
+                                            CompYZ[#1] * CompYX[#1],   CompYZ[#1] * CompYY[#1],   CompYZ[#1] * CompYZ[#1]];
+  // zy
+  C_Mat_lambda~{3}~{2}[] = lambda * Tensor[ CompZX[$1#1] * CompYX[#1],   CompZX[#1] * CompYY[#1],   CompZX[#1] * CompYZ[#1],
+                                            CompZY[#1] * CompYX[#1],   CompZY[#1] * CompYY[#1],   CompZY[#1] * CompYZ[#1],
+                                            CompZZ[#1] * CompYX[#1],   CompZZ[#1] * CompYY[#1],   CompZZ[#1] * CompYZ[#1]];
+  // xz
+  C_Mat_lambda~{1}~{3}[] = lambda * Tensor[ CompXX[$1#1] * CompZX[#1],   CompXX[#1] * CompZY[#1],   CompXX[#1] * CompZZ[#1],
+                                            CompXY[#1] * CompZX[#1],   CompXY[#1] * CompZY[#1],   CompXY[#1] * CompZZ[#1],
+                                            CompXZ[#1] * CompZX[#1],   CompXZ[#1] * CompZY[#1],   CompXZ[#1] * CompZZ[#1]];
+  // yz
+  C_Mat_lambda~{2}~{3}[] = lambda * Tensor[ CompYX[$1#1] * CompZX[#1],   CompYX[#1] * CompZY[#1],   CompYX[#1] * CompZZ[#1],
+                                            CompYY[#1] * CompZX[#1],   CompYY[#1] * CompZY[#1],   CompYY[#1] * CompZZ[#1],
+                                            CompYZ[#1] * CompZX[#1],   CompYZ[#1] * CompZY[#1],   CompYZ[#1] * CompZZ[#1]];
+  // zz
+  C_Mat_lambda~{3}~{3}[] = lambda * Tensor[ CompZX[$1#1] * CompZX[#1],   CompZX[#1] * CompZY[#1],   CompZX[#1] * CompZZ[#1],
+                                            CompZY[#1] * CompZX[#1],   CompZY[#1] * CompZY[#1],   CompZY[#1] * CompZZ[#1],
+                                            CompZZ[#1] * CompZX[#1],   CompZZ[#1] * CompZY[#1],   CompZZ[#1] * CompZZ[#1]];
+  
+  // 5.2. Total material stiffness : C_Mat_ij = C_Mat_lambda_ij +  C_Mat_mu_ij
+  //==========================================================================
+  For i In {1:3}
+    For j In {1:3}
+    C_Mat~{i}~{j}[] =  C_Mat_lambda~{i}~{j}[$1#1]  + Transpose[C_Mat_mu~{i}~{j}[#1]] ;
+    EndFor
+  EndFor
+
+}
+
+Function {
+  // 5.3 Total stiffness : C_Mat_ij +  C_Geo_ij
+  //===========================================
+  For i In {1:3} 
+    For j In {1:3}
+      C_Tot~{i}~{j}[] =  C_Mat~{i}~{j}[$1] + C_Geo~{i}~{j}[$1];      
+    EndFor
+  EndFor
+}
+
+Constraint {
+  { Name Displacement_1 ; // x
+    Case {
+      { Region Left; Value 0.0 ; }
+    }
+  }
+  { Name Displacement_2 ; // y
+    Case {
+      { Region Left; Value 0.0 ; }
+    }
+  }
+  { Name Displacement_3 ; // z
+    Case {
+      { Region Left; Value 0.0 ; }
+    }
+  }
+}
+
+
+Include "Jacobian_Lib.pro"
+Include "Integration_Lib.pro"
+Include "Total_Lagrangian.pro"
+
+Resolution {
+  { Name Mec_SMP ;
+    System {
+      { Name Sys_Mec ; NameOfFormulation Total_Lagrangian ; }
+    }
+    Operation {
+      SetGlobalSolverOptions["-petsc_prealloc 400"];
+      CreateDir["res/"];
+      Evaluate[ $syscount = 0 ];
+      InitSolution[Sys_Mec];
+      Evaluate[ $var_DTime = period/(num_steps)];
+      Evaluate[ $var_Dt_Max = period/(num_steps)];
+        
+      TimeLoopTheta[t_0, 1.0 * t_end, dtime[], theta_value]{ //TimeLoopNewmark[t_0, t_end, dtime[], beta, gamma] {
+
+        Generate[Sys_Mec]; Solve[Sys_Mec]; Evaluate[ $syscount = $syscount + 1 ];
+        Generate[Sys_Mec]; GetResidual[Sys_Mec, $res0]; Evaluate[ $res = $res0, $iter = 0 ];
+        Print[{$iter, $res, $res / $res0}, Format "Residual %03g: abs %14.12e rel %14.12e"];
+        While[$res > tol_abs[] && $res / $res0 > tol_rel[] && $res / $res0 <= 1 && $iter < iter_max]{
+          Solve[Sys_Mec]; Evaluate[ $syscount = $syscount + 1 ];
+          Generate[Sys_Mec]; GetResidual[Sys_Mec, $res]; Evaluate[ $iter = $iter + 1 ];
+          Print[{$iter, $res, $res / $res0}, Format "Residual %03g: abs %14.12e rel %14.12e"];
+        }
+        Test[ ($iter < iter_max) && ($res / $res0 <= 1) ]{
+          SaveSolution[Sys_Mec];
+          Test[ GetNumberRunTime[visu]{"Input/Solver/Visu"} ]{
+            PostOperation[Mec];
+          }
+          Test[ ($iter < iter_max / 10) && ($DTime < dt_max[]) ]{
+            Evaluate[ $dt_new = Min[$DTime * 2.0, dt_max[]] ];
+            Print[{$dt_new}, Format "*** Fast convergence: increasing time step to %g"];
+            SetDTime[$dt_new];
+          }
+        }
+        {
+          Evaluate[ $dt_new = $DTime/2.0 ];
+          Print[{$iter, $dt_new}, Format "*** Non convergence (iter %g): recomputing with reduced step %g"];
+          SetTime[$Time - $DTime];
+          SetTimeStep[$TimeStep - 1];
+          RemoveLastSolution[Sys_Mec];
+          SetDTime[$dt_new];
+          Evaluate[ $var_DTime = $dt_new ];
+        }
+      }
+      Print[{$syscount}, Format "Total number of linear systems solved: %g"];
+    }
+  }
+}
+
+PostOperation {
+  { Name Mec ; NameOfPostProcessing Total_Lagrangian; 
+    Operation {
+      Print[ u, OnElementsOf Domain_Mecha, File "res/u.pos"] ;
+
+      Print[ PK2_xx, OnElementsOf Domain_Mecha, File "res/PK2_xx.pos"] ;
+      Print[ PK2_xy, OnElementsOf Domain_Mecha, File "res/PK2_xy.pos"] ;
+      Print[ PK2_xz, OnElementsOf Domain_Mecha, File "res/PK2_xz.pos"] ;
+      Print[ PK2_yx, OnElementsOf Domain_Mecha, File "res/PK2_yx.pos"] ;
+      Print[ PK2_yy, OnElementsOf Domain_Mecha, File "res/PK2_yy.pos"] ;
+      Print[ PK2_yz, OnElementsOf Domain_Mecha, File "res/PK2_yz.pos"] ;
+      Print[ PK2_zx, OnElementsOf Domain_Mecha, File "res/PK2_zx.pos"] ;
+      Print[ PK2_zy, OnElementsOf Domain_Mecha, File "res/PK2_zy.pos"] ;
+      Print[ PK2_zz, OnElementsOf Domain_Mecha, File "res/PK2_zz.pos"] ;
+
+      Print[ F_tot_xx, OnElementsOf Domain_Mecha, File "res/F_tot_xx.pos"] ;
+      Print[ F_tot_xy, OnElementsOf Domain_Mecha, File "res/F_tot_xy.pos"] ;
+      Print[ F_tot_xz, OnElementsOf Domain_Mecha, File "res/F_tot_xz.pos"] ;
+      Print[ F_tot_yx, OnElementsOf Domain_Mecha, File "res/F_tot_yx.pos"] ;
+      Print[ F_tot_yy, OnElementsOf Domain_Mecha, File "res/F_tot_yy.pos"] ;
+      Print[ F_tot_yz, OnElementsOf Domain_Mecha, File "res/F_tot_yz.pos"] ;
+      Print[ F_tot_zx, OnElementsOf Domain_Mecha, File "res/F_tot_zx.pos"] ;
+      Print[ F_tot_zy, OnElementsOf Domain_Mecha, File "res/F_tot_zy.pos"] ;
+      Print[ F_tot_zz, OnElementsOf Domain_Mecha, File "res/F_tot_zz.pos"] ;
+    }
+  }
+}
+
+DefineConstant[
+  R_ = {"Mec_SMP", Name "GetDP/1ResolutionChoices", Visible 0},
+  C_ = {"-solve -bin -v 3 -v2", Name "GetDP/9ComputeCommand", Visible 0},
+  P_ = { "", Name "GetDP/2PostOperationChoices", Visible 0}
+];