Added a general flux reconstruction method, unfortunately Raviart Thomas...

Added a general flux reconstruction method, unfortunately Raviart Thomas elements are not implemented in all of the dimensions and degrees we need
Signed-off-by: Santiago Ospina <santiago.ospina@iup.uni-heidelberg.de>

dune/dorie/interface/richards_simulation.cc

 #include "richards_simulation.hh"
+#include "../solver/util_flux_reconstruction.hh"
 
 namespace Dune{
 namespace Dorie{
@@ -213,18 +214,32 @@ void RichardsSimulation<Traits>::write_data (RF time)
 	if constexpr (Traits::write_output)
 	{
 		update_adapters();
+
+	std::shared_ptr<Dune::Dorie::FluxReconstruction<typename Traits::GV,
+		typename Traits::RangeField,
+		Traits::GV::dimension,
+		1,
+		Traits::GridGeometryType>> lift = std::make_shared<
+								Dune::Dorie::FluxReconstruction<typename Traits::GV,
+								typename Traits::RangeField,
+								Traits::GV::dimension,
+								1,
+								Traits::GridGeometryType>>(udgf);
+
 		if (inifile.get<bool>("output.vertexData")) {
 			vtkwriter->addVertexData(udgf,"head");
 			vtkwriter->addVertexData(fluxdgf,"flux");
 			vtkwriter->addVertexData(condgf,"K_0");
 			vtkwriter->addVertexData(waterdgf,"theta_w");
 			vtkwriter->addVertexData(satdgf,"Theta");
+			// vtkwriter->addVertexData(lift,"lift");
 		} else {
 			vtkwriter->addCellData(udgf,"head");
 			vtkwriter->addCellData(fluxdgf,"flux");
 			vtkwriter->addCellData(condgf,"K_0");
 			vtkwriter->addCellData(waterdgf,"theta_w");
 			vtkwriter->addCellData(satdgf,"Theta");
+			// vtkwriter->addCellData(lift,"lift");
 		}
 		try{
 			vtkwriter->write(time,output_type);

dune/dorie/solver/util_flux_reconstruction.hh

+// -*- tab-width: 4; indent-tabs-mode: nil -*-
+#ifndef DUNE_DORIE_FLUX_RECONSTRUCTION_HH
+#define DUNE_DORIE_FLUX_RECONSTRUCTION_HH
+
+#include <dune/pdelab/finiteelementmap/raviartthomasfem.hh>
+
+namespace Dune {
+namespace Dorie {
+
+/**
+ * @brief      Class for global lifting.
+ *
+ * @tparam     <unnamed>  { description }
+ */
+template<class GF>
+class GlobalLifting 
+  : public Dune::PDELab::LocalOperatorDefaultFlags
+{
+public:
+  GlobalLifting(std::shared_ptr<GF> grid_function)
+    : _grid_function(grid_function)
+  {}
+
+  enum {doAlphaSkeleton = true};
+  enum {doSkeletonTwoSided = true}; // FIXME, my programer was too lazy to do the two sides at once
+
+  /*-----------------------------------------------------------------------*//**
+   * @brief      Skeleton integral depending on test and ansatz functions. Each
+   *             face is only visited once since this method is symmetric.
+   *
+   * @param[in]  ig      The intersection entity of the grid (inside + outside
+   *                     entities)
+   * @param[in]  lfsu_i  The inside ansatz local function space
+   * @param[in]  x_i     The coefficients of the lfsu_i
+   * @param[in]  lfsv_i  The inside test local function space
+   * @param[in]  lfsu_o  The outside ansatz local function space
+   * @param[in]  x_o     The coefficients of the lfsu_o
+   * @param[in]  lfsv_o  The outside test local function space
+   * @param      r_i     The view of the residual vector w.r.t lfsu_i
+   * @param      r_o     The view of the residual vector w.r.t lfsu_o
+   *
+   * @tparam     IG      The type for ig
+   * @tparam     LFSU    The type for lfsu_i and lfsu_o
+   * @tparam     X       The type for x_i and x_o
+   * @tparam     LFSV    The type for lfsv_i and lfsv_o
+   * @tparam     R       The type for r_i and r_o
+   */
+  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
+  void alpha_skeleton(const IG& ig, 
+      const LFSU& lfsu_i, const X& x_i, const LFSV& lfsv_i,
+      const LFSU& lfsu_o, const X& x_o, const LFSV& lfsv_o, 
+      R& r_i, R& r_o) const
+  {
+
+    using RF = typename LFSU::Traits::FiniteElementType::
+          Traits::LocalBasisType::Traits::RangeFieldType;
+
+    const int intorder = 1; //FIXME
+    auto face_geo = ig.geometry();
+
+    assert(lfsu_i.size()==lfsv_i.size());
+    assert(lfsu_o.size()==lfsv_o.size());
+
+    int n = lfsu_i.size();
+    Dune::DynamicMatrix<RF> A(n,n);
+    Dune::DynamicVector<RF> b(n);
+
+    // loop over quadrature points and integrate normal flux
+    for (const auto& it : quadratureRule(face_geo,intorder)) {
+      for (std::size_t i = 0; i < lfsv_i.size(); ++i) {
+        for (std::size_t j = 0; j < lfsv_o.size(); ++j) {
+          if (lfsv_i.localIndex(i) == lfsv_o.localIndex(j)) {
+            const auto average = 0.5*(x_i(lfsv_i,i) + x_o(lfsv_o,j)); // FIXME: average of a vector field!
+            const auto normal = ig.unitOuterNormal(it.position());
+            typename GF::Traits::RangeType flux;
+            typename GF::Traits::DomainType  x_local = ig.geometryInInside().global(it.position());
+            _grid_function->evaluate(ig.inside(),x_local,flux);
+            b[i] += average * flux * normal ; // FIXME eval theta = jump in flux
+          }
+        }
+        for (std::size_t j = 0; j < lfsu_i.size(); ++j)
+          A[j][i] = x_i(lfsu_i,j)*x_i(lfsv_i,i); // Check order of i and j access
+      }
+    }
+    Dune::DynamicVector<RF> x(n);
+    A.solve(b,x);
+    for (std::size_t i = 0; i < lfsu_i.size(); ++i)
+      r_i.accumulate(lfsu_i,i,x[i]);
+  }
+private:
+  std::shared_ptr<GF> _grid_function;
+};
+
+template<class GV, class RF, int d, int k, Dune::GeometryType::BasicType gt>
+class FluxReconstruction : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<GV,RF,GV::dimension,
+                                                  Dune::FieldVector<RF,GV::dimension>>,FluxReconstruction<GV,RF,d,k,gt>>
+{
+public:
+  using Traits = Dune::PDELab::GridFunctionTraits<GV,RF,GV::dimension,
+                                                  Dune::FieldVector<RF,GV::dimension>>;
+private:
+  static constexpr int dim = d;
+  using DF = typename Traits::DomainFieldType;
+
+  using FEM = Dune::PDELab::RaviartThomasLocalFiniteElementMap<GV,DF,RF,k,gt>;
+  using CON = Dune::PDELab::NoConstraints;
+  using VBE = Dune::PDELab::ISTL::VectorBackend<>;
+  using GFS = Dune::PDELab::GridFunctionSpace<GV,FEM,CON,VBE>;
+  using X = Dune::PDELab::Backend::Vector<GFS,RF>;
+  using DGF = Dune::PDELab::DiscreteGridFunction<GFS,X>;
+
+public:
+  template<class GF>
+  FluxReconstruction(std::shared_ptr<GF> grid_function)
+    : _gv(grid_function->getGridView())
+    , _fem(_gv)
+    , _gfs(_gv,_fem)
+    , _r(_gfs,0.0)
+  {
+    using MBE = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
+    MBE mbe(5);
+
+    using LOP = GlobalLifting<GF>;
+    LOP lop(grid_function);
+    using GO = Dune::PDELab::GridOperator<GFS,GFS,LOP,MBE,DF,RF,RF>;
+    GO go(_gfs,_gfs,lop,mbe);
+    X x(_gfs,1.0);
+    using LS = Dune::PDELab::ISTLBackend_SEQ_BCGS_AMG_SSOR<GO>;
+    LS ls;
+    go.residual(x,_r);
+    _dfg = std::make_unique<DGF>(_gfs,_r);
+  }
+
+  void evaluate ( const typename Traits::ElementType& e,
+                  const typename Traits::DomainType& x,
+                        typename Traits::RangeType& y) const
+  {
+    // _dfg->evaluate(e,x,y);
+  }
+
+  const typename Traits::GridViewType& getGridView() const
+  {
+    return _gv;
+  }
+
+private:
+  GV _gv;
+  FEM _fem;
+  GFS _gfs;
+  X _r;
+  std::unique_ptr<DGF> _dfg;
+};
+
+} // namespace Dorie
+} // namespace Dune
+
+
+#endif // DUNE_DORIE_FLUX_RECONSTRUCTION_HH

dune/dorie/solver/util_writer.hh

@@ -44,7 +44,8 @@ public:
     static_assert(std::is_same<typename GF::Traits::GridViewType,GV>::value, 
       "GridFunctionVTKSequenceWriter only works with only one GridView type.");
 
-    auto vtk_gf_ptr = Dune::PDELab::makeVTKGridFunctionAdapter(gf_ptr,name);
+
+    auto vtk_gf_ptr = std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<GF> >(gf_ptr, name);
     Dune::VTKSequenceWriter<GV>::addVertexData(vtk_gf_ptr);
   }
 
@@ -62,7 +63,7 @@ public:
     static_assert(std::is_same<typename GF::Traits::GridViewType,GV>::value, 
       "GridFunctionVTKSequenceWriter only works with only one GridView type.");
     
-    auto vtk_gf_ptr = Dune::PDELab::makeVTKGridFunctionAdapter(gf_ptr,name);
+    auto vtk_gf_ptr = std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<GF> >(gf_ptr, name);
     Dune::VTKSequenceWriter<GV>::addCellData(vtk_gf_ptr);
   }
 

python/dorie/dorie/cli/cmds.py

@@ -9,8 +9,8 @@ except NameError:
     pass
 
 # paths set by cmake
-DORIEDIR = "/Users/lriedel/dune/dorie/build-cmake"
-MPIEXEC = "/usr/local/bin/mpiexec"
+DORIEDIR = "/home/saospina/Codes/DUNE_INSTALL/Debug/dorie"
+MPIEXEC = "/usr/bin/mpiexec"
 MPIEXEC_NUMPROC_FLAG = "-n"
 MPIEXEC_PREFLAG = ""
 MPIEXEC_POSTFLAGS = ""