[Getfem-commits] r4911 - in /trunk/getfem: interface/src/ interface/tests/matlab/ src/ src/getfem/

[Yves . Renard]

Author: renard
Date: Thu Mar 26 14:16:45 2015
New Revision: 4911

URL: http://svn.gna.org/viewcvs/getfem?rev=4911&view=rev
Log:
Mindlin-Reissner plate model with MITC element

Added:
    trunk/getfem/interface/tests/matlab/demo_Mindlin_Reissner_plate.m
Modified:
    trunk/getfem/interface/src/gf_model_set.cc
    trunk/getfem/interface/tests/matlab/Makefile.am
    trunk/getfem/src/getfem/getfem_linearized_plates.h
    trunk/getfem/src/getfem_generic_assembly.cc
    trunk/getfem/src/getfem_linearized_plates.cc

Modified: trunk/getfem/interface/src/gf_model_set.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/src/gf_model_set.cc?rev=4911&r1=4910&r2=4911&view=diff
==============================================================================
--- trunk/getfem/interface/src/gf_model_set.cc  (original)
+++ trunk/getfem/interface/src/gf_model_set.cc  Thu Mar 26 14:16:45 2015
@@ -1755,9 +1755,6 @@
        out.pop().from_integer(int(ind));
        );
     
-       
-
-
     /address@hidden ind = ('add Kirchhoff-Love plate brick', @tmim mim, @str 
varname, @str dataname_D, @str dataname_nu [, @int region])
       Add a bilaplacian brick on the variable
       `varname` and on the mesh region `region`.
@@ -1927,6 +1924,55 @@
        );
 
 
+    /address@hidden ind = ('add Reisner-Mindlin plate brick', @tmim mim, @tmim 
mim_reduced, @str varname_U, @str varname_theta , @str param_E, @str param_nu, 
@str param_epsilon, @str param_kappa [,@int variant [, @int region]])
+      Add a term corresponding to the classical Reissner-Mindlin plate
+      model for which `U` is the transverse displacement,
+      `Theta` the rotation of
+      fibers normal to the midplane, 'param_E' the Young Modulus,
+      `param_nu` the poisson ratio,
+      `param_epsilon` the plate thickness,
+      `param_kappa` the shear correction factor. Note that since this brick
+      uses the high level generic assembly language, the parameter can
+      be regular expression of this language.
+      There are three variants.
+      `variant = 0` corresponds to the an
+      unreduced formulation and in that case only the integration
+      method `mim` is used. Practically this variant is not usable since
+      it is subject to a strong locking phenomenon.
+      `variant = 1` corresponds to a reduced integration where `mim` is
+      used for the rotation term and `mim_reduced` for the transverse
+      shear term. `variant = 2` (default) corresponds to the projection onto
+      a rotated RT0 element of the transverse shear term. For the moment, this
+      is adapted to quadrilateral only (because it is not sufficient to
+      remove the locking phenomenon on triangle elements). Note also that if
+      you use high order elements, the projection on RT0 will reduce the order
+      of the approximation.
+      Returns the brick index in the model.
+      @*/
+     sub_command
+        ("add Mindlin Reissner plate brick", 7, 9, 0, 1,
+         getfemint_mesh_im *gfi_mim = in.pop().to_getfemint_mesh_im();
+         getfemint_mesh_im *gfi_mim_reduced = in.pop().to_getfemint_mesh_im();
+         std::string varname_U = in.pop().to_string();
+         std::string varname_theta = in.pop().to_string();
+         std::string param_E = in.pop().to_string();
+         std::string param_nu = in.pop().to_string();
+         std::string param_epsilon = in.pop().to_string();
+         std::string param_kapa = in.pop().to_string();
+        size_type variant = size_type(2);
+         if (in.remaining()) variant = in.pop().to_integer();
+         size_type region = size_type(-1);
+         if (in.remaining()) region = in.pop().to_integer();
+         size_type ind = add_Mindlin_Reissner_plate_brick
+         (md->model(), gfi_mim->mesh_im(), gfi_mim_reduced->mesh_im(),
+          varname_U, varname_theta, param_E, param_nu, param_epsilon,
+          param_kapa, variant, region);
+         workspace().set_dependance(md, gfi_mim);
+         out.pop().from_integer(int(ind));
+         );
+       
+
+
     /address@hidden ind = ('add mass brick', @tmim mim, @str varname[, @str 
dataname_rho[, @int region]])
       Add mass term to the model relatively to the variable `varname`.
       If specified, the data `dataname_rho` should contain the

