Merge branch 'feature/master_grid' into 'master'

Feature/master grid

See merge request !4

dune/multimesh/mmgridfactory.hh

@@ -51,9 +51,8 @@ namespace Dune
     template <class... Args>
     explicit GridFactory (std::size_t n, Args&&... args)
     {
-      gridFactories_.reserve(n);
       for (std::size_t i = 0; i < n; ++i)
-        gridFactories_.emplace_back(std::forward<Args>(args)...);
+        gridFactories_.emplace_back(new GridFactory<HostGrid>{args...});
     }
 
     // initialize at least 1 grid
@@ -68,7 +67,7 @@ namespace Dune
     virtual void insertVertex (const GlobalCoordinate& pos) override
     {
       for (auto& gridFactory : gridFactories_)
-        gridFactory.insertVertex(pos);
+        gridFactory->insertVertex(pos);
     }
 
     /// \brief Insert an element into the coarse grid
@@ -80,7 +79,7 @@ namespace Dune
                                 const std::vector<unsigned int>& vertices) override
     {
       for (auto& gridFactory : gridFactories_)
-        gridFactory.insertElement(type, vertices);
+        gridFactory->insertElement(type, vertices);
     }
 
 
@@ -104,7 +103,7 @@ namespace Dune
       Hybrid::ifElse(Std::is_detected<HasInsertElement, GridFactory<HostGrid>>{},
       [&](auto id) {
         for (auto& gridFactory : gridFactories_)
-          id(gridFactory).insertElement(type, vertices, param);
+          id(gridFactory)->insertElement(type, vertices, param);
       });
     }
 
@@ -117,7 +116,7 @@ namespace Dune
         const std::vector<unsigned int>& vertices) override
     {
       for (auto& gridFactory : gridFactories_)
-        gridFactory.insertBoundarySegment(vertices);
+        gridFactory->insertBoundarySegment(vertices);
     }
 
 
@@ -138,7 +137,7 @@ namespace Dune
       Hybrid::ifElse(Std::is_detected<HasInsertBoundarySegment, GridFactory<HostGrid>>{},
       [&](auto id) {
         for (auto& gridFactory : gridFactories_)
-          id(gridFactory).insertBoundarySegment(vertices, boundarySegment);
+          id(gridFactory)->insertBoundarySegment(vertices, boundarySegment);
       });
     }
 
@@ -152,12 +151,12 @@ namespace Dune
     {
       Grid* multimesh = new Grid{};
       for (auto& gridFactory : gridFactories_)
-        multimesh->grids_.emplace_back(gridFactory.createGrid());
+        multimesh->grids_.emplace_back(gridFactory->createGrid());
       return multimesh;
     }
 
   private:
-    std::vector<GridFactory<HostGrid> > gridFactories_;
+    std::vector<std::unique_ptr<GridFactory<HostGrid>>> gridFactories_;
   };
 
 } // end namespace Dune

dune/multimesh/mmgridview.hh

