Mixed function interpolation between different mesh

nikoleeaa · April 4, 2025, 7:40pm

I am trying to interpolate data from a mixed finite element function space on one mesh to another function space on a finer mesh. However, my code fails when attempting to interpolate using f1.sub(0).interpolate(f2.sub(0), nmm_interpolation_data=...), while a similar setup works fine for scalar function spaces.

Below is a minimal reproducible example: (dolfinx version is 0.7.2)

from mpi4py import MPI
from dolfinx import fem
from dolfinx.mesh import create_rectangle
from dolfinx.io import XDMFFile
import ufl


domain1 = create_rectangle(MPI.COMM_WORLD, [(0,0),(12.5e-3,12.5e-3)], (10,10))
domain2 = create_rectangle(MPI.COMM_WORLD, [(0,0),(12.5e-3,12.5e-3)], (200,200))

e1 = ufl.FiniteElement("Lagrange", domain1.ufl_cell(), 1)
element_V1 = ufl.MixedElement([e1, e1])
V1 = fem.FunctionSpace(domain1, element_V1)

e2 = ufl.FiniteElement("Lagrange", domain2.ufl_cell(), 1)
element_V2 = ufl.MixedElement([e2, e2])
V2 = fem.FunctionSpace(domain2, element_V2)

f1 = fem.Function(V1)
f2 = fem.Function(V2)

f1.sub(0).interpolate(lambda x: x[0]*2+x[1])

f2.sub(0).interpolate(f1.sub(0).collapse(), nmm_interpolation_data=fem.create_nonmatching_meshes_interpolation_data(
        f1.sub(0).collapse().function_space.mesh._cpp_object,
        f1.sub(0).collapse().function_space.element,
        f2.sub(0).collapse().function_space.mesh._cpp_object,
        padding = 1e-14))

with XDMFFile(MPI.COMM_WORLD, "data2D.xdmf", "w") as ufile_xdmf:
            ufile_xdmf.write_mesh(domain2)
            ufile_xdmf.write_function(f2.sub(0))

There are the following error prompts：

Loguru caught a signal: SIGSEGV
Stack trace:
26      0x5624c9f0bf25 _start + 37
25      0x7fd317e67e40 __libc_start_main + 128
24      0x7fd317e67d90 /lib/x86_64-linux-gnu/libc.so.6(+0x29d90) [0x7fd317e67d90]
23      0x5624c9f0c02d Py_BytesMain + 45
22      0x5624c9f35a5e Py_RunMain + 702
21      0x5624c9f42f83 _PyRun_AnyFileObject + 67
20      0x5624c9f43338 _PyRun_SimpleFileObject + 424
19      0x5624c9f43e55 /bin/python3(+0x25fe55) [0x5624c9f43e55]
18      0x5624c9f3d9cb /bin/python3(+0x2599cb) [0x5624c9f3d9cb]
17      0x5624c9f44108 /bin/python3(+0x260108) [0x5624c9f44108]
16      0x5624c9f19256 PyEval_EvalCode + 134
15      0x5624c9e239c6 /bin/python3(+0x13f9c6) [0x5624c9e239c6]
14      0x5624c9e2853c _PyEval_EvalFrameDefault + 6540
13      0x5624c9e4c7f1 /bin/python3(+0x1687f1) [0x5624c9e4c7f1]
12      0x5624c9e2726d _PyEval_EvalFrameDefault + 1725
11      0x5624c9e3e9fc _PyFunction_Vectorcall + 124
10      0x5624c9e295d7 _PyEval_EvalFrameDefault + 10791
9       0x5624c9e3e9fc _PyFunction_Vectorcall + 124
8       0x5624c9e2ccfa _PyEval_EvalFrameDefault + 24906
7       0x5624c9e4cacb /bin/python3(+0x168acb) [0x5624c9e4cacb]
6       0x5624c9e34a7b _PyObject_MakeTpCall + 603
5       0x5624c9e3e10e /bin/python3(+0x15a10e) [0x5624c9e3e10e]
4       0x7fd30c0d55dc /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0xb25dc) [0x7fd30c0d55dc]
3       0x7fd30c19afc4 /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0x177fc4) [0x7fd30c19afc4]
2       0x7fd30c17eb2c /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0x15bb2c) [0x7fd30c17eb2c]
1       0x7fd30c34e567 /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0x32b567) [0x7fd30c34e567]
0       0x7fd317e80520 /lib/x86_64-linux-gnu/libc.so.6(+0x42520) [0x7fd317e80520]
2025-04-04 19:33:29.024 (   0.297s) [main            ]                       :0     FATL| Signal: SIGSEGV
Segmentation fault (core dumped)

