# Using the GMSH Python API to generate complex meshes

In this tutorial, we will use the gmsh API to generate complex meshes. We will in this tutorial learn how to make the 3D mesh used in the DFG 3D laminar benchmark. The first part of this tutorial can be completed with the dokken92/pygmsh docker images, as described in the pygmsh tutorial.

For the second and third part of the tutorial, dolfinx is required. You can obtain a jupyter-lab image with dolfinx/lab and a normal docker image with dolfinx/dolfinx.

This tutorial can be downloaded as a Python-file or as a Jupyter notebook.

## 1. How to create a 3D mesh with physical tags with the GMSH python API.

We start by import the gmsh module, and initializing the Python API. NOTE: This has to be done in order to run any function in the GMSH API.

import warnings
warnings.filterwarnings("ignore")
import gmsh
gmsh.initialize()


The next step is to create the rectangular channel of the benchmark. In GMSH, there are two kernels for geometry computations; the built_in kernel ( gmsh.model.geo), and the OpenCascade kernel (gmsh.model.opencascade). In this tutorial, we will use the the occ kernel, as it is better suited. Other demos of the usage of the gmsh python API can be found in their GitLab repository.

gmsh.model.add("DFG 3D")
L, B, H, r = 2.5, 0.41, 0.41, 0.05
channel = gmsh.model.occ.addBox(0, 0, 0, L, B, H)


The next step is to create the cylinder, which is done in the following way

cylinder = gmsh.model.occ.addCylinder(0.5, 0,0.2,0, B, 0, r)


where the first three arguments describes the (x,y,z) coordinate of the center of the first circular face. The next three arguments describes the axis of height of the cylinder, in this case, it is (0,0.41,0). The final parameter is the radius of the cylinder. We have now created two geometrical objects, that each can be meshed with GMSH. However, we are only interested in the fluid volume in the channel, which whould be the channel excluding the sphere. We use the GMSH command BooleanDifference for this

fluid = gmsh.model.occ.cut([(3, channel)], [(3, cylinder)])


where the first argument [(3, channel)]is a list of tuples, where the first argument is the geometrical dimension of the entity (Point=0, Line=1, Surface=2, Volume=3). and channel is a unique integer identifying the channel. Similarly, the second argument is the list of tuples of entities we would like to exclude from the newly created fluid volume. The next step is to tag physical entities, such as the fluid volume, and inlets, outlets, channel walls and obstacle walls. We start by finding the volumes, which after the cut-operation is only the fluid volume. We could have kept the other volumes by supply keyword arguments to the cutoperation. See the GMSH Python API for more information. We also need to syncronize the CAD module before tagging entities.

gmsh.model.occ.synchronize()
volumes = gmsh.model.getEntities(dim=3)
assert(volumes == fluid)
fluid_marker = 11
gmsh.model.setPhysicalName(volumes, fluid_marker, "Fluid volume")


For the surfaces, we start by finding all surfaces, and then compute the geometrical center such that we can indentify which are inlets, outlets, walls and the obstacle. As the walls will consist of multiple surfaces, and the obstacle is circular, we need to find all entites before addin the physical group.

import numpy as np

surfaces = gmsh.model.occ.getEntities(dim=2)
inlet_marker, outlet_marker, wall_marker, obstacle_marker = 1, 3, 5, 7
walls = []
obstacles = []
for surface in surfaces:
com = gmsh.model.occ.getCenterOfMass(surface, surface)
if np.allclose(com, [0, B/2, H/2]):
inlet = surface
gmsh.model.setPhysicalName(surface, inlet_marker, "Fluid inlet")
elif np.allclose(com, [L, B/2, H/2]):
gmsh.model.setPhysicalName(surface, outlet_marker, "Fluid outlet")
elif np.isclose(com, 0) or np.isclose(com, B) or np.isclose(com, H) or np.isclose(com,0):
walls.append(surface)
else:
obstacles.append(surface)
gmsh.model.setPhysicalName(2, wall_marker, "Walls")
gmsh.model.setPhysicalName(2, obstacle_marker, "Obstacle")


The final step is to set mesh resolutions. We will use GMSH Fields to do this. One can alternatively set mesh resolutions at points with the command gmsh.model.occ.mesh.setSize. We start by specifying a distance field from the obstacle surface

gmsh.model.mesh.field.add("Distance", 1)
gmsh.model.mesh.field.setNumbers(1, "FacesList", obstacles)


The next step is to use a threshold function vary the resolution from these surfaces in the following way:

LcMax -                  /--------
/
LcMin -o---------/
|         |       |
Point    DistMin DistMax

resolution = r/10
gmsh.model.mesh.field.setNumber(2, "IField", 1)
gmsh.model.mesh.field.setNumber(2, "LcMin", resolution)
gmsh.model.mesh.field.setNumber(2, "LcMax", 20*resolution)
gmsh.model.mesh.field.setNumber(2, "DistMin", 0.5*r)
gmsh.model.mesh.field.setNumber(2, "DistMax", r)


We add a similar threshold at the inlet to ensure fully resolved inlet flow

gmsh.model.mesh.field.add("Distance", 3)
gmsh.model.mesh.field.setNumbers(3, "FacesList", [inlet])
gmsh.model.mesh.field.setNumber(4, "IField", 3)
gmsh.model.mesh.field.setNumber(4, "LcMin", 5*resolution)
gmsh.model.mesh.field.setNumber(4, "LcMax", 10*resolution)
gmsh.model.mesh.field.setNumber(4, "DistMin", 0.1)
gmsh.model.mesh.field.setNumber(4, "DistMax", 0.5)


We combine these fields by using the minimum field

gmsh.model.mesh.field.add("Min", 5)
gmsh.model.mesh.field.setNumbers(5, "FieldsList", [2, 4])
gmsh.model.mesh.field.setAsBackgroundMesh(5)