@@ -223,13 +223,13 @@ namespace Dune
   class MultiMeshLeafGridView
   {
   public:
-    using Traits = MultiMeshLevelGridViewTraits<GridImp>;
+    using Traits = MultiMeshLeafGridViewTraits<GridImp>;
 
     /// The MultiMesh GridType
     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;
@@ -255,7 +255,7 @@ namespace Dune
 
     const HostGrid& grid (std::size_t i) const
     {
-      return multiMesh_->grid(i);
+      return multiMesh_[i];
     }
 
     /// Obtain the level-indexSet

dune/multimesh/mmhierarchiciterator.hh

@@ -22,8 +22,8 @@ namespace Dune
             typename HostGrid::template Codim<0>::Entity>
   {
   private:
-    template <int codim, PartitionIteratorType pitype, class HostGrid_>
-    friend class MultiMeshIterator;
+    template <PartitionIteratorType pitype, class HG, class D>
+    friend class MultiMeshIteratorBase;
 
     using HostEntity = typename HostGrid::template Codim<0>::Entity;
     using EntityTest = std::function<bool(HostEntity)>;

dune/multimesh/mmiteratorbase.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_MULTIMESH_ITERATOR_BASE_HH
+#define DUNE_MULTIMESH_ITERATOR_BASE_HH
+
+#include <numeric>
+#include <stack>
+#include <type_traits>
+
+#include <dune/common/std/type_traits.hh>
+#include <dune/grid/common/gridenums.hh>
+
+#include "mmentity.hh"
+#include "mmhierarchiciterator.hh"
+#include "mmiteratorinterface.hh"
+
+namespace Dune
+{
+  template <PartitionIteratorType pitype, class HostGrid, class Derived>
+  class MultiMeshIteratorBase
+      : public MultiMeshIteratorInterface<HostGrid,Derived>
+  {
+    using Super = MultiMeshIteratorInterface<HostGrid,Derived>;
+
+  protected:
+    // LevelIterator to the equivalent entity in the host grid
+    using HostGridLevelIterator =
+      typename HostGrid::template Codim<0>::template Partition<pitype>::LevelIterator;
+
+    template <class GV>
+    using IsGridView = Std::disjunction<
+      std::is_same<std::decay_t<GV>,typename HostGrid::LeafGridView>,
+      std::is_same<std::decay_t<GV>,typename HostGrid::LevelGridView> >;
+
+    using HostEntity = typename HostGrid::template Codim<0>::Entity;
+
+    using EntityTest = typename MultiMeshHierarchicIterator<HostGrid>::EntityTest;
+    using EntityStack = typename MultiMeshHierarchicIterator<HostGrid>::EntityStack;
+
+  public:
+
+    /// Constructor which creates the begin iterator
+    /**
+     * \param  multiMesh  Pointer to grid instance
+     */
+    template <class GridImp>
+    MultiMeshIteratorBase (tag::begin_iterator, const GridImp* multiMesh)
+      : maxLevel_(multiMesh->size())
+    {
+      for (std::size_t i = 0; i < multiMesh->size(); ++i) {
+        macroIterators_.push_back(multiMesh->levelGridView(i,0).template begin<0,pitype>());
+        macroEndIterators_.push_back(multiMesh->levelGridView(i,0).template end<0,pitype>());
+      }
+
+      // 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());
+      }
+    }
+
+    /// Constructor which creates the end iterator
+    /**
+     * \param  multiMesh  Pointer to grid instance
+     */
+    template <class GridImp>
+    MultiMeshIteratorBase (tag::end_iterator, const GridImp* multiMesh)
+    {
+      for (std::size_t i = 0; i < multiMesh->size(); ++i)
+        macroIterators_.push_back(multiMesh->levelGridView(i,0).template end<0,pitype>());
+    }
+
+    /// Construct an iterator from n gridViews of possibly different type
+    template <class... GridViews,
+      std::enable_if_t<Std::conjunction<IsGridView<GridViews>...>::value, int> = 0>
+    MultiMeshIteratorBase (tag::begin_iterator, GridViews const&... gridViews)
+      : maxLevel_{gridViews.grid().maxLevel()...}
+      , macroIterators_{gridViews.grid().levelGridView(0).template begin<0,pitype>()...}
+      , macroEndIterators_{gridViews.grid().levelGridView(0).template end<0,pitype>()...}
+      , entityStacks_{EntityStack{gridViews.grid().maxLevel()}...}
+    {}
+
+    template <class... GridViews,
+      std::enable_if_t<Std::conjunction<IsGridView<GridViews>...>::value, int> = 0>
+    MultiMeshIteratorBase (tag::end_iterator, GridViews const&... gridViews)
+      : macroIterators_{gridViews.grid().levelGridView(0).template end<0,pitype>()...}
+    {}
+
+#if DUNE_HAVE_CXX_VARIANT
+    /// Construct an iterator from a vector of GridViews, where each GridView can be
+    /// either a leaf- or a level-gridView, thus we need to use a variant.
+    template <class... GV>
+    MultiMeshIteratorBase (tag::begin_iterator, const std::vector<std::variant<GV...>>& gridViews)
+    {
+      entityStacks_.reserve(gridViews.size());
+      maxLevel_.reserve(gridViews.size());
+
+      for (auto const& gvs : gridViews) {
+        macroIterators_.push_back(std::visit([](auto const& gv) {
+          return gv.grid().levelGridView(0).template begin<0,pitype>(); }, gvs));
+        macroEndIterators_.push_back(std::visit([](auto const& gv) {
+          return gv.grid().levelGridView(0).template end<0,pitype>(); }, gvs));
+        maxLevel_.push_back(std::visit([](auto const& gv) {
+          return gv.grid().maxLevel(); }, gvs));
+        entityStacks_.emplace_back(std::visit([](auto const& gv) {
+          return gv.grid().maxLevel(); }, gvs));
+      }
+    }
+
+    template <class... GV>
+    MultiMeshIteratorBase (tag::end_iterator, const std::vector<std::variant<GV...>>& gridViews)
+    {
+      for (auto const& gvs : gridViews) {
+        macroIterators_.push_back(std::visit([](auto const& gv) {
+          return gv.grid().levelGridView(0).template end<0,pitype>(); }, gvs));
+      }
+    }
+#endif
+
+    /// equality
+    bool equals_impl (const MultiMeshIteratorBase& that) const
+    {
+      return macroIterators_ == that.macroIterators_;
+    }
+
+    // got to next entity in grid i
+    int increment (std::size_t i)
+    {
+      auto& entityStack = entityStacks_[i];
+      auto& macroIt = macroIterators_[i];
+      auto const& macroEnd = macroEndIterators_[i];
+
+      // 1. go up in tree or to next entity on current level until we can go down again
+      while (!entityStack.empty()) {
+        auto& top = entityStack.top();
+        ++top.it;
+        if (top.it == top.end) {
+          entityStack.pop();
+        } else {
+          break;
+        }
+      }
+
+      // 2. if entityStack is empty, go to next macroElement
+      if (entityStack.empty()) {
+        ++macroIt;
+        if (macroIt == macroEnd)
+          return 0;
+      }
+
+      // 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]);
+      }
+
+      return entityStack.size();
+      // assert(contains_[i](child) && "No valid child element found in gridView");
+    }
+    using Super::increment;
+
+    /// Return true, if all stacks with size > stack[i].size are finished
+    bool incrementAllowed (std::size_t i) const
+    {
+      std::size_t size = entityStacks_[i].size();
+      return std::accumulate(entityStacks_.begin(), entityStacks_.end(), true,
+        [i,size](bool allowed, auto const& entityStack) {
+          return allowed && (entityStack.size() <= size || entityStack.finished(size));
+        });
+    }
+
+    // got to first leaf entity on grid i
+    int initialIncrement (std::size_t i)
+    {
+      auto& entityStack = entityStacks_[i];
+      auto& macroIt = macroIterators_[i];
+      auto const& macroEnd = macroEndIterators_[i];
+
+      assert(entityStack.empty());
+      if (macroIt == macroEnd)
+        return 0;
+
+      // 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]);
+      }
+
+      return entityStack.size();
+      // assert(contains_[i](child) && "No valid child element found in gridView");
+    }
+    using Super::initialIncrement;
+
+    HostEntity dereference (std::size_t i) const
+    {
+      if (entityStacks_[i].empty()) {
+        assert(macroIterators_[i] != macroEndIterators_[i]);
+        return *macroIterators_[i];
+      } else {
+        assert(entityStacks_[i].top().it != entityStacks_[i].top().end);
+        return *entityStacks_[i].top().it;
+      }
+    }
+    using Super::dereference;
+
+    std::size_t size () const
+    {
+      return entityStacks_.size();
+    }
+
+  private:
+    std::vector<int> maxLevel_;
+
+    std::vector<HostGridLevelIterator> macroIterators_;
+    std::vector<HostGridLevelIterator> macroEndIterators_;
+    std::vector<EntityStack> entityStacks_;
+  };
+
+}  // end namespace Dune
+
+#endif // DUNE_MULTIMESH_ITERATOR_BASE_HH

