diff --git a/examples/cahn_hilliard.cc b/examples/cahn_hilliard.cc
index a26e9ccf229ace431993e1284318162341a0eb57..fe37bdd7cfc013741b564b23efe12afe8e399df7 100644
--- a/examples/cahn_hilliard.cc
+++ b/examples/cahn_hilliard.cc
@@ -6,6 +6,8 @@
#include <amdis/GridFunctions.hpp>
#include <amdis/Marker.hpp>
+#include <dune/grid/albertagrid.hh>
+
using namespace AMDiS;
using Grid = Dune::AlbertaGrid<GRIDDIM, WORLDDIM>;
diff --git a/examples/init/heat.dat.2d b/examples/init/heat.dat.2d
index 31836835ef599d98b3335ba67d7c6e9c64980fcc..e938b21c65fa0df47152df6fb850de5525c2a96c 100644
--- a/examples/init/heat.dat.2d
+++ b/examples/init/heat.dat.2d
@@ -1,9 +1,8 @@
heatMesh->macro file name: ./macro/macro.stand.2d
heatMesh->global refinements: 4
heatMesh->min corner: 0 0
-heatMesh->max corner: 1 1
+heatMesh->max corner: 2 2
heatMesh->num cells: 2 2
-heatMesh->dimension: 2 2
heat->mesh: heatMesh
heat->names: u
diff --git a/examples/init/stokes.dat.2d b/examples/init/stokes.dat.2d
index 54542825dba336eba2af96e3b3371f548ff8bad6..ded25d60812bc6e5f3eac56f709fb43b1fa9afe3 100644
--- a/examples/init/stokes.dat.2d
+++ b/examples/init/stokes.dat.2d
@@ -1,5 +1,5 @@
stokesMesh->global refinements: 0
-stokesMesh->dimension: 1.0 1.0
+stokesMesh->max corner: 1.0 1.0
stokesMesh->num cells: 4 4
stokes->mesh: stokesMesh
diff --git a/src/amdis/BoundaryManager.hpp b/src/amdis/BoundaryManager.hpp
index 9e9191eb25e45f19a42b3c13910e905a68f90598..1b4358396b23021e0d7a6a5957a4e63fecb0d931 100644
--- a/src/amdis/BoundaryManager.hpp
+++ b/src/amdis/BoundaryManager.hpp
@@ -154,6 +154,13 @@ namespace AMDiS
}
}
+ template <class I>
+ void setBoundaryIds(std::vector<I> const& ids)
+ {
+ test_exit(ids.size() == boundaryIds_.size(), "Number of boundary IDs does not match!");
+ std::copy(ids.begin(), ids.end(), boundaryIds_.begin());
+ }
+
private:
std::shared_ptr<Grid> grid_;
using Super::boundaryIds_;
diff --git a/src/amdis/CMakeLists.txt b/src/amdis/CMakeLists.txt
index ceaf0b69d0750e114d7f43c85b0d91e509c49f36..1237bbc906cd4a0c5806d782675d552daeb8afd2 100644
--- a/src/amdis/CMakeLists.txt
+++ b/src/amdis/CMakeLists.txt
@@ -50,7 +50,7 @@ install(FILES
LocalOperators.hpp
Marker.hpp
Marker.inc.hpp
- Mesh.hpp
+ MeshCreator.hpp
Operations.hpp
OperatorList.hpp
Output.hpp
diff --git a/src/amdis/Mesh.hpp b/src/amdis/Mesh.hpp
deleted file mode 100644
index d4524bf3b98ca98eaea1e30d50f649c1cf134e77..0000000000000000000000000000000000000000
--- a/src/amdis/Mesh.hpp
+++ /dev/null
@@ -1,183 +0,0 @@
-#pragma once
-
-#include <array>
-#include <memory>
-#include <string>
-
-#include <dune/common/filledarray.hh>
-#include <dune/common/fvector.hh>
-
-#include <dune/grid/albertagrid.hh>
-#include <dune/grid/uggrid.hh>
-#include <dune/grid/yaspgrid.hh>
-
-#include <dune/grid/albertagrid/albertareader.hh>
-#include <dune/grid/io/file/gmshreader.hh>
-
-#include <amdis/Initfile.hpp>
-#include <amdis/Output.hpp>
-#include <amdis/common/Filesystem.hpp>
-
-namespace AMDiS
-{
- namespace tag
- {
- struct albertagrid {};
- struct uggrid {};
- struct yaspgrid {};
- struct unknowngrid {};
-
- } // end namespace tag
-
- namespace Impl
- {
- template <class Grid>
- struct MeshTag
- {
- using type = tag::unknowngrid;
- };
-
- // specialization for some grid types from DUNE
-#if HAVE_ALBERTA
- template <int dim, int dimworld>
- struct MeshTag<Dune::AlbertaGrid<dim, dimworld>>
- {
- using type = tag::albertagrid;
- };
-#endif
-
-#if HAVE_UG
- template <int dim>
- struct MeshTag<Dune::UGGrid<dim>>
- {
- using type = tag::uggrid;
- };
-#endif
-
- template <int dim, class Coordinates>
- struct MeshTag<Dune::YaspGrid<dim, Coordinates>>
- {
- using type = tag::yaspgrid;
- };
-
- } // end namespace Impl
-
- template <class Grid>
- using MeshTag_t = typename Impl::MeshTag<Grid>::type;
-
-
- /// A creator class for meshes. Each mesh needs different way of initialization
- template <class Grid>
- struct MeshCreator
- {
- static std::unique_ptr<Grid> create(std::string const& meshName)
- {
- error_exit("Creator not yet implemented for this mesh type.");
- return {};
- }
- };
-
-#if HAVE_ALBERTA
- template <int dim, int dimworld>
- struct MeshCreator<Dune::AlbertaGrid<dim, dimworld>>
- {
- using Grid = Dune::AlbertaGrid<dim, dimworld>;
-
- static std::unique_ptr<Grid> create(std::string const& meshName)
- {
- std::string macro_filename = "";
- Parameters::get(meshName + "->macro file name", macro_filename);
-
- // TODO: if filename_extension is ".2d" or ".3d" read it directly from file
- // otherwise use a factory method
-
- return std::make_unique<Grid>(macro_filename);
- }
- };
-#endif
-
-
-#if HAVE_UG
- template <int dim>
- struct MeshCreator<Dune::UGGrid<dim>>
- {
- using Grid = Dune::UGGrid<dim>;
-
- static std::unique_ptr<Grid> create(std::string const& meshName)
- {
-
- std::string filename = "";
- Parameters::get(meshName + "->macro file name", filename);
-
- test_exit(!filename.empty(),
- "Construction of UGGrid without filename not yet implemented!");
-
- filesystem::path fn(filename);
- auto ext = fn.extension();
-
- if (ext == ".msh") {
- Dune::GmshReader<Grid> reader;
- return std::unique_ptr<Grid>{reader.read(filename)};
- }
-// #if HAVE_ALBERTA
-// else if (ext == ".1d" || ext == ".2d" || ext == ".3d") {
-// Dune::Hybrid::ifElse(bool_t<ALBERTA_DIM == Grid::dimensionworld>{},
-// [&](auto id) {
-// Dune::GridFactory<Grid> factory;
-// Dune::AlbertaReader<Grid> reader;
-// id(reader).readGrid(filename, id(factory));
-// return std::unique_ptr<Grid>{factory.createGrid()};
-// });
-
-// warning("WORLDDIM must be given in Alberta flags.");
-// }
-// #endif
-
- error_exit("No way to construct UG-Grid found");
- return {};
- }
- };
-#endif
-
- template <int dim, class T>
- struct MeshCreator<Dune::YaspGrid<dim, Dune::EquidistantCoordinates<T,dim>>>
- {
- using Grid = Dune::YaspGrid<dim, Dune::EquidistantCoordinates<T,dim>>;
-
- static std::unique_ptr<Grid> create(std::string const& meshName)
- {
- Dune::FieldVector<T, dim> L; L = 1.0; // extension of the domain
- Parameters::get(meshName + "->dimension", L);
-
- auto s = Dune::filledArray<std::size_t(dim)>(1); // number of cells on coarse mesh in each direction
- Parameters::get(meshName + "->num cells", s);
-
- // TODO: add more parameters for yasp-grid (see constructor)
-
- return std::make_unique<Grid>(L, s);
- }
- };
-
-
- template <int dim, class T>
- struct MeshCreator<Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<T, dim>>>
- {
- using Grid = Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<T, dim>>;
-
- static std::unique_ptr<Grid> create(std::string const& meshName)
- {
- Dune::FieldVector<T, dim> lowerleft; lowerleft = 0.0; // Lower left corner of the domain
- Dune::FieldVector<T, dim> upperright; upperright = 1.0; // Upper right corner of the domain
- Parameters::get(meshName + "->min corner", lowerleft);
- Parameters::get(meshName + "->max corner", upperright);
-
- auto s = Dune::filledArray<std::size_t(dim)>(1); // number of cells on coarse mesh in each direction
- Parameters::get(meshName + "->num cells", s);
-
- // TODO: add more parameters for yasp-grid (see constructor)
-
- return std::make_unique<Grid>(lowerleft, upperright, s);
- }
- };
-
-} // end namespace AMDiS
diff --git a/src/amdis/MeshCreator.hpp b/src/amdis/MeshCreator.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2dc00df3b33a6e4581ddf7419f7134cc739cc36a
--- /dev/null
+++ b/src/amdis/MeshCreator.hpp
@@ -0,0 +1,261 @@
+#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>
+
+#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>
+
+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;
+
+ /// Construct a new MeshCreator
+ /**
+ * \param name The name of the mesh used in the initifile
+ **/
+ 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
+ *
+ * Otherwise tries to create a grid depending on the grid type.
+ **/
+ std::unique_ptr<Grid> create() const
+ {
+ auto filename = Parameters::get<std::string>(name_ + "->macro file name");
+ auto structured = Parameters::get<std::string>(name_ + "->structured");
+ if (filename) {
+ // read a macro file
+ return create_unstructured_grid(filename.value());
+ } else if (structured) {
+ // create structured grids
+ if (structured.value() == "cube") {
+ return create_cube_grid();
+ } else if (structured && structured.value() == "simplex") {
+ return create_simplex_grid();
+ } else {
+ error_exit("Unknown structured grid type. Must be either `cube` or `simplex` in parameter [meshName]->structured.");
+ return {};
+ }
+ } else {
+ // decide by inspecting the grid type how to create the grid
+ return create_by_gridtype<Grid>(Dune::PriorityTag<42>{});
+ }
+ }
+
+ /// Create a structured cube grid
+ std::unique_ptr<Grid> create_cube_grid() const
+ {
+ return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells)
+ {
+ return Dune::StructuredGridFactory<Grid>::createCubeGrid(lower, upper, numCells);
+ });
+ }
+
+ /// Create a structured simplex grid
+ std::unique_ptr<Grid> create_simplex_grid() const
+ {
+ return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells)
+ {
+ return MacroGridFactory<Grid>::createSimplexGrid(lower, upper, numCells);
+ });
+ }
+
+ /// 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:
+ // use the structured grid factory to create the grid
+ template <class Factory>
+ std::unique_ptr<Grid> create_structured_grid(Factory factory) const
+ {
+ // 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
+ auto numCells = Dune::filledArray<std::size_t(dimension),unsigned int>(1);
+
+ 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
+ std::unique_ptr<Grid> create_unstructured_grid(std::string const& filename) const
+ {
+ filesystem::path fn(filename);
+ auto ext = fn.extension();
+
+ if (ext == ".msh") {
+ return read_gmsh_file<Grid>(filename, Dune::PriorityTag<42>{});
+ }
+ else if (ext == ".1d" || ext == ".2d" || ext == ".3d" || ext == ".amc") {
+ return read_alberta_file<Grid>(filename, Dune::PriorityTag<42>{});
+ }
+ else {
+ error_exit("Can not read grid file. Unknown file extension.");
+ 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
+ template <class GridType = Grid,
+ REQUIRES(Dune::Std::is_detected<SupportsGmshReader, GridType>::value)>
+ std::unique_ptr<GridType> read_gmsh_file(std::string const& filename, Dune::PriorityTag<2>) const
+ {
+ Dune::GmshReader<GridType> reader;
+ return std::unique_ptr<GridType>{reader.read(filename, boundaryIds_, elementIds_)};
+ }
+
+ // fallback if GmshReader cannot be used
+ template <class GridType = Grid>
+ 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
+ template <class GridType = Grid,
+ 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
+ template <class GridType = Grid,
+ 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;
+ Dune::AlbertaReader<GridType> reader;
+ reader.readGrid(filename, factory);
+ return std::unique_ptr<GridType>{factory.createGrid()};
+ }
+
+ // error if WORLDDIM not the same as Grid::dimensionworld
+ template <class GridType = Grid,
+ REQUIRES(GridType::dimensionworld != DIM_OF_WORLD)>
+ std::unique_ptr<GridType> read_alberta_file(std::string const& filename, Dune::PriorityTag<1>) const
+ {
+ error_exit("Alberta (and AlbertaReader) require WORLDDIM == Grid::dimensionworld. Change the cmake parameters!");
+ 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
+ {
+ error_exit("Alberta (and AlbertaReader) not available. Set AlbertaFlags to your executable in cmake!");
+ return {};
+ }
+#endif
+
+
+#if HAVE_ALBERTA
+ // albertagrid -> simplex
+ template <class GridType = Grid,
+ REQUIRES(Dune::Std::is_detected<IsAlbertaGrid, GridType>::value)>
+ std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<3>) const
+ {
+ return create_simplex_grid();
+ }
+#endif
+
+ // yasp grid -> cube
+ template <class GridType = Grid,
+ class = typename GridType::YGrid>
+ std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<2>) const
+ {
+ return create_cube_grid();
+ }
+
+#if HAVE_DUNE_SPGRID
+ // spgrid -> cube
+ template <class GridType = Grid,
+ class = typename GridType::ReferenceCube,
+ class = typename GridType::MultiIndex>
+ std::unique_ptr<GridType> create_by_gridtype(Dune::PriorityTag<1>) const
+ {
+ return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells)
+ {
+ std::array<int,dimworld> cells;
+ std::copy(std::begin(numCells), std::end(numCells), std::begin(cells));
+ typename GridType::MultiIndex multiIndex(cells);
+ return std::make_unique<GridType>(lower, upper, multiIndex);
+ });
+ }
+#endif
+
+ // final fallback
+ template <class GridType = Grid>
+ 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
diff --git a/src/amdis/ProblemStat.hpp b/src/amdis/ProblemStat.hpp
index 005cf68eeb8afa62ed9e4213293fe8222d377395..faf6a34219f968efe7d481dfd0001d3695f97a1e 100644
--- a/src/amdis/ProblemStat.hpp
+++ b/src/amdis/ProblemStat.hpp
@@ -22,7 +22,7 @@
#include <amdis/LinearAlgebra.hpp>
#include <amdis/OperatorList.hpp>
#include <amdis/Marker.hpp>
-#include <amdis/Mesh.hpp>
+#include <amdis/MeshCreator.hpp>
#include <amdis/PeriodicBC.hpp>
#include <amdis/ProblemStatBase.hpp>
#include <amdis/ProblemStatTraits.hpp>
diff --git a/src/amdis/ProblemStat.inc.hpp b/src/amdis/ProblemStat.inc.hpp
index 0a942639d2b8038814a3c898a271384b12087f2d..34d552a6fe77401f218dbcf0df33163d3f6a4666 100644
--- a/src/amdis/ProblemStat.inc.hpp
+++ b/src/amdis/ProblemStat.inc.hpp
@@ -126,8 +126,11 @@ template <class Traits>
void ProblemStat<Traits>::createGrid()
{
Parameters::get(name_ + "->mesh", gridName_);
- grid_ = MeshCreator<Grid>::create(gridName_);
+ MeshCreator<Grid> creator(gridName_);
+ grid_ = creator.create();
boundaryManager_ = std::make_shared<BoundaryManager<Grid>>(grid_);
+ if (!creator.boundaryIds().empty())
+ boundaryManager_->setBoundaryIds(creator.boundaryIds());
msg("Create grid:");
msg("#elements = {}" , grid_->size(0));
diff --git a/src/amdis/common/Concepts.hpp b/src/amdis/common/Concepts.hpp
index fc4aba7ad538be05eb359d8f4c6a73029e0bb8df..3a43c784c6badb7a470ff8b6b5a2577e4a4ad2e0 100644
--- a/src/amdis/common/Concepts.hpp
+++ b/src/amdis/common/Concepts.hpp
@@ -53,6 +53,17 @@ namespace AMDiS
struct IsReferenceWrapper<std::reference_wrapper<T>>
: std::true_type {};
+
+ template <class, class = void>
+ struct IsDefined
+ : std::false_type {};
+
+ template <class T>
+ struct IsDefined<T, std::enable_if_t<std::is_object<T>::value &&
+ !std::is_pointer<T>::value &&
+ (sizeof(T) > 0)> >
+ : std::true_type {};
+
} // end namespace Traits
namespace Concepts
@@ -117,6 +128,9 @@ namespace AMDiS
template <class MI>
constexpr bool MultiIndex = models<Definition::MultiIndex(MI)>;
+ template <class T>
+ constexpr bool Defined = Traits::IsDefined<T>::value;
+
/** @} **/
} // end namespace Concepts
diff --git a/src/amdis/utility/CMakeLists.txt b/src/amdis/utility/CMakeLists.txt
index f9300e39e6bb22f3525dfc4a04c93363f0e6e03c..870230af8d6cc5e42061d2bffb6e135918017f11 100644
--- a/src/amdis/utility/CMakeLists.txt
+++ b/src/amdis/utility/CMakeLists.txt
@@ -2,5 +2,6 @@ install(FILES
AssembleOperators.hpp
LocalBasisCache.hpp
LocalToGlobalAdapter.hpp
+ MacroGridFactory.hpp
QuadratureFactory.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/utility)
diff --git a/src/amdis/utility/MacroGridFactory.hpp b/src/amdis/utility/MacroGridFactory.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..26b0eae6acb09391bd5cf01564d04417f10e98b1
--- /dev/null
+++ b/src/amdis/utility/MacroGridFactory.hpp
@@ -0,0 +1,150 @@
+#pragma once
+
+#include <dune/common/typeutilities.hh>
+#include <dune/common/version.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+
+#if ! DUNE_VERSION_GT(DUNE_GRID,2,6)
+ #include <dune/grid/yaspgrid.hh>
+#endif
+
+namespace Dune
+{
+ template <class GridType>
+ struct CubedWrapper : public GridType {};
+
+
+ template <class GridType>
+ class GridFactory<CubedWrapper<GridType>>
+ : public GridFactoryInterface<GridType>
+ {
+ using Self = GridFactory;
+ using ctype = typename GridType::ctype;
+
+ enum { dim = GridType::dimension };
+ enum { dimworld = GridType::dimensionworld };
+
+ public:
+ template <class... Args, disableCopyMove<Self, Args...> = 0>
+ GridFactory (Args&&... args)
+ : factory_(std::make_shared<GridFactory<GridType>>(std::forward<Args>(args)...))
+ {}
+
+ GridFactory (GridFactory<GridType>& factory)
+ : factory_(stackobject_to_shared_ptr(factory))
+ {}
+
+ /// \brief Insert a vertex into the coarse grid
+ void insertVertex (const FieldVector<ctype,dimworld>& pos) override
+ {
+ factory_->insertVertex(pos);
+ }
+
+ /// \brief Insert simplex elements into the coarse grid
+ /**
+ * Creates a simplex subdividion of the cube element.
+ *
+ * \param type The GeometryType of the box grid
+ * \param vertices Indices of the cube corners
+ **/
+ void insertElement (const GeometryType& type,
+ const std::vector<unsigned int>& vertices) override
+ {
+ // triangulation of reference cube
+ static const auto reference_cubes = std::make_tuple(
+ std::array<std::array<int,2>, 1>{std::array<int,2>{0,1}},
+ std::array<std::array<int,3>, 2>{std::array<int,3>{3,0,1}, std::array<int,3>{0,3,2}},
+ std::array<std::array<int,4>, 6>{std::array<int,4>{0,7,3,1}, std::array<int,4>{0,7,5,1},
+ std::array<int,4>{0,7,5,4}, std::array<int,4>{0,7,3,2},
+ std::array<int,4>{0,7,6,2}, std::array<int,4>{0,7,6,4}} );
+
+ assert(type == GeometryTypes::cube(dim));
+
+ auto const& simplices = std::get<dim-1>(reference_cubes);
+ thread_local std::vector<unsigned int> corners(dim+1);
+ for (auto const& simplex : simplices) {
+ for (std::size_t i = 0; i < simplex.size(); ++i)
+ corners[i] = vertices[simplex[i]];
+
+ factory_->insertElement(GeometryTypes::simplex(dim), corners);
+ }
+ }
+
+ /// \brief insert a boundary segment
+ // TODO: maybe split boundary segment in simplices
+ void insertBoundarySegment (const std::vector<unsigned int>& vertices) override
+ {
+ factory_->insertBoundarySegment(vertices);
+ }
+
+ /// \brief Finalize grid creation and hand over the grid
+#if DUNE_VERSION_GT(DUNE_GRID,2,6)
+ ToUniquePtr<GridType> createGrid () override
+#else
+ GridType* createGrid () override
+#endif
+ {
+ return factory_->createGrid();
+ }
+
+ private:
+ std::shared_ptr<GridFactory<GridType>> factory_;
+ };
+
+} // end namespace Dune
+
+
+namespace AMDiS
+{
+ template <class GridType>
+ class MacroGridFactory
+ {
+ using ctype = typename GridType::ctype;
+
+ enum { dim = GridType::dimension };
+ enum { dimworld = GridType::dimensionworld };
+
+ public:
+#if DUNE_VERSION_GT(DUNE_GRID,2,6)
+ /// \brief insert structured simplex grid into grid factory
+ static std::unique_ptr<GridType> createSimplexGrid (Dune::GridFactory<GridType>& originalFactory,
+ const Dune::FieldVector<ctype,dimworld>& lowerLeft,
+ const Dune::FieldVector<ctype,dimworld>& upperRight,
+ const std::array<unsigned int,dim>& numElements)
+ {
+ Dune::GridFactory<Dune::CubedWrapper<GridType>> factory(originalFactory);
+ Dune::StructuredGridFactory<Dune::CubedWrapper<GridType>>::createCubeGrid(factory, lowerLeft, upperRight, numElements);
+ return std::unique_ptr<GridType>(factory.createGrid());
+ }
+#endif
+
+ /// \brief Create a structured simplex grid
+ static std::unique_ptr<GridType> createSimplexGrid (const Dune::FieldVector<ctype,dimworld>& lowerLeft,
+ const Dune::FieldVector<ctype,dimworld>& upperRight,
+ const std::array<unsigned int,dim>& numElements)
+ {
+ Dune::GridFactory<Dune::CubedWrapper<GridType>> factory;
+#if DUNE_VERSION_GT(DUNE_GRID,2,6)
+ Dune::StructuredGridFactory<Dune::CubedWrapper<GridType>>::createCubeGrid(factory, lowerLeft, upperRight, numElements);
+#else
+ // fallback implementation using temporary YaspGrid
+ using TempGrid = Dune::YaspGrid<dim,Dune::EquidistantOffsetCoordinates<double,dim>>;
+ auto grid = Dune::StructuredGridFactory<TempGrid>::createCubeGrid(lowerLeft, upperRight, numElements);
+ for (auto const& v : vertices(grid->leafGridView()))
+ factory.insertVertex(v.geometry().corner(0));
+
+ auto const& indexSet = grid->leafIndexSet();
+ for (auto const& e : elements(grid->leafGridView()))
+ {
+ thread_local std::vector<unsigned int> vertices;
+ vertices.resize(e.subEntities(dim));
+ for (unsigned int i = 0; i < e.subEntities(dim); ++i)
+ vertices[i] = indexSet.subIndex(e,i,dim);
+ factory.insertElement(Dune::GeometryTypes::cube(dim), vertices);
+ }
+#endif
+ return std::unique_ptr<GridType>(factory.createGrid());
+ }
+ };
+
+} // end namespace AMDiS
diff --git a/test/MarkerTest.cpp b/test/MarkerTest.cpp
index a120ce9a6cae6719a2e91321d3656b1cd4f5a825..6ae27b0def0ca4b167c257ce9bcf5428ebca2811 100644
--- a/test/MarkerTest.cpp
+++ b/test/MarkerTest.cpp
@@ -1,6 +1,8 @@
#include <amdis/AMDiS.hpp>
#include <amdis/ProblemStat.hpp>
+#include <dune/grid/uggrid.hh>
+
#include "Tests.hpp"
using namespace AMDiS;