partition-type and partition-iterators added. A mapper to map local indices to...

partition-type and partition-iterators added. A mapper to map local indices to continuouse range added

dune/dec/Dec.cpp

@@ -22,27 +22,14 @@
 
 namespace Dec {
 
-#ifdef DEC_HAS_PETSC
-decpde::decpde(int& argc, char**& argv)
-{
-  PetscInitialize(&argc, &argv, NULL, NULL);
-  std::srand(std::time(0));
-}
-#else
 decpde::decpde(int& argc, char**& argv)
+  : mpihelper_(Dune::MPIHelper::instance(argc, argv))
 {
-  Dune::MPIHelper::instance(argc, argv);
   std::srand(std::time(0));
 }
-#endif
 
 decpde::~decpde()
 {
-#ifdef DEC_HAS_PETSC
-  PetscFinalize();
-#else
-//   MPI_Finalize();
-#endif
 }
 
 } // end namespace Dec

dune/dec/Dec.hpp

@@ -11,6 +11,8 @@
 #define DEC_DOW 3
 #endif
 
+#include <dune/common/parallel/mpihelper.hh>
+
 namespace Dec
 {
   /// type of floating point values
@@ -30,6 +32,8 @@ namespace Dec
     /// Destructor. Calls \ref PetscFinalize in case that PETSc backend is enabled.
     ~decpde();
 
+    Dune::MPIHelper& mpihelper_;
+
   public:
     /// Create an singleton instance of the class \ref decpde
     static decpde& init(int& argc, char**& argv)
@@ -37,5 +41,10 @@ namespace Dec
       static decpde instance(argc, argv);
       return instance;
     }
+
+    Dune::MPIHelper& mpihelper()
+    {
+      return mpihelper_;
+    }
   };
 } // end namespace Dec

dune/dec/PartitionTypeMapper.hpp

+#pragma once
+
+#include <vector>
+
+#include <dune/grid/common/gridenums.hh>
+#include <dune/grid/common/mapper.hh>
+#include <dune/grid/common/partitionset.hh>
+
+#include <dune/dec/DecGrid.hpp>
+
+namespace Dec
+{
+  template <class GridView, int codim, Dune::PartitionIteratorType pit>
+  class PartitionTypeMapper
+      : public Dune::Mapper<typename GridView::Grid, PartitionTypeMapper<GridView, codim, pit>, typename HalfEdge::IndexType>
+  {
+  public:
+
+    using IndexType = typename HalfEdge::IndexType;
+
+  public:
+
+    /// Constructor, stores a copy of `gv`
+    PartitionTypeMapper(GridView const& gv)
+      : gridView_(gv)
+    {
+      update();
+    }
+
+
+    /// Map entity to array index
+    template <class Entity>
+    IndexType index(Entity const& e) const
+    {
+      assert( Entity::codimension == codim );
+      assert_msg( Dune::partitionSet<pit>().contains(e.partitionType()), "PartitionType ",e.partitionType()," is not in ",pit );
+
+      auto const& indexSet = gridView_.indexSet();
+      IndexType idx = indexSet.index(e);
+
+      return indexMap_[idx];
+    }
+
+
+    /// Return total number of entities in the entity set managed by the mapper.
+    int size() const
+    {
+      return size_;
+    }
+
+
+    /// Returns true if the entity is contained in the index set and at the same time the array index is returned.
+    template <class Entity>
+    bool contains(Entity const& e, int& result) const
+    {
+      if (Entity::codimension != codim)
+        return false;
+      if (!Dune::partitionSet<pit>().contains(e.partitionType()))
+        return false;
+
+      result = int(index(e));
+      return true;
+    }
+
+
+    /// Reinitialize mapper after grid has been modified.
+    void update()
+    {
+      auto const& indexSet = gridView_.indexSet();
+      indexMap_.resize(gridView_.size(codim));
+
+      IndexType idx = 0;
+      for (auto const& e : entities(gridView_, Codim_<codim>, Dune::partitionSet<pit>()))
+        indexMap_[indexSet.index(e)] = idx++;
+
+      size_ = int(idx);
+    }
+
+
+  private:
+
+    GridView gridView_;
+
+    // maps indices of PartitionIteratorType to continuouse range
+    std::vector<IndexType> indexMap_;
+
+    // number of indices in continuouse range
+    int size_;
+  };
+
+} // end namespace Dec

