Fix interpolation onto submesh with maps requiring dof transformations

[jorgensd]

Two bug fixes:

Interpolation onto submesh with same element type

Currently, interpolation from space A to space B doesn't use permutation info if element_A == element_B.
However, this is not correct when working with submeshes, as the cell_permutation_data might be different between the meshes.

The initial issue was discovered by @afasil9 when working on TEAM30 (https://github.com/Wells-Group/TEAM30/).

This pull request add a test that fails on 2,3 and 4 processors on main and a fix to interpolate.

Interpolation onto same function space (and mesh), with subset of cells

If the number of cells (cells0.size()) passed into u.interpolate(v, cells0=cells0)
matches the number of cells on the process (including ghosts), we currently assume that cells0=np.arange(num_cells_on_process), which is not True, as cells0 can contain duplicate cells.

This is fixed by adding a progressive set of checks to determine if this is True.

[garth-wells]

Should there be changes in the core code too? I only see a new test in the PR.

[jorgensd]

Should there be changes in the core code too? I only see a new test in the PR.

Yes, just waiting for the CI to fail before I push the fix.

[jorgensd]

Should there be changes in the core code too? I only see a new test in the PR.

Fix pushed.

[garth-wells]

Great to have this fixed. Suggested change of order for the if statement.

[garth-wells]

Looking further up in the function, we have

auto mesh = u1.function_space()->mesh();
...
auto cell_map0 = mesh->topology()->index_map(mesh->topology()->dim());

This looks a bit wrong (in a way related to what this PR fixes). We have mesh from u1, but then use mesh to extract cell_map0.

I'd suggest changing

auto mesh = u1.function_space()->mesh();

to

auto mesh1 = u1.function_space()->mesh();

and checking for foo0 vs foo1 consistency in the function.

This would also make the fix in this PR clearer with:

 if ((mesh1 == u0.function_space()->mesh()) and (element1 == element0 or *element1 == *element0))

[jorgensd]

Looking further up in the function, we have

auto mesh = u1.function_space()->mesh();
...
auto cell_map0 = mesh->topology()->index_map(mesh->topology()->dim());

This looks a bit wrong (in a way related to what this PR fixes). We have mesh from u1, but then use mesh to extract cell_map0.

I'd suggest changing

auto mesh = u1.function_space()->mesh();

to

auto mesh1 = u1.function_space()->mesh();

and checking for foo0 vs foo1 consistency in the function.

This would also make the fix in this PR clearer with:

 if ((mesh1 == u0.function_space()->mesh()) and (element1 == element0 or *element1 == *element0))

Clarification pushed.
There is still the underlying assumption

  auto mesh1 = u1.function_space()->mesh();
  assert(mesh1);
  assert(cells0.size() == cells1.size());

  auto cell_map1 = mesh1->topology()->index_map(mesh1->topology()->dim());
  assert(cell_map1);
  std::size_t num_cells1 = cell_map1->size_local() + cell_map1->num_ghosts();
  if (u1.function_space() == u0.function_space()
      and cells0.size() == num_cells1)

which in theory could be wrong, as cells0 = {0, 1, 3, 3, 2, 2} and num_cells1 = 6; would assume that cells0 = {0,1,2,3,4,5};.
It is not a very likely use-case, but an awful one to debug. One would have to add the check

auto mismatch = std::mismatch(cells0.begin(),cells0.end(), cells1.begin());
if (u1.function_space() == u0.function_space()
    and cells0.size() == num_cells1 and mismatch.first()==cells0.end())

Stand alone example of proper check (tested in godbolt):

#include<vector>
#include<algorithm>
#include<cstdint>
#include<iostream>
#include<iterator>
#include<ranges> 
int main ()
{

    std::vector<std::int32_t> a = {0,1,2,3,5,4};
    std::vector<std::int32_t> c = {0,1,2,4,4,5};
    [[maybe_unused]] auto [first, last] = std::mismatch(a.begin(), a.end(), c.begin());
    [[maybe_unused]] auto [first_iota, last_iota] = std::mismatch(a.begin(), a.end(), std::ranges::iota_view{0, (int)a.size()}.begin());
    std::cout <<  std::distance(a.begin(), first) << "\n";
    std::cout << std::distance(a.begin(), first_iota) << "\n";
    return 1;
}

which returns

3
4

[jorgensd]

@garth-wells I've added a progressive check that tests if:

1. if same function space
    2. If (1.) check if cells0.size() == num_cells_local + num_ghost_cells
        3. if (2.) check that cells0 == range(num_cells_local + num_ghost_cells) (where cells0 can be unsorted)

if 3. is true 
   do direct copy of the underlying array
else:
   Make further checks on the equality of elements, map type, same mesh, etc.

A check to test this has been added to the code.