Read OpenFoam results in dolfinx

dolfinx
11

Many thanks ! Final process looks like :

  1. Read the .xmf file using meshio
  2. Write the mesh again thanks to meshio, pruning out the z elements
  3. Read the mesh in dolfinx
  4. Use it to create FiniteElements of order 1 and associated FunctionSpace
  5. Jiggle around the data indexing
  6. Set the Function using vector directly
  7. Interpolate the lower order Function onto the higher order one
  8. Write the Functions into a mixed FunctionSpace and save

Save 8 is required for the rest of my code, but this whole process still feels a tad cumbersome (especially steps 2, 4). Nonetheless, it works, thanks again !

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("front3.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)
U_b,p_b = ufl.split(q)
U_b,p_b=U,p
viewer = pet.Viewer().createMPIIO("../cases/nozzle/baseflow/dat_complex/baseflow_S=0.000.dat", 'w', COMM_WORLD)
q.vector.view(viewer)