dune/dec/halfedgegrid/Entity.hpp

@@ -55,7 +55,7 @@ namespace Dec
   public:
 
     /// Partition type of this entity
-    Dune::PartitionType partitionType() const { return Dune::InteriorEntity; }
+    Dune::PartitionType partitionType() const { return halfEdge().partitionType_; }
 
     /// Returns the index of this entity
     IndexType index() const

dune/dec/halfedgegrid/GridView.hpp

@@ -30,6 +30,7 @@ namespace Dec
     {
       using Entity = HalfEdgeEntity<cd, dimension, Grid>;
 
+      template <Dune::PartitionIteratorType pit = Dune::All_Partition>
       struct Iterator
           : public std::iterator<std::random_access_iterator_tag, Entity>
       {
@@ -37,57 +38,60 @@ namespace Dec
         using difference_type = std::ptrdiff_t;
         using Reference = Entity;
 
-        Iterator(HalfEdgeGridView const& gv, IndexType index)
+        template <Dune::PartitionIteratorType p>
+        using pit_t = std::integral_constant<Dune::PartitionIteratorType, p>;
+
+        Iterator(HalfEdgeGridView const& gv, IndexType index, IndexType end)
           : gv_(gv)
           , index_(index)
-        {}
+          , end_(end)
+        {
+          if (pit != Dune::All_Partition) {
+            while (index_ != end_ && !Dune::partitionSet<pit>().contains(halfEdge().partitionType_))
+              ++index_;
+          }
+        }
 
-        Iterator& operator++()    { ++index_; return *this; }
+        Iterator& operator++()    { increment(pit_t<pit>{}); return *this; }
         Iterator  operator++(int) { Iterator tmp(*this); ++(*this); return tmp; }
-        Iterator& operator--()    { --index_; return *this; }
-        Iterator  operator--(int) { Iterator tmp(*this); --(*this); return tmp; }
-
-        Iterator& operator+=(difference_type rhs) { index_ += rhs; return *this; }
-        Iterator& operator-=(difference_type rhs) { index_ -= rhs; return *this; }
-
-        difference_type operator-(Iterator const& rhs) const { return index_ - rhs.index_; }
-
-        Iterator operator+(difference_type rhs) const { return {index_ + rhs}; }
-        Iterator operator-(difference_type rhs) const { return {index_ - rhs}; }
-
-        friend Iterator operator+(difference_type lhs, const Iterator& rhs) { return {lhs + rhs.index_}; }
-        friend Iterator operator-(difference_type lhs, const Iterator& rhs) { return {lhs - rhs.index_}; }
 
         bool operator==(Iterator const& other) const { return index_ == other.index_; }
         bool operator!=(Iterator const& other) const { return !(*this == other); }
-        bool operator> (Iterator const& other) const { return index_ > other.index_; }
-        bool operator< (Iterator const& other) const { return index_ < other.index_; }
-        bool operator>=(Iterator const& other) const { return index_ >= other.index_; }
-        bool operator<=(Iterator const& other) const { return index_ <= other.index_; }
 
         std::shared_ptr<Entity> operator->() const
         {
-          return std::make_shared<Entity>(gv_.indexSet(), gv_.coordinates(), halfEdge());
+          return std::make_shared<Entity>(gv_.indexSet(), gv_.coordinates(), he());
+        }
+        Entity operator*() const { return {gv_.indexSet(), gv_.coordinates(), he()}; }
+
+      private:
+
+        void increment(pit_t<Dune::All_Partition>)
+        {
+          ++index_;
         }
-        Entity operator*() const { return {gv_.indexSet(), gv_.coordinates(), halfEdge()}; }
-        Entity operator[](difference_type idx)
+
+        template <Dune::PartitionIteratorType p>
+        void increment(pit_t<p>)
         {
-          index_ += idx;
-          return {gv_.indexSet(), gv_.coordinates(), halfEdge()};
+          ++index_;
+          while (index_ != end_ && !Dune::partitionSet<p>().contains(halfEdge().partitionType_))
+            ++index_;
         }
 
-      private:
-        IndexType halfEdge() const { return gv_.indexSet().template halfEdgeOfIndex<cd>(index_); }
+        IndexType he() const { return gv_.indexSet().template halfEdgeOfIndex<cd>(index_); }
+        HalfEdge const& halfEdge() const { return gv_.indexSet().halfEdge(he()); }
 
       private:
         HalfEdgeGridView gv_;
         IndexType index_;
+        IndexType end_;
       };
 
-      template <Dune::PartitionIteratorType /*pit*/>
+      template <Dune::PartitionIteratorType pit>
       struct Partition
       {
-        using Iterator = Codim::Iterator;
+        using Iterator = Codim::Iterator<pit>;
       };
     };
 