dune/multimesh/mmiteratorinterface.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_MULTIMESH_ITERATOR_INTERFACE_HH
+#define DUNE_MULTIMESH_ITERATOR_INTERFACE_HH
+
+#include <dune/common/iteratorfacades.hh>
+#include <dune/grid/common/exceptions.hh>
+
+#include "mmentity.hh"
+
+namespace Dune
+{
+  namespace tag
+  {
+    struct begin_iterator {};
+    struct end_iterator {};
+  }
+
+  class DummyIterator
+  {
+  public:
+    template <class... Args>
+    DummyIterator(Args&&... args)
+    {
+      DUNE_THROW(NotImplemented, "Not implemented for entity of codim != 0!");
+    }
+
+    int operator*() const { return 0; }
+    DummyIterator& operator++() { return *this; }
+    DummyIterator operator++(int) { return *this; }
+    bool operator!=(const DummyIterator&) const { return false; }
+  };
+
+
+  /** \brief Iterator over all entities of a given codimension and level of a grid.
+   *  \ingroup MultiMesh
+   *
+   *  Implemented for codim 0 entities only
+   */
+  template <class HostGrid, class Derived>
+  class MultiMeshIteratorInterface
+      : public ForwardIteratorFacade<Derived, MultiEntity<HostGrid>, const MultiEntity<HostGrid>&>
+  {
+    using HostEntity = typename HostGrid::template Codim<0>::Entity;
+
+  public:
+
+    void initialIncrement ()
+    {
+      derived().initialIncrement_impl();
+    }
+
+    /// prefix increment
+    void increment ()
+    {
+      derived().increment_impl();
+    }
+
+    /// dereferencing
+    const MultiEntity<HostGrid>& dereference () const
+    {
+      return derived().dereference_impl();
+    }
+
+    /// equality
+    template <class Iterator>
+    bool equals (const Iterator& that) const
+    {
+      return derived().equals_impl(that);
+    }
+
+    bool levelReached (std::size_t i, const HostEntity& entity) const
+    {
+      return derived().levelReached_impl(i, entity);
+    }
+
+  private:
+    const Derived& derived() const { return static_cast<const Derived&>(*this); }
+    Derived& derived() { return static_cast<Derived&>(*this); }
+  };
+
+}  // end namespace Dune
+
+#endif // DUNE_MULTIMESH_ITERATOR_INTERFACE_HH

