AtPoints Operators - Add "foreign" field to support mapping between distinct element sets with consistent point sets

[zatkins-dev]

Motivation

The particular use case is that I need to project a quantity from different mesh elements to the same set of points. The points here are "shared" between the two sides - they have the same quadrature weights and areas, etc., but they are grouped differently and have different reference coordinates for each set of elements.

After interpolating from each set of FEM nodes to points, I need to be able to evaluate a QFunction which takes the pair of values at each point as an input to compute some physics.

Finally, the resulting values are interp-transposed back to one of the element sets and finally added into the FEM L-vector.

Current Approach

Currently, our at-point operators require that all element restrictions refer to the same elements, so that we can group points by their enclosing elements. As a result, this operation requires 2 "setup" operators (one for each set of mesh elements) that just interp to points, then 2 "output" operators (again, one for each set of mesh elements) to compute the physics and interp-transpose back to the mesh nodes.

This approach is particularly problematic for assembly (let's focus on diagonal), since our assembly routines only affect active outputs with the current active inputs. As such, all of the coupling between the two sets of mesh elements is lost. See https://gitlab.com/micromorph/ratel/-/merge_requests/1213 for an example of this approach.

Proposed Solution

Introduce CeedOperatorAtPointsSetForeignField() which, in addition to the arguments to CeedOperatorSetField(), takes the reference coordinates and restriction for the layout of the other elements.

In implementations, we would have to first interpolate from the foreign fields to an intermediate point vector, then use an internally created at-points element restriction to extract the values needed for the QFunction.

Additionally, for diagonal assembly, we would have to include the foreign field's elements in the loop over input elements.

Here's a hand drawn diagram of the process as I foresee it.

rough_foreign_field_diagram.pdf

[zatkins-dev]

@jeremylt what do you think about this allowing different per-field point coordinate and source elements idea? It would certainly be less efficient than our current approach, but IMO it's the cleanest way to do what I'm trying to do

[jeremylt]

I'm pondering. I want to make sure we are not violating any implicit assumptions for the user

[jeremylt]

I think one of the biggest potential issues is the fact that the parallelization strategy for gen will struggle to be compatible with this proposal. We can likely get it to work, but it will have a lot of dead threads during some parts of the kernel execution if there is high disparity in the number or size of the elements.

Also, we would need to selectively disable some of the consistency checking at the interface level (in interface/ceed-operator.c).

I think really clear diagramming and descriptions that we would put into the documentation would be a good first step in helping specify the conceptual model and starting to see what functions we would need to add to the API. (Good start above) Perhaps on the user side just a CeedOperatorSetForeignField(op, field_name, rstr_foreign, rstr_native, basis, vec), but lots of work internally to make such a think work? (spitballing here, but we should do this in the order conceptual model -> API -> impl, I think, since this is a big changeup in assumptions)

I do think this only really makes sense 'AtPoints', which gives us more room to breathe vs the current conceptual model.