Merge branch 'feature/interpolation' into 'master'

Feature/interpolation

See merge request !5

.gitlab-ci.yml

+---
+before_script:
+  - source ~/toolchain
+  - export CMAKE_FLAGS="-DCMAKE_C_COMPILER='$CC' -DCMAKE_CXX_COMPILER='$CXX'"
+
+dune-2.6 gcc-6-14:
+  image: registry.dune-project.org/docker/ci/dune-pdelab-deps:2.6-debian-9-gcc-6-14
+  script: duneci-standard-test
+
+dune-2.6 gcc-8-17:
+  image: registry.dune-project.org/docker/ci/dune-pdelab-deps:2.6-debian-10-gcc-8-17
+  script: duneci-standard-test
+
+dune-2.6 clang-6-17:
+  image: registry.dune-project.org/docker/ci/dune-pdelab-deps:2.6-debian-10-clang-6-libcpp-17
+  script: duneci-standard-test
+
+dune-git gcc-7-14:
+  image: registry.dune-project.org/docker/ci/dune-pdelab-deps:git-debian-10-gcc-7-14
+  script: duneci-standard-test
+
+dune-git gcc-8-17:
+  image: registry.dune-project.org/docker/ci/dune-pdelab-deps:git-debian-10-gcc-8-noassert-17
+  script: duneci-standard-test
+
+dune-git clang-6-17:
+  image: registry.dune-project.org/docker/ci/dune-pdelab-deps:git-debian-10-clang-6-libcpp-17
+  script: duneci-standard-test
\ No newline at end of file

dune.module

@@ -7,5 +7,5 @@ Module: dune-multimesh
 Version: 0.1
 Maintainer: simon.praetorius@tu-dresden.de
 #depending on
-Depends: dune-common dune-grid
-Suggests: dune-uggrid dune-alugrid dune-foamgrid dune-functions dune-geometry dune-istl
+Depends: dune-grid
+Suggests: dune-alugrid dune-foamgrid dune-functions dune-istl dune-spgrid dune-uggrid

dune/multimesh/mmentity.hh

@@ -19,20 +19,21 @@ namespace Dune
   private:
     using ctype = typename HostGrid::ctype;
 
+    /// Entity in the host grid
+    using HostGridEntity = typename HostGrid::Traits::template Codim<0>::Entity;
+
   public:
-    enum { dimension = HostGrid::dimension };
+    enum { codimension = 0 };
+    enum { dimension = HostGridEntity::dimension };
+    enum { mydimension = HostGridEntity::mydimension };
 
     /// The type of a local geometry
     using LocalGeometry = MultiMeshLocalGeometry<dimension, dimension, HostGrid>;
 
-    using EntitySeed = MultiMeshEntitySeed<0, HostGrid>;
-
-    /// Entity in the host grid
-    using HostGridEntity = typename HostGrid::Traits::template Codim<0>::Entity;
-
-    /// Containertype
-    using Super = std::vector<HostGridEntity>;
+    /// Seed representing the vector of entities in the MultiMesh
+    using EntitySeed = MultiEntitySeed<0, HostGrid>;
 
+  public:
     /// Constructor from std::vector
     using std::vector<HostGridEntity>::vector;
 
@@ -62,6 +63,39 @@ namespace Dune
         return localGeometry(source, target).global(sourceLocal);
     }
 
+    /// Return the partitionType of any entity in the multiEntity
+    // [[ expects: all partitionType are the same ]]
+    PartitionType partitionType () const
+    {
+      return max().partitionType();
+    }
+
+    /// Return the geometry of the entity with maximal level
+    typename HostGridEntity::Geometry geometry () const
+    {
+      return max().geometry();
+    }
+
+    /// Return the type of the entity with maximal level
+    GeometryType type () const
+    {
+      return max().type();
+    }
+
+    /// The entities are always regular, since irregular entities are not allowed in multimesh
+    bool isRegular() const { return true; }
+
+    /// Return the maximal level
+    int level () const
+    {
+      return max().level();
+    }
+
+    /// Return whether the entity with minimal level has boundary intersections
+    bool hasBoundaryIntersections () const
+    {
+      return min().hasBoundaryIntersections();
+    }
 
     /// \brief Return the entity seed which contains sufficient information
     /// to generate the entity again and uses as little memory as possible.
