MeshCreator.hpp 11.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#pragma once

#include <algorithm>
#include <array>
#include <memory>
#include <string>

#include <dune/common/filledarray.hh>
#include <dune/common/fvector.hh>
#include <dune/common/typeutilities.hh>

#if HAVE_ALBERTA
#include <dune/grid/albertagrid/albertareader.hh>
#endif
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/grid/utility/structuredgridfactory.hh>

18
#include <amdis/AdaptiveGrid.hpp>
19
20
21
22
23
24
25
#include <amdis/Initfile.hpp>
#include <amdis/Output.hpp>
#include <amdis/common/Concepts.hpp>
#include <amdis/common/Filesystem.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/utility/MacroGridFactory.hpp>

26
27
28
29
30
31
32
33
34
35
namespace Dune
{
  // forward declarations
  template <int dim, class Coordinates> class YaspGrid;
  template <class ct, int dim> class EquidistantCoordinates;
  template <class ct, int dim> class EquidistantOffsetCoordinates;
  template <class ct, int dim> class TensorProductCoordinates;
}


36
37
38
39
40
41
42
43
44
45
namespace AMDiS
{
  /// A creator class for dune grids.
  template <class Grid>
  struct MeshCreator
  {
    enum { dimension = Grid::dimension };
    enum { dimworld = Grid::dimensionworld };

    using ctype = typename Grid::ctype;
Praetorius, Simon's avatar
Praetorius, Simon committed
46
    using HostGrid = typename AdaptiveGrid_t<Grid>::HostGrid;
47
48
49

    /// Construct a new MeshCreator
    /**
Praetorius, Simon's avatar
Praetorius, Simon committed
50
     * \param name  The name of the mesh used in the initfile
51
52
53
54
55
56
57
58
59
60
61
62
63
     **/
    MeshCreator(std::string const& name)
      : name_(name)
    {}

    /// Create a new grid by inspecting the intifile parameter group `[meshName]`
    /**
     * Reads first the parameter `[meshName]->macro file name` and if set
     * - reads the grid from file
     *
     * Otherwise reads the parameter `[meshName]->structured` and if set, creates:
     * - cube: a structured cube grid
     * - simplex: a structured simplex grid
64
     * using a StructuredGridFactory
65
66
67
     *
     * Otherwise tries to create a grid depending on the grid type.
     **/
Praetorius, Simon's avatar
Praetorius, Simon committed
68
    std::shared_ptr<Grid> create() const
69
70
71
    {
      auto filename = Parameters::get<std::string>(name_ + "->macro file name");
      auto structured = Parameters::get<std::string>(name_ + "->structured");
Praetorius, Simon's avatar
Praetorius, Simon committed
72
      std::unique_ptr<HostGrid> gridPtr;
73
74
      if (filename) {
        // read a macro file
75
        gridPtr = create_unstructured_grid(filename.value());
76
77
78
      } else if (structured) {
        // create structured grids
        if (structured.value() == "cube") {
79
          gridPtr = create_cube_grid();
80
        } else if (structured && structured.value() == "simplex") {
81
          gridPtr = create_simplex_grid();
82
83
84
85
86
        } else {
          error_exit("Unknown structured grid type. Must be either `cube` or `simplex` in parameter [meshName]->structured.");
        }
      } else {
        // decide by inspecting the grid type how to create the grid
Praetorius, Simon's avatar
Praetorius, Simon committed
87
        gridPtr = create_by_gridtype<HostGrid>(Dune::PriorityTag<42>{});
88
      }
89
90
91
92
93
94
95
96
97

      // Perform initial refinement and load balance if requested in the initfile
      auto globalRefinements = Parameters::get<int>(name_ + "->global refinements");
      if (globalRefinements)
        gridPtr->globalRefine(globalRefinements.value());
      auto loadBalance = Parameters::get<bool>(name_ + "->load balance");
      if (loadBalance && loadBalance.value())
        gridPtr->loadBalance();

Praetorius, Simon's avatar
Praetorius, Simon committed
98
      return construct(std::move(gridPtr));
99
100
101
    }

    /// Create a structured cube grid
Praetorius, Simon's avatar
Praetorius, Simon committed
102
    std::unique_ptr<HostGrid> create_cube_grid() const
103
104
105
    {
      return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells)
      {
Praetorius, Simon's avatar
Praetorius, Simon committed
106
        return Dune::StructuredGridFactory<HostGrid>::createCubeGrid(lower, upper, numCells);
107
108
109
110
      });
    }

    /// Create a structured simplex grid
Praetorius, Simon's avatar
Praetorius, Simon committed
111
    std::unique_ptr<HostGrid> create_simplex_grid() const
112
113
114
    {
      return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells)
      {
Praetorius, Simon's avatar
Praetorius, Simon committed
115
        return Dune::StructuredGridFactory<HostGrid>::createSimplexGrid(lower, upper, numCells);
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
      });
    }

    /// Return the filled vector of boundary ids. NOTE: not all creators support reading this.
    std::vector<int> const& boundaryIds() const
    {
      return boundaryIds_;
    }

    /// Return the filled vector of element ids. NOTE: not all creators support reading this.
    std::vector<int> const& elementIds() const
    {
      return elementIds_;
    }

  private:
