Point evaluation raises Not In Domain for point in domain

[taupalosaurus]

This issue contains an example where the Function.at() function raises a PointNotInDomainError error for a point that is obviously in the domain. This bug seems to depend very much on the machine/platform, so I hope it can be reproduced on machines other than mine. I give two versions here, one that works on my mac (os x 10.9) and one that works on a 64-bit linux server (Ubuntu 4.8.2-19ubuntu1).

The domain is a 2D unit square, meshed with triangles. It is stored as a gmsh mesh (the file is attached at the end of the issue). The minimal code below replicates the problem, comment uncomment the appropriate lines if you want the mac/linux64 version

from firedrake import *
mesh = Mesh("mesh_mac.msh")
#mesh = Mesh("mesh_linux.msh")
V = FunctionSpace(mesh, "CG", 1)
u = Function(V)
u.interpolate(Expression("x[0]"))
crd = [0.8003202852202465, 0.5373200327560730]
crd = [ 0.7670007215677767, 0.5308980767672131] # linux
try :
    u.at(crd)
except PointNotInDomainError :
    print "####  Vertex not in domain: %1.3f %1.3f" % (crd[0], crd[1])

Attached (at the end) is another source file, which does the same as above, but also computes the barycentric coordinates of the point for each cell of the vertex and prints them if they are all three positive (i.e. the point is inside this cell). We notice that one coordinate is null, which means the point is on an edge of the triangle.
Copying here Lawrence's email:

So what you're saying is that the point is right on the boundary between two cells? I guess there are two ways the point location can fail. Either the point is not found in any bounding box. Or, the point is found by libspatialindex in a bounding box. But that bucket does not actually contain the cell in question, so then the linear search for point location may fail due to floating point rounding.

It seems likely that we are in the second case, and the fact that it is so sensible to the machine/architecture speaks in favor of this possibility.

bug.zip

[miklos1]

Thanks for the report, I can reproduce the problem.

[miklos1]

libspatialindex finds that the point is within the bounding box of four cells. Physical to reference space transformation yields:

(0.333333, 0.333333) -> (0.797195, 0.531502) <=> (0.80032, 0.53732)
FINAL: (0.962841, -0.569417)
(0.333333, 0.333333) -> (0.803104, 0.53286) <=> (0.80032, 0.53732)
FINAL: (0.49871, 0.50129)
(0.333333, 0.333333) -> (0.79878, 0.538044) <=> (0.80032, 0.53732)
FINAL: (0.50129, -1.28231e-14)
(0.333333, 0.333333) -> (0.803378, 0.541148) <=> (0.80032, 0.53732)
FINAL: (-0.635631, 1.24134)
####  Vertex not in domain: 0.800 0.537

As we can see, the point is almost at the midpoint of an edge in two triangles, slightly outside in both cases.

[miklos1]

Increasing the limit to the number of Newton iterations seems to help in this case:

(0.333333, 0.333333) -> (0.797195, 0.531502) <=> (0.80032, 0.53732)
(0.962841, -0.569417) -> (0.80032, 0.53732) <=> (0.80032, 0.53732)
FINAL: (0.962841, -0.569417)
(0.333333, 0.333333) -> (0.803104, 0.53286) <=> (0.80032, 0.53732)
(0.49871, 0.50129) -> (0.80032, 0.53732) <=> (0.80032, 0.53732)
FINAL: (0.49871, 0.50129)

Increasing the tolerance from 1e-14 to 1e-13 seems to also help:

(0.333333, 0.333333) -> (0.797195, 0.531502) <=> (0.80032, 0.53732)
FINAL: (0.962841, -0.569417)
(0.333333, 0.333333) -> (0.803104, 0.53286) <=> (0.80032, 0.53732)
FINAL: (0.49871, 0.50129)

[miklos1]

The short-term solution could be to tune parameters. There are three tunable parameters in firedrake/pointquery_utils.py, at lines 179, 259 and 260. I let @taupalosaurus propose values that work for him, since he's the only one who've complained so far.

[taupalosaurus]

Thanks Miklós for working this out.
I'm not sure which values would be good, I'm going to run some tests to see what would suit me. But changing thresholds does not look like a good long term solution, right ? You can always find a point that will be too close to an edge for a certain threshold... Which raises a question: why do these points fail, and not points located on mesh vertices ?

[miklos1]

There are several things one can do.

Tune parameters

This won't make the thing any more robust though, and all parameter setting is just a different compromise.

Immersed manifolds

Currently, we have a function that returns a logical value whether the point is in the cell or not. Instead, that function could return non-negative real, describing the shortest distance between the point and a point of the cell. This is exactly zero if the the point is clearly inside the cell. Then we could choose the cell which the point is closest to, given that the distance is below a certain threshold.

Pro:

• Support for immersed manifolds
• More robust cell selection in case of ambiguity.

Con:

• For a given cell, it doesn't make the in/out question any more robust, and still needs a threshold.

Floating-point aware algorithms to decide if point is in cell

The nice thing about the current approach is that it is very general. One just need the tabulation of the coordinate element and the matrix inversion snippets to carry out the physical to reference coordinates transformation, then whether the inside check is simply done using the reference coordinates. We will need the reference coordinates any way to evaluate the field there.

However, as we have seen, this isn't very robust with floating-point arithmetics. A point near a facet may appear to be outside both cells incident to that facet - a logical impossibility for a continuous coordinate space. For a non-immersed, affine simplex, any facet defines hyperplane. One can check for each facet if the vertex that isn't incident to the facet and the evaluation point are in the same half-space. If this is the case for every facet, then the point is inside the cell. This way the "which side of a facet a point is on" question will be evaluated by exactly the same numerical operations, irrespective of which cell the facet was approached from. So if a point is near a facet, only the following outcomes are possible, even considering floating-point arithmetics:

