MultilinearGeometry: method local fails when global point is not inside the element.
Every now and then one needs to check whether a global point is inside an element or not. The way to do this in DUNE is the following:
auto local = geometry.local( globalPoint );
auto refElem = referenceElement( geometry );
if( refElem.checkInside( local ) )
{
// do something
}This fails with MultilinearGeometry because in this case the Newton algorithm inside the local method won't converge and there is not max iteration criterion implemented.
Activity
added bug discussion labels
mentioned in commit cda92064
mentioned in merge request !125 (merged)
mentioned in commit b71d15ba
Is this issue possibly related to dune-grid#84 (closed) where @felix.mueller pointed out that the Hierarchic search in grids may fail due to not reached tolerance in
checkInside()?@robert.kloefkorn Do you have any proposal on how to fix this in the interface?
@simon.praetorius: Most likely the same problem I would say. I stumbled upon this in a similar application.
-
@carsten.graeser: Actually I don't have a good solution besides introducing a max iteration count that breaks the while loop. That worked in my case but introduces a free parameter, i.e. whats a good maxIter number.
@robert.kloefkorn I'm not sure if this really solves the problem. Furthermore I think that the proposed use case will still fail in general. For an affine geometry with world-dim>dim
local()will probably compute the local coordinate of the normal projection which may be inside although the global coordinate is far away.Before trying to fix this, it should be clear what the precise semantics of
local()is.-
@carsten.graeser: But we do agree that calling local on a 2d simplex with an point outside of the element should at least terminate, right? I didn't say that my fix was the most brilliant possible, other and better ideas very welcome here.
For the 2,3 case you pointed out you should add a test case.
I agree that maybe the semantics of checkInside + local should be discussed, that's why I labled this issue "discussion".
@robert.kloefkorn We can even agree that it should always terminate. I also don't object against your patch, because it fixes one problem. But it will not make your example work in general.
However, writing a test for the 2,3 case is not possible unless it is clear what the precise semantics is. I'm also not convinced that this should work at all. But that's part of the discussion.
I don't think it was ever properly documented anywhere, I remember that the interpretation once was...
I project to local coordinates. Now it might be, that the point is projected onto the element, although it was significantly away in normal direction. We can them account for this by transforming back into global coordinates and checking the distance to the initial point.
-
Ok, lets do it like this: We fix the easy cases, such as dim=dimworld=2 and dim=dimworld=3. We mark this as might not work in dim=2,dimworld=3 case and we put the general work flow up for discussion for the next meeting.
It is still unclear to me why the algorithm does not terminate for the 2d simplex case. My patch is just a quick-fix.
There are two issues here (1) a non terminating Newton loop (2) semantics of local methods for points outside the physical element
Concerning 1: an infinite loop is always bad in my opinion so I don't understand @carsten.graeser 's comment that stopping this from happening is not something we can agree on. Of course what the correct solution is, is a different issue but the minimum I would think is an exception
Concerning 2: the only implementation of the hierarchical search seems to be to use the local method on the 'geometryInFather' - at least I can't really think of a better way. So allowing 'local(x)' to be called at least on local geometries with x outside the element seems a necessity - or does someone know of a different way of doing this?
How about changing 'local' to return a std::optional? Or to introduce a second 'safeLocal' method returning an optional. We could then add an exception in the standard 'local' method and catch that in the 'safe' version. The check for the dimWorld>dimDomain case could also then be added. I would like to directly change 'local' but that would change things for people who used the return value also for points 'outside' e.g. to figure out which child to look at in a hierarchic travels. Of course that would be using undefined behaviour I guess...
changed milestone to %DUNE 2.7.0
mentioned in commit 38586f59
changed milestone to %DUNE 2.7.0
