DG Method From Scratch (in Python) - Gmsh to Numpy
Extracting information from a Gmsh mesh
last update: Feb 6, 2025
Article Motivation
- Create a python class called
Mesh3dthat extracts important information about a Gmsh mesh.- Store mesh connectivity information about cells (tetrahedra) and faces in numpy arrays.
- Visualize the
Mesh3dobject.
Introduction
My goal is to create a python class that lets me take any 3D tetrahedral .msh file from Gmsh, and extract all the relevant information for a discontinuous Galerkin (DG) finite element problem.
Gmsh is an open source 3D finite element mesh generator ran by Christophe Geuzaine out of the Université de Liège in Belgium [1]. It does a lot of great stuff, but what I'm mainly interested in is its ability to generate a 3D unstructured tetrahedral mesh out of common 3D files (like .stl files for example). This will allow me to take MRI data of the human brain and use to produce a 3D volumetric mesh to use in a DG simulation.
How do you turn MRI data into a .msh file?
Typically MRI data comes in the from of a DICOM (.dcm) file. There are a number of ways to turn this into a .stl file (often used in CAD and 3D printing). Gmsh can turn the .stl file into a .msh file.
For now I'm just going to focus on a 3D tetrahedral mesh.
Extracting information from the mesh
Our goal is to create a python class called Mesh that can convert the data in a .msh file from Gmsh into numpy arrays that we can use in a DG wave simulation.
Naming conventions for meshes can vary from person to person. In this cite we'll stick to the following.
A mesh has cells. 3D cells (tetrahedron in our case) have faces, edges, and vertices. All of these are examples of entities.
An element refers to a finite element which can have quadrature points called nodes.
See the posts about meshes and finite elements for more information.
Unfortunately Gmsh uses element and node terminology for meshes, so you might see them interchanged at times.
Really we have three problems to solve. The first is to figure out how to extract information from a .msh file, the second is to decide what information to extract, and the third is to come up with a way to structure the data in numpy arrays so the finite element calculations will be easy to do.
We can solve the first problem using using Gmsh's python API that can return a bunch of useful information. For example:
gmsh.model.mesh.getElementProperties()- tells us the general properties of our meshgmsh.model.mesh.getNodesByElementType()- gets the coordinates of all the verticesgmsh.model.mesh.getElementEdgeNodes()- get the two vertex indices for each edge
The second problem is specific to the discontinuous Galerkin method. We'll need to be able to find each cells neighbors, and which faces are adjacent.
What information do we want to extract?
dimension- dimension of the meshnum_vertices= number of vertices in the meshnum_cells= number of tetrahedral cells in the meshvertexCoordinates= x, y, and z coordinates of all verticesedgeVertices= numpy array of vertex indices for edgescell2vertices= numpy array of vertex indices for cellscell2cells= numpy array of adjacent cell indices for cellscell2faces= numpy array of which faces each cell touchesjacobians= Jacobian matrices from the mappings from reference celldeterminants= determinant of the Jacobian matrices
The third problem can be solved by a clever scheme used in Hesthaven and Warburton's book, Nodal Discontinuous Galerkin Methods [2].
cell2verticeswill have a size of(number of cells, 4)since each tetrahedron has 4 vertices.cell2cellswill have a size of(number of cells, 4)since each cell touches four other cells. This matrix is initiated such that the values in each row are the index of that row. That means is a connection is not found,cells2cellswill report that the cell is connected to itself.cell2faceswill have a size of(number of cells, 4)since each cell touches four faces. This matrix is initialized with every row as[0,1,2,3].
For example, if cell2cells[5] is [29, 12, 5, 40] and cell2faces[5] is [2, 0, 2, 1] that means the following:
- face 0 of cell 5 is connected to face 2 of cell 29
- face 1 of cell 5 is connected to face 0 of cell 12
- face 2 of cell 5 is connected to itself (boundary cell)
- face 3 of cell 5 is connected to face 1 of cell 40
This implies there is a standardized way to keep track of the local face number for each cell. This turns out to be inherent in the way the code is written.
Using Gmsh
Gmsh is a very complex tool, and we are only going to use a small subset of it. Gmsh has a full CAD GUI and can even be used as a post-processor like Paraview.
We only care about a few functions from the python API.
What is a tag in Gmsh?
A tag is just a global identification model for the entity. It is strictly positive and, when combined with the dimension, uniquely defines any entity. From the documentation:
- each point must possess a unique tag;
- each curve must possess a unique tag;
- each surface must possess a unique tag;
- each volume must possess a unique tag.
The Mesh3d python class
class Mesh3d: def __init__(self, msh_file): gmsh.initialize() gmsh.open(msh_file) self.num_vertices = 0 self.num_elements = 0 self.vertexCoordinates = [] self.element2vertices = [] self.element2elements = [] self.element2faces = [] self.faceNormals = [] self.jacobians = {} self.determinants = {} self._extract_mesh_info() self._build_connectivityMatricies() self._compute_jacobians() gmsh.finalize() def _extract_mesh_info(self): # get vertex information ntags, coords, _ = gmsh.model.mesh.getNodes(4) self.num_vertices= len(ntags) self.vertexCoordinates = coords.reshape(-1, 3) # get element information # get all the nodes from tetrahedrons (elementType = 4) nodeTags, _, _ = gmsh.model.mesh.getNodesByElementType(4) self.num_elements = int(len(nodeTags)/4) self.elements2vertices = nodeTags.reshape(-1, 4).astype(int) def _build_connectivityMatricies(self): """tetrahedral face connect algorithm from Toby Isaac""" num_faces = 4 K = self.num_elements EtoV = self.elements2vertices num_vertices = self.num_vertices # create list of all faces faceVertices = np.vstack((EtoV[:, [0, 1, 2]], EtoV[:, [0, 1, 3]], EtoV[:, [1, 2, 3]], EtoV[:, [0, 2, 3]])) # sort each row from low to high for hash algorithm faceVertices = np.sort(faceVertices, axis=1) - 1 # unique hash for each set of three faces by their vertex numbers faceHashes = faceVertices[:, 0] * num_vertices * num_vertices + \ faceVertices[:, 1] * num_vertices + \ faceVertices[:, 2] + 1 # vertex id from 1 - num_vertices * num_elements vertex_ids = np.arange(1, num_faces*K+1) # set up default element to element and element to faces connectivity EtoE = np.tile(np.arange(1, K+1)[:, np.newaxis], (1, num_faces)) EtoF = np.tile(np.arange(1, num_faces+1), (K, 1)) # build a master matrix (mappingTable) that we will solve by # sorting by one column to create the connectivity matricies mappingTable = np.column_stack((faceHashes, vertex_ids, np.ravel(EtoE, order='F'), np.ravel(EtoF, order='F'))) # Now we sort by global face number. sorted_mapTable= np.array(sorted(mappingTable, key=lambda x: (x[0], x[1]))) # find matches in the sorted face list matches = np.where(sorted_mapTable[:-1, 0] == sorted_mapTable[1:, 0])[0] # make links reflexive matchL = np.vstack((sorted_mapTable[matches], sorted_mapTable[matches + 1])) matchR = np.vstack((sorted_mapTable[matches + 1], sorted_mapTable[matches])) # insert matches EtoE_tmp = np.ravel(EtoE, order='F') - 1 EtoF_tmp = np.ravel(EtoF, order='F') - 1 EtoE_tmp[matchL[:, 1] - 1] = (matchR[:, 2] - 1) EtoF_tmp[matchL[:, 1] - 1] = (matchR[:, 3] - 1) EtoE = EtoE_tmp.reshape(EtoE.shape, order='F') EtoF = EtoF_tmp.reshape(EtoF.shape, order='F') self.element2elements = EtoE self.element2faces = EtoF def _compute_jacobians(self): # get local coordinates of the verticies in the # reference tetrahedron name, dim, order, numNodes, localCoords, _ = gmsh.model.mesh.getElementProperties(4) jacobians, determinants, _ = gmsh.model.mesh.getJacobians(4, localCoords) self.jacobians = jacobians.reshape(-1, 3, 3) self.determinants = determinants pass
Plotting in Matplotlib
Now that we have our Mesh3d class, lets test it out. In the Mesh3d class we made a bunch of connectivity matrices. Lets visualize some of these to get a feel for the information we need in a discontinuous Galerkin algorithm.
Lets start by building a function that can take in any cell's index and output a 3D visualization of that cell and its neighbors.
We'll put this function in a file called plot-adjacent-cells.py that contains the following:
import pyvista as pv import numpy as np def plot_mesh(mesh, highlight_cells=None): # Create PyVista UnstructuredGrid cells = np.hstack([np.full((mesh.num_cells, 1), 4), mesh.cell2vertices]).flatten() cell_types = np.full(mesh.num_cells, pv.CellType.TETRA) # Tetrahedral elements grid = pv.UnstructuredGrid(cells, cell_types, mesh.vertexCoordinates) plotter = pv.Plotter() # Plot mesh as wireframe plotter.add_mesh(grid, style='wireframe', color='black') # Highlight specific cells if provided if highlight_cells is not None: # First cell in yellow first_cell = highlight_cells[0] highlight_grid = pv.UnstructuredGrid( np.hstack([[4], mesh.cell2vertices[first_cell]]).flatten(), [pv.CellType.TETRA], mesh.vertexCoordinates ) plotter.add_mesh(highlight_grid, color='#ebcb8b', opacity=1.0) # Remaining cells in blue for cell in highlight_cells[1:]: highlight_grid = pv.UnstructuredGrid( np.hstack([[4], mesh.cell2vertices[cell]]).flatten(), [pv.CellType.TETRA], mesh.vertexCoordinates ) plotter.add_mesh(highlight_grid, color='5e81ac', opacity=0.5) plotter.show() if __name__ == "__main__": import gmsh import sys from mesh import Mesh3d if len(sys.argv) > 1: mesh_file = sys.argv[1] else: gmsh.initialize() gmsh.model.add("simple") gmsh.model.occ.addBox(0, 0, 0, 1, 1, 1) #gmsh.option.setNumber("Mesh.MeshSizeMin", 2.0) gmsh.option.setNumber("Mesh.MeshSizeMax", 0.1) gmsh.model.occ.synchronize() gmsh.model.mesh.generate(3) gmsh.write("default.msh") mesh_file = "default.msh" gmsh.finalize() mesh = Mesh3d(mesh_file) cell_number = 1604 # the cell to visualize adjacent_cells = mesh.cell2cells[cell_number] all_cells = np.insert(adjacent_cells, 0, cell_number) all_cells, indices = np.unique(all_cells, return_index=True) all_cells = all_cells[np.argsort(indices)] # Restore original order plot_mesh(mesh, all_cells)
You can change the cell_number variable at the bottom to visualize other cells.
This produces the following.
Interactive 3D mesh Containing cell number 1604 and the immediate 4 neighbors. Click and drag to move around the mesh.
Conclusion
In this article we've created a python class, Mesh3d that allows us to take any Gmsh .msh file that contains unstructured tetrahedral mesh data and extract the information we need for a discontinuous Galerkin finite element problem. We've stored this data in numpy arrays that we call connectivity matrices. We've also created a 3D plotting function using pyvista that takes in an element and outputs a interactive 3D viewer that plots the cell and its immediate neighbors.