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 #include +#include + using namespace AMDiS; using Grid = Dune::AlbertaGrid; 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 + void setBoundaryIds(std::vector 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_; 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 -#include -#include - -#include -#include - -#include -#include -#include - -#include -#include - -#include -#include -#include - -namespace AMDiS -{ - namespace tag - { - struct albertagrid {}; - struct uggrid {}; - struct yaspgrid {}; - struct unknowngrid {}; - - } // end namespace tag - - namespace Impl - { - template - struct MeshTag - { - using type = tag::unknowngrid; - }; - - // specialization for some grid types from DUNE -#if HAVE_ALBERTA - template - struct MeshTag> - { - using type = tag::albertagrid; - }; -#endif - -#if HAVE_UG - template - struct MeshTag> - { - using type = tag::uggrid; - }; -#endif - - template - struct MeshTag> - { - using type = tag::yaspgrid; - }; - - } // end namespace Impl - - template - using MeshTag_t = typename Impl::MeshTag::type; - - - /// A creator class for meshes. Each mesh needs different way of initialization - template - struct MeshCreator - { - static std::unique_ptr create(std::string const& meshName) - { - error_exit("Creator not yet implemented for this mesh type."); - return {}; - } - }; - -#if HAVE_ALBERTA - template - struct MeshCreator> - { - using Grid = Dune::AlbertaGrid; - - static std::unique_ptr 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(macro_filename); - } - }; -#endif - - -#if HAVE_UG - template - struct MeshCreator> - { - using Grid = Dune::UGGrid; - - static std::unique_ptr 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 reader; - return std::unique_ptr{reader.read(filename)}; - } -// #if HAVE_ALBERTA -// else if (ext == ".1d" || ext == ".2d" || ext == ".3d") { -// Dune::Hybrid::ifElse(bool_t{}, -// [&](auto id) { -// Dune::GridFactory factory; -// Dune::AlbertaReader reader; -// id(reader).readGrid(filename, id(factory)); -// return std::unique_ptr{factory.createGrid()}; -// }); - -// warning("WORLDDIM must be given in Alberta flags."); -// } -// #endif - - error_exit("No way to construct UG-Grid found"); - return {}; - } - }; -#endif - - template - struct MeshCreator>> - { - using Grid = Dune::YaspGrid>; - - static std::unique_ptr create(std::string const& meshName) - { - Dune::FieldVector L; L = 1.0; // extension of the domain - Parameters::get(meshName + "->dimension", L); - - auto s = Dune::filledArray(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(L, s); - } - }; - - - template - struct MeshCreator>> - { - using Grid = Dune::YaspGrid>; - - static std::unique_ptr create(std::string const& meshName) - { - Dune::FieldVector lowerleft; lowerleft = 0.0; // Lower left corner of the domain - Dune::FieldVector 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(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(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 +#include +#include +#include + +#include +#include +#include + +#if HAVE_ALBERTA +#include +#endif +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace AMDiS +{ + /// A creator class for dune grids. + template + 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 create() const + { + auto filename = Parameters::get(name_ + "->macro file name"); + auto structured = Parameters::get(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(Dune::PriorityTag<42>{}); + } + } + + /// Create a structured cube grid + std::unique_ptr create_cube_grid() const + { + return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells) + { + return Dune::StructuredGridFactory::createCubeGrid(lower, upper, numCells); + }); + } + + /// Create a structured simplex grid + std::unique_ptr create_simplex_grid() const + { + return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells) + { + return MacroGridFactory::createSimplexGrid(lower, upper, numCells); + }); + } + + /// Return the filled vector of boundary ids. NOTE: not all creators support reading this. + std::vector const& boundaryIds() const + { + return boundaryIds_; + } + + /// Return the filled vector of element ids. NOTE: not all creators support reading this. + std::vector const& elementIds() const + { + return elementIds_; + } + + private: + // use the structured grid factory to create the grid + template + std::unique_ptr create_structured_grid(Factory factory) const + { + // Lower left corner of the domain + Dune::FieldVector lower(0); + // Upper right corner of the domain + Dune::FieldVector upper(1); + // number of blocks in each coordinate direction + auto numCells = Dune::filledArray(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 create_unstructured_grid(std::string const& filename) const + { + filesystem::path fn(filename); + auto ext = fn.extension(); + + if (ext == ".msh") { + return read_gmsh_file(filename, Dune::PriorityTag<42>{}); + } + else if (ext == ".1d" || ext == ".2d" || ext == ".3d" || ext == ".amc") { + return read_alberta_file(filename, Dune::PriorityTag<42>{}); + } + else { + error_exit("Can not read grid file. Unknown file extension."); + return {}; + } + } + + template + using SupportsGmshReader = decltype(std::declval>().insertBoundarySegment( + std::declval>(), + std::declval >>()) ); + + // use GmshReader if GridFactory supports insertBoundarySegments + template ::value)> + std::unique_ptr read_gmsh_file(std::string const& filename, Dune::PriorityTag<2>) const + { + Dune::GmshReader reader; + return std::unique_ptr{reader.read(filename, boundaryIds_, elementIds_)}; + } + + // fallback if GmshReader cannot be used + template + std::unique_ptr 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 + using IsAlbertaGrid = decltype(std::declval().alberta2dune(0,0)); + + // construct the albertagrid directly from a filename + template ::value)> + std::unique_ptr read_alberta_file(std::string const& filename, Dune::PriorityTag<3>) const + { + return std::make_unique(filename); + } + + // use a gridfactory and the generic AlbertaReader + template + std::unique_ptr read_alberta_file(std::string const& filename, Dune::PriorityTag<2>) const + { + Dune::GridFactory factory; + Dune::AlbertaReader reader; + reader.readGrid(filename, factory); + return std::unique_ptr{factory.createGrid()}; + } + + // error if WORLDDIM not the same as Grid::dimensionworld + template + std::unique_ptr 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 + std::unique_ptr 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 ::value)> + std::unique_ptr create_by_gridtype(Dune::PriorityTag<3>) const + { + return create_simplex_grid(); + } +#endif + + // yasp grid -> cube + template + std::unique_ptr create_by_gridtype(Dune::PriorityTag<2>) const + { + return create_cube_grid(); + } + +#if HAVE_DUNE_SPGRID + // spgrid -> cube + template + std::unique_ptr create_by_gridtype(Dune::PriorityTag<1>) const + { + return create_structured_grid([](auto&& lower, auto&& upper, auto&& numCells) + { + std::array cells; + std::copy(std::begin(numCells), std::end(numCells), std::begin(cells)); + typename GridType::MultiIndex multiIndex(cells); + return std::make_unique(lower, upper, multiIndex); + }); + } +#endif + + // final fallback + template + std::unique_ptr 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 boundaryIds_; + mutable std::vector 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 #include #include -#include +#include #include #include #include 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 void ProblemStat::createGrid() { Parameters::get(name_ + "->mesh", gridName_); - grid_ = MeshCreator::create(gridName_); + MeshCreator creator(gridName_); + grid_ = creator.create(); boundaryManager_ = std::make_shared>(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::true_type {}; + + template + struct IsDefined + : std::false_type {}; + + template + struct IsDefined::value && + !std::is_pointer::value && + (sizeof(T) > 0)> > + : std::true_type {}; + } // end namespace Traits namespace Concepts @@ -117,6 +128,9 @@ namespace AMDiS template constexpr bool MultiIndex = models; + template + constexpr bool Defined = Traits::IsDefined::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 +#include +#include + +#if ! DUNE_VERSION_GT(DUNE_GRID,2,6) + #include +#endif + +namespace Dune +{ + template + struct CubedWrapper : public GridType {}; + + + template + class GridFactory> + : public GridFactoryInterface + { + using Self = GridFactory; + using ctype = typename GridType::ctype; + + enum { dim = GridType::dimension }; + enum { dimworld = GridType::dimensionworld }; + + public: + template = 0> + GridFactory (Args&&... args) + : factory_(std::make_shared>(std::forward(args)...)) + {} + + GridFactory (GridFactory& factory) + : factory_(stackobject_to_shared_ptr(factory)) + {} + + /// \brief Insert a vertex into the coarse grid + void insertVertex (const FieldVector& 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& vertices) override + { + // triangulation of reference cube + static const auto reference_cubes = std::make_tuple( + std::array, 1>{std::array{0,1}}, + std::array, 2>{std::array{3,0,1}, std::array{0,3,2}}, + std::array, 6>{std::array{0,7,3,1}, std::array{0,7,5,1}, + std::array{0,7,5,4}, std::array{0,7,3,2}, + std::array{0,7,6,2}, std::array{0,7,6,4}} ); + + assert(type == GeometryTypes::cube(dim)); + + auto const& simplices = std::get(reference_cubes); + thread_local std::vector 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& vertices) override + { + factory_->insertBoundarySegment(vertices); + } + + /// \brief Finalize grid creation and hand over the grid +#if DUNE_VERSION_GT(DUNE_GRID,2,6) + ToUniquePtr createGrid () override +#else + GridType* createGrid () override +#endif + { + return factory_->createGrid(); + } + + private: + std::shared_ptr> factory_; + }; + +} // end namespace Dune + + +namespace AMDiS +{ + template + 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 createSimplexGrid (Dune::GridFactory& originalFactory, + const Dune::FieldVector& lowerLeft, + const Dune::FieldVector& upperRight, + const std::array& numElements) + { + Dune::GridFactory> factory(originalFactory); + Dune::StructuredGridFactory>::createCubeGrid(factory, lowerLeft, upperRight, numElements); + return std::unique_ptr(factory.createGrid()); + } +#endif + + /// \brief Create a structured simplex grid + static std::unique_ptr createSimplexGrid (const Dune::FieldVector& lowerLeft, + const Dune::FieldVector& upperRight, + const std::array& numElements) + { + Dune::GridFactory> factory; +#if DUNE_VERSION_GT(DUNE_GRID,2,6) + Dune::StructuredGridFactory>::createCubeGrid(factory, lowerLeft, upperRight, numElements); +#else + // fallback implementation using temporary YaspGrid + using TempGrid = Dune::YaspGrid>; + auto grid = Dune::StructuredGridFactory::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 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(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 #include +#include + #include "Tests.hpp" using namespace AMDiS;