Modified: trunk/getfem/interface/tests/matlab/Makefile.am
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/Makefile.am?rev=4911&r1=4910&r2=4911&view=diff
==============================================================================
--- trunk/getfem/interface/tests/matlab/Makefile.am     (original)
+++ trunk/getfem/interface/tests/matlab/Makefile.am     Thu Mar 26 14:16:45 2015
@@ -25,6 +25,7 @@
        check_plasticity.m \
        demo_elasticity.m \
        demo_bilaplacian.m \
+       demo_Mindlin_Reissner_plate.m \
        demo_laplacian.m \
        demo_periodic_laplacian.m \
        demo_laplacian1D.m \

Added: trunk/getfem/interface/tests/matlab/demo_Mindlin_Reissner_plate.m
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/demo_Mindlin_Reissner_plate.m?rev=4911&view=auto
==============================================================================
--- trunk/getfem/interface/tests/matlab/demo_Mindlin_Reissner_plate.m   (added)
+++ trunk/getfem/interface/tests/matlab/demo_Mindlin_Reissner_plate.m   Thu Mar 
26 14:16:45 2015
@@ -0,0 +1,114 @@
+% Copyright (C) 20015-2015 FABRE Mathieu, SECK Mamadou, DALLERIT Valentin, 
Yves Renard
+%
+% This file is a part of GETFEM++
+%
+% Getfem++  is  free software;  you  can  redistribute  it  and/or modify it
+% under  the  terms  of the  GNU  Lesser General Public License as published
+% by  the  Free Software Foundation;  either version 3 of the License,  or
+% (at your option) any later version along with the GCC Runtime Library
+% Exception either version 3.1 or (at your option) any later version.
+% This program  is  distributed  in  the  hope  that it will be useful,  but
+% WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+% or  FITNESS  FOR  A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+% License and GCC Runtime Library Exception for more details.
+% You  should  have received a copy of the GNU Lesser General Public License
+% along  with  this program;  if not, write to the Free Software Foundation,
+% Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.
+
+% Simple supported Mindlin-Reissner plate
+
+
+Emodulus = 1;          % Young Modulus
+nu       = 0.5;        % Poisson Coefficient
+epsilon  = 0.001;      % Plate thickness
+kappa     = 5/6;       % Shear correction factor
+f = -10*epsilon^3;     % Prescribed force on the top of the plate
+
+dirichlet_version = 1; % 0 = simplification, 1 = with multipliers, 2 = 
penalization
+quadrangles = true;    % Locking free only on quadrangle for the moment
+K = 1;                 % Degree of the finite element method
+variant = 1;           % 0 : not reduced, 1 : with reduced integration, 2 : 
MITC reduction
+with_Mindlin_brick = true; % Uses the Reissner-Mindlin predefined brick or not
+
+plot_mesh = false;
+draw_solution = true;
+
+% trace on;
+gf_workspace('clear all');
+NX = 20;
+if (quadrangles)
+  m = gf_mesh('cartesian',[0:1/NX:1],[0:1/NX:1]);
+else
+  
m=gf_mesh('import','structured',sprintf('GT="GT_PK(2,1)";SIZES=[1,1];NOISED=0;NSUBDIV=[%d,%d];',
 NX, NX));
