Commit d580cccc authored by Praetorius, Simon's avatar Praetorius, Simon

added traversal over iregular elements

parent 31f111d0
......@@ -62,6 +62,32 @@ namespace Dune
return localGeometry(source, target).global(sourceLocal);
}
/// Return the geometry of the entity with maximal level
typename HostGridEntity::Geometry geometry () const
{
return max().geometry();
}
/// Return the type of the entity with maximal level
GeometryType type () const
{
return max().type();
}
/// The entities are always regular, since irregular entities are not allowed in multimesh
bool isRegular() const { return true; }
/// Return the maximal level
int level () const
{
return max().level();
}
/// Return whether the entity with minimal level has boundary intersections
bool hasBoundaryIntersections () const
{
return min().hasBoundaryIntersections();
}
/// \brief Return the entity seed which contains sufficient information
/// to generate the entity again and uses as little memory as possible.
......
......@@ -16,10 +16,10 @@ namespace Dune
protected:
/// Entity type of the hostgrid
typedef typename GridImp::HostGridType::Traits::template Codim<codim>::Entity HostEntity;
typedef typename GridImp::HostGrid::Traits::template Codim<codim>::Entity HostEntity;
/// EntitySeed type of the hostgrid
typedef typename GridImp::HostGridType::Traits::template Codim<codim>::EntitySeed HostEntitySeed;
typedef typename GridImp::HostGrid::Traits::template Codim<codim>::EntitySeed HostEntitySeed;
public:
enum { codimension = codim };
......@@ -61,6 +61,59 @@ namespace Dune
std::size_t gridIdx_;
};
template <class HG>
class MultiEntity;
// An entity seed for a vector of entities
template <int codim, class GridImp>
class MultiEntitySeed
{
protected:
using HostGrid = typename GridImp::HostGrid;
/// Entity type of the hostgrid
using HostEntity = typename HostGrid::Traits::template Codim<codim>::Entity;
/// EntitySeed type of the hostgrid
using HostEntitySeed = typename HostGrid::Traits::template Codim<codim>::EntitySeed;
public:
enum { codimension = codim };
/// Construct an empty (i.e. isValid() == false) seed.
MultiEntitySeed()
{}
/// \brief Create EntitySeed from hostgrid Entity
/**
* We call hostEntity.seed() directly in the constructor
* of MultiMeshEntitySeed to allow for return value optimization.
*/
MultiEntitySeed (const MultiEntity<HostGrid>& multiEntity)
{
// for (auto const& e : multiEntity)
// hostEntitySeeds_.emplace_back(e.seed());
}
/// Get stored HostEntitySeed
const std::vector<HostEntitySeed>& hostEntitySeeds () const
{
return hostEntitySeeds_;
}
/// Check whether it is safe to create an Entity from this Seed
bool isValid () const
{
return std::all_of(hostEntitySeeds_.begin(), hostEntitySeeds_.end(),
[](auto const& seed) { return seed.isValid(); });
}
private:
std::vector<HostEntitySeed> hostEntitySeeds_;
};
} // end namespace Dune
#endif // DUNE_MULTIMESH_ENTITYSEED_HH
......@@ -84,7 +84,7 @@ namespace Dune
const HostGrid& grid (std::size_t i) const
{
return multiMesh_->grid(i);
return (*multiMesh_)[i];
}
/// Obtain the level-indexSet
......@@ -255,7 +255,7 @@ namespace Dune
const HostGrid& grid (std::size_t i) const
{
return multiMesh_[i];
return (*multiMesh_)[i];
}
/// Obtain the level-indexSet
......
......@@ -55,8 +55,8 @@ namespace Dune
// go to first leaf entity on all grids
entityStacks_.reserve(multiMesh->size());
for (std::size_t i = 0; i < multiMesh->size(); ++i) {
maxLevel_[i] = (*multiMesh)[i].maxLevel();
entityStacks_.emplace_back((*multiMesh)[i].maxLevel());
maxLevel_[i] = multiMesh->maxLevel(i);
entityStacks_.emplace_back(multiMesh->maxLevel(i));
}
}
......@@ -152,13 +152,18 @@ namespace Dune
// 3. go down in tree until leaf entity
auto child = dereference(i);
for (; !this->levelReached(i, child); child = dereference(i)) {
assert(child.isRegular() && "No irregular elements allowed in multi-mesh traversal");
entityStack.emplace(child);
assert(int(entityStack.size()) <= maxLevel_[i]);
}
// 4. go up in tree again to the first regular entity, since
// irregular element can not be traversed in a multi-mesh sense
while (!child.isRegular() && !entityStack.empty()) {
entityStack.pop();
child = dereference(i);
}
return entityStack.size();
// assert(contains_[i](child) && "No valid child element found in gridView");
}
using Super::increment;
......@@ -172,7 +177,7 @@ namespace Dune
});
}
// got to first leaf entity on grid i
// go to first entity on grid i on the desired level
int initialIncrement (std::size_t i)
{
auto& entityStack = entityStacks_[i];
......@@ -186,16 +191,22 @@ namespace Dune
// 1. go down in tree until desired level is reached
auto child = dereference(i);
for (; !this->levelReached(i, child); child = dereference(i)) {
assert(child.isRegular() && "No irregular elements allowed in multi-mesh traversal");
entityStack.emplace(child);
assert(int(entityStack.size()) <= maxLevel_[i]);
}
// 2. go up in tree again to the first regular entity, since
// irregular element can not be traversed in a multi-mesh sense
while (!child.isRegular() && !entityStack.empty()) {
entityStack.pop();
child = dereference(i);
}
return entityStack.size();
// assert(contains_[i](child) && "No valid child element found in gridView");
}
using Super::initialIncrement;
// Return the current entity the hierarchic iterator or macro iterator pointing to
HostEntity dereference (std::size_t i) const
{
if (entityStacks_[i].empty()) {
......
......@@ -97,6 +97,8 @@ namespace Dune
using Super = GridDefaultImplementation<HG::dimension, HG::dimensionworld,
typename HG::ctype, MultiMeshFamily<HG> >;
using Self = MultiMesh;
public:
using HostGrid = HG;
......@@ -127,7 +129,7 @@ namespace Dune
MultiMesh () = default;
/// Return the number of grids handled by this MultiMesh
std::size_t size() const
std::size_t size () const
{
return grids_.size();
}
......@@ -157,6 +159,11 @@ namespace Dune
[](int level, auto const& grid) { return std::max(level, grid->maxLevel()); });
}
/// Return maximum level in the `idx`th grid
int maxLevel (std::size_t idx) const
{
return grids_[idx]->maxLevel();
}
/// Iterator to first entity of given codim on level
template <int codim, PartitionIteratorType PiType = All_Partition>
......@@ -303,12 +310,23 @@ namespace Dune
/// Create Entity from EntitySeed
template <class EntitySeed>
typename Traits::template Codim<EntitySeed::codimension>::Entity
entity(const EntitySeed& seed) const
entity (const EntitySeed& seed) const
{
std::size_t gridIdx = this->getRealImplementation(seed).gridIndex();
return grids_[gridIdx]->entity(this->getRealImplementation(seed).hostEntitySeed());
}
/// Create a MultiEntity from a MultiEntitySeed
template <int codim>
MultiEntity<HostGrid> entity (const MultiEntitySeed<codim, Self>& seed) const
{
assert(seed.isValid());
MultiEntity<HostGrid> multiEntity;
auto const& hostEntitySeeds = seed.hostEntitySeeds();
for (std::size_t i = 0; i < hostEntitySeeds.size(); ++i)
multiEntity.emplace_back(grids_[i]->entity(hostEntitySeeds[i]));
return multiEntity;
}
/** @name Grid Refinement Methods */
/** @{ */
......@@ -320,6 +338,24 @@ namespace Dune
grid->globalRefine(refCount);
}
/// Mark all entities in a multiEntitity for refinement or coarsening
bool mark (int refCount, const MultiEntity<HostGrid>& entities)
{
bool b = false;
for (std::size_t i = 0; i < grids_.size(); ++i)
b = grids_[i]->mark(refCount, entities[i]) || b;
return b;
}
/// Return the marks of the the entities in a multiEntity
std::vector<int> getMark (const MultiEntity<HostGrid>& entities) const
{
std::vector<int> marks(grids_.size());
for (std::size_t i = 0; i < grids_.size(); ++i)
marks[i] = grids_[i]->getMark(entities[i]);
return marks;
}
/// Returns true, if at least one entity is marked for adaption in any of the
/// handled grids
bool preAdapt()
......@@ -344,29 +380,39 @@ namespace Dune
/** @} */
/** @name Grid Parallelization Methods */
/** @{ */
/// Size of the overlap on the leaf level
unsigned int overlapSize (int codim) const
{
return 0u; // TODO: implement
return std::accumulate(grids_.begin(), grids_.end(), 0u,
[codim](unsigned int overlap, auto const& grid) {
return std::max(overlap, grid->overlapSize(codim)); });
}
/// Size of the ghost cell layer on the leaf level
unsigned int ghostSize (int codim) const
{
return 0u; // TODO: implement
return std::accumulate(grids_.begin(), grids_.end(), 0u,
[codim](unsigned int ghost, auto const& grid) {
return std::max(ghost, grid->ghostSize(codim)); });
}
/// Size of the overlap on a given level
unsigned int overlapSize (int level, int codim) const
{
return 0u; // TODO: implement
return std::accumulate(grids_.begin(), grids_.end(), 0u,
[level,codim](unsigned int overlap, auto const& grid) {
return std::max(overlap, grid->overlapSize(level,codim)); });
}
/// Size of the ghost cell layer on a given level
unsigned int ghostSize (int level, int codim) const
{
return 0u; // TODO: implement
return std::accumulate(grids_.begin(), grids_.end(), 0u,
[level,codim](unsigned int ghost, auto const& grid) {
return std::max(ghost, grid->ghostSize(level,codim)); });
}
/// dummy collective communication
......@@ -376,6 +422,22 @@ namespace Dune
return comm(0);
}
/// Rebalance all grids the same way based on an average load calculation.
bool loadBalance ()
{
assert(false && "Needs to be implemented.");
return false;
}
/// Rebalance all grids the same way based on an average load calculation.
template<class DataHandle>
bool loadBalance (DataHandle& data)
{
assert(false && "Needs to be implemented.");
return false;
}
/** @} */
public: // implementation details
......
......@@ -29,6 +29,7 @@
#include <dune/geometry/quadraturerules.hh>
#include <dune/grid/albertagrid.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/albertagrid/albertareader.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>
......@@ -40,6 +41,8 @@
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include "interpolation.hh"
using namespace Dune;
using namespace Dune::Indices;
......@@ -59,126 +62,17 @@ template <class Grids, class SignedDistFct>
void markElements(Grids& grids, SignedDistFct&& signedDistFct, double eps = 0.1)
{
for (auto const& elems : elements(grids.leafGridView())) {
auto geo = elems[0].geometry();
int level = elems[0].level();
auto geo = elems.geometry();
int level = elems.level();
int newLevel = 0;
for (int j = 0; j < geo.corners(); ++j)
newLevel = std::max(newLevel, std::abs(signedDistFct(geo.corner(j))) < eps ? 12 : 0);
newLevel = std::max(newLevel, std::abs(signedDistFct(geo.corner(j))) < eps ? 8 : 0);
int m = level < newLevel ? 1 : -1;
grids[0].mark(m, elems[0]);
grids[1].mark(m, elems[1]);
grids.mark(m, elems);
}
}
template <class Element, class GlobalCoordinate>
auto findLeafEntity(Element const& father, GlobalCoordinate const& global)
{
const int childLevel = father.level()+1;
// loop over all child Entities
const auto end = father.hend(childLevel);
for (auto it = father.hbegin(childLevel); it != end; ++it)
{
Element child = *it;
auto geo = child.geometry();
auto local = geo.local(global);
if (referenceElement(geo).checkInside(local))
{
// return if we found the leaf, else search through the child entites
if (child.isLeaf())
return std::make_pair(std::move(child), local);
else
return findLeafEntity(child, global);
}
}
assert(father.isLeaf());
return std::make_pair(father, father.geometry().local(global));
}
template <class Grids, class OldBasis, class NewBasis>
void mmInterpolate(Grids const& grids,
VectorType const& oldCoeff, OldBasis const& oldBasis,
VectorType& newCoeff, NewBasis const& newBasis)
{
using namespace Dune::Indices;
using LocalCoordinate = typename Grids::Traits::template Codim<0>::Geometry::LocalCoordinate;
auto oldLocalView = oldBasis.localView();
auto newLocalView = newBasis.localView();
newCoeff.resize(newBasis.dimension());
for (const auto& e : master_leaf_elements(grids, 1))
{
auto const& oldEntity = e[0];
auto const& newEntity = e[1]; // always a leaf entity
assert(newEntity.isLeaf());
newLocalView.bind(newEntity);
auto newGeo = newEntity.geometry();
auto const& node = newLocalView.tree();
auto const& fe = node.finiteElement();
std::size_t newFeSize = fe.size();
std::vector<double> newDOFs(newFeSize);
// entity not changed
if (newEntity.level() == oldEntity.level() && oldEntity.isLeaf()) {
oldLocalView.bind(oldEntity);
std::size_t oldFeSize = oldLocalView.tree().finiteElement().size();
assert(newFeSize == oldFeSize);
for (std::size_t i = 0; i < newFeSize; ++i)
newDOFs[i] = oldCoeff[oldLocalView.index(oldLocalView.tree().localIndex(i))];
}
// entity was coarsened
if (newEntity.level() == oldEntity.level() && !oldEntity.isLeaf()) {
auto evalLeaf = [&](LocalCoordinate const& x) {
thread_local std::vector<FieldVector<double,1>> shapeValues;
auto leafLocal = findLeafEntity(oldEntity, newGeo.global(x));
oldLocalView.bind(leafLocal.first);
auto const& node = oldLocalView.tree();
auto const& fe = node.finiteElement();
fe.localBasis().evaluateFunction(leafLocal.second, shapeValues);
double y = 0;
for (std::size_t i = 0; i < shapeValues.size(); ++i)
y += shapeValues[i] * oldCoeff[oldLocalView.index(node.localIndex(i))];
return y;
};
auto const& node = newLocalView.tree();
auto const& fe = node.finiteElement();
using FFC = Functions::FunctionFromCallable<double(LocalCoordinate), decltype(evalLeaf)>;
fe.localInterpolation().interpolate(FFC(evalLeaf), newDOFs);
}
// entity was refined
else if (newEntity.level() >= oldEntity.level()) {
oldLocalView.bind(oldEntity);
auto oldGeo = oldEntity.geometry();
auto evalOld = [&](LocalCoordinate const& x)
{
thread_local std::vector<FieldVector<double,1>> shapeValues;
auto const& fe = oldLocalView.tree().finiteElement();
fe.localBasis().evaluateFunction(oldGeo.local(newGeo.global(x)), shapeValues);
double y = 0;
for (std::size_t i = 0; i < shapeValues.size(); ++i)
y += shapeValues[i] * oldCoeff[oldLocalView.index(node.localIndex(i))];
return y;
};
using FFC = Functions::FunctionFromCallable<double(LocalCoordinate), decltype(evalOld)>;
fe.localInterpolation().interpolate(FFC(evalOld), newDOFs);
}
for (std::size_t i = 0; i < newFeSize; ++i)
newCoeff[newLocalView.index(node.localIndex(i))] = newDOFs[i];
}
}
int main(int argc, char** argv)
{
......@@ -188,7 +82,8 @@ int main(int argc, char** argv)
std::string filename = argv[1];
// using HostGrid = Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>;
using HostGrid = Dune::AlbertaGrid<2,2>;
// using HostGrid = Dune::AlbertaGrid<2,2>;
using HostGrid = Dune::UGGrid<2>;
using WorldVector = typename HostGrid::template Codim<0>::Geometry::GlobalCoordinate;
using Grid = MultiMesh<HostGrid>;
AlbertaReader<Grid> albertaReader;
......@@ -197,7 +92,7 @@ int main(int argc, char** argv)
std::unique_ptr<Grid> gridPtr(factory.createGrid());
auto& grids = *gridPtr;
grids.globalRefine(6);
grids.globalRefine(3);
double eps = 0.1;
WorldVector shift(0.0);
......@@ -231,6 +126,8 @@ int main(int argc, char** argv)
double timeInterpolate = 0;
double timeWriteFile = 0;
Interpolation<Grid, VectorType> interpolation(grids);
// assemble a sequence of systems for finer and finer grids
for (int i = 1; i < 20; ++i) {
std::cout << i << "\n";
......@@ -254,7 +151,7 @@ int main(int argc, char** argv)
timeBasisUpdate += t.elapsed();
t.reset();
mmInterpolate(grids, oldPhase, basis0, phase, basis1); // interpolate from old grid to new grid
interpolation(oldPhase, basis0, phase, basis1); // interpolate from old grid to new grid
timeInterpolate += t.elapsed();
t.reset();
......
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