Adds a fake local function space and a corrected assembler.

Signed-off-by: Santiago Ospina <santiago.ospina@iup.uni-heidelberg.de>

dune/dorie/solver/flux_reconstruction/localengeine.hh

@@ -10,6 +10,8 @@
 #include <dune/common/dynvector.hh>
 #include <dune/common/dynmatrix.hh>
 
+#include "raviart_thomas_operator/localfunctionspace.hh"
+
 namespace Dune{
 namespace Dorie{
 
@@ -657,8 +659,8 @@ private:
   const std::shared_ptr<const GFSVSkeleton> pgfsv_skeleton;
 
   /* local function spaces */
-  typedef Dune::PDELab::LocalFunctionSpace<GFSVVolume, LocalTestSpaceTag> LFSVVolume;
-  typedef Dune::PDELab::LocalFunctionSpace<GFSVSkeleton, LocalTestSpaceTag> LFSVSkeleton;
+  typedef Dune::Dorie::FakeLocalFunctionSpace<GFSVVolume, LocalTestSpaceTag> LFSVVolume;
+  typedef Dune::Dorie::FakeLocalFunctionSpace<GFSVSkeleton, LocalTestSpaceTag> LFSVSkeleton;
 
   // local function spaces in local cell
   LFSVVolume lfsv_volume;

dune/dorie/solver/flux_reconstruction/raviart_thomas_operator/assembler.hh

+#ifndef DUNE_DORIE_GRIDOPERATOR_DEFAULT_ASSEMBLER_HH
+#define DUNE_DORIE_GRIDOPERATOR_DEFAULT_ASSEMBLER_HH
+
+/* This is a modification of the DefaultAssembler in PDELab. It binds and 
+ * unbinds inside local function spaces with the intersection. Notice that to 
+ * merge this into PDELab, it is necessary to remove the inside bind/unbind 
+ * methods in all the local assamblers! They are not used currently, but if we 
+ * add this changes, local function spaces will be binded several times for the
+ * same entity. */
+
+#include <dune/common/typetraits.hh>
+#include <dune/pdelab/gridoperator/common/assemblerutilities.hh>
+#include <dune/pdelab/gridfunctionspace/localfunctionspace.hh>
+#include <dune/pdelab/gridfunctionspace/lfsindexcache.hh>
+#include <dune/pdelab/common/geometrywrapper.hh>
+#include <dune/pdelab/common/intersectiontype.hh>
+
+namespace Dune{
+  namespace Dorie{
+
+    namespace { using namespace Dune::PDELab;}
+    /**
+       \brief The assembler for standard DUNE grid
+
+       * \tparam GFSU GridFunctionSpace for ansatz functions
+       * \tparam GFSV GridFunctionSpace for test functions
+       */
+
+    template<typename GFSU, typename GFSV, typename CU, typename CV>
+    class DefaultAssembler {
+    public:
+
+      //! Types related to current grid view
+      //! @{
+      using EntitySet = typename GFSU::Traits::EntitySet;
+      using Element = typename EntitySet::Element;
+      using Intersection = typename EntitySet::Intersection;
+      //! @}
+
+      //! Grid function spaces
+      //! @{
+      typedef GFSU TrialGridFunctionSpace;
+      typedef GFSV TestGridFunctionSpace;
+      //! @}
+
+      //! Size type as used in grid function space
+      typedef typename GFSU::Traits::SizeType SizeType;
+
+      //! Static check on whether this is a Galerkin method
+      static const bool isGalerkinMethod = std::is_same<GFSU,GFSV>::value;
+
+      DefaultAssembler (const GFSU& gfsu_, const GFSV& gfsv_, const CU& cu_, const CV& cv_)
+        : gfsu(gfsu_)
+        , gfsv(gfsv_)
+        , cu(cu_)
+        , cv(cv_)
+        , lfsu(gfsu_)
+        , lfsv(gfsv_)
+        , lfsun(gfsu_)
+        , lfsvn(gfsv_)
+      { }
+
+      DefaultAssembler (const GFSU& gfsu_, const GFSV& gfsv_)
+        : gfsu(gfsu_)
+        , gfsv(gfsv_)
+        , cu()
+        , cv()
+        , lfsu(gfsu_)
+        , lfsv(gfsv_)
+        , lfsun(gfsu_)
+        , lfsvn(gfsv_)
+      { }
+
+      //! Get the trial grid function space
+      const GFSU& trialGridFunctionSpace() const
+      {
+        return gfsu;
+      }
+
+      //! Get the test grid function space
+      const GFSV& testGridFunctionSpace() const
+      {
+        return gfsv;
+      }
+
+      // Assembler (const GFSU& gfsu_, const GFSV& gfsv_)
+      //   : gfsu(gfsu_), gfsv(gfsv_), lfsu(gfsu_), lfsv(gfsv_),
+      //     lfsun(gfsu_), lfsvn(gfsv_),
+      //     sub_triangulation(ST(gfsu_.gridview(),Dune::PDELab::NoSubTriangulationImp()))
+      // { }
+
+      template<class LocalAssemblerEngine>
+      void assemble(LocalAssemblerEngine & assembler_engine) const
+      {
+        typedef LFSIndexCache<LFSU,CU> LFSUCache;
+
+        typedef LFSIndexCache<LFSV,CV> LFSVCache;
+
+        const bool needs_constraints_caching = assembler_engine.needsConstraintsCaching(cu,cv);
+
+        LFSUCache lfsu_cache(lfsu,cu,needs_constraints_caching);
+        LFSVCache lfsv_cache(lfsv,cv,needs_constraints_caching);
+        LFSUCache lfsun_cache(lfsun,cu,needs_constraints_caching);
+        LFSVCache lfsvn_cache(lfsvn,cv,needs_constraints_caching);
+
+        // Notify assembler engine about oncoming assembly
+        assembler_engine.preAssembly();
+
+        // Extract integration requirements from the local assembler
+        const bool require_uv_skeleton = assembler_engine.requireUVSkeleton();
+        const bool require_v_skeleton = assembler_engine.requireVSkeleton();
+        const bool require_uv_boundary = assembler_engine.requireUVBoundary();
+        const bool require_v_boundary = assembler_engine.requireVBoundary();
+        const bool require_uv_processor = assembler_engine.requireUVBoundary();
+        const bool require_v_processor = assembler_engine.requireVBoundary();
+        const bool require_uv_post_skeleton = assembler_engine.requireUVVolumePostSkeleton();
+        const bool require_v_post_skeleton = assembler_engine.requireVVolumePostSkeleton();
+        const bool require_skeleton_two_sided = assembler_engine.requireSkeletonTwoSided();
+
+        auto entity_set = gfsu.entitySet();
+        auto& index_set = entity_set.indexSet();
+
+        // Traverse grid view
+        for (const auto& element : elements(entity_set))
+          {
+            // Compute unique id
+            auto ids = index_set.uniqueIndex(element);
+
+            ElementGeometry<Element> eg(element);
+
+            if(assembler_engine.assembleCell(eg))
+              continue;
+
+            // Bind local test function space to element
+            lfsv.bind( element );
+            lfsv_cache.update();
+
+            // Notify assembler engine about bind
+            assembler_engine.onBindLFSV(eg,lfsv_cache);
+
+            // Volume integration
+            assembler_engine.assembleVVolume(eg,lfsv_cache);
+
+            // Bind local trial function space to element
+            lfsu.bind( element );
+            lfsu_cache.update();
+
+            // Notify assembler engine about bind
+            assembler_engine.onBindLFSUV(eg,lfsu_cache,lfsv_cache);
+
+            // Load coefficients of local functions
+            assembler_engine.loadCoefficientsLFSUInside(lfsu_cache);
+
+            // Volume integration
+            assembler_engine.assembleUVVolume(eg,lfsu_cache,lfsv_cache);
+
+            // Skip if no intersection iterator is needed
+            if (require_uv_skeleton || require_v_skeleton ||
+                require_uv_boundary || require_v_boundary ||
+                require_uv_processor || require_v_processor)
+              {
+                // Traverse intersections
+                unsigned int intersection_index = 0;
+                for(const auto& intersection : intersections(entity_set,element))
+                  {
+
+                    IntersectionGeometry<Intersection> ig(intersection,intersection_index);
+
+                    auto intersection_data = classifyIntersection(entity_set,intersection);
+                    auto intersection_type = std::get<0>(intersection_data);
+                    auto& outside_element = std::get<1>(intersection_data);
+
+                    switch (intersection_type)
+                      {
+                      case IntersectionType::skeleton:
+                        // the specific ordering of the if-statements in the old code caused periodic
+                        // boundary intersection to be handled the same as skeleton intersections
+                      case IntersectionType::periodic:
+                        if (require_uv_skeleton || require_v_skeleton)
+                          {
+                            // compute unique id for neighbor
+
+                            auto idn = index_set.uniqueIndex(outside_element);
+
+                            // Visit face if id is bigger
+                            bool visit_face = ids > idn || require_skeleton_two_sided;
+
+                            // unique vist of intersection
+                            if (visit_face)
+                              {
+                                // Bind local test space to neighbor element
+                                lfsvn.bind(outside_element);
+                                lfsvn_cache.update();
+
+                                // Notify assembler engine about binds
+                                assembler_engine.onBindLFSVInside(ig,lfsv_cache);
+
+                                // Notify assembler engine about binds
+                                assembler_engine.onBindLFSVOutside(ig,lfsv_cache,lfsvn_cache);
+
+                                // Skeleton integration
+                                assembler_engine.assembleVSkeleton(ig,lfsv_cache,lfsvn_cache);
+
+                                if(require_uv_skeleton){
+
+                                  // Bind local trial space to neighbor element
+                                  lfsun.bind(outside_element);
+                                  lfsun_cache.update();
+
+                                  // Notify assembler engine about binds
+                                  assembler_engine.onBindLFSUVInside(ig,
+                                                                      lfsu_cache,lfsv_cache);
+
+                                  // Notify assembler engine about binds
+                                  assembler_engine.onBindLFSUVOutside(ig,
+                                                                      lfsu_cache,lfsv_cache,
+                                                                      lfsun_cache,lfsvn_cache);
+
+                                  // Load coefficients of local functions
+                                  assembler_engine.loadCoefficientsLFSUOutside(lfsun_cache);
+
+                                  // Skeleton integration
+                                  assembler_engine.assembleUVSkeleton(ig,lfsu_cache,lfsv_cache,lfsun_cache,lfsvn_cache);
+
+                                  // Notify assembler engine about unbinds
+                                  assembler_engine.onUnbindLFSUVOutside(ig,
+                                                                        lfsu_cache,lfsv_cache,
+                                                                        lfsun_cache,lfsvn_cache);
+
+                                  // Notify assembler engine about unbinds
+                                  assembler_engine.onUnbindLFSUVInside(ig,
+                                                                      lfsu_cache,lfsv_cache);
+
+                                }
+
+                                // Notify assembler engine about unbinds
+                                assembler_engine.onUnbindLFSVOutside(ig,lfsv_cache,lfsvn_cache);
+
+                                // Notify assembler engine about unbinds
+                                assembler_engine.onUnbindLFSVOutside(ig,lfsv_cache,lfsvn_cache);
+
+                              }
+                          }
+                        break;
+
+                      case IntersectionType::boundary:
+                        if(require_uv_boundary || require_v_boundary )
+                          {
+                            // Notify assembler engine about binds
+                            assembler_engine.onBindLFSVInside(ig,lfsv_cache);
+
+                            // Boundary integration
+                            assembler_engine.assembleVBoundary(ig,lfsv_cache);
+
+                            if(require_uv_boundary){
+
+                              // Notify assembler engine about binds
+                              assembler_engine.onBindLFSUVInside(ig,
+                                                                  lfsu_cache,lfsv_cache);
+
+                              // Boundary integration
+                              assembler_engine.assembleUVBoundary(ig,lfsu_cache,lfsv_cache);
+
+                              // Notify assembler engine about unbinds
+                              assembler_engine.onUnbindLFSUVInside(ig,
+                                                                   lfsu_cache,lfsv_cache);
+
+                            }
+
+                            // Notify assembler engine about un binds
+                            assembler_engine.onUnbindLFSVInside(ig,lfsv_cache);
+
+                          }
+                        break;
+
+                      case IntersectionType::processor:
+                        if(require_uv_processor || require_v_processor )
+                          {
+
+                            // Processor integration
+                            assembler_engine.assembleVProcessor(ig,lfsv_cache);
+
+                            if(require_uv_processor){
+                              // Processor integration
+                              assembler_engine.assembleUVProcessor(ig,lfsu_cache,lfsv_cache);
+                            }
+                          }
+                        break;
+                      } // switch
+
+                    ++intersection_index;
+                  } // iit
+              } // do skeleton
+
+            if(require_uv_post_skeleton || require_v_post_skeleton){
+              // Volume integration
+              assembler_engine.assembleVVolumePostSkeleton(eg,lfsv_cache);
+
+              if(require_uv_post_skeleton){
+                // Volume integration
+                assembler_engine.assembleUVVolumePostSkeleton(eg,lfsu_cache,lfsv_cache);
+              }
+            }
+
+            // Notify assembler engine about unbinds
+            assembler_engine.onUnbindLFSUV(eg,lfsu_cache,lfsv_cache);
+
+            // Notify assembler engine about unbinds
+            assembler_engine.onUnbindLFSV(eg,lfsv_cache);
+
+          } // it
+
+        // Notify assembler engine that assembly is finished
+        assembler_engine.postAssembly(gfsu,gfsv);
+
+      }
+
+    private:
+
+      /* global function spaces */
+      const GFSU& gfsu;
+      const GFSV& gfsv;
+
+      typename std::conditional<
+        std::is_same<CU,EmptyTransformation>::value,
+        const CU,
+        const CU&
+        >::type cu;
+      typename std::conditional<
+        std::is_same<CV,EmptyTransformation>::value,
+        const CV,
+        const CV&
+        >::type cv;
+
+      /* local function spaces */
+      typedef LocalFunctionSpace<GFSU, TrialSpaceTag> LFSU;
+      typedef LocalFunctionSpace<GFSV, TestSpaceTag> LFSV;
+      // local function spaces in local cell
+      mutable LFSU lfsu;
+      mutable LFSV lfsv;
+      // local function spaces in neighbor
+      mutable LFSU lfsun;
+      mutable LFSV lfsvn;
+
+    };
+
+  }
+}
+#endif // DUNE_DORIE_GRIDOPERATOR_DEFAULT_ASSEMBLER_HH

dune/dorie/solver/flux_reconstruction/raviart_thomas_operator/localfunctionspace.hh

+#ifndef DUNE_DORIE_LOCAL_FUNCTION_SPACE_HH
+#define DUNE_DORIE_LOCAL_FUNCTION_SPACE_HH
+
+#include <dune/localfunctions/common/interfaceswitch.hh>
+
+#include <dune/pdelab/gridfunctionspace/localfunctionspace.hh>
+
+namespace Dune{
+namespace Dorie{
+
+/*  This is a special local function space for the raviart thomas engine. 
+ *  In particular, it is thought to work with one grid function space, that is, 
+ *  without composite or vector spaces so that we avoid anything of the type 
+ *  tree and ordering nightmare. It is simply not needed here. Indeed we don't 
+ *  even need the dof indices since the engine does a local solution and there 
+ *  is no need to map to entities. In any case, the reason of this object is 
+ *  because the user is expecting a local function space as an argument of the 
+ *  local operator (even though no many really know how a local function space 
+ *  looks like ¬_¬).
+ */
+
+  template <class GFS, class TAG=Dune::PDELab::AnySpaceTag>
+  class FakeLocalFunctionSpace
+  {
+    static_assert(not GFS::Traits::isComposite,
+      "This local function space does not work with "
+      "composited or power grid function spaces");
+
+    using DOFIndex = typename Dune::PDELab::gfs_to_lfs<GFS>::DOFIndex;
+
+  public:
+
+    using Traits = typename Dune::PDELab::LocalFunctionSpace<GFS>::Traits;
+
+    using FESwitch = Dune::FiniteElementInterfaceSwitch<
+                                  typename Traits::FiniteElementType>;
+
+    FakeLocalFunctionSpace(std::shared_ptr<const GFS> gfs) 
+      : pgfs(gfs)
+    {}
+
+    FakeLocalFunctionSpace(const GFS & gfs) 
+      : pgfs(stackobject_to_shared_ptr(gfs)) 
+    {}
+
+    FakeLocalFunctionSpace(const FakeLocalFunctionSpace & lfs) {}
+
+    typename Traits::IndexContainer::size_type 
+    localIndex(typename Traits::IndexContainer::size_type index) const
+    {
+      return index;
+    }
+
+    typename Traits::IndexContainer::size_type size () const
+    {
+      return n;
+    }
+
+    const typename Traits::DOFIndex& 
+    dofIndex(typename Traits::IndexContainer::size_type index) const
+    {
+      DUNE_THROW(Dune::NotImplemented,
+        "This local function space never constructs dof indices");
+      return Traits::DOFIndex();
+    }
+
+    // This is indeed instantiated as an intersection. Since we don't care about 
+    // the dof intex, a template is is ok to forward it to the finite element 
+    // map. In other cases one whould have to hard code the type.
+    template<class Entity>
+    void bind(Entity& e) 
+    {
+      FESwitch::setStore(pfe, pgfs->finiteElementMap().find(e));
+      n = FESwitch::basis(*pfe).size();
+    }
+
+    //! get finite element
+    const typename Traits::FiniteElementType& finiteElement () const
+    {
+      assert(pfe);
+      return *pfe;
+    }
+
+  private:
+    std::shared_ptr<const GFS> pgfs;
+    typename Traits::IndexContainer::size_type n;
+    typename FESwitch::Store pfe;
+  };
+
+} // namespace Dorie
+} // namespace Dune
+
+#endif // DUNE_DORIE_LOCAL_FUNCTION_SPACE_HH

dune/dorie/solver/flux_reconstruction/rt_projection.hh

@@ -4,6 +4,7 @@
 #include "localengeine.hh"
 #include "skeleton_fem.hh"
 #include "volume_fem.hh"
+#include "raviart_thomas_operator/assembler.hh"
 
 #include <dune/pdelab/finiteelementmap/raviartthomasfem.hh>
 
@@ -61,6 +62,8 @@ class RaviartThomasFluxReconstruction
     using SkeletonFEM = typename Dune::Dorie::SkeletonFEM<DF,RF,dim,order,gt>::type;
     using GFSVSkeleton = Dune::PDELab::GridFunctionSpace<ES,SkeletonFEM>;
     using LocalRaviartThomasEngine = Dune::Dorie::LocalRaviartThomasAssemblerEngine<LA,GFSVVolume,GFSVSkeleton>;
+    // using Assember = typename GOP::Assembler;
+    using Assember = Dune::Dorie::DefaultAssembler<GFSU,GFSW,Dune::PDELab::EmptyTransformation,Dune::PDELab::EmptyTransformation>;
 
 public:
 
@@ -99,7 +102,12 @@ public:
         MBE mbe(0);
         GOP gop(_gfsu,gfsw,local_operator,mbe);
 
-        const auto& global_assembler_gop = gop.assembler();
+        // TODO: send our DefaultAssamber to PDELab so that it can be retrived 
+        // from the grid operator:
+        // const auto& global_assembler_gop = gop.assembler();
+
+        Assember global_assembler_gop(_gfsu,gfsw);
+
         const auto& local_assembler_gop = gop.localAssembler();
 
 

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 = ""