[Getfem-commits] (no subject)

[Yves Renard]

branch: master
commit 92c168bc8cd62790417849f8f61acf235cc5868d
Author: Yves Renard <address@hidden>
Date:   Thu May 9 14:36:15 2019 +0200

    Adding a local P0 projection transformation
---
 interface/src/gf_model_set.cc         |  14 ++++-
 src/getfem/getfem_linearized_plates.h |   6 ++
 src/getfem_linearized_plates.cc       | 105 +++++++++++++++++++++++++++++++++-
 3 files changed, 122 insertions(+), 3 deletions(-)

diff --git a/interface/src/gf_model_set.cc b/interface/src/gf_model_set.cc
index 67e05e4..4b9bf6d 100644
--- a/interface/src/gf_model_set.cc
+++ b/interface/src/gf_model_set.cc
@@ -424,13 +424,25 @@ void gf_model_set(getfemint::mexargs_in& m_in,
        );
 
     /address@hidden ('add elementary rotated RT0 projection', @str transname)
-      Experimental method ... @*/
+      Add the elementary transformation corresponding to the projection
+      on rotated RT0 element for two-dimensional elements to the model.
+      The name is the name given to the elementary transformation. @*/
     sub_command
       ("add elementary rotated RT0 projection", 1, 1, 0, 0,
        std::string transname = in.pop().to_string();
        add_2D_rotated_RT0_projection(*md, transname);
        );
 
+    /address@hidden ('add elementary P0 projection', @str transname)
+      Add the elementary transformation corresponding to the projection
+      P0 element.
+      The name is the name given to the elementary transformation. @*/
+    sub_command
+      ("add P0 projection", 1, 1, 0, 0,
+       std::string transname = in.pop().to_string();
+       add_P0_projection(*md, transname);
+       );
+
     /address@hidden ('add interpolate transformation from expression', @str 
transname, @tmesh source_mesh, @tmesh target_mesh, @str expr)
       Add a transformation to the model from mesh `source_mesh` to mesh
       `target_mesh` given by the expression `expr` which corresponds to a
