Sampling at a point in space

[majosm]

@dshtey2 is looking for a way to evaluate DOFArray-based fields at given spatial positions. This issue is intended for discussion of that functionality: what code might already exist for this, what needs to be added, etc.

cc @inducer

[dshtey2]

Here is some additional information on my idea for how to accomplish this. I am looking for a routine that can succinctly accomplish this, or one that can perform an even simpler procedure:

1. For a given sample point or points, determine which element(s) contain(s) that or those points (through some convex hull checking function, such as taking the sums of distances to all nodal points in an element and minimizing over all elements)
2. After finding the desired element, reconstruct the solution in that element from its basis functions, and evaluate the solution at the desired sample point(s).

[inducer]

Good question! Various parts of this might exist in a way that we might be able to factor out. First of all, @xywei might already have this functionality in a form that can be reused. @xywei, could you comment?

• Finding elements is probably the hardest bit. It looks like https://github.com/xywei/volumential/blob/main/volumential/interpolation.py solves this via geometric queries from https://documen.tician.de/boxtree/. That's fast if you need many points. https://github.com/inducer/pytools/blob/main/pytools/spatial_btree.py could be used for smaller point counts.
• Next, you need to find the reference coordinates of each real-space point with respect to the element mapping. This routine already does that and, AFAICT, only needs a slight generalization.
• Next, use the to-modal connection to compute a modal representation of the function you're looking to evaluate.
• Last, you can symbolic expressions for the basis functions from modepy, as in this test. You could easily toss them into a loopy kernel to produce the function evaluations.

All of the above is for single-rank. Multi-rank is a bit more complicated, but many building blocks can be reused in that case.

[majosm]

@inducer Do you have a sense of how well-behaved the basis functions are when evaluated outside of the unit element coordinate range?

I had suggested to @dshtey2 first using the elements' bounding boxes to find those that are "close" to the query point, then computing the unit coordinates of the point relative to them, and then using those unit coordinates to determine whether the query point is inside each element or not. This is sort of working, but we're running into some intermittent convergence issues when the point is outside of the element. For example, for a uniform 2D mesh on (0,1)x(0,1) and a query point (0, 0), the bounding box search finds 2 candidate elements (the lower and upper triangle halves of the first "box"). The coordinate computation succeeds for the lower triangle, but not for the upper. Interestingly, if you move the query point slightly off-axis in the y direction (either way), they both converge. Any idea if this is expected?

If necessary we could do an alternative test to find out if the point lies inside the element, but I only know how to do that for non-curved elements -- so I think we'd lose some generality there.

[inducer]

While I don't think that in general there's much of an expectation that evaluating functions outside of the reference element would yield meaningful results, in your case (given purely affine mappings), you should be able to evaluate basically anywhere. (For context, in #228, I observed what seems like a similar non-uniqueness issue, but that was for a genuinely curvi mesh.)

When you say "the coordinate computation does not succeed", what does that mean? Gauss-Newton doesn't converge? Something else happens?

[majosm]

While I don't think that in general there's much of an expectation that evaluating functions outside of the reference element would yield meaningful results, in your case (given purely affine mappings), you should be able to evaluate basically anywhere. (For context, in #228, I observed what seems like a similar non-uniqueness issue, but that was for a genuinely curvi mesh.)

Is there maybe a different way we should be approaching this?

When you say "the coordinate computation does not succeed", what does that mean? Gauss-Newton doesn't converge?

Yep. (Even after cranking up the iteration count.)

For a query point (0, 0) I get (each column is a candidate element; the second column is the one that fails):

src_unit_query_points=array([
       [-1., 21.],
       [-1.,  1.]])

For (0, 1e-5) I get:

src_unit_query_points=array([
       [-1.      ,  1.      ],
       [-0.999994,  0.999994]])

and for (0, -1e-5) I get:

src_unit_query_points=array([
       [-1.      ,  1.      ],
       [-1.000006,  1.000006]])

[inducer]

Huh, interesting. That should be fixable. Could you put code to reproduce this somewhere? (Maybe a temporary PR)

[majosm]

Huh, interesting. That should be fixable. Could you put code to reproduce this somewhere? (Maybe a temporary PR)

The code is in mirgecom at the moment, in @dshtey2's branch. I tweaked his test into a pytest test and put it up at illinois-ceesd/mirgecom#479. (It can also be run manually with pytest test_sampling.py.)