Cells argument of create_interpolation_data and interpolate_nonmatching

General
1

Concerning non-matchng mesh interpolation in the main branch, I only get expected behaviour if the cells argument of create_interpolation_data and interpolate_nonmatching is the same.

In particular, I can’t create the interpolation data for all the cells and then only interpolate on some cells.

Is this expected behaviour? Here is a MWE:

from dolfinx import io, fem, mesh
import ufl
from mpi4py import MPI
import dolfinx.fem.petsc
import numpy as np

comm = MPI.COMM_WORLD
rank = comm.Get_rank()

def interpolate(sending_func,
                receiving_func_mid,
                receiving_func_all,
                ):
    '''
    Interpolate sending_func to receiving_func,
    each comming from separate meshes
    '''
    rmesh = receiving_func_mid.function_space.mesh
    rtopology = rmesh.topology
    cmap = rtopology.index_map(rtopology.dim)
    num_cells = cmap.size_local + cmap.num_ghosts
    middle_nodes = mesh.locate_entities(rmesh,0,lambda x : np.isclose(x[0],0.5))
    rtopology.create_connectivity(0,2)
    middle_cells = mesh.compute_incident_entities(rtopology,middle_nodes,0,rtopology.dim)
    all_cells = np.arange(num_cells,dtype=np.int32)
    nmmid_mid = fem.create_interpolation_data(
                                 receiving_func_mid.function_space,
                                 sending_func.function_space,
                                 middle_cells,
                                 padding=0,)
    nmmid_all = fem.create_interpolation_data(
                                 receiving_func_mid.function_space,
                                 sending_func.function_space,
                                 all_cells,
                                 padding=0,)
    receiving_func_mid.interpolate_nonmatching(sending_func, middle_cells, interpolation_data=nmmid_mid)
    receiving_func_all.interpolate_nonmatching(sending_func, middle_cells, interpolation_data=nmmid_all)

def write_post(mesh1,mesh2,fs1,fs2):
    writer1, writer2 = io.VTKFile(mesh1.comm, f"out/result1.vtk",'w'),io.VTKFile(mesh2.comm, f"out/result2.vtk",'w')
    writer1.write_function(fs1)
    writer2.write_function(fs2)
    writer1.close()
    writer2.close()

def main():
    # Mesh and problems
    points_side = 16
    mesh1 = mesh.create_unit_square(MPI.COMM_WORLD, points_side, points_side, mesh.CellType.triangle)
    mesh2 = mesh.create_unit_square(MPI.COMM_WORLD, points_side, points_side, mesh.CellType.quadrilateral)

    (V1, V2) = (fem.functionspace(mesh1, ("Lagrange", 1)),fem.functionspace(mesh2, ("Lagrange", 1)))
    (u1, u2a, u2b)   = (fem.Function(V1),fem.Function(V2,name="mid_cells"),fem.Function(V2,name="all_cells"))
    # Solve Poisson problem on mesh1
    x = ufl.SpatialCoordinate(mesh1)
    u1.interpolate(fem.Expression(ufl.sin(x[1]),V1.element.interpolation_points()))
    # Interpolate from mesh1 to mesh2
    interpolate(u1,u2a,u2b)
    # Write post
    write_post(mesh1,mesh2,[u1],[u2a,u2b])

if __name__=="__main__":
    main()

Creating the interpolation data over the middle cells and interpolating to the middle cells:

image
image1265×950 43.4 KB

Creating the interpolation data over all cells and interpolating to the middle cells:
image
image1256×932 59.9 KB

2

Looking at dolfinx/python/test/unit/fem/test_interpolation.py at ab79530e1ae60a646002609b510e8392cc09e213 · FEniCS/dolfinx · GitHub, I don’t think we a test covering this case.

My understanding is that the cells argument in the data creation and interpolation call should be the same. If that is true, at the very least dolfinx/python/dolfinx/fem/function.py at ab79530e1ae60a646002609b510e8392cc09e213 · FEniCS/dolfinx · GitHub should be updated to include this information. Maybe wait a couple of days to give other developers time to confirm (or deny) that my understanding is correct, and if that is confirmed feel free to open a PR in the dolfinx repo with a suggestion on how you would write this documentation improvement.

1 Like