diff --git a/src/getfem/getfem_linearized_plates.h 
b/src/getfem/getfem_linearized_plates.h
index d8c90fc..17de751 100644
--- a/src/getfem/getfem_linearized_plates.h
+++ b/src/getfem/getfem_linearized_plates.h
@@ -48,6 +48,12 @@ namespace getfem {
   */
   void add_2D_rotated_RT0_projection(model &md, std::string name);
 
+  /** Add the elementary transformation corresponding to the projection
+      on P0 element.
+      The name is the name given to the elementary transformation.
+  */
+  void add_P0_projection(model &md, std::string name);
+
 
   /** Add a term corresponding to the classical Reissner-Mindlin plate
       model for which `u3` is the transverse displacement,
diff --git a/src/getfem_linearized_plates.cc b/src/getfem_linearized_plates.cc
index d7fb4c3..92080fe 100644
--- a/src/getfem_linearized_plates.cc
+++ b/src/getfem_linearized_plates.cc
@@ -189,8 +189,6 @@ namespace getfem {
     }
   };
 
-
-
   void add_2D_rotated_RT0_projection(model &md, std::string name) {
     pelementary_transformation
       p = std::make_shared<_2D_Rotated_RT0_projection_transformation>();
@@ -199,6 +197,109 @@ namespace getfem {
 
 
 
+  // Can be simplified ...
+  class _P0_projection_transformation
+    : public virtual_elementary_transformation {
+
+  public:
+
+    virtual void give_transformation(const mesh_fem &mf, size_type cv,
+                                     base_matrix &M) const{
+
+      THREAD_SAFE_STATIC base_matrix M_old;
+      THREAD_SAFE_STATIC pfem pf_old = nullptr;
+        
+      // Obtaining the fem descriptors
+      pfem pf1 = mf.fem_of_element(cv);
+      size_type N = 2;
+      // GMM_ASSERT1(pf1->dim() == 2, "This projection is only defined "
+      //             "for two-dimensional elements");
+      size_type qmult =  N / pf1->target_dim();
+      
+      bool simplex = false;
+      if (pf1->ref_convex(cv) == bgeot::simplex_of_reference(dim_type(N))) {
+        simplex = true;
+      } else if (pf1->ref_convex(cv)
+                 == bgeot::parallelepiped_of_reference(dim_type(N))) {
+        simplex = false;
+      } else {
+        GMM_ASSERT1(false, "Cannot adapt the method for such an element.");
+      }
+
+      if (pf1 == pf_old && pf1->is_equivalent() && M.size() == M_old.size()) {
+        gmm::copy(M_old, M);
+        return;
+      }
+
+      std::stringstream fem_desc;
+      fem_desc << "FEM_" << (simplex ? "P0":"Q0") << "(" << N << ")";
+      pfem pf2 = fem_descriptor(fem_desc.str());
+
+      // Obtaining a convenient integration method
+      size_type degree = pf1->estimated_degree() + pf2->estimated_degree();
+      bgeot::pgeometric_trans pgt = mf.linked_mesh().trans_of_convex(cv);
+      papprox_integration pim
+        = classical_approx_im(pgt, dim_type(degree))->approx_method();
+
+      // Computation of mass matrices
+      size_type ndof1 = pf1->nb_dof(cv) * qmult;
+      size_type ndof2 = pf2->nb_dof(0);
+      base_matrix M1(ndof1, ndof1), M2(ndof2, ndof2), B(ndof1, ndof2);
+      base_matrix aux0(ndof1, ndof1), aux1(ndof1, ndof2), aux2(ndof1, ndof2);
+      base_matrix aux3(ndof2, ndof2);
+
+      
+      base_matrix G;
+      bgeot::vectors_to_base_matrix(G, mf.linked_mesh().points_of_convex(cv));
+      fem_interpolation_context ctx1(pgt, pf1, base_node(N), G, cv);
+      fem_interpolation_context ctx2(pgt, pf2, base_node(N), G, cv);
+
+      base_tensor t1, t2;
+      base_matrix tv1, tv2;
+        
+      for (size_type i = 0; i < pim->nb_points_on_convex(); ++i) {
+
+        scalar_type coeff = pim->coeff(i); // Mult by ctx.J() not useful here
+        ctx1.set_xref(pim->point(i));
+        ctx2.set_xref(pim->point(i));    
+        pf1->real_base_value(ctx1, t1);
+        vectorize_base_tensor(t1, tv1, ndof1, pf1->target_dim(), N);
+        pf2->real_base_value(ctx2, t2);
+        vectorize_base_tensor(t2, tv2, ndof2, pf2->target_dim(), N);
+        for (size_type j = 0; j < ndof2; ++j) std::swap(tv2(j,0), tv2(j,1));
+       
+        gmm::mult(tv1, gmm::transposed(tv1), aux0);
+        gmm::add(gmm::scaled(aux0, coeff), M1);
+        gmm::mult(tv2, gmm::transposed(tv2), aux3);
+        gmm::add(gmm::scaled(aux3, coeff), M2);
+        gmm::mult(tv1, gmm::transposed(tv2), aux1);
+        gmm::add(gmm::scaled(aux1, coeff), B);
+      }
+      
+      
+      // Computation of M
+      gmm::lu_inverse(M1);
+      gmm::lu_inverse(M2);
+      gmm::mult(M1, B, aux1);
+      gmm::mult(aux1, M2, aux2);
+      GMM_ASSERT1(gmm::mat_nrows(M) == ndof1,
+                  "Element not convenient for projection");
+      gmm::mult(aux2, gmm::transposed(B), M);
+      gmm::clean(M, 1E-15);
+      M_old = M; pf_old = pf1;
+    }
+  };
+
+
+
+
+  void add_P0_projection(model &md, std::string name) {
+    pelementary_transformation
+      p = std::make_shared<_P0_projection_transformation>();
+    md.add_elementary_transformation(name, p);
+  }
+
+
 
 
   // RT0 projection in any dimension. Unused for the moment.