• Both cells (incident to the facet) claim that the point is exactly on the facet.
• One cell claims the point is inside, the other claims it is outside.
• One cell claims the point is outside, the other claims it is inside.

In particular, the following outcomes are not possible:

• Both cells claim they don't contain the cell.
• Both cells claim they contain the cell.

DMPlex also has a similar algorithm for quadrilateral cells. Can we do a similar thing for triangular prism and hexes? How can we generalise this for higher-order coordinate spaces?

Pro:

• Robust in/out decisions. Threshold might only be useful to allow a tolerance around the domain boundary.

Con:

• Requires cell-specific algorithms.
• Difficult to generalise the approach to all meshes supported by Firedrake.
• Likely more expensive, since we still need to the physical to reference space transformation to do the evaluation.

More advanced arithmetics

If floating-point errors are a problem, one could be exact using rational numbers. Unfortunately, immersed manifolds need a square root to calculate the pseudo-inverse, an operation that isn't closed on rational numbers. So one can try arbitrary-precision floating-point arithmetics, but then the question arises: what precision to use? Interval arithmetics can tell if the precision isn't high enough.

Pro:

• "Unbreaks" the general algorithm we have already implemented.

Con:

• These advanced arithmetics are slow as their implementation isn't fully in hardware.
• These advanced arithmetics require additional libraries.
• Their integration to the existing code is not trivial at all.

And finally the answers to your questions:

But changing thresholds does not look like a good long term solution, right ?

Correct.

why do these points fail, and not points located on mesh vertices ?

I guess it is because floating-point errors make it random whether a point is found to be in or out, and vertices have more incident cells than facets, consequently it is less likely that a point near a vertex will be outside every cell that that vertex is incident to.

[taupalosaurus]

Would you want to use your better approach all the time, or only when the first approach fails ?
I don't know well how Firedrake works, so I may be missing something, but since the bug is not so frequent, the latter seems reasonable, and in that case, the better approach is the robust one evaluating on which side of each facet the point is. Why would it not be generic to all meshes ? Can't we always loop over the edges/facets of a cell ?

Your method "Immersed manifolds" is interesting when you have a point that is indeed slightly out of the domain (which happens all the time in mesh adaptation when you generate a new mesh)

Also, there should soon be (or maybe it's even already there) some localization mechanism in dmplex for unstructured meshes, which we might want to use, and may bring some answers.

[miklos1]

Why would it not be generic to all meshes?

Not all facets are planar. For example:

• The vertices of a quadrilateral may not be coplanar.
• For higher-order cells, the facets are higher-order as well, so they are generally not planar.
• For immersed manifolds, facets geometrically have co-dimension greater than 1.

[miklos1]

I mean, the idea to iterate over all facets and check on which side of the facet a point is should work in almost all cases. The question is rather how can we generalise for all supported cell types the algorithm to determine on which side of a facet a point is.

[taupalosaurus]

Okay, so how do you deal with non-planar cell facets for now ?

Also, what do people use point evaluation for ? (I guess knowing how people use a function help coming up with a good strategy).
My current usage, is to interpolate a function from a mesh adapted at time t^n to a mesh adapted at time t^n+1 (for each vertex of the latter mesh, I evaluate the function on the former mesh). It seems that pragmatic is likely to create vertices or or very close to edges that trigger the bug (the bigger the meshes, the more such vertices). I will not always use this method, so other usages should be considered...
For now, raising the epsilon from firedrake/pointquery_utils.py, line 179 works, but what is the harm of doing that ?

[wence-]

FWIW, one could imagine using the current (non-robust) scheme, but adding a "nearest cell" evaluation mode. I imagine, although do not know, that if your numerical scheme is sensitive to interpolation errors due to floating point round-off, it will probably also be highly sensitive to mesh-adaptivity (or numerous other issues) as well, so "do we care?".

A "nearest cell" evaluation mode would also be useful in the general case: imagine that I want to approximately represent a field on some coarse mesh on a finer one, where the boundary facets do not coincide. It is likely that I will have fine vertices (say) that are outside the domain of the coarse mesh.

[taupalosaurus]

I agree. If we use consistent interpolation, the precision is not great anyway and we probably don't care (it might need to be confirmed, notably if very small elements are involved ?). And we do need a mechanism for points slightly out of the domain: this happens all the time when you adapt a boundary and try to recover the curvature - which we are working on in Pragmatic.

[rbpiccinini]

Hello, I'm also having this "not in domain" error for a function evaluation of a point inside an extruded mesh.
I'll try to make a reproducible example to attach here.

[vabhi6]

I am getting a Point Not In Domain error for even a simple example like this. The code runs fine without errors if I remove the call to at(). I would appreciate it if I could get some help on how I can debug this.

from firedrake import *
import numpy as np

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "CG", 1)
v = TestFunction(V)
u = TrialFunction(V)
a = inner(grad(u), grad(v)) * dx
L = Constant(0) * v * dx
bc1 = DirichletBC(V, 0.0, 1)
bc2 = DirichletBC(V, 1.0, 2)
sol = Function(V)
problem1 = LinearVariationalProblem(a, L, sol, bcs=[bc1, bc2], constant_jacobian=True)
solver1 = LinearVariationalSolver(problem1, solver_parameters={"ksp_type": "cg", "pc_type": "jacobi"})
solver1.solve()
crd=[0.525,0.525]
print(sol.at(crd))

[wence-]

If you perturb the point by a small amount do you still get an error?

[vabhi6]

If you perturb the point by a small amount do you still get an error?

I tried crd=[0.5256,0.5051] and still got the error:

domain <Mesh #1> does not contain point [0.5256 0.5051]