@@ -69,9 +103,15 @@ namespace Dune
      * The MultiMeshEntitySeed contains the HostGridEntitySeed and the index of
      * the grid it is extracted from.
      **/
-    EntitySeed seed (std::size_t entity_i) const
+    MultiMeshEntitySeed<0,HostGrid> seed (std::size_t idx) const
+    {
+      return {(*this)[idx], idx};
+    }
+
+    /// Return an entity seed to restore the multi-entity.
+    EntitySeed seed () const
     {
-      return {(*this)[entity_i], entity_i};
+      return {*this};
     }
 
     /// Return the entity with maximal level

dune/multimesh/mmentityseed.hh

@@ -10,16 +10,16 @@ namespace Dune
    * \ingroup MultiMesh
    *
    */
-  template <int codim, class GridImp>
+  template <int codim, class HostGrid>
   class MultiMeshEntitySeed
   {
   protected:
 
     /// Entity type of the hostgrid
-    typedef typename GridImp::HostGridType::Traits::template Codim<codim>::Entity HostEntity;
+    using HostEntity = typename HostGrid::Traits::template Codim<codim>::Entity;
 
     /// EntitySeed type of the hostgrid
-    typedef typename GridImp::HostGridType::Traits::template Codim<codim>::EntitySeed HostEntitySeed;
+    using HostEntitySeed = typename HostGrid::Traits::template Codim<codim>::EntitySeed;
 
   public:
     enum { codimension = codim };
@@ -61,6 +61,55 @@ namespace Dune
     std::size_t gridIdx_;
   };
 
+
+  template <class HG>
+  class MultiEntity;
+
+  // An entity seed for a vector of entities
+  template <int codim, class HostGrid>
+  class MultiEntitySeed
+  {
+  protected:
+
+    /// Entity type of the hostgrid
+    using HostEntity = typename HostGrid::Traits::template Codim<codim>::Entity;
+
+    /// EntitySeed type of the hostgrid
+    using HostEntitySeed = typename HostGrid::Traits::template Codim<codim>::EntitySeed;
+
+  public:
+    enum { codimension = codim };
+
+    /// Construct an empty (i.e. isValid() == false) seed.
+    MultiEntitySeed()
+    {}
+
+    /// \brief Create EntitySeed from hostgrid Entity
+    template <class MultiEntity>
+    MultiEntitySeed (const MultiEntity& multiEntity)
+    {
+      for (auto const& e : multiEntity)
+        hostEntitySeeds_.emplace_back(e.seed());
+    }
+
+    /// Get stored HostEntitySeed
+    const std::vector<HostEntitySeed>& hostEntitySeeds () const
+    {
+      return hostEntitySeeds_;
+    }
+
+    /// Check whether it is safe to create an Entity from this Seed
+    bool isValid () const
+    {
+      return hostEntitySeeds_.empty() ? false :
+        std::all_of(hostEntitySeeds_.begin(), hostEntitySeeds_.end(),
+          [](auto const& seed) { return seed.isValid(); });
+    }
+
+  private:
+    std::vector<HostEntitySeed> hostEntitySeeds_;
+  };
+
 } // end namespace Dune
 
 #endif  // DUNE_MULTIMESH_ENTITYSEED_HH

dune/multimesh/mmgridfactory.hh

@@ -8,6 +8,8 @@
 #include <vector>
 
 #include <dune/common/hybridutilities.hh>
+#include <dune/common/to_unique_ptr.hh>
+#include <dune/common/version.hh>
 #include <dune/common/std/type_traits.hh>
 #include <dune/grid/common/gridfactory.hh>
 
@@ -147,6 +149,7 @@ namespace Dune
      * Create n copies of the grid and store it in unique_pointers
      * inside the multimesh grid.
      **/
+#if DUNE_VERSION_LT(DUNE_GRID,2,7)
     virtual Grid* createGrid () override
     {
       Grid* multimesh = new Grid{};
@@ -154,6 +157,15 @@ namespace Dune
         multimesh->grids_.emplace_back(gridFactory->createGrid());
       return multimesh;
     }
