Meshio does not read Gmsh physical groups

mesh
8

… and, because this might get important: once created a gmsh file, you have an interest for the import of the gmsh mesh into FEniCS. Below you find two definitions, how I import 2D or 3D meshes into FEniCS (this was derived, e.g., from here and so on).

First you use convert_mesh_xdmf and then load_mesh_to_fenics. It is written such that you can use MPI (see, e.g., here).

Greetings … .

import os
import numpy as np
import meshio
from fenics import *


# _________________________________________________________________ Definitions

# Convert the 'gmsh' specific mesh into a XDMF mesh. This is done for MPI
#
# i             => process number (MPI)
# file_mesh     => gmsh file
# dir_data      => directory of the mesh files
# dimension     => either "2D" or "3D"
#
def convert_mesh_xdmf(i, file_mesh, dir_data, dimension):

    mesh = meshio.read(file_mesh)

    string_nr = "{:04d}".format(i)

    tetra_cells    = None # For 3D only
    tetra_data     = None # For 3D only
    triangle_cells = None # For 2D and 3D
    triangle_data  = None # For 2D and 3D
    line_cells     = None # For 2D only
    line_data      = None # For 2D only

    for key in mesh.cell_data_dict["gmsh:physical"].keys():
        if key == "triangle":
            triangle_data = mesh.cell_data_dict["gmsh:physical"][key]
        elif key == "tetra":
            tetra_data = mesh.cell_data_dict["gmsh:physical"][key]
        elif key == "line":
            line_data = mesh.cell_data_dict["gmsh:physical"][key]

    for cell in mesh.cells:
        if cell.type == "tetra":
            if tetra_cells is None:
                tetra_cells = cell.data
            else:
                tetra_cells = np.vstack([tetra_cells, cell.data])
        elif cell.type == "triangle":
            if triangle_cells is None:
                triangle_cells = cell.data
            else:
                triangle_cells = np.vstack([triangle_cells, cell.data])
        elif cell.type == "line":
            if line_cells is None:
                line_cells = cell.data
            else:
                line_cells = np.vstack([line_cells, cell.data])

    line_mesh  = meshio.Mesh(points=mesh.points,
                             cells={"line": line_cells},
                             cell_data={"name_to_read":[line_data]})
    tetra_mesh = meshio.Mesh(points=mesh.points,
                             cells={"tetra": tetra_cells},
                             cell_data={"name_to_read":[tetra_data]})

    triangle_mesh =meshio.Mesh(points=mesh.points,
                               cells=[("triangle", triangle_cells)],
                               cell_data={"name_to_read":[triangle_data]})

    file_msh_subd = os.path.join(dir_data, "Data_"+string_nr+"_msh_subd.xdmf")
    file_bound    = os.path.join(dir_data, "Data_"+string_nr+"_bound.xdmf")

    if dimension == "2D":
        meshio.write(file_msh_subd, triangle_mesh)
        meshio.write(file_bound,    line_mesh)
    if dimension == "3D":
        meshio.write(file_msh_subd, tetra_mesh)
        meshio.write(file_bound,    triangle_mesh)


# Import the mesh, boundaries and subdomains into FEniCS.
#
# i             => process number (MPI)
# dir_data      => directory of the mesh files
#
def load_mesh_to_fenics(i, dir_data):

    string_nr = "{:04d}".format(i)

    file_msh_subd = os.path.join(dir_data, "Data_"+string_nr+"_msh_subd.xdmf")
    file_bound    = os.path.join(dir_data, "Data_"+string_nr+"_bound.xdmf")

    # Open the latter files and obtain the mesh, subdomains, boundaries,
    # dx and ds.
    mesh = Mesh(MPI.comm_self)
    with XDMFFile(MPI.comm_self, file_msh_subd) as infile:
        infile.read(mesh)

    mvc = MeshValueCollection("size_t", mesh, 2)
    with XDMFFile(MPI.comm_self, file_bound) as infile:
        infile.read(mvc, "name_to_read")
    boundaries = cpp.mesh.MeshFunctionSizet(mesh, mvc)

    mvc2 = MeshValueCollection("size_t", mesh, 3)
    with XDMFFile(MPI.comm_self, file_msh_subd) as infile:
        infile.read(mvc2, "name_to_read")
    subdomains = cpp.mesh.MeshFunctionSizet(mesh, mvc2)

    dx = Measure("dx", domain=mesh, subdomain_data=subdomains)
    ds = Measure("ds", domain=mesh, subdomain_data=boundaries)
    dS = Measure("dS", domain=mesh, subdomain_data=boundaries)

    # We note some statistical information into the parameter class
    n_cells     = mesh.num_cells()
    n_facets    = mesh.num_faces()
    n_vertices  = mesh.num_vertices()
    n_edges     = mesh.num_edges()

    # We print out this info on the terminal for the case that we decide to
    # eventually stop the calculations.
    print("Mesh information")
    print("================")
    print("No. of cells   : ", n_cells)
    print("No. of faces   : ", n_facets)
    print("No. of vertices: ", n_vertices)
    print("No. of edges   : ", n_edges)

    return n_cells, n_facets, n_vertices, n_edges, mesh, subdomains, boundaries, dx, ds, dS