AMDiS  0.3
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 #include <amdis/utility/MacroGridFactory.hpp>
37 
38 namespace Dune
39 {
40  // forward declarations
41  template <int dim, class Coordinates> class YaspGrid;
42  template <class ct, int dim> class EquidistantCoordinates;
43  template <class ct, int dim> class EquidistantOffsetCoordinates;
44  template <class ct, int dim> class TensorProductCoordinates;
45 }
46 
47 
48 namespace AMDiS
49 {
51  template <class G>
52  struct MeshCreator
53  {
54  enum { dimension = G::dimension };
55  enum { dimworld = G::dimensionworld };
56 
57  using Grid = AdaptiveGrid_t<G>;
58  using HostGrid = typename Grid::HostGrid;
59 
60  using ctype = typename Grid::ctype;
61 
63 
66  MeshCreator(std::string const& name)
67  : name_(name)
68  {}
69 
71  static std::shared_ptr<Grid> create(std::string name)
72  {
73  return MeshCreator{name}.create();
74  }
75 
77 
88  std::unique_ptr<Grid> create() const
89  {
90  auto filename = Parameters::get<std::string>(name_ + "->macro file name");
91  auto structured = Parameters::get<std::string>(name_ + "->structured");
92  std::unique_ptr<HostGrid> gridPtr;
93  if (filename) {
94  // read a macro file
95  gridPtr = create_unstructured_grid(filename.value());
96  } else if (structured) {
97  // create structured grids
98  if (structured.value() == "cube") {
99  gridPtr = create_cube_grid();
100  } else if (structured && structured.value() == "simplex") {
101  gridPtr = create_simplex_grid();
102  } else {
103  error_exit("Unknown structured grid type. Must be either `cube` or `simplex` in parameter [meshName]->structured.");
104  }
105  } else {
106  // decide by inspecting the grid type how to create the grid
107  gridPtr = create_by_gridtype<HostGrid>(Dune::PriorityTag<42>{});
108  }
109 
110  // Perform initial refinement and load balance if requested in the initfile
111  auto globalRefinements = Parameters::get<int>(name_ + "->global refinements");
112  if (globalRefinements)
113  gridPtr->globalRefine(globalRefinements.value());
114  auto loadBalance = Parameters::get<bool>(name_ + "->load balance");
115  if (loadBalance && loadBalance.value())
116  gridPtr->loadBalance();
117 
118  return construct(std::move(gridPtr));
119  }
120 
122  std::unique_ptr<HostGrid> create_cube_grid() const
123  {
124  return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells)
125  {
126  return Dune::StructuredGridFactory<HostGrid>::createCubeGrid(lower, upper, numCells);
127  });
128  }
129 
131  std::unique_ptr<HostGrid> create_simplex_grid() const
132  {
133  return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells)
134  {
135  return Dune::StructuredGridFactory<HostGrid>::createSimplexGrid(lower, upper, numCells);
136  });
137  }
138 
140  std::vector<int> const& boundaryIds() const
141  {
142  return boundaryIds_;
143  }
144 
146  std::vector<int> const& elementIds() const
147  {
148  return elementIds_;
149  }
150 
151  private:
152  static std::unique_ptr<Grid> construct(std::unique_ptr<Grid> hostGrid)
153  {
154  return std::move(hostGrid);
155  }
156 
157  template <class HG>
158  static std::unique_ptr<Grid> construct(std::unique_ptr<HG> hostGrid)
159  {
160  return std::make_unique<Grid>(std::move(hostGrid));
161  }
162 
163  // use the structured grid factory to create the grid
164  template <class Size = unsigned int, class Factory>
165  std::unique_ptr<HostGrid> create_structured_grid(Factory factory) const
166  {
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  }
179 
180  // read a filename from `[meshName]->macro file name` and determine from the extension the fileformat
181  std::unique_ptr<HostGrid> create_unstructured_grid(std::string const& filename) const
182  {
183  filesystem::path fn(filename);
184  auto ext = fn.extension();
185 
186  if (ext == ".msh") {
187 #if HAVE_DUNE_GMSH4
188  if (Dune::Gmsh4::fileVersion(filename)[0] >= 4)
189  return Dune::Gmsh4Reader<HostGrid, Dune::Gmsh4::LagrangeGridCreator<HostGrid>>::createGridFromFile(filename);
190  else
191 #endif
192  return read_gmsh_file<HostGrid>(filename, Dune::PriorityTag<42>{});
193  }
194 #if HAVE_DUNE_VTK
195  else if (ext == ".vtu") {
196  return Dune::VtkReader<HostGrid, Dune::Vtk::LagrangeGridCreator<HostGrid>>::createGridFromFile(filename);
197  }
198 #endif
199  else if (ext == ".1d" || ext == ".2d" || ext == ".3d" || ext == ".amc") {
200  return read_alberta_file<HostGrid>(filename, Dune::PriorityTag<42>{});
201  }
202  else {
203  error_exit("Cannot read grid file. Unknown file extension.");
204  return {};
205  }
206  }
207 
208  template <class GridType>
209  using SupportsGmshReader = decltype(std::declval<Dune::GridFactory<GridType>>().insertBoundarySegment(
210  std::declval<std::vector<unsigned int>>(),
211  std::declval<std::shared_ptr<Dune::BoundarySegment<GridType::dimension, GridType::dimensionworld> >>()) );
212 
213  template <class GridType, class LC>
214  using SupportsInsertionIndex = decltype(std::declval<Dune::GridFactory<GridType>>().insertionIndex(std::declval<LC>()));
215 
216  // use GmshReader if GridFactory supports insertBoundarySegments
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  Dune::GridFactory<GridType> factory;
223  std::vector<int> boundaryIds, elementIds;
224  reader.read(factory, filename, boundaryIds, elementIds);
225 
226  using HasInsertionIndexEntity
227  = Dune::Std::is_detected<SupportsInsertionIndex, GridType, typename GridType::template Codim<0>::Entity>;
228  using HasInsertionIndexIntersection
229  = Dune::Std::is_detected<SupportsInsertionIndex, GridType, typename GridType::LeafIntersection>;
230 
231  auto gridPtr = factory.createGrid();
232  if ((boundaryIds.empty() && elementIds.empty()) ||
233  (!HasInsertionIndexEntity::value && !HasInsertionIndexIntersection::value))
234  return std::unique_ptr<GridType>(std::move(gridPtr));
235 
236  // map boundaryIds and elementIds read from file to grid indexing.
237  if (!boundaryIds.empty() && HasInsertionIndexIntersection::value)
238  boundaryIds_.resize(gridPtr->numBoundarySegments());
239  if (!elementIds.empty() && HasInsertionIndexEntity::value)
240  elementIds_.resize(gridPtr->size(0));
241 
242  auto const& indexSet = gridPtr->leafIndexSet();
243  for (auto const& e : elements(gridPtr->leafGridView())) {
244  if constexpr(HasInsertionIndexEntity::value) {
245  if (!elementIds.empty())
246  elementIds_[indexSet.index(e)] = elementIds[factory.insertionIndex(e)];
247  }
248 
249  if (boundaryIds.empty() || !e.hasBoundaryIntersections())
250  continue;
251 
252  for (auto const& it : intersections(gridPtr->leafGridView(), e)) {
253  if constexpr(HasInsertionIndexIntersection::value) {
254  if (it.boundary())
255  boundaryIds_[it.boundarySegmentIndex()]
256  = boundaryIds[factory.insertionIndex(it)];
257  }
258  }
259  }
260 
261  return std::unique_ptr<GridType>(std::move(gridPtr));
262  }
263 
264  // fallback if GmshReader cannot be used
265  template <class GridType = HostGrid>
266  std::unique_ptr<GridType> read_gmsh_file(std::string const&, Dune::PriorityTag<0>) const
267  {
268  error_exit("Gmsh reader not supported for this grid type!");
269  return {};
270  }
271 
272  // read from Alberta file
273 
274 #if HAVE_ALBERTA
275  template <class GridType>
276  using IsAlbertaGrid = decltype(std::declval<GridType>().alberta2dune(0,0));
277 
278  // construct the albertagrid directly from a filename
279  template <class GridType = HostGrid,
280  REQUIRES(GridType::dimensionworld == DIM_OF_WORLD),
281  REQUIRES(Dune::Std::is_detected<IsAlbertaGrid, GridType>::value)>
282  std::unique_ptr<GridType> read_alberta_file(std::string const& filename, Dune::PriorityTag<3>) const
283  {
284  return std::make_unique<GridType>(filename);
285  }
286 
287  // use a gridfactory and the generic AlbertaReader
288  template <class GridType = HostGrid,
289  REQUIRES(GridType::dimensionworld == DIM_OF_WORLD)>
290  std::unique_ptr<GridType> read_alberta_file(std::string const& filename, Dune::PriorityTag<2>) const
291  {
292  Dune::GridFactory<GridType> factory;
293  if (Environment::mpiRank() == 0) {
294  Dune::AlbertaReader<GridType> reader;
295  reader.readGrid(filename, factory);
296  }
297  return std::unique_ptr<GridType>{factory.createGrid()};
298  }
299 
300  // error if WORLDDIM not the same as Grid::dimensionworld
301  template <class GridType = HostGrid,
302  REQUIRES(GridType::dimensionworld != DIM_OF_WORLD)>
303  std::unique_ptr<GridType> read_alberta_file(std::string const& filename, Dune::PriorityTag<1>) const
304  {
305  error_exit("AlbertaGrid (and AlbertaReader) require WORLDDIM == Grid::dimensionworld. Change the cmake parameters!");
306  return {};
307  }
308 #else
309  // fallback if no Alberta is available
310  template <class GridType>
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  return {};
315  }
316 #endif
317 
318 
319 #if HAVE_ALBERTA
320  // albertagrid -> simplex
321  template <class GridType = HostGrid,
322  REQUIRES(Dune::Std::is_detected<IsAlbertaGrid, GridType>::value)>
323  std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<3>) const
324  {
325  return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells)
326  {
327  return MacroGridFactory<GridType>::createSimplexGrid(lower, upper, numCells);
328  });
329  }
330 #endif
331 
332  // yasp grid -> cube
333  template <class GridType = HostGrid,
334  class = typename GridType::YGrid>
335  std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<2>) const
336  {
337  int overlap = 0;
338  Parameters::get(name_ + "->overlap", overlap);
339  std::string periodic(dimension, '0');
340  Parameters::get(name_ + "->periodic", periodic); // e.g. 000 01 111
341 
342  return create_yaspgrid(Types<GridType>{}, overlap, std::bitset<dimension>(periodic));
343  }
344 
345  template <int dim, class ct>
346  std::unique_ptr<HostGrid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::EquidistantCoordinates<ct,dim>>>,
347  int overlap, std::bitset<dimension> const& per) const
348  {
349  return create_structured_grid<int>([&](auto&& /*lower*/, auto&& upper, std::array<int,dimension> const& numCells)
350  {
351  return std::make_unique<HostGrid>(upper, numCells, per, overlap);
352  });
353  }
354 
355  template <int dim, class ct>
356  std::unique_ptr<HostGrid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::EquidistantOffsetCoordinates<ct,dim>>>,
357  int overlap, std::bitset<dimension> const& per) const
358  {
359  return create_structured_grid<int>([&](auto&& lower, auto&& upper, std::array<int,dimension> const& numCells)
360  {
361  return std::make_unique<HostGrid>(lower, upper, numCells, per, overlap);
362  });
363  }
364 
365  template <int dim, class ct>
366  std::unique_ptr<HostGrid> create_yaspgrid(Types<Dune::YaspGrid<dim,Dune::TensorProductCoordinates<ct,dim>>>,
367  int, std::bitset<dimension> const&) const
368  {
369  error_exit("MeshCreator cannot create YaspGrid with TensorProductCoordinates.");
370  return {};
371  }
372 
373 
374 #if HAVE_DUNE_SPGRID
375  // spgrid -> cube
376  template <class GridType = HostGrid,
377  class = typename GridType::ReferenceCube,
378  class = typename GridType::MultiIndex>
379  std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<1>) const
380  {
381  return create_structured_grid<int>([](auto&& lower, auto&& upper, std::array<int,dimension> const& numCells)
382  {
383  typename GridType::MultiIndex multiIndex(numCells);
384  return std::make_unique<GridType>(lower, upper, multiIndex);
385  });
386  }
387 #endif
388 
389  // final fallback
390  template <class GridType = HostGrid>
391  std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<0>) const
392  {
393  error_exit("Don't know how to create the grid.");
394  return {};
395  }
396 
397  private:
398  std::string name_;
399  mutable std::vector<int> boundaryIds_;
400  mutable std::vector<int> elementIds_;
401  };
402 
403 
404 } // end namespace AMDiS
constexpr bool MultiIndex
A multi-index type.
Definition: Concepts.hpp:150
path extension() const
Returns the file extension path component.
Definition: Filesystem.cpp:83
Definition: AdaptiveGrid.hpp:373
A variadic type list.
Definition: TypeTraits.hpp:88
void error_exit(std::string const &str, Args &&... args)
print a message and exit
Definition: Output.hpp:142
typename Impl::AdaptiveGridImpl< HostGrid >::type AdaptiveGrid_t
Definition: AdaptiveGrid.hpp:367
Contains all classes needed for solving linear and non linear equation systems.
Definition: AdaptBase.hpp:6
std::unique_ptr< HostGrid > create_simplex_grid() const
Create a structured simplex grid.
Definition: MeshCreator.hpp:131
std::vector< int > const & elementIds() const
Return the filled vector of element ids. NOTE: not all creators support reading this.
Definition: MeshCreator.hpp:146
A creator class for dune grids.
Definition: MeshCreator.hpp:52
static std::shared_ptr< Grid > create(std::string name)
Static create mthod. See create()
Definition: MeshCreator.hpp:71
std::unique_ptr< HostGrid > create_cube_grid() const
Create a structured cube grid.
Definition: MeshCreator.hpp:122
Definition: Filesystem.hpp:11
MeshCreator(std::string const &name)
Construct a new MeshCreator.
Definition: MeshCreator.hpp:66
std::unique_ptr< Grid > create() const
Create a new grid by inspecting the intifile parameter group [meshName]
Definition: MeshCreator.hpp:88
std::vector< int > const & boundaryIds() const
Return the filled vector of boundary ids. NOTE: not all creators support reading this.
Definition: MeshCreator.hpp:140
Definition: MacroGridFactory.hpp:105