[Getfem-commits] r4916 - in /trunk/getfem: doc/sphinx/source/ doc/sphinx/source/userdoc/ interface/src/ interface/tests/matlab/

[Yves . Renard]

Author: renard
Date: Fri Mar 27 14:26:07 2015
New Revision: 4916

URL: http://svn.gna.org/viewcvs/getfem?rev=4916&view=rev
Log:
Minor fixes and documentation for Mindlin plates

Added:
    trunk/getfem/doc/sphinx/source/biblio.bib
    trunk/getfem/doc/sphinx/source/userdoc/model_Mindlin_plate.rst
Modified:
    trunk/getfem/doc/sphinx/source/biblio.rst
    trunk/getfem/doc/sphinx/source/userdoc/gasm_high.rst
    trunk/getfem/doc/sphinx/source/userdoc/model.rst
    trunk/getfem/interface/src/gf_model_set.cc
    trunk/getfem/interface/tests/matlab/demo_Mindlin_Reissner_plate.m

Added: trunk/getfem/doc/sphinx/source/biblio.bib
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/doc/sphinx/source/biblio.bib?rev=4916&view=auto
==============================================================================
--- trunk/getfem/doc/sphinx/source/biblio.bib   (added)
+++ trunk/getfem/doc/sphinx/source/biblio.bib   Fri Mar 27 14:26:07 2015
@@ -0,0 +1,11 @@
+
+
+
address@hidden,
+  author =      " K.J. Bathe, E.N. Dvorkin",
+  title =       A four-node plate bending element based on Mindlin-Reissner 
plate theory and a mixed interpolation,
+  year =        "1985",
+  journal =     "Internat. J. Numer. Methods Engrg.",
+  volume =       "21",
+  pages =       "367-383",
+}

Modified: trunk/getfem/doc/sphinx/source/biblio.rst
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/doc/sphinx/source/biblio.rst?rev=4916&r1=4915&r2=4916&view=diff
==============================================================================
--- trunk/getfem/doc/sphinx/source/biblio.rst   (original)
+++ trunk/getfem/doc/sphinx/source/biblio.rst   Fri Mar 27 14:26:07 2015
@@ -13,6 +13,13 @@
 
 .. [bank1983] R.E. Bank, A.H. Sherman, A. Weiser,
    *Refinement algorithms and data structures for regular local mesh 
refinement*, in Scientific Computing IMACS, Amsterdam, North-Holland, pp 3-17, 
(1983).
+
+.. [ba-dv1985] K.J. Bathe, E.N. Dvorkin,
+   *A four-node plate bending element based on Mindlin-Reissner plate theory 
and a mixed interpolation*, Internat. J. Numer. Methods Engrg., 21 (1985), 
367-383.
+
+.. [br-ba-fo1989] F. Brezzi, K.J. Bathe, M. Fortin,
+   *Mixed-interpolated element for Reissner-Mindlin plates*, Internat. J. 
Numer. Methods Engrg., 28 (1989), 1787-1801.
+
 
 .. [ca-re-so1994] D. Calvetti, L. Reichel and D.C. Sorensen.
    *An implicitely restarted Lanczos method for large symmetric eigenvalue 
problems*. Electronic Transaction on Numerical Analysis}. 2:1-21, (1994).
@@ -31,10 +38,16 @@
    *The Finite Element Method Displayed*, J. Wiley & Sons, New York, (1984).
 
 .. [Dh-Go-Ku2003] A. Dhooge, W. Govaerts and Y. A. Kuznetsov,
-   *MATCONT: A MATLAB Package for Numerical Bifurcation Analysis of ODEs*, ACM 
Trans. Math. Software 31 (2003), 141-164.
+   *MATCONT: A MATLAB Package for Numerical Bifurcation Analysis of ODEs*,
+   ACM Trans. Math. Software 31 (2003), 141-164.
+
+.. [Duan2014] H. Duan,
+   *A finite element method for Reissner-Mindlin plates*,
+   Math. Comp., 83:286 (2014), 701-733.
 
 .. [Fa-Pa2003] F. Facchinei and J.-S. Pang,
