convert vtk to dolfinx

read_vtk.py

import numpy as np
import dolfinx
import pyvista as pv
import basix
import ufl
from dolfinx.mesh import create_mesh
from mpi4py import MPI

xrng = np.arange(0, 20, 5, dtype=np.float32)
yrng = np.arange(0, 20, 5, dtype=np.float32)
zrng = np.arange(0, 20, 5, dtype=np.float32)
x, y, z = np.meshgrid(xrng, yrng, zrng, indexing="ij")
grid = pv.StructuredGrid(x, y, z)

# add cell data
grid.cell_data["mean"] = np.arange(grid.n_cells)


# save to vtk file
grid.save("original.vtk")


# grid = pv.read("neutron_dose.vtk")

num_cells = grid.GetNumberOfCells()

# Extract connectivity information
cell_connectivity = []
points = []
ordering = [0, 1, 3, 2, 4, 5, 7, 6]

for i in range(num_cells):
    cell = grid.GetCell(i)  # Get the i-th cell
    point_ids = [cell.GetPointId(j) for j in ordering]  # Extract connectivity
    cell_connectivity.append(point_ids)

degree = 1  # Set polynomial degree
cell = ufl.Cell("hexahedron")
mesh_element = basix.ufl.element("Lagrange", cell.cellname(), degree, shape=(3,))

# Create dolfinx Mesh
mesh_ufl = ufl.Mesh(mesh_element)
dolfinx_mesh = create_mesh(MPI.COMM_WORLD, cell_connectivity, grid.points, mesh_ufl)


function_space = dolfinx.fem.functionspace(dolfinx_mesh, ("DG", 0))
u = dolfinx.fem.Function(function_space)

u.x.array[:] = grid.cell_data["mean"][dolfinx_mesh.topology.original_cell_index]

import dolfinx.io

writer = dolfinx.io.VTXWriter(MPI.COMM_WORLD, "u_as_dg0.bp", u, "BP5")
writer.write(t=0)

# to properly visualise in paraview we need to project the function to a continuous space
V_dg_1 = dolfinx.fem.functionspace(dolfinx_mesh, ("DG", 1))
u_dg_1 = dolfinx.fem.Function(V_dg_1)
u_dg_1.interpolate(u)
writer = dolfinx.io.VTXWriter(MPI.COMM_WORLD, "u_as_dg1.bp", u_dg_1, "BP5")
writer.write(t=0)