Create a function from cell data

dolfinx
1

Hello,

I am trying to port some fenics code to fenicsx. Currently trying to define a function from mesh cell data. Let’s say I have a XDMF mesh on which a cell data field f is defined. I basically want to load this field into a dolfinx Function object.

In fenics, I would load the cell data into a numpy array with meshio, then define a UserExpression with a custom eval_cell function, which would get the right value in the numpy array with the proper cell index. For efficiency, I would implement the UserExpression in a separate C++ file, an then compile it with a dolfin cpp_expression.

However there seems to be no way to do that in dolfinx. Firstly, the use of expressions with custom eval function seems to be deprecated (Equivalent for Expression in dolfinx), and from what I’ve seen from examples, the most common way to define a custom expression is to define a class with an eval function, and then interpolate it on a dolfinx function. However the interpolate function must take coordinates as input, but in my case I have no way to map coordinates to values since I’m working with cell data.

Below a simplified version of the fenics code I’m trying to port to fenicsx

mesh = dolfin.Mesh()

mesh_xdmf = dolfin.XDMFFile(mesh_file)
mesh_xdmf.read(mesh)

mesh_meshio = meshio.read(mesh_file)

f = mesh_meshio.cell_data["f"][0]

cpp_file = "my_expression.cpp"
cpp_code = dolfin.compile_cpp_code(open(cpp_file, 'r').read())
my_expression = dolfin.CompiledExpression(cpp_code.MyExpression, f)
my_function = dolfin.interpolate(my_expression, V) # V is a function space defined earlier

So my problem can be summarized with the two following questions:

  • What’s the proper way in fenicsx to load cell data from a XDMF file into a Function object?
  • What’s the new way of doing this efficiently (what’s the replacement in dolfinx for CompiledExpression)?

Thanks

2

This depends on what your expression is supposed to do. If it is used to define a material parameter, you can have a look at: Defining subdomains for different materials — FEniCSx tutorial
where input cell data from xdmf is converted to a DG 0 function:


Q = dolfinx.FunctionSpace(mesh, ("DG", 0))
kappa = dolfinx.Function(Q)
with kappa.vector.localForm() as loc:
    bottom_cells = ct.indices[ct.values==bottom_marker]
    loc.setValues(bottom_cells, np.full(len(bottom_cells), 1))
    top_cells = ct.indices[ct.values==top_marker]
    loc.setValues(top_cells, np.full(len(top_cells), 0.1))
3

Thanks, you code snippet is what I was looking for. However I’m still running into an issue.

My dolfinx code now looks like something like this:

mesh_XDMF = dolfinx.io.XDMFFile(MPI.COMM_WORLD, mesh_path, "r")
mesh = mesh_XDMF.read_mesh(name='Grid')

DG0 = dolfinx.FunctionSpace(mesh, ("DG", 0))
f = dolfinx.Function(DG0, name="f")

# read f cell data with meshio
mesh_meshio = meshio.read(mesh_path)
f_meshio = mesh_meshio.cell_data['f'][0]

with f.vector.localForm() as loc:
    loc.setValues(np.array(range(loc.local_size), dtype=np.int32), f_meshio)

Note that I am unable to use the read_meshtags as the f function contains float values (it indeed represents a material parameter, but which is not a location marker, think of it as a local geometric feature that has its values between 0 and 1).
Ideally, there would be a read_function which could be used here, that would mirror the write_function function. That’s why I’m still using meshio to retrieve the float data.

However this code does not work. Once the f vector is set, the values are assigned to wrong cells. I guess it has to do with some kind of internal reordering of the cells that is not compatible with the ordering of the meshio mesh.

Do you have any idea of what I can do here?

Thanks

4

Currently you cannot read in float’s / double’s with XDMFFile. You are right that dolfin-x reorders the cells (The reason for this is data locality).

I would suggest one of the following:

  1. Create an issue at Issues · FEniCS/dolfinx · GitHub
  2. Rewrite your mesh data such that it uses integers, and map the integers to material data after reading in the data, as I suggested above.