Skip to content
Snippets Groups Projects
Commit 818b7a19 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Merge branch 'feature/mesh_creator' into 'master'

rewritten mesh creator to support more general creation

See merge request !47
parents bbf80751 000e35b6
No related branches found
No related tags found
1 merge request!47rewritten mesh creator to support more general creation
......@@ -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>;
......
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
......
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
......
......@@ -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_;
......
......@@ -50,7 +50,7 @@ install(FILES
LocalOperators.hpp
Marker.hpp
Marker.inc.hpp
Mesh.hpp
MeshCreator.hpp
Operations.hpp
OperatorList.hpp
Output.hpp
......
#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
#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
......@@ -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>
......
......@@ -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));
......
......@@ -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
......
......@@ -2,5 +2,6 @@ install(FILES
AssembleOperators.hpp
LocalBasisCache.hpp
LocalToGlobalAdapter.hpp
MacroGridFactory.hpp
QuadratureFactory.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/utility)
#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
#include <amdis/AMDiS.hpp>
#include <amdis/ProblemStat.hpp>
#include <dune/grid/uggrid.hh>
#include "Tests.hpp"
using namespace AMDiS;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment