read_from_msh is broken because gmsh does not return ValueError's

[mildblimp]

Trying to read a .msh file results in an error, because gmsh does not throw a ValueError here. Instead, it prints Error : Gmsh has not been initialized, without throwing an error. This breaks the program, and causes me to have a really hard time importing a mesh made with gmsh into FEniCSx

Additionally, this line should say _gmsh.model.getCurrent().

I think this change in gmsh is also relevant https://gitlab.onelab.info/gmsh/gmsh/-/commit/e79c83ac4289fd45019f137fa2a50fb409d98763.

[jorgensd]

I would like to note that an up to date interface for reading msh directly with FEniCSx is found at: https://jsdokken.com/src/tutorial_gmsh.html (Although it does not recognize that the ValueError is no longer thrown). To circumvent the value error, you could simply call import gmsh; gmsh.initialize(). My source for this code (as github is a mirror) is down for maintenance, I will publish the modified code here:

# Copyright (C) 2021 Jørgen Schartum Dokken
#
# This file is part of DOLFINX_MPC
#
# SPDX-License-Identifier:    MIT

import gmsh as _gmsh
import numpy
from dolfinx.cpp.graph import AdjacencyList_int32
from dolfinx.cpp.io import distribute_entity_data, perm_gmsh
from dolfinx.cpp.mesh import cell_entity_type, create_meshtags, to_type
from dolfinx.io import (extract_gmsh_geometry,
                        extract_gmsh_topology_and_markers, ufl_mesh_from_gmsh)
from dolfinx.mesh import create_mesh
from mpi4py import MPI as _MPI

def read_from_msh(filename: str, cell_data=False, facet_data=False, gdim=None):
    """
    Reads a mesh from a msh-file and returns the dolfin-x mesh.
    Input:
        filename: Name of msh file
        cell_data: Boolean, True of a mesh tag for cell data should be returned
                   (Default: False)
        facet_data: Boolean, True if a mesh tag for facet data should be
                    returned (Default: False)
        gdim: Geometrical dimension of problem (Default: 3)
    """
    if _MPI.COMM_WORLD.rank == 0:
        # Check if gmsh is already initialized
        _gmsh.initialize()
        current_model = _gmsh.model.getCurrent()
        _gmsh.model.add("Mesh from file")
        _gmsh.merge(filename)
        _gmsh.model.setCurrent("Mesh from file")
    output = gmsh_model_to_mesh(_gmsh.model, cell_data=cell_data, facet_data=facet_data, gdim=gdim)
    if _MPI.COMM_WORLD.rank == 0:
        if current_model is None:
            _gmsh.finalize()
        else:
            _gmsh.model.setCurrent(current_model)
    return output