If I use scalar function spaces instead of MixedElement as below, the interpolation works without issues:

from mpi4py import MPI
from dolfinx import fem
from dolfinx.mesh import create_rectangle
from dolfinx.io import XDMFFile

domain = create_rectangle(MPI.COMM_WORLD, [(0,0),(12.5e-3,12.5e-3)], (500,500))
domain1 = create_rectangle(MPI.COMM_WORLD, [(0,0),(12.5e-3,12.5e-3)], (200,200))

V = fem.functionspace(domain, ("Lagrange", 1))
V1 = fem.functionspace(domain1, ("Lagrange", 1))

f = fem.Function(V)
f1 = fem.Function(V1)

f.interpolate(lambda x: x[0]*2+x[1])
f1.interpolate(f, nmm_interpolation_data=fem.create_nonmatching_meshes_interpolation_data(
        f1.function_space.mesh._cpp_object,
        f1.function_space.element,
        f.function_space.mesh._cpp_object,
        padding = 1e-14))

with XDMFFile(MPI.COMM_WORLD, "data2D.xdmf", "w") as ufile_xdmf:
            ufile_xdmf.write_mesh(domain1)
            ufile_xdmf.write_function(f1)

Any guidance would be greatly appreciated.

dokken · April 4, 2025, 7:45pm

you are using an old version of DOLFINx (v0.7.2). Is the same error present on a newer version of DOLFINx?
If your «mixed element» is just copies of a scalar element, why not use a VectorElement with two components?

nikoleeaa:

element_V1 = ufl.MixedElement([e1, e1])
V1 = fem.FunctionSpace(domain1, element_V1)

e2 = ufl.FiniteElement("Lagrange", domain2.ufl_cell(), 1)
element_V2 = ufl.MixedElement([e2, e2])
V2 = fem.FunctionSpace(domain2, element_V2)

nikoleeaa · April 4, 2025, 7:51pm

Thank you so much for your quick response, @dokken!

I haven’t tried running it on the latest version of DOLFINx yet.

Will this improve? I’ll give it a try.

dokken · April 4, 2025, 8:06pm

VectorElement has a very different implementation under the hood, and should be easier to use for interpolation.

nikoleeaa · April 4, 2025, 8:07pm

from mpi4py import MPI
from dolfinx import fem
from dolfinx.mesh import create_rectangle
from dolfinx.io import XDMFFile


domain1 = create_rectangle(MPI.COMM_WORLD, [(0,0),(12.5e-3,12.5e-3)], (10,10))
domain2 = create_rectangle(MPI.COMM_WORLD, [(0,0),(12.5e-3,12.5e-3)], (200,200))

V1 = fem.FunctionSpace(domain1, ("Lagrange", 1, (2,))) 
V2 = fem.FunctionSpace(domain2, ("Lagrange", 1, (2,))) 

f1 = fem.Function(V1)
f2 = fem.Function(V2)

f1.sub(0).interpolate(lambda x: x[0]*2+x[1])

f2.sub(0).interpolate(f1.sub(0).collapse(), nmm_interpolation_data=fem.create_nonmatching_meshes_interpolation_data(
        f1.sub(0).collapse().function_space.mesh._cpp_object,
        f1.sub(0).collapse().function_space.element,
        f2.sub(0).collapse().function_space.mesh._cpp_object,
        padding = 1e-14))

with XDMFFile(MPI.COMM_WORLD, "data2D.xdmf", "w") as ufile_xdmf:
            ufile_xdmf.write_mesh(domain2)
            ufile_xdmf.write_function(f2.sub(0))

This seems to generate the same error message：

