Merge branch 'feature/parallel_grid' into 'master'

Feature/parallel grid

See merge request !9

.gitignore

@@ -8,13 +8,18 @@
 *.parh
 *.tarh
 
+*.opts
+
 # ignore build directories
 build*/
 install*/
 output/
+html/
+latex/
 
 error.txt
 *.out
 *.out.*
 *.tmp
 *.*-swp
+*.log

dune.module

@@ -6,4 +6,4 @@ Module: dune-dec
 Version: 0.1
 Maintainer: simon.praetorius@tu-dresden.de
 Depends: dune-common dune-geometry dune-grid dune-localfunctions dune-functions
-Suggests: dune-foamgrid dune-uggrid
\ No newline at end of file
+Suggests: dune-foamgrid dune-uggrid dune-alugrid
\ No newline at end of file

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/EntityData.hpp

+#pragma once
+
+#include <array>
+#include <vector>
+
+#include <dune/grid/common/gridenums.hh>
+
+#include <dune/dec/common/ConceptsBase.hpp>
+#include <dune/dec/common/StaticLoops.hpp>
+
+#include <dune/dec/parallel/ParallelDofMapper.hpp>
+
+namespace Dec
+{
+  // storage for various data per entity. TODO: maybe activate/deactive some data with flags
+  template <int dim>
+  class EntityData
+  {
+  public:
+
+    template <class BaseGridView>
+    EntityData(BaseGridView const& gv)
+      : mapper_(constructTuple<ParallelDofMapper>(gv))
+    {
+      init(gv);
+    }
+
+    /// Returns the partitionType of the entity `e`
+    template <class Entity,
+      REQUIRES(Entity::codimension <= dim) >
+    Dune::PartitionType partitionType(Entity const& e) const
+    {
+      return partitionType_[Entity::codimension][e.index()];
+    }
+
+    /// Return the partitionType of the entity of codim `cd` with index `i`
+    template <int cd>
+    Dune::PartitionType partitionType(std::size_t i) const
+    {
+      return partitionType_[cd][i];
+    }
+
+
+    /// Returns whether entity `e` is owned by the current rank or not
+    template <class Entity,
+      REQUIRES(Entity::codimension <= dim) >
+    bool owns(Entity const& e) const
+    {
+      return std::get<Entity::codimension>(mapper_).owns(e);
+    }
+
+    /// Returns whether entity with codim `cd` and index `i` is owned by the current rank or not
+    template <int cd>
+    bool owns(std::size_t i) const
+    {
+      return std::get<cd>(mapper_).owns(i);
+    }
+
+    std::size_t memory() const
+    {
+      std::size_t mem = 0;
+
+      for (auto const& p : partitionType_)
+        mem += p.capacity()*sizeof(Dune::PartitionType) + sizeof(p);
+
+      forEach(range_<0,dim+1>, [&](auto const cd) {
+        mem += std::get<cd>(mapper_).memory();
+      });
+
+      return mem;
+    }
+
+
+  private:
+
+    template <class BaseGridView>
+    void init(BaseGridView const& gv);
+
+
+  private:
+
+    std::array<std::vector<Dune::PartitionType>, dim+1> partitionType_;
+    ParallelDofMapper mapper_;
+  };
+
+
+
+  template <int dim>
+    template <class BaseGridView>
+  void EntityData<dim>::init(BaseGridView const& gv)
+  {
+    for (int cd = 0; cd <= dim; ++cd)
+      partitionType_[cd].resize(gv.size(cd));
+
+    auto const& indexSet = gv.indexSet();
+    for (auto const& e : elements(gv)) {
+      partitionType_[0][indexSet.index(e)] = e.partitionType();
+      forEach(range_<1,dim+1>, [&](auto const cd) {
+        for (unsigned int i = 0; i < e.subEntities(cd); ++i)
+          partitionType_[cd][indexSet.subIndex(e,i,cd)] = e.template subEntity<cd>(i).partitionType();
+      });
+    }
+  }
+
+} // end namespace Dec

dune/dec/GridTransfer.hpp

@@ -83,6 +83,16 @@ namespace Dec
 //     }
 
 
+    std::size_t memory() const
+    {
+      std::size_t mem = 0;
+      mem += indexmap_.capacity()*sizeof(std::size_t) + sizeof(indexmap_);
+      mem += mat1_.memory();
+      mem += mat2_.memory();
+
+      return mem;
+    }
+
   private: // implementation of restriction strategies
 
     template <class Vector1, class Vector2>

