structured-gridfactory with factory argument

dune/multimesh/mmgridfactory.hh

@@ -41,7 +41,7 @@ namespace Dune
     using GlobalCoordinate = FieldVector<ctype, dimensionworld>;
 
     using BoundarySegment = Dune::BoundarySegment<dimension, dimensionworld>;
-    using ElementParametrization = VirtualFunction<LocalCoordinate,GlobalCoordinate>
+    using ElementParametrization = VirtualFunction<LocalCoordinate,GlobalCoordinate>;
 
   public:
 
@@ -55,6 +55,11 @@ namespace Dune
       : gridFactories_(n, GridFactory<HostGrid>{std::forward<Args>(args)...})
     {}
 
+    // initialize at least 1 grid
+    GridFactory ()
+      : GridFactory{1}
+    {}
+
     /// \brief insert a vertex into the macro grid
     /**
      *  \param[in]  pos  position of the vertex (in world coordinates)
@@ -74,7 +79,7 @@ namespace Dune
                                 const std::vector<unsigned int>& vertices)
     {
       for (auto& gridFactory : gridFactories_)
-        gridFactory.insertElement(vertices);
+        gridFactory.insertElement(type, vertices);
     }
 
     /// \brief Insert a parametrized element into the coarse grid
@@ -83,13 +88,13 @@ namespace Dune
      *   \param[in] vertices The vertices of the new element, using the DUNE numbering
      *   \param[in] elementParametrization A function prescribing the shape of this element
      **/
-    virtual void insertElement (const GeometryType& type,
-                                const std::vector<unsigned int>& vertices,
-                                const std::shared_ptr<ElementParametrization>& elementParametrization)
-    {
-      for (auto& gridFactory : gridFactories_)
-        gridFactory.insertElement(type, vertices, elementParametrization);
-    }
+    // virtual void insertElement (const GeometryType& type,
+    //                             const std::vector<unsigned int>& vertices,
+    //                             const std::shared_ptr<ElementParametrization>& elementParametrization)
+    // {
+    //   for (auto& gridFactory : gridFactories_)
+    //     gridFactory.insertElement(type, vertices, elementParametrization);
+    // }
 
     /// \brief insert a boundary segment into the macro grid
     /**
@@ -124,7 +129,7 @@ namespace Dune
     {
       Grid* multimesh = new Grid{};
       for (auto& gridFactory : gridFactories_)
-        multimesh.grids_.emplace_back(gridFactory.createGrid());
+        multimesh->grids_.emplace_back(gridFactory.createGrid());
       return multimesh;
     }
 

dune/multimesh/multimesh.hh

@@ -123,23 +123,32 @@ namespace Dune
      * \param n The number of host grids to handle by the MultiMesh
      */
     template <class... Args>
-    explicit MultiMesh(std::size_t n, Args&&... args)
+    explicit MultiMesh (std::size_t n, Args&&... args)
     {
       grids_.reserve(n);
       for (std::size_t i = 0; i < n; ++i)
         grids_.emplace_back(std::make_unique<HostGrid>(std::forward<Args>(args)...));
     }
 
+    // initialize 0 grids
+    MultiMesh () = default;
+
+    std::size_t numGrids() const
+    {
+      return grids_.size();
+    }
 
     /// Returns the i'th grid managed by this MultiMesh
     HostGridType& grid (std::size_t i)
     {
+      assert(i < grids_.size());
       return *grids_[i];
     }
 
     /// Returns the i'th grid managed by this MultiMesh
     HostGridType const& grid (std::size_t i) const
     {
+      assert(i < grids_.size());
       return *grids_[i];
     }
 

