From 30d765d1d4489a79796c7ebdf53aed0659277c08 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Thu, 15 Mar 2018 21:55:27 +0100
Subject: [PATCH] up
---
Elasticity/wrench2D.pro | 191 ++++++++++++++++++----------------------
1 file changed, 88 insertions(+), 103 deletions(-)
diff --git a/Elasticity/wrench2D.pro b/Elasticity/wrench2D.pro
index ce97783..5751251 100644
--- a/Elasticity/wrench2D.pro
+++ b/Elasticity/wrench2D.pro
@@ -13,41 +13,34 @@
/* Linear elasticity with GetDP:
- GetDP has a peculiar way to deal with linear elasticity.
- Instead of vector field "u = Vector[ ux, uy, uz ]",
- the displacement field is regarded as two (2D case) or 3 (3D case) scalar fields.
- Unlike conventional formulations, GetDP's formulation is written
- in terms of the gradient "Grad u" of the displacement field,
- which is a non-symmetric tensor, and the needed symmetrization
- (to define the strain tensor and relate it to the stress tensor)
- is done through the constitutive relationship (Hooke law).
- The reason for this unusual formulation is to be able to use also for elastic problems
- the powerful geometrical and homological kernel of GetDP,
- which relies on the operators Grad, Curl and Div.
-
- The "Grad u" formulation entails a small increase of assembly work
- but makes in counterpart lots of geometrical features
- implemented in GetDP (change of coordinates, function spaces, etc...)
- applicable to elastic problems out-of-the-box,
- since the scalar fields { ux, uy, uz } have exactly the same geometrical properties
- as, e.g. a scalar electrip potential or a temperature field.
-*/
-
+ GetDP has a peculiar way to deal with linear elasticity. Instead of a vector
+ field "u = Vector[ ux, uy, uz ]", the displacement field is regarded as two
+ (2D case) or 3 (3D case) scalar fields. Unlike conventional formulations
+ then, GetDP's formulation is written in terms of the gradient "Grad u" of the
+ displacement field, which is a non-symmetric tensor, and the needed
+ symmetrization (to define the strain tensor and relate it to the stress
+ tensor) is done through the constitutive relationship (Hooke law). The
+ reason for this unusual formulation is to be able to use also for elastic
+ problems the powerful geometrical and homological kernel of GetDP, which
+ relies on the operators Grad, Curl and Div.
+
+ The "Grad u" formulation entails a small increase of assembly work but makes
+ in counterpart lots of geometrical features implemented in GetDP (change of
+ coordinates, function spaces, etc...) applicable to elastic problems
+ out-of-the-box, since the scalar fields { ux, uy, uz } have exactly the same
+ geometrical properties as, e.g. a scalar electric potential or a temperature
+ field. */
Include "wrench2D_common.pro";
-Young = 1e9 *
- DefineNumber[ 200, Name "Material/Young modulus [GPa]"];
-Poisson =
- DefineNumber[ 0.3, Name "Material/Poisson coefficient []"];
-AppliedForce =
- DefineNumber[ 100, Name "Material/Applied force [N]"];
+Young = 1e9 * DefineNumber[ 200, Name "Material/Young modulus [GPa]"];
+Poisson = DefineNumber[ 0.3, Name "Material/Poisson coefficient []"];
+AppliedForce = DefineNumber[ 100, Name "Material/Applied force [N]"];
// Approximation of the maximum deflection by an analytical model:
// Deflection = PL^3/(3EI) with I = Width^3*Thickness/12
-Deflection =
- DefineNumber[4*AppliedForce*((LLength-0.018)/Width)^3/(Young*Thickness)*1e3,
- Name "Solution/Deflection (analytical) [mm]", ReadOnly 1];
+Deflection = DefineNumber[4*AppliedForce*((LLength-0.018)/Width)^3/(Young*Thickness)*1e3,
+ Name "Solution/Deflection (analytical) [mm]", ReadOnly 1];
Group {
// Physical regions: Give explicit labels to the regions defined in the .geo file
@@ -56,17 +49,16 @@ Group {
Force = Region[ 3 ];
// Abstract regions:
- Vol_Elast_Mec = Region[ { Wrench } ];
- Vol_Force_Mec = Region[ { Wrench } ];
- Sur_Clamp_Mec = Region[ { Grip } ];
- Sur_Force_Mec = Region[ { Force } ];
-/*
- Signification of the abstract regions:
- Vol_Elast_Mec Elastic domain
- Vol_Force_Mec Region with imposed volumic force
- Sur_Force_Mec Surface with imposed surface traction
- Sur_Clamp_Mec Surface with imposed zero displacements (all components)
-*/
+ Vol_Mec = Region[ Wrench ];
+ Vol_Force_Mec = Region[ Wrench ];
+ Sur_Clamp_Mec = Region[ Grip ];
+ Sur_Force_Mec = Region[ Force ];
+
+ /* Signification of the abstract regions:
+ - Vol_Mec : Elastic domain
+ - Vol_Force_Mec : Region with imposed volumic force
+ - Sur_Force_Mec : Surface with imposed surface traction
+ - Sur_Clamp_Mec : Surface with imposed zero displacements (all components) */
}
Function {
@@ -84,18 +76,17 @@ Function {
pressure_y[] = -AppliedForce/(SurfaceArea[]*Thickness); // downward vertical force
}
-
-/* Hooke law
+/* Hooke's law
The material law
sigma_ij = C_ijkl epsilon_ij
- is represented in 2D by 4 2x2 tensors C_ij[], i,j=1,2
- depending on the Lamé coefficients of the isotropic linear material,
+ is represented in 2D by 4 2x2 tensors C_ij[], i,j=1,2, depending on the Lamé
+ coefficients of the isotropic linear material,
- lambda = E[]*nu[]/(1.+nu[])/(1.-2.*nu[]);
- mu = E[]/2./(1.+nu[]);
+ lambda = E[]*nu[]/(1.+nu[])/(1.-2.*nu[])
+ mu = E[]/2./(1.+nu[])
as follows
@@ -103,8 +94,7 @@ Function {
EPD: a[] = lambda + 2 mu b[] = mu c[] = lambda
3D: a[] = lambda + 2 mu b[] = mu c[] = lambda
- respectively for the 2D plane strain (EPD), 2D plane stress (EPS) and 3D cases.
-*/
+ respectively for the 2D plane strain (EPD), 2D plane stress (EPS) and 3D cases. */
Function {
Flag_EPC = 1;
@@ -138,45 +128,40 @@ Constraint {
}
}
-/* As explained above, the displacement field is discretized
- as two scalar fields "ux" and "uy", which are the spatial components
- of the vector field "u" in a fixed Cartesian coordinate system.
+/* As explained above, the displacement field is discretized as two scalar
+ fields "ux" and "uy", which are the spatial components of the vector field
+ "u" in a fixed Cartesian coordinate system.
Boundary conditions like
ux = ... ;
uy = ... ;
- translate naturally into Dirichlet constraints
- on the scalar "ux" and "uy" FunctionSpaces.
- Conditions like, however,
+ translate naturally into Dirichlet constraints on the scalar "ux" and "uy"
+ FunctionSpaces. More exotic conditions like
u . n = ux Cos [th] + uy Sin [th] = ... ;
- are less naturally accounted for within the "Grad u" formulation.
- Fortunately, they are rather uncommon.
+ are less naturally accounted for within the "Grad u" formulation; but they
+ could be easily implemented with a Lagrange multplier.
- Finite element shape (triangles or quadrangles) makes no difference
- in the definition of the FunctionSpaces. The appropriate shape functions to be used
+ Finite element shape (triangles or quadrangles) makes no difference in the
+ definition of the FunctionSpaces. The appropriate shape functions to be used
are determined by GetDP at a much lower level on basis of the information
contained in the *.msh file.
- Second order elements, on the other hand, are implemented in the hierarchical fashion
- by adding to the first order node-based shape functions
- a set of second order edge-based functions
- to complete a basis for 2d order polynomials on the reference element.
- */
-
+ Second order elements, on the other hand, are implemented in the hierarchical
+ fashion by adding to the first order node-based shape functions a set of
+ second order edge-based functions to complete a basis for 2d order
+ polynomials on the reference element. */
// Domain of definition of the "ux" and "uy" FunctionSpaces
Group {
- Dom_H_u_Mec = Region[ { Vol_Elast_Mec,
- Sur_Force_Mec,
- Sur_Clamp_Mec} ];
+ Dom_H_u_Mec = Region[ { Vol_Mec, Sur_Force_Mec, Sur_Clamp_Mec} ];
}
-Flag_Degree =
- DefineNumber[ 0, Name "Geometry/Use degree 2 (hierarch.)", Choices{0,1}, Visible 1];
+Flag_Degree = DefineNumber[ 0, Name "Geometry/Use degree 2 (hierarch.)",
+ Choices{0,1}, Visible 1];
FE_Degree = ( Flag_Degree == 0 ) ? 1 : 2; // Convert flag value into polynomial degree
FunctionSpace {
@@ -233,36 +218,30 @@ Jacobian {
}
/* Adapt the number of Gauss points to the polynomial degree of the finite elements
-is as simple as this: */
-If (FE_Degree == 1)
+ is as simple as this: */
Integration {
{ Name Gauss_v;
Case {
- {Type Gauss;
+ If (FE_Degree == 1)
+ { Type Gauss;
Case {
- { GeoElement Line ; NumberOfPoints 3; }
+ { GeoElement Line ; NumberOfPoints 3; }
{ GeoElement Triangle ; NumberOfPoints 3; }
{ GeoElement Quadrangle ; NumberOfPoints 4; }
}
}
- }
- }
-}
-ElseIf (FE_Degree == 2)
-Integration {
- { Name Gauss_v;
- Case {
- {Type Gauss;
+ Else
+ { Type Gauss;
Case {
{ GeoElement Line ; NumberOfPoints 5; }
{ GeoElement Triangle ; NumberOfPoints 7; }
{ GeoElement Quadrangle ; NumberOfPoints 7; }
}
}
+ EndIf
}
}
}
-EndIf
Formulation {
{ Name Elast_u ; Type FemEquation ;
@@ -272,13 +251,13 @@ Formulation {
}
Equation {
Integral { [ -C_xx[] * Dof{d ux}, {d ux} ] ;
- In Vol_Elast_Mec ; Jacobian Vol ; Integration Gauss_v ; }
+ In Vol_Mec ; Jacobian Vol ; Integration Gauss_v ; }
Integral { [ -C_xy[] * Dof{d uy}, {d ux} ] ;
- In Vol_Elast_Mec ; Jacobian Vol ; Integration Gauss_v ; }
+ In Vol_Mec ; Jacobian Vol ; Integration Gauss_v ; }
Integral { [ -C_yx[] * Dof{d ux}, {d uy} ] ;
- In Vol_Elast_Mec ; Jacobian Vol ; Integration Gauss_v ; }
+ In Vol_Mec ; Jacobian Vol ; Integration Gauss_v ; }
Integral { [ -C_yy[] * Dof{d uy}, {d uy} ] ;
- In Vol_Elast_Mec ; Jacobian Vol ; Integration Gauss_v ; }
+ In Vol_Mec ; Jacobian Vol ; Integration Gauss_v ; }
Integral { [ force_x[] , {ux} ];
In Vol_Force_Mec ; Jacobian Vol ; Integration Gauss_v ; }
@@ -310,29 +289,36 @@ Resolution {
PostProcessing {
{ Name Elast_u ; NameOfFormulation Elast_u ;
PostQuantity {
- { Name u ; Value { Term { [ Vector[ {ux}, {uy}, 0 ]];
- In Vol_Elast_Mec ; Jacobian Vol ; } } }
- { Name uy ; Value { Term { [ 1e3*{uy} ];
- In Vol_Elast_Mec ; Jacobian Vol ; } } }
- { Name sig_xx ; Value { Term {
- [ CompX[ C_xx[]*{d ux} + C_xy[]*{d uy} ] ];
- In Vol_Elast_Mec ; Jacobian Vol ; } } }
- { Name sig_xy ; Value { Term {
- [ CompY[ C_xx[]*{d ux} + C_xy[]*{d uy} ] ];
- In Vol_Elast_Mec ; Jacobian Vol ; } } }
- { Name sig_yy ; Value { Term {
- [ CompY [ C_yx[]*{d ux} + C_yy[]*{d uy} ] ];
- In Vol_Elast_Mec ; Jacobian Vol ; } } }
+ { Name u ; Value {
+ Term { [ Vector[ {ux}, {uy}, 0 ]]; In Vol_Mec ; Jacobian Vol ; }
+ }
+ }
+ { Name uy ; Value {
+ Term { [ 1e3*{uy} ]; In Vol_Mec ; Jacobian Vol ; }
+ }
+ }
+ { Name sig_xx ; Value {
+ Term { [ CompX[ C_xx[]*{d ux} + C_xy[]*{d uy} ] ];
+ In Vol_Mec ; Jacobian Vol ; }
+ }
+ }
+ { Name sig_xy ; Value {
+ Term { [ CompY[ C_xx[]*{d ux} + C_xy[]*{d uy} ] ];
+ In Vol_Mec ; Jacobian Vol ; }
+ }
+ }
+ { Name sig_yy ; Value {
+ Term { [ CompY [ C_yx[]*{d ux} + C_yy[]*{d uy} ] ];
+ In Vol_Mec ; Jacobian Vol ; }
+ }
+ }
}
}
}
-
-
PostOperation {
{ Name pos; NameOfPostProcessing Elast_u;
Operation {
-
If(FE_Degree == 1)
Print[ sig_xx, OnElementsOf Wrench, File "sigxx.pos" ];
Print[ u, OnElementsOf Wrench, File "u.pos" ];
@@ -356,7 +342,6 @@ PostOperation {
}
}
-
// Tell Gmsh which GetDP commands to execute when running the model.
DefineConstant[
R_ = {"Elast_u", Name "GetDP/1ResolutionChoices", Visible 0},
--
GitLab