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

[Martin Nolte]

Shouldn't move the discussion to FlySpray?

On 11/08/2012 11:25 PM, Carsten Gräser wrote:
> 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
>
> _______________________________________________
> Dune mailing list
> Dune at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune

-- 
Dr. Martin Nolte <nolte at mathematik.uni-freiburg.de>

Universität Freiburg                                   phone: +49-761-203-5630
Abteilung für angewandte Mathematik                    fax:   +49-761-203-5632
Hermann-Herder-Straße 10
79104 Freiburg, Germany