Praetorius, Simon's avatar
Praetorius, Simon committed
132
133
134
135
136
137
138
139
140
141
142
    static std::shared_ptr<Grid> construct(std::unique_ptr<Grid> hostGrid)
    {
      return std::move(hostGrid);
    }

    template <class HG>
    static std::shared_ptr<Grid> construct(std::unique_ptr<HG> hostGrid)
    {
      return std::make_shared<Grid>(std::move(hostGrid));
    }

143
    // use the structured grid factory to create the grid
144
    template <class Size = unsigned int, class Factory>
Praetorius, Simon's avatar
Praetorius, Simon committed
145
    std::unique_ptr<HostGrid> create_structured_grid(Factory factory) const
146
147
148
149
150
151
    {
      // Lower left corner of the domain
      Dune::FieldVector<ctype,int(dimworld)> lower(0);
      // Upper right corner of the domain
      Dune::FieldVector<ctype,int(dimworld)> upper(1);
      // number of blocks in each coordinate direction
152
      auto numCells = Dune::filledArray<std::size_t(dimension),Size>(1);
153
154
155
156
157
158
159
160

      Parameters::get(name_ + "->min corner", lower);
      Parameters::get(name_ + "->max corner", upper);
      Parameters::get(name_ + "->num cells", numCells);
      return factory(lower, upper, numCells);
    }

    // read a filename from `[meshName]->macro file name` and determine from the extension the fileformat
Praetorius, Simon's avatar
Praetorius, Simon committed
161
    std::unique_ptr<HostGrid> create_unstructured_grid(std::string const& filename) const
162
163
164
165
166
    {
      filesystem::path fn(filename);
      auto ext = fn.extension();

      if (ext == ".msh") {
Praetorius, Simon's avatar
Praetorius, Simon committed
167
        return read_gmsh_file<HostGrid>(filename, Dune::PriorityTag<42>{});
168
169
      }
      else if (ext == ".1d" || ext == ".2d" || ext == ".3d" || ext == ".amc") {
Praetorius, Simon's avatar
Praetorius, Simon committed
170
        return read_alberta_file<HostGrid>(filename, Dune::PriorityTag<42>{});
171
172
      }
      else {
173
        error_exit("Cannot read grid file. Unknown file extension.");
174
175
176
177
178
179
180
181
182
183
        return {};
      }
    }

    template <class GridType>
    using SupportsGmshReader = decltype(std::declval<Dune::GridFactory<GridType>>().insertBoundarySegment(
      std::declval<std::vector<unsigned int>>(),
      std::declval<std::shared_ptr<Dune::BoundarySegment<GridType::dimension, GridType::dimensionworld> >>()) );

    // use GmshReader if GridFactory supports insertBoundarySegments
Praetorius, Simon's avatar
Praetorius, Simon committed
184
    template <class GridType = HostGrid,
185
      REQUIRES(Dune::Std::is_detected<SupportsGmshReader, GridType>::value)>
186
    std::unique_ptr<GridType> read_gmsh_file(std::string const& filename, Dune::PriorityTag<1>) const
187
188
    {
      Dune::GmshReader<GridType> reader;
189
      return std::unique_ptr<GridType>{reader.read(filename)}; // , boundaryIds_, elementIds_)};
190
191
192
    }

    // fallback if GmshReader cannot be used
Praetorius, Simon's avatar
Praetorius, Simon committed
193
    template <class GridType = HostGrid>
194
195
196
197
198
199
200
201
202
203
204
205
206
    std::unique_ptr<GridType> read_gmsh_file(std::string const& filename, Dune::PriorityTag<0>) const
    {
      error_exit("Gmsh reader not supported for this grid type!");
      return {};
    }

    // read from Alberta file

#if HAVE_ALBERTA
    template <class GridType>
    using IsAlbertaGrid = decltype(std::declval<GridType>().alberta2dune(0,0));

    // construct the albertagrid directly from a filename
Praetorius, Simon's avatar
Praetorius, Simon committed
207
    template <class GridType = HostGrid,
208
209
210
211
212
213
214
215
      REQUIRES(GridType::dimensionworld == DIM_OF_WORLD),
      REQUIRES(Dune::Std::is_detected<IsAlbertaGrid, GridType>::value)>
    std::unique_ptr<GridType> read_alberta_file(std::string const& filename, Dune::PriorityTag<3>) const
    {
      return std::make_unique<GridType>(filename);
    }

    // use a gridfactory and the generic AlbertaReader
Praetorius, Simon's avatar
Praetorius, Simon committed
216
    template <class GridType = HostGrid,
217
218
219
220
      REQUIRES(GridType::dimensionworld == DIM_OF_WORLD)>
    std::unique_ptr<GridType> read_alberta_file(std::string const& filename, Dune::PriorityTag<2>) const
    {
      Dune::GridFactory<GridType> factory;
221
222
223
224
      if (Environment::mpiRank() == 0) {
        Dune::AlbertaReader<GridType> reader;
        reader.readGrid(filename, factory);
      }
225
226
227
228
      return std::unique_ptr<GridType>{factory.createGrid()};
    }

    // error if WORLDDIM not the same as Grid::dimensionworld