dune/multimesh/utility/structuredgridfactory.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_MULTIMESH_STRUCTURED_GRIDFACTORY_HH
+#define DUNE_MULTIMESH_STRUCTURED_GRIDFACTORY_HH
+
+/** \file
+    \brief A class to construct structured cube and simplex grids using the grid factory
+ */
+
+#include <algorithm>
+#include <array>
+#include <cstddef>
+#include <cstdlib>
+#include <memory>
+
+#include <dune/common/classname.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/parallel/mpihelper.hh>
+
+#include <dune/grid/common/gridfactory.hh>
+#include <dune/grid/utility/multiindex.hh>
+
+namespace Dune
+{
+  /** \brief Construct structured cube and simplex grids in unstructured grid managers
+   */
+  template <class GridType>
+  class StructuredGridFactoryWrapper
+  {
+    typedef typename GridType::ctype ctype;
+
+    static const int dim = GridType::dimension;
+
+    static const int dimworld = GridType::dimensionworld;
+
+    /** \brief Insert a structured set of vertices into the factory */
+    static void insertVertices(GridFactory<GridType>& factory,
+                               const FieldVector<ctype,dimworld>& lowerLeft,
+                               const FieldVector<ctype,dimworld>& upperRight,
+                               const std::array<unsigned int,dim>& vertices)
+    {
+      FactoryUtilities::MultiIndex<dim> index(vertices);
+
+      // Compute the total number of vertices to be created
+      int numVertices = index.cycle();
+
+      // Create vertices
+      for (int i=0; i<numVertices; i++, ++index) {
+
+        // scale the multiindex to obtain a world position
+        FieldVector<double,dimworld> pos(0);
+        for (int j=0; j<dim; j++)
+          pos[j] = lowerLeft[j] + index[j] * (upperRight[j]-lowerLeft[j])/(vertices[j]-1);
+        for (int j=dim; j<dimworld; j++)
+          pos[j] = lowerLeft[j];
+
+        factory.insertVertex(pos);
+
+      }
+
+    }
+
+    // Compute the index offsets needed to move to the adjacent vertices
+    // in the different coordinate directions
+    static std::array<unsigned int, dim> computeUnitOffsets(const std::array<unsigned int,dim>& vertices)
+    {
+      std::array<unsigned int, dim> unitOffsets;
+      if (dim>0)        // paranoia
+        unitOffsets[0] = 1;
+
+      for (int i=1; i<dim; i++)
+        unitOffsets[i] = unitOffsets[i-1] * vertices[i-1];
+
+      return unitOffsets;
+    }
+
+  public:
+
+    /** \brief Create a structured cube grid
+
+        If the grid dimension is less than the world dimension, the coefficients (dim+1,...,dimworld) in
+        the vertex coordinates are set to the corresponding values of the lowerLeft input argument.
+
+        \param lowerLeft Lower left corner of the grid
+        \param upperRight Upper right corner of the grid
+        \param elements Number of elements in each coordinate direction
+     */
+    static std::unique_ptr<GridType> createCubeGrid(GridFactory<GridType>& factory,
+        const FieldVector<ctype,dimworld>& lowerLeft,
+        const FieldVector<ctype,dimworld>& upperRight,
+        const std::array<unsigned int,dim>& elements)
+    {
+      if (MPIHelper::getCollectiveCommunication().rank() == 0) {
+        // Insert uniformly spaced vertices
+        std::array<unsigned int,dim> vertices = elements;
+        for( size_t i = 0; i < vertices.size(); ++i )
+          vertices[i]++;
+
+        // Insert vertices for structured grid into the factory
+        insertVertices(factory, lowerLeft, upperRight, vertices);
+
+        // Compute the index offsets needed to move to the adjacent
+        // vertices in the different coordinate directions
+        std::array<unsigned int, dim> unitOffsets =
+          computeUnitOffsets(vertices);
+
+        // Compute an element template (the cube at (0,...,0).  All
+        // other cubes are constructed by moving this template around
+        unsigned int nCorners = 1<<dim;
+
+        std::vector<unsigned int> cornersTemplate(nCorners,0);
+
+        for (size_t i=0; i<nCorners; i++)
+          for (int j=0; j<dim; j++)
+            if ( i & (1<<j) )
+              cornersTemplate[i] += unitOffsets[j];
+
+        // Insert elements
+        FactoryUtilities::MultiIndex<dim> index(elements);
+
+        // Compute the total number of elementss to be created
+        int numElements = index.cycle();
+
+        for (int i=0; i<numElements; i++, ++index) {
+
+          // 'base' is the index of the lower left element corner
+          unsigned int base = 0;
+          for (int j=0; j<dim; j++)
+            base += index[j] * unitOffsets[j];
+
+          // insert new element
+          std::vector<unsigned int> corners = cornersTemplate;
+          for (size_t j=0; j<corners.size(); j++)
+            corners[j] += base;
+
+          factory.insertElement(GeometryTypes::cube(dim), corners);
+
+        }
+
+      }       // if(rank == 0)
+
+      // Create the grid and hand it to the calling method
+      return std::unique_ptr<GridType>(factory.createGrid());
+    }
+
+    static std::unique_ptr<GridType> createCubeGrid(const FieldVector<ctype,dimworld>& lowerLeft,
+        const FieldVector<ctype,dimworld>& upperRight,
+        const std::array<unsigned int,dim>& elements)
+    {
+      GridFactory<GridType> factory;
+      return createCubeGrid(factory, lowerLeft, upperRight, elements);
+    }
+
+    /** \brief Create a structured simplex grid
+
+        This works in all dimensions.  The Coxeter-Freudenthal-Kuhn triangulation is
+        used, which splits each cube into dim! (i.e., dim faculty) simplices.  See Allgower and Georg,
+        'Numerical Path Following' for a description.
+
+        If the grid dimension is less than the world dimension, the coefficients (dim+1,...,dimworld) in
+        the vertex coordinates are set to the corresponding values of the lowerLeft input argument.
+
+        \param lowerLeft Lower left corner of the grid
+        \param upperRight Upper right corner of the grid
+        \param elements Number of elements in each coordinate direction
+     */
+    static std::unique_ptr<GridType> createSimplexGrid(GridFactory<GridType>& factory,
+        const FieldVector<ctype,dimworld>& lowerLeft,
+        const FieldVector<ctype,dimworld>& upperRight,
+        const std::array<unsigned int,dim>& elements)
+    {
+      if (MPIHelper::getCollectiveCommunication().rank() == 0) {
+        // Insert uniformly spaced vertices
+        std::array<unsigned int,dim> vertices = elements;
+        for (std::size_t i=0; i<vertices.size(); i++)
+          vertices[i]++;
+
+        insertVertices(factory, lowerLeft, upperRight, vertices);
+
+        // Compute the index offsets needed to move to the adjacent
+        // vertices in the different coordinate directions
+        std::array<unsigned int, dim> unitOffsets =
+          computeUnitOffsets(vertices);
+
+        // Loop over all "cubes", and split up each cube into dim!
+        // (factorial) simplices
+        FactoryUtilities::MultiIndex<dim> elementsIndex(elements);
+        size_t cycle = elementsIndex.cycle();
+
+        for (size_t i=0; i<cycle; ++elementsIndex, i++) {
+
+          // 'base' is the index of the lower left element corner
+          unsigned int base = 0;
+          for (int j=0; j<dim; j++)
+            base += elementsIndex[j] * unitOffsets[j];
+
+          // each permutation of the unit vectors gives a simplex.
+          std::vector<unsigned int> permutation(dim);
+          for (int j=0; j<dim; j++)
+            permutation[j] = j;
+
+          do {
+
+            // Make a simplex
+            std::vector<unsigned int> corners(dim+1);
+            corners[0] = base;
+
+            for (int j=0; j<dim; j++)
+              corners[j+1] =
+                corners[j] + unitOffsets[permutation[j]];
+
+            factory.insertElement(GeometryTypes::simplex(dim), corners);
+
+          } while (std::next_permutation(permutation.begin(),
+                                         permutation.end()));
+
+        }
+
+      }       // if(rank == 0)
+
+      // Create the grid and hand it to the calling method
+      return std::unique_ptr<GridType>(factory.createGrid());
+    }
+
+    // for default-constructable gridFactories
+    static std::unique_ptr<GridType> createSimplexGrid(const FieldVector<ctype,dimworld>& lowerLeft,
+        const FieldVector<ctype,dimworld>& upperRight,
+        const std::array<unsigned int,dim>& elements)
+    {
+      GridFactory<GridType> factory;
+      return createCubeGrid(factory, lowerLeft, upperRight, elements);
+    }
+  };
+
+}  // end namespace Dune
+
+#endif // DUNE_MULTIMESH_STRUCTURED_GRIDFACTORY_HH

