structuredgridbuilder.hh 7.97 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
// -*- 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

#include <algorithm>
#include <array>
#include <cstddef>
#include <cstdlib>
#include <memory>

#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
{
Praetorius, Simon's avatar
Praetorius, Simon committed
20
  /// \brief Construct structured cube and simplex grids in unstructured grid managers
21
  template <class GridType>
Praetorius, Simon's avatar
Praetorius, Simon committed
22
  class StructuredGridBuilder
23
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
24
    using ctype = typename GridType::ctype;
25 26 27 28 29

    static const int dim = GridType::dimension;
    static const int dimworld = GridType::dimensionworld;

  public:
Praetorius, Simon's avatar
Praetorius, Simon committed
30 31 32 33 34 35 36 37 38 39
    /// \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
     **/
Praetorius, Simon's avatar
Praetorius, Simon committed
40
    static void createCubeGrid (
Praetorius, Simon's avatar
Praetorius, Simon committed
41
        GridFactory<GridType>& factory,
42 43 44 45 46 47 48
        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;
Praetorius, Simon's avatar
Praetorius, Simon committed
49
        for (std::size_t i = 0; i < vertices.size(); ++i)
50 51 52 53 54 55 56
          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
Praetorius, Simon's avatar
Praetorius, Simon committed
57
        std::array<unsigned int,dim> unitOffsets =
58 59 60 61 62 63 64 65
          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);

Praetorius, Simon's avatar
Praetorius, Simon committed
66 67 68
        for (std::size_t i = 0; i < nCorners; ++i)
          for (int j = 0; j < dim; ++j)
            if (i & (1<<j))
69 70 71 72 73 74 75 76
              cornersTemplate[i] += unitOffsets[j];

        // Insert elements
        FactoryUtilities::MultiIndex<dim> index(elements);

        // Compute the total number of elementss to be created
        int numElements = index.cycle();

Praetorius, Simon's avatar
Praetorius, Simon committed
77
        for (int i = 0; i < numElements; ++i, ++index) {
78 79
          // 'base' is the index of the lower left element corner
          unsigned int base = 0;
Praetorius, Simon's avatar
Praetorius, Simon committed
80
          for (int j = 0; j < dim; ++j)
81 82 83 84
            base += index[j] * unitOffsets[j];

          // insert new element
          std::vector<unsigned int> corners = cornersTemplate;
Praetorius, Simon's avatar
Praetorius, Simon committed
85
          for (std::size_t j = 0; j < corners.size(); ++j)
86 87 88 89
            corners[j] += base;

          factory.insertElement(GeometryTypes::cube(dim), corners);
        }
Praetorius, Simon's avatar
Praetorius, Simon committed
90
      }
91 92
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
93 94
    static std::unique_ptr<GridType> createCubeGrid (
        const FieldVector<ctype,dimworld>& lowerLeft,
95 96 97 98
        const FieldVector<ctype,dimworld>& upperRight,
        const std::array<unsigned int,dim>& elements)
    {
      GridFactory<GridType> factory;
Praetorius, Simon's avatar
Praetorius, Simon committed
99 100
      createCubeGrid(factory, lowerLeft, upperRight, elements);
      return std::unique_ptr<GridType>(factory.createGrid());
101 102
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
103 104 105 106 107 108 109 110 111 112 113 114 115 116
    /// \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
     **/
117
    static void createSimplexGrid (
Praetorius, Simon's avatar
Praetorius, Simon committed
118
        GridFactory<GridType>& factory,
119 120 121 122 123 124 125
        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;
Praetorius, Simon's avatar
Praetorius, Simon committed
126
        for (std::size_t i = 0; i < vertices.size(); ++i)
127 128 129 130 131 132
          vertices[i]++;

        insertVertices(factory, lowerLeft, upperRight, vertices);

        // Compute the index offsets needed to move to the adjacent
        // vertices in the different coordinate directions
Praetorius, Simon's avatar
Praetorius, Simon committed
133
        std::array<unsigned int,dim> unitOffsets =
134 135 136 137 138
          computeUnitOffsets(vertices);

        // Loop over all "cubes", and split up each cube into dim!
        // (factorial) simplices
        FactoryUtilities::MultiIndex<dim> elementsIndex(elements);
Praetorius, Simon's avatar
Praetorius, Simon committed
139
        std::size_t cycle = elementsIndex.cycle();
140

Praetorius, Simon's avatar
Praetorius, Simon committed
141
        for (std::size_t i = 0; i < cycle; ++elementsIndex, ++i) {
142 143
          // 'base' is the index of the lower left element corner
          unsigned int base = 0;
Praetorius, Simon's avatar
Praetorius, Simon committed
144
          for (int j = 0; j < dim; ++j)
145 146 147 148
            base += elementsIndex[j] * unitOffsets[j];

          // each permutation of the unit vectors gives a simplex.
          std::vector<unsigned int> permutation(dim);
Praetorius, Simon's avatar
Praetorius, Simon committed
149
          for (int j = 0; j < dim; ++j)
150 151 152 153 154 155 156
            permutation[j] = j;

          do {
            // Make a simplex
            std::vector<unsigned int> corners(dim+1);
            corners[0] = base;

Praetorius, Simon's avatar
Praetorius, Simon committed
157
            for (int j = 0; j < dim; ++j)
158 159 160 161 162 163 164 165
              corners[j+1] =
                corners[j] + unitOffsets[permutation[j]];

            factory.insertElement(GeometryTypes::simplex(dim), corners);

          } while (std::next_permutation(permutation.begin(),
                                         permutation.end()));
        }
Praetorius, Simon's avatar
Praetorius, Simon committed
166
      }
167 168
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
169 170
    static std::unique_ptr<GridType> createSimplexGrid (
        const FieldVector<ctype,dimworld>& lowerLeft,
171 172 173 174
        const FieldVector<ctype,dimworld>& upperRight,
        const std::array<unsigned int,dim>& elements)
    {
      GridFactory<GridType> factory;
175 176
      createSimplexGrid(factory, lowerLeft, upperRight, elements);
      return std::unique_ptr<GridType>(factory.createGrid());
177
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218

  private:
    /// 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);
      }
    }

    /// \brief 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;
    }
219 220 221 222 223
  };

}  // end namespace Dune

#endif // DUNE_MULTIMESH_STRUCTURED_GRIDFACTORY_HH