[Getfem-commits] r5470 - in /trunk/getfem: doc/sphinx/source/userdoc/ src/

[Yves . Renard]

Author: renard
Date: Mon Nov  7 14:19:14 2016
New Revision: 5470

URL: http://svn.gna.org/viewcvs/getfem?rev=5470&view=rev
Log:
fix on incomplete Q2 geometric transformation

Modified:
    trunk/getfem/doc/sphinx/source/userdoc/model_contact_friction.rst
    trunk/getfem/src/bgeot_convex_ref.cc
    trunk/getfem/src/getfem_mesh_slicers.cc

Modified: trunk/getfem/doc/sphinx/source/userdoc/model_contact_friction.rst
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/doc/sphinx/source/userdoc/model_contact_friction.rst?rev=5470&r1=5469&r2=5470&view=diff
==============================================================================
--- trunk/getfem/doc/sphinx/source/userdoc/model_contact_friction.rst   
(original)
+++ trunk/getfem/doc/sphinx/source/userdoc/model_contact_friction.rst   Mon Nov 
 7 14:19:14 2016
@@ -84,13 +84,7 @@
 where :math:`\alpha_i` is a parameter which can be added for the 
homogenization of the augmentation parameter, :math:`(B_T U)_i` denotes here 
the sub-vector of indices from :math:`(d-1)(i-1)+1` to :math:`(d-1)i` for the 
sake of simplicity and the sliding velocity :math:`B_T \dot{U}` have been 
discretized into :math:`\frac{(B_T U - B_T U^{0})}{\Delta t}` with 
:math:`U^{0}` the displacement at the previous time step. Note that of course 
another discretization of the sliding velocity is possible and that the time 
step :math:`\Delta t` do not appear in the expression of the friction condition 
since it does not influence the direction of the sliding velocity.
 
 
-In that case, the homogenization coefficient :math:`\alpha_i` can be taken
-
-.. math::
-
-  \alpha_i = \frac{\int_{\Gamma_c} \varphi_i d\Gamma}{\ell}
-
-where :math:`\Gamma_c` is the contact boundary, :math:`\varphi_i` is the 
displacement shape function corresponding to the node :math:`a_i` and 
:math:`\ell` is a characteristic length, for instance the radius of the domain. 
In this way, the augmentation parameter :math:`r` can be expressed in 
:math:`N/m^2` and chosen closed to the Young modulus of the elastic body. Note 
that the solution is not very sensitive to the value of the augmentation 
parameter.
+In that case, the homogenization coefficient :math:`\alpha_i` can be taken 
proportional to :math:`h^{d-2}` (:math:`h` being the diameter of the element). 
In this way, the augmentation parameter :math:`r` can be expressed in 
:math:`N/m^2` and chosen closed to the Young modulus of the elastic body. Note 
that the solution is not sensitive to the value of the augmentation parameter.
 
 
 Weak nodal contact condition
@@ -134,13 +128,7 @@
 
   -\frac{1}{r\alpha_i}(\lambda_T^i - P_{{\mathscr B}(-{\mathscr F}P_{]-\infty, 
0]}(\lambda_N^i - \alpha_i r ((B_N U)_i - \text{gap}_i)))}(\lambda_T^i - 
\alpha_i r (B_T U - B_T U^{0})_i)) = 0, ~~ i = 1..N_c,
 
-except that now :math:`\lambda_N^i` and :math:`\lambda_T^i` are force 
densities, and a good value for :math:`\alpha_i` is now
-
-.. math::
-
-  \alpha_i = \frac{1}{\ell \int_{\Gamma_c}\psi_i},
-
-where :math:`\psi_i` is the shape function of the multiplier for the node 
:math:`a_i`. In that case, the augmentation parameter :math:`r` can still be 
chosen close to the Young modulus of the elastic body.
+except that now :math:`\lambda_N^i` and :math:`\lambda_T^i` are force 
densities, and :math:`\alpha_i` has to be now chosen proportional to 
:math:`1/h^d` such that the augmentation parameter :math:`r` can still be 
chosen close to the Young modulus of the elastic body.
 
 
 Note that without additional stabilization technique (see [HI-RE2010]_) an 
