How to load and inspect a gmsh mesh (*.msh file)

Question

I want to programmatically inspect *.msh files (get a number of nodes and elements, possibly also edges and faces). How can I do this with either gmsh or pygmsh module? All tutorials I have found so far focus mostly on mesh generation. I can not find any function to read the mesh, not to mention its inspection. Do I have to resort to meshio?

Answer 1

Personally, I tend to export meshes in CalculiX's .inp format by setting Mesh.Format = 39; in my .geo file.

The following code is used to read that .inp file, and extract nodes and (C3D20; 2nd order hexaeder elements) elements.

def read_input(ifn):
    """
    Read an Abaqus INP file, read its sections.
    Return the section headings and the lines.
    """
    with open(ifn) as inf:
        lines = [ln.strip() for ln in inf.readlines()]
    # Remove comments
    lines = [ln for ln in lines if not ln.startswith("**")]
    # Find section headers
    headings = [(ln[1:], n) for n, ln in enumerate(lines) if ln.startswith("*")]
    # Filter the headings so that every heading has a start-of-data and
    # end-of-data index.
    headings.append(("end", -1))
    ln = [h[1] for h in headings]
    headings = [
        (name, start + 1, end) for (name, start), end in zip(headings[:-1], ln[1:])
    ]
    return headings, lines

def retrieve_nodes(headings, lines):
    """
    Extract the nodes out of lines.
    Return a dict of nodes, indexed by the node number.
    A node is a 3-tuple of coordinate strings.
    The node coordinates are *not* converted to floats, so as to not lose precision.

    Arguments:
        headings: list of (name, start, end) tuples.
        lines: list of lines.

    Returns:
        A dict mapping node numbers to (x,y,z)-tuples.
    """
    nodes = {}
    for h in headings:
        if h[0].lower().startswith("node"):
            for ln in lines[h[1] : h[2]]:
                idx, x, y, z = ln.split(",")
                nodes[int(idx)] = (x.strip(), y.strip(), z.strip())
            # Assuming there is only one NODE section.
            break
    return nodes

def retrieve_C3D20(headings, lines):
    """
    Extract C3D20 element sets from the file contents.

    Arguments:
        headings: list of (name, start, end) tuples.
        lines: list of lines.

    Returns:
        1) A dict mapping element numbers to lists of node numbers.
           Note that this is not a separately defined ELSET, but the built-in one.
        2) A dict mapping volume names to lists of element numbers.
        3) A dict mapping node numbers to lists of element numbers.
        4) A dict mapping set name to sets of node numbers.
    """
    all_elements = {}
    element_sets = {}
    all_nreverse = {}
    all_setnodes = {}
    n = 1
    for h in headings:
        name = h[0]
        lname = name.lower()
        if not lname.startswith("element"):
            continue
        if "type=C3D20".lower() not in lname:
            continue
        setname = [s.strip()[6:] for s in name.split(",") if "elset=" in s.lower()]
        if not setname:
            setname = f"unknown{n}"
            n += 1
        else:
            setname = setname[0]
        (
            elements,
            setnodes,
            nreverse,
        ) = read_elements(lines, h[1], h[2], setname)
        element_sets[setname] = elements
        all_elements.update(elements)
        for k, v in nreverse.items():
            if k in all_nreverse:
                all_nreverse[k] += v
            else:
                all_nreverse[k] = v
        all_setnodes[setname] = setnodes
    return all_elements, element_sets, all_nreverse, all_setnodes

This is an excerpt from my gmsh-inp-filter program.

Edit

The file format is defined here

Writing a parser for the $Nodes and $Elements sections doesn't seem too hard.

Answer 2

You can use pygmsh and meshio. Load the file with meshio first and then convert it to a pygmsh Geometry object. If your file has the gmsh:physical and gmsh:geometrical you can define the geometry and then create a mesh.

Here is an example (assumes a simple shape in the mesh file)

import pygmsh
import meshio

mesh = meshio.read("file.msh")
geom = pygmsh.built_in.Geometry()

# Define the nodes
nodes = []
for i in range(mesh.points.shape[0]):
    nodes.append(geom.add_point(mesh.points[i], mesh.point_data["gmsh:physical"][i]))

# Define the elements
element_type = list(mesh.cells.keys())[0] # assume only one element type in the mesh
elements = []
for i in range(mesh.cells[element_type].shape[0]):
    element_points = []
    for j in range(mesh.cells[element_type].shape[1]):
        node_index = mesh.cells[element_type][i][j]
        element_points.append(nodes[node_index])
    elements.append(geom.add_polygon(element_points, mesh.cell_data["gmsh:physical"][i]))

# Define the edges
edges = []
if "line" in mesh.cells:
    for i in range(mesh.cells["line"].shape[0]):
        edge_points = []
        for j in range(mesh.cells["line"].shape[1]):
            node_index = mesh.cells["line"][i][j]
            edge_points.append(nodes[node_index])
        edges.append(geom.add_line(edge_points))

# Define the faces
faces = []
if "quad" in mesh.cells:
    for i in range(mesh.cells["quad"].shape[0]):
        face_points = []
        for j in range(mesh.cells["quad"].shape[1]):
            node_index = mesh.cells["quad"][i][j]
            face_points.append(nodes[node_index])
        faces.append(geom.add_polygon(face_points))

# Generate the mesh
mesh = geom.generate_mesh()

num_nodes = len(nodes)
num_elements = len(elements)
num_edges = len(edges)
num_faces = len(faces)