+#else
+    virtual ToUniquePtr<Grid> createGrid () override
+    {
+      auto multimesh = makeToUnique<Grid>();
+      for (auto& gridFactory : gridFactories_)
+        multimesh->grids_.emplace_back(gridFactory->createGrid());
+      return std::move(multimesh);
+    }
+#endif
 
   private:
     std::vector<std::unique_ptr<GridFactory<HostGrid>>> gridFactories_;

dune/multimesh/mmgridview.hh

@@ -56,7 +56,7 @@ namespace Dune
     using Grid = typename Traits::Grid;
 
     /// Type of the corresponding GridType hosted by the MultiMesh
-    using HostGrid = typename GridImp::HostGridType;
+    using HostGrid = typename GridImp::HostGrid;
 
     using IndexSet = typename Traits::IndexSet;
     using Intersection = typename Traits::Intersection;
@@ -84,7 +84,7 @@ namespace Dune
 
     const HostGrid& grid (std::size_t i) const
     {
-      return multiMesh_->grid(i);
+      return (*multiMesh_)[i];
     }
 
     /// Obtain the level-indexSet
@@ -255,7 +255,7 @@ namespace Dune
 
     const HostGrid& grid (std::size_t i) const
     {
-      return multiMesh_[i];
+      return (*multiMesh_)[i];
     }
 
     /// Obtain the level-indexSet

dune/multimesh/mmhierarchiciterator.hh

@@ -5,10 +5,11 @@
 
 #include <numeric>
 #include <stack>
-#include <variant>
 
 #include <dune/common/iteratorfacades.hh>
+#include <dune/common/iteratorrange.hh>
 #include <dune/common/std/type_traits.hh>
+#include <dune/grid/common/gridenums.hh>
 
 namespace Dune
 {
@@ -30,7 +31,8 @@ namespace Dune
 
     struct EntityStackEntry
     {
-      template <class Entity>
+      template <class Entity,
+        std::enable_if_t<std::is_same<std::decay_t<Entity>, HostEntity>::value, int> = 0>
       explicit EntityStackEntry (Entity&& entity)
         : it(entity.hbegin(entity.level()+1))
         , end(entity.hend(entity.level()+1))
@@ -117,11 +119,17 @@ namespace Dune
       // 3. go down in tree until leaf entity
       auto child = dereference();
       for (; !(contains_(child) || child.isLeaf()); child = dereference()) {
-        assert(child.isRegular() && "No irregular elements allowed in multi-mesh traversal");
         entityStack_.emplace(child);
         assert( entityStack_.size() <= maxLevel_ );
       }
 
+      // 4. go up in tree again to the first regular entity, since
+      // irregular element can not be traversed in a multi-mesh sense
+      while (!child.isRegular() && !entityStack_.empty()) {
+        entityStack_.pop();
+        child = dereference();
+      }
+
       assert(contains_(child) && "No valid child element found in gridView");
     }
 
@@ -150,11 +158,17 @@ namespace Dune
       // 1. go down in tree until leaf entity
       auto child = dereference();
       for (; !(contains_(child) || child.isLeaf()); child = dereference()) {
-        assert(child.isRegular() && "No irregular elements allowed in multi-mesh traversal");
         entityStack_.emplace(child);
         assert( entityStack_.size() <= maxLevel_ );
       }
 
+      // 2. go up in tree again to the first regular entity, since
+      // irregular element can not be traversed in a multi-mesh sense
+      while (!child.isRegular() && !entityStack_.empty()) {
+        entityStack_.pop();
+        child = dereference();
+      }
+
       assert(contains_(child) && "No valid child element found in gridView");
     }
 
@@ -167,7 +181,7 @@ namespace Dune
   };
 
   template <class Entity, class GridView>
