simplify

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
+  Sur_Clamp = Region[{1}];
+  Vol_Mec = 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.8e8;
-    force_ext[] = Vector[0., force_amplitude, 0.0];
-}
-
-//Include "smp_stent_KMat_Elasticity.pro";
-
-
-
-
-
-
-
-
-
-
-
-
-
+  E  = 771.0e6; // Young modulus
+  nu = 0.29; // Poisson ratio
+  lambda = E/(2.0 * (1.0 + nu)); // first Lame parameter
+  mu = (nu * E)/((1.0 + nu) * (1.0 - 2 * nu)); // second Lame parameter
 
+  force_ext[] = Vector[0., 1.8e8, 0.0]; // volume force
 
+  tol_abs[] = 1e-8; // absolute tolerance for Newton-Raphson solver
+  tol_rel[] = 1e-8; // relative tolerance for Newton-Rasphon solver
+  iter_max = 50; // maximum number of iterations for Newton-Raphon solver
+}
 
 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
-  //=================================
+  // The nonlinear right hand side P_i: (i = 1 for x, 2 for y, 3 for z)
+
   //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
-  //==========================================================
+  // The stiffness matrix: C_Tot_ij = C_Mat_ij +  C_Geo_ij
+
+  // 1. Geometric stiffness: C_Geo_ij
+  //    (see (4.69) in "P. Wriggers, Nonlinear finite element methods, 2008")
 
-  // 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
-  //==========================================================================
+  // 2 Material stiffness: C_Mat_ij
+
+  // 2.1. Material stiffness - the part resulting from "mu" : C_Mat_mu_ij
+
+  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] ];
+  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] ];
+  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] ];
+  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] ];
+  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] ];
+  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] ];
+  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] ];
+  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] ];
+  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] ];
+
+  // 2.2. Material stiffness - the part resulting from "lambda" : C_Mat_lambda_ij
+
+  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]];
+  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]];
+  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]];
+  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]];
+  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]];
+  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]];
+  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]];
+  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]];
+  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]];
+
+  // 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]] ;
+      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
-  //===========================================
+  // 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];
@@ -276,25 +201,159 @@ Function {
 Constraint {
   { Name Displacement_1 ; // x
     Case {
-      { Region Left; Value 0.0 ; }
+      { Region Sur_Clamp; Value 0.0 ; }
     }
   }
   { Name Displacement_2 ; // y
     Case {
-      { Region Left; Value 0.0 ; }
+      { Region Sur_Clamp; Value 0.0 ; }
     }
   }
   { Name Displacement_3 ; // z
     Case {
-      { Region Left; Value 0.0 ; }
+      { Region Sur_Clamp; Value 0.0 ; }
+    }
+  }
+}
+
+Jacobian {
+  { Name JVol ;
+    Case {
+      { Region All ; Jacobian Vol ; }
+    }
+  }
+  { Name JSur ;
+    Case {
+      { Region All ; Jacobian Sur ; }
+    }
+  }
+  { Name JLin ;
+    Case {
+      { Region All ; Jacobian Lin ; }
+    }
+  }
+}
+
+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 ; } }
+      }
+    }
+  }
+}
+
+FunctionSpace {
+  For i In {1:3}
+    { Name H_u~{i} ; Type Form0 ;
+      BasisFunction {
+        { Name sn~{i} ; NameOfCoef un~{i} ; Function BF_Node ;
+          Support Region[ { Vol_Mec} ] ; 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 {
+      For i In {1:3}
+        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 Vol_Mec ; Jacobian JVol ; 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 Vol_Mec ; Jacobian JVol ; Integration GradGrad ; }
+        EndFor
+        Galerkin { [ P~{i}[ TensorDiag[1., 1., 1.] +
+                            TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ] , {d u~{i}} ] ;
+          In Vol_Mec; Jacobian JVol; Integration GradGrad; }
+      EndFor
+
+      Galerkin { [ - CompX[ -force_ext[] ] , {u~{1}} ] ;
+        In Vol_Mec; Jacobian JVol; Integration GradGrad; }
+      Galerkin { [ - CompY[ -force_ext[] ] , {u~{2}} ] ;
+        In Vol_Mec; Jacobian JVol; Integration GradGrad; }
+      Galerkin { [ - CompZ[ -force_ext[] ] , {u~{3}} ] ;
+        In Vol_Mec; Jacobian JVol; 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 Vol_Mec ; Jacobian JVol ; } } }
+      { Name um ; Value{ Term{ [ Norm[ Vector[ {u~{1}}, {u~{2}}, {u~{3}} ] ] ] ;
+            In Vol_Mec  ; Jacobian JVol; } } }
+      { Name ux  ; Value{ Term{ [ {u~{1}} ] ;
+            In Vol_Mec  ; Jacobian JVol; } } }
+      { Name uy  ; Value{ Term{ [ {u~{2}} ] ;
+            In Vol_Mec  ; Jacobian JVol; } } }
+      { Name uz  ; Value{ Term{ [ {u~{3}} ] ;
+            In Vol_Mec  ; Jacobian JVol; } } }
+
+      { 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 Vol_Mec  ; Jacobian JVol; } } }
+      { 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 Vol_Mec  ; Jacobian JVol; } } }
+      { 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 Vol_Mec  ; Jacobian JVol; } } }
+      { 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 Vol_Mec  ; Jacobian JVol; } } }
+      { 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 Vol_Mec  ; Jacobian JVol; } } }
+      { 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 Vol_Mec  ; Jacobian JVol; } } }
+      { 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 Vol_Mec  ; Jacobian JVol; } } }
+      { 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 Vol_Mec  ; Jacobian JVol; } } }
+      { 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 Vol_Mec  ; Jacobian JVol; } } }
+
+      { 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 Vol_Mec  ; Jacobian JVol; } } }
+      { 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 Vol_Mec  ; Jacobian JVol; } } }
+      { 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 Vol_Mec  ; Jacobian JVol; } } }
+      { 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 Vol_Mec  ; Jacobian JVol; } } }
+      { 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 Vol_Mec  ; Jacobian JVol; } } }
+      { 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 Vol_Mec  ; Jacobian JVol; } } }
+      { 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 Vol_Mec  ; Jacobian JVol; } } }
+      { 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 Vol_Mec  ; Jacobian JVol; } } }
+      { 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 Vol_Mec  ; Jacobian JVol; } } }
+    }
+  }
+}
 