inf-sup condition have to be satisfied between the finite element of the 
displacement and the one for the multipliers. This means in particular that the 
finite element for the multiplier have to be "less rich" than the one for the 
displacement.
@@ -190,12 +178,12 @@
 .. math::
 
   \left\{\begin{array}{l}
-  a(u^h, v^h) + \displaystyle \int_{\Gamma_c} \lambda^h \cdot v^h d\Gamma = 
l(v^h) ~~~~ \forall v^h \in V^h, \\
+  a(u^h, v^h) + \displaystyle \int_{\Gamma_c} \lambda^h \cdot v^h d\Gamma = 
\ell(v^h) ~~~~ \forall v^h \in V^h, \\
   \displaystyle -\frac{1}{r}\int_{\Gamma_c} (\lambda^h_N + (\lambda^h_N - 
r(u^h_N-gap))_-)\mu^h_N d\Gamma \\
   ~~~~~~~~~~\displaystyle -\frac{1}{r}\int_{\Gamma_c} (\lambda^h_T 
-P_{B(\rho)}(\lambda^h_T - r\alpha(u^h_T-w^h_T)))\cdot \mu^h_T d\Gamma = 0 ~~~~ 
\forall \mu^h \in W^h,
   \end{array}\right.
 
-where :math:`a(\cdot, \cdot)` and :math:`l(v)` represent the remaining parts 
of the problem in  :math:`u`, for instance linear elasticity and 
:math:`\rho={\mathscr F}(\lambda^h_N - r(u^h_N-gap))_-`. In order to write a 
Newton iteration, one has to derive the tangent system. It can be written, 
reporting only the contact and friction terms and not the right hand side:
+where :math:`a(\cdot, \cdot)` and :math:`\ell(v)` represent the remaining 
parts of the problem in  :math:`u`, for instance linear elasticity and 
:math:`\rho={\mathscr F}(\lambda^h_N - r(u^h_N-gap))_-`. Note that in this 
case, the mathematical analysis leads to choose a value for the augmentation 
parameter of the kind :math:`r = r_0 / r` with :math:`r_0` having the dimension 
of a elasticity modulus (a classical choice is the value of Young's modulus). 
In order to write a Newton iteration, one has to derive the tangent system. It 
can be written, reporting only the contact and friction terms and not the right 
hand side:
 
 .. math::
 
@@ -218,7 +206,7 @@
 
   \left\{\begin{array}{l}
   a(u^h, v^h) + \displaystyle \int_{\Gamma_c} (\lambda^h_N - r(u^h_N-gap))_- 
v^h_N d\Gamma \\
-  ~~~~~~ - \displaystyle \int_{\Gamma_c} P_{B(\rho)}(\lambda^h_T - 
r\alpha(u^h_T-w^h_T)))\cdot v^h_T d\Gamma = l(v^h) ~~~~ \forall v^h \in V^h, \\
+  ~~~~~~ - \displaystyle \int_{\Gamma_c} P_{B(\rho)}(\lambda^h_T - 
r\alpha(u^h_T-w^h_T)))\cdot v^h_T d\Gamma = \ell(v^h) ~~~~ \forall v^h \in V^h, 
\\
   \displaystyle -\frac{1}{r}\int_{\Gamma_c} (\lambda^h_N + (\lambda^h_N - 
r(u^h_N-gap))_-)\mu^h_N d\Gamma \\
   ~~~~~~~~~~\displaystyle -\frac{1}{r}\int_{\Gamma_c} (\lambda^h_T 
-P_{B(\rho)}(\lambda^h_T - r\alpha(u^h_T-w^h_T)))\cdot \mu^h_T d\Gamma = 0 ~~~~ 
\forall \mu^h \in W^h,
   \end{array}\right.
@@ -249,7 +237,7 @@
 
   \left\{\begin{array}{l}
   a(u^h, v^h) + \displaystyle \int_{\Gamma_c} r(u^h_N-gap)_+ v^h_N d\Gamma \\
-  ~~~~~~ + \displaystyle \int_{\Gamma_c} P_{B(\mathscr F 
r(u^h_N-gap)_+)}(r\alpha(u^h_T-w^h_T))\cdot v^h_T d\Gamma = l(v^h) ~~~~ \forall 
v^h \in V^h,
+  ~~~~~~ + \displaystyle \int_{\Gamma_c} P_{B(\mathscr F 
r(u^h_N-gap)_+)}(r\alpha(u^h_T-w^h_T))\cdot v^h_T d\Gamma = \ell(v^h) ~~~~ 
\forall v^h \in V^h,
   \end{array}\right.
 
 and the tangent system:

Modified: trunk/getfem/src/bgeot_convex_ref.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/bgeot_convex_ref.cc?rev=5470&r1=5469&r2=5470&view=diff
==============================================================================
--- trunk/getfem/src/bgeot_convex_ref.cc        (original)
+++ trunk/getfem/src/bgeot_convex_ref.cc        Mon Nov  7 14:19:14 2016
@@ -214,29 +214,16 @@
   /* By Yao Koutsawa  <address@hidden> 2012-12-10                  */
 
   class Q2_incomplete_of_ref_ : public convex_of_reference {
+    pconvex_ref pllref;
   public :
-    scalar_type is_in(const base_node& pt) const {
-      // return a negative or null number if pt is in the convex
-      GMM_ASSERT1(pt.size() == cvs->dim(),
-                  "Q2_incomplete_of_ref_::is_in: Dimension does not match");
-      scalar_type e = -1.0, r = (pt.size() > 0) ? -pt[0] : 0.0;
-      base_node::const_iterator it = pt.begin(), ite = pt.end();
-      for (; it != ite; e += *it, ++it) r = std::max(r, -(*it));
-      return std::max(r, e);
-    }
-    scalar_type is_in_face(short_type f, const base_node& pt) const {
-      // return a null number if pt is in the face of the convex
-      // negative if the point is on the side of the face where the element is
-      GMM_ASSERT1(pt.size() == cvs->dim(), "Q2_incomplete_of_ref_::is_in_face: 
"
-                  "Dimension does not match");
-      if (f > 0) return -pt[f-1];
-      scalar_type e = -1.0;
-      base_node::const_iterator it = pt.begin(), ite = pt.end();
-      for (; it != ite; e += *it, ++it) {};
-      return e / sqrt(scalar_type(pt.size()));
-    }
+    scalar_type is_in(const base_node& pt) const
+    { return pllref->is_in(pt); }
+    scalar_type is_in_face(short_type f, const base_node& pt) const
+    { return pllref->is_in_face(f, pt); }
     
     Q2_incomplete_of_ref_(dim_type nc) {
+      GMM_ASSERT1(nc == 2 || nc == 3, "Sorry exist only in dimension 2 or 3");
+      pllref = parallelepiped_of_reference(nc);
       cvs = Q2_incomplete_structure(nc);
       convex<base_node>::points().resize(cvs->nb_points());
       normals_.resize(nc == 2 ? 4: 6);

Modified: trunk/getfem/src/getfem_mesh_slicers.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_mesh_slicers.cc?rev=5470&r1=5469&r2=5470&view=diff
==============================================================================
--- trunk/getfem/src/getfem_mesh_slicers.cc     (original)
+++ trunk/getfem/src/getfem_mesh_slicers.cc     Mon Nov  7 14:19:14 2016
@@ -655,11 +655,13 @@
     exec_(&n, 0, cvlst);
   }
 
-  void mesh_slicer::exec(const std::vector<short_type> &nrefine, const 
mesh_region& cvlst) {
+  void mesh_slicer::exec(const std::vector<short_type> &nrefine,
+                        const mesh_region& cvlst) {
     exec_(&nrefine[0], 1, cvlst);
   }
 
-  static bool check_orient(size_type cv, bgeot::pgeometric_trans pgt, const 
mesh& m) {
+  static bool check_orient(size_type cv, bgeot::pgeometric_trans pgt,
+                          const mesh& m) {
     if (pgt->dim() == m.dim() && m.dim()>=2) { /* no orient check for 
                                                   convexes of lower dim */
       base_matrix G; bgeot::vectors_to_base_matrix(G,m.points_of_convex(cv));
@@ -677,7 +679,8 @@
   }
 
 #if OLD_MESH_SLICER
-  void mesh_slicer::exec_(const short_type *pnrefine, int nref_stride, const 
mesh_region& cvlst) {
+  void mesh_slicer::exec_(const short_type *pnrefine, int nref_stride,
+                         const mesh_region& cvlst) {
     std::vector<base_node> cvm_pts;
     const bgeot::basic_mesh *cvm = 0;
     const bgeot::mesh_structure *cvms = 0;
@@ -829,16 +832,14 @@
       bool revert_orientation = check_orient(cv, pgt,m);
 
       /* update structure-dependent data */
-      
       /* TODO : fix levelset handling when slicing faces .. */
-      if (prev_cvr != cvr || nrefine != prev_nrefine || discont || 
prev_discont) {
+      if (prev_cvr != cvr || nrefine != prev_nrefine
+         || discont || prev_discont) {
         if (discont) {
-          //if (prev_cv != it.cv())
           cvm = &refined_simplex_mesh_for_convex_cut_by_level_set
             (mls->mesh_of_convex(cv), unsigned(nrefine));
         } else {
-          cvm = bgeot::refined_simplex_mesh_for_convex(cvr,
-                                                       short_type(nrefine));
+          cvm = 
bgeot::refined_simplex_mesh_for_convex(cvr,short_type(nrefine));
         }
         cvm_pts.resize(cvm->points().card());
         std::copy(cvm->points().begin(), cvm->points().end(), cvm_pts.begin());
@@ -848,15 +849,14 @@
       }
       if (face < dim_type(-1)) {
         if (!discont) {
-          cvms = bgeot::refined_simplex_mesh_for_convex_faces(cvr, 
short_type(nrefine))[face].get();
+          cvms = bgeot::refined_simplex_mesh_for_convex_faces
+           (cvr, short_type(nrefine))[face].get();
         } else {
           cvms = &refined_simplex_mesh_for_convex_faces_cut_by_level_set(face);
         }
       } else {
         cvms = cvm; 
       }
-
-
 
       /* apply the initial geometric transformation */
       std::vector<size_type> ptsid(cvm_pts.size()); std::fill(ptsid.begin(), 
ptsid.end(), size_type(-1));
@@ -880,7 +880,6 @@
           }
           G /= scalar_type(simplexes[snum].inodes.size());
         }
-          
 
         for (std::vector<size_type>::iterator itp = 
                simplexes[snum].inodes.begin();