From: Roy Stogner <roystgnr@ic...>  20060807 20:18:39

On Mon, 7 Aug 2006, Tim Kr=F6ger wrote: > On Fri, 4 Aug 2006, Roy Stogner wrote: > >> The inverse map iteration definitely won't converge for some elements >> if the physical point is outside the element, but that shouldn't cause >> a problem for Elem::constains_point() anymore; if its call to >> inverse_map() doesn't converge, the result is a distant point that >> will make sure on_reference_element() is false. > > No: If the nonlinear solver doesn't converge, it will call error() and = thus=20 > crash. At least, this call to error() should be removed  perhaps by=20 > giving the solver class an optional argument. That's exactly why the inverse_map function has the "bool secure" argument. If secure is false (as it is when contains_point() calls inverse_map() ), then Newton isn't expected to converge and inverse_map just exits after 10 iterations. Of course, inverse_map() is *supposed* to exit by returning an outofbounds point, not just whatever its last iterate was. Derek, didn't we fix this? Or did you send me a patch I forgot to commit to CVS? > The other problem that I see when leaving the linear search in is that = when=20 > I will later implement the possibiliy to evaluate the MeshFunction outs= ide=20 > the domain covered by the grid, the linear search will slow this down=20 > considerabely. That is undesirable. > What about adding some heuristic safety margin around the bounding box = in=20 > the case of quadratic mappings? Well, not really satisfying ... As long as the margin is a proveable upper bound, I'll be happy with it. Trying to figure out what would be an upper bound on a quadratic Hex makes my head spin, though. > Perhaps better: Enable the linear search only if quadratic mappings are= =20 > involved. That sounds good too, but deciding whether or not quadratic mappings are involved is easier said than done. I'd go so far as to say that secondorder geometric elements with linear mappings is one of the most common usage patterns in libMesh, since we need the extra nodes in the secondorder elements to store degree of freedom numbers for any p>1 finite element. > Any suggestions? No good ones. All I can think to do is warn about the problem. When you add an interface for outofmesh MeshFunction extensions, let that interface turn off (well, replace) the linear search, and add an assert() or #ifdef DEBUG block in that code path which verifies that it doesn't see any quadratic mappings.  Roy 