Praetorius, Simon's avatar
Praetorius, Simon committed
229
    template <class GridType = HostGrid,
230
231
232
      REQUIRES(GridType::dimensionworld != DIM_OF_WORLD)>
    std::unique_ptr<GridType> read_alberta_file(std::string const& filename, Dune::PriorityTag<1>) const
    {
233
      error_exit("AlbertaGrid (and AlbertaReader) require WORLDDIM == Grid::dimensionworld. Change the cmake parameters!");
234
235
236
237
238
239
240
      return {};
    }
#else
    // fallback if no Alberta is available
    template <class GridType>
    std::unique_ptr<GridType> read_alberta_file(std::string const& filename, Dune::PriorityTag<0>) const
    {
241
      error_exit("AlbertaGrid (and AlbertaReader) not available. Set AlbertaFlags to your executable in cmake!");
242
243
244
245
246
247
248
      return {};
    }
#endif


#if HAVE_ALBERTA
    // albertagrid -> simplex
Praetorius, Simon's avatar
Praetorius, Simon committed
249
    template <class GridType = HostGrid,
250
251
252
      REQUIRES(Dune::Std::is_detected<IsAlbertaGrid, GridType>::value)>
    std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<3>) const
    {
253
254
      return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells)
      {
Praetorius, Simon's avatar
Praetorius, Simon committed
255
        return MacroGridFactory<GridType>::createSimplexGrid(lower, upper, numCells);
256
      });
257
258
259
260
    }
#endif

    // yasp grid -> cube
Praetorius, Simon's avatar
Praetorius, Simon committed
261
    template <class GridType = HostGrid,
262
263
264
      class = typename GridType::YGrid>
    std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<2>) const
    {
265
266
267
268
269
270
271
272
273
      int overlap = 0;
      Parameters::get(name_ + "->overlap", overlap);
      std::string periodic(dimension, '0');
      Parameters::get(name_ + "->periodic", periodic); // e.g. 000 01 111

      return create_yaspgrid(Types<GridType>{}, overlap, std::bitset<dimension>(periodic));
    }

    template <int dim, class ct>
Praetorius, Simon's avatar
Praetorius, Simon committed
274
    std::unique_ptr<HostGrid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::EquidistantCoordinates<ct,dim>>>,
275
276
277
278
                                          int overlap, std::bitset<dimension> const& per) const
    {
      return create_structured_grid<int>([&](auto&& /*lower*/, auto&& upper, std::array<int,dimension> const& numCells)
      {
Praetorius, Simon's avatar
Praetorius, Simon committed
279
        return std::make_unique<HostGrid>(upper, numCells, per, overlap);
280
281
282
283
      });
    }

    template <int dim, class ct>
Praetorius, Simon's avatar
Praetorius, Simon committed
284
    std::unique_ptr<HostGrid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::EquidistantOffsetCoordinates<ct,dim>>>,
285
286
287
288
                                          int overlap, std::bitset<dimension> const& per) const
    {
      return create_structured_grid<int>([&](auto&& lower, auto&& upper, std::array<int,dimension> const& numCells)
      {
Praetorius, Simon's avatar
Praetorius, Simon committed
289
        return std::make_unique<HostGrid>(lower, upper, numCells, per, overlap);
290
      });
291
292
    }

293
    template <int dim, class ct>
Praetorius, Simon's avatar
Praetorius, Simon committed
294
    std::unique_ptr<HostGrid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::TensorProductCoordinates<ct,dim>>>,
295
296
297
298
299
300
301
                                          int overlap, std::bitset<dimension> const& per) const
    {
      error_exit("MeshCreator cannot create YaspGrid with TensorProductCoordinates.");
      return {};
    }


302
303
#if HAVE_DUNE_SPGRID
    // spgrid -> cube
Praetorius, Simon's avatar
Praetorius, Simon committed
304
    template <class GridType = HostGrid,
305
306
307
308
      class = typename GridType::ReferenceCube,
      class = typename GridType::MultiIndex>
    std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<1>) const
    {
309
      return create_structured_grid<int>([](auto&& lower, auto&& upper, std::array<int,dimension> const& numCells)
310
      {
311
        typename GridType::MultiIndex multiIndex(numCells);
312
313
314
315
316
317
        return std::make_unique<GridType>(lower, upper, multiIndex);
      });
    }
#endif

    // final fallback
Praetorius, Simon's avatar
Praetorius, Simon committed
318
    template <class GridType = HostGrid>
319
320
321
322
323
324
325
326
327
328
329
330
331
332
    std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<0>) const
    {
      error_exit("Don't know how to create the grid.");
      return {};
    }

  private:
    std::string name_;
    mutable std::vector<int> boundaryIds_;
    mutable std::vector<int> elementIds_;
  };


} // end namespace AMDiS