Added local solution of the raviart thomas local engine. In theory the engine...

Added local solution of the raviart thomas local engine. In theory the engine should be already working.
Signed-off-by: Santiago Ospina <santiago.ospina@iup.uni-heidelberg.de>

dune/dorie/solver/flux_reconstruction/localengeine.hh

@@ -7,6 +7,9 @@
 #include <dune/pdelab/constraints/common/constraints.hh>
 #include <dune/pdelab/localoperator/callswitch.hh>
 
+#include <dune/common/dynvector.hh>
+#include <dune/common/dynmatrix.hh>
+
 namespace Dune{
 namespace Dorie{
 
@@ -165,6 +168,7 @@ public:
     // Check that local matrix will be conforming
     assert(lfsu_cache.size() == lfsw_cache.size());
     assert(lfsu_cache.size() == (lfsv_volume.size() + lfsv_skeleton.size() * eg.entity().subEntities(1)));
+    assert(lfsv_skeleton.size() == lfsv_skeleton_n.size());
 
     global_sl_view.bind(lfsu_cache);
     xl.resize(lfsu_cache.size());
@@ -174,21 +178,34 @@ public:
   void onBindLFSV(const EG & eg, const LFSWC & lfsw_cache)
   {
     global_pl_view.bind(lfsw_cache);
-    pl.assign(lfsw_cache.size(),0.0);
+    pl.resize(lfsw_cache.size());
 
     lfsv_volume.bind( eg.entity() );
-    lfsv_skeleton.bind( eg.entity() );
 
     rl.assign(lfsv_volume.size(),0.0);
 
-    // TODO local matrix = 0
+    // variables for the local problem
+    auto n = lfsv_volume.size() + lfsv_skeleton.size() * eg.entity().subEntities(1);
+
+    // local system to be solved on each element. Matrix has block form:
+    // A_{eleme,eleme} A_{eleme,face1} A_{eleme,face2} ... A_{eleme,facen}
+    // A_{face1,eleme} A_{face1,face1} A_{face1,face2} ... A_{face1,facen}
+    // :                                                   :              
+    // A_{facen,eleme} A_{facen,face1} A_{facen,face2} ... A_{facen,facen}
+  
+    r_vec.resize(n,0.0);
+    m_matrix.resize(n,n,0.0);
   }
 
   template<typename IG, typename LFSUC, typename LFSWC>
   void onBindLFSUVInside(const IG & ig, const LFSUC & lfsu_cache, const LFSWC & lfsw_cache)
   {
+    lfsv_skeleton.bind( ig.inside() );
     global_sl_view.bind(lfsu_cache);
-    xl.resize(lfsu_cache.size());
+    xl.assign(lfsu_cache.size(),0.0);
+
+    global_pl_view.bind(lfsw_cache);
+    pl.resize(lfsw_cache.size());
   }
 
   template<typename IG, typename LFSUC, typename LFSWC>
@@ -196,14 +213,18 @@ public:
                           const LFSUC & lfsu_s_cache, const LFSWC & lfsw_s_cache,
                           const LFSUC & lfsu_n_cache, const LFSWC & lfsw_n_cache)
   {
+    lfsv_skeleton_n.bind( ig.outside() );
     global_sn_view.bind(lfsu_n_cache);
-    xn.resize(lfsu_n_cache.size());
+    xn.assign(lfsu_n_cache.size(),0.0);
+
+    global_pn_view.bind(lfsw_n_cache);
+    pn.resize(lfsw_n_cache.size());
   }
 
   template<typename IG, typename LFSWC>
   void onBindLFSVInside(const IG & ig, const LFSWC & lfsw_cache)
   {
-    global_pl_view.bind(lfsw_cache);
+    // global_rl_view.bind(lfsw_cache);
     rl.assign(lfsw_cache.size(),0.0);
   }
 
@@ -212,39 +233,22 @@ public:
                          const LFSWC & lfsw_s_cache,
                          const LFSWC & lfsw_n_cache)
   {
-    global_pn_view.bind(lfsw_n_cache);
+    // global_rn_view.bind(lfsw_n_cache);
     rn.assign(lfsw_n_cache.size(),0.0);
   }
 
   //! @}
 
