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