MeshCreator.hpp 14.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
#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
15

16
17
18
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/grid/utility/structuredgridfactory.hh>

19
20
21
22
23
24
25
26
27
28
#if HAVE_DUNE_VTK
#include <dune/vtk/vtkreader.hh>
#include <dune/vtk/gridcreators/lagrangegridcreator.hh>
#endif

#if HAVE_DUNE_GMSH4
#include <dune/gmsh4/gmsh4reader.hh>
#include <dune/gmsh4/gridcreators/lagrangegridcreator.hh>
#endif

29
#include <amdis/AdaptiveGrid.hpp>
30
31
32
33
34
35
36
#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>

37
38
39
40
41
42
43
44
45
46
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;
}


47
48
49
namespace AMDiS
{
  /// A creator class for dune grids.
50
  template <class G>
51
52
  struct MeshCreator
  {
53
54
55
56
57
    enum { dimension = G::dimension };
    enum { dimworld = G::dimensionworld };

    using Grid = AdaptiveGrid_t<G>;
    using HostGrid = typename Grid::HostGrid;
58
59
60
61
62

    using ctype = typename Grid::ctype;

    /// Construct a new MeshCreator
    /**
Praetorius, Simon's avatar
Praetorius, Simon committed
63
     * \param name  The name of the mesh used in the initfile
64
65
66
67
68
     **/
    MeshCreator(std::string const& name)
      : name_(name)
    {}

69
70
71
72
73
74
    /// Static create mthod. See \ref create()
    static std::shared_ptr<Grid> create(std::string name)
    {
      return MeshCreator{name}.create();
    }

75
76
77
78
79
80
81
82
    /// 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
83
     * using a StructuredGridFactory
84
85
86
     *
     * Otherwise tries to create a grid depending on the grid type.
     **/
87
    std::unique_ptr<Grid> create() const
88
89
90
    {
      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
91
      std::unique_ptr<HostGrid> gridPtr;
92
93
      if (filename) {
        // read a macro file
94
        gridPtr = create_unstructured_grid(filename.value());
95
96
97
      } else if (structured) {
        // create structured grids
        if (structured.value() == "cube") {
98
          gridPtr = create_cube_grid();
99
        } else if (structured && structured.value() == "simplex") {
100
          gridPtr = create_simplex_grid();
101
102
103
104
105
        } 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
106
        gridPtr = create_by_gridtype<HostGrid>(Dune::PriorityTag<42>{});
107
      }
108
109
110
111
112
113
114
115
116

      // 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
117
      return construct(std::move(gridPtr));
118
119
120
    }

    /// Create a structured cube grid
Praetorius, Simon's avatar
Praetorius, Simon committed
121
    std::unique_ptr<HostGrid> create_cube_grid() const
122
123
124
    {
      return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells)
      {
Praetorius, Simon's avatar
Praetorius, Simon committed
125
        return Dune::StructuredGridFactory<HostGrid>::createCubeGrid(lower, upper, numCells);
126
127
128
129
      });
    }

    /// Create a structured simplex grid
Praetorius, Simon's avatar
Praetorius, Simon committed
130
    std::unique_ptr<HostGrid> create_simplex_grid() const
131
132
133
    {
      return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells)
      {
Praetorius, Simon's avatar
Praetorius, Simon committed
134
        return Dune::StructuredGridFactory<HostGrid>::createSimplexGrid(lower, upper, numCells);
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
      });
    }

    /// 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:
151
    static std::unique_ptr<Grid> construct(std::unique_ptr<Grid> hostGrid)
Praetorius, Simon's avatar
Praetorius, Simon committed
152
153
154
155
156
    {
      return std::move(hostGrid);
    }

    template <class HG>
157
    static std::unique_ptr<Grid> construct(std::unique_ptr<HG> hostGrid)
Praetorius, Simon's avatar
Praetorius, Simon committed
158
    {
159
      return std::make_unique<Grid>(std::move(hostGrid));
Praetorius, Simon's avatar
Praetorius, Simon committed
160
161
    }

162
    // use the structured grid factory to create the grid
163
    template <class Size = unsigned int, class Factory>
Praetorius, Simon's avatar
Praetorius, Simon committed
164
    std::unique_ptr<HostGrid> create_structured_grid(Factory factory) const
