libMesh
Function Elem::contains_point may return incorrect results
I was hitting an error which is caused by Elem::contains_point function. The mesh I was using contains a HEX 3D element, from which I build a side element in type of QUAD4. The four vertices of this QUAD4 is
-0.086602540378444323, -0.049999999999999753, 0;
-0.085949231939090279, -0.04999999999999899, 0;
-0.085785904829253024, -0.050471484754155152, 0;
-0.086112559048929102, -0.050282890852493239, 0;
which are showed in the picture:
As you can see, the element is actually a triangle, but with an extra vertex on one of its side.
In this case, calling Elem::contains_point with the last vertex returns a false although it should return true. @roystgnr suggested that
We do an inverse_map() as the first step of contains_point() and if that doesn't converge we return false. That (element)'s got to be singular.
I think the QUAD4 element is valid even though it is indeed a triangle shape. Fixing this may not be a trivial task though.
The “element is valid”? In what sense? The element map is not invertible at the degenerate point (and possibly along the entire degenerate edge) so any code that calls inverse_map there can potentially fail.
Yes to that but its shape functions are still defined and mass matrix is not singular (I suppose), so we still can find valid solutions. In this sense, I call it valid.
Fixing this may not be a trivial task though.
IMHO fixing this one very-specific case would be trivial: if we added the simple "make sure we're on the correct side of each line segment" override to Quad4::contains_point(), it would work for the barely-degenerate case here (though it would still fail for the concave cases epsilon away!) and it would be faster for other users too.
However, I don't think this would be useful for your final mesh (as soon as the element here got extruded into a Hex8 there'd no longer be a simple contains_point() override possible and you'd probably break again there), or for most users (once you proceed from here to an FE solve all bets are off; IIRC "the mapping Jacobian and its inverse are bounded" are postulates in convergence theorems for most formulations in H1 and up).
If I could figure out a way to support all degenerate elements, I'd think that was worth doing, even if the next step from there is "but you're on your own if you try to solve on these things", but if we support one type of degenerate element and still break on the other dozen categories that's just giving people false hope.
We could cut the element into simplex and only implement the method for simplex, which should not allow degenerated cases.
This only works for the planar Quad4 case, but yes that could be a possible optimization of Quad4::contains_point().