Merge branch 'multientity' into 'master'

Multientity

See merge request !1

dune/multimesh/CMakeLists.txt

 #install headers
 install(FILES
+  mmentity.hh
   mmentityseed.hh
+  mmgeometry.hh
   mmgridfactory.hh
   mmgridview.hh
   mmiterator.hh

dune/multimesh/mmentity.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_MULTIMESH_ENTITY_HH
+#define DUNE_MULTIMESH_ENTITY_HH
+
+#include <vector>
+
+#include <dune/grid/common/geometry.hh>
+
+#include "mmentityseed.hh"
+#include "mmgeometry.hh"
+
+namespace Dune
+{
+  template <class GridImp>
+  class MultiEntity
+      : public std::vector<typename GridImp::HostGridType::Traits::template Codim<0>::Entity>
+  {
+  private:
+    using ctype = typename GridImp::ctype;
+
+  public:
+    enum { dimension = GridImp::dimension };
+
+    /// The type of a local geometry
+    using LocalGeometry = Dune::Geometry<dimension, dimension, const GridImp, MultiMeshLocalGeometry>;
+
+    using EntitySeed = typename GridImp::Traits::template Codim<0>::EntitySeed;
+
+    /// Entity in the host grid
+    using HostGridEntity = typename GridImp::HostGridType::Traits::template Codim<0>::Entity;
+
+    /// Containertype
+    using Super = std::vector<HostGridEntity>;
+
+    /// Constructor from std::vector
+    using Super::vector;
+
+    /// Return a local geometry of source in target
+    static LocalGeometry localGeometry (HostGridEntity const& source, HostGridEntity const& target)
+    {
+      using LocalGeometryImp = MultiMeshLocalGeometry<dimension, dimension, const GridImp>;
+      return LocalGeometry{LocalGeometryImp{source, target}};
+    }
+
+    /// Return a local geometry of source_i'th entity in target_t'th entity
+    LocalGeometry localGeometry (std::size_t source_i, std::size_t target_i) const
+    {
+      using LocalGeometryImp = MultiMeshLocalGeometry<dimension, dimension, const GridImp>;
+      return LocalGeometry{LocalGeometryImp{(*this)[source_i], (*this)[target_i]}};
+    }
+
+    /// \brief Return the entity seed which contains sufficient information
+    /// to generate the entity again and uses as little memory as possible.
+    /**
+     * The MultiMeshEntitySeed contains the HostGridEntitySeed and the index of
+     * the grid it is extracted from.
+     **/
+    EntitySeed seed (std::size_t entity_i) const
+    {
+      using EntitySeedImp = MultiMeshEntitySeed<0, const GridImp>;
+      return EntitySeedImp{(*this)[entity_i], entity_i};
+    }
+
+    /// Return the entity with maximal level
+    HostGridEntity const& max() const
+    {
+      // TODO: cache max-element
+      auto it_max = std::max_element(this->begin(), this->end(),
+        [](auto const& e1, auto const& e2) { return e1.level() < e2.level(); });
+      return *it_max;
+    }
+
+    /// Return the entity with minimal level
+    HostGridEntity const& min() const
+    {
+      // TODO: cache min-element
+      auto it_min = std::min_element(this->begin(), this->end(),
+        [](auto const& e1, auto const& e2) { return e1.level() < e2.level(); });
+      return *it_min;
+    }
+  };
+
+}  // end namespace Dune
+
+#endif // DUNE_MULTIMESH_ENTITY_HH

dune/multimesh/mmgeometry.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_MULTIMESH_LOCALGEOMETRY_HH
+#define DUNE_MULTIMESH_LOCALGEOMETRY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/typetraits.hh>
+#include <dune/grid/common/geometry.hh>
+
+namespace Dune
+{
+  template <int mydim, int coorddim, class GridImp>
+  class MultiMeshLocalGeometry
+      : public GeometryDefaultImplementation<mydim, coorddim, GridImp, MultiMeshLocalGeometry>
+  {
+  private:
+    using ctype = typename GridImp::ctype;
+
+  public:
+    // The codimension of this entitypointer wrt the host grid
+    enum { codimension = GridImp::HostGridType::dimension - mydim };
+    enum { dimensionworld = GridImp::HostGridType::dimensionworld };
+
+    // select appropriate hostgrid geometry via typeswitch
+    using HostLocalGeometry = typename GridImp::HostGridType::Traits::template Codim<codimension>::LocalGeometry;
+
+    /// entity in the host grid
+    using HostGridEntity = typename GridImp::HostGridType::Traits::template Codim<codimension>::Entity;
+
+    //! type of jacobian transposed
+    using JacobianInverseTransposed = typename HostLocalGeometry::JacobianInverseTransposed;
+    using JacobianTransposed = typename HostLocalGeometry::JacobianTransposed;
+
+
+    /// Constructor from host entities
+    MultiMeshLocalGeometry (const HostGridEntity& sourceEntity, const HostGridEntity& targetEntity)
+    {
+      hostLocalGeometries_.reserve(sourceEntity.level() - targetEntity.level());
+      auto entity = sourceEntity;
+      for (int l = sourceEntity.level(); l > targetEntity.level(); --l) {
+        assert( entity.hasFather() );
+        hostLocalGeometries_.push_back(entity.geometryInFather());
+        if (l > targetEntity.level()+1)
+          entity = entity.father();
+      }
+    }
+
+
+    /// Return the element type identifier
+    GeometryType type () const
+    {
+      return hostLocalGeometries_[0].type();
+    }
+
+    /// Return wether we have an affine mapping
+    bool affine () const
+    {
+      return std::all_of(hostLocalGeometries_.begin(), hostLocalGeometries_.end(),
+        [](auto const& localGeo) { return localGeo.affine(); });
+    }
+
+    /// Returns true if the point is in the current element
+    bool checkInside (const FieldVector<ctype, mydim>& local) const
+    {
+      return hostLocalGeometries_[0].checkInside(local);
+    }
+
+    /// Return the number of corners of this element. Corners are numbered 0...n-1
+    int corners () const
+    {
+      return hostLocalGeometries_[0].corners();
+    }
+
+    /// Access to coordinates of corners. Index is the number of the corner
+    const FieldVector<ctype, coorddim> corner (int i) const
+    {
+      FieldVector<ctype, coorddim> c = hostLocalGeometries_[0].corner(i);
+      for (std::size_t j = 1; j < hostLocalGeometries_.size(); ++j)
+        c = hostLocalGeometries_[j].global(c);
+      return c;
+    }
+
+
+    /// Maps a local coordinate within reference element to global coordinate in element
+    FieldVector<ctype, coorddim> global (const FieldVector<ctype, mydim>& local) const
+    {
+      auto c = hostLocalGeometries_[0].global(local);
+      for (std::size_t j = 1; j < hostLocalGeometries_.size(); ++j)
+        c = hostLocalGeometries_[j].global(c);
+      return c;
+    }
+
+    /// Maps a global coordinate within the element to a local coordinate in its reference element
+    FieldVector<ctype, mydim> local (const FieldVector<ctype, coorddim>& global) const
+    {
+      auto c = hostLocalGeometries_.back().local(global);
+      std::size_t n = hostLocalGeometries_.size()-1;
+      for (std::size_t j = 1; j < hostLocalGeometries_.size(); ++j)
+        c = hostLocalGeometries_[n-j].local(c);
+      return c;
+    }
+
+    /// Return the factor appearing in the integral transformation formula
+    ctype integrationElement (const FieldVector<ctype, mydim>& local) const
+    {
+      ctype dx = 1;
+      auto c = local;
+      for (std::size_t j = 0; j < hostLocalGeometries_.size(); ++j) {
+        dx *= hostLocalGeometries_[j].integrationElement(c);
+        if (j < hostLocalGeometries_.size()-1)
+          c = hostLocalGeometries_[j].global(local);
+      }
+      return dx;
+    }
+
+    /// The Jacobian matrix of the mapping from the reference element to this element
+    JacobianInverseTransposed jacobianInverseTransposed (const FieldVector<ctype, mydim>& local) const
+    {
+      auto J = hostLocalGeometries_[0].jacobianInverseTransposed(local);
+      auto c = hostLocalGeometries_[0].global(local);
+      for (std::size_t j = 1; j < hostLocalGeometries_.size(); ++j) {
+        J.leftmultiply(hostLocalGeometries_[j].jacobianInverseTransposed(c));
+        if (j < hostLocalGeometries_.size()-1)
+          c = hostLocalGeometries_[j].global(local);
+      }
+      return J;
+    }
+
+    /// Return the transposed of the Jacobian
+    JacobianTransposed jacobianTransposed (const FieldVector<ctype, mydim>& local) const
+    {
+      auto J = hostLocalGeometries_[0].jacobianTransposed(local);
+      auto c = hostLocalGeometries_[0].global(local);
+      for (std::size_t j = 1; j < hostLocalGeometries_.size(); ++j) {
+        J.rightmultiply(hostLocalGeometries_[j].jacobianTransposed(c));
+        if (j < hostLocalGeometries_.size()-1)
+          c = hostLocalGeometries_[j].global(local);
+      }
+      return J;
+    }
+
+  private:
+    std::vector<HostLocalGeometry> hostLocalGeometries_;
+  };
+
+}  // end namespace Dune
+
+#endif // DUNE_MULTIMESH_LOCALGEOMETRY_HH

dune/multimesh/mmiterator.hh

@@ -10,9 +10,7 @@
 #include <dune/grid/common/exceptions.hh>
 #include <dune/grid/common/gridenums.hh>
 
-/** \file
- * \brief The MultiMeshIterator class
- */
+#include "mmentity.hh"
 
 namespace Dune
 {
@@ -40,8 +38,8 @@ namespace Dune
   template <PartitionIteratorType pitype, class GridImp>
   class MultiMeshIterator<0, pitype, GridImp>
       : public ForwardIteratorFacade<MultiMeshIterator<0,pitype,GridImp>,
-            std::vector<typename GridImp::Traits::template Codim<0>::Entity>,
-            std::vector<typename GridImp::Traits::template Codim<0>::Entity> >
+            MultiEntity<GridImp>,
+            MultiEntity<GridImp> const&>
   {
   private:
     // LevelIterator to the equivalent entity in the host grid
@@ -148,7 +146,7 @@ namespace Dune
     }
 
     /// dereferencing
-    std::vector<HostEntity> dereference () const
+    MultiEntity<GridImp> const& dereference () const
     {
       // update entries in multiEntity that have changed
       for (std::size_t i = 0; i < entityStacks_.size(); ++i) {
@@ -253,7 +251,7 @@ namespace Dune
     std::vector<bool> incrementAllowed_;
     int level_ = -1;
 
-    mutable std::vector<HostEntity> multiEntity_;
+    mutable MultiEntity<GridImp> multiEntity_;
   };
 
 }  // end namespace Dune

dune/multimesh/multimesh.hh

@@ -13,6 +13,7 @@
 #include <dune/grid/common/gridfactory.hh>
 
 // The components of the MultiMesh interface
+#include "mmentity.hh"
 #include "mmentityseed.hh"
 #include "mmgridview.hh"
 #include "mmiterator.hh"
@@ -82,22 +83,25 @@ namespace Dune
       : public GridDefaultImplementation<HostGrid::dimension, HostGrid::dimensionworld,
           typename HostGrid::ctype, MultiMeshFamily<HostGrid::dimension, HostGrid> >
   {
-    template <int codim, PartitionIteratorType pitype, class GridImp_>
+    template <int codim, PartitionIteratorType pitype, class GridImp>
     friend class MultiMeshIterator;
 
-    template <class GridImp_>
+    template <class GridImp>
     friend class MultiMeshLevelGridView;
 
-    template <class GridImp_>
+    template <class GridImp>
     friend class MultiMeshLeafGridView;
 
+    template <class GridImp>
+    friend class MultiEntity;
+
     friend class GridFactory<MultiMesh<HostGrid> >;
 
-    using HostGridType = HostGrid;
     using Super = GridDefaultImplementation<HostGrid::dimension, HostGrid::dimensionworld,
       typename HostGrid::ctype, MultiMeshFamily<HostGrid::dimension, HostGrid> >;
 
   public:
+    using HostGridType = HostGrid;
 
     /// Type of the used GridFamily for this grid
     using GridFamily = MultiMeshFamily<HostGrid::dimension,HostGrid>;
@@ -222,7 +226,7 @@ namespace Dune
     /// Access to the GlobalIdSet
     const typename HostGrid::GlobalIdSet& globalIdSet (std::size_t idx) const
     {
-      return grids_[idx].globalIdSet();
+      return grids_[idx]->globalIdSet();
     }
 
     /// Return the globalIdSet for grid 0
@@ -236,7 +240,7 @@ namespace Dune
     /// Access to the LocalIdSet
     const typename HostGrid::LocalIdSet& localIdSet (std::size_t idx) const
     {
-      return grids_[idx].localIdSet();
+      return grids_[idx]->localIdSet();
     }
 
     /// Return the localIdSet for grid 0
@@ -253,7 +257,7 @@ namespace Dune
       if (level < 0 || level > grids_[idx].maxLevel()) {
         DUNE_THROW(GridError, "levelIndexSet of nonexisting level " << level << " requested!");
       }
-      return grids_[idx].levelIndexSet(level);
+      return grids_[idx]->levelIndexSet(level);
     }
 
     /// Return the levelIndexSet for grid 0
@@ -267,7 +271,7 @@ namespace Dune
     /// Access to the LeafIndexSet
     const typename HostGrid::LeafIndexSet& leafIndexSet(std::size_t idx) const
     {
-      return grids_[idx].leafIndexSet();
+      return grids_[idx]->leafIndexSet();
     }
 
     /// Return the leafIndexSet for grid 0
@@ -280,7 +284,7 @@ namespace Dune
     /// View for a grid level for All_Partition of the idx'th grid
     typename HostGrid::LevelGridView levelGridView(std::size_t idx, int level) const
     {
-      return grids_[idx].levelGridView(level);
+      return grids_[idx]->levelGridView(level);
     }
 
     /// Return \ref MultiMeshLevelGridView, i.e. a view on the combined
@@ -290,7 +294,7 @@ namespace Dune
     /// View for the leaf grid for All_Partition of the idx'th grid
     typename HostGrid::LeafGridView leafGridView(std::size_t idx) const
     {
-      return grids_[idx].leafGridView();
+      return grids_[idx]->leafGridView();
     }
 
     /// Return \ref MultiMeshLeafGridView, i.e. a view on the combined

src/localrefinement.cc

@@ -9,10 +9,16 @@
 #include <dune/common/exceptions.hh> // We use exceptions
 #include <dune/common/timer.hh>
 
+#if HAVE_ALBERTA
 #include <dune/grid/albertagrid.hh>
+#endif
 #include <dune/grid/yaspgrid.hh>
+#if HAVE_UG
 #include <dune/grid/uggrid.hh>
+#endif
+#if HAVE_DUNE_ALUGRID
 #include <dune/alugrid/grid.hh>
+#endif
 
 #include <dune/multimesh/multimesh.hh>
 #include <dune/multimesh/mmgridfactory.hh>
@@ -55,6 +61,7 @@ int main(int argc, char** argv)
 {
   MPIHelper::instance(argc, argv);
 
+#if HAVE_DUNE_ALUGRID
   FieldVector<double,2> lower_left = {-1.5, -1.5};
   FieldVector<double,2> bbox = {1.5, 1.5};
   std::array<unsigned int,2> num_elements = {2, 2};
@@ -88,4 +95,5 @@ int main(int argc, char** argv)
 
   printGrid2(grid[0]);
   printGrid(grid);
+#endif
 }

src/testiterator.cc

@@ -9,10 +9,16 @@
 #include <dune/common/exceptions.hh> // We use exceptions
 #include <dune/common/timer.hh>
 
+#if HAVE_ALBERTA
 #include <dune/grid/albertagrid.hh>
+#endif
 #include <dune/grid/yaspgrid.hh>
+#if HAVE_UG
 #include <dune/grid/uggrid.hh>
+#endif
+#if HAVE_DUNE_ALUGRID
 #include <dune/alugrid/grid.hh>
+#endif
 
 #include <dune/multimesh/multimesh.hh>
 #include <dune/multimesh/mmgridfactory.hh>
@@ -28,10 +34,10 @@ void printGrid(Grid const& grid)
   volatile std::size_t n = 0;
   Dune::Timer t;
   for (auto const& entities : elements(grid.leafGridView())) {
-    n += grid[0].leafIndexSet().index(entities[0]);
+    n += grid.leafIndexSet(0).index(entities[0]);
     std::cout << "indices = [";
     for (std::size_t i = 0; i < entities.size(); ++i) {
-      std::cout << (i > 0 ? ", " : "") << grid[i].leafIndexSet().index(entities[i]);
+      std::cout << (i > 0 ? ", " : "") << grid.leafIndexSet(i).index(entities[i]);
     }
     std::cout << "]\n";
   }
@@ -55,15 +61,7 @@ int main(int argc, char** argv)
 {
   MPIHelper::instance(argc, argv);
 
-#if GRID == 1
-
-  FieldVector<double,3> lower_left = {-1.5, -1.5, -1.5};
-  FieldVector<double,3> bbox = {1.5, 1.5, 1.5};
-  std::array<int,3> num_elements = {2, 2, 2};
-  using HostGrid = YaspGrid<3, EquidistantOffsetCoordinates<double,3>>;
-  MultiMesh<HostGrid> grid(3, lower_left, bbox, num_elements);
-
-#elif GRID == 2
+#if GRID == 2 && HAVE_DUNE_ALUGRID
 
   FieldVector<double,2> lower_left = {-1.5, -1.5};
   FieldVector<double,2> bbox = {1.5, 1.5};
@@ -76,7 +74,7 @@ int main(int argc, char** argv)
   auto gridPtr = Factory::createSimplexGrid(gridFactory, lower_left, bbox, num_elements);
   auto& grid = *gridPtr;
 
-#elif GRID == 3
+#elif GRID == 3 && HAVE_ALBERTA
 
   assert(argc > 1);
   std::string filename = argv[1];
@@ -85,6 +83,14 @@ int main(int argc, char** argv)
   grid[0].globalRefine(8);
   grid[1].globalRefine(8);
 
+#else
+
+  FieldVector<double,3> lower_left = {-1.5, -1.5, -1.5};
+  FieldVector<double,3> bbox = {1.5, 1.5, 1.5};
+  std::array<int,3> num_elements = {2, 2, 2};
+  using HostGrid = YaspGrid<3, EquidistantOffsetCoordinates<double,3>>;
+  MultiMesh<HostGrid> grid(3, lower_left, bbox, num_elements);
+
 #endif
 
   std::cout << "number of grids = " << grid.size() << "\n";

src/uggrid.cc

@@ -9,7 +9,9 @@
 #include <dune/common/exceptions.hh> // We use exceptions
 #include <dune/common/timer.hh>
 
+#if HAVE_UG
 #include <dune/grid/uggrid.hh>
+#endif
 #include <dune/multimesh/multimesh.hh>
 #include <dune/multimesh/mmgridfactory.hh>
 #include <dune/multimesh/utility/structuredgridbuilder.hh>
@@ -20,6 +22,7 @@ int main(int argc, char** argv)
 {
   MPIHelper::instance(argc, argv);
 
+#if HAVE_UG
   FieldVector<double,2> lower_left = {-1.5, -1.5};
   FieldVector<double,2> bbox = {1.5, 1.5};
   std::array<unsigned int,2> num_elements = {2, 2};
@@ -28,4 +31,5 @@ int main(int argc, char** argv)
   GridFactory<MultiMesh<HostGrid> > factory(1);
 
   auto gridPtr = Factory::createSimplexGrid(factory, lower_left, bbox, num_elements);
+#endif
 }

test/CMakeLists.txt

 dune_add_test(SOURCES testvolume.cc)
 dune_add_test(SOURCES testnumentities.cc)
+dune_add_test(SOURCES testtransform.cc)
\ No newline at end of file

test/testtransform.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#ifdef HAVE_CONFIG_H
+# include "config.h"
+#endif
+
+#include <functional>
+#include <iostream>
+#include <numeric>
+
+#include <dune/common/filledarray.hh>
+#include <dune/common/test/testsuite.hh>
+#include <dune/common/parallel/mpihelper.hh>
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dune/multimesh/multimesh.hh>
+
+using namespace Dune;
+
+template <std::size_t dim, class Test>
+void test_dim(Test& test)
+{
+  FieldVector<double,dim> lower; lower = -1.5;
+  FieldVector<double,dim> upper; upper =  1.5;
+  auto num_elements = filledArray<dim>(2);
+
+  using HostGrid = YaspGrid<dim, EquidistantOffsetCoordinates<double,dim>>;
+  MultiMesh<HostGrid> grid(3, lower, upper, num_elements);
+
+  grid[0].globalRefine(1);
+  grid[1].globalRefine(2);
+  grid[2].globalRefine(3);
+
+  for (auto const& entities : elements(grid.leafGridView())) {
+    auto trafo21 = entities.localGeometry(entities[2], entities[1]);
+    auto trafo20 = entities.localGeometry(entities[2], entities[0]);
+    auto trafo10 = entities.localGeometry(entities[1], entities[0]);
+
+    auto refElem = referenceElement(entities[2].geometry());
+    auto c = refElem.position(0,0);
+
+    // check whether transform.global returns a point inside the entity
+    test.check( refElem.checkInside(trafo21.global(c)) );
+    test.check( refElem.checkInside(trafo20.global(c)) );
+    test.check( refElem.checkInside(trafo10.global(c)) );
+
+    // check whether local is inverse of global
+    test.check( trafo21.local(trafo21.global(c)) == c );
+    test.check( trafo20.local(trafo20.global(c)) == c );
+    test.check( trafo10.local(trafo10.global(c)) == c );
+  }
+}
+
+int main(int argc, char** argv)
+{
+  MPIHelper::instance(argc, argv);
+
+  Dune::TestSuite test;
+  test_dim<1>(test);
+  test_dim<2>(test);
+  test_dim<3>(test);
+  return test.exit();
+}

test/testvolume.cc

@@ -41,10 +41,7 @@ bool test_dim(Test& test)
   // calculate volume by summing up the entity volumes of the smalles leaf entities
   double volume_leaf = 0.0;
   for (auto const& entities : elements(grid.leafGridView())) {
-    auto it_small = std::max_element(entities.begin(), entities.end(),
-      [](auto const& e1, auto const& e2) { return e1.level() < e2.level(); });
-
-    auto geo_small = it_small->geometry();
+    auto geo_small = entities.max().geometry();
     volume_leaf += geo_small.volume();
   }
   std::cout << "volume_leaf(elements<" << dim << ">) = " << volume_leaf << "\n";