dune/multimesh/multigridview.hh

@@ -16,8 +16,8 @@
 #include <dune/grid/common/capabilities.hh>
 #include <dune/grid/common/gridview.hh>
 
-#include "mmiterator.hh"
 #include "multiindexset.hh"
+#include "multiiterator.hh"
 
 namespace Dune
 {
@@ -50,7 +50,7 @@ namespace Dune
       template <PartitionIteratorType pit>
       struct Partition
       {
-        using Iterator = MultiMeshIterator<cd,pit,HostGrid>;
+        using Iterator = MultiIterator<cd,pit,HostGrid>;
       };
     };
 
@@ -144,7 +144,7 @@ namespace Dune
     typename Codim<cd>::template Partition<pit>::Iterator begin () const
     {
       static_assert(cd == 0, "Implemented for codim == 0 only");
-      return MultiMeshIterator<cd,pit,HostGrid>(gridViews_);
+      return MultiIterator<cd,pit,HostGrid>(tag::begin_iterator{}, gridViews_);
     }
 
     /// Obtain end iterator for this view
@@ -152,7 +152,7 @@ namespace Dune
     typename Codim<cd>::template Partition<pit>::Iterator end () const
     {
       static_assert(cd == 0, "Implemented for codim == 0 only");
-      return MultiMeshIterator<cd,pit,HostGrid>(gridViews_, true);
+      return MultiIterator<cd,pit,HostGrid>(tag::end_iterator{}, gridViews_);
     }
 
     /// Obtain begin intersection iterator with respect to this view