165
166
167
168
169
170
    {
      // 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
171
      auto numCells = Dune::filledArray<std::size_t(dimension),Size>(1);
172
173
174
175
176
177
178
179

      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
180
    std::unique_ptr<HostGrid> create_unstructured_grid(std::string const& filename) const
181
182
183
184
185
    {
      filesystem::path fn(filename);
      auto ext = fn.extension();

      if (ext == ".msh") {
186
187
188
189
190
#if HAVE_DUNE_GMSH4
        if (Dune::Gmsh4::fileVersion(filename)[0] >= 4)
          return Dune::Gmsh4Reader<HostGrid, Dune::Gmsh4::LagrangeGridCreator<HostGrid>>::createGridFromFile(filename);
        else
#else
Praetorius, Simon's avatar
Praetorius, Simon committed
191
        return read_gmsh_file<HostGrid>(filename, Dune::PriorityTag<42>{});
192
193
#endif
      }
194
#if HAVE_DUNE_VTK
195
196
      else if (ext == ".vtu") {
        return Dune::VtkReader<HostGrid, Dune::Vtk::LagrangeGridCreator<HostGrid>>::createGridFromFile(filename);
197
      }
198
#endif
199
      else if (ext == ".1d" || ext == ".2d" || ext == ".3d" || ext == ".amc") {
Praetorius, Simon's avatar
Praetorius, Simon committed
200
        return read_alberta_file<HostGrid>(filename, Dune::PriorityTag<42>{});
201
202
      }
      else {
203
        error_exit("Cannot read grid file. Unknown file extension.");
204
205
206
207
208
209
210
211
212
        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> >>()) );

213
214
215
    template <class GridType, class LC>
    using SupportsInsertionIndex = decltype(std::declval<Dune::GridFactory<GridType>>().insertionIndex(std::declval<LC>()));

216
    // use GmshReader if GridFactory supports insertBoundarySegments
Praetorius, Simon's avatar
Praetorius, Simon committed
217
    template <class GridType = HostGrid,
218
      REQUIRES(Dune::Std::is_detected<SupportsGmshReader, GridType>::value)>
219
    std::unique_ptr<GridType> read_gmsh_file(std::string const& filename, Dune::PriorityTag<1>) const
220
221
    {
      Dune::GmshReader<GridType> reader;
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
      Dune::GridFactory<GridType> factory;
      std::vector<int> boundaryIds, elementIds;
      reader.read(factory, filename, boundaryIds, elementIds);

      using HasInsertionIndexEntity
        = Dune::Std::is_detected<SupportsInsertionIndex, GridType, typename GridType::template Codim<0>::Entity>;
      using HasInsertionIndexIntersection
        = Dune::Std::is_detected<SupportsInsertionIndex, GridType, typename GridType::LeafIntersection>;

      auto gridPtr = factory.createGrid();
      if ((boundaryIds.empty() && elementIds.empty()) ||
          (!HasInsertionIndexEntity::value && !HasInsertionIndexIntersection::value))
        return std::unique_ptr<GridType>(std::move(gridPtr));

      // map boundaryIds and elementIds read from file to grid indexing.
      if (!boundaryIds.empty() && HasInsertionIndexIntersection::value)
        boundaryIds_.resize(gridPtr->numBoundarySegments());
      if (!elementIds.empty() && HasInsertionIndexEntity::value)
        elementIds_.resize(gridPtr->size(0));

      auto const& indexSet = gridPtr->leafIndexSet();
      for (auto const& e : elements(gridPtr->leafGridView())) {
        if constexpr(HasInsertionIndexEntity::value) {
          if (!elementIds.empty())
            elementIds_[indexSet.index(e)] = elementIds[factory.insertionIndex(e)];
        }

        if (boundaryIds.empty() || !e.hasBoundaryIntersections())
          continue;

        for (auto const& it : intersections(gridPtr->leafGridView(), e)) {
          if constexpr(HasInsertionIndexIntersection::value) {
            if (it.boundary())
              boundaryIds_[it.boundarySegmentIndex()]
                = boundaryIds[factory.insertionIndex(it)];
          }
        }
      }

      return std::unique_ptr<GridType>(std::move(gridPtr));
262
263
264
    }

    // fallback if GmshReader cannot be used
Praetorius, Simon's avatar
Praetorius, Simon committed
265
    template <class GridType = HostGrid>
Praetorius, Simon's avatar
Praetorius, Simon committed
266
    std::unique_ptr<GridType> read_gmsh_file(std::string const&, Dune::PriorityTag<0>) const
267
268
269
270
271
272
273
274
275
276
277
278
    {
      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
279
    template <class GridType = HostGrid,
280
281
282
283
284
285
286
287
      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
288
    template <class GridType = HostGrid,
289
290
291
292
      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;
293
294
295
296
      if (Environment::mpiRank() == 0) {
        Dune::AlbertaReader<GridType> reader;
        reader.readGrid(filename, factory);
      }
297
298
299
300
      return std::unique_ptr<GridType>{factory.createGrid()};
    }

    // error if WORLDDIM not the same as Grid::dimensionworld
Praetorius, Simon's avatar
Praetorius, Simon committed
301
    template <class GridType = HostGrid,
302
303
304
      REQUIRES(GridType::dimensionworld != DIM_OF_WORLD)>
    std::unique_ptr<GridType> read_alberta_file(std::string const& filename, Dune::PriorityTag<1>) const
    {
305
      error_exit("AlbertaGrid (and AlbertaReader) require WORLDDIM == Grid::dimensionworld. Change the cmake parameters!");
306
307
308
309
310
      return {};
    }
#else
    // fallback if no Alberta is available
    template <class GridType>
Praetorius, Simon's avatar
Praetorius, Simon committed
311
    std::unique_ptr<GridType> read_alberta_file(std::string const&, Dune::PriorityTag<0>) const
312
    {
313
      error_exit("AlbertaGrid (and AlbertaReader) not available. Set AlbertaFlags to your executable in cmake!");
314
315
316
317
318
319
320
      return {};
    }
#endif


#if HAVE_ALBERTA
    // albertagrid -> simplex
Praetorius, Simon's avatar
Praetorius, Simon committed
321
    template <class GridType = HostGrid,
322
323
324
      REQUIRES(Dune::Std::is_detected<IsAlbertaGrid, GridType>::value)>
    std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<3>) const
    {
325
326
      return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells)
      {
Praetorius, Simon's avatar
Praetorius, Simon committed
327
        return MacroGridFactory<GridType>::createSimplexGrid(lower, upper, numCells);
328
      });