-  //! Called when the local function space is about to be rebound or
-  //! discarded
-  //! @{
-  template<typename EG, typename LFSWC>
-  void onUnbindLFSV(const EG & eg, const LFSWC & lfsw_cache)
+  template<typename EG, typename LFSUC, typename LFSVC>
+  void onUnbindLFSUV(const EG & eg, 
+                     const LFSUC & lfsu_cache, 
+                     const LFSVC & lfsv_cache)
   {
-    //FIXME
-    // global_pl_view.add(rl);
-    // global_pl_view.commit();
-  }
+    auto& lfsu = lfsu_cache.localFunctionSpace();
 
-  template<typename IG, typename LFSWC>
-  void onUnbindLFSVInside(const IG & ig, const LFSWC & lfsw_cache)
-  {
-    //FIXME
-    // global_pl_view.add(rl);
-    // global_pl_view.commit();
-  }
-
-  template<typename IG, typename LFSWC>
-  void onUnbindLFSVOutside(const IG & ig,
-                           const LFSWC & lfsw_s_cache,
-                           const LFSWC & lfsw_n_cache)
-  {
-    //FIXME
-    // global_pn_view.add(rn);
-    // global_pn_view.commit();
+    for (unsigned int i = 0; i < lfsu.size(); ++i)
+      xl(lfsu,i) = x_vec[i];
+    global_sl_view.add(xl);
   }
   //! @}
 
@@ -258,7 +262,7 @@ public:
   template<typename LFSUC>
   void loadCoefficientsLFSUOutside(const LFSUC & lfsu_n_cache)
   {
-    global_sn_view.read(xn);
+    global_pn_view.read(pn);
   }
   template<typename LFSUC>
   void loadCoefficientsLFSUCoupling(const LFSUC & lfsu_c_cache)
@@ -298,9 +302,69 @@ public:
   template<typename EG, typename LFSUC, typename LFSWC>
   void assembleUVVolume(const EG & eg, const LFSUC & lfsu_cache, const LFSWC & lfsw_cache)
   {
+    auto& lfsw = lfsw_cache.localFunctionSpace();
+
     rl_view.setWeight(local_assembler.weight());
     Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaVolume>::
-      alpha_volume(lop,eg,lfsw_cache.localFunctionSpace(),pl,lfsv_volume,rl_view);
+      alpha_volume(lop,eg,lfsw,pl,lfsv_volume,rl_view);
+
+    // Load residual view into the local residual vector
+    for (unsigned int i = 0; i < lfsv_volume.size(); ++i)
+      r_vec[i] = rl(lfsv_volume,i);
+
+    auto& lfsu = lfsu_cache.localFunctionSpace();
+
+    using LFSU = typename LFSUC::LocalFunctionSpace;
+    using FESwitchTrial = Dune::FiniteElementInterfaceSwitch<typename LFSU::Traits::FiniteElementType>;
+    using BasisSwitchTrial = Dune::BasisInterfaceSwitch<typename FESwitchTrial::Basis>;
+    using LFSURange = typename BasisSwitchTrial::Range;
+
+    using FESwitchTest = Dune::FiniteElementInterfaceSwitch<typename LFSVVolume::Traits::FiniteElementType>;
+    using BasisSwitchTest = Dune::BasisInterfaceSwitch<typename FESwitchTest::Basis>;
+    using LFSVVolumeRange = typename BasisSwitchTest::Range;
+
+    // Calculate the mass matrix
+    int intorderadd = 2; // FIXME!
+    int quadrature_factor = 2; // FIXME!
+
+    auto gt = eg.geometry();
+
+    const int intorder = intorderadd+quadrature_factor*lfsu.finiteElement().localBasis().order();
+
+    for (const auto& it : quadratureRule(gt,intorder))
+    {
+     // evaluate position in element local and global coordinates
+      const auto p_local = it.position();
+
+      auto BTransposed = gt.jacobianTransposed(p_local);
+      auto detB = BTransposed.determinant();
+
+      std::vector<LFSURange> phiu(lfsu.size());
+      std::vector<LFSVVolumeRange> phiv(lfsv_volume.size());
+
+      // evaluate basis functions
+      lfsu.finiteElement().localBasis().evaluateFunction(p_local,phiu);
+      lfsv_volume.finiteElement().localBasis().evaluateFunction(p_local,phiv);
+
+      for (unsigned int i=0; i<lfsu.size(); i++)
+      { // Piola transformation without 1/det
+        LFSURange y(phiu[i]);
+        BTransposed.mtv(y,phiu[i]);
+      }
+
+      for (unsigned int i=0; i<lfsv_volume.size(); i++)
+      { // Piola transformation without 1/det
+        LFSVVolumeRange y(phiv[i]);
+        BTransposed.mtv(y,phiv[i]);
+      }
+
+      // integration factor
+      auto factor = it.weight() * gt.integrationElement(p_local);
+
+      for (unsigned int i = 0; i<lfsv_volume.size(); i++)
+        for (unsigned int j = 0; i<lfsu.size(); j++)
+          m_matrix[i][j] = (phiu[j]*phiv[i])/detB * factor/detB;
+    }
   }
 
   template<typename EG, typename LFSWC>
