[Dune] [Dune-devel] Checking whether a given point is contained in a given element

[Carsten Gräser]

Am 08.11.2012 16:11, schrieb Oliver Sander:
> Hi Martin,
>> well, the restriction is only necessary for the local method on
>> pyramid geometries (e.g., simplices, four-sided pyramid). Assuming
>> that you will always use CachedMultiLinearGeometry, it will never be
>> called for simplices, so you can safely remove the assertion. In the
>> case of a 3d pyramid, however, the trouble will persist.
>>
>> Please note, that the old generic geometry implementation behaved as
>> described above with the assertion commented out. So this is not a new
>> bug, but only a new assertion.
> That is strange, because my code used to work with the old implementation.
> But it is good to know that in the eyes of the implementor the two codes
> should implement the same algorithm.  That way maybe I'll find the time to
> put them both side-by-side and see where they differ.
> 
>>
>> While the mapping might be globally well-defined (I did not check
>> this), it becomes singular in the tip for a four-sided pyramid. This
>> causes trouble with the Newton scheme used in local, because the
>> derivative can no longer be computed. Maybe we need to resort to some
>> bisection scheme in this case, but that might become extremely slow.
> That may be a way out.  If we get a good automatic switching between Newton
> and bisection then the run-time penalty may not even be all that big.
> 
> Slightly off-topic: actually I see more problems with the Newton solver
> lurking.
> Currently the algorithm terminates when the norm of the correction drops
> below a fixed number.  That will fail when considering very large elements,
> because then there may be Newton corrections that are above the threshold,
> but too small to be added to the current iterate.  Hence the algorithm will
> iterate forever.
I can't see a problem here. By the affine covariance of the Newton
method it is completely invariant under linear transformations
of the image space like, e.g., scaling the element. I'm not
sure if this is also true for the implemented algorithm which
is in fact Gauß-Newton and not Newton for dim!=dimworld.

However, the iterates are still local coordinates and should
thus be O(1) in general. IMHO that's a place where controling the
absolute error does really make sense.

Best,
Carsten

>> I am not sure, whether the problem can be resolved satisfactorily.
>> Maybe we should revive the method checkInside on the geometry
>> returning true if and only if a point in global coordinates lies
>> within the image of the reference mapping (up to a tolerance). If the
>> method returns true, it should be safe to call local. My feeling is
>> that this method is simpler to implement that a local method that has
>> to work for arbitrary points.
> That would mean another interface method for a fairly exotic use case.
> I'd rather improve the 'local' method.
> 
>>
>> A last question: What is the exact semantics of local, if the point
>> does not lie inside the geometric realization of an entity? For
>> surface grids, you currently have the following (unintuitive)
>> behavior: Assume that a point lies slightly "above" the actual
>> geometric realization. Now, local will return the local coordinates of
>> its orthogonal projection and checkInside will return true despite the
>> fact that the point actually lies outside the geometric realization.
> Yes, that question was implied in my initial mail: what are the exact
> semantics of 'local'?
> Personally, I think the best is what's implemented and documented in
> multilineargeometry.
> It's easy to understand and easy to implement.
> 
> best,
> Oliver
>>
>> Best,
>>
>> Martin
>>
>>
>>
>> On 11/08/2012 02:29 PM, Oliver Sander wrote:
>>> Am 08.11.2012 13:20, schrieb Dedner, Andreas:
>>>> I think that the local method should not abort if the point is
>>>> inside the
>>>> reference element
>>>> or not - and one of the reasons is that you need as Oliver said you
>>>> need the
>>>> local method
>>>> to verify if a point is inside. Also I can think of applications
>>>> where I
>>>> would like to
>>>> represent for example the barycenter of neighbors in the local
>>>> coordinate
>>>> system
>>>> of an element - requiring me to use its local method...
>>>> I would thus suggest to remove that check (and think its
>>>> mathematically fine
>>>> because these
>>>> transformations are global mappings in all the cases I can think of...)
>>> I certainly wouldn't object to that solution.  Strangely enough, though,
>>> after removing all assertions from MultiLinearGeometry, the code still
>>> crashed due to infs and nans appearing (only for simplices, now that I
>>> think about it).  Does the MLG code rely on the points being inside the
>>> reference element?
>>> best,
>>> Oliver
>>>
>>>>
>>>> Andreas
>>>> ________________________________________
>>>> From: dune-devel-bounces+a.s.dedner=warwick.ac.uk at dune-project.org
>>>> [dune-devel-bounces+a.s.dedner=warwick.ac.uk at dune-project.org] on
>>>> behalf of
>>>> Oliver Sander [sander at igpm.rwth-aachen.de]
>>>> Sent: 08 November 2012 12:09
>>>> To: dune-devel at dune-project.org
>>>> Subject: [Dune-devel] Checking whether a given point is contained in a
>>>> given    element
>>>>
>>>> Dear Dune,
>>>> I have recently run into the following issue: I have a given point
>>>> in world
>>>> space and an element.  I want to decide whether the element contains
>>>> the point.
>>>> So far, I have solved this by computing the coordinates of the point
>>>> wrt to the element (using Geometry::local(point)), and then calling
>>>> checkInside on the appropriate reference element.  That used to work
>>>> nicely.
>>>> Unfortunately, some grids have recently changed to MultiLinearGeometry
>>>> for their geometry implementations.  MLG has fairly tight
>>>> preconditions:
>>>> in particular, its 'local' method aborts with an assertion failure
>>>> if the
>>>> input argument is not contained in the element.
>>>> That does not seem like an unreasonable restriction to make.  After
>>>> all,
>>>> strictly speaking, local element coordinates only exist on the element.
>>>> Unfortunately, now my old trick doesn't work anymore.  And I can't
>>>> think of another elegant why to decide whether a given point is
>>>> contained
>>>> in a given element.
>>>>
>>>> Any ideas?
>>>>
>>>> cheers,
>>>> Oliver