def gmsh_model_to_mesh(model, cell_data=False, facet_data=False, gdim=None):
    """
    Given a GMSH model, create a DOLFIN-X mesh and MeshTags.
        model: The GMSH model
        cell_data: Boolean, True of a mesh tag for cell data should be returned
                   (Default: False)
        facet_data: Boolean, True if a mesh tag for facet data should be
                    returned (Default: False)
        gdim: Geometrical dimension of problem (Default: 3)
    """

    if gdim is None:
        gdim = 3

    if _MPI.COMM_WORLD.rank == 0:
        # Get mesh geometry
        x = extract_gmsh_geometry(model)

        # Get mesh topology for each element
        topologies = extract_gmsh_topology_and_markers(model)

        # Get information about each cell type from the msh files
        num_cell_types = len(topologies.keys())
        cell_information = {}
        cell_dimensions = numpy.zeros(num_cell_types, dtype=numpy.int32)
        for i, element in enumerate(topologies.keys()):
            properties = model.mesh.getElementProperties(element)
            name, dim, order, num_nodes, local_coords, _ = properties
            cell_information[i] = {"id": element, "dim": dim,
                                   "num_nodes": num_nodes}
            cell_dimensions[i] = dim

        # Sort elements by ascending dimension
        perm_sort = numpy.argsort(cell_dimensions)
        # Broadcast cell type data and geometric dimension
        cell_id = cell_information[perm_sort[-1]]["id"]
        tdim = cell_information[perm_sort[-1]]["dim"]
        num_nodes = cell_information[perm_sort[-1]]["num_nodes"]
        cell_id, num_nodes = _MPI.COMM_WORLD.bcast([cell_id, num_nodes], root=0)

        # Check for facet data and broadcast if found
        if facet_data:
            if tdim - 1 in cell_dimensions:
                num_facet_nodes = _MPI.COMM_WORLD.bcast(
                    cell_information[perm_sort[-2]]["num_nodes"], root=0)
                gmsh_facet_id = cell_information[perm_sort[-2]]["id"]
                marked_facets = numpy.asarray(topologies[gmsh_facet_id]["topology"], dtype=numpy.int64)
                facet_values = numpy.asarray(topologies[gmsh_facet_id]["cell_data"], dtype=numpy.int32)
            else:
                raise ValueError("No facet data found in file.")

        cells = numpy.asarray(topologies[cell_id]["topology"], dtype=numpy.int64)
        cell_values = numpy.asarray(topologies[cell_id]["cell_data"], dtype=numpy.int32)

    else:
        cell_id, num_nodes = _MPI.COMM_WORLD.bcast([None, None], root=0)
        cells, x = numpy.empty([0, num_nodes], dtype=numpy.int32), numpy.empty([0, gdim])
        cell_values = numpy.empty((0,), dtype=numpy.int32)
        if facet_data:
            num_facet_nodes = _MPI.COMM_WORLD.bcast(None, root=0)
            marked_facets = numpy.empty((0, num_facet_nodes), dtype=numpy.int32)
            facet_values = numpy.empty((0,), dtype=numpy.int32)

    # Create distributed mesh
    ufl_domain = ufl_mesh_from_gmsh(cell_id, gdim)
    gmsh_cell_perm = perm_gmsh(to_type(str(ufl_domain.ufl_cell())), num_nodes)
    cells = cells[:, gmsh_cell_perm]
    mesh = create_mesh(_MPI.COMM_WORLD, cells, x[:, :gdim], ufl_domain)
    # Create MeshTags for cells
    if cell_data:
        local_entities, local_values = distribute_entity_data(
            mesh, mesh.topology.dim, cells, cell_values)
        mesh.topology.create_connectivity(mesh.topology.dim, 0)
        adj = AdjacencyList_int32(local_entities)
        ct = create_meshtags(mesh, mesh.topology.dim, adj, local_values)
        ct.name = "Cell tags"

    # Create MeshTags for facets
    if facet_data:
        # Permute facets from MSH to Dolfin-X ordering
        # FIXME: This does not work for prism meshes
        facet_type = cell_entity_type(to_type(str(ufl_domain.ufl_cell())),
                                      mesh.topology.dim - 1, 0)
        gmsh_facet_perm = perm_gmsh(facet_type, num_facet_nodes)
        marked_facets = marked_facets[:, gmsh_facet_perm]

        local_entities, local_values = distribute_entity_data(
            mesh, mesh.topology.dim - 1, marked_facets, facet_values)
        mesh.topology.create_connectivity(
            mesh.topology.dim - 1, mesh.topology.dim)
        adj = AdjacencyList_int32(local_entities)
        ft = create_meshtags(mesh, mesh.topology.dim - 1, adj, local_values)
        ft.name = "Facet tags"

    if cell_data and facet_data:
        return mesh, ct, ft
    elif cell_data and not facet_data:
        return mesh, ct
    elif not cell_data and facet_data:
        return mesh, ft
    else:
        return mesh

Note that the recommended way of using GMSH with DOLFINx is to create the mesh, convert it to DOLFINx and then save is as an XDMF file, as it will then be much quicker to read in if you want to run the program multiple times (especially in parallel for large meshes).

[jorgensd]

You could also use meshio, as explained in: https://jsdokken.com/src/pygmsh_tutorial.html#second to convert meshes to a DOLFINx compatible format.

[mildblimp]

You could also use meshio, as explained in: https://jsdokken.com/src/pygmsh_tutorial.html#second to convert meshes to a DOLFINx compatible format.

I've done this initially, after first having tried to use the xdmf format directly (but this gave me an error that cell types were mixed). But in following the tutorial you linked it is unclear to me how you eventually import the mesh into DOLFINx.

[mildblimp]

I will try to use the tutorial from your first comment with the updated script, and see if that resolves it. Thanks for your help so far!

[jorgensd]

Read_from_msh has now been removed from dolfinx_mpc, as it is part of the DOLFINx library now.