Commit 695e7aba authored by Praetorius, Simon's avatar Praetorius, Simon

Multientity

parent bf547e21
#install headers
install(FILES
mmentity.hh
mmentityseed.hh
mmgeometry.hh
mmgridfactory.hh
mmgridview.hh
mmiterator.hh
......
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_MULTIMESH_ENTITY_HH
#define DUNE_MULTIMESH_ENTITY_HH
#include <vector>
#include <dune/grid/common/geometry.hh>
#include "mmentityseed.hh"
#include "mmgeometry.hh"
namespace Dune
{
template <class GridImp>
class MultiEntity
: public std::vector<typename GridImp::HostGridType::Traits::template Codim<0>::Entity>
{
private:
using ctype = typename GridImp::ctype;
public:
enum { dimension = GridImp::dimension };
/// The type of a local geometry
using LocalGeometry = Dune::Geometry<dimension, dimension, const GridImp, MultiMeshLocalGeometry>;
using EntitySeed = typename GridImp::Traits::template Codim<0>::EntitySeed;
/// Entity in the host grid
using HostGridEntity = typename GridImp::HostGridType::Traits::template Codim<0>::Entity;
/// Containertype
using Super = std::vector<HostGridEntity>;
/// Constructor from std::vector
using Super::vector;
/// Return a local geometry of source in target
static LocalGeometry localGeometry (HostGridEntity const& source, HostGridEntity const& target)
{
using LocalGeometryImp = MultiMeshLocalGeometry<dimension, dimension, const GridImp>;
return LocalGeometry{LocalGeometryImp{source, target}};
}
/// Return a local geometry of source_i'th entity in target_t'th entity
LocalGeometry localGeometry (std::size_t source_i, std::size_t target_i) const
{
using LocalGeometryImp = MultiMeshLocalGeometry<dimension, dimension, const GridImp>;
return LocalGeometry{LocalGeometryImp{(*this)[source_i], (*this)[target_i]}};
}
/// \brief Return the entity seed which contains sufficient information
/// to generate the entity again and uses as little memory as possible.
/**
* The MultiMeshEntitySeed contains the HostGridEntitySeed and the index of
* the grid it is extracted from.
**/
EntitySeed seed (std::size_t entity_i) const
{
using EntitySeedImp = MultiMeshEntitySeed<0, const GridImp>;
return EntitySeedImp{(*this)[entity_i], entity_i};
}
/// Return the entity with maximal level
HostGridEntity const& max() const
{
// TODO: cache max-element
auto it_max = std::max_element(this->begin(), this->end(),
[](auto const& e1, auto const& e2) { return e1.level() < e2.level(); });
return *it_max;
}
/// Return the entity with minimal level
HostGridEntity const& min() const
{
// TODO: cache min-element
auto it_min = std::min_element(this->begin(), this->end(),
[](auto const& e1, auto const& e2) { return e1.level() < e2.level(); });
return *it_min;
}
};
} // end namespace Dune
#endif // DUNE_MULTIMESH_ENTITY_HH
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_MULTIMESH_LOCALGEOMETRY_HH
#define DUNE_MULTIMESH_LOCALGEOMETRY_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/typetraits.hh>
#include <dune/grid/common/geometry.hh>
namespace Dune
{
template <int mydim, int coorddim, class GridImp>
class MultiMeshLocalGeometry
: public GeometryDefaultImplementation<mydim, coorddim, GridImp, MultiMeshLocalGeometry>
{
private:
using ctype = typename GridImp::ctype;
public:
// The codimension of this entitypointer wrt the host grid
enum { codimension = GridImp::HostGridType::dimension - mydim };
enum { dimensionworld = GridImp::HostGridType::dimensionworld };
// select appropriate hostgrid geometry via typeswitch
using HostLocalGeometry = typename GridImp::HostGridType::Traits::template Codim<codimension>::LocalGeometry;
/// entity in the host grid
using HostGridEntity = typename GridImp::HostGridType::Traits::template Codim<codimension>::Entity;
//! type of jacobian transposed
using JacobianInverseTransposed = typename HostLocalGeometry::JacobianInverseTransposed;
using JacobianTransposed = typename HostLocalGeometry::JacobianTransposed;
/// Constructor from host entities
MultiMeshLocalGeometry (const HostGridEntity& sourceEntity, const HostGridEntity& targetEntity)
{
hostLocalGeometries_.reserve(sourceEntity.level() - targetEntity.level());
auto entity = sourceEntity;
for (int l = sourceEntity.level(); l > targetEntity.level(); --l) {
assert( entity.hasFather() );
hostLocalGeometries_.push_back(entity.geometryInFather());
if (l > targetEntity.level()+1)
entity = entity.father();
}
}
/// Return the element type identifier
GeometryType type () const
{
return hostLocalGeometries_[0].type();
}
/// Return wether we have an affine mapping
bool affine () const
{
return std::all_of(hostLocalGeometries_.begin(), hostLocalGeometries_.end(),
[](auto const& localGeo) { return localGeo.affine(); });
}
/// Returns true if the point is in the current element
bool checkInside (const FieldVector<ctype, mydim>& local) const
{
return hostLocalGeometries_[0].checkInside(local);
}
/// Return the number of corners of this element. Corners are numbered 0...n-1
int corners () const
{
return hostLocalGeometries_[0].corners();
}
/// Access to coordinates of corners. Index is the number of the corner
const FieldVector<ctype, coorddim> corner (int i) const
{
FieldVector<ctype, coorddim> c = hostLocalGeometries_[0].corner(i);
for (std::size_t j = 1; j < hostLocalGeometries_.size(); ++j)
c = hostLocalGeometries_[j].global(c);
return c;
}
/// Maps a local coordinate within reference element to global coordinate in element
FieldVector<ctype, coorddim> global (const FieldVector<ctype, mydim>& local) const
{
auto c = hostLocalGeometries_[0].global(local);
for (std::size_t j = 1; j < hostLocalGeometries_.size(); ++j)
c = hostLocalGeometries_[j].global(c);
return c;
}
/// Maps a global coordinate within the element to a local coordinate in its reference element
FieldVector<ctype, mydim> local (const FieldVector<ctype, coorddim>& global) const
{
auto c = hostLocalGeometries_.back().local(global);
std::size_t n = hostLocalGeometries_.size()-1;
for (std::size_t j = 1; j < hostLocalGeometries_.size(); ++j)
c = hostLocalGeometries_[n-j].local(c);
return c;
}
/// Return the factor appearing in the integral transformation formula
ctype integrationElement (const FieldVector<ctype, mydim>& local) const
{
ctype dx = 1;
auto c = local;
for (std::size_t j = 0; j < hostLocalGeometries_.size(); ++j) {
dx *= hostLocalGeometries_[j].integrationElement(c);
if (j < hostLocalGeometries_.size()-1)
c = hostLocalGeometries_[j].global(local);
}
return dx;
}
/// The Jacobian matrix of the mapping from the reference element to this element
JacobianInverseTransposed jacobianInverseTransposed (const FieldVector<ctype, mydim>& local) const
{
auto J = hostLocalGeometries_[0].jacobianInverseTransposed(local);
auto c = hostLocalGeometries_[0].global(local);
for (std::size_t j = 1; j < hostLocalGeometries_.size(); ++j) {
J.leftmultiply(hostLocalGeometries_[j].jacobianInverseTransposed(c));
if (j < hostLocalGeometries_.size()-1)
c = hostLocalGeometries_[j].global(local);
}
return J;
}
/// Return the transposed of the Jacobian
JacobianTransposed jacobianTransposed (const FieldVector<ctype, mydim>& local) const
{
auto J = hostLocalGeometries_[0].jacobianTransposed(local);
auto c = hostLocalGeometries_[0].global(local);
for (std::size_t j = 1; j < hostLocalGeometries_.size(); ++j) {
J.rightmultiply(hostLocalGeometries_[j].jacobianTransposed(c));
if (j < hostLocalGeometries_.size()-1)
c = hostLocalGeometries_[j].global(local);
}
return J;
}
private:
std::vector<HostLocalGeometry> hostLocalGeometries_;
};
} // end namespace Dune
#endif // DUNE_MULTIMESH_LOCALGEOMETRY_HH
......@@ -10,9 +10,7 @@
#include <dune/grid/common/exceptions.hh>
#include <dune/grid/common/gridenums.hh>
/** \file
* \brief The MultiMeshIterator class
*/
#include "mmentity.hh"
namespace Dune
{
......@@ -40,8 +38,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> >
MultiEntity<GridImp>,
MultiEntity<GridImp> const&>
{
private:
// LevelIterator to the equivalent entity in the host grid
......@@ -148,7 +146,7 @@ namespace Dune
}
/// dereferencing
std::vector<HostEntity> dereference () const
MultiEntity<GridImp> const& dereference () const
{
// update entries in multiEntity that have changed
for (std::size_t i = 0; i < entityStacks_.size(); ++i) {
......@@ -253,7 +251,7 @@ namespace Dune
std::vector<bool> incrementAllowed_;
int level_ = -1;
mutable std::vector<HostEntity> multiEntity_;
mutable MultiEntity<GridImp> multiEntity_;
};
} // end namespace Dune
......
......@@ -13,6 +13,7 @@
#include <dune/grid/common/gridfactory.hh>
// The components of the MultiMesh interface
#include "mmentity.hh"
#include "mmentityseed.hh"
#include "mmgridview.hh"
#include "mmiterator.hh"
......@@ -82,22 +83,25 @@ namespace Dune
: public GridDefaultImplementation<HostGrid::dimension, HostGrid::dimensionworld,
typename HostGrid::ctype, MultiMeshFamily<HostGrid::dimension, HostGrid> >
{
template <int codim, PartitionIteratorType pitype, class GridImp_>
template <int codim, PartitionIteratorType pitype, class GridImp>
friend class MultiMeshIterator;
template <class GridImp_>
template <class GridImp>
friend class MultiMeshLevelGridView;
template <class GridImp_>
template <class GridImp>
friend class MultiMeshLeafGridView;
template <class GridImp>
friend class MultiEntity;
friend class GridFactory<MultiMesh<HostGrid> >;
using HostGridType = HostGrid;
using Super = GridDefaultImplementation<HostGrid::dimension, HostGrid::dimensionworld,
typename HostGrid::ctype, MultiMeshFamily<HostGrid::dimension, HostGrid> >;
public:
using HostGridType = HostGrid;
/// Type of the used GridFamily for this grid
using GridFamily = MultiMeshFamily<HostGrid::dimension,HostGrid>;
......@@ -222,7 +226,7 @@ namespace Dune
/// Access to the GlobalIdSet
const typename HostGrid::GlobalIdSet& globalIdSet (std::size_t idx) const
{
return grids_[idx].globalIdSet();
return grids_[idx]->globalIdSet();
}
/// Return the globalIdSet for grid 0
......@@ -236,7 +240,7 @@ namespace Dune
/// Access to the LocalIdSet
const typename HostGrid::LocalIdSet& localIdSet (std::size_t idx) const
{
return grids_[idx].localIdSet();
return grids_[idx]->localIdSet();
}
/// Return the localIdSet for grid 0
......@@ -253,7 +257,7 @@ namespace Dune
if (level < 0 || level > grids_[idx].maxLevel()) {
DUNE_THROW(GridError, "levelIndexSet of nonexisting level " << level << " requested!");
}
return grids_[idx].levelIndexSet(level);
return grids_[idx]->levelIndexSet(level);
}
/// Return the levelIndexSet for grid 0
......@@ -267,7 +271,7 @@ namespace Dune
/// Access to the LeafIndexSet
const typename HostGrid::LeafIndexSet& leafIndexSet(std::size_t idx) const
{
return grids_[idx].leafIndexSet();
return grids_[idx]->leafIndexSet();
}
/// Return the leafIndexSet for grid 0
......@@ -280,7 +284,7 @@ namespace Dune
/// View for a grid level for All_Partition of the idx'th grid
typename HostGrid::LevelGridView levelGridView(std::size_t idx, int level) const
{
return grids_[idx].levelGridView(level);
return grids_[idx]->levelGridView(level);
}
/// Return \ref MultiMeshLevelGridView, i.e. a view on the combined
......@@ -290,7 +294,7 @@ namespace Dune
/// View for the leaf grid for All_Partition of the idx'th grid
typename HostGrid::LeafGridView leafGridView(std::size_t idx) const
{
return grids_[idx].leafGridView();
return grids_[idx]->leafGridView();
}
/// Return \ref MultiMeshLeafGridView, i.e. a view on the combined
......
......@@ -9,10 +9,16 @@
#include <dune/common/exceptions.hh> // We use exceptions
#include <dune/common/timer.hh>
#if HAVE_ALBERTA
#include <dune/grid/albertagrid.hh>
#endif
#include <dune/grid/yaspgrid.hh>
#if HAVE_UG
#include <dune/grid/uggrid.hh>
#endif
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#endif
#include <dune/multimesh/multimesh.hh>
#include <dune/multimesh/mmgridfactory.hh>
......@@ -55,6 +61,7 @@ int main(int argc, char** argv)
{
MPIHelper::instance(argc, argv);
#if HAVE_DUNE_ALUGRID
FieldVector<double,2> lower_left = {-1.5, -1.5};
FieldVector<double,2> bbox = {1.5, 1.5};
std::array<unsigned int,2> num_elements = {2, 2};
......@@ -88,4 +95,5 @@ int main(int argc, char** argv)
printGrid2(grid[0]);
printGrid(grid);
#endif
}
......@@ -9,10 +9,16 @@
#include <dune/common/exceptions.hh> // We use exceptions
#include <dune/common/timer.hh>
#if HAVE_ALBERTA
#include <dune/grid/albertagrid.hh>
#endif
#include <dune/grid/yaspgrid.hh>
#if HAVE_UG
#include <dune/grid/uggrid.hh>
#endif
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#endif
#include <dune/multimesh/multimesh.hh>
#include <dune/multimesh/mmgridfactory.hh>
......@@ -28,10 +34,10 @@ void printGrid(Grid const& grid)
volatile std::size_t n = 0;
Dune::Timer t;
for (auto const& entities : elements(grid.leafGridView())) {
n += grid[0].leafIndexSet().index(entities[0]);
n += grid.leafIndexSet(0).index(entities[0]);
std::cout << "indices = [";
for (std::size_t i = 0; i < entities.size(); ++i) {
std::cout << (i > 0 ? ", " : "") << grid[i].leafIndexSet().index(entities[i]);
std::cout << (i > 0 ? ", " : "") << grid.leafIndexSet(i).index(entities[i]);
}
std::cout << "]\n";
}
......@@ -55,15 +61,7 @@ int main(int argc, char** argv)
{
MPIHelper::instance(argc, argv);
#if GRID == 1
FieldVector<double,3> lower_left = {-1.5, -1.5, -1.5};
FieldVector<double,3> bbox = {1.5, 1.5, 1.5};
std::array<int,3> num_elements = {2, 2, 2};
using HostGrid = YaspGrid<3, EquidistantOffsetCoordinates<double,3>>;
MultiMesh<HostGrid> grid(3, lower_left, bbox, num_elements);
#elif GRID == 2
#if GRID == 2 && HAVE_DUNE_ALUGRID
FieldVector<double,2> lower_left = {-1.5, -1.5};
FieldVector<double,2> bbox = {1.5, 1.5};
......@@ -76,7 +74,7 @@ int main(int argc, char** argv)
auto gridPtr = Factory::createSimplexGrid(gridFactory, lower_left, bbox, num_elements);
auto& grid = *gridPtr;
#elif GRID == 3
#elif GRID == 3 && HAVE_ALBERTA
assert(argc > 1);
std::string filename = argv[1];
......@@ -85,6 +83,14 @@ int main(int argc, char** argv)
grid[0].globalRefine(8);
grid[1].globalRefine(8);
#else
FieldVector<double,3> lower_left = {-1.5, -1.5, -1.5};
FieldVector<double,3> bbox = {1.5, 1.5, 1.5};
std::array<int,3> num_elements = {2, 2, 2};
using HostGrid = YaspGrid<3, EquidistantOffsetCoordinates<double,3>>;
MultiMesh<HostGrid> grid(3, lower_left, bbox, num_elements);
#endif
std::cout << "number of grids = " << grid.size() << "\n";
......
......@@ -9,7 +9,9 @@
#include <dune/common/exceptions.hh> // We use exceptions
#include <dune/common/timer.hh>
#if HAVE_UG
#include <dune/grid/uggrid.hh>
#endif
#include <dune/multimesh/multimesh.hh>
#include <dune/multimesh/mmgridfactory.hh>
#include <dune/multimesh/utility/structuredgridbuilder.hh>
......@@ -20,6 +22,7 @@ int main(int argc, char** argv)
{
MPIHelper::instance(argc, argv);
#if HAVE_UG
FieldVector<double,2> lower_left = {-1.5, -1.5};
FieldVector<double,2> bbox = {1.5, 1.5};
std::array<unsigned int,2> num_elements = {2, 2};
......@@ -28,4 +31,5 @@ int main(int argc, char** argv)
GridFactory<MultiMesh<HostGrid> > factory(1);
auto gridPtr = Factory::createSimplexGrid(factory, lower_left, bbox, num_elements);
#endif
}
dune_add_test(SOURCES testvolume.cc)
dune_add_test(SOURCES testnumentities.cc)
\ No newline at end of file
dune_add_test(SOURCES testnumentities.cc)
dune_add_test(SOURCES testtransform.cc)
\ No newline at end of file
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifdef HAVE_CONFIG_H
# include "config.h"
#endif
#include <functional>
#include <iostream>
#include <numeric>
#include <dune/common/filledarray.hh>
#include <dune/common/test/testsuite.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/multimesh/multimesh.hh>
using namespace Dune;
template <std::size_t dim, class Test>
void test_dim(Test& test)
{
FieldVector<double,dim> lower; lower = -1.5;
FieldVector<double,dim> upper; upper = 1.5;
auto num_elements = filledArray<dim>(2);
using HostGrid = YaspGrid<dim, EquidistantOffsetCoordinates<double,dim>>;
MultiMesh<HostGrid> grid(3, lower, upper, num_elements);
grid[0].globalRefine(1);
grid[1].globalRefine(2);
grid[2].globalRefine(3);
for (auto const& entities : elements(grid.leafGridView())) {
auto trafo21 = entities.localGeometry(entities[2], entities[1]);
auto trafo20 = entities.localGeometry(entities[2], entities[0]);
auto trafo10 = entities.localGeometry(entities[1], entities[0]);
auto refElem = referenceElement(entities[2].geometry());
auto c = refElem.position(0,0);
// check whether transform.global returns a point inside the entity
test.check( refElem.checkInside(trafo21.global(c)) );
test.check( refElem.checkInside(trafo20.global(c)) );
test.check( refElem.checkInside(trafo10.global(c)) );
// check whether local is inverse of global
test.check( trafo21.local(trafo21.global(c)) == c );
test.check( trafo20.local(trafo20.global(c)) == c );
test.check( trafo10.local(trafo10.global(c)) == c );
}
}
int main(int argc, char** argv)
{
MPIHelper::instance(argc, argv);
Dune::TestSuite test;
test_dim<1>(test);
test_dim<2>(test);
test_dim<3>(test);
return test.exit();
}
......@@ -41,10 +41,7 @@ bool test_dim(Test& test)
// calculate volume by summing up the entity volumes of the smalles leaf entities
double volume_leaf = 0.0;
for (auto const& entities : elements(grid.leafGridView())) {
auto it_small = std::max_element(entities.begin(), entities.end(),
[](auto const& e1, auto const& e2) { return e1.level() < e2.level(); });
auto geo_small = it_small->geometry();
auto geo_small = entities.max().geometry();
volume_leaf += geo_small.volume();
}
std::cout << "volume_leaf(elements<" << dim << ">) = " << volume_leaf << "\n";
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment