# Mesh generation with Gmsh

Copyright (C) 2020-2023 Garth N. Wells and Jørgen S. Dokken

Download sources

- `Python script`

- `Jupyter notebook`

This demo shows how to create meshes using the Gmsh Python interface.

The Gmsh module is required for this demo.

    from mpi4py import MPI

    import gmsh  # type: ignore

    from dolfinx.io import XDMFFile
    from dolfinx.io import gmsh as gmshio

## Gmsh model builders

The following functions add Gmsh meshes to a ‘model’.

    def gmsh_sphere(model: gmsh.model, name: str) -> gmsh.model:
        """Create a Gmsh model of a sphere and tag sub entitites
        from all co-dimensions (peaks, ridges, facets and cells).

        Args:
            model: Gmsh model to add the mesh to.
            name: Name (identifier) of the mesh to add.

        Returns:
            Gmsh model with a sphere mesh added.

        """
        model.add(name)
        model.setCurrent(name)
        sphere = model.occ.addSphere(0, 0, 0, 1, tag=1)

        # Synchronize OpenCascade representation with gmsh model
        model.occ.synchronize()

        # Add physical tag for sphere
        model.add_physical_group(dim=3, tags=[sphere], tag=1)

        # Embed all sub-entities from the GMSH model into the sphere and tag
        # them
        for dim in [0, 1, 2]:
            entities = model.getEntities(dim)
            entity_ids = [entity[1] for entity in entities]
            model.mesh.embed(dim, entity_ids, 3, sphere)
            model.add_physical_group(dim=dim, tags=entity_ids, tag=dim)

        # Generate the mesh
        model.mesh.generate(dim=3)
        return model

    def gmsh_sphere_minus_box(model: gmsh.model, name: str) -> gmsh.model:
        """Create a Gmsh model of a sphere with a box from the sphere removed.

        Args:
            model: Gmsh model to add the mesh to.
            name: Name (identifier) of the mesh to add.

        Returns:
            Gmsh model with a sphere mesh added.
        """
        model.add(name)
        model.setCurrent(name)

        sphere_dim_tags = model.occ.addSphere(0, 0, 0, 1)
        box_dim_tags = model.occ.addBox(0, 0, 0, 1, 1, 1)
        model_dim_tags = model.occ.cut([(3, sphere_dim_tags)], [(3, box_dim_tags)])
        model.occ.synchronize()

        # Add physical tag 1 for exterior surfaces
        boundary = model.getBoundary(model_dim_tags[0], oriented=False)
        boundary_ids = [b[1] for b in boundary]
        model.addPhysicalGroup(2, boundary_ids, tag=1)
        model.setPhysicalName(2, 1, "Sphere surface")

        # Add physical tag 2 for the volume
        volume_entities = [model[1] for model in model.getEntities(3)]
        model.addPhysicalGroup(3, volume_entities, tag=2)
        model.setPhysicalName(3, 2, "Sphere volume")

        model.mesh.generate(dim=3)
        return model

    def gmsh_ring(model: gmsh.model, name: str) -> gmsh.model:
        """Create a Gmsh model of a ring-type geometry using hexahedral cells.

        Args:
            model: Gmsh model to add the mesh to.
            name: Name (identifier) of the mesh to add.

        Returns:
            Gmsh model with a sphere mesh added.
        """
        model.add(name)
        model.setCurrent(name)

        # Recombine tetrahedra to hexahedra
        gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
        gmsh.option.setNumber("Mesh.RecombineAll", 2)
        gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 1)

        circle = model.occ.addDisk(0, 0, 0, 1, 1)
        circle_inner = model.occ.addDisk(0, 0, 0, 0.5, 0.5)
        cut = model.occ.cut([(2, circle)], [(2, circle_inner)])[0]
        extruded_geometry = model.occ.extrude(cut, 0, 0, 0.5, numElements=[5], recombine=True)
        model.occ.synchronize()

        model.addPhysicalGroup(2, [cut[0][1]], tag=1)
        model.setPhysicalName(2, 1, "2D cylinder")
        boundary_entities = model.getEntities(2)
        other_boundary_entities = []
        for entity in boundary_entities:
            if entity != cut[0][1]:
                other_boundary_entities.append(entity[1])
        model.addPhysicalGroup(2, other_boundary_entities, tag=3)
        model.setPhysicalName(2, 3, "Remaining boundaries")

        model.mesh.generate(3)
        model.mesh.setOrder(2)
        volume_entities = []
        for entity in extruded_geometry:
            if entity[0] == 3:
                volume_entities.append(entity[1])
        model.addPhysicalGroup(3, volume_entities, tag=1)
        model.setPhysicalName(3, 1, "Mesh volume")
        return model

