Commit 482fa372 authored by Müller, Felix's avatar Müller, Felix
Browse files

Change implementation to use geometryInFather

parent 98999a15
......@@ -50,7 +50,8 @@ preAdapt(B const& basis, C const& coeff, bool mightCoarsen)
{
GridView gv = basis.gridView();
LocalView lv = basis.localView();
auto const& idSet = gv.grid().localIdSet();
auto const& grid = gv.grid();
auto const& idSet = grid.localIdSet();
Traversal::forEachLeafNode(lv.tree(), [&](auto const& node, auto const& tp) {
nodeDataTransfer_[tp].preAdaptInit(lv, coeff, node);
......@@ -73,9 +74,12 @@ preAdapt(B const& basis, C const& coeff, bool mightCoarsen)
return;
// Interpolate from possibly vanishing elements
auto maxLevel = gv.grid().maxLevel();
auto maxLevel = grid.maxLevel();
using std::sqrt;
typename Grid::ctype const checkInsideTolerance = sqrt(std::numeric_limits<typename Grid::ctype>::epsilon());
using Seed = typename Grid::template Codim<0>::EntitySeed;
std::vector<Seed> seeds(maxLevel+1);
for (auto const& e : entitySet(basis))
{
auto father = e;
......@@ -91,7 +95,10 @@ preAdapt(B const& basis, C const& coeff, bool mightCoarsen)
bool init = true; // init flag for first call on new father element
bool restrictLocalCompleted = false;
auto hItEnd = father.hend(maxLevel);
for (auto hIt = father.hbegin(maxLevel); hIt != hItEnd; ++hIt) {
for (auto hIt = father.hbegin(maxLevel); hIt != hItEnd; ++hIt)
{
seeds[hIt->level()] = hIt->seed(); // Save element in hierarchy to access geometryInFather at each step later
if (!hIt->isLeaf())
continue;
......@@ -111,12 +118,21 @@ preAdapt(B const& basis, C const& coeff, bool mightCoarsen)
// Transfers input father-local point x into child-local point y
// Returns false if x is not inside the child
auto xInChild = [&](LocalCoordinate const& x) -> BoolCoordPair {
LocalCoordinate local = childGeo.local(geo.global(x));
// TODO(FM): Using an implementation detail as workaround for insufficient
// tolerance, see https://gitlab.dune-project.org/core/dune-grid/issues/84
bool isInside = Dune::Geo::Impl::checkInside(refTypeId, Geometry::mydimension, local, checkInsideTolerance);
return BoolCoordPair(isInside, std::move(local));
LocalCoordinate y = x;
for (std::size_t i = father.level()+1; i <= child.level(); ++i)
{
auto currentElement = grid.entity(seeds[i]);
auto const& geoInFather = currentElement.geometryInFather();
y = geoInFather.local(y);
// TODO(FM): Using an implementation detail as workaround for insufficient
// tolerance, see https://gitlab.dune-project.org/core/dune-grid/issues/84
bool isInside = Dune::Geo::Impl::checkInside(Dune::referenceElement(geoInFather).type().id(), Geometry::mydimension, y, checkInsideTolerance);
if (!isInside)
return BoolCoordPair(false, std::move(y));
}
return BoolCoordPair(true, std::move(y));
};
// Cache result of xInChild for subsequent calls from other basis nodes
// TODO(FM): Disable for single-node basis
ChildCache childCache;
auto xInChildCached = [&](LocalCoordinate const& x) -> BoolCoordPair {
......@@ -201,10 +217,16 @@ void InterpolationDataTransfer<B,C>::adapt(B const& basis, C& coeff)
lv.bind(child);
// coordinate transform from child to father element
auto xInFather = [&fatherGeo, childGeo = child.geometry()]
(LocalCoordinate const& x) -> LocalCoordinate
auto xInFather = [&](LocalCoordinate const& x) -> LocalCoordinate
{
return fatherGeo.local(childGeo.global(x));
auto y = x;
auto currentElement = child;
while (currentElement.level() != father.level())
{
y = currentElement.geometryInFather().global(y);
currentElement = currentElement.father();
}
return y;
};
Traversal::forEachLeafNode(lv.tree(), [&](auto const& /*node*/, auto const& tp) {
......
Supports Markdown
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