[Getfem-commits] (no subject)

[Yves Renard]

branch: master
commit 172ffff70336a48798daa591efc8e2474dfdd3a9
Author: Yves Renard <address@hidden>
Date:   Sat Aug 31 16:18:10 2019 +0200

    minor fixes
---
 doc/sphinx/source/userdoc/hho.rst             |   2 +-
 interface/tests/python/demo_elasticity_HHO.py |  14 +-
 src/getfem_HHO.cc                             | 236 ++++++++++++++++++++++++++
 3 files changed, 244 insertions(+), 8 deletions(-)

diff --git a/doc/sphinx/source/userdoc/hho.rst 
b/doc/sphinx/source/userdoc/hho.rst
index 758d313..5f277c1 100644
--- a/doc/sphinx/source/userdoc/hho.rst
+++ b/doc/sphinx/source/userdoc/hho.rst
@@ -13,7 +13,7 @@ Tools for HHO (Hybrid High-Order) methods
 
 HHO method are hybrid methods in the sense that they have both degrees of 
freedom located on the element of a mesh and on the faces of the elements which 
represent separated approximations. HHO method are primal methods in the sense 
that both the degree of freedom in the element and on the faces represent the 
main unknown of the problem (no lagrange multipliers is introduced). The 
interest of these methods, first developped in  [Di-Er2015]_, [Di-Er2017]_ is 
their accuracy and their great [...]
 
-HHO methods can be applied to arbitrary shape elements. However, the 
implementation in |gf| is for the moment limited to standard elements : 
simplices, quadrilaterals, hexahedrons, ... Moreover this implementation is 
still experimental and not pretending to optimality. For the moment, no tool to 
make an auyomatic condensation of internal node is proposed and te method wotk 
only for simplices.
+HHO methods can be applied to arbitrary shape elements. However, the 
implementation in |gf| is for the moment limited to standard elements : 
simplices, quadrilaterals, hexahedrons, ... Moreover this implementation is 
still experimental and not pretending to optimality. For the moment, there is 
no tool to make an automatic condensation of internal dofs.
 
 HHO elements
 ------------
diff --git a/interface/tests/python/demo_elasticity_HHO.py 
b/interface/tests/python/demo_elasticity_HHO.py
index aaf6b07..b439056 100644
--- a/interface/tests/python/demo_elasticity_HHO.py
+++ b/interface/tests/python/demo_elasticity_HHO.py
@@ -65,10 +65,10 @@ mfrhs = gf.MeshFem(m, 1)
 
 if (using_HHO):
   if (use_quad):
-    mfu.set_fem(gf.Fem('FEM_HHO(FEM_QUAD_IPK(2,1),FEM_SIMPLEX_CIPK(1,1))'))
+    mfu.set_fem(gf.Fem('FEM_HHO(FEM_QUAD_IPK(2,2),FEM_SIMPLEX_CIPK(1,2))'))
     # 
mfu.set_fem(gf.Fem('FEM_HHO(FEM_QK_DISCONTINUOUS(2,2,0.1),FEM_SIMPLEX_CIPK(1,2))'))
     # mfu.set_fem(gf.Fem('FEM_HHO(FEM_QK(2,2),FEM_PK(1,2))'))
-    mfur.set_fem(gf.Fem('FEM_QUAD_IPK(2,2)'))
+    mfur.set_fem(gf.Fem('FEM_QUAD_IPK(2,3)'))
     # mfur.set_fem(gf.Fem('FEM_QK(2,3)'))
   else:
     mfu.set_fem(gf.Fem('FEM_HHO(FEM_SIMPLEX_IPK(2,2),FEM_SIMPLEX_CIPK(1,2))'))
@@ -82,9 +82,9 @@ else:
     mfur.set_fem(gf.Fem('FEM_PK(2,2)'))
     
 if (use_quad):
-  mfgu.set_fem(gf.Fem('FEM_QUAD_IPK(2,1)'))
-  # mfgu.set_fem(gf.Fem('FEM_QK(2,1)'))
-  mfrhs.set_fem(gf.Fem('FEM_QK(2,2)'))
+  mfgu.set_fem(gf.Fem('FEM_QUAD_IPK(2,2)'))
+  # mfgu.set_fem(gf.Fem('FEM_QK(2,2)'))
+  mfrhs.set_fem(gf.Fem('FEM_QK(2,3)'))
 else:
   mfgu.set_fem(gf.Fem('FEM_PK(2,2)'))
   mfrhs.set_fem(gf.Fem('FEM_PK(2,3)'))
@@ -158,7 +158,7 @@ if (using_HHO):
   md.add_linear_term(mim, 'clambda*Trace(HHO_Grad_u)*Trace(HHO_Grad_Test_u)'
                      +    '+ 2*cmu*Sym(HHO_Grad_u):Sym(HHO_Grad_Test_u)')
   # Stabilization term
