structuredgridbuilder.hh 8.06 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
    /// \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
     **/
119
    static void createSimplexGrid (
Praetorius, Simon's avatar
Praetorius, Simon committed
120
        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
    }

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

  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;
    }
221
222
223
224
225
  };

}  // end namespace Dune

#endif // DUNE_MULTIMESH_STRUCTURED_GRIDFACTORY_HH