AMDiS  2.10
The Adaptive Multi-Dimensional Simulation Toolbox
MeshCreator.hpp
1 #pragma once
2 
3 #include <algorithm>
4 #include <array>
5 #include <memory>
6 #include <string>
7 
8 #include <dune/common/filledarray.hh>
9 #include <dune/common/fvector.hh>
10 #include <dune/common/typeutilities.hh>
11 
12 #if HAVE_ALBERTA
13 #include <dune/grid/albertagrid/albertareader.hh>
14 #endif
15 
16 #include <dune/grid/io/file/gmshreader.hh>
17 #include <dune/grid/utility/structuredgridfactory.hh>
18 
19 #if HAVE_DUNE_VTK
20 #include <dune/vtk/vtkreader.hh>
21 #include <dune/vtk/gridcreators/lagrangegridcreator.hh>
22 #endif
23 
24 #if HAVE_DUNE_GMSH4
25 #include <dune/gmsh4/gmsh4reader.hh>
26 #include <dune/gmsh4/gridcreators/lagrangegridcreator.hh>
27 #include <dune/gmsh4/utility/version.hh>
28 #endif
29 
30 #include <amdis/AdaptiveGrid.hpp>
31 #include <amdis/Initfile.hpp>
32 #include <amdis/Output.hpp>
33 #include <amdis/common/Concepts.hpp>
34 #include <amdis/common/Filesystem.hpp>
35 #include <amdis/common/TypeTraits.hpp>
36 
37 namespace Dune
38 {
39  // forward declarations
40  template <int dim, class Coordinates> class YaspGrid;
41  template <class ct, int dim> class EquidistantCoordinates;
42  template <class ct, int dim> class EquidistantOffsetCoordinates;
43  template <class ct, int dim> class TensorProductCoordinates;
44 }
45 
46 
47 namespace AMDiS
48 {
50  template <class G>
51  struct MeshCreator
52  {
53  enum { dimension = G::dimension };
54  enum { dimworld = G::dimensionworld };
55 
56  using Grid = AdaptiveGrid_t<G>;
57  using HostGrid = typename Grid::HostGrid;
58 
59  using ctype = typename Grid::ctype;
60 
62 
65  MeshCreator(std::string const& name)
66  : name_(name)
67  {}
68 
70  static std::shared_ptr<Grid> create(std::string name)
71  {
72  return MeshCreator{name}.create();
73  }
74 
76 
87  std::unique_ptr<Grid> create() const
88  {
89  auto filename = Parameters::get<std::string>(name_ + "->macro file name");
90  auto structured = Parameters::get<std::string>(name_ + "->structured");
91  std::unique_ptr<HostGrid> gridPtr;
92  if (filename) {
93  // read a macro file
94  gridPtr = create_unstructured_grid(filename.value());
95  } else if (structured) {
96  // create structured grids
97  if (structured.value() == "cube") {
98  gridPtr = create_cube_grid();
99  } else if (structured.value() == "simplex") {
100  gridPtr = create_simplex_grid();
101  } else {
102  error_exit("Unknown structured grid type. Must be either `cube` or `simplex` in parameter [meshName]->structured.");
103  }
104  } else {
105  // decide by inspecting the grid type how to create the grid
106  gridPtr = create_by_gridtype<HostGrid>(Dune::PriorityTag<42>{});
107  }
108 
109  // Perform initial refinement and load balance if requested in the initfile
110  auto globalRefinements = Parameters::get<int>(name_ + "->global refinements");
111  if (globalRefinements)
112  gridPtr->globalRefine(globalRefinements.value());
113  auto loadBalance = Parameters::get<bool>(name_ + "->load balance");
114  if (loadBalance && loadBalance.value())
115  gridPtr->loadBalance();
116 
117  return construct(std::move(gridPtr));
118  }
119 
121  std::unique_ptr<HostGrid> create_cube_grid() const
122  {
123  return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells)
124  {
125  return Dune::StructuredGridFactory<HostGrid>::createCubeGrid(lower, upper, numCells);
126  });
127  }
128 
130  std::unique_ptr<HostGrid> create_simplex_grid() const
131  {
132  return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells)
133  {
134  return Dune::StructuredGridFactory<HostGrid>::createSimplexGrid(lower, upper, numCells);
135  });
136  }
137 
139  std::vector<int> const& boundaryIds() const
140  {
141  return boundaryIds_;
142  }
143 
145  std::vector<int> const& elementIds() const
146  {
147  return elementIds_;
148  }
149 
150  private:
151  static std::unique_ptr<Grid> construct(std::unique_ptr<Grid> hostGrid)
152  {
153  return std::move(hostGrid);
154  }
155 
156  template <class HG>
157  static std::unique_ptr<Grid> construct(std::unique_ptr<HG> hostGrid)
158  {
159  return std::make_unique<Grid>(std::move(hostGrid));
160  }
161 
162  // use the structured grid factory to create the grid
163  template <class Size = unsigned int, class Factory>
164  std::unique_ptr<HostGrid> create_structured_grid(Factory factory) const
165  {
166  if constexpr (int(Grid::dimension) == int(Grid::dimensionworld)) {
167  // Lower left corner of the domain
168  Dune::FieldVector<ctype,int(dimworld)> lower(0);
169  // Upper right corner of the domain
170  Dune::FieldVector<ctype,int(dimworld)> upper(1);
171  // number of blocks in each coordinate direction
172  auto numCells = Dune::filledArray<std::size_t(dimension),Size>(1);
173 
174  Parameters::get(name_ + "->min corner", lower);
175  Parameters::get(name_ + "->max corner", upper);
176  Parameters::get(name_ + "->num cells", numCells);
177  return factory(lower, upper, numCells);
178  } else {
179  error_exit("Structured grids can be created for dim == dow only.");
180  return nullptr;
181  }
182  }
183 
184  // read a filename from `[meshName]->macro file name` and determine from the extension the fileformat
185  std::unique_ptr<HostGrid> create_unstructured_grid(std::string const& filename) const
186  {
187  filesystem::path fn(filename);
188  auto ext = fn.extension();
189 
190  if (ext == ".msh") {
191 #if HAVE_DUNE_GMSH4
192  if (Dune::Gmsh4::fileVersion(filename)[0] >= 4)
193  return Dune::Gmsh4Reader<HostGrid, Dune::Gmsh4::LagrangeGridCreator<HostGrid>>::createGridFromFile(filename);
194  else
195 #endif
196  return read_gmsh_file<HostGrid>(filename, Dune::PriorityTag<42>{});
197  }
198 #if HAVE_DUNE_VTK
199  else if (ext == ".vtu") {
200  return Dune::VtkReader<HostGrid, Dune::Vtk::LagrangeGridCreator<HostGrid>>::createGridFromFile(filename);
201  }
202 #endif
203  else if (ext == ".1d" || ext == ".2d" || ext == ".3d" || ext == ".amc") {
204  return read_alberta_file<HostGrid>(filename, Dune::PriorityTag<42>{});
205  }
206  else {
207  error_exit("Cannot read grid file. Unknown file extension.");
208  return {};
209  }
210  }
211 
212  template <class GridType>
213  using SupportsGmshReader = decltype(std::declval<Dune::GridFactory<GridType>>().insertBoundarySegment(
214  std::declval<std::vector<unsigned int>>(),
215  std::declval<std::shared_ptr<Dune::BoundarySegment<GridType::dimension, GridType::dimensionworld> >>()) );
216 
217  template <class GridType, class LC>
218  using SupportsInsertionIndex = decltype(std::declval<Dune::GridFactory<GridType>>().insertionIndex(std::declval<LC>()));
219 
220  // use GmshReader if GridFactory supports insertBoundarySegments
221  template <class GridType = HostGrid,
222  REQUIRES(Dune::Std::is_detected<SupportsGmshReader, GridType>::value)>
223  std::unique_ptr<GridType> read_gmsh_file(std::string const& filename, Dune::PriorityTag<1>) const
224  {
225  Dune::GridFactory<GridType> factory;
226  std::vector<int> boundaryIds, elementIds;
227  Dune::GmshReader<GridType>::read(factory, filename, boundaryIds, elementIds);
228 
229  using HasInsertionIndexEntity
230  = Dune::Std::is_detected<SupportsInsertionIndex, GridType, typename GridType::template Codim<0>::Entity>;
231  using HasInsertionIndexIntersection
232  = Dune::Std::is_detected<SupportsInsertionIndex, GridType, typename GridType::LeafIntersection>;
233 
234  auto gridPtr = factory.createGrid();
235  if ((boundaryIds.empty() && elementIds.empty()) ||
236  (!HasInsertionIndexEntity::value && !HasInsertionIndexIntersection::value))
237  return std::unique_ptr<GridType>(std::move(gridPtr));
238 
239  // map boundaryIds and elementIds read from file to grid indexing.
240  if (!boundaryIds.empty() && HasInsertionIndexIntersection::value)
241  boundaryIds_.resize(gridPtr->numBoundarySegments());
242  if (!elementIds.empty() && HasInsertionIndexEntity::value)
243  elementIds_.resize(gridPtr->size(0));
244 
245  auto const& indexSet = gridPtr->leafIndexSet();
246  for (auto const& e : elements(gridPtr->leafGridView())) {
247  if constexpr(HasInsertionIndexEntity::value) {
248  if (!elementIds.empty())
249  elementIds_[indexSet.index(e)] = elementIds[factory.insertionIndex(e)];
250  }
251 
252  if (boundaryIds.empty() || !e.hasBoundaryIntersections())
253  continue;
254 
255  for (auto const& it : intersections(gridPtr->leafGridView(), e)) {
256  if constexpr(HasInsertionIndexIntersection::value) {
257  if (it.boundary())
258  boundaryIds_[it.boundarySegmentIndex()]
259  = boundaryIds[factory.insertionIndex(it)];
260  }
261  }
262  }
263 
264  return std::unique_ptr<GridType>(std::move(gridPtr));
265  }
266 
267  // fallback if GmshReader cannot be used
268  template <class GridType = HostGrid>
269  std::unique_ptr<GridType> read_gmsh_file(std::string const&, Dune::PriorityTag<0>) const
270  {
271  error_exit("Gmsh reader not supported for this grid type!");
272  return {};
273  }
274 
275  // read from Alberta file
276 
277 #if HAVE_ALBERTA
278  template <class GridType>
279  using IsAlbertaGrid = decltype(std::declval<GridType>().alberta2dune(0,0));
280 
281  // construct the albertagrid directly from a filename
282  template <class GridType = HostGrid,
283  REQUIRES(GridType::dimensionworld == DIM_OF_WORLD),
284  REQUIRES(Dune::Std::is_detected<IsAlbertaGrid, GridType>::value)>
285  std::unique_ptr<GridType> read_alberta_file(std::string const& filename, Dune::PriorityTag<3>) const
286  {
287  return std::make_unique<GridType>(filename);
288  }
289 
290  // use a gridfactory and the generic AlbertaReader
291  template <class GridType = HostGrid,
292  REQUIRES(GridType::dimensionworld == DIM_OF_WORLD)>
293  std::unique_ptr<GridType> read_alberta_file(std::string const& filename, Dune::PriorityTag<2>) const
294  {
295  Dune::GridFactory<GridType> factory;
296  if (Environment::mpiRank() == 0) {
297  Dune::AlbertaReader<GridType> reader;
298  reader.readGrid(filename, factory);
299  }
300  return std::unique_ptr<GridType>{factory.createGrid()};
301  }
302 
303  // error if WORLDDIM not the same as Grid::dimensionworld
304  template <class GridType = HostGrid,
305  REQUIRES(GridType::dimensionworld != DIM_OF_WORLD)>
306  std::unique_ptr<GridType> read_alberta_file(std::string const& filename, Dune::PriorityTag<1>) const
307  {
308  error_exit("AlbertaGrid (and AlbertaReader) require WORLDDIM == Grid::dimensionworld. Change the cmake parameters!");
309  return {};
310  }
311 #else
312  // fallback if no Alberta is available
313  template <class GridType>
314  std::unique_ptr<GridType> read_alberta_file(std::string const&, Dune::PriorityTag<0>) const
315  {
316  error_exit("AlbertaGrid (and AlbertaReader) not available. Set AlbertaFlags to your executable in cmake!");
317  return {};
318  }
319 #endif
320 
321  // yasp grid -> cube
322  template <class GridType = HostGrid,
323  class = typename GridType::YGrid>
324  std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<2>) const
325  {
326  int overlap = 0;
327  Parameters::get(name_ + "->overlap", overlap);
328  std::string periodic(dimension, '0');
329  Parameters::get(name_ + "->periodic", periodic); // e.g. 000 01 111
330 
331  return create_yaspgrid(Types<GridType>{}, overlap, std::bitset<dimension>(periodic));
332  }
333 
334  template <int dim, class ct>
335  std::unique_ptr<HostGrid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::EquidistantCoordinates<ct,dim>>>,
336  int overlap, std::bitset<dimension> const& per) const
337  {
338  return create_structured_grid<int>([&](auto&& /*lower*/, auto&& upper, std::array<int,dimension> const& numCells)
339  {
340  return std::make_unique<HostGrid>(upper, numCells, per, overlap);
341  });
342  }
343 
344  template <int dim, class ct>
345  std::unique_ptr<HostGrid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::EquidistantOffsetCoordinates<ct,dim>>>,
346  int overlap, std::bitset<dimension> const& per) const
347  {
348  return create_structured_grid<int>([&](auto&& lower, auto&& upper, std::array<int,dimension> const& numCells)
349  {
350  return std::make_unique<HostGrid>(lower, upper, numCells, per, overlap);
351  });
352  }
353 
354  template <int dim, class ct>
355  std::unique_ptr<HostGrid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::TensorProductCoordinates<ct,dim>>>,
356  int, std::bitset<dimension> const&) const
357  {
358  error_exit("MeshCreator cannot create YaspGrid with TensorProductCoordinates.");
359  return {};
360  }
361 
362 
363 #if HAVE_DUNE_SPGRID
364  // spgrid -> cube
365  template <class GridType = HostGrid,
366  class = typename GridType::ReferenceCube,
367  class = typename GridType::MultiIndex>
368  std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<1>) const
369  {
370  return create_structured_grid<int>([](auto&& lower, auto&& upper, std::array<int,dimension> const& numCells)
371  {
372  typename GridType::MultiIndex multiIndex(numCells);
373  return std::make_unique<GridType>(lower, upper, multiIndex);
374  });
375  }
376 #endif
377 
378  // final fallback
379  template <class GridType = HostGrid>
380  std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<0>) const
381  {
382  error_exit("Don't know how to create the grid.");
383  return {};
384  }
385 
386  private:
387  std::string name_;
388  mutable std::vector<int> boundaryIds_;
389  mutable std::vector<int> elementIds_;
390  };
391 
392 
393 } // end namespace AMDiS
constexpr bool MultiIndex
A multi-index type.
Definition: Concepts.hpp:149
path extension() const
Returns the file extension path component.
Definition: Filesystem.cpp:83
Definition: AdaptiveGrid.hpp:393
A variadic type list.
Definition: TypeTraits.hpp:146
Definition: AdaptBase.hpp:6
std::unique_ptr< HostGrid > create_simplex_grid() const
Create a structured simplex grid.
Definition: MeshCreator.hpp:130
std::vector< int > const & elementIds() const
Return the filled vector of element ids. NOTE: not all creators support reading this.
Definition: MeshCreator.hpp:145
A creator class for dune grids.
Definition: MeshCreator.hpp:51
static std::shared_ptr< Grid > create(std::string name)
Static create mthod. See create()
Definition: MeshCreator.hpp:70
std::unique_ptr< HostGrid > create_cube_grid() const
Create a structured cube grid.
Definition: MeshCreator.hpp:121
Definition: Filesystem.hpp:11
MeshCreator(std::string const &name)
Construct a new MeshCreator.
Definition: MeshCreator.hpp:65
std::unique_ptr< Grid > create() const
Create a new grid by inspecting the intifile parameter group [meshName]
Definition: MeshCreator.hpp:87
std::vector< int > const & boundaryIds() const
Return the filled vector of boundary ids. NOTE: not all creators support reading this.
Definition: MeshCreator.hpp:139