-Include "Jacobian_Lib.pro"
-Include "Integration_Lib.pro"
-Include "Total_Lagrangian.pro"
 
 Resolution {
   { Name Mec_SMP ;
@@ -304,43 +363,17 @@ Resolution {
     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 ];
+      Generate[Sys_Mec]; Solve[Sys_Mec];
+      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"];
-        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"];
+      SaveSolution[Sys_Mec];
+      PostOperation[Mec];
     }
   }
 }
@@ -348,33 +381,33 @@ Resolution {
 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"] ;
+      Print[ u, OnElementsOf Vol_Mec, File "res/u.pos"] ;
+
+      Print[ PK2_xx, OnElementsOf Vol_Mec, File "res/PK2_xx.pos"] ;
+      Print[ PK2_xy, OnElementsOf Vol_Mec, File "res/PK2_xy.pos"] ;
+      Print[ PK2_xz, OnElementsOf Vol_Mec, File "res/PK2_xz.pos"] ;
+      Print[ PK2_yx, OnElementsOf Vol_Mec, File "res/PK2_yx.pos"] ;
+      Print[ PK2_yy, OnElementsOf Vol_Mec, File "res/PK2_yy.pos"] ;
+      Print[ PK2_yz, OnElementsOf Vol_Mec, File "res/PK2_yz.pos"] ;
+      Print[ PK2_zx, OnElementsOf Vol_Mec, File "res/PK2_zx.pos"] ;
+      Print[ PK2_zy, OnElementsOf Vol_Mec, File "res/PK2_zy.pos"] ;
+      Print[ PK2_zz, OnElementsOf Vol_Mec, File "res/PK2_zz.pos"] ;
+
+      Print[ F_tot_xx, OnElementsOf Vol_Mec, File "res/F_tot_xx.pos"] ;
+      Print[ F_tot_xy, OnElementsOf Vol_Mec, File "res/F_tot_xy.pos"] ;
+      Print[ F_tot_xz, OnElementsOf Vol_Mec, File "res/F_tot_xz.pos"] ;
+      Print[ F_tot_yx, OnElementsOf Vol_Mec, File "res/F_tot_yx.pos"] ;
+      Print[ F_tot_yy, OnElementsOf Vol_Mec, File "res/F_tot_yy.pos"] ;
+      Print[ F_tot_yz, OnElementsOf Vol_Mec, File "res/F_tot_yz.pos"] ;
+      Print[ F_tot_zx, OnElementsOf Vol_Mec, File "res/F_tot_zx.pos"] ;
+      Print[ F_tot_zy, OnElementsOf Vol_Mec, File "res/F_tot_zy.pos"] ;
+      Print[ F_tot_zz, OnElementsOf Vol_Mec, File "res/F_tot_zz.pos"] ;
     }
   }
 }
 
 DefineConstant[
   R_ = {"Mec_SMP", Name "GetDP/1ResolutionChoices", Visible 0},
-  C_ = {"-solve -bin -v 3 -v2 -pos", Name "GetDP/9ComputeCommand", Visible 0},
+  C_ = {"-solve -bin -v 3 -v2", Name "GetDP/9ComputeCommand", Visible 0},
   P_ = { "", Name "GetDP/2PostOperationChoices", Visible 0}
 ];