How to get facet_region.xml from a .msh file by using meshio

Hello everyone,

I’m now trying to convert my .msh file to .xml files, all works fine except for the facet_region.xml. It seems that it transforms only the lines defined as physical lines in .geo file and eliminate all edges of mesh.

this is my inital .geo file:

d = 0.015;
R = 0.25;
Point(1) = {0., 0., 0., d};
Point(2) = {0.4-R, 0., 0., d};
Point(3) = {0.4, 0., 0., d};
Point(4) = {0.4+R, 0., 0., d};
Point(5) = {1., 0., 0., d};
Point(6) = {1., 0.5, 0., d};
Point(7) = {0.6+R, 0.5, 0., d};
Point(8) = {0.6, 0.5, 0., d};
Point(9) = {0.6-R, 0.5, 0., d};
Point(10) = {0., 0.5, 0., d};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 5};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 9};
Line(9) = {9, 10};
Line(10) = {10, 1};
Circle(11) = {4,3,2};
Circle(12) = {9,8,7};
Line Loop(1) = {4, 5, 6, -12, 9, 10, 1, -11};
Plane Surface(1) = {1};
Line Loop(2) = {11, 2, 3};
Plane Surface(2) = {2};
Line Loop(3) = {12, 7, 8};
Plane Surface(3) = {3};
Physical Line(1) = {10};
Physical Line(2) = {5};
Physical Line(3) = {11, 12};
Physical Surface(1) = {1};
Physical Surface(2) = {2,3};

my purpose is to give labels for the facets of the mesh, either exterior ones (labelled 1 and 2 for the left and right boundary respectively) or labelled 3 for the interface between the two materials, and all the others are labelled by default value 0.

but by runding this code:

def collect_segments(mesh):
    segment_data = {}
    segment_types = {"line", "line3"}  # Assuming 'line' and 'line3' are your segment types

    if 'gmsh:physical' in mesh.cell_data:
        # Iterate over each cell block and the associated physical tags
        for cell_block, physical_tags in zip(mesh.cells, mesh.cell_data['gmsh:physical']):
            if cell_block.type in segment_types:
                segment_data[cell_block.type] = dict()
                for tag, cell in zip(physical_tags,
                    segment_data[cell_block.type][tuple(sorted(cell))] = tag
    return segment_data

i’ve got 174 facets (segments) in total, while the real number should be 9156. (it ignores all lines with label 0 )

To solve this problem, i’ve tried a new approach:

def collect_segments(mesh):
    segments = set()
    for cell_block in mesh.cells:
        for cell in
            # Check all two-node combinations within each cell
            if len(cell) >= 2:
                for i in range(len(cell)):
                    for j in range(i + 1, len(cell)):
                        # Create a sorted tuple to avoid duplicate edges being counted
                        edge = tuple(sorted([cell[i], cell[j]]))
                        # print(edge)
    return segments

Fortunately, it gives out all 9156 elements, but since this approach is not Mesh Dependency, it seems to have a problem with the ordering of the tuple of the segments. and the resulting facet_region.xml file can not be correctly read by fenics.

Is there any efficient approach to solve this problem?

Why do you want to mark all facets in the conversion?

I would simply read in the current facet markers, and then find those not marked by the xml file inside dolfin, instead of adding it to the xml file.

Please also note that the xml format is long deprecated, and one should use the xdmf file format, as for instance explained here:

thanks dokken, yes by using xdmf format i’ve solved my problem.

the reason why i want all facet to be marked comes from the fenics example:

Cohesive zone modelling of debonding and bulk fracture

in which, all facets are marked and stocked in .xml file, and in this example, if we simply read the current facet markers, the code will run with error. So i was just wondering how he did it and making it a try.

Unfortunately, i didn’t find a solution to do so, but what you recommand “find those not marked by the xml file inside dolfin, instead of adding it to the xml file” may a new way to solve the problem.

I’ve solved this problem by using xdmf format

as explained in:

Mesh generation and conversion with GMSH and PYGMSH

the facet.xdmf file can be visualised in paraview and it results in:

which only takes into account the physical lines instead of all edges of mesh.

So i’ve corrected the initial conversion code by:

import meshio
import numpy as np

def extract_edges_with_tags(mesh):
    # Dictionary to store edges with tags
    edge_tags = {}

    # First, handle explicitly defined "line" cells
    if "line" in mesh.cells_dict and "line" in mesh.cell_data_dict["gmsh:physical"]:
        line_cells = mesh.cells_dict["line"]
        line_tags = mesh.cell_data_dict["gmsh:physical"]["line"]
        for edge, tag in zip(line_cells, line_tags):
            # Sort the edge tuple to prevent directional issues
            sorted_edge = tuple(sorted(edge))
            edge_tags[sorted_edge] = tag

    # Now, handle other cell types to extract their edges
    for cell_type, cells in mesh.cells_dict.items():
        if cell_type != "line":
            #print("check 2")
            for cell in cells:
                # This assumes triangular cells; adjust for other types as necessary
                for start, end in zip(cell, np.roll(cell, -1)):
                    sorted_edge = tuple(sorted([start, end]))
                    # Only assign default if not already tagged
                    if sorted_edge not in edge_tags:
                        edge_tags[sorted_edge] = 0  # Default tag

    # Convert dictionary to two lists for meshio
    edges, tags = zip(*edge_tags.items()) if edge_tags else ([], [])
    return np.array(edges), np.array(tags)

# Load the mesh
mesh ="plate_incl.msh")

# Extract all unique edges with physical tags or default
edges, tags = extract_edges_with_tags(mesh)

# Create mesh data structure with these edges
facet_mesh = meshio.Mesh(
    cells=[("line", edges)],
    cell_data={"gmsh:physical": [tags]}

# Write to an XDMF file
meshio.write("facet.xdmf", facet_mesh, file_format="xdmf")

and it results in:

by using this result, the former mentioned fenics example works fine.