libMesh / libmesh

libMesh github repository
http://libmesh.github.io
GNU Lesser General Public License v2.1
654 stars 287 forks source link

contains_point() fails on edges #2938

Open BalticPinguin opened 3 years ago

BalticPinguin commented 3 years ago

Recently, I get from time to time assertions from find_point_neighbor() that points are not contained in elements which are found by the PointLocatorTree. This happens when I try to map a point on an elements edge (which seems quite likely for me just for numerical reasons, given that we use Newton-iteration).

My idea is to write some function similar to FEAbstract::on_reference_element but which moves points from slightly outside to the boarder of elements and add this into the Newton-iteration in FEMap::inverse_map, hoping not to destroy some slow convergence with points that belong just a bit outside of that element?

Since this is a bit of work, I'd like to discuss this in advance; other possible solutions are very welcome. Best regards

jwpeterson commented 3 years ago

If I understand correctly, you are saying that PointLocatorTree returns an Elem when you search for a Point, but then when you later try to inverse_map() that Point, it fails? There are a couple of different tolerances on PointLocator (set_contains_point_tol() and set_close_to_point_tol()) so one thing you could check is that they are not somehow too loose, and are returning Elems which the Point really is not in.

Another comment is that the inverse_map() function itself does some additional "sanity checking" that you can possibly turn off (see the secure and extra_checks parameters below)

static Point inverse_map (const unsigned int dim,
                            const Elem * elem,
                            const Point & p,
                            const Real tolerance = TOLERANCE,
                            const bool secure = true,
                            const bool extra_checks = true);

in your case, that would then allow you to accept some points just barely outside of elements. I'm not sure how one would do this reliably (or if one would want to do this) in general, though.