Write/read mesh and impose boundary condition

mesh
1

Hello there,

I’m new to dolfinx and stuck in the meshing process. I’m trying to construct a 2D notched square mesh, write it into a xdmf file and then read it. A simplified code is here. I can open the mesh file in paraview but it went wrong when setting the boundary. The error message is

RuntimeError: Cannot use mesh::locate_entities_boundary (boundary) for cells.

The mesh is written by

msh, cell_markers, facet_markers = gmshio.model_to_mesh(model, MPI.COMM_SELF, 0)
msh.name = "notched_sample"
cell_markers.name = f"{msh.name}_cells"
facet_markers.name = f"{msh.name}_facets"

with XDMFFile(msh.comm, f"mesh.xdmf", "w") as file:
    file.write_mesh(msh)
    file.write_meshtags(cell_markers)
    msh.topology.create_connectivity(msh.topology.dim - 1, msh.topology.dim)
    file.write_meshtags(facet_markers)

and read by

infile = XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r")
domain = infile.read_mesh(name="notched_sample")
infile.close()

Could you take a look and help me out?

2

First of all, in your mwe you are creating the mesh in two different ways. In general I would suggest to make sure you keep minimal examples to a minimal, removing all code that isn’t needed.

In your particular case, you are calling:

which should take in the geometrical dimension as input, i.e.
msh, cell_markers, facet_markers = gmshio.model_to_mesh(model, MPI.COMM_SELF, 0, gdim=2)

Finally, you are calling:
bottom_facets = mesh.locate_entities_boundary(domain, domain.geometry.dim-1, bottom)
which would be fine if domain.topology.dim == domain.geometry.dim.
In general I would suggest doing
bottom_facets = mesh.locate_entities_boundary(domain, domain.topology.dim-1, bottom)
which would make your previous code run (even if the geometrical dimension was 3 instead of 2).

3

Thanks very much! It works now. I deleted unnecessary code to keep it minimal.