-  md.add_linear_term(mim, 'cmu*HHO_Stab_u.HHO_Stab_Test_u', ALL_FACES)
+  md.add_linear_term(mim, '20*cmu*HHO_Stab_u.HHO_Stab_Test_u', ALL_FACES)
 else:
   md.add_isotropic_linearized_elasticity_brick(mim, 'u', 'clambda', 'cmu')
 
@@ -201,6 +201,6 @@ print('You can view the solution with (for example):')
 print('gmsh elasticity.pos')
 
 
-if (H1error > 0.013):
+if (H1error > 0.011):
     print('Error too large !')
     exit(1)
diff --git a/src/getfem_HHO.cc b/src/getfem_HHO.cc
index 52f08f0..0d451ea 100644
--- a/src/getfem_HHO.cc
+++ b/src/getfem_HHO.cc
@@ -644,6 +644,241 @@ namespace getfem {
     md.add_elementary_transformation(name, p);
   }
 
+#if 1 //  Old single mef version
+
+class _HHO_stabilization
+    : public virtual_elementary_transformation {
+
+  public:
+
+    virtual void give_transformation(const mesh_fem &mf1, const mesh_fem &mf2,
+                                     size_type cv, base_matrix &M) const {
+      // The reconstructed variable "S" is described on mf2 and computed by
+      // S(v) = P_{\dT}(v_{dT} - D(v)  - P_T(v_T - D(v)))
+      // where P__{\dT} et P_T are L2 projections on the boundary and on the
+      // interior of T on the corresponding discrete spaces.
+      // Note that P_{\dT}(v_{dT}) = v_{dT} and P_T(v_T) = v_T and D is
+      // the reconstructed value on P^{k+1} given by the formula:
+      //   \int_T Grad(D).Grad(w) =   \int_T Grad(v_T).Grad(w)
+      //                            + \int_{dT}(v_{dT} - v_T).(Grad(w).n)
+      // with the constraint
+      //   \int_T D = \int_T v_T
+      // where "w" is the test function arbitrary in mf2, "v_T" is the field
+      // inside the element whose gradient is to be reconstructed,
+      // "v_{dT}" is the field on the boundary of T and "n" is the outward
+      // unit normal.
+      // The implemented formula is
+      // S(v) = v_{dT} - P_{\dT}D(v) - P_{\dT}(v_T) + P_{\dT}(P_T(D(v)))
+      // by the mean of the projection matrix from P^{k+1} to the original 
space
+      // and the projection matrix from interior space to the boundary space.
+      // As it is built, S(v) is zero on interior dofs.
+      
+      GMM_ASSERT1(&mf1 == &mf2, "The HHO stabilization transformation is "
+                  "only defined on the HHO space to itself");
+
+      // Obtaining the fem descriptors
+      pfem pf1 = mf1.fem_of_element(cv);
+      short_type degree = pf1->estimated_degree();
+      bgeot::pgeometric_trans pgt = mf1.linked_mesh().trans_of_convex(cv);
+      pfem pf2 = classical_fem(pgt, short_type(degree + 1)); // Should be
+                                         // changed for an interior PK method
+      pfem pfi = interior_fem_of_hho_method(pf1);
+
+      papprox_integration pim
+        = classical_approx_im(pgt, dim_type(2*degree+2))->approx_method();
+
+      base_matrix G;
+      bgeot::vectors_to_base_matrix(G, mf1.linked_mesh().points_of_convex(cv));
+
+      bgeot::pgeotrans_precomp pgp
+        = HHO_pgp_pool(pgt, pim->pintegration_points());
+      pfem_precomp pfp1 = HHO_pfp_pool(pf1, pim->pintegration_points());
+      pfem_precomp pfp2 = HHO_pfp_pool(pf2, pim->pintegration_points());
+      pfem_precomp pfpi = HHO_pfp_pool(pfi, pim->pintegration_points());
+      
+      fem_interpolation_context ctx1(pgp, pfp1, 0, G, cv);
+      fem_interpolation_context ctx2(pgp, pfp2, 0, G, cv);
+      fem_interpolation_context ctxi(pgp, pfpi, 0, G, cv);
+
+      size_type Q = mf1.get_qdim(), N = mf1.linked_mesh().dim();
+      base_vector un(N);
+      size_type qmult1 =  Q / pf1->target_dim();
+      size_type ndof1 = pf1->nb_dof(cv) * qmult1;
+      size_type qmult2 =  Q / pf2->target_dim();
+      size_type ndof2 = pf2->nb_dof(cv) * qmult2;
+      size_type qmulti =  Q / pfi->target_dim();
+      size_type ndofi = pfi->nb_dof(cv) * qmulti;
+
+      
+      base_tensor t1, t2, ti, tv1, tv2, t1p, t2p;
+      base_matrix tv1p, tv2p, tvi;
+      base_matrix M1(ndof2, ndof1), M2(ndof2, ndof2), M2inv(ndof2, ndof2);
+      base_matrix M3(Q, ndof1), M4(Q, ndof2);
+      base_matrix aux1(ndof2, ndof1), aux2(ndof2, ndof2);
+      base_matrix M7(ndof1, ndof1), M7inv(ndof1, ndof1), M8(ndof1, ndof2);
+      base_matrix M9(ndof1, ndof1), MD(ndof2, ndof1);
+      scalar_type area(0);
+
+      // Integrals on the element : \int_T Grad(D).Grad(w) (M2)
+      //                            \int_T Grad(v_T).Grad(w) (M1)
+      //                            \int_T D (M4)  and \int_T v_T (M3)
+      for (size_type ipt = 0; ipt < pim->nb_points_on_convex(); ++ipt) {
+        ctx1.set_ii(ipt); ctx2.set_ii(ipt);
+        scalar_type coeff = pim->coeff(ipt) * ctx1.J();
+        area += coeff;
+        
+        ctx1.grad_base_value(t1);
+        vectorize_grad_base_tensor(t1, tv1, ndof1, pf1->target_dim(), Q);
+
+        ctx1.base_value(t1p);
+        vectorize_base_tensor(t1p, tv1p, ndof1, pf1->target_dim(), Q);
+
+        ctx2.grad_base_value(t2);
+        vectorize_grad_base_tensor(t2, tv2, ndof2, pf2->target_dim(), Q);
+
+        ctx2.base_value(t2p);
+        vectorize_base_tensor(t2p, tv2p, ndof2, pf2->target_dim(), Q);
+
+        for (size_type i = 0; i < ndof2; ++i) // To be optimized
+          for (size_type j = 0; j < ndof2; ++j)
+            for (size_type k = 0; k < Q*N; ++k)
+              M2(j, i) += coeff * tv2.as_vector()[i+k*ndof2]
+                                * tv2.as_vector()[j+k*ndof2];
+
+        for (size_type i = 0; i < ndof2; ++i)
+          for (size_type k = 0; k < Q; ++k)
+            M4(k,  i) += coeff * tv2p(i, k);
+              
+        for (size_type i = 0; i < ndof1; ++i) // To be optimized
+          for (size_type j = 0; j < ndof2; ++j)
+            for (size_type k = 0; k < Q*N; ++k)
+              M1(j, i) += coeff * tv1.as_vector()[i+k*ndof1]
+                                * tv2.as_vector()[j+k*ndof2];
+
+        for (size_type i = 0; i < ndof1; ++i)
+          for (size_type k = 0; k < Q; ++k)
+            M3(k,  i) += coeff * tv1p(i, k);
+
+        for (size_type i = 0; i < ndof1; ++i) // To be optimized
+          for (size_type j = 0; j < ndof1; ++j)
+            for (size_type k = 0; k < Q; ++k)
+              M7(i, j) += coeff * tv1p(i, k) * tv1p(j, k);
+        
+        for (size_type i = 0; i < ndof1; ++i) // To be optimized
+          for (size_type j = 0; j < ndof2; ++j)
+            for (size_type k = 0; k < Q; ++k)
+              M8(i, j) += coeff * tv1p(i, k) * tv2p(j, k);
+
+      }
+
+      // Integrals on the faces : \int_{dT}(v_{dT} - v_T).(Grad(w).n) (M1)
+      for (short_type ifc = 0; ifc < pgt->structure()->nb_faces(); ++ifc) {
+        ctx1.set_face_num(ifc); ctx2.set_face_num(ifc); ctxi.set_face_num(ifc);
+        size_type first_ind = pim->ind_first_point_on_face(ifc);
+        for (size_type ipt = 0; ipt < pim->nb_points_on_face(ifc); ++ipt) {
+          ctx1.set_ii(first_ind+ipt);
+          ctx2.set_ii(first_ind+ipt);
+          ctxi.set_ii(first_ind+ipt);
+          scalar_type coeff = pim->coeff(first_ind+ipt) * ctx1.J();
+          gmm::mult(ctx1.B(), pgt->normals()[ifc], un);
+          scalar_type normun = gmm::vect_norm2(un);
+          
+          ctx2.grad_base_value(t2);
+          vectorize_grad_base_tensor(t2, tv2, ndof2, pf2->target_dim(), Q);
+
+          ctx2.base_value(t2p);
+          vectorize_base_tensor(t2p, tv2p, ndof2, pf2->target_dim(), Q);
+
+          ctx1.base_value(t1);
+          vectorize_base_tensor(t1, tv1p, ndof1, pf1->target_dim(), Q);
+          
+          ctxi.base_value(ti);
+          vectorize_base_tensor(ti, tvi, ndofi, pfi->target_dim(), Q);
+
+
+          for (size_type i = 0; i < ndof1; ++i) // To be optimized
+            for (size_type j = 0; j < ndof2; ++j)
+              for (size_type k1 = 0; k1 < Q; ++k1) {
+                scalar_type b(0), a = coeff *
+                  (tv1p(i, k1) - (i < ndofi ? tvi(i, k1) : 0.));
+                for (size_type k2 = 0; k2 < N; ++k2)
+                  b += a * tv2.as_vector()[j+(k1 + k2*Q)*ndof2] * un[k2];
+                M1(j, i) += b;
+              }
+
+          for (size_type i = 0; i < ndof1; ++i) // To be optimized
+            for (size_type j = 0; j < ndof1; ++j)
+              for (size_type k = 0; k < Q; ++k)
+                M7(i, j) += coeff * normun * tv1p(i,k) * tv1p(j, k);
+
+          for (size_type i = 0; i < ndof1; ++i) // To be optimized
+            for (size_type j = 0; j < ndof2; ++j)
+              for (size_type k = 0; k < Q; ++k)
+                M8(i, j) += coeff * normun * tv1p(i,k) * tv2p(j, k);
+
+          for (size_type i = 0; i < ndof1; ++i) // To be optimized
+            for (size_type j = 0; j < ndofi; ++j)
+              for (size_type k = 0; k < Q; ++k)
+                M9(i, j) += coeff * normun * tv1p(i,k) * tvi(j, k); 
+        }
+      }
+
+      // Add the constraint with penalization
+      scalar_type coeff_p = pow(area, -1. - 2./scalar_type(N));
+      gmm::mult(gmm::transposed(M4), M4, aux2);
+      gmm::add (gmm::scaled(aux2, coeff_p), M2);
+      gmm::mult(gmm::transposed(M4), M3, aux1);
+      gmm::add (gmm::scaled(aux1, coeff_p), M1);
+
+      if (pf2->target_dim() == 1 && Q > 1) {
+        gmm::sub_slice I(0, ndof2/Q, Q);
+        gmm::lu_inverse(gmm::sub_matrix(M2, I, I));
+        for (size_type i = 0; i < Q; ++i) {
+          gmm::sub_slice I2(i, ndof2/Q, Q);
+          gmm::copy(gmm::sub_matrix(M2, I, I), gmm::sub_matrix(M2inv, I2, I2));
+        }
+      } else 
+        { gmm::copy(M2, M2inv); gmm::lu_inverse(M2inv); }
+      
+      if (pf1->target_dim() == 1 && Q > 1) {
+        gmm::sub_slice I(0, ndof1/Q, Q);
+        gmm::lu_inverse(gmm::sub_matrix(M7, I, I));
+        for (size_type i = 0; i < Q; ++i) {
+          gmm::sub_slice I2(i, ndof1/Q, Q);
+          gmm::copy(gmm::sub_matrix(M7, I, I), gmm::sub_matrix(M7inv, I2, I2));
+        }
+      } else
+        { gmm::copy(M7, M7inv); gmm::lu_inverse(M7inv); }
+      
+      gmm::mult(M2inv, M1, MD);
+      gmm::clean(MD, gmm::vect_norminf(MD.as_vector()) * 1E-13);
+
+      // S  = (I - inv(M7)*M9)(I - inv(M7)*M8*MD)
+      base_matrix MPB(ndof1, ndof1);
+      gmm::mult(M7inv, M9, MPB);
+      gmm::copy(gmm::identity_matrix(), M9);
+      gmm::add(gmm::scaled(MPB, scalar_type(-1)), M9);
+
+      base_matrix MPC(ndof1, ndof1), MPD(ndof1, ndof1);
+      gmm::mult(M8, MD, MPC);
+      gmm::mult(M7inv, MPC, MPD);
+      gmm::copy(gmm::identity_matrix(), M7);
+      gmm::add(gmm::scaled(MPD, scalar_type(-1)), M7);
+
+      gmm::mult(M9, M7, M);
+      gmm::clean(M, 1E-13);
+    }
+  };
+
+  void add_HHO_stabilization(model &md, std::string name) {
+    pelementary_transformation
+      p = std::make_shared<_HHO_stabilization>();
+    md.add_elementary_transformation(name, p);
+  }
+
+
+#else
+  
 
   class _HHO_stabilization
     : public virtual_elementary_transformation {
@@ -899,6 +1134,7 @@ namespace getfem {
     md.add_elementary_transformation(name, p);
   }
 
+#endif
 
   class _HHO_symmetrized_stabilization
     : public virtual_elementary_transformation {