[Flux Reconstruction]

* Added test for conforming jumps in normal direction of the skeleton grid.
Signed-off-by: Santiago Ospina <santiago.ospina@iup.uni-heidelberg.de>

.gitlab-ci.yml

@@ -101,7 +101,8 @@ build:unit-tests:
   <<: *build-tests
   script:
     - CMAKE_FLAGS="$CMAKE_FLAGS
-                   -DCMAKE_BUILD_TYPE=Debug"
+                   -DCMAKE_BUILD_TYPE=Debug
+                   -DDORIE_CHECK_FLUX_RECONSTRUCTION_JUMPS=true"
       $DUNECONTROL --only=dorie configure
     - $DUNECONTROL --only=dorie make $MAKE_FLAGS build_unit_tests
 

dune/dorie/common/flux_reconstruction/check_normal_jump_skeleton.hh

+#ifndef DUNE_DORIE_CHECK_NORMAL_JUMP_SKELETON_HH
+#define DUNE_DORIE_CHECK_NORMAL_JUMP_SKELETON_HH
+
+namespace Dune {
+namespace Dorie {
+
+template<class GF>
+double check_normal_jump_skeleton(const GF& gf, int int_order)
+{
+  using Range = typename GF::Traits::RangeType;
+  constexpr int dimDomain = GF::Traits::dimDomain;
+  constexpr int dimRange = GF::Traits::dimRange;
+  static_assert(dimDomain==dimRange,
+    "Function 'check_normal_jump_skeleton' is only defined for vector fields "
+    "with same dimension as the domain");
+
+  auto& gv = gf.getGridView();
+  auto& is = gv.indexSet();
+
+  double norm2;
+
+  // Traverse grid view
+  for (const auto& element : elements(gv))
+  {
+    // Traverse intersections
+    unsigned int intersection_index = 0;
+    for(const auto& intersection : intersections(gv,element))
+    {
+      // Only for skeleton
+      if (intersection.neighbor())
+      {
+        // get entities
+        const auto& entity_f = intersection;
+        const auto& entity_i = intersection.inside();
+        const auto& entity_o = intersection.outside();
+
+        auto id_i = is.index(entity_i);
+        auto id_o = is.index(entity_o);
+
+        // Compute face only once
+        if (id_o > id_i)
+        {
+          // Get element geometry
+          auto geo_f = entity_f.geometry();
+          auto geo_i = entity_i.geometry();
+          auto geo_o = entity_o.geometry();
+          auto geo_in_i = entity_f.geometryInInside();
+          auto geo_in_o = entity_f.geometryInOutside();
+
+          // Loop over quadrature points and integrate normal flux
+          for (const auto& it : quadratureRule(geo_f,int_order))
+          {
+            const auto& position_f = it.position();
+
+            // position of quadrature point in local coordinates of elements
+            const auto position_i = geo_in_i.global(position_f);
+            const auto position_o = geo_in_o.global(position_f);
+
+            // face normal vector
+            const auto normal_f = entity_f.unitOuterNormal(position_f);
+            Range v_i, v_o;
+            gf.evaluate(entity_i,position_i,v_i);
+            gf.evaluate(entity_o,position_o,v_o);
+
+            // integration factor
+            auto factor = it.weight() * geo_f.integrationElement(position_f);
+
+            auto jump = normal_f*(v_i-v_o);
+            norm2 += jump*jump*factor;
+          }
+        }
+      }
+    }
+  }
+  return sqrt(norm2);
+}
+
+} // namespace Dorie
+} // namespace Dune
+
+#endif // DUNE_DORIE_CHECK_NORMAL_JUMP_SKELETON_HH
\ No newline at end of file

dune/dorie/common/flux_reconstruction/raviart_thomas.hh

@@ -10,13 +10,16 @@
 #include <dune/dorie/common/finite_element/raviart_thomas_volume_fem.hh>
 #include <dune/dorie/common/finite_element/raviart_thomas_fem.hh>
 
+
+#include <dune/dorie/common/flux_reconstruction/check_normal_jump_skeleton.hh>
+
 namespace Dune{
 namespace Dorie{
 
 /**
  * @brief      Helper for the raviart thomas flux reconstrucion.
  * @details    This helper exports a bool that says whether the flux
- *             reconstruction mtehod is available for the convination of the
+ *             reconstruction method is available for the combination of the
  *             used template parameters.
  *
  * @tparam     GO     Grid operator type.
@@ -162,6 +165,14 @@ public:
         local_raviart_thomas_engine.setPrescription(p);
         local_raviart_thomas_engine.setSolution(_x);
         global_assembler_gop.assemble(local_raviart_thomas_engine);
+
+#ifdef DORIE_CHECK_FLUX_RECONSTRUCTION_JUMPS
+        if (check_jumps(std::numeric_limits<double>::epsilon()))
+            DUNE_THROW(Dune::InvalidStateException, 
+                        "Flux reconstruction is non-conforming " <<
+                        "(i.e. has jumps between interfaces). " <<
+                        "This is only expected for non-conforming grids")
+#endif        
     }
 
     // evaluate velocity
@@ -176,12 +187,33 @@ public:
     void setTime (Time time_)
     {}
 
-  //! get a reference to the GridView
+  /*-----------------------------------------------------------------------*//**
+   * @brief      get a reference to the GridView
+   *
+   * @return     The grid view.
+   */
   inline const typename Traits::GridViewType& getGridView () const
   {
     return _gv;
   }
 
+  /*-----------------------------------------------------------------------*//**
+   * @brief      Check that jumps in skeleton are below some value.
+   * @details    Tolerance is compared against the
+   *             @f$ L^2(\Omega)\f$ norm of the normal jumps in the skeleton.
+   *
+   * @param[in]  tol   The tolerance.
+   *
+   * @return     { description_of_the_return_value }
+   */
+  bool check_jumps(double tol)
+  {
+    int quadrature_factor = 2; // FIXME!
+    int int_order = _int_order_add+quadrature_factor*order;
+    auto jumps_norm = check_normal_jump_skeleton(*this,int_order);
+    return Dune::FloatCmp::le(jumps_norm,tol);
+  }
+
 private:
     GO& _go;
     GV _gv;
@@ -190,6 +222,7 @@ private:
     Range _x;
     Dune::PDELab::DiscreteGridFunctionPiola<GFSU,Range> _dgfp;
     const int _int_order_add;
+    double jumps_norm;
 };
 
 } // namespace Dorie

python/dorie/dorie/cli/cmds.py

@@ -9,7 +9,7 @@ except NameError:
     pass
 
 # paths set by cmake
-DORIEDIR = "/home/saospina/Codes/DUNE_INSTALL/Debug/dorie"
+DORIEDIR = "/home/saospina/Codes/DUNE_INSTALL/Release/dorie"
 MPIEXEC = "/usr/bin/mpiexec"
 MPIEXEC_NUMPROC_FLAG = "-n"
 MPIEXEC_PREFLAG = ""