Object oriented "Projection" implementation

[javijv4]

Currently, svMultiPhysics has the option to link DOFs from different meshes using <Add_Projection>. The current implementation is used in two main applications

1. FSI problems (see tests/cases/fsi/pipe_3d for example) to match the inner-wall DOFs from the structure problem to the lumen wall DOFs from the fluid.
2. EP problem to match nodal values from one mesh to another. The main application is to match 1D nodes from a Purkinje representation to a 3D myocardium (see tests/cases/cep/cylinder_purkinje_1d3d).

Both applications match nodal positions between meshes and assign the same DOFs to both meshes. I believe the fluid and solid meshes are always conforming in the FSI case, so the nodal position match is direct. In EP, however, the 1D mesh nodes will not match the 3D mesh nodal positions, so the code finds the nearest-neighbor node, which changes the geometry (see #499).

The current implementation exists mostly in read_msh.cpp where the main set_projector() function calls the match_faces() function, where the matching of nodes occurs, and the shared DOFs are assigned.

Related desired features

• Multi-point constraint, MPC (Implement multipoint constraint for EP simulations #499): This is a direct alternative to the current implementation. Instead of snapping nodes to match nodal positions, it constrains the background mesh solution to match the value at the inner-node position from the projected mesh. This can be achieved by using a Lagrange Multiplier that can be solved for before the system solve and used to modify the RHS (again see Implement multipoint constraint for EP simulations #499). I have this implemented in this branch. This implementation should work with any physics, although I've only tested it with scalar equations (EP, Heat, Darcy).
• Purkinje-myocardium junctions (PMJ): This is a more EP-specific implementation. It couples the Purkinje end node with myocardial elements within a certain range through a series of resistors (see below, taken from this paper)
  .
  In EP, the amount of current you inject into the myocardium determines whether depolarization is triggered, making this approach more appropriate because it is not mesh-dependent (like the others).

These two more complex approaches require calling "projector" related functions during the solve, so not everything is contained in the presolve stage.

Proposed implementation

Create a Projector class (not sure about the name) that handles the setup of these projections (like the match_faces() function) and contains all the required functions that can be called during the solve if needed.
This branch has a (mostly vibe-coded) object-oriented implementation of my old MPC implementation. I can extend it to the other projection types (and clean up the implementation), but I have a few questions.

Questions

1. Should this class handle the Face projection as well (as used by FSI)? This face projection is basically also doing a node-by-node matching, so I think it is analogous and not complicated at all. This would mean that the proposed class would handle

• Face-to-face projection (by nearest neighbor)
• Node-to-face projection (by nearest neighbor)
• MPC
• PMJ

2. Who should own this class? In the OO implementation of the MPC, com_mod owns the Projector object, so it is easy to access from anywhere, but I am not sure this is the best option.
3. Lastly, before everything above, do you think it is worth creating such a class that encompasses the old projection methods, or should I focus on adding OO implementations for MPC and PMJ only?

[michelebucelli]

Thanks for describing this @javijv4!

I think that all the applications you describe share a requirement for spatial queries (e.g. find the mesh node closest to a given point, or find the point in the mesh - including element interiors - closest to a given point, or find the mesh nodes that are within a certain distance of a given point). Similar tasks are implemented for URIS (uris_find_closest_face_centroid and find_closest_element_centroid from uris.cpp), remeshing (interp in remesh.cpp), and possibly a few other places.

I suggest that that is the first thing that can be abstracted. This could be done through a class, say PointLocator, that is constructed from mesh information, stores that in a query-efficient format, and exposes a few query methods (e.g. PointLocator::find_nearest_neighbor, PointLocator::find_containing_element, and so on). This class would be a drop-in replacement for the places in the code where nearest-neighbor search is implemented (according to Claude, this is currently happening in 6 independent places).

(As a bonus: those nearest-neighbor searches are mostly linear searches right now, which is possibly quite inefficient. Once this is centralized in a single class, it becomes easier to optimize this operation through e.g. spatial partitioning trees.)

Having done that, we could look into the similarities (if any) between the different kinds of couplings that rely on this. I think that the four use cases that you describe look different enough to warrant different classes or functions, if nothing else to clearly communicate projection/interpolation semantics. However, a lot of the logic can probably be abstracted away to the spatial query helper.

(Partially related, I think that for the EP application matching the nodes of the 1D and 3D meshes leads to an ill-posed problem, but I've explained this better here.)

[javijv4]

Sounds good! I think it makes sense to implement first a class that deals with the spatial search and query first. I will open an issue about it soon.

[michelebucelli]

For efficient spatial queries, in the past I have used R-trees from boost::geometry, as wrapped by deal.II. The library ArborX offers a similar functionality, but I have little experience with it. I think it is designed for parallel queries and various forms of distributed hardware, while boost is maybe a bit more basic in that respect.

Anyways, these optimizations can come after the basic infrastructure for spatial queries is established.

[ktbolt]

I agree with first creating classes for spatial search. Both the uris and remeshing code use linear search for point location.

Kd-trees and Octrees are typically used for 3D point location. We might could try using VTK for this, maybe vtkOctreePointLocator.