src/testiterator.cc

@@ -10,8 +10,14 @@
 #include <dune/common/timer.hh>
 
 #include <dune/grid/yaspgrid.hh>
+#include <dune/grid/uggrid.hh>
+#include <dune/alugrid/grid.hh>
 
 #include <dune/multimesh/multimesh.hh>
+#include <dune/multimesh/mmgridfactory.hh>
+#include <dune/multimesh/utility/structuredgridfactory.hh>
+
+#define GRID 3
 
 using namespace Dune;
 
@@ -42,21 +48,38 @@ 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 = {16, 16, 16};
-
-#if GRID == 1
   using HostGrid = YaspGrid<3, EquidistantOffsetCoordinates<double,3>>;
   MultiMesh<HostGrid> grid(2, lower_left, bbox, num_elements);
 #elif GRID == 2
-  using UGGrid<3>;
-  StructuredGridFactory<MultiMesh<HostGrid> > gridFactory(2);
-  StructuredGridFactory
+  FieldVector<double,3> lower_left = {-1.5, -1.5, -1.5};
+  FieldVector<double,3> bbox = {1.5, 1.5, 1.5};
+  std::array<unsigned int,3> num_elements = {2, 2, 2};
+  using HostGrid = Dune::ALUGrid<3, 3, Dune::simplex, Dune::conforming>;
+
+  using Factory = StructuredGridFactoryWrapper<MultiMesh<HostGrid> >;
+  GridFactory<MultiMesh<HostGrid> > gridFactory(2);
+
+  auto gridPtr = Factory::createSimplexGrid(gridFactory, lower_left, bbox, num_elements);
+  auto& grid = *gridPtr;
+#elif GRID == 3
+  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};
+  using HostGrid = Dune::UGGrid<2>;
 
-  MultiMesh<HostGrid> grid(2);
+  using Factory = StructuredGridFactoryWrapper<MultiMesh<HostGrid> >;
+  GridFactory<MultiMesh<HostGrid> > gridFactory(2);
+
+  auto gridPtr = Factory::createSimplexGrid(gridFactory, lower_left, bbox, num_elements);
+  auto& grid = *gridPtr;
 #endif
 
+  std::cout << "number of grids = " << grid.numGrids() << "\n";
+
   printGrid2(grid.grid(0));
   printGrid(grid);
   grid.grid(0).globalRefine(1);