Re: [Libmesh-devel] FE::inverse_map() not converging...

[Roy S. <roy...@ic...>]

On Wed, 21 Jun 2006, Kirk, Benjamin (JSC-EG) wrote:

> Which type of elements are these?

They're QUAD4 elements.

> For Quads the zero point corresponds to the centroid of the element,
> which should be as good a guess as any.  For triangles, the reference
> element is an isosceles triangle with the right angle at the origin,
> hence the zero point is at a node.  Probably not the best place,
> especially if that angle on the physical element approaches (or exceeds)
> 180 degrees.

Starting at element centroids (or approximate element centroids) would
probably be more robust, but we've looked at things in detail and
found two more potential problems:

1.  Elem::contains_point passes a tolerance of 1e-4 to inverse_map()
but then uses the default tolerance of 1e-6 for on_reference_element()

This means that it's slightly more likely for contains_point to return
a false positive (not a big problem), but it also means it's possible
for contains_point to return a false negative (which can break
programs).

2.  Because Elem::contains_point uses inverse_map() and because it
gets called even for points which are far away from the physical
elements (in the point locator functions is where Derek's code tripped
it), if the physical elements are badly distorted then the Newton
iteration won't converge.

For now Derek's just going to silence the Newton warning messages and
return a grossly invalid point when Newton fails to converge.  Is that
safe to do in general, though?  It seems like there are many
situations where inverse_map gets called and we expect to get a valid
answer, and it would be nice to keep asserting convergence in DEBUG
mode for such calls.  We should avoid changing the declaration of
inverse_map() if possible, so there's no way for inverse_map() itself
to distinguish between the "must converge" and "may not converge"
cases, but on the other hand inverse_map() get called dosens of times
by the library and it would be tedious to add asserts() after every
such call.
---
Roy