## DOLFINx mesh creation and file output

The following function creates a DOLFINx mesh from a Gmsh model, and cell and facets tags. The mesh and the tags are written to an XDMF file for visualisation, e.g. using ParaView.

    def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mode: str):
        """Create a DOLFINx from a Gmsh model and output to file.

        Args:
            comm: MPI communicator top create the mesh on.
            model: Gmsh model.
            name: Name (identifier) of the mesh to add.
            filename: XDMF filename.
            mode: XDMF file mode. "w" (write) or "a" (append).
        """
        mesh_data = gmshio.model_to_mesh(model, comm, rank=0)
        mesh_data.mesh.name = name
        if mesh_data.cell_tags is not None:
            mesh_data.cell_tags.name = f"{name}_cells"
        if mesh_data.facet_tags is not None:
            mesh_data.facet_tags.name = f"{name}_facets"
        if mesh_data.ridge_tags is not None:
            mesh_data.ridge_tags.name = f"{name}_ridges"
        if mesh_data.peak_tags is not None:
            mesh_data.peak_tags.name = f"{name}_peaks"
        with XDMFFile(mesh_data.mesh.comm, filename, mode) as file:
            mesh_data.mesh.topology.create_connectivity(2, 3)
            mesh_data.mesh.topology.create_connectivity(1, 3)
            mesh_data.mesh.topology.create_connectivity(0, 3)
            file.write_mesh(mesh_data.mesh)
            if mesh_data.cell_tags is not None:
                file.write_meshtags(
                    mesh_data.cell_tags,
                    mesh_data.mesh.geometry,
                    geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{name}']/Geometry",
                )
            if mesh_data.facet_tags is not None:
                file.write_meshtags(
                    mesh_data.facet_tags,
                    mesh_data.mesh.geometry,
                    geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{name}']/Geometry",
                )
            if mesh_data.ridge_tags is not None:
                file.write_meshtags(
                    mesh_data.ridge_tags,
                    mesh_data.mesh.geometry,
                    geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{name}']/Geometry",
                )
            if mesh_data.peak_tags is not None:
                file.write_meshtags(
                    mesh_data.peak_tags,
                    mesh_data.mesh.geometry,
                    geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{name}']/Geometry",
                )

## Generate meshes

Create a Gmsh model and set the verbosity level.

    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0)

    # Create model
    model = gmsh.model()

First, we create a Gmsh model of a sphere using tetrahedral cells (linear geometry), then create independent meshes on each MPI rank and write each mesh to an XDMF file. The MPI rank is appended to the filename since the meshes are not distributed.

    model = gmsh_sphere(model, "Sphere")
    model.setCurrent("Sphere")
    create_mesh(MPI.COMM_SELF, model, "sphere", f"out_gmsh/mesh_rank_{MPI.COMM_WORLD.rank}.xdmf", "w")

Next, we create a Gmsh model of a sphere with a box removed and using tetrahedral cells (linear geometry), then create a distributed mesh. The distributed mesh is written to file. The write option `"w"` is passed to create a new XDMF file.

    model = gmsh_sphere_minus_box(model, "Sphere minus box")
    model.setCurrent("Sphere minus box")
    create_mesh(MPI.COMM_WORLD, model, "ball_d1", "out_gmsh/mesh.xdmf", "w")

For the mesh of the sphere with a box remove, we can increase the degree of the geometry representation to 2 (quadratic geometry representation). The higher-order distributed mesh is appended to the XDMF file.

    model.mesh.generate(3)
    gmsh.option.setNumber("General.Terminal", 1)
    model.mesh.setOrder(2)
    gmsh.option.setNumber("General.Terminal", 0)
    create_mesh(MPI.COMM_WORLD, model, "ball_d2", "out_gmsh/mesh.xdmf", "a")

Finally, we create a distributed mesh using hexahedral cells of geometric degree 2, and append the mesh to the XDMF file.

    model = gmsh_ring(model, "Hexahedral mesh")
    model.setCurrent("Hexahedral mesh")
    create_mesh(MPI.COMM_WORLD, model, "hex_d2", "out_gmsh/mesh.xdmf", "a")

The generated meshes can be visualised using ParaView.