Add functionalities to cache data when interpolating on nonmatching meshes

[massimiliano-leoni]

After merging #2245 dolfinx can interpolate discrete functions defined on different meshes. This is, however, quite costly.
This PR adds a class that performs this interpolation caching as much data as possible so that interpolating again between the same two functions, when their meshes haven't moved, is much cheaper.

Some timings for the same example as in the discussion of #2245 [mesh with 750'000 cells]

serial
Interpolating       7.311494
Interpolating again 0.572593

n = 6 processes
Interpolating       1.772976
Interpolating again 0.164332

[jhale]

Would it be possible to implement this not using a class? For example, the interpolate method could return a map containing the cache which could then be passed as an optional argument. std::optional might help.

[jhale]

Another (perhaps better) option would be to bring the expensive call to determine_point_ownership outside the function.

[IgorBaratta]

Another (perhaps better) option would be to bring the expensive call to determine_point_ownership outside the function.

There's also opportunity of improvement in eval, which is quite expensive (evaluating basis at interpolation points).
Maybe, for starters, we could use a more explicit way of caching data, with the functions that are already exist (in a demo for example).

[massimiliano-leoni]

@jhale I guess it is possible to avoid using a class. Out of ignorance curiosity, what is the matter with using a class?

[jhale]

@massimiliano-leoni A core design principles in DOLFINx is to avoid implicit state (as much as possible) and to prefer functional approaches as much as possible. I think what you want to achieve could be more straightforwardly achieved by making the expensive call to determine_point_ownership explicit.

[massimiliano-leoni]

@jhale I understand that it's a design decision, but I was wondering for my own knowledge what the reason behind the decision is!

[garth-wells]

@jhale I understand that it's a design decision, but I was wondering for my own knowledge what the reason behind the decision is!

Avoiding, where possible, objects with (mutable) state leads to simpler code (not need to check state), avoids inconsistent state, avoids hidden memory 'leaks' (not technically leaks, but objects that might hold a lot of memory and which isn't obvious to the user), and thread safety.

This is why we prefer explicit approaches first. Sometimes it might then be appropriate to build a convenience class/function on top of the low-level implementation.

For this particular case, it seems simplest to me that the user calls geometry::determine_point_ownership and then uses the data. This give more flexibility, e.g. a case where the caller knows the cells and just needs to determine the ranks.

[massimiliano-leoni]

@garth-wells Thanks for the explanation. I implemented it the way you and @jhale suggested. Is it better now?

[jhale]

Pretty much I think, just a few style things I left in the review. Another thing we do is try to make expensive calls explicit, so I don't think we would include a code path where the data is not pre-computed in this low-level interface.

[jhale]

I think this is close. How does it all look from the Python interface?

[massimiliano-leoni]

I think this is close. How does it all look from the Python interface?

I think I updated it accordingly. At least, the CI is satisfied. I suggest someone quickly goes through my additions because I'm not very familiar with pybind11 and writing Python wrappers so I might have overlooked something.