-   *Finite-Dimensional Variational Inequalities and Complementarity Problems, 
Vol. II*, Springer Series in Operations Research, Springer, New York, 2003.
+   *Finite-Dimensional Variational Inequalities and Complementarity Problems, 
Vol. II*,
+   Springer Series in Operations Research, Springer, New York, 2003.
 
 .. [Georg2001] K. Georg,
    *Matrix-free numerical continuation and bifurcation*, Numer. Funct. Anal. 
Optimization 22 (2001), 303-320.
@@ -63,8 +76,11 @@
 .. [Li-Re2015hal] T. |ligursky| and Y. Renard,
    *A Method of Piecewise-Smooth Numerical Branching.* Submitted (2015). 
<`hal-01113564 <https://hal.archives-ouvertes.fr/hal-01113564>`_>
 
+.. [Mi-Zh2002] P. Ming and Z. Shi,
+   *Optimal L2 error bounds for MITC3 type element.* Numer. Math. 91 (2002), 
77-91.
+
 .. [Xfem] N. Moes, J. Dolbow and T. Belytschko,
-   *A finite element method for crack growth without remeshing*, Int. J. Num. 
Meth. Engng. 46 (1999), 131-150.
+   *A finite element method for crack growth without remeshing*, Internat. J. 
Numer. Methods Engrg., 46 (1999), 131-150.
 
 .. [nedelec1991] J.-C. Nedelec.
    *Notions sur les techniques d'elements finis*, Ellipses, SMAI, 
Mathematiques & Applications no 7, (1991).
@@ -76,7 +92,7 @@
    *Domain decomposition algorithms for the p-version finite element method 
for elliptic problems*, Luca F. Pavarino. PhD thesis, Courant Institute of 
Mathematical Sciences}. 1992.
 
 .. [remacle2002] J-F. Remacle, M. Shephard,
-   *An algorithm oriented database*,  Int. J. Num. Meth. Engng. 58 (2003), 
349-374.
+   *An algorithm oriented database*,  Internat. J. Numer. Methods Engrg., 58 
(2003), 349-374.
 
 .. [so-se-do2004] P. |solin|, K. Segeth, I. |dolezel| ,
    *Higher-Order Finite Element Methods*, Chapman and Hall/CRC, Studies in 
advanced mathematics, 2004.
@@ -93,3 +109,9 @@
 .. |ligursky| unicode:: Ligursk U+00FD
 .. |solin| unicode:: U+0160 ol U+00ED n 
    :rtrim:
+
+
+.. For recent sphinx version, see 
http://sphinxcontrib-bibtex.readthedocs.org/en/latest/usage.html
+   .. bibliography:: biblio.bib
+   :style: unsrt
+   :all:

Modified: trunk/getfem/doc/sphinx/source/userdoc/gasm_high.rst
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/doc/sphinx/source/userdoc/gasm_high.rst?rev=4916&r1=4915&r2=4916&view=diff
==============================================================================
--- trunk/getfem/doc/sphinx/source/userdoc/gasm_high.rst        (original)
+++ trunk/getfem/doc/sphinx/source/userdoc/gasm_high.rst        Fri Mar 27 
14:26:07 2015
@@ -663,6 +663,8 @@
 
 **CAUTION**: You have to think that when some variables are used in the 
transformation, the computation of the tangent system takes into account these 
dependence. However, the second derivative of a transformation with respect to 
a variable used has not been implemented. Thus, such a transformation is not 
allowed in the definition of a potential since it cannot be derived twice.
 
+.. _ud-gasm-high-elem-trans:
+
 Elementary transformations
 **************************
 

Modified: trunk/getfem/doc/sphinx/source/userdoc/model.rst
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/doc/sphinx/source/userdoc/model.rst?rev=4916&r1=4915&r2=4916&view=diff
==============================================================================
--- trunk/getfem/doc/sphinx/source/userdoc/model.rst    (original)
+++ trunk/getfem/doc/sphinx/source/userdoc/model.rst    Fri Mar 27 14:26:07 2015
@@ -53,6 +53,7 @@
    model_linear_elasticity
    model_mass
    model_bilaplacian
+   model_Mindlin_plate
    model_time_integration
    model_nonlinear_elasticity
    model_contact_friction