-  inline auto childs(const Entity& entity, const GridView& gridView)
+  inline auto childs (const Entity& entity, const GridView& gridView)
   {
     using Iterator = MultiMeshHierarchicIterator<typename GridView::Grid>;
     return IteratorRange<Iterator>{ Iterator{entity, gridView}, Iterator{true} };

dune/multimesh/mmiteratorbase.hh

@@ -6,6 +6,9 @@
 #include <numeric>
 #include <stack>
 #include <type_traits>
+#if DUNE_HAVE_CXX_VARIANT
+#include <variant>
+#endif
 
 #include <dune/common/std/type_traits.hh>
 #include <dune/grid/common/gridenums.hh>
@@ -55,8 +58,8 @@ namespace Dune
       // go to first leaf entity on all grids
       entityStacks_.reserve(multiMesh->size());
       for (std::size_t i = 0; i < multiMesh->size(); ++i) {
-        maxLevel_[i] = (*multiMesh)[i].maxLevel();
-        entityStacks_.emplace_back((*multiMesh)[i].maxLevel());
+        maxLevel_[i] = multiMesh->maxLevel(i);
+        entityStacks_.emplace_back(multiMesh->maxLevel(i));
       }
     }
 
@@ -152,13 +155,18 @@ namespace Dune
       // 3. go down in tree until leaf entity
       auto child = dereference(i);
       for (; !this->levelReached(i, child); child = dereference(i)) {
-        assert(child.isRegular() && "No irregular elements allowed in multi-mesh traversal");
         entityStack.emplace(child);
         assert(int(entityStack.size()) <= maxLevel_[i]);
       }
 
+      // 4. go up in tree again to the first regular entity, since
+      // irregular element can not be traversed in a multi-mesh sense
+      while (!child.isRegular() && !entityStack.empty()) {
+        entityStack.pop();
+        child = dereference(i);
+      }
+
       return entityStack.size();
-      // assert(contains_[i](child) && "No valid child element found in gridView");
     }
     using Super::increment;
 
