# Interpolating between meshes - Python Binding **URL:** \ **Category:** General\ **Created:** [April 13, 2023, 10:39am UTC](https://fenicsproject.discourse.group/t/interpolating-between-meshes-python-binding/10922 "2023-04-13T10:39:42Z")\ **Posts on this page:** 16\ **Page:** 1
**Author:** ![Khuldoon\_Usman](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/khuldoon_usman/32/4884_2.png) [@Khuldoon\_Usman](https://fenicsproject.discourse.group/u/Khuldoon_Usman)\ **Post date:** [April 13, 2023, 10:39am UTC](https://fenicsproject.discourse.group/t/interpolating-between-meshes-python-binding/10922/1 "2023-04-13T10:39:42Z")
Hello All, I am trying to interpolate/project the a function from one mesh to another mesh but somehow the use of projection or just using interpolate is not working. Then I found this solution: > > > This pull request supersedes #1598 . > This pull request implements the same algo…rithm using the MPI Neighbourhood communication paradigm instead of the MPI One-sided communication one. The paradigm switch is joint work with @IgorBaratta . But this has no Python Bindings created as of yet. How to use the above solution in Python? Best, Khuldoon ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [April 13, 2023, 11:06am UTC](https://fenicsproject.discourse.group/t/interpolating-between-meshes-python-binding/10922/2 "2023-04-13T11:06:38Z")
See [dolfinx/test\_function.py at 160ed13eb476df99944072aec70bd46a6fcb9bcf · FEniCS/dolfinx · GitHub](https://github.com/FEniCS/dolfinx/blob/160ed13eb476df99944072aec70bd46a6fcb9bcf/python/test/unit/fem/test_function.py#L199) ---
**Author:** ![Khuldoon\_Usman](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/khuldoon_usman/32/4884_2.png) [@Khuldoon\_Usman](https://fenicsproject.discourse.group/u/Khuldoon_Usman)\ **Post date:** [April 13, 2023, 11:15am UTC](https://fenicsproject.discourse.group/t/interpolating-between-meshes-python-binding/10922/3 "2023-04-13T11:15:20Z")
```auto from dolfinx.fem import create_nonmatching_meshes_interpolation_data ``` ```auto ImportError: cannot import name 'create_nonmatching_meshes_interpolation_data' from 'dolfinx.fem' (/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/ __init__.py) ``` Note: I am running dolfinx via the container ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [April 13, 2023, 12:14pm UTC](https://fenicsproject.discourse.group/t/interpolating-between-meshes-python-binding/10922/4 "2023-04-13T12:14:20Z")
What version of DOLFINx are you running? What is the output of `python3 -c "import dolfinx;print(dolfinx. __version__ )"` ---
**Author:** ![Khuldoon\_Usman](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/khuldoon_usman/32/4884_2.png) [@Khuldoon\_Usman](https://fenicsproject.discourse.group/u/Khuldoon_Usman)\ **Post date:** [April 13, 2023, 12:16pm UTC](https://fenicsproject.discourse.group/t/interpolating-between-meshes-python-binding/10922/5 "2023-04-13T12:16:19Z")
> [@dokken](#): > > python3 -c "import dolfinx;print(dolfinx. **version** ) I have version 0.6.0.0 although I know that there is a higher version out (0.7.0.0). However, I do not know if there is a docker image out there of this version. ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [April 13, 2023, 1:26pm UTC](https://fenicsproject.discourse.group/t/interpolating-between-meshes-python-binding/10922/6 "2023-04-13T13:26:47Z")
This is the v0.6.0 way of doing it [dolfinx/test\_function.py at v0.6.0 · FEniCS/dolfinx · GitHub](https://github.com/FEniCS/dolfinx/blob/v0.6.0/python/test/unit/fem/test_function.py#L198) ---
**Author:** ![Khuldoon\_Usman](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/khuldoon_usman/32/4884_2.png) [@Khuldoon\_Usman](https://fenicsproject.discourse.group/u/Khuldoon_Usman)\ **Post date:** [April 14, 2023, 8:44am UTC](https://fenicsproject.discourse.group/t/interpolating-between-meshes-python-binding/10922/7 "2023-04-14T08:44:21Z")
Hello Dokken, I tried this but that gives me the wrong result. I will now share the full code and try to explain my approach. I have extracted Point Data from VTU file and mesh conversion via Meshio and then interpolated the Point Data into a Fenics function. Now I want to interpolate this function (called `point_function` in the code) into a seperate mesh and for this I created another function in this mesh (called `T`). ```auto ! pip3 install pandas ! pip3 install meshio ! pip3 install h5py ! pip install pyvista ! pip install panel import numpy as np import pyvista import panel import meshio import ufl import dolfinx as dfx import numpy as np from dolfinx.io import XDMFFile from dolfinx import plot from mpi4py.MPI import COMM_WORLD # use meshio to read vtk mesh filename = "new.vtu" vtk_mesh = meshio.vtu.read(filename) # cell data field as named in vtu (see paraview) fieldname = "Temperature" # use meshio to convert to and write a dolfinx mesh dim = 3 ps=vtk_mesh.points[:,:dim] cs = vtk_mesh.get_cells_type("tetra") dolfinx_mesh = meshio.Mesh(points=ps, cells={"tetra": cs}) meshio.write("test.xdmf", dolfinx_mesh) # read mesh with dolfinx comm = COMM_WORLD with XDMFFile(comm, "test.xdmf", "r") as file: dolfinx_mesh = file.read_mesh(name="Grid") # create the necessary grid functions cell_element=ufl.FiniteElement("DG", dolfinx_mesh.ufl_cell(), 0) point_element=ufl.FiniteElement("Lagrange", dolfinx_mesh.ufl_cell(), 1) cell_function_space=dfx.fem.FunctionSpace(dolfinx_mesh, cell_element) point_function_space=dfx.fem.FunctionSpace(dolfinx_mesh, point_element) cell_function = dfx.fem.Function(cell_function_space) point_function = dfx.fem.Function(point_function_space) point_function.name = fieldname # Global indexes owned locally point_ids = dolfinx_mesh.geometry.input_global_indices #cell_ids = dolfinx_mesh.geometry.input_cell_indices # Map VTK data directy onto dolfinx vectors #cell_function.x.array[:] = vtk_mesh.cell_data[fieldname][ids] point_function.x.array[:] = vtk_mesh.point_data[fieldname][point_ids] # create new mesh test_mesh = dfx.mesh.create_box(COMM_WORLD, [np.array([0.4, -0.125, -0.3]), np.array([0.9, 0.125, 0.0])], [50, 25, 25], cell_type=dfx.mesh.CellType.tetrahedron) tdim = test_mesh.topology.dim fdim = tdim - 1 test_mesh.topology.create_connectivity(fdim, tdim) boundary_facets = dfx.mesh.exterior_facet_indices(test_mesh.topology) pyvista.set_jupyter_backend("panel") pyvista.start_xvfb() topology, cell_types, geometry = plot.create_vtk_mesh(test_mesh, tdim) grid = pyvista.UnstructuredGrid(topology, cell_types, geometry) plotter = pyvista.Plotter() plotter.add_mesh(grid, show_edges=True) plotter.view_xy() if not pyvista.OFF_SCREEN: plotter.show() else: figure = plotter.screenshot("fundamentals_mesh.png") #create new function in new mesh V0 = ufl.FiniteElement("Lagrange", test_mesh.ufl_cell(), 1) V = dfx.fem.FunctionSpace(test_mesh, V0) T = dfx.fem.Function(V) # this gives an error from dolfinx.fem import create_nonmatching_meshes_interpolation_data # this is what dokken suggests for interpolating function on old mesh in new mesh fo 0.6.0 dolfinx T.interpolate(point_function) T.x.scatter_forward() u_topology, u_cell_types, u_geometry = plot.create_vtk_mesh(V) u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry) u_grid.point_data["T"] = T.x.array.real u_grid.set_active_scalars("T") u_plotter = pyvista.Plotter() u_plotter.add_mesh(u_grid, show_edges=True, cmap="rainbow") u_plotter.view_xy() if not pyvista.OFF_SCREEN: # result is not the same u_plotter.show() p_topology, p_cell_types, p_geometry = plot.create_vtk_mesh(point_function_space) p_grid = pyvista.UnstructuredGrid(p_topology, p_cell_types, p_geometry) p_grid.point_data["T"] = point_function.x.array.real p_grid.set_active_scalars("T") p_plotter = pyvista.Plotter() p_plotter.add_mesh(p_grid, show_edges=True, cmap="rainbow") p_plotter.view_xy() if not pyvista.OFF_SCREEN: # this is what it should look like p_plotter.show() ``` Original Function: ![Bild2](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/f/f47efecb96498736d980955c44395cfba2ba9954.jpeg) Function (on new mesh) into which it should be interpolated - currently not being properly interpolated: ![Bild1](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/e/ec9f6156a68d98a4ce693667bb7b753938ef8868.jpeg) One thing I should have mentioned before is that the original mesh is unstructured and the new one is structured. Does that cause an issue? ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [April 14, 2023, 9:37am UTC](https://fenicsproject.discourse.group/t/interpolating-between-meshes-python-binding/10922/8 "2023-04-14T09:37:27Z")
> [@Khuldoon\_Usman](#): > > Now I want to interpolate this function (called `point_function` in the code) Have you verified that `point_function` has the correct results, because I doubt that > [@Khuldoon\_Usman](#): > > `point_function.x.array[:] = vtk_mesh.point_data[fieldname][point_ids]` is correct, as you are assuming a 1-1 mapping between the mesh vertices of the input mesh in meshio and the degrees of freedom in the DOLFINx mesh (read by XDMF). ---
**Author:** ![Khuldoon\_Usman](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/khuldoon_usman/32/4884_2.png) [@Khuldoon\_Usman](https://fenicsproject.discourse.group/u/Khuldoon_Usman)\ **Post date:** [April 14, 2023, 1:30pm UTC](https://fenicsproject.discourse.group/t/interpolating-between-meshes-python-binding/10922/9 "2023-04-14T13:30:59Z")
Well, I am a little confused by what this means. The plot of point function seems to be valid and we checked this by saving the meshio mesh once again as an XDMF (with point function included) and compared it with the original VTU file in Paraview. The mesh on the other hand seems to be dest ```auto # Sanity check (use paraview to compare the data). with XDMFFile(COMM_WORLD, "sanity_check.xdmf", "w") as xdmf: xdmf.write_mesh(dolfinx_mesh) xdmf.write_function(point_function) ``` A comparison file is attached. And the imported meshio mesh also looks a little bit jadded. But still, how would this affect the interpolation on the new mesh? ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/0/094bdeaf3643abe0eb3dbadad3dee08df26ba953.jpeg) ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [April 14, 2023, 1:54pm UTC](https://fenicsproject.discourse.group/t/interpolating-between-meshes-python-binding/10922/10 "2023-04-14T13:54:54Z")
> [@Khuldoon\_Usman](#): > > ```auto > # Sanity check (use paraview to compare the data). > with XDMFFile(COMM_WORLD, "sanity_check.xdmf", "w") as xdmf: > xdmf.write_mesh(dolfinx_mesh) > xdmf.write_function(point_function) > > ``` > > A comparison file is attached. And the imported meshio mesh also looks a little bit jadded. But still, how would this affect the interpolation on the new mesh? First of all, what you are doing works in serial on this given mesh, and is not guaranteed to work on every mesh. Secondly, it looks like your mesh does not properly contain the mesh you are interpolating from. The current interpolation operator does not do extrapolation, i.e. Any degree of freedom in the mesh you are Interpolating to whose interpolation point is not in the original mesh, will be set to 0. Which is most likely why you see the behavior that the function jumps up and down. ---
**Author:** ![Khuldoon\_Usman](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/khuldoon_usman/32/4884_2.png) [@Khuldoon\_Usman](https://fenicsproject.discourse.group/u/Khuldoon_Usman)\ **Post date:** [April 14, 2023, 2:22pm UTC](https://fenicsproject.discourse.group/t/interpolating-between-meshes-python-binding/10922/11 "2023-04-14T14:22:12Z")
I see. This definitely makes a lot of sense, so then how do you suggest to extract point data and mesh data from a VTU and interpolate onto a dolfinx mesh as well as a dolfinx function. The original point data comes from a “box clip” where only a small bounding box is extracted from the original VTU file which is the reason for the extrapolation issue. My original VTU file only has constant cell data and thus through the “point data to cell data” filter I got the point data and then clipped it into a smaller box before importing it into the above script. ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [April 14, 2023, 3:05pm UTC](https://fenicsproject.discourse.group/t/interpolating-between-meshes-python-binding/10922/12 "2023-04-14T15:05:58Z")
First of all, i would definitely use cell data, and not a point to cell data filter (as this changes the variable from DG0 to CG1). Secondly, why have you clipped the original mesh. ---
**Author:** ![Khuldoon\_Usman](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/khuldoon_usman/32/4884_2.png) [@Khuldoon\_Usman](https://fenicsproject.discourse.group/u/Khuldoon_Usman)\ **Post date:** [April 17, 2023, 9:41am UTC](https://fenicsproject.discourse.group/t/interpolating-between-meshes-python-binding/10922/13 "2023-04-17T09:41:04Z")
All right, if I use cell data, then how does one interpolate that cell data to get the point data? Because I would like to use this data for another variation problem where this data will become a source term on a new mesh (hence the interpolation on another mesh). I clipped the mesh because I only wanted use this part of the mesh as my domain for my problem. Hence I used the paraview “Clip” filter to clip a small box and used it for import into fenics and thus subsequent projection/interpolation onto another mesh. This is also why I wanted to use point data (filtered from cell data) so that clipping any arbitrary part of the original domain would stell give me point values at the borders of the newly clipped domain. ---
**Author:** ![Khuldoon\_Usman](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/khuldoon_usman/32/4884_2.png) [@Khuldoon\_Usman](https://fenicsproject.discourse.group/u/Khuldoon_Usman)\ **Post date:** [April 17, 2023, 1:08pm UTC](https://fenicsproject.discourse.group/t/interpolating-between-meshes-python-binding/10922/14 "2023-04-17T13:08:17Z")
Also, I think I figure out the problem. The original VTK mesh had both 2D triangle elements as well as 3D tetrahedral elements and thus it is a mixed mesh. However since XDMF cannot handle mixed meshes, the meshio mesh was created using only tetra elements from the VTK file which therefore gives this jagged look in PyVista. Would it be possible to somehow handle the mixed mesh and import into a dolfinx function? ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [April 17, 2023, 1:23pm UTC](https://fenicsproject.discourse.group/t/interpolating-between-meshes-python-binding/10922/15 "2023-04-17T13:23:31Z")
No, that is not currently supported. There is some initial work on this: > > > This adds some rudimentary support for mixed topology/geometry. > At present, mu…ch of the functionality is at the 'stub' level, i.e. replacing: \`CellType\`, \`ElementDofLayout\` and \`CoordinateElement\` with lists instead of a single value for a \`Mesh\`. > The main new feature is the \`entity\_group\_offsets\` which defines the blocks of entities (i.e. cells or facets) which exist. > > Cells are kept in blocks in the topology, e.g. triangles, then quads, followed by ghost triangles and ghost quads. The offsets point to the start/end point of each block. But the IO from XDMFFile does not support it. Is there a reason your mesh is mixed-element? ---
**Author:** ![Khuldoon\_Usman](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/khuldoon_usman/32/4884_2.png) [@Khuldoon\_Usman](https://fenicsproject.discourse.group/u/Khuldoon_Usman)\ **Post date:** [April 17, 2023, 1:39pm UTC](https://fenicsproject.discourse.group/t/interpolating-between-meshes-python-binding/10922/16 "2023-04-17T13:39:21Z")
Well because the mesh is joint between two seperate meshes and thus their interconnectivity required surface elements in between