Added: trunk/getfem/doc/sphinx/source/userdoc/model_Mindlin_plate.rst
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/doc/sphinx/source/userdoc/model_Mindlin_plate.rst?rev=4916&view=auto
==============================================================================
--- trunk/getfem/doc/sphinx/source/userdoc/model_Mindlin_plate.rst      (added)
+++ trunk/getfem/doc/sphinx/source/userdoc/model_Mindlin_plate.rst      Fri Mar 
27 14:26:07 2015
@@ -0,0 +1,115 @@
+.. $Id: model.rst 3655 2010-07-17 20:42:08Z renard $
+
+.. include:: ../replaces.txt
+
+.. highlightlang:: c++
+
+.. index:: models, model bricks
+
+.. _ud-model-mindlin-plate:
+
+Mindlin-Reissner plate model
+----------------------------
+
+This brick implements the classical Mindlin-Reissner bending model for 
isotropic plates.
+
+The Mindlin-Reissner plate model
+++++++++++++++++++++++++++++++++
+
+Let :math:`\Omega \subset \R^2` be the reference configuration of the 
mid-plane of a plate of thickness :math:`\epsilon`.
+
+The weak formulation of the Mindlin-Reissner model for isotropic material can 
be written as follows for :math:`u_3` the transverse displacement and 
:math:`\theta` the rotation of fibers normal to the mid-plane:
+
+.. math::
+
+  & \int_{\Omega} D \epsilon^3\left((1-v)\gamma(\theta):\gamma(\psi) + \nu 
\mbox{div}(\theta)\mbox{div}(\psi)\right) dx \\
+  & ~~~~~~~~~~~~~~ + \int_{\Omega}G\epsilon (\nabla u_3 - \theta)\cdot(\nabla 
v_3 - \psi)dx = \int_{\Omega} F_3v_3 + M.\psi dx,
+
+for all admissible test functions :math:`v_3 : \Omega \rightarrow \R,$ $\psi : 
\Omega \rightarrow \R^2` and where:
+
+.. math::
+  & D = \Frac{E}{12(1-\nu^2)}, ~~ G = \Frac{E\kappa}{2(1+\nu)}, \\
+  & \gamma(\theta) = (\nabla \theta + \nabla \theta^T)/2, \\
+  & F_3 = \int_{-\epsilon/2}^{\epsilon/2} f_3dx_3 + g_3^+ + g_3^-, \\
+  & M_{\alpha} = \epsilon(g^+_{\alpha} - g^-_{\alpha})/2 +  
\int_{-\epsilon/2}^{\epsilon/2} x_3 f_{\alpha}dx_3, \alpha \in \{1, 2\},
+
+:math:`f` being a volumic force applied inside the three dimensional plate, 
:math:`g^+` and :math:`g^-` a force applied on the top and bottom surface of 
the plate, :math:`E` Young's modulus, :math:`\nu` Poisson's ratio and 
:math:`\kappa` the shear correction factor (usually set to 5/6).
+
+The classical boundary conditions are the following:
+
+  - Simple support :  a dirichlet condition on :math:`u_3`.
+
+  - Clamped support : a dirichlet condition on both :math:`u_3` and 
:math:`\theta`.
+
+  - Prescribed transverse force : boundary source term on :math:`u_3`.
+
+  - Prescribed moment : boundary source term on :math:`\theta`.
+
+
+An important issue of this model is that it is subjected to the so called 
shear locking so that a direct Galerkin procedure do not give a satisfactory 
approximation. There is several ways to circumvent the shear locking problem : 
reduced integration, projection of the transverse shear term, mixed methods. 
The two first method are proposed.
+
+
+Reduced integration of the transverse shear term
+================================================
+
+This strategy consists simply to use a lower order integration method to 
numerically compute the term
+
+.. math::
+
+  \int_{\Omega}G\epsilon (\nabla u_3 - \theta)\cdot(\nabla v_3 - \psi)dx
+
+This strategy is working properly at least when both the rotation and the 
transverse displacement is approximated with Q1 quadrilateral element with a 
degree one reduced integration method (the so-called QUAD4 element).
+
+Projection of the transverse shear term
+=======================================
+
+Another strategy comes from the MITC elements (Mixed Interpolation of 
Tensorial Components) which correspond to a re-interpretation in terms of 
projection of some mixed methods. The most popular element of this type is the 
MITC4 which correspond to the quadrilateral element Q1 with a projection of the 
transverse shear term on a rotated Raviart-Thomas element of lowest degree 
(RT0) (see [ba-dv1985]_, [br-ba-fo1989]_). This means that the transverse shear 
term becomes
+
+.. math::
+
+  \int_{\Omega}G\epsilon P^h(\nabla u_3 - \theta)\cdot P^h(\nabla v_3 - \psi)dx
+
+where :math:`P^h(T)` is the elementwize :math:`L^2`-projection onto the 
rotated RT0 space.  For the moment, the only projection implemented is the 
previous described one (projection on rotated RT0 space for quadrilateral 
element). Higher degree elements and triangular elements can be found in the 
litterature (see [Mi-Zh2002]_, [br-ba-fo1989]_, [Duan2014]_) and will be under 
consideration for a future implementation. Note also that since 
:math:`P^h(\nabla u_3) = \nabla u_3`, the term reduces to
+
+.. math::
+
+  \int_{\Omega}G\epsilon (\nabla u_3 - P^h(\theta))\cdot(\nabla v_3 - 
P^h(\psi))dx
+
+The principle of the definition of an elementary projection is explained if 
the description of the high-level generica assembly language (see 
:ref:`ud-gasm-high-elem-trans`) and an example can be found in the file 
:file:`src/getfem_linearized_plate.cc`.
+
+
+
+Add a Mindlin-Reissner plate model brick to a model
++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+The following function defined in :file:`getfem/getfem_linearized_plate.h` 
allows to add a Mindlin-Reissner plate model term to a transverse displacement 
`u3` and a rotation `theta`::
+
+  size_type add_Mindlin_Reissner_plate_brick
+  (model, mim, mim_reduced, name_u3, name_theta, param_E,
+   param_nu, param_epsilon, param_kappa, variant = 2, region)
+
+where `name_u3` is name of the variable which represents the transverse 
displacmenent, `name_theta` the variable which represents the rotation, 
'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.
+
+
+
+
+
+The projection on rotated RTO element can be added to a model thanks to the 
following function::
+
+  void add_2D_rotated_RT0_projection(model, transname);
+

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=4916&r1=4915&r2=4916&view=diff
==============================================================================
--- trunk/getfem/interface/src/gf_model_set.cc  (original)
+++ trunk/getfem/interface/src/gf_model_set.cc  Fri Mar 27 14:26:07 2015
@@ -1924,10 +1924,10 @@
        );
 
 
