From de5d952ec367601b0f482d0fcc6c732f14234990 Mon Sep 17 00:00:00 2001 From: Christophe Geuzaine <cgeuzaine@ulg.ac.be> Date: Sat, 30 Jun 2018 09:33:29 +0200 Subject: [PATCH] simplify --- HyperElasticity/beam.pro | 615 +++++++++++++++++++++------------------ 1 file changed, 324 insertions(+), 291 deletions(-) diff --git a/HyperElasticity/beam.pro b/HyperElasticity/beam.pro index 2afda84..30568cc 100644 --- a/HyperElasticity/beam.pro +++ b/HyperElasticity/beam.pro @@ -1,271 +1,196 @@ -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} ]; -- GitLab