Hierarchic search failing for points very close but outside an element/ checkInside tolerance insufficient
Summary
The precision of geometry.local(global)
depends on the element size. HierarchicSearch
uses referenceElement.checkInside(local)
which uses a fixed tolerance. On higher refinement levels of the grid this check fails and HierarchicSearch
triggers an exception.
This is not strictly a dune-grid issue but also relates to dune-geometry.
Description
The current implementation of referenceElement.checkInside()
[dune/geometry/referenceelementimplementation.hh] uses a fixed tolerance to check if a local coordinate is included in the reference element.
HierarchicSearch::hFindEntity()
[dune/grid/utility/hierarchicsearch.hh] recursively transforms a given global point global
into local coordinates using local = element.geometry().local(global)
and then proceeds to check if the result is within an element using checkInside(local)
.
If a starting point of the recursion is close to an element (at level 0) by a distance eps
but outside, it can still be considered inside, but may not be inside any child element on a finer grid level. This triggers an "Unexpected internal Error: none of the children of the entity ... contains coordinate ...".
In the case of affine geometries this is due to the observation that the magnitude |M|
of the coefficients of the transformation matrix M
used in the call to global/local()
is dependent on the size of the element in global coordinates. A global point x+eps
with inaccuracy eps
is transformed into Mx+epsNew
with epsNew ~ |M|*eps
. As a result the fixed tolerance of checkInside()
may be sufficient when transforming a global point into a coarse element, but insufficient when transforming into a finer one.
I ran into the issue with when interpolating data during grid adaption. I tried to find a child containing each of the Lagrange nodes of a father element and at a certain refinement level the process failed due to not finding any (even though the nodes are all inside the father element and so one of its children must contain them). The code used for finding the children was very similar to that of HierarchicSearch
.
Possible Solution
A possible solution to the problem consists of the following two parts:
- Add a function to the reference element
bool checkInside( const Coordinate &local, const ctype tolerance )
{
return Impl::template checkInside< ctype, dim >( type().id(), dim, local, tolerance );
}
- Replace calls to
checkInside(local)
withcheckInside(local, tol)
and compute a tolerance within each step of the hierarchic search, e.g.
...
constexpr double baseTol = 1e-14;
...
const auto& vol = geo.volume();
const auto tol = baseTol / std::pow(vol, 1.0/dimw);
if (referenceElement( geo ).checkInside(local, tol))
...
Reproducing the bug
You can find small tests showing the problem in https://gitlab.mn.tu-dresden.de/s5272799/bug-checkinside.
src/bug-checkinside_short
prints the required tolerance to find an entity containing a test point outside the coarse element with a given distance as well as the tolerance computed by the method above. It also performs a classic hierarchic search and the proposed alternative.
Output on my machine with default parameters:
===== UGGrid =====
Level required tolerance adaptive tolerance
0 9.992007e-15 1.414214e-14
1 1.998401e-14 2.828427e-14
2 3.996803e-14 5.656854e-14
3 7.993606e-14 1.131371e-13
4 1.598721e-13 2.262742e-13
5 3.197442e-13 4.525483e-13
6 6.394885e-13 9.050967e-13
7 1.278977e-12 1.810193e-12
8 2.557954e-12 3.620387e-12
Actual tolerance used: 1.421085e-14
Trying to find leaf entity containing the test point ... - Failed:
Exception [hFindEntity:/usr/local/include/dune/grid/utility/hierarchicsearch.hh:109]: {Dune::HierarchicSearch<Dune::UGGrid<2>, Dune::IndexSet<Dune::UGGrid<2> const, Dune::UGGridLeafIndexSet<Dune::UGGrid<2> const>, unsigned int, std::vector<Dune::GeometryType, std::allocator<Dune::GeometryType> > > > const} Unexpected internal Error: none of the children of the entity {level=0 partition=interior center=(0.666667 0.333333) corners=[(0 0) (1 0) (1 1)]} contains coordinate (1 0.125). Children are: [{level=1 partition=interior center=(0.666667 0.333333) corners=[(0.5 0) (1 0.5) (0.5 0.5)]} {level=1 partition=interior center=(0.833333 0.666667) corners=[(1 0.5) (1 1) (0.5 0.5)]} {level=1 partition=interior center=(0.833333 0.166667) corners=[(0.5 0) (1 0) (1 0.5)]} {level=1 partition=interior center=(0.333333 0.166667) corners=[(0 0) (0.5 0) (0.5 0.5)]}].
Trying to find leaf entity containing the test point with variable tolerance search ... - Success!
...
You can see the hierarchic search failing on reaching level 1, while the alternative implementation succeeds.
src/real-situation
computes the tolerance requirement to find the nodes of a dune-functions Lagrange basis of order 3 of a leaf grid entity inside its father entity.
Output:
Required tolerance for nodal basis points of child to be found within father element:
0 0: 0.000000e+00
3.333333e-01 0.000000e+00: 9.464651e-15
6.666667e-01 0.000000e+00: 4.718448e-15
1.000000e+00 0.000000e+00: -0.000000e+00
0.000000e+00 3.333333e-01: 0.000000e+00
3.333333e-01 3.333333e-01: 0.000000e+00
6.666667e-01 3.333333e-01: 0.000000e+00
0.000000e+00 6.666667e-01: 0.000000e+00
3.333333e-01 6.666667e-01: 0.000000e+00
0.000000e+00 1.000000e+00: 0.000000e+00
The 2nd and 3rd node are transformed to points outside the reference element. You can see here, that transformation to points outside occurs in real examples.