structuredgridbuilder.hh 8.14 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 40 41
    /// \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,
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 93 94 95

      // Create the grid and hand it to the calling method
      return std::unique_ptr<GridType>(factory.createGrid());
    }

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

Praetorius, Simon's avatar
Praetorius, Simon committed
105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
    /// \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,
121 122 123 124 125 126 127
        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
128
        for (std::size_t i = 0; i < vertices.size(); ++i)
129 130 131 132 133 134
          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
135
        std::array<unsigned int,dim> unitOffsets =
136 137 138 139 140
          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
141
        std::size_t cycle = elementsIndex.cycle();
142

Praetorius, Simon's avatar
Praetorius, Simon committed
143
        for (std::size_t i = 0; i < cycle; ++elementsIndex, ++i) {
144 145
          // 'base' is the index of the lower left element corner
          unsigned int base = 0;
Praetorius, Simon's avatar
Praetorius, Simon committed
146
          for (int j = 0; j < dim; ++j)
147 148 149 150
            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
151
          for (int j = 0; j < dim; ++j)
152 153 154 155 156 157 158
            permutation[j] = j;

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

Praetorius, Simon's avatar
Praetorius, Simon committed
159
            for (int j = 0; j < dim; ++j)
160 161 162 163 164 165 166 167
              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
168
      }
169 170 171 172 173

      // Create the grid and hand it to the calling method
      return std::unique_ptr<GridType>(factory.createGrid());
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
174 175
    static std::unique_ptr<GridType> createSimplexGrid (
        const FieldVector<ctype,dimworld>& lowerLeft,
176 177 178 179 180 181
        const FieldVector<ctype,dimworld>& upperRight,
        const std::array<unsigned int,dim>& elements)
    {
      GridFactory<GridType> factory;
      return createCubeGrid(factory, lowerLeft, upperRight, elements);
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
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 219 220 221 222

  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;
    }
223 224 225 226 227
  };

}  // end namespace Dune

#endif // DUNE_MULTIMESH_STRUCTURED_GRIDFACTORY_HH