+end
+
+% Create a mesh_fem of for a 2 dimension vector field
+mftheta = gf_mesh_fem(m,2);
+mfu = gf_mesh_fem(m,1);
+% Assign the QK or PK fem to all convexes of the mesh_fem, and define an
+% integration method
+if (quadrangles)
+  gf_mesh_fem_set(mftheta,'fem',gf_fem(sprintf('FEM_QK(2,%d)', K)));
+  gf_mesh_fem_set(mfu,'fem',gf_fem(sprintf('FEM_QK(2,%d)', K)));
+  mim = gf_mesh_im(m, gf_integ('IM_GAUSS_PARALLELEPIPED(2,6)'));
+  mim_reduced = gf_mesh_im(m, gf_integ('IM_GAUSS_PARALLELEPIPED(2,1)'));
+else
+  gf_mesh_fem_set(mftheta,'fem',gf_fem(sprintf('FEM_PK(2,%d)', K)));
+  gf_mesh_fem_set(mfu,'fem',gf_fem(sprintf('FEM_PK(2,%d)', K)));
+  mim = gf_mesh_im(m, gf_integ('IM_TRIANGLE(6)'));
+  mim_reduced = gf_mesh_im(m, gf_integ('IM_TRIANGLE(1)'));
+end
+
+% detect the border of the mesh
+border = gf_mesh_get(m,'outer faces');
+% mark it as boundary #1
+gf_mesh_set(m, 'boundary', 1, border);
+if (plot_mesh)
+  gf_plot_mesh(m, 'regions', [1]); % the boundary edges appears in red
+  pause(1);
+end
+
+md=gf_model('real');
+gf_model_set(md, 'add elementary rotated RT0 projection', 'RT0_projection');
+gf_model_set(md, 'add fem variable', 'u', mfu);
+gf_model_set(md, 'add fem variable', 'theta', mftheta);
+gf_model_set(md, 'add initialized data', 'E', Emodulus);
+gf_model_set(md, 'add initialized data', 'nu', nu);
+gf_model_set(md, 'add initialized data', 'epsilon', epsilon);
+gf_model_set(md, 'add initialized data', 'kappa', kappa);
+
+
+if (with_Mindlin_brick)
+  gf_model_set(md, 'add Mindlin Reissner plate brick', mim, mim_reduced, 'u', 
'theta', 'E', 'nu', 'epsilon', 'kappa', variant);
+else
+  gf_model_set(md, 'add linear generic assembly brick', mim, 
'(E*epsilon*epsilon*epsilon*(1-nu)/(48 * (1 - nu*nu))) * 
((Grad_theta+Grad_theta''):(Grad_Test_theta+Grad_Test_theta''))');
+  gf_model_set(md, 'add linear generic assembly brick', mim, 
'(E*epsilon*epsilon*epsilon*nu/(12 * (1 - nu*nu))) * 
(Trace(Grad_theta)*Trace(Grad_Test_theta))');
+  if (variant == 0)
+    gf_model_set(md, 'add linear generic assembly brick', mim, 
'(E*kappa*epsilon/(1 + nu)) * ((Grad_u + theta).Grad_Test_u) + 
(E*kappa*epsilon/(1 + nu)) * ((Grad_u + theta).Test_theta)');
+  elseif (variant == 1)
+    gf_model_set(md, 'add linear generic assembly brick', mim_reduced, 
'(E*kappa*epsilon/(1 + nu)) * ((Grad_u + theta).Grad_Test_u) + 
(E*kappa*epsilon/(1 + nu)) * ((Grad_u + theta).Test_theta)');
+  else
+    gf_model_set(md, 'add linear generic assembly brick', mim, 
'(E*kappa*epsilon/(1 + nu)) * ((Grad_u + 
Elementary_transformation(theta,RT0_projection)).Grad_Test_u) + 
(E*kappa*epsilon/(1 + nu)) * ((Grad_u + Elementary_transformation(theta, 
RT0_projection)).(Elementary_transformation(Test_theta, RT0_projection)))');  
+  end
+end
+
+gf_model_set(md, 'add initialized data', 'VolumicData', f);
+
+gf_model_set(md, 'add source term brick', mim, 'u', 'VolumicData');
+gf_model_set(md, 'add initialized data', 'DirichletData', 0);
+switch (dirichlet_version)
+  case 0,
+    gf_model_set(md, 'add Dirichlet condition with simplification', 'u', 1, 
'DirichletData');   
+  case 1, 
+    gf_model_set(md, 'add Dirichlet condition with multipliers', mim, 'u', 
mfu, 1, 'DirichletData');
+  case 2,
+    gf_model_set(md, 'add Dirichlet condition with penalization', mim, 'u', r, 
1, 'DirichletData');
+end
+gf_model_get(md, 'solve');
+U = gf_model_get(md, 'variable', 'u');
+
+if (draw)
+  gf_plot(mfu,U,'mesh','off', 'zplot', 'on'); 
+  colorbar; title('computed solution');
+end
+

Modified: trunk/getfem/src/getfem/getfem_linearized_plates.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_linearized_plates.h?rev=4911&r1=4910&r2=4911&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_linearized_plates.h  (original)
+++ trunk/getfem/src/getfem/getfem_linearized_plates.h  Thu Mar 26 14:16:45 2015
@@ -49,6 +49,45 @@
   */
   void add_2D_rotated_RT0_projection(model &md, std::string name);
 
+
+  /** Add a term corresponding to the classical Reissner-Mindlin plate
+      model for which `U` is the transverse displacement,
+      `Theta` the rotation of
+      fibers normal to the midplane, 'param_E' the Young Modulus,
+      `param_nu` the poisson ratio,
+      `param_epsilon` the plate thickness,
+      `param_kappa` the shear correction factor. Note that since this brick
+      uses the high level generic assembly language, the parameter can
+      be regular expression of this language.
+      There are three variants.
+      `variant = 0` corresponds to the an
+      unreduced formulation and in that case only the integration
+      method `mim` is used. Practically this variant is not usable since
+      it is subject to a strong locking phenomenon.
+      `variant = 1` corresponds to a reduced integration where `mim` is
+      used for the rotation term and `mim_reduced` for the transverse
+      shear term. `variant = 2` (default) corresponds to the projection onto
+      a rotated RT0 element of the transverse shear term. For the moment, this
+      is adapted to quadrilateral only (because it is not sufficient to
+      remove the locking phenomenon on triangle elements). Note also that if
+      you use high order elements, the projection on RT0 will reduce the order
+      of the approximation.
+      Returns the brick index in the model.
+   */
+  size_type add_Mindlin_Reissner_plate_brick
+  (model &md, const mesh_im &mim, const mesh_im &mim_reduced,
+   const std::string &U,
+   const std::string &Theta, const std::string &param_E,
+   const std::string &param_nu, const std::string &param_epsilon,
+   const std::string &param_kappa, size_type variant = size_type(2), 
+   size_type region = size_type(-1));
+
+
+
+
+
+
+
   /* ******************************************************************** */
   /*            Linear plate specific assembly procedures.                */
   /* ******************************************************************** */

Modified: trunk/getfem/src/getfem_generic_assembly.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_generic_assembly.cc?rev=4911&r1=4910&r2=4911&view=diff
==============================================================================
--- trunk/getfem/src/getfem_generic_assembly.cc (original)
+++ trunk/getfem/src/getfem_generic_assembly.cc Thu Mar 26 14:16:45 2015
@@ -232,14 +232,17 @@
 
   static void ga_throw_error_msg(const std::string &expr, size_type pos,
                                  const std::string &msg) {
+    int lenght_before = 40, lenght_after = 40;
     if (expr.size()) {
-      int first = std::max(0, int(pos)-40);
-      int last = std::min(int(pos)+20, int(expr.size()));
-      if (last - first < 60)
-      first = std::max(0, int(pos)-40-(60-last+first));
-      if (last - first < 60)
-        last = std::min(int(pos)+20+(60-last+first),int(expr.size()));
-
+      int first = std::max(0, int(pos)-lenght_before);
+      int last = std::min(int(pos)+lenght_after, int(expr.size()));
+      if (last - first < lenght_before+lenght_after)
+      first = std::max(0, int(pos)-lenght_before
+                       -(lenght_before+lenght_after-last+first));
+      if (last - first < lenght_before+lenght_after)
+        last = std::min(int(pos)+lenght_after
+                        +(lenght_before+lenght_after-last+first),
+                        int(expr.size()));
       if (first > 0) cerr << "...";
       cerr << expr.substr(first, last-first);
       if (last < int(expr.size())) cerr << "...";
@@ -6215,8 +6218,8 @@
             { test = 1; name = name.substr(5); }
 
           if (!(workspace.variable_exists(name)))
-            ga_throw_error(expr, pnode->pos,
-                           "Unknown variable, function, operator or data");
+            ga_throw_error(expr, pnode->pos, "Unknown variable, function, "
+                           "operator or data " + name);
 
           if (pnode->der1)
             ga_throw_error(expr, pnode->pos, "Derivative is for functions or "

Modified: trunk/getfem/src/getfem_linearized_plates.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_linearized_plates.cc?rev=4911&r1=4910&r2=4911&view=diff
==============================================================================
--- trunk/getfem/src/getfem_linearized_plates.cc        (original)
+++ trunk/getfem/src/getfem_linearized_plates.cc        Thu Mar 26 14:16:45 2015
@@ -1,7 +1,9 @@
 /* -*- c++ -*- (enables emacs c++ mode) */
 /*===========================================================================
  
- Copyright (C) 2004-2015 Yves Renard
+ Copyright (C) 2004-2015 Yves Renard, Jeremie Lasry, Mathieu Fabre
+                         SECK Mamadou, DALLERIT Valentin
+ 
  
  This file is a part of GETFEM++
  
@@ -26,6 +28,75 @@
 
 
 namespace getfem {
+
+  size_type add_Mindlin_Reissner_plate_brick(model &md,
+                                             const mesh_im & mim,
+                                             const mesh_im & mim_red,
+                                             const std::string &U,
+                                             const std::string &Theta,
+                                             const std::string &param_E,
+                                             const std::string &param_nu,
+                                             const std::string &param_epsilon,
+                                             const std::string &param_kappa,
+                                             size_type variant,
+                                             size_type region) {
+    
+    
+    std::string test_U = "Test_" + sup_previous_and_dot_to_varname(U);
+    std::string test_Theta = "Test_" + sup_previous_and_dot_to_varname(Theta);
+    std::string proj_Theta = (variant == 2) ?
+      ("Elementary_transformation("+Theta+",_2D_rotated_RT0_projection__434)")
+      : Theta;
+    std::string proj_test_Theta = (variant == 2) ?
+      ("Elementary_transformation("+test_Theta
+       +",_2D_rotated_RT0_projection__434)") : test_Theta;
+    
+    std::string D = "(("+param_E+")*pow("+param_epsilon+
+      ",3))/(12*(1-sqr("+param_nu+")))";
+    std::string G = "(("+param_E+")*("+param_epsilon+"))*("+param_kappa+
+      ")/(2*(1+("+param_nu+")))";
+    std::string E_theta = "(Grad_" + Theta + "+(Grad_" + Theta + ")')/2";
+    std::string E_test_Theta="(Grad_"+test_Theta+"+(Grad_"+test_Theta+")')/2";
+    
+    std::string expr_left =
+      D+"*(( 1-("+param_nu+"))*("+E_theta+"):("+E_test_Theta+")+("+param_nu+
+      ")*Trace("+E_theta+")*Trace("+E_test_Theta+"))";
+    
+    std::string expr_right = 
+      "("+G+")*(Grad_"+U+"-"+proj_Theta+").Grad_"+test_U+
+      "-("+G+")*(Grad_"+U+"-"+proj_Theta+")."+proj_test_Theta;
+    
+    switch(variant) {
+    case 0: // Without reduction
+      return add_linear_generic_assembly_brick
+        (md, mim, expr_left+"+"+expr_right, region, false, false,
+         "Reissner-Mindlin plate model brick");
+    case 1: // With reduced integration
+      add_linear_generic_assembly_brick
+        (md, mim, expr_left, region, false, false,
+         "Reissner-Mindlin plate model brick, rotation term");
+      return add_linear_generic_assembly_brick
+        (md, mim_red, expr_right, region, false, false,
+         "Reissner-Mindlin plate model brick, transverse shear term");
+    case 2: // Variant with projection on rotated RT0
+      add_2D_rotated_RT0_projection(md, "_2D_rotated_RT0_projection__434");
+      return add_linear_generic_assembly_brick
+        (md, mim, expr_left+"+"+expr_right, region, false, false,
+         "Reissner-Mindlin plate model brick");
+      break;
+    default: GMM_ASSERT1(false, "Invalid variant for Reissner-Mindlin brick.");
+    }
+    return size_type(-1);
+  }
+  
+
+
+
+
+
+
+
+
 
 
   // For the moment, only projection onto rotated RT0 element in dimension 2