@@ -127,18 +131,30 @@ namespace Dec
       return grid_.size(level_, type);
     }
 
+    //// Return size of the overlap region for a given codim on the grid view.
+    int overlapSize(int codim) const
+    {
+      return grid_.grid_.levelGridView(level_).overlapSize(codim);
+    }
+
+    /// Return size of the ghost region for a given codim on the grid view.
+    int ghostSize(int codim) const
+    {
+      return grid_.grid_.levelGridView(level_).ghostSize(codim);
+    }
+
     /// Obtain begin iterator for this view.
-    template <int cd, Dune::PartitionIteratorType pitype = Dune::All_Partition>
-    typename Codim<cd>::Iterator begin() const
+    template <int cd, Dune::PartitionIteratorType pit = Dune::All_Partition>
+    typename Codim<cd>::template Partition<pit>::Iterator begin() const
     {
-      return {*this, 0u};
+      return {*this, 0u, size(cd)};
     }
 
     /// Obtain end iterator for this view.
-    template <int cd, Dune::PartitionIteratorType pitype = Dune::All_Partition>
-    typename Codim<cd>::Iterator end() const
+    template <int cd, Dune::PartitionIteratorType pit = Dune::All_Partition>
+    typename Codim<cd>::template Partition<pit>::Iterator end() const
     {
-      return {*this, size(cd)};
+      return {*this, size(cd), size(cd)};
     }
 
     int level() const { return level_; }

dune/dec/halfedgegrid/HalfEdge.hpp

@@ -74,6 +74,8 @@ namespace Dec
     IndexType vertex_ = invalid; // target of the half edge
     IndexType edge_ = invalid;   // edge that the half edge is part of
     IndexType face_ = invalid;   // left element of the half edge
+
+    Dune::PartitionType partitionType_ = Dune::InteriorEntity;
   };
 
 } // end namespace Dec

dune/dec/halfedgegrid/IndexSet.hpp

@@ -358,6 +358,8 @@ namespace Dec
         h.edge_ = edge;
         h.face_ = face;
 