@@ -318,6 +382,59 @@ public:
                      lfsw_s_cache.localFunctionSpace(),pl,lfsv_skeleton,
                      lfsw_n_cache.localFunctionSpace(),pn,lfsv_skeleton_n,
                      rl_view,rn_view);
+
+    unsigned int offset = lfsv_volume.size() + lfsv_skeleton.size()*ig.intersectionIndex();
+    
+    // Load residual view into the local residual vector
+    for (unsigned int i = 0; i < lfsv_skeleton.size(); ++i)
+      r_vec[i+offset] = rl(lfsv_skeleton,i) + rn(lfsv_skeleton_n,i);
+
+    auto& lfsu_s = lfsu_s_cache.localFunctionSpace();
+
+    using LFSU = typename LFSUC::LocalFunctionSpace;
+    using FESwitchTrial = Dune::FiniteElementInterfaceSwitch<typename LFSU::Traits::FiniteElementType>;
+    using BasisSwitchTrial = Dune::BasisInterfaceSwitch<typename FESwitchTrial::Basis>;
+    using LFSURange = typename BasisSwitchTrial::Range;
+
+    using FESwitchTest = Dune::FiniteElementInterfaceSwitch<typename LFSVSkeleton::Traits::FiniteElementType>;
+    using BasisSwitchTest = Dune::BasisInterfaceSwitch<typename FESwitchTest::Basis>;
+    using LFSVSkeletonRange = typename BasisSwitchTest::Range;
+
+    // Calculate the mass matrix
+    int intorderadd = 2; // FIXME!
+    int quadrature_factor = 2; // FIXME!
+
+    // get polynomial degree
+    const int order_s  = lfsu_s.finiteElement().localBasis().order();
+    const int intorder = intorderadd + quadrature_factor * order_s;
+
+    auto gtface = ig.geometry();
+
+    for (const auto& it : quadratureRule(gtface,intorder))
+    {
+      // position of quadrature point in local coordinates of elements
+      const auto p_local_s = ig.geometryInInside().global(it.position());
+
+      Dune::FieldVector<SolutionElement,IG::Entity::mydimension> y;
+      auto n_F = ig.unitOuterNormal(it.position());
+      auto BTransposed = ig.inside().geometry().jacobianTransposed(p_local_s);
+      auto detB = BTransposed.determinant();
+      BTransposed.mv(n_F,y);
+
+      std::vector<LFSURange> phiu_s(lfsu_s.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(it.position(),phiv_s);
+
+      // integration factor
+      const auto factor = it.weight() * ig.geometry().integrationElement(it.position());
+
+      for (unsigned int i = 0; i<lfsv_skeleton.size(); i++)
+        for (unsigned int j = 0; i<lfsu_s.size(); j++)
+          m_matrix[i+offset][j] = (y*phiu_s[j])/detB * phiv_s[i+offset] * factor;
+    }
   }
 
   template<typename IG, typename LFSWC>
@@ -330,6 +447,58 @@ public:
     rl_view.setWeight(local_assembler.weight());
     Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaBoundary>::
       alpha_boundary(lop,ig,lfsw_s_cache.localFunctionSpace(),pl,lfsv_skeleton,rl_view);