dune/dec/Operator.hpp

@@ -2,8 +2,7 @@
 
 #include <limits>
 
-#include "utility/FunctorVector.hpp"
-#include "Mapper.hpp"
+#include <dune/dec/utility/FunctorVector.hpp>
 
 namespace Dec
 {
@@ -172,7 +171,7 @@ namespace Dec
     template <class Inserter>
     void assemble_impl(Inserter& A, float_type factor) const
     {
-      for (auto const& e : entities(gv_, Dim_<N>))
+      for (auto const& e : entities(gv_, Dim_<N>, Dune::Partitions::InteriorBorder{}))
         assembleRow(A, e, factor);
     }
 
@@ -200,7 +199,6 @@ namespace Dec
   protected: // Member variables:
 
     GridView gv_;
-
   };
 
   /** @} */

dune/dec/common/Common.hpp

@@ -95,13 +95,6 @@ namespace Dec
   template <class T>
   using Store_t = typename aux::StoreType<Decay_t<T>>::type;
 
-  /// Variadic type list
-  template <class... Ts>
-  struct Types {};
-
-  template <class... Ts>
-  using Types_t = Types<Decay_t<Ts>...>;
-
   /// Identity type wrapper, represents the type itself
   template <class T>
   struct Type
@@ -109,6 +102,23 @@ namespace Dec
     using type = T;
   };
 
+  /// Variadic type list
+  template <class... Ts>
+  struct Types
+  {
+    template <bool...> struct _bools {};
+
+    template <class T>
+    static constexpr bool contains(Type<T> = {})
+    {
+      return not std::is_same< _bools<false, std::is_same<T,Ts>::value...>,
+                               _bools<std::is_same<T,Ts>::value..., false> >::value;
+    }
+  };
+
+  template <class... Ts>
+  using Types_t = Types<Decay_t<Ts>...>;
+
   template <class T>
   using owner = T;
 

dune/dec/common/Joiner.hpp

+#pragma once
+
+#include <string>
+#include <sstream>
+
+#include <dune/dec/common/Output.hpp>
+
+namespace Dec
+{
+  struct Joiner
+  {
+    Joiner(const char* separator, std::size_t len)
+      : separator_{separator, len}
+    {}
+
+    template <class Container>
+    std::string operator()(Container const& container) const
+    {
+      if (container.empty())
+        return "";
+
+      auto it = container.begin();
+      std::stringstream ss; ss << *it;
+      for (++it; it != container.end(); ++it)
+        ss << separator_ << *it;
+
+      return ss.str();
+    }
+
+    std::string separator_;
+  };
+
+  Joiner operator ""_join(const char* separator, std::size_t len)
+  {
+    return {separator, len};
+  }
+
+} // end namespace Dec

dune/dec/common/Operations.hpp

@@ -13,6 +13,36 @@ namespace Dec
      *  @{
      **/
 
+    /// Functor returning the argument directly
+    struct id
+    {
+      template <class Iter>
+      decltype(auto) operator()(Iter&& it) const { return std::forward<Iter>(it); }
+    };
+
+    /// Functor representing a pre-increment operator
+    struct increment
+    {
+      template <class Iter, class... EndIter>
+      void operator()(Iter& it, EndIter&&...) const { ++it; }
+    };
+
+    /// Functor representing a dereferencing operator
+    struct dereferencing
+    {
+      template <class Iter>
+      auto operator()(Iter it) const { return *it; }
+    };
+
+    struct assign
+    {
+      template <class T, class S>
+      constexpr void operator()(T const& from, S& to) const
+      {
+        to = from;
+      }
+    };
+
     /// Operation that represents A+B
     struct plus
     {
@@ -23,6 +53,15 @@ namespace Dec
       }
     };
 
+    struct plus_assign
+    {
+      template <class T, class S>
+      constexpr void operator()(T const& from, S& to) const
+      {
+        to += from;
+      }
+    };
+
     /// Operation that represents A-B
     struct minus
     {
@@ -33,6 +72,15 @@ namespace Dec
       }
     };
 
