From a07c3eb660ba46c0552b6fdb3b48e3578bb71864 Mon Sep 17 00:00:00 2001
From: Inno <niyinno@gmail.com>
Date: Fri, 29 Jun 2018 17:14:18 +0200
Subject: [PATCH] first version
---
HyperElasticity/Integration_Lib.pro | 37 +++
HyperElasticity/Jacobian_Lib.pro | 59 +++++
HyperElasticity/Total_Lagrangian.pro | 104 ++++++++
HyperElasticity/beam.pro | 380 +++++++++++++++++++++++++++
4 files changed, 580 insertions(+)
create mode 100644 HyperElasticity/Integration_Lib.pro
create mode 100644 HyperElasticity/Jacobian_Lib.pro
create mode 100644 HyperElasticity/Total_Lagrangian.pro
create mode 100644 HyperElasticity/beam.pro
diff --git a/HyperElasticity/Integration_Lib.pro b/HyperElasticity/Integration_Lib.pro
new file mode 100644
index 0000000..99599cf
--- /dev/null
+++ b/HyperElasticity/Integration_Lib.pro
@@ -0,0 +1,37 @@
+/*
+ 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 ; } }
+ }
+ }
+ }
+}
+
+/* --------------------------------------------------------------------------*/
diff --git a/HyperElasticity/Jacobian_Lib.pro b/HyperElasticity/Jacobian_Lib.pro
new file mode 100644
index 0000000..71f7b08
--- /dev/null
+++ b/HyperElasticity/Jacobian_Lib.pro
@@ -0,0 +1,59 @@
+/*
+ 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 ; }
+ }
+ }
+}
+
+
+
+/* --------------------------------------------------------------------------*/
diff --git a/HyperElasticity/Total_Lagrangian.pro b/HyperElasticity/Total_Lagrangian.pro
new file mode 100644
index 0000000..5f42170
--- /dev/null
+++ b/HyperElasticity/Total_Lagrangian.pro
@@ -0,0 +1,104 @@
+//=====================================================================
+
+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; } } }
+ }
+ }
+}
diff --git a/HyperElasticity/beam.pro b/HyperElasticity/beam.pro
new file mode 100644
index 0000000..00fab21
--- /dev/null
+++ b/HyperElasticity/beam.pro
@@ -0,0 +1,380 @@
+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}
+];
--
GitLab