Commit 4423c539 authored by Praetorius, Simon's avatar Praetorius, Simon

gridfactory cleaned up

parent 670041db
......@@ -9,8 +9,11 @@
*/
#include <memory>
#include <type_traits>
#include <vector>
#include <dune/common/hybridutilities.hh>
#include <dune/common/std/type_traits.hh>
#include <dune/grid/common/gridfactory.hh>
#include "multimesh.hh"
......@@ -41,86 +44,111 @@ namespace Dune
using GlobalCoordinate = FieldVector<ctype, dimensionworld>;
using BoundarySegment = Dune::BoundarySegment<dimension, dimensionworld>;
using ElementParametrization = VirtualFunction<LocalCoordinate,GlobalCoordinate>;
using ElementParametrization = VirtualFunction<LocalCoordinate, GlobalCoordinate>;
public:
/// \brief Constructor
/**
* \param[in] n The number of grids to construct
* \param[in] args Additional parameters passed to the constructor of the hostgridfactory
* \param[in] args Additional parameters passed to the constructor of the
* hostgridfactory
**/
template <class... Args>
explicit GridFactory (std::size_t n, Args&&... args)
: gridFactories_(n, GridFactory<HostGrid>{std::forward<Args>(args)...})
{}
{
gridFactories_.reserve(n);
for (std::size_t i = 0; i < n; ++i)
gridFactories_.emplace_back(std::forward<Args>(args)...);
}
// initialize at least 1 grid
GridFactory ()
: GridFactory{1}
{}
/// \brief insert a vertex into the macro grid
/// \brief Insert a vertex into the macro grid
/**
* \param[in] pos position of the vertex (in world coordinates)
* \param[in] pos Position of the vertex (in world coordinates)
**/
virtual void insertVertex (const GlobalCoordinate& pos)
virtual void insertVertex (const GlobalCoordinate& pos) override
{
for (auto& gridFactory : gridFactories_)
gridFactory.insertVertex(pos);
}
/// \brief insert an element into the coarse grid
/// \brief Insert an element into the coarse grid
/**
* \param[in] type GeometryType of the new element
* \param[in] vertices indices of the element vertices
* \param[in] vertices Indices of the element vertices
**/
virtual void insertElement (const GeometryType& type,
const std::vector<unsigned int>& vertices)
const std::vector<unsigned int>& vertices) override
{
for (auto& gridFactory : gridFactories_)
gridFactory.insertElement(type, vertices);
}
template <class GF>
using HasInsertElement = decltype( std::declval<GF>().insertElement(
std::declval<GeometryType>(),
std::declval<std::vector<unsigned int>>(),
std::declval<std::shared_ptr<ElementParametrization>>()) );
/// \brief Insert a parametrized element into the coarse grid
/**
* \param[in] type The GeometryType of the new element
* \param[in] vertices The vertices of the new element, using the DUNE numbering
* \param[in] elementParametrization A function prescribing the shape of this element
* \param[in] type The GeometryType of the new element
* \param[in] vertices The vertices of the new element
* \param[in] param A function prescribing the shape of this element
**/
// virtual void insertElement (const GeometryType& type,
// const std::vector<unsigned int>& vertices,
// const std::shared_ptr<ElementParametrization>& elementParametrization)
// {
// for (auto& gridFactory : gridFactories_)
// gridFactory.insertElement(type, vertices, elementParametrization);
// }
/// \brief insert a boundary segment into the macro grid
virtual void insertElement (
const GeometryType& type,
const std::vector<unsigned int>& vertices,
const std::shared_ptr<ElementParametrization>& param) override
{
Hybrid::ifElse(Std::is_detected<HasInsertElement, GridFactory<HostGrid>>{},
[&](auto id) {
for (auto& gridFactory : gridFactories_)
id(gridFactory).insertElement(type, vertices, param);
});
}
/// \brief Insert a boundary segment into the macro grid
/**
* Only influences the ordering of the boundary segments
* \param[in] vertices vertex indices of boundary face
* \param[in] vertices Vertex indices of boundary face
**/
virtual void insertBoundarySegment (const std::vector<unsigned int>& vertices) override
virtual void insertBoundarySegment (
const std::vector<unsigned int>& vertices) override
{
for (auto& gridFactory : gridFactories_)
gridFactory.insertBoundarySegment(vertices);
}
/** \brief insert a shaped boundary segment into the macro grid
template <class GF>
using HasInsertBoundarySegment = decltype( std::declval<GF>().insertBoundarySegment(
std::declval<std::vector<unsigned int>>(),
std::declval<std::shared_ptr<BoundarySegment>>()) );
/** \brief Insert a shaped boundary segment into the macro grid
*
* \param[in] vertices vertex indices of boundary face
* \param[in] boundarySegment geometric realization of shaped boundary
* \param[in] vertices Vertex indices of boundary face
* \param[in] boundarySegment Geometric realization of shaped boundary
*/
virtual void insertBoundarySegment (const std::vector<unsigned int>& vertices,
const std::shared_ptr<BoundarySegment>& boundarySegment) override
virtual void insertBoundarySegment (
const std::vector<unsigned int>& vertices,
const std::shared_ptr<BoundarySegment>& boundarySegment) override
{
for (auto& gridFactory : gridFactories_)
gridFactory.insertBoundarySegment(vertices, boundarySegment);
Hybrid::ifElse(Std::is_detected<HasInsertBoundarySegment, GridFactory<HostGrid>>{},
[&](auto id) {
for (auto& gridFactory : gridFactories_)
id(gridFactory).insertBoundarySegment(vertices, boundarySegment);
});
}
/// \brief finalize grid creation and hand over the grid
/// \brief Finalize grid creation and hand over the grid
/**
* Create n copies of the grid and store it in unique_pointers
* inside the multimesh grid.
......
......@@ -40,8 +40,8 @@ namespace Dune
template <PartitionIteratorType pitype, class GridImp>
class MultiMeshIterator<0, pitype, GridImp>
: public ForwardIteratorFacade<MultiMeshIterator<0,pitype,GridImp>,
std::vector<typename GridImp::Traits::template Codim<0>::Entity>,
std::vector<typename GridImp::Traits::template Codim<0>::Entity> >
std::vector<typename GridImp::Traits::template Codim<0>::Entity>,
std::vector<typename GridImp::Traits::template Codim<0>::Entity> >
{
private:
// LevelIterator to the equivalent entity in the host grid
......@@ -53,7 +53,7 @@ namespace Dune
struct EntityStackEntry
{
template <class Entity>
explicit EntityStackEntry(Entity&& entity)
explicit EntityStackEntry (Entity&& entity)
: it(entity.hbegin(entity.level()+1))
, end(entity.hend(entity.level()+1))
{}
......@@ -68,21 +68,18 @@ namespace Dune
using Super = std::stack<EntityStackEntry, std::vector<EntityStackEntry>>;
public:
explicit EntityStack(const GridImp* multiMesh)
explicit EntityStack (const GridImp* multiMesh)
{
c.reserve(multiMesh->maxLevel());
}
// return true of all levels >= l are finished, i.e. hierarchic iterators it == end
bool finished(std::size_t l = 0) const
// return true if all levels >= l are finished, i.e. hierarchic iterators it == end
bool finished (std::size_t l = 0) const
{
bool f = true;
while (f && l < c.size()) {
auto tmp = c[l].it;
++tmp;
f = f && tmp == c[l].end;
++l;
}
for (auto tmp = c[l].it; f && l < c.size(); ++l,
tmp = c[l].it)
f = f && (++tmp) == c[l].end;
return f;
}
......@@ -91,7 +88,6 @@ namespace Dune
};
public:
/// Constructor. Stores a pointer to the grid
explicit MultiMeshIterator (const GridImp* multiMesh, int level)
: multiMesh_(multiMesh)
......@@ -124,10 +120,12 @@ namespace Dune
}
}
/// Constructor which create the leaf-iterator
MultiMeshIterator (const GridImp* multiMesh)
: MultiMeshIterator{multiMesh, -1}
{}
/// Constructor which create the end leaf-iterator
MultiMeshIterator (const GridImp* multiMesh, bool endDummy)
: MultiMeshIterator{multiMesh, -1, endDummy}
{}
......@@ -160,11 +158,8 @@ namespace Dune
return macroIterators_ == that.macroIterators_;
}
protected:
// got to next entity in grid i
// iterator always points to a leaf entity
void increment (std::size_t i)
{
auto& entityStack = entityStacks_[i];
......@@ -191,9 +186,8 @@ namespace Dune
// 3. go down in tree until leaf entity
for (auto child = dereference(i); !entityTest(child);
child = dereference(i)) {
child = dereference(i))
entityStack.emplace(child);
}
}
/// Return true, if all stacks with size > stack[i].size are finished
......@@ -219,9 +213,8 @@ namespace Dune
// 2. go down in tree until leaf entity
for (auto child = dereference(i); !entityTest(child);
child = dereference(i)) {
child = dereference(i))
entityStack.emplace(child);
}
}
HostEntity dereference (std::size_t i) const
......@@ -232,13 +225,12 @@ namespace Dune
return *entityStacks_[i].top().it;
}
bool entityTest(HostEntity const& entity)
bool entityTest (HostEntity const& entity)
{
return entity.level() == level_ || entity.isLeaf() || !entity.isRegular();
}
private:
const GridImp* multiMesh_;
std::vector<HostGridLevelIterator> macroIterators_;
......
......@@ -13,8 +13,6 @@
#include <cstdlib>
#include <memory>
#include <dune/common/classname.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dune/common/parallel/mpihelper.hh>
......@@ -23,70 +21,28 @@
namespace Dune
{
/** \brief Construct structured cube and simplex grids in unstructured grid managers
*/
/// \brief Construct structured cube and simplex grids in unstructured grid managers
template <class GridType>
class StructuredGridFactoryWrapper
class StructuredGridBuilder
{
typedef typename GridType::ctype ctype;
using ctype = typename GridType::ctype;
static const int dim = GridType::dimension;
static const int dimworld = GridType::dimensionworld;
/** \brief Insert a structured set of vertices into the factory */
static void insertVertices(GridFactory<GridType>& factory,
const FieldVector<ctype,dimworld>& lowerLeft,
const FieldVector<ctype,dimworld>& upperRight,
const std::array<unsigned int,dim>& vertices)
{
FactoryUtilities::MultiIndex<dim> index(vertices);
// Compute the total number of vertices to be created
int numVertices = index.cycle();
// Create vertices
for (int i=0; i<numVertices; i++, ++index) {
// scale the multiindex to obtain a world position
FieldVector<double,dimworld> pos(0);
for (int j=0; j<dim; j++)
pos[j] = lowerLeft[j] + index[j] * (upperRight[j]-lowerLeft[j])/(vertices[j]-1);
for (int j=dim; j<dimworld; j++)
pos[j] = lowerLeft[j];
factory.insertVertex(pos);
}
}
// Compute the index offsets needed to move to the adjacent vertices
// in the different coordinate directions
static std::array<unsigned int, dim> computeUnitOffsets(const std::array<unsigned int,dim>& vertices)
{
std::array<unsigned int, dim> unitOffsets;
if (dim>0) // paranoia
unitOffsets[0] = 1;
for (int i=1; i<dim; i++)
unitOffsets[i] = unitOffsets[i-1] * vertices[i-1];
return unitOffsets;
}
public:
/** \brief Create a structured cube grid
If the grid dimension is less than the world dimension, the coefficients (dim+1,...,dimworld) in
the vertex coordinates are set to the corresponding values of the lowerLeft input argument.
\param lowerLeft Lower left corner of the grid
\param upperRight Upper right corner of the grid
\param elements Number of elements in each coordinate direction
*/
static std::unique_ptr<GridType> createCubeGrid(GridFactory<GridType>& factory,
/// \brief Create a structured cube grid
/**
* If the grid dimension is less than the world dimension, the coefficients
* (dim+1,...,dimworld) in the vertex coordinates are set to the corresponding
* values of the lowerLeft input argument.
*
* \param lowerLeft Lower left corner of the grid
* \param upperRight Upper right corner of the grid
* \param elements Number of elements in each coordinate direction
**/
static std::unique_ptr<GridType> createCubeGrid (
GridFactory<GridType>& factory,
const FieldVector<ctype,dimworld>& lowerLeft,
const FieldVector<ctype,dimworld>& upperRight,
const std::array<unsigned int,dim>& elements)
......@@ -94,7 +50,7 @@ namespace Dune
if (MPIHelper::getCollectiveCommunication().rank() == 0) {
// Insert uniformly spaced vertices
std::array<unsigned int,dim> vertices = elements;
for( size_t i = 0; i < vertices.size(); ++i )
for (std::size_t i = 0; i < vertices.size(); ++i)
vertices[i]++;
// Insert vertices for structured grid into the factory
......@@ -102,7 +58,7 @@ namespace Dune
// Compute the index offsets needed to move to the adjacent
// vertices in the different coordinate directions
std::array<unsigned int, dim> unitOffsets =
std::array<unsigned int,dim> unitOffsets =
computeUnitOffsets(vertices);
// Compute an element template (the cube at (0,...,0). All
......@@ -111,9 +67,9 @@ namespace Dune
std::vector<unsigned int> cornersTemplate(nCorners,0);
for (size_t i=0; i<nCorners; i++)
for (int j=0; j<dim; j++)
if ( i & (1<<j) )
for (std::size_t i = 0; i < nCorners; ++i)
for (int j = 0; j < dim; ++j)
if (i & (1<<j))
cornersTemplate[i] += unitOffsets[j];
// Insert elements
......@@ -122,29 +78,27 @@ namespace Dune
// Compute the total number of elementss to be created
int numElements = index.cycle();
for (int i=0; i<numElements; i++, ++index) {
for (int i = 0; i < numElements; ++i, ++index) {
// 'base' is the index of the lower left element corner
unsigned int base = 0;
for (int j=0; j<dim; j++)
for (int j = 0; j < dim; ++j)
base += index[j] * unitOffsets[j];
// insert new element
std::vector<unsigned int> corners = cornersTemplate;
for (size_t j=0; j<corners.size(); j++)
for (std::size_t j = 0; j < corners.size(); ++j)
corners[j] += base;
factory.insertElement(GeometryTypes::cube(dim), corners);
}
} // if(rank == 0)
}
// Create the grid and hand it to the calling method
return std::unique_ptr<GridType>(factory.createGrid());
}
static std::unique_ptr<GridType> createCubeGrid(const FieldVector<ctype,dimworld>& lowerLeft,
static std::unique_ptr<GridType> createCubeGrid (
const FieldVector<ctype,dimworld>& lowerLeft,
const FieldVector<ctype,dimworld>& upperRight,
const std::array<unsigned int,dim>& elements)
{
......@@ -152,20 +106,22 @@ namespace Dune
return createCubeGrid(factory, lowerLeft, upperRight, elements);
}
/** \brief Create a structured simplex grid
This works in all dimensions. The Coxeter-Freudenthal-Kuhn triangulation is
used, which splits each cube into dim! (i.e., dim faculty) simplices. See Allgower and Georg,
'Numerical Path Following' for a description.
If the grid dimension is less than the world dimension, the coefficients (dim+1,...,dimworld) in
the vertex coordinates are set to the corresponding values of the lowerLeft input argument.
\param lowerLeft Lower left corner of the grid
\param upperRight Upper right corner of the grid
\param elements Number of elements in each coordinate direction
*/
static std::unique_ptr<GridType> createSimplexGrid(GridFactory<GridType>& factory,
/// \brief Create a structured simplex grid
/**
* This works in all dimensions. The Coxeter-Freudenthal-Kuhn triangulation
* is used, which splits each cube into dim! (i.e., dim faculty) simplices.
* See Allgower and Georg, 'Numerical Path Following' for a description.
*
* If the grid dimension is less than the world dimension, the coefficients
* (dim+1,...,dimworld) in the vertex coordinates are set to the corresponding
* values of the lowerLeft input argument.
*
* \param lowerLeft Lower left corner of the grid
* \param upperRight Upper right corner of the grid
* \param elements Number of elements in each coordinate direction
**/
static std::unique_ptr<GridType> createSimplexGrid (
GridFactory<GridType>& factory,
const FieldVector<ctype,dimworld>& lowerLeft,
const FieldVector<ctype,dimworld>& upperRight,
const std::array<unsigned int,dim>& elements)
......@@ -173,40 +129,38 @@ namespace Dune
if (MPIHelper::getCollectiveCommunication().rank() == 0) {
// Insert uniformly spaced vertices
std::array<unsigned int,dim> vertices = elements;
for (std::size_t i=0; i<vertices.size(); i++)
for (std::size_t i = 0; i < vertices.size(); ++i)
vertices[i]++;
insertVertices(factory, lowerLeft, upperRight, vertices);
// Compute the index offsets needed to move to the adjacent
// vertices in the different coordinate directions
std::array<unsigned int, dim> unitOffsets =
std::array<unsigned int,dim> unitOffsets =
computeUnitOffsets(vertices);
// Loop over all "cubes", and split up each cube into dim!
// (factorial) simplices
FactoryUtilities::MultiIndex<dim> elementsIndex(elements);
size_t cycle = elementsIndex.cycle();
for (size_t i=0; i<cycle; ++elementsIndex, i++) {
std::size_t cycle = elementsIndex.cycle();
for (std::size_t i = 0; i < cycle; ++elementsIndex, ++i) {
// 'base' is the index of the lower left element corner
unsigned int base = 0;
for (int j=0; j<dim; j++)
for (int j = 0; j < dim; ++j)
base += elementsIndex[j] * unitOffsets[j];
// each permutation of the unit vectors gives a simplex.
std::vector<unsigned int> permutation(dim);
for (int j=0; j<dim; j++)
for (int j = 0; j < dim; ++j)
permutation[j] = j;
do {
// Make a simplex
std::vector<unsigned int> corners(dim+1);
corners[0] = base;
for (int j=0; j<dim; j++)
for (int j = 0; j < dim; ++j)
corners[j+1] =
corners[j] + unitOffsets[permutation[j]];
......@@ -214,23 +168,62 @@ namespace Dune
} while (std::next_permutation(permutation.begin(),
permutation.end()));
}
} // if(rank == 0)
}
// Create the grid and hand it to the calling method
return std::unique_ptr<GridType>(factory.createGrid());
}
// for default-constructable gridFactories
static std::unique_ptr<GridType> createSimplexGrid(const FieldVector<ctype,dimworld>& lowerLeft,
static std::unique_ptr<GridType> createSimplexGrid (
const FieldVector<ctype,dimworld>& lowerLeft,
const FieldVector<ctype,dimworld>& upperRight,
const std::array<unsigned int,dim>& elements)
{
GridFactory<GridType> factory;
return createCubeGrid(factory, lowerLeft, upperRight, elements);
}
private:
/// Insert a structured set of vertices into the factory
static void insertVertices (
GridFactory<GridType>& factory,
const FieldVector<ctype,dimworld>& lowerLeft,
const FieldVector<ctype,dimworld>& upperRight,
const std::array<unsigned int,dim>& vertices)
{
FactoryUtilities::MultiIndex<dim> index(vertices);
// Compute the total number of vertices to be created
int numVertices = index.cycle();
// Create vertices
for (int i = 0; i < numVertices; ++i, ++index) {
// scale the multiindex to obtain a world position
FieldVector<double,dimworld> pos(0);
for (int j = 0; j < dim; ++j)
pos[j] = lowerLeft[j] + index[j] * (upperRight[j]-lowerLeft[j])/(vertices[j]-1);
for (int j = dim; j < dimworld; ++j)
pos[j] = lowerLeft[j];
factory.insertVertex(pos);
}
}
/// \brief Compute the index offsets needed to move to the adjacent vertices
/// in the different coordinate directions
static std::array<unsigned int,dim> computeUnitOffsets (
const std::array<unsigned int,dim>& vertices)
{
std::array<unsigned int,dim> unitOffsets;
if (dim > 0) // paranoia
unitOffsets[0] = 1;
for (int i = 1; i < dim; ++i)
unitOffsets[i] = unitOffsets[i-1] * vertices[i-1];
return unitOffsets;
}
};
} // end namespace Dune
......
add_executable("testiterator" testiterator.cc)
target_link_dune_default_libraries("testiterator")
add_dune_alberta_flags(GRIDDIM 3 WORLDDIM 3 "testiterator")
add_executable("uggrid" uggrid.cc)
target_link_dune_default_libraries("uggrid")
DIM: 2
DIM_OF_WORLD: 2
number of elements: 4
number of vertices: 5
element vertices:
0 1 4
1 2 4
2 3 4
3 0 4
element boundaries:
0 0 1
0 0 2
0 0 3
0 0 4
vertex coordinates:
0.0 0.0
1.0 0.0
1.0 1.0
0.0 1.0
0.5 0.5
element neighbours:
1 3 -1
2 0 -1
3 1 -1
0 2 -1
DIM: 3
DIM_OF_WORLD: 3