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

[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 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
>>>
>>>
>>> _______________________________________________
>>> Dune-devel mailing list
>>> Dune-devel at dune-project.org
>>> http://lists.dune-project.org/mailman/listinfo/dune-devel
>>>
>>>
>>
>>
>> _______________________________________________
>> Dune-devel mailing list
>> Dune-devel at dune-project.org
>> http://lists.dune-project.org/mailman/listinfo/dune-devel
>