+
+    unsigned int offset = lfsv_volume.size() + lfsv_skeleton.size()*ig.intersectionIndex();
+    // Load residual view into the local residual vector
+    for (unsigned int i = 0; i < lfsv_skeleton.size(); ++i)
+      r_vec[i+offset] = rl(lfsv_skeleton,i);
+
+    auto& lfsu_s = lfsu_s_cache.localFunctionSpace();
+
+    using LFSU = typename LFSUC::LocalFunctionSpace;
+    using FESwitchTrial = Dune::FiniteElementInterfaceSwitch<typename LFSU::Traits::FiniteElementType>;
+    using BasisSwitchTrial = Dune::BasisInterfaceSwitch<typename FESwitchTrial::Basis>;
+    using LFSURange = typename BasisSwitchTrial::Range;
+
+    using FESwitchTest = Dune::FiniteElementInterfaceSwitch<typename LFSVSkeleton::Traits::FiniteElementType>;
+    using BasisSwitchTest = Dune::BasisInterfaceSwitch<typename FESwitchTest::Basis>;
+    using LFSVSkeletonRange = typename BasisSwitchTest::Range;
+
+    // Calculate the mass matrix
+    int intorderadd = 2; // FIXME!
+    int quadrature_factor = 2; // FIXME!
+
+    // get polynomial degree
+    const int order_s  = lfsu_s.finiteElement().localBasis().order();
+    const int intorder = intorderadd + quadrature_factor * order_s;
+
+    auto gtface = ig.geometry();
+
+    for (const auto& it : quadratureRule(gtface,intorder))
+    {
+      // position of quadrature point in local coordinates of elements
+      const auto p_local_s = ig.geometryInInside().global(it.position());
+
+      Dune::FieldVector<SolutionElement,IG::Entity::mydimension> y;
+      auto n_F = ig.unitOuterNormal(it.position());
+      auto BTransposed = ig.inside().geometry().jacobianTransposed(p_local_s);
+      auto detB = BTransposed.determinant();
+      BTransposed.mv(n_F,y);
+
+      std::vector<LFSURange> phiu_s(lfsu_s.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(it.position(),phiv_s);
+
+      // integration factor
+      const auto factor = it.weight() * ig.geometry().integrationElement(it.position());
+
+      for (unsigned int i = 0; i<lfsv_skeleton.size(); i++)
+        for (unsigned int j = 0; i<lfsu_s.size(); j++)
+          m_matrix[i+offset][j] = (y*phiu_s[j])/detB * phiv_s[i+offset] * factor;
+    }
   }
 
   template<typename IG, typename LFSWC>
@@ -360,6 +529,9 @@ public:
     rl_view.setWeight(local_assembler.weight());
     Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaVolumePostSkeleton>::
       alpha_volume_post_skeleton(lop,eg,lfsw_cache.localFunctionSpace(),pl,lfsv_volume,rl_view);
+
+    x_vec.resize(lfsu_cache.size(),0.0);
+    m_matrix.solve(x_vec,r_vec);
   }
 
   template<typename EG, typename LFSWC>
@@ -393,7 +565,7 @@ private:
   typedef Dune::PDELab::LocalVector<ResidualElement, LocalTestSpaceTag> ResidualVector;
   typedef Dune::PDELab::LocalVector<PrescriptionElement, LocalTestSpaceTag> PrescriptionVector;
 
-  typedef Dune::PDELab::LocalMatrix<MassMatrixElement> MassMatrix;
+  // typedef Dune::PDELab::LocalMatrix<MassMatrixElement> MassMatrix;
 
   //! Inside local coefficients
   SolutionVector xl;
@@ -426,6 +598,14 @@ private:
   // local function spaces in neighbor
   LFSVSkeleton lfsv_skeleton_n;
 
+  typedef Dune::DynamicVector<ResidualElement> LocalResidualVector;
+  typedef Dune::DynamicVector<ResidualElement> LocalSolutionVector;
+  typedef Dune::DynamicMatrix<ResidualElement> LocalMassMatrix;
+
+  LocalResidualVector r_vec;
+  LocalSolutionVector x_vec;
+  LocalMassMatrix     m_matrix;
+
 }; // End of class LocalRaviartThomasAssemblerEngine
 
 } // namespace Dorie

dune/dorie/solver/operator_DG.hh

@@ -287,7 +287,7 @@ public:
         if constexpr (dimRangeTest == dim) {
           // (we assume lfsv = RTk)
           lfsv.finiteElement().localBasis().evaluateFunction(p_local,gradphiv);
-          // Piola transformation without 1/det
+          // Piola transformation
           auto BTransposed = gt.jacobianTransposed(p_local);
           auto detB = BTransposed.determinant();
           for (unsigned int i=0; i<lfsv.size(); i++) {