Do local FEM on all RT elements

dune/dorie/common/flux_reconstruction/raviart_thomas/local_engine.hh

@@ -300,22 +300,21 @@ public:
   {
     auto& lfsu = lfsu_cache.localFunctionSpace();
 
-    const int order  = lfsu.finiteElement().localBasis().order();
-
     // solution of the local problem.
-    if (order > 1)
-    {
+    // const int order  = lfsu.finiteElement().localBasis().order();
+    // if (order > 1)
+    // {
       x_vec.resize(lfsu_cache.size(),0.0);
       m_matrix.solve(x_vec,r_vec);
 
       for (unsigned int i = 0; i < lfsu.size(); ++i)
         xl(lfsu,i) = x_vec[i];
       global_sl_view.write(xl);
-    } else {
-      for (unsigned int i = 0; i < lfsu.size(); ++i)
-        xl(lfsu,i) = r_vec[i];
-      global_sl_view.write(xl);
-    }
+    // } else {
+    //   for (unsigned int i = 0; i < lfsu.size(); ++i)
+    //     xl(lfsu,i) = r_vec[i];
+    //   global_sl_view.write(xl);
+    // }
 
   }
   //! @}
@@ -395,8 +394,8 @@ public:
     for (unsigned int i = 0; i < lfsv_volume.size(); ++i)
       r_vec[i] += rl(lfsv_volume,i);
 
-    const int order = lfsu.finiteElement().localBasis().order();
-    if (order > 1)
+    // const int order = lfsu.finiteElement().localBasis().order();
+    // if (order > 1)
       assembleVolume(eg, lfsu);
   }
 
@@ -448,8 +447,8 @@ public:
     }
 
     auto& lfsu = lfsu_s_cache.localFunctionSpace();
-    const int order = lfsu.finiteElement().localBasis().order();
-    if (order > 1)
+    // const int order = lfsu.finiteElement().localBasis().order();
+    // if (order > 1)
       assembleSkeleton(ig, lfsu);
   }
 
@@ -483,8 +482,8 @@ public:
     }
 
     auto& lfsu = lfsu_s_cache.localFunctionSpace();
-    const int order = lfsu.finiteElement().localBasis().order();
-    if (order > 1)
+    // const int order = lfsu.finiteElement().localBasis().order();
+    // if (order > 1)
       assembleSkeleton(ig, lfsu);
   }
 
@@ -528,6 +527,7 @@ public:
   template<class EG, class LFSU>
   void assembleVolume(const EG& eg, const LFSU& lfsu)
   {
+    auto& local_keys_u = lfsu.finiteElement().localCoefficients();
     const int order = lfsu.finiteElement().localBasis().order();
 
     using FESwitchTrial = Dune::FiniteElementInterfaceSwitch<
@@ -577,18 +577,21 @@ public:
       auto factor = it.weight() * gt.integrationElement(p_local);
 
       for (unsigned int i = 0; i < lfsv_volume.size(); i++)
-        for (unsigned int j = 0; j < lfsu.size(); j++)
-          m_matrix[i][j] += (phiu[j] * phiv[i]) * factor / (detB * detB);
+          for (unsigned int j = 0; j < lfsu.size(); j++)
+            if (local_keys_u.localKey(j).codim() == 0)
+              m_matrix[i][j] += (phiu[j] * phiv[i]) * factor / (detB * detB);
     }
   }
 
   template<class IG, class LFSU>
-  void assembleSkeleton(const IG& ig, const LFSU& lfsu_s)
+  void assembleSkeleton(const IG& ig, const LFSU& lfsu)
   {
+    auto& local_keys_u = lfsu.finiteElement().localCoefficients();
+
     unsigned int offset =
       lfsv_volume.size()
       + lfsv_skeleton.size() * ig.intersection().indexInInside();
-    const int order_s = lfsu_s.finiteElement().localBasis().order();
+    const int order_s = lfsu.finiteElement().localBasis().order();
 
     using FESwitchTrial = Dune::FiniteElementInterfaceSwitch<
       typename LFSU::Traits::FiniteElementType>;
@@ -612,32 +615,32 @@ public:
       const auto p_local_s = ig.geometryInInside().global(it.position());
 
       using LocalCoordinate = typename IG::LocalGeometry::LocalCoordinate;
-      LocalCoordinate sub_p_local_s;
+      LocalCoordinate sub_p_local;
+
+      auto sub_entity_id = ig.intersection().indexInInside();
 
       if (ig.intersection().conforming())
       {
-        sub_p_local_s = it.position();
+        sub_p_local = it.position();
       }
       else
       {
-        auto sub_entity_s_id = ig.intersection().indexInInside();
-        auto sub_entity_s = ig.inside().template subEntity<1>(sub_entity_s_id);
-        sub_p_local_s =
-          sub_entity_s.geometry().local(gtface.global(it.position()));
+        auto sub_entity = ig.inside().template subEntity<1>(sub_entity_id);
+        sub_p_local = sub_entity.geometry().local(gtface.global(it.position()));
       }
 
-      std::vector<LFSURange> phiu_s(lfsu_s.size());
+      std::vector<LFSURange> phiu_s(lfsu.size());
       std::vector<LFSVSkeletonRange> phiv_s(lfsv_skeleton.size());
 
       // evaluate basis functions
-      lfsu_s.finiteElement().localBasis().evaluateFunction(p_local_s, phiu_s);
-      lfsv_skeleton.finiteElement().localBasis().evaluateFunction(sub_p_local_s,
+      lfsu.finiteElement().localBasis().evaluateFunction(p_local_s, phiu_s);
+      lfsv_skeleton.finiteElement().localBasis().evaluateFunction(sub_p_local,
                                                                   phiv_s);
 
       auto n_F = ig.unitOuterNormal(it.position());
       auto BTransposed = ig.inside().geometry().jacobianTransposed(p_local_s);
       auto detB = BTransposed.determinant();
-      for (unsigned int i = 0; i < lfsu_s.size(); i++)
+      for (unsigned int i = 0; i < lfsu.size(); i++)
       {
         LFSURange y_s{ phiu_s[i] };
         BTransposed.mtv(y_s, phiu_s[i]);
@@ -648,9 +651,11 @@ public:
         it.weight() * ig.geometry().integrationElement(it.position());
 
       for (unsigned int i = 0; i < lfsv_skeleton.size(); i++)
-        for (unsigned int j = 0; j < lfsu_s.size(); j++)
-          m_matrix[i + offset][j] +=
-            (n_F * phiu_s[j]) * phiv_s[i] * factor / detB;
+          for (unsigned int j = 0; j < lfsu.size(); j++)
+            // if (local_keys_u.localKey(j).codim() == 1
+                // and local_keys_u.localKey(j).subEntity() == sub_entity_id)
+              m_matrix[i + offset][j] +=
+                (n_F * phiu_s[j]) * phiv_s[i] * factor / detB;
     }
   }