Breaking Changes: Fenics and Meshio>=4.0.0

Package versions this was tested with (2020-12-22):
gmsh 4.7.1
fenics 2019.1.0:latest
meshio 4.3.7
pygmsh 7.1.5

Check your versions with the code provided here.

Meshio is a great piece of software if you’re into converting meshes, or for that matter, if you’re using Fenics for finite elements calculations. However, the problem I experienced with both packages is their rapid development, consequently breaking changes and sometimes a lack of up-to-date documentation.

So, the Fenics community decided to embrace XDMF as mesh format of choice (good idea) but as usual with new stuff it’s WIP. So, there’s an ongoing thread on the discussion group which is extensive, but frankly hard to follow. After a day of fiddeling, I realised that a lot of example code discussed in the thread broke with Meshio version 4.0.0 and newer. It’s a subtle change in Meshio, which makes it non-obvious.

So, with no futher ado, here are code blocks that

1. Generate a simple mesh using Gmsh
   1. A geo file
   2. A code block using Gmsh’s python API
2. Convert msh-file to two XDMF files (mesh and boundary markers) for Fenics to digest
3. Import of the XDMF files into Fenics and a minimal use-case

Generate a simple mesh using Gmsh

The three methods shown below are equivalent. Which one to choose depends on your workflow. It should be noted that getting the Gmsh Python API is not hard but still a bit tricky depending on context.

The goal is to build a rectangle of 5 mm × 2.5 mm with two physical groups (the boundary and the inner area).

Let’s start with the most static but straight-forward aproach: The geo-file straight from Gmsh’s UI. We’ll use the OpenCASCADE kernel (because, why not?):

SetFactory("OpenCASCADE");
Rectangle(1) = {-1.3, 0.5, 0, 5e-3, 2.5e-3, 0};
Physical Surface(2) = {1};
Physical Curve(1) = {1,2,3,4};

Save this file and generate a .msh file using the following command (we’re forcing a 2D mesh with -2 and gmsh format version 4.1 using -format msh41 here):

gmsh -2 -format msh41 test.geo test.msh

And last but not least using Gmsh’s Python API it looks like this:

import gmsh
import numpy as np

gmsh.initialize()
gmsh.model.add("test")
gmsh.logger.start()

r1 = gmsh.model.occ.addRectangle(0, 0, 0, 5e-3, 2.5e-3)
gmsh.model.occ.synchronize()

pg1 = gmsh.model.addPhysicalGroup(2, [1], tag=2) # inner surface
pg2 = gmsh.model.addPhysicalGroup(1, [1, 2, 3, 4], tag=1) # outer wall

gmsh.model.mesh.generate(dim=2)

gmsh.write("test.msh")

This example is closely based on an example in gmsh’s docs (t16.geo, t16.py). I would not say I’d know what gmsh.initialize, gmsh.model.add, gmsh.logger.start and gmsh.model.occ.synchronize actually do.

Convert msh-file to two XDMF files (mesh and boundary markers) for Fenics to digest

As I understand it, Fenics currently requires separate files for mesh and – let’s say – mesh annotations (like boundary names/numbers etc.). So, what the discussion came down to is converting the mesh as obtained from gmsh into two xdmf-files: test_mesh.xdmf and test_boundaries.xdmf.

The code below is shamelessly stolen mostly from dokken‘s post and Frankenstein-copy-pasted together with other contributions from forum users.

import meshio

# Let's introduce some symbolic names:
infile_mesh = "test.msh"
outfile_mesh = "test_mesh.xdmf"
outfile_boundary = "test_boundaries.xdmf"

# read input from infile
inmsh = meshio.read(infile_mesh)

# delete third (obj=2) column (axis=1), this strips the z-component
# outpoints = np.delete(arr=inmsh.points, obj=2, axis=1), create (two dimensional!) triangle mesh file
outmsh = meshio.Mesh(points=outpoints,
                      cells=[('triangle', inmsh.get_cells_type("triangle"))],
                      cell_data={'Subdomain': [inmsh.cell_data_dict['gmsh:physical']['triangle']]},
                      field_data=inmsh.field_data)

# write mesh to file
meshio.write(outfile_mesh, outmsh)

# create (two dimensional!) boundary data file
outboundary = meshio.Mesh(points=outpoints,
                           cells=[('line', inmsh.get_cells_type('line') )],
                           cell_data={'Boundary': [inmsh.cell_data_dict['gmsh:physical']['line']]},
                           field_data=inmsh.field_data)

