up

Elasticity/wrench2D.pro

@@ -13,40 +13,33 @@
 
 /* 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,
+Deflection = DefineNumber[4*AppliedForce*((LLength-0.018)/Width)^3/(Young*Thickness)*1e3,
                           Name "Solution/Deflection (analytical) [mm]", ReadOnly 1];
 
 Group {
@@ -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 {
@@ -234,10 +219,10 @@ Jacobian {
 
 /* Adapt the number of Gauss points to the polynomial degree of the finite elements
    is as simple as this: */
-If (FE_Degree == 1)
 Integration {
   { Name Gauss_v;
     Case {
+      If (FE_Degree == 1)
       { Type Gauss;
         Case {
           { GeoElement Line       ; NumberOfPoints  3; }
@@ -245,13 +230,7 @@ Integration {
           { GeoElement Quadrangle ; NumberOfPoints  4; }
         }
       }
-    }
-  }
-}
-ElseIf (FE_Degree == 2)
-Integration {
-  { Name Gauss_v;
-    Case {
+      Else
       { Type Gauss;
         Case {
 	  { GeoElement Line       ; NumberOfPoints  5; }
@@ -259,10 +238,10 @@ Integration {
           { 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},