Before meshing the model, we need to use the syncronize command

gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(3)


We can write the mesh to msh to be visualized with gmsh with the following command

gmsh.write("mesh3D.msh")


The mesh can be visualized in the GMSH GUI. The figure above visualized the marked facets of the 3D geometry. In the next tutorials, we will learn two diffferent methods of loading this mesh into DOLFINx

## 2. How to load this mesh into DOLFINx without IO

In his tutorial, we will go through the steps of how to load the gmsh geometry above into dolfin-X, without using an intermediate file for storing mesh data.

### Short tutorial

from gmsh_helpers import gmsh_model_to_mesh

mesh, cell_tags, facet_tags = gmsh_model_to_mesh(gmsh.model, cell_data=True, facet_data=True, gdim=3)


This function creates a mesh on processor 0 with GMSH and distributes the mesh data in a dolfinx.Mesh for parallel usage. The flags cell_data and facet_data are booleans that indicates that you would like to extract cell and facet markers from the gmsh model. The last flag gdim indicates the geometrical dimension of your mesh, and should be set to 2 if you want to have a 2D geometry.

### Long tutorial

If you want to learn what the gmsh_model_to_mesh function is actualy doing, the rest of this section will go through it step by step. Note that the long tutorial assumes that you are running in serial, and does therefore not require the MPI-communication found in gmsh_model_to_mesh. We start by using some convenience functions from dolfin-x to extract the mesh geometry (the nodes of the mesh) and the mesh topology (the cell connectivities) for the mesh.

from dolfinx.io import extract_gmsh_geometry, extract_gmsh_topology_and_markers

x = extract_gmsh_geometry(gmsh.model)
topologies = extract_gmsh_topology_and_markers(gmsh.model)


The mesh geometry is a (number of points, 3) array of all the coordinates of the nodes in the mesh. The topologies is a dictionary, where the key is the unique GMSH cell identifier. For each cell type, there is a sub-dictionary containing the mesh topology, an array (number_of_cells, number_of_nodes_per_cell) array containing an integer referring to a row (coordinate) in the mesh geometry, and a 1D array (cell_data) with mesh markers for each cell in the topology.

As an MSH-file can contain meshes for multiple topological dimensions (0=vertices, 1=lines, 2=surfaces, 3=volumes), we have to determine which of the cells has to highest topological dimension. We do this with the following snippet

import numpy

# 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 = gmsh.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)


We extract the topology of the cell with the highest topological dimension from topologies, and create the corresponding ufl.domain.Mesh for the given celltype

from dolfinx.io import ufl_mesh_from_gmsh

cell_id = cell_information[perm_sort[-1]]["id"]
cells = numpy.asarray(topologies[cell_id]["topology"], dtype=numpy.int64)
ufl_domain = ufl_mesh_from_gmsh(cell_id, 3)


As the GMSH model has the cell topology ordered as specified in the MSH format, we have to permute the topology to the FIAT format. The permuation is done using the perm_gmsh function from dolfin-X.

from dolfinx.cpp.io import perm_gmsh
from dolfinx.cpp.mesh import to_type

num_nodes = cell_information[perm_sort[-1]]["num_nodes"]
gmsh_cell_perm = perm_gmsh(to_type(str(ufl_domain.ufl_cell())), num_nodes)
cells = cells[:, gmsh_cell_perm]


The final step is to create the mesh from the topology and geometry

from dolfinx.mesh import create_mesh
from mpi4py import MPI

mesh = create_mesh(MPI.COMM_WORLD, cells, x, ufl_domain)


As the meshes can contain markers for the cells or any sub entity, the next snippets show how to extract this info to GMSH into dolfinx.MeshTags.

from dolfinx.cpp.graph import AdjacencyList_int32

from dolfinx.cpp.io import distribute_entity_data
from dolfinx.cpp.mesh import cell_entity_type, create_meshtags

# Create MeshTags for cell data
cell_values = numpy.asarray(topologies[cell_id]["cell_data"], dtype=numpy.int32)
local_entities, local_values = distribute_entity_data(mesh, mesh.topology.dim, cells, cell_values)
mesh.topology.create_connectivity(mesh.topology.dim, 0)
ct = create_meshtags(mesh, mesh.topology.dim, adj, local_values)
ct.name = "Cell tags"

# Create MeshTags for facets
# 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_id = cell_information[perm_sort[-2]]["id"]
num_facet_nodes = cell_information[perm_sort[-2]]["num_nodes"]
gmsh_facet_perm = perm_gmsh(facet_type, num_facet_nodes)
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)
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)
ft = create_meshtags(mesh, mesh.topology.dim - 1, adj, local_values)
ft.name = "Facet tags"

# Output DOLFINx meshes to file
from dolfinx.io import XDMFFile

with XDMFFile(MPI.COMM_WORLD, "mesh_out.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_meshtags(ft)
xdmf.write_meshtags(ct)


## 3. How to load msh files into dolfin-X

In the previous tutorial, we learnt how to load a gmsh python model into dolfin-X. In this section, we will learn how to load an “msh” file into dolfin-X. We will do this by using the gmsh_model_to_mesh function explained in detail in the previous section. We load the read_from_msh-function from gmsh_helpers.py and use it in the following way

from gmsh_helpers import read_from_msh

mesh, cell_tags, facet_tags = read_from_msh("mesh3D.msh", cell_data=True, facet_data=True, gdim=3)


What this function does, is that it uses the gmsh.merge command to create a gmsh model of the msh file and then in turn calls the gmsh_model_to_mesh function. The read_from_msh function also handles MPI communication and gmsh initialization/finalization.

gmsh.finalize()
gmsh.initialize()