Initializing mesh and function from a file (read sequentially)

I/O
1

I require to read mesh and initial conditions from a file which is generated by another code (outside of FEniCSX). If I write them in a xdmf format, to read mesh is trivial, however, reading function to initialize fields is challenge. If I correctly understood the challenge is attributed to reading in parallel. If so, my question is that how I can read mesh and function in serial, and then FENICSX takes care of domain decomposition and parallel run (any minimal script or pseudo script is helpful)?

2

I would suggest to read:

You have not specified what the input field is tied to though, is it tied to the mesh nodes or the cell?

3

Hi dokken,

Thanks for your suggestion. The data is defined on vertexes of simple CG Hex elements (8-DoFs per element).

The target is initialize a function with a distance function that is computed using VTK library (distance field is computed from surface triangulated geometry; STL file).

I am curios, whether it is possible to do this job without using adios4dolfinx module?

4

You would have to implement the precise code that is provided by adios4dolfinx for communicating the data from an input process to the one owning the degree of freedom.

5

Thanks for your feedback. I just followed your first recipe and use the same scripte you refereed, i.e. this one. I did not modify the code but change input stuff based on my data. However, running the script I face the following error. I did not realize if there are something specific to a element types in this code to be modified. Could you please help me, whether I missed something etc (I put my python script together with data in a zip archive on this link if I can help to realize the problem):

Traceback (most recent call last):
File “/home/rtavakoli/Downloads/tesio/test_init.py”, line 59, in
uh = read_mesh_data(infile, mesh, data)
File “/home/rtavakoli/Downloads/tesio/test_init.py”, line 32, in read_mesh_data
element = basix.ufl.element(
File “/root/miniconda3/envs/dolfinx-checkpoint/lib/python3.10/site-packages/basix/ufl.py”, line 1671, in element
e = _basix.create_element(family, cell, degree, lagrange_variant, dpc_variant, discontinuous)
TypeError: create_element(): incompatible function arguments. The following argument types are supported:
1. (family_name: basix._basixcpp.ElementFamily, cell_name: basix._basixcpp.CellType, degree: int, lagrange_variant: basix._basixcpp.LagrangeVariant = <LagrangeVariant.unset: -1>, dpc_variant: basix._basixcpp.DPCVariant = <DPCVariant.unset: -1>, discontinuous: bool = False, dof_ordering: List[int] = ) → basix._basixcpp.FiniteElement

Invoked with: <ElementFamily.P: 1>, <CellType.hexahedron: 5>, 1, 1, <DPCVariant.unset: -1>, False
Traceback (most recent call last):
File “h5py/_objects.pyx”, line 201, in h5py._objects.ObjectID.dealloc
RuntimeError: Can’t decrement id ref count (can’t close file, there are objects still open)
Exception ignored in: ‘h5py._objects.ObjectID.dealloc’
Traceback (most recent call last):
File “h5py/_objects.pyx”, line 201, in h5py._objects.ObjectID.dealloc
RuntimeError: Can’t decrement id ref count (can’t close file, there are objects still open)

Many thanks for your consideration

6

I would need to know what version of dolfinx you are running. From

it looks like you are running through conda? so you are probably using v0.7.x?

I would appreciate if you could put your code here, instead of in a zip file (for the sake of transparency and accessibility in the future).
What I do assume in the link I referenced is that the data is associated with the vertices of the mesh, and that one has to map these to an appropriate 1st order Lagrange space.

7

Yes I used conda (conda install -n dolfinx-checkpoint -c conda-forge fenics-dolfinx mpich) and installed ADIOS4DOLFINx v0.7.3

And data are attributed to vertices, in fact number of data are identical to number of vertices, is it an issue?

The code is here (identical to your former comment but changing input data; however, it is not possible to attach data to this forum; I put the contents of xdmf file here too)

from pathlib import Path

import adios4dolfinx
import dolfinx
import numpy as np
from mpi4py import MPI

filename = "scaffold.xdmf"
encoding=dolfinx.io.XDMFFile.Encoding.HDF5
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filename, "r", encoding=encoding) as file:
        mesh = file.read_mesh(name="Grid")
        
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", driver="mpio", comm=mesh.comm)
    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)
    assert len(mesh.geometry.cmaps) == 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 doflayout 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("scaffold.h5")


for data in ["data2", "data3"]:
    uh = read_mesh_data(infile, mesh, data)
    with dolfinx.io.XDMFFile(mesh.comm, f"{data}.xdmf", "w") as xdmf:
        xdmf.write_mesh(mesh)
        xdmf.write_function(uh)

scaffold.xdmf file:

<Xdmf Version="3.0"><Domain><Grid Name="Grid"><Geometry GeometryType="XYZ"><DataItem DataType="Float" Dimensions="238328 3" Format="HDF" Precision="8">scaffold.h5:/data0</DataItem></Geometry><Topology TopologyType="Hexahedron" NumberOfElements="226981" NodesPerElement="8"><DataItem DataType="Int" Dimensions="226981 8" Format="HDF" Precision="4">scaffold.h5:/data1</DataItem></Topology><Attribute Name="signed-dist" AttributeType="Scalar" Center="Node"><DataItem DataType="Float" Dimensions="238328 1" Format="HDF" Precision="8">scaffold.h5:/data2</DataItem></Attribute><Attribute Name="level-set" AttributeType="Scalar" Center="Node"><DataItem DataType="Float" Dimensions="238328 1" Format="HDF" Precision="8">scaffold.h5:/data3</DataItem></Attribute></Grid></Domain></Xdmf>

Thanks

8

The same question. I use dolfinx v0.7.3. And I save the mesh and its function (tied to cells) to “.xdmf” and “.h5" file through the method dolfinx.io.XDMFFile.write_mesh and dolfinx.io.XDMFFile.write_function. Now I want to read the mesh data and the function to initialize a dolfinx.fem.Function object. What should I do with it？