dune/multimesh/multiiterator.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_MULTI_ITERATOR_HH
+#define DUNE_MULTI_ITERATOR_HH
+
+#include <type_traits>
+#include <variant>
+#include <vector>
+
+#include <dune/common/iteratorfacades.hh>
+#include <dune/common/std/type_traits.hh>
+#include <dune/grid/common/gridenums.hh>
+
+#include "mmentity.hh"
+#include "mmhierarchiciterator.hh"
+#include "mmiteratorbase.hh"
+
+namespace Dune
+{
+  template <int codim, PartitionIteratorType pitype, class HostGrid>
+  class MultiIterator : public DummyIterator
+  {
+    using DummyIterator::DummyIterator;
+  };
+
+  template <PartitionIteratorType pitype, class HostGrid>
+  class MultiIterator<0,pitype,HostGrid>
+      : public MultiMeshIteratorBase<pitype,HostGrid,MultiIterator<0,pitype,HostGrid>>
+  {
+    using Self = MultiIterator;
+    using Super = MultiMeshIteratorBase<pitype,HostGrid,Self>;
+    using HostEntity = typename HostGrid::template Codim<0>::Entity;
+    using EntityTest = typename Super::EntityTest;
+
+    template <class GV>
+    using IsGridView = Std::disjunction<
+      std::is_same<std::decay_t<GV>,typename HostGrid::LeafGridView>,
+      std::is_same<std::decay_t<GV>,typename HostGrid::LevelGridView> >;
+
+  public:
+    template <class... GridViews,
+      std::enable_if_t<Std::conjunction<IsGridView<GridViews>...>::value, int> = 0>
+    MultiIterator (tag::begin_iterator, GridViews const&... gridViews)
+      : Super{tag::begin_iterator{}, gridViews...}
+      , incrementAllowed_(sizeof...(GridViews), true)
+      , multiEntity_(sizeof...(GridViews))
+    {
+      contains_.reserve(sizeof...(GridViews));
+      Hybrid::forEach(std::forward_as_tuple(gridViews...), [this](auto const& gv) {
+        contains_.emplace_back([gv](const HostEntity& entity) { return gv.contains(entity); });
+      });
+
+      this->initialIncrement();
+    }
+
+#if DUNE_HAVE_CXX_VARIANT
+    /// Construct an iterator from a vector of GridViews, where each GridView can be
+    /// either a leaf- or a level-gridView, thus we need to use a variant.
+    template <class... GV>
+    MultiIterator (tag::begin_iterator, const std::vector<std::variant<GV...>>& gridViews)
+      : Super{tag::begin_iterator{}, gridViews}
+      , incrementAllowed_(gridViews.size(), true)
+      , multiEntity_(gridViews.size())
+    {
+      contains_.reserve(gridViews.size());
+
+      for (auto const& gvs : gridViews) {
+        contains_.emplace_back(std::visit([](auto const& gv) {
+          return EntityTest{[gv](const HostEntity& entity) { return gv.contains(entity); }}; }, gvs));
+      }
+
+      this->initialIncrement();
+    }
+#endif
+
+    // End iterator
+    template <class... Args>
+    MultiIterator (tag::end_iterator, Args&&... args)
+      : Super{tag::end_iterator{}, std::forward<Args>(args)...}
+    {}
+
+    void initialIncrement_impl ()
+    {
+      for (std::size_t i = 0; i < this->size(); ++i)
+        this->initialIncrement(i);
+    }
+
+    void increment_impl ()
+    {
+      for (std::size_t i = 0; i < this->size(); ++i)
+        incrementAllowed_[i] = this->incrementAllowed(i);
+
+      for (std::size_t i = 0; i < this->size(); ++i) {
+        if (incrementAllowed_[i])
+          this->increment(i);
+      }
+    }
+
+    const MultiEntity<HostGrid>& dereference_impl () const
+    {
+      // update entries in multiEntity that have changed
+      for (std::size_t i = 0; i <  this->size(); ++i) {
+        if (incrementAllowed_[i])
+          multiEntity_[i] = this->dereference(i);
+      }
+      return multiEntity_;
+    }
+
+    bool levelReached_impl (std::size_t i, const HostEntity& entity) const
+    {
+      return contains_[i](entity) || entity.isLeaf();
+    }
+
+  private:
+    std::vector<bool> incrementAllowed_;
+    std::vector<EntityTest> contains_;
+    mutable MultiEntity<HostGrid> multiEntity_;
+  };
+
+
+  template <class... GridViews>
+  inline auto multi_elements(GridViews const&... gridViews)
+  {
+    using GridView0 = std::tuple_element_t<0,std::tuple<GridViews...>>;
+    using Iterator = MultiIterator<0,All_Partition,typename GridView0::Grid>;
+    using Range = IteratorRange<Iterator>;
+    return Range{ Iterator{tag::begin_iterator{}, gridViews...},
+                  Iterator{tag::end_iterator{}, gridViews...} };
+  }
+
+}  // end namespace Dune
+
+#endif // DUNE_MULTI_ITERATOR_HH

dune/multimesh/multimesh.hh

@@ -24,7 +24,7 @@ namespace Dune
   template <class HostGrid>
   class MultiMesh;
 
-  template <int dim, class HostGrid>
+  template <class HostGrid>
   struct MultiMeshFamily
   {
     struct Traits : HostGrid::GridFamily::Traits
@@ -46,11 +46,11 @@ namespace Dune
         {
           /// The type of the iterator over the level entities of this codim on this partition.
           // using LevelIterator = Dune::EntityIterator<cd, const Grid, MultiMeshIterator<cd, pitype, const Grid> >;
-          using LevelIterator = MultiMeshIterator<cd, pitype, HostGrid>;
+          using LevelIterator = MultiMeshLevelIterator<cd, pitype, HostGrid>;
 
           /// The type of the iterator over the leaf entities of this codim on this partition.
           // using LeafIterator = Dune::EntityIterator<cd, const Grid, MultiMeshIterator<cd, pitype, const Grid> >;
-          using LeafIterator = MultiMeshIterator<cd, pitype, HostGrid>;
+          using LeafIterator = MultiMeshLeafIterator<cd, pitype, HostGrid>;
         };
 
         /// The type of the iterator over all leaf entities of this codim.
@@ -78,14 +78,11 @@ namespace Dune
    *
    * \tparam  HostGrid  The host grid type wrapped by the MultiMesh
    */
-  template <class HostGrid>
+  template <class HG>
   class MultiMesh
-      : public GridDefaultImplementation<HostGrid::dimension, HostGrid::dimensionworld,
-          typename HostGrid::ctype, MultiMeshFamily<HostGrid::dimension, HostGrid> >
+      : public GridDefaultImplementation<HG::dimension, HG::dimensionworld,
+          typename HG::ctype, MultiMeshFamily<HG> >
   {
-    template <int codim, PartitionIteratorType pitype, class GridImp>
-    friend class MultiMeshIterator;
-
     template <class GridImp>
     friend class MultiMeshLevelGridView;