Importing domain/boundary markers in fenicsx

mesh
1

Dear all,

I have created a mesh with different physical groups in gmsh, converted in to xdmf and loaded it in dolfinx (see code below).

However, I wasn’t able to get the physical groups in dolfinx. I’ve found several tutorials on this topic, but all seem outdated?

Version 1. As suggested here (resp. in this code, line 143), I tried using read_information_int, but this method seems to have been removed.

Version 2. In this tutorial, it is suggested to use the find-method on MeshTags, which also seem to have been removed.

Version 3. In this post using boundary_markers = dolfinx.mesh.MeshFunction("size_t", mesh, mesh.topology().dim() - 1) is suggested, but MeshFunctions seem to also have been removed. Which, simultaneously, makes msh2xdmf obsolete, I think?

So, what to do? How can I access the physical groups in the current dolfinx release (0.4)?
Thank you in advance.

Here’s the code:

# load_mesh.py
import meshio

import dolfinx
from dolfinx.io import XDMFFile
from mpi4py import MPI

with XDMFFile(MPI.COMM_WORLD, "mesh_vol.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")
    ct = xdmf.read_meshtags(mesh, name="Grid")
    ## Version 1.
    # tag_info = xdmf.read_information_int()

mesh.topology.create_connectivity(mesh.topology.dim-1, 0)

with XDMFFile(MPI.COMM_WORLD, "mesh_surf.xdmf", "r") as xdmf:
    ft = xdmf.read_meshtags(mesh, name="Grid")

## Version 2.
# ft.find

## Version 3
# boundary_markers = dolfinx.mesh.cpp.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)

# convert.py
import meshio

from dolfinx.io import XDMFFile
from mpi4py import MPI

### Convert from msh to XDMF

def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh

mesh = meshio.read("mesh.msh")

mesh_surf = create_mesh(mesh, "triangle")
meshio.write("mesh_surf.xdmf", mesh_surf)

mesh_vol = create_mesh(mesh, "tetra")
meshio.write("mesh_vol.xdmf", mesh_vol)

//+ mesh.geo
SetFactory("OpenCASCADE");
Box(1) = {0, 0, 0, 1, 1, 1};
//+
Physical Surface("s1") = {6};
//+
Physical Surface("s2") = {1};
//+
Physical Surface("s3") = {5};
//+
Physical Surface("s4") = {2};
//+
Physical Surface("s5") = {4};
//+
Physical Surface("s6") = {3};
//+
Physical Volume("dom") = {1};
2

Since you are working with v0.4.0, while the tutorial is using v0.4.2.0, I would suggest you have a look at: https://github.com/jorgensd/dolfinx-tutorial/blob/v0.4.0/chapter3/subdomains.ipynb
In particular: bottom_cells = ct.indices[ct.values==bottom_marker]

1 Like
3

Thanks a bunch, mate. For you answer as well as all the tutorials.

I shouldn’t have dropped the minor version, apparently I’m using 0.4.1 though, at least that’s what condaforge is saying. Your solution still works.

If I would also save the labels into the xdmf-file, would there be a way to read them into fenicsx and access them directly via strings? Or would I still need to use the integer-indices? That is, there is no equivalent to Version 1 anymore, right?

4

It has been a very long time (if ever) that there has been support for labelling with strings. The reason for this is that a string can have arbitrary length, and therefore be incredibly inefficient to save/load.

The solution in version one was removed a long time ago, as MeshFunctions are inefficient for storing information about any sub-entity (vertices, edges, facets), as one almost never need to tag all of them.

5

Sounds reasonable. Thanks again.