Loguru caught a signal: SIGSEGV
Stack trace:
26      0x55c5ef671f25 _start + 37
25      0x7fdb1b70ee40 __libc_start_main + 128
24      0x7fdb1b70ed90 /lib/x86_64-linux-gnu/libc.so.6(+0x29d90) [0x7fdb1b70ed90]
23      0x55c5ef67202d Py_BytesMain + 45
22      0x55c5ef69ba5e Py_RunMain + 702
21      0x55c5ef6a8f83 _PyRun_AnyFileObject + 67
20      0x55c5ef6a9338 _PyRun_SimpleFileObject + 424
19      0x55c5ef6a9e55 /bin/python3(+0x25fe55) [0x55c5ef6a9e55]
18      0x55c5ef6a39cb /bin/python3(+0x2599cb) [0x55c5ef6a39cb]
17      0x55c5ef6aa108 /bin/python3(+0x260108) [0x55c5ef6aa108]
16      0x55c5ef67f256 PyEval_EvalCode + 134
15      0x55c5ef5899c6 /bin/python3(+0x13f9c6) [0x55c5ef5899c6]
14      0x55c5ef58e53c _PyEval_EvalFrameDefault + 6540
13      0x55c5ef5b27f1 /bin/python3(+0x1687f1) [0x55c5ef5b27f1]
12      0x55c5ef58d26d _PyEval_EvalFrameDefault + 1725
11      0x55c5ef5a49fc _PyFunction_Vectorcall + 124
10      0x55c5ef58f5d7 _PyEval_EvalFrameDefault + 10791
9       0x55c5ef5a49fc _PyFunction_Vectorcall + 124
8       0x55c5ef592cfa _PyEval_EvalFrameDefault + 24906
7       0x55c5ef5b2acb /bin/python3(+0x168acb) [0x55c5ef5b2acb]
6       0x55c5ef59aa7b _PyObject_MakeTpCall + 603
5       0x55c5ef5a410e /bin/python3(+0x15a10e) [0x55c5ef5a410e]
4       0x7fdb0f97d5dc /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0xb25dc) [0x7fdb0f97d5dc]
3       0x7fdb0fa42fc4 /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0x177fc4) [0x7fdb0fa42fc4]
2       0x7fdb0fa26b2c /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0x15bb2c) [0x7fdb0fa26b2c]
1       0x7fdb0fbf6567 /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0x32b567) [0x7fdb0fbf6567]
0       0x7fdb1b727520 /lib/x86_64-linux-gnu/libc.so.6(+0x42520) [0x7fdb1b727520]
2025-04-04 20:07:07.941 (   0.285s) [main            ]                       :0     FATL| Signal: SIGSEGV
Segmentation fault (core dumped)

I also tried to replace f2. sub(0).interpolate() with the following:

f2.interpolate(f1, nmm_interpolation_data=fem.create_nonmatching_meshes_interpolation_data(
        f1.function_space.mesh._cpp_object,
        f1.function_space.element,
        f2.function_space.mesh._cpp_object,
        padding = 1e-14))

But it still doesn’t help.

nikoleeaa · April 4, 2025, 8:27pm

Dear dokken, I’m currently using dolfinx v0.7.2, and my development environment is not yet set up for the newer version. I’m wondering if you’d be willing to run the snippet on a newer version of DOLFINx just to confirm whether the issue still exists?

If it works on the latest version, I’ll prioritize migrating my code accordingly. Your help would mean a lot — really appreciate your time!

dokken · April 4, 2025, 8:33pm

Looking closer at your example: Are you trying to only interpolate a sub-space onto the other mesh?

dokken · April 4, 2025, 8:44pm

The following runs nicely for me on v0.9.0:

from mpi4py import MPI
from dolfinx import fem
from dolfinx.mesh import create_rectangle
from dolfinx.io import XDMFFile
import numpy as np


domain1 = create_rectangle(MPI.COMM_WORLD, [(0, 0), (12.5e-3, 12.5e-3)], (10, 10))
domain2 = create_rectangle(MPI.COMM_WORLD, [(0, 0), (12.5e-3, 12.5e-3)], (200, 200))

V1 = fem.functionspace(domain1, ("Lagrange", 1, (2,)))
V2 = fem.functionspace(domain2, ("Lagrange", 1, (2,)))

f1 = fem.Function(V1)
f2 = fem.Function(V2)

f1.sub(0).interpolate(lambda x: x[0] * 2 + x[1])


f10 = f1.sub(0).collapse()


V20, sub2_to_V2 = V2.sub(0).collapse()
f20 = fem.Function(V20)

c2_map = V20.mesh.topology.index_map(V20.mesh.topology.dim)
cells = np.arange(c2_map.size_local + c2_map.num_ghosts, dtype=np.int32)
nid = fem.create_interpolation_data(
    f20.function_space,
    f10.function_space,
    cells=cells,
    padding=1e-14,
)

f20.interpolate_nonmatching(
    f10,
    cells=cells,
    interpolation_data=nid,
)


f2.x.array[sub2_to_V2] = f20.x.array[:]
f2.x.scatter_forward()

with XDMFFile(MPI.COMM_WORLD, "data2D.xdmf", "w") as ufile_xdmf:
    ufile_xdmf.write_mesh(domain2)
    ufile_xdmf.write_function(f2.sub(0).collapse())

To me it seems like you have the order of the input arguments to nonmatching interpolation data the wrong way around, considering:

github.com/FEniCS/dolfinx

cpp/dolfinx/fem/interpolate.h