-    /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]])
+    /address@hidden ind = ('add Reisner-Mindlin plate brick', @tmim mim, @tmim 
mim_reduced, @str varname_u3, @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
+      model for which `varname_u3` is the transverse displacement,
+      `varname_theta` the rotation of
       fibers normal to the midplane, 'param_E' the Young Modulus,
       `param_nu` the poisson ratio,
       `param_epsilon` the plate thickness,

Modified: 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=4916&r1=4915&r2=4916&view=diff
==============================================================================
--- trunk/getfem/interface/tests/matlab/demo_Mindlin_Reissner_plate.m   
(original)
+++ trunk/getfem/interface/tests/matlab/demo_Mindlin_Reissner_plate.m   Fri Mar 
27 14:26:07 2015
@@ -22,20 +22,20 @@
 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
+f = -5*epsilon^3;      % Prescribed force on the top of the plate
 
-dirichlet_version = 1; % 0 = simplification, 1 = with multipliers, 2 = 
penalization
+variant = 2;           % 0 : not reduced, 1 : with reduced integration, 2 : 
MITC reduction
 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
+dirichlet_version = 1; % 0 = simplification, 1 = with multipliers, 2 = 
penalization
 
 plot_mesh = false;
 draw_solution = true;
 
 % trace on;
 gf_workspace('clear all');
-NX = 20;
+NX = 80;
 if (quadrangles)
   m = gf_mesh('cartesian',[0:1/NX:1],[0:1/NX:1]);
 else