@@ -167,12 +175,12 @@ namespace Dune
     {
       std::size_t size = entityStacks_[i].size();
       return std::accumulate(entityStacks_.begin(), entityStacks_.end(), true,
-        [i,size](bool allowed, auto const& entityStack) {
+        [size](bool allowed, auto const& entityStack) {
           return allowed && (entityStack.size() <= size || entityStack.finished(size));
         });
     }
 
-    // got to first leaf entity on grid i
+    // go to first entity on grid i on the desired level
     int initialIncrement (std::size_t i)
     {
       auto& entityStack = entityStacks_[i];
@@ -186,16 +194,22 @@ namespace Dune
       // 1. go down in tree until desired level is reached
       auto child = dereference(i);
       for (; !this->levelReached(i, child); child = dereference(i)) {
-        assert(child.isRegular() && "No irregular elements allowed in multi-mesh traversal");
         entityStack.emplace(child);
         assert(int(entityStack.size()) <= maxLevel_[i]);
       }
 
+      // 2. go up in tree again to the first regular entity, since
+      // irregular element can not be traversed in a multi-mesh sense
+      while (!child.isRegular() && !entityStack.empty()) {
+        entityStack.pop();
+        child = dereference(i);
+      }
+
       return entityStack.size();
-      // assert(contains_[i](child) && "No valid child element found in gridView");
     }
     using Super::initialIncrement;
 
+    // Return the current entity the hierarchic iterator or macro iterator pointing to
     HostEntity dereference (std::size_t i) const
     {
       if (entityStacks_[i].empty()) {

dune/multimesh/multimesh.hh

@@ -97,6 +97,8 @@ namespace Dune
     using Super = GridDefaultImplementation<HG::dimension, HG::dimensionworld,
       typename HG::ctype, MultiMeshFamily<HG> >;
 
+    using Self = MultiMesh;
+
   public:
     using HostGrid = HG;
 
@@ -109,6 +111,8 @@ namespace Dune
     /// The type used to store coordinates, inherited from the HostGrid
     using ctype = typename HostGrid::ctype;
 
+    enum { dimension = HostGrid::dimension };
+    enum { dimensionworld = HostGrid::dimensionworld };
 
     /// \brief Constructor, stores n copies of the hostgrid
     /**
@@ -127,7 +131,7 @@ namespace Dune
     MultiMesh () = default;
 
     /// Return the number of grids handled by this MultiMesh
-    std::size_t size() const
+    std::size_t size () const
     {
       return grids_.size();
     }
@@ -157,6 +161,11 @@ namespace Dune
         [](int level, auto const& grid) { return std::max(level, grid->maxLevel()); });
     }
 
+    /// Return maximum level in the `idx`th grid
+    int maxLevel (std::size_t idx) const
+    {
+      return grids_[idx]->maxLevel();
+    }
 
     /// Iterator to first entity of given codim on level
     template <int codim, PartitionIteratorType PiType = All_Partition>
@@ -303,12 +312,23 @@ namespace Dune
     /// Create Entity from EntitySeed
     template <class EntitySeed>
     typename Traits::template Codim<EntitySeed::codimension>::Entity
-    entity(const EntitySeed& seed) const
+    entity (const EntitySeed& seed) const
     {
       std::size_t gridIdx = this->getRealImplementation(seed).gridIndex();
       return grids_[gridIdx]->entity(this->getRealImplementation(seed).hostEntitySeed());
     }
 
+    /// Create a MultiEntity from a MultiEntitySeed
+    template <int codim>
+    MultiEntity<HostGrid> entity (const MultiEntitySeed<codim, Self>& seed) const
+    {
+      assert(seed.isValid());
+      MultiEntity<HostGrid> multiEntity;
+      auto const& hostEntitySeeds = seed.hostEntitySeeds();
+      for (std::size_t i = 0; i < hostEntitySeeds.size(); ++i)
+        multiEntity.emplace_back(grids_[i]->entity(hostEntitySeeds[i]));
+      return multiEntity;
+    }
 
     /** @name Grid Refinement Methods */
     /** @{ */
@@ -320,6 +340,24 @@ namespace Dune
         grid->globalRefine(refCount);
     }
 
+    /// Mark all entities in a multiEntitity for refinement or coarsening
+    bool mark (int refCount, const MultiEntity<HostGrid>& entities)
+    {
+      bool b = false;
+      for (std::size_t i = 0; i < grids_.size(); ++i)
+        b = grids_[i]->mark(refCount, entities[i]) || b;
+      return b;
+    }
+
+    /// Return the marks of the the entities in a multiEntity
+    std::vector<int> getMark (const MultiEntity<HostGrid>& entities) const
+    {
+      std::vector<int> marks(grids_.size());
+      for (std::size_t i = 0; i < grids_.size(); ++i)
+        marks[i] = grids_[i]->getMark(entities[i]);
+      return marks;
+    }
+
     /// Returns true, if at least one entity is marked for adaption in any of the
     /// handled grids
     bool preAdapt()
@@ -344,29 +382,39 @@ namespace Dune
 
     /** @} */
 
+    /** @name Grid Parallelization Methods */
+    /** @{ */
 
     /// Size of the overlap on the leaf level
     unsigned int overlapSize (int codim) const
     {
-      return 0u; // TODO: implement
+      return std::accumulate(grids_.begin(), grids_.end(), 0u,
+        [codim](unsigned int overlap, auto const& grid) {
+          return std::max(overlap, grid->overlapSize(codim)); });
     }
 
     /// Size of the ghost cell layer on the leaf level
     unsigned int ghostSize (int codim) const
     {
-      return 0u; // TODO: implement
+      return std::accumulate(grids_.begin(), grids_.end(), 0u,
+        [codim](unsigned int ghost, auto const& grid) {
+          return std::max(ghost, grid->ghostSize(codim)); });
     }
 
     /// Size of the overlap on a given level
     unsigned int overlapSize (int level, int codim) const
     {
-      return 0u; // TODO: implement
+      return std::accumulate(grids_.begin(), grids_.end(), 0u,
+        [level,codim](unsigned int overlap, auto const& grid) {
+          return std::max(overlap, grid->overlapSize(level,codim)); });
     }
 
     /// Size of the ghost cell layer on a given level
     unsigned int ghostSize (int level, int codim) const
     {
-      return 0u; // TODO: implement
+      return std::accumulate(grids_.begin(), grids_.end(), 0u,
+        [level,codim](unsigned int ghost, auto const& grid) {
+          return std::max(ghost, grid->ghostSize(level,codim)); });
     }
 
     /// dummy collective communication
@@ -376,6 +424,22 @@ namespace Dune
       return comm(0);
     }
 
+    /// Rebalance all grids the same way based on an average load calculation.
+    bool loadBalance ()
+    {
+      assert(false && "Needs to be implemented.");
+      return false;
+    }
+
+    /// Rebalance all grids the same way based on an average load calculation.
+    template<class DataHandle>
+    bool loadBalance (DataHandle& data)
+    {
+      assert(false && "Needs to be implemented.");
+      return false;
+    }
+
+    /** @} */
 
   public: // implementation details
 
@@ -413,7 +477,7 @@ namespace Dune
     template <class HostGrid, int codim>
     struct hasEntityIterator<MultiMesh<HostGrid>, codim>
     {
-      static const bool v = hasEntityIterator<HostGrid, codim>::v;
+      static const bool v = (codim == 0);
     };
 
     /** \brief has conforming level grids when host grid has

src/CMakeLists.txt

-add_executable("testiterator" testiterator.cc)
-target_link_dune_default_libraries("testiterator")
-if (HAVE_ALBERTA)
-  add_dune_alberta_flags(GRIDDIM 3 WORLDDIM 3 "testiterator")
-endif (HAVE_ALBERTA)
 
 add_executable("localrefinement" localrefinement.cc)
 target_link_dune_default_libraries("localrefinement")
@@ -10,42 +5,47 @@ target_link_dune_default_libraries("localrefinement")
 add_executable("uggrid" uggrid.cc)
 target_link_dune_default_libraries("uggrid")
 
-add_executable("multigridview" multigridview.cc)
-target_link_dune_default_libraries("multigridview")
+add_executable("profiling" profiling.cc)
+target_link_dune_default_libraries("profiling")
+
+add_executable("benchmark" benchmark.cc)
+target_link_dune_default_libraries("benchmark")
 
-add_executable("hierarchiciterator" hierarchiciterator.cc)
-target_link_dune_default_libraries("hierarchiciterator")
+add_executable("testiterator" testiterator.cc)
+target_link_dune_default_libraries("testiterator")
 
-add_executable("interpolate" interpolate.cc)
-target_link_dune_default_libraries("interpolate")
 if (HAVE_ALBERTA)
-  add_dune_alberta_flags("interpolate")
-endif ()
+  add_dune_alberta_flags(GRIDDIM 3 WORLDDIM 3 "testiterator")
+  add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 "benchmark")
+endif (HAVE_ALBERTA)
 
-find_package(MTL PATHS /opt/development/mtl4)
-if (MTL_FOUND)
-  set(CXX_ELEVEN_FEATURE_LIST "MOVE" "AUTO" "RANGEDFOR" "INITLIST" "STATICASSERT" "DEFAULTIMPL")
-  set(MTL_COMPILE_DEFINITIONS "")
-  foreach(feature ${CXX_ELEVEN_FEATURE_LIST})
-    list(APPEND MTL_COMPILE_DEFINITIONS "MTL_WITH_${feature}")
-  endforeach()
 
-  if (HAVE_UMFPACK OR ENABLE_SUITESPARSE OR SuiteSparse_FOUND)
-    list(APPEND MTL_COMPILE_DEFINITIONS "MTL_HAS_UMFPACK")
-  endif ()
+if (DUNE_HAVE_CXX_VARIANT)
+  add_executable("multigridview" multigridview.cc)
+  target_link_dune_default_libraries("multigridview")
 
-  set(MTL_TARGETS "")
-  list(APPEND MTL_TARGETS "phasefield" "phasefield2" "phasefield3" "phasefield4")
+  add_executable("hierarchiciterator" hierarchiciterator.cc)
+  target_link_dune_default_libraries("hierarchiciterator")
+endif (DUNE_HAVE_CXX_VARIANT)
 
-  foreach(target ${MTL_TARGETS})
+
+if (dune-functions_FOUND AND dune-alugrid_FOUND AND HAVE_ALBERTA)
+  add_executable("interpolate" interpolate.cc)
+  target_link_dune_default_libraries("interpolate")
+  add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 "interpolate")
+
+  add_executable("taylorhood" taylorhood.cc)
+  target_link_dune_default_libraries("taylorhood")
+  add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 "taylorhood")
+
+  find_package(MTL QUIET PATHS /opt/development/mtl4)
+  if (MTL_FOUND)