+        h.partitionType_ = elem.template subEntity<1>(i).partitionType();
+
         // set opposite half_edge, if edge is assigned an half-edge already.
         if (edges_[edge] != invalid) {
           HalfEdge& h_opp = half_edges_[edges_[edge]];

examples/CMakeLists.txt

@@ -5,6 +5,8 @@ set(DEC_DOW 3)
 
 add_compile_options(-Wall -Wextra -Wno-unused-parameter -Wno-deprecated-declarations)
 
+add_dec_executable(domain_decomposition ALBERTA domain_decomposition.cpp)
+target_compile_definitions(domain_decomposition PRIVATE DEC_HAS_MPI=1)
 add_dec_executable(geometry ALBERTA geometry.cpp)
 add_dec_executable(laplace ALBERTA laplace.cpp)
 add_dec_executable(laplace_operator ALBERTA laplace_operator.cpp)
@@ -15,7 +17,7 @@ add_dec_executable(transfer ALBERTA transfer.cpp)
 add_dec_executable(vecellipt ALBERTA vecellipt.cpp)
 # add_dec_executable(weighted_triangulation ALBERTA weighted_triangulation.cpp)
 add_dec_executable(weighted_triangulation2 ALBERTA weighted_triangulation2.cpp)
-add_dependencies(examples geometry laplace laplace_operator simple_grid heat weighted_triangulation2 orientation transfer)
+add_dependencies(examples geometry laplace laplace_operator simple_grid heat weighted_triangulation2 orientation transfer domain_decomposition)
 
 add_dec_executable(helmholtz ALBERTA helmholtz.cpp)
 add_dec_executable(consistency ALBERTA consistency.cpp)

examples/domain_decomposition.cpp

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/dec/DecGrid.hpp>
+#include <dune/dec/PartitionTypeMapper.hpp>
+#include <dune/dec/common/Output.hpp>
+
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/grid/utility/parmetisgridpartitioner.hh>
+#include <dune/grid/uggrid.hh>
+
+
+using namespace Dec;
+
+
+int main(int argc, char** argv)
+{
+  auto& mpihelper = Dune::MPIHelper::instance(argc, argv);
+
+  // Create ug grid from structured grid
+  std::array<unsigned int, 2> n = {{3, 2}};
+  Dune::FieldVector<double, 2> lower = {{0, 0}};
+  Dune::FieldVector<double, 2> upper = {{1, 1}};
+
+  using GridBase = Dune::UGGrid<2>;
+  std::shared_ptr<GridBase> gridBase = Dune::StructuredGridFactory<GridBase>::createSimplexGrid(lower, upper, n);
+
+  // Create initial partitioning using ParMETIS
+  std::vector<unsigned> part(Dune::ParMetisGridPartitioner<GridBase::LeafGridView>::partition(gridBase->leafGridView(), mpihelper));
+
+  // Transfer partitioning from ParMETIS to our grid
+  gridBase->loadBalance(part, 0);
+
+  using Grid = DecGrid<GridBase>;
+  Grid grid(*gridBase);
+
+  using GridView = Grid::LeafGridView;
+  const GridView gv = grid.leafGridView();
+
+  msg("All:");
+  auto const& indexSet = gv.indexSet();
+  for (auto const& e : elements(gv)) {
+    msg("element ", indexSet.index(e));
+  }
+
+  msg("----------------------------------");
+
+  msg("Interior:");
+  for (auto const& e : elements(gv, Dune::Partitions::Interior{})) {
+    msg("element ", indexSet.index(e));
+  }
+
+  msg("Ghost:");
+  for (auto const& e : elements(gv, Dune::Partitions::Ghost{})) {
+    msg("element ", indexSet.index(e));
+  }
+
+  msg("----------------------------------");
+
+  msg("All edges:");
+  for (auto const& e : edges(gv)) {
+    msg("element ", indexSet.index(e));
+  }
+
+  msg("----------------------------------");
+
+  msg("Interior edges:");
+  PartitionTypeMapper<GridView,1,Dune::Interior_Partition> mapper0(gv);
+  for (auto const& e : edges(gv, Dune::Partitions::Interior{})) {
+    msg("element ", mapper0.index(e));
+  }
+
+  msg("InteriorBorder edges:");
+  PartitionTypeMapper<GridView,1,Dune::InteriorBorder_Partition> mapper1(gv);
+  for (auto const& e : edges(gv, Dune::Partitions::InteriorBorder{})) {
+    msg("element ", mapper1.index(e));
+  }
+
+  msg("Ghost edges:");
+  PartitionTypeMapper<GridView,1,Dune::Ghost_Partition> mapper2(gv);
+  for (auto const& e : edges(gv, Dune::Partitions::Ghost{})) {
+    msg("element ", mapper2.index(e));
+  }
+
+  // NOTE: index is not continuouse any more!
+
+}

external/json11/json11.cpp

@@ -227,9 +227,9 @@ struct Statics {
     const std::shared_ptr<JsonValue> null = make_shared<JsonNull>();
     const std::shared_ptr<JsonValue> t = make_shared<JsonBoolean>(true);
     const std::shared_ptr<JsonValue> f = make_shared<JsonBoolean>(false);
-    const string empty_string;
-    const vector<Json> empty_vector;
-    const map<string, Json> empty_map;
+    const string empty_string = "";
+    const vector<Json> empty_vector = {};
+    const map<string, Json> empty_map = {};
     Statics() {}
 };