Quick update, my naive assignement to a mixed element was wrong. Building on this post I can now present a working code :
import meshio, ufl
import numpy as np
import dolfinx as dfx
from dolfinx.io import XDMFFile
from mpi4py.MPI import COMM_WORLD
from petsc4py import PETSc as pet
# Step 1 : read mesh and point data
in_mesh = meshio.read("input.xmf")
# Step 2 : write it out again
cells = mesh.get_cells_type("quad")
out_mesh = meshio.Mesh(points=in_mesh.points[:,:2], cells={"quad": cells})
meshio.write("mesh.xdmf", out_mesh)
# Step 3 : read it again in dolfinx
with XDMFFile(COMM_WORLD, "mesh.xdmf", "r") as file:
mesh = file.read_mesh(name="Grid")
# Step 4 : create FiniteElement, FunctionSpace & Functions
FE_vector_1=ufl.VectorElement("Lagrange",mesh.ufl_cell(),1,3)
FE_vector_2=ufl.VectorElement("Lagrange",mesh.ufl_cell(),2,3)
FE_scalar =ufl.FiniteElement("Lagrange",mesh.ufl_cell(),1)
V_Space=dfx.FunctionSpace(mesh,FE_vector_1)
U_Space=dfx.FunctionSpace(mesh,FE_vector_2)
p_Space=dfx.FunctionSpace(mesh,FE_scalar)
Space =dfx.FunctionSpace(mesh,FE_vector_2*FE_scalar)
V = dfx.Function(V_Space)
U = dfx.Function(U_Space)
p = dfx.Function(p_Space)
# Step 5 : jiggle indexing
idcs = np.argsort(mesh.geometry.input_global_indices).astype('int32')
vec_idcs = np.repeat(3*idcs,3)
vec_idcs[1::3]+=1
vec_idcs[2::3]+=2
# Step 6 : map one onto the other
V.vector[vec_idcs] = in_mesh.point_data['U'].flatten()
p.vector[idcs] = in_mesh.point_data['p']
# Step 7: interpolation
U.interpolate(V)
# Step 8 : write result as mixed
q = dfx.Function(Space)
_, map_U = Space.sub(0).collapse(collapsed_dofs=True)
_, map_p = Space.sub(1).collapse(collapsed_dofs=True)
q.vector[map_U]=U.vector
q.vector[map_p]=p.vector
viewer = pet.Viewer().createMPIIO("../cases/nozzle/baseflow/dat_complex/baseflow_S=0.000.dat", 'w', COMM_WORLD)
q.vector.view(viewer)
Also notice OpenFOAM v2112 & Paraview 5.10 output includes a entry that causes meshio to fail. If this line is removed in the .xmf file produced by Paraview however, above code runs smoothly.