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:
But this has no Python Bindings created as of yet. How to use the above solution in Python?
Best,
Khuldoon
from dolfinx.fem import create_nonmatching_meshes_interpolation_data
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
What version of DOLFINx are you running?
What is the output of 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.
This is the v0.6.0 way of doing it dolfinx/test_function.py at v0.6.0 · FEniCS/dolfinx · GitHub
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).
! 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:
Function (on new mesh) into which it should be interpolated - currently not being properly interpolated:
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?
Have you verified that point_function has the correct results, because I doubt that
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).
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
# 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.
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.
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.
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.
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?
No, that is not currently supported. There is some initial work on this:
But the IO from XDMFFile does not support it.
Is there a reason your mesh is mixed-element?
Well because the mesh is joint between two seperate meshes and thus their interconnectivity required surface elements in between
