Not Properly Loading Solution Data From .xdmf File

I/O
1

I am trying to load in a .xdmf and .h5 file that have solution data. My code runs without errors but when I print out the velocity field it is only zeros. I have looked at the solution in Paraview, and my solver worked, the solution is nonzero. How do I correctly read in a function space/solution from a .xdmf file?

Below is my code.

import numpy as np

from skimage import io
import skimage as sk
import scipy.ndimage as ndimage
from rdp import rdp
import os

import dolfinx
import dolfinx.io
import dolfinx.fem
from mpi4py import MPI
import sys
import h5py

comm = MPI.COMM_WORLD
img_fname = sys.argv[1] # File name of input image
solname_xdmf = sys.argv[2] # File name of solution (.xdmf file)
solname_h5 = sys.argv[3] # File name of solution (.h5 file)

def read_sol_file(comm, sol_fname_xdmf, sol_fname_h5):
    with dolfinx.io.XDMFFile(comm, sol_fname_xdmf, "r") as infile:
        # Read the mesh
        mesh = infile.read_mesh(name="mesh")
        mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)

        # Create a function space (this needs to match the space used in the saved file)
        V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))

        # Create a function in the space
        u = dolfinx.fem.Function(V)

        # --- Load function values from HDF5 file ---
    with h5py.File(sol_fname_h5, "r") as h5f:
        # print("Available datasets in HDF5 file:", list(h5f.keys())) # Use to see datasets in the h5 file
        # print("Datasets inside 'Function':", list(group.keys()))  # Check available datasets
        data = h5f['Function']['Velocity']
        u.x.array[:] = data

        return[u, V, mesh]
    
u, V, mesh = read_sol_file(comm, solname_xdmf, solname_h5)

print(u.x.array)

I inspected the .xdmf file with vim, and it contains information about the velocity which, however I do not know how to load it in.

XdmfInspect
XdmfInspect1112×471 20.2 KB

I am using FEniCSx version 0.9, installed on WSL using conda.

2

What is the output of print(data)?

3

The output of print(data) is:
<HDF5 group "/Function/Velocity" (1 members)>

4

You need to actually access the data in that array.
Please note that the proposed method is not going to work in general, and you would have to consider something along the lines of:

5

Thank you for your help, the code almost works. When I run the code with I get an error on line 44 (assert len(mesh.geometry.cmaps) == 1 ), the error is: AttributeError: 'Geometry' object has no attribute 'cmaps'. Did you mean: 'cmap'?

When I change the line to read assert len(mesh.geometry.cmap) == 1 , it then gives the error TypeError: object of type 'CoordinateElement' has no len()

print(mesh.geometry.cmap) gives <dolfinx.fem.element.CoordinateElement object at 0x7f556658d550>

Here is the code I am using

import numpy as np

from rdp import rdp
import os

import dolfinx
import dolfinx.io
import dolfinx.fem
from mpi4py import MPI
import sys
import h5py

comm = MPI.COMM_WORLD
img_fname = sys.argv[1] # File name of input image
solname_xdmf = sys.argv[2] # File name of solution (.xdmf file)
solname_h5 = sys.argv[3] # File name of solution (.h5 file)

from pathlib import Path
import adios4dolfinx

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, solname_xdmf, "r") as xdmf:
    mesh = xdmf.read_mesh(name="mesh")

import basix
import h5py


def read_mesh_data(file_path:Path, mesh: dolfinx.mesh.Mesh, data_path:str):
    assert file_path.is_file(), f"File {file_path} does not exist"
    infile = h5py.File(file_path, "r")
    num_nodes_global = mesh.geometry.index_map().size_global
    assert data_path in infile.keys(), f"Data {data_path} does not exist"
    dataset = infile[data_path]
    shape = dataset.shape
    assert shape[0] == num_nodes_global, f"Got data of shape {shape}, expected {num_nodes_global, shape[1]}"
    dtype = dataset.dtype
    # Read data locally on each process
    local_input_range = adios4dolfinx.comm_helpers.compute_local_range(mesh.comm, num_nodes_global)    
    local_input_data = dataset[local_input_range[0]:local_input_range[1]]

    # Create appropriate function space (based on coordinate map)
    print(mesh.geometry.cmap)
    assert len(mesh.geometry.cmap) == 1, "Mixed cell-type meshes not supported"
    element = basix.ufl.element(
        basix.ElementFamily.P,
        mesh.topology.cell_name(),
        mesh.geometry.cmaps[0].degree,
        mesh.geometry.cmaps[0].variant,
        shape=(shape[1],),
        gdim=mesh.geometry.dim)

    # Assumption: Same dof layout for geometry and function space, cannot test in python    
    V = dolfinx.fem.FunctionSpace(mesh, element)
    uh = dolfinx.fem.Function(V, name=data_path)
    # Assume that mesh is first order for now
    assert mesh.geometry.cmaps[0].degree == 1, "Only linear meshes supported"
    x_dofmap = mesh.geometry.dofmap
    igi = np.array(mesh.geometry.input_global_indices, dtype=np.int64)
    global_geom_input = igi[x_dofmap]
    global_geom_owner = adios4dolfinx.utils.index_owner(mesh.comm, global_geom_input.reshape(-1), num_nodes_global)
    for i in range(shape[1]):
        arr_i = adios4dolfinx.comm_helpers.send_dofs_and_recv_values(global_geom_input.reshape(-1), global_geom_owner, mesh.comm, local_input_data[:,i], local_input_range[0])
        dof_pos = x_dofmap.reshape(-1)*shape[1]+i
        uh.x.array[dof_pos] = arr_i
    infile.close()
    return uh

infile = Path(solname_h5)

uh = read_mesh_data(infile, mesh, "Function/Velocity/0")
print(uh)
6

Just remove this assert, and replace,

with

  mesh.geometry.cmap.degree,
        mesh.geometry.cmap.variant,

(also in the following part of the code, where it uses mesh.geometry.cmaps[0]. This is an API change since I made this workaround.