329
330
331
332
    }
#endif

    // yasp grid -> cube
Praetorius, Simon's avatar
Praetorius, Simon committed
333
    template <class GridType = HostGrid,
334
335
336
      class = typename GridType::YGrid>
    std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<2>) const
    {
337
338
339
340
341
342
343
344
345
      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
346
    std::unique_ptr<HostGrid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::EquidistantCoordinates<ct,dim>>>,
347
348
349
350
                                          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
351
        return std::make_unique<HostGrid>(upper, numCells, per, overlap);
352
353
354
355
      });
    }

    template <int dim, class ct>
Praetorius, Simon's avatar
Praetorius, Simon committed
356
    std::unique_ptr<HostGrid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::EquidistantOffsetCoordinates<ct,dim>>>,
357
358
359
360
                                          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
361
        return std::make_unique<HostGrid>(lower, upper, numCells, per, overlap);
362
      });
363
364
    }

365
    template <int dim, class ct>
Praetorius, Simon's avatar
Praetorius, Simon committed
366
    std::unique_ptr<HostGrid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::TensorProductCoordinates<ct,dim>>>,
Praetorius, Simon's avatar
Praetorius, Simon committed
367
                                          int, std::bitset<dimension> const&) const
368
369
370
371
372
373
    {
      error_exit("MeshCreator cannot create YaspGrid with TensorProductCoordinates.");
      return {};
    }


374
375
#if HAVE_DUNE_SPGRID
    // spgrid -> cube
Praetorius, Simon's avatar
Praetorius, Simon committed
376
    template <class GridType = HostGrid,
377
378
379
380
      class = typename GridType::ReferenceCube,
      class = typename GridType::MultiIndex>
    std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<1>) const
    {
381
      return create_structured_grid<int>([](auto&& lower, auto&& upper, std::array<int,dimension> const& numCells)
382
      {
383
        typename GridType::MultiIndex multiIndex(numCells);
384
385
386
387
388
389
        return std::make_unique<GridType>(lower, upper, multiIndex);
      });
    }
#endif

    // final fallback
Praetorius, Simon's avatar
Praetorius, Simon committed
390
    template <class GridType = HostGrid>
391
392
393
394
395
396
397
398
399
400
401
402
403
404
    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