# write boundary data to file
meshio.write(filename=outfile_boundary, mesh=outboundary)

So, after running this code snipped, you should see two files: test_mesh.xdmf and test_boundaries.xdmf

Import of the XDMF files into Fenics and a minimal use-case

Well, minimal may be something different. But since I’m using this approach for RF-related eigenvalue problems (EVP), this is as close to minimal as it gets:

import dolfin as do
import numpy as np
import matplotlib.pyplot as plt # for plotting mesh/results

# Import Mesh
mesh = do.Mesh()
with do.XDMFFile(outfile_mesh) as meshfile, \
        do.XDMFFile(outfile_boundary) as boundaryfile:
    meshfile.read(mesh)
    mvc = do.MeshValueCollection("size_t", mesh, 2)
    boundaryfile.read(mvc, "Boundary")
    outerwall = do.MeshFunction("size_t", mesh, mvc)

# Uncomment, if you want to see the mesh
# do.plot(mesh); plt.show()

# Generate Function Space
FE = do.FiniteElement("RTE", mesh.ufl_cell(), 1)
FS = do.FunctionSpace(mesh, FE)

# Use markers for boundary conditions
bcs = [
     do.DirichletBC(FS, do.Constant((0.0, 0.0)), outerwall, 1),
     ]

# Trial and Test functions
E = do.TrialFunction(FS)
EE = do.TestFunction(FS)

# Helmholtz EVP (Ae = -k_co^2B*e)
a = do.inner(do.curl(E), do.curl(EE))*do.dx
b = do.inner(E, EE)*do.dx

# For EVP make use of PETSc and SLEPc
dummy = E[0]*do.dx
A = do.PETScMatrix()
B = do.PETScMatrix()

# Assemble System
do.assemble_system(a, dummy, bcs, A_tensor=A)
do.assemble_system(b, dummy, bcs, A_tensor=B)
# Apply Boundaries
[bc.apply(B) for bc in bcs]
[bc.apply(A) for bc in bcs]

# Let SLEPc solve that generalised EVP
solver = do.SLEPcEigenSolver(A, B)
solver.parameters["solver"] = "krylov-schur"
solver.parameters["tolerance"] = 1e-16
solver.parameters["problem_type"] = "gen_hermitian"
solver.parameters["spectrum"] = "target magnitude"
solver.parameters["spectral_transform"] = "shift-and-invert"
solver.parameters["spectral_shift"] = -(2.np.pi/10e-3)*2

neigs = 10
solver.solve(neigs)

print(f"Found {solver.get_number_converged():d} solutions.")

# Return the computed eigenvalues in a sorted array
computed_eigenvalues = []

for i in range(min(neigs, solver.get_number_converged())):
    r, _, fieldRe, fieldIm = solver.get_eigenpair(i) # ignore the imaginary part
    f = do.Function(FS)
    f.vector()[:] = fieldRe

    if np.abs(r) > 1.1:
    # With r = k_co^2, find k_co = sqrt(r)
    k_co = np.sqrt(r)
    do.plot(f); plt.title(f"{gamma:f}"); plt.show()
    computed_eigenvalues.append(r)

print("All computed eigenvalues:")
print(np.sort(np.array(computed_eigenvalues)))

This gives me a 3 modes for this 5 mm by 2.5 mm waveguide and a bunch of spurious modes. That are filtered out because their k_co^2 is 1 or smaller.

Conclusion

This used to be easier (let’s say back in 2016). I’m fairly sure this is not the most efficient way (both programmatically and numerically) and I’m sure as well that there are better ways to solve this problem. For instance, instead of a generic triangular mesh a more structured triangular mesh or even a quad mesh would produce better solutions. Also, a higher order finite element should produce “smoother” solutions.

I also wanted to add a code block for pygmsh, but it turned out that there are a few catches with physical labels, how they’re numbered and labelled, and how to handle them in Fenics that I could not figure out. So, there’s a separate post on that here.

I am happy and thankful for any feedback!

How to check your package versions

import fenics
import meshio
import gmsh
import pygmsh

print("fenics:", fenics.__version__,
      "meshio:", meshio.__version__,
      "gmsh:", gmsh.__version__,
      "pygmsh:", pygmsh.__version__)