v0.7.2


      
          /// @brief Generate data needed to interpolate discrete functions across
          /// different meshes.
          /// @param[in] geometry0 Mesh geometry of the space to interpolate into
          /// @param[in] element0 Element of the space to interpolate into
          /// @param[in] mesh1 Mesh of the function to interpolate from
          /// @param[in] cells Indices of the cells in the destination mesh on
          /// which to interpolate. Should be the same as the list used when
          /// calling fem::interpolation_coords.
          /// @param[in] padding Absolute padding of bounding boxes of all entities on
          /// `mesh1`. This is used avoid floating point issues when an interpolation
          /// point from `mesh0` is on the surface of a cell in `mesh1`. This parameter
          /// can also be used for extrapolation, i.e. if cells in `mesh0` is not
          /// overlapped by `mesh1`.
          ///
          /// @note Setting the `padding` to a large value will increase runtime of this
          /// function, as one has to determine what entity is closest if there is no
          /// intersection.

github.com/FEniCS/dolfinx

python/dolfinx/wrappers/fem.cpp

v0.7.2


      
          "create_nonmatching_meshes_interpolation_data",
          [](const dolfinx::mesh::Mesh<T>& mesh0,
             const dolfinx::fem::FiniteElement<T>& element0,
             const dolfinx::mesh::Mesh<T>& mesh1, T padding)
          {
            int tdim = mesh0.topology()->dim();
            auto cell_map = mesh0.topology()->index_map(tdim);
            assert(cell_map);
            std::int32_t num_cells
                = cell_map->size_local() + cell_map->num_ghosts();
            std::vector<std::int32_t> cells(num_cells, 0);
            std::iota(cells.begin(), cells.end(), 0);
            return dolfinx::fem::create_nonmatching_meshes_interpolation_data(
                mesh0.geometry(), element0, mesh1,
                std::span(cells.data(), cells.size()), padding);
          },
          py::arg("mesh0"), py::arg("element0"), py::arg("mesh1"),
          py::arg("padding"));

nikoleeaa · April 5, 2025, 5:17am

Thank you so much! I carefully compared your working code in v0.9.0 with mine and tried adapting it to match the v0.7.2 API.

As you pointed out, it turns out the issue was caused by a very basic mistake — I had mixed up the order of the input arguments when calling create_nonmatching_meshes_interpolation_data.

Here’s the correct code that works with dolfinx v0.7.2：

from mpi4py import MPI
from dolfinx import fem
from dolfinx.mesh import create_rectangle
from dolfinx.io import XDMFFile

domain1 = create_rectangle(MPI.COMM_WORLD, [(0,0),(12.5e-3,12.5e-3)], (10,10))
domain2 = create_rectangle(MPI.COMM_WORLD, [(0,0),(12.5e-3,12.5e-3)], (200,200))

V1 = fem.FunctionSpace(domain1, ("Lagrange", 1, (2,))) 
V2 = fem.FunctionSpace(domain2, ("Lagrange", 1, (2,))) 

f1 = fem.Function(V1)
f2 = fem.Function(V2)

f1.sub(0).interpolate(lambda x: x[0]*2+x[1])

f2.sub(0).interpolate(f1.sub(0), nmm_interpolation_data=fem.create_nonmatching_meshes_interpolation_data(
        f2.sub(0).function_space.mesh._cpp_object,
        f2.sub(0).function_space.element,
        f1.sub(0).function_space.mesh._cpp_object,
        padding = 1e-14))

with XDMFFile(MPI.COMM_WORLD, "data2D.xdmf", "w") as ufile_xdmf:
            ufile_xdmf.write_mesh(domain2)
            ufile_xdmf.write_function(f2.sub(0))

As you can see, I confused the positions of f1 and f2 in the previous error code. Rookie error!

I found one more thing, when working with subfunctions, calling .sub(i) with or without .collapse() seems to give the same .function_space. That is to say, the above code can run well regardless of whether using .collapse(). Also, if vector space is not used, using the MixedElement space directly also works fine.

nikoleeaa · April 5, 2025, 5:22am

Yes, exactly — this was just a minimal example to demonstrate the issue, so it may look a bit artificial or oversimplified.

I really appreciate your thoughtful observation and kind reminder — thank you for taking the time to look into it so carefully, @dokken!

Topic		Replies	Views
Interpolation between different meshes in 0.7.1 General	1	235	November 6, 2023
Problem interpolating mixed function dolfinx dolfinx	11	2523	April 15, 2024
Change the mesh of a function dolfinx	10	1019	November 8, 2022
Interpolation on different grids General	3	39	September 15, 2024
Projection of DG0 Function over Non-Matching Meshes dolfinx interpolation	1	335	March 4, 2024

Mixed function interpolation between different mesh

Related topics