+    struct minus_assign
+    {
+      template <class T, class S>
+      constexpr void operator()(T const& from, S& to) const
+      {
+        to -= from;
+      }
+    };
+
     /// Operation that represents A*B
     struct times
     {
@@ -43,6 +91,15 @@ namespace Dec
       }
     };
 
+    struct times_assign
+    {
+      template <class T, class S>
+      constexpr void operator()(T const& from, S& to) const
+      {
+        to *= from;
+      }
+    };
+
     /// Operation that represents A/B
     struct divides
     {
@@ -53,6 +110,15 @@ namespace Dec
       }
     };
 
+    struct divides_assign
+    {
+      template <class T, class S>
+      constexpr void operator()(T const& from, S& to) const
+      {
+        to /= from;
+      }
+    };
+
     /// Operation that represents A-B
     struct negate
     {

dune/dec/common/Output.hpp

@@ -5,6 +5,8 @@
 #include <sstream>
 #include <string>
 
+#include <dune/dec/common/Common.hpp>
+
 /**
  * \def DEC_NO_THROW
  * \brief The preprocessor constant sets whether to use c-asserts (if defined) or
@@ -40,6 +42,15 @@
 
 namespace Dec
 {
+  namespace tag
+  {
+    struct own_rank {};
+    struct all_ranks {};
+
+    template <class OStream> OStream& operator<<(OStream& out, own_rank) { return out; }
+    template <class OStream> OStream& operator<<(OStream& out, all_ranks) { return out; }
+  }
+
   namespace aux
   {
     template <class OStream>
@@ -67,9 +78,12 @@ namespace Dec
       int num_ranks = -1;
       MPI_Comm_rank(MPI_COMM_WORLD, &rank);
       MPI_Comm_size(MPI_COMM_WORLD, &num_ranks);
-      if (num_ranks > 1 && rank == 0) {
-        concatImpl(out, "[0] ",std::forward<Args>(args)...);
-      } else if (num_ranks == 1) {
+
+      bool print = rank == 0 || Types<Decay_t<Args>...>::contains(Type<tag::all_ranks>{});
+      if (num_ranks > 1) {
+        if (print)
+          concatImpl(out, "[",rank,"] ",std::forward<Args>(args)...);
+      } else {
         concatImpl(out, std::forward<Args>(args)...);
       }
 #else

dune/dec/halfedgegrid/Entity.hpp

@@ -55,7 +55,12 @@ namespace Dec
   public:
 
     /// Partition type of this entity
-    Dune::PartitionType partitionType() const { return Dune::InteriorEntity; }
+    Dune::PartitionType partitionType() const
+    {
+      return indexSet_->data_.template partitionType<codimension>(index());
+    }
+
+    //TODO: add level() and isLeaf()
 
     /// Returns the index of this entity
     IndexType index() const
@@ -201,6 +206,10 @@ namespace Dec
     friend auto edges(HalfEdgeEntity const& entity) { return entity.edges(); }
     friend auto faces(HalfEdgeEntity const& entity) { return entity.faces(); }
 
+    std::size_t memory() const
+    {
+      return sizeof(IndexType) + sizeof(indexSet_) + sizeof(coordinates_);
+    }
 
   // private:
   public:

dune/dec/halfedgegrid/Geometry.hpp

@@ -207,7 +207,7 @@ namespace Dec
       auto sqr_e02 = unary_dot(e02);
 
       auto D2 = std::abs(sqr_e01*sqr_e02 - sqr(dot(e01, e02)));
-      assert_msg( D2 > tol, "D2:=",D2," < tol:=",tol );
+      assert_msg( D2 > tol_, "D2:=",D2," < tol:=",tol_ );
 
       auto a1 = -sqr_e02 * dot(e01, e21) / (2*D2);
       auto a2 =  sqr_e01 * dot(e02, e21) / (2*D2);
@@ -258,6 +258,10 @@ namespace Dec
       float_type dual_vol = 0;
       for (auto const& f : faces(entity)) // faces adjacent to the edge
       {
+        // ignore contributions from ghost entities
+        if (f.partitionType() == Dune::GhostEntity)
+          continue;
+
         auto edge = center(f, Dim_<2>) - c;
         dual_vol += edge.two_norm();
       }
@@ -275,6 +279,11 @@ namespace Dec
     }
 
 
+    std::size_t memory() const
+    {
+      return sizeof(float_type) + sizeof(duneGeometry_);
+    }
+
   private: // implementation details:
 
     template <class EntityType, class Container, int cd0, int cd>
@@ -287,6 +296,10 @@ namespace Dec
       float_type dual_vol = 0;
       for (auto const& e : entities(entity, Codim_<cd-1>))
       {
+        // ignore contributions from ghost entities
+        if (e.partitionType() == Dune::GhostEntity)
+          continue;
+
         auto sub_vol = fill_dual_volume(e, corners, n+1, Codim_<cd0>, Codim_<cd-1>);
         dual_vol += sub_vol;
       }
@@ -332,7 +345,7 @@ namespace Dec
     std::vector<GlobalCoordinate> const& coordinates_; //< a reference to the vector of all grid coordinates
 
     ReferenceElement const& refElem_ = ReferenceElements::simplex(); //< the corresponding reference element
-    float_type const tol = std::numeric_limits<float_type>::epsilon();
+    float_type const tol_ = std::numeric_limits<float_type>::epsilon();
 
     mutable std::shared_ptr<DuneGeometry> duneGeometry_; //< a geometry pointer, initialized only if needed, by using \ref duneGeometry().
   };

dune/dec/halfedgegrid/GeometryCache.hpp

@@ -16,6 +16,7 @@ namespace Dec
   {
     template <class> friend class HalfEdgeGrid;
     friend class WeightedTriangulation;
+    template <class GridView, class GeometryCache> friend class GhostSynchHandle;
 
   public:
 
@@ -57,6 +58,9 @@ namespace Dec
       init_center(indexSet);
       init_volume(indexSet);
       init_dual_volume(indexSet);
+
+      // TODO: synchronize dual_volume on ghost entities
+      // TODO: sum up dual_volumes on border entities
     }
 
 
@@ -122,6 +126,21 @@ namespace Dec
     }
 
 
+    std::size_t memory() const
+    {
+      std::size_t mem = 0;
+      for (auto const& c : centers_)
+        mem += c.capacity()*sizeof(GlobalCoordinate);
+
+      for (auto const& v : volumes_)
+        mem += v.capacity()*sizeof(float_type);
+
+      for (auto const& dv : dual_volumes_)
+        mem += dv.capacity()*sizeof(float_type);
+
+      return mem + sizeof(centers_) + sizeof(volumes_) + sizeof(dual_volumes_);
+    }
+
   private:
 
     // Calculate circumcenter for all entities and store in `this->centers_`.

dune/dec/halfedgegrid/GeometrySynchronization.hpp

+#pragma once
+
+#include <dune/grid/common/datahandleif.hh>
+#include <dune/grid/common/gridenums.hh>
+#include <dune/grid/common/mapper.hh>
+#include <dune/grid/common/partitionset.hh>
+#include <dune/common/binaryfunctions.hh>
+#include <dune/common/parallel/mpihelper.hh>
+
+namespace Dec
+{
+  template <class GridView, class GeometryCache>
+  class GhostSynchHandle
+      : public Dune::CommDataHandleIF<GhostSynchHandle<GridView, GeometryCache>, double>
+  {
+    using IndexType = typename HalfEdge::IndexType;
+
+  public:
+    GhostSynchHandle(GridView const& gv, GeometryCache& cache)
+      : gv_(gv)
+      , cache_(cache)
+    {}
+
+    bool contains(int /*dim*/, int codim) const
+    {
+      return codim > 0;
+    }
+
+    bool fixedSize(int /*dim*/, int /*codim*/) const
+    {
+      return true;
+    }
+
+    template <class BaseEntity>
+    std::size_t size(BaseEntity const& /*e_*/) const
+    {
+      return 1;
+    }
+
+    template <class MessageBuffer, class BaseEntity>
+    void gather(MessageBuffer& buff, BaseEntity const& e_) const
+    {
+      auto e = gv_.entity(e_);
+      buff.write(gv_.dual_volume(e));
+//       msg("send entity<",BaseEntity::codimension,">(",e.index(),").dual_volume = ",gv_.dual_volume(e));
+    }
+
+    template <class MessageBuffer, class BaseEntity>
+    void scatter(MessageBuffer& buff, BaseEntity const& e_, std::size_t /*n*/)
+    {
+      auto e = gv_.entity(e_);
+      float_type dual_volume = 0.0;
+      buff.read(dual_volume);