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
So, I thought, let’s do the article on Meshio>=4.0.0 and Fenics and show how interchangeable the gmsh itself, gmsh Python API, and pygmsh are. Well, I seem to not get it, but at least I did not manage to show that.
Using pygmsh with physical labels and Fenics is a bit unclear to me. I got it to work eventually, so here are the code blocks.
Generate the mesh using pygmsh
import pygmsh
# OpenCascade in pygmsh seems not to support extraction of lines from a rectangle (... to use with physical labels).
# So, let's use the geo kernel:
with pygmsh.geo.Geometry() as geom:
r1 = geom.add_rectangle(0., 5e-3, 0., 2.5e-3, z=0.)
geom.add_physical(r1.lines, label="1")
geom.add_physical(r1.surface, label="2")
mesh = geom.generate_mesh(dim=2)
# We'll use gmsh format version 2.2 here, as there's a problem
# with writing nodes in the format version 4.1 here, that I cannot figure out
mesh.write("test.msh", file_format="gmsh22")
So, here’s the first oddity I would not get my head around: There seems to be no easy way to access the boundaries of the rectangle generated with the OpenCASCADE kernel. In gmsh’s API they were there as the first four one-dimensional items (although without the tutorial file there would have been no way I could have guessed that).
The second problem was writing the generated mesh to a gmsh format version 4.1, which resulted in an error message I could not quite track back:
>>> mesh.write("test.msh", file_format="gmsh") # That's gmsh41
---------------------------------------------------------------------------
WriteError Traceback (most recent call last)
in
13 # with writing nodes in the format version 4.1 here, that I cannot
14 # figure out
---> 15 mesh.write("test.msh", file_format="gmsh")
/usr/local/lib/python3.6/dist-packages/meshio/_mesh.py in write(self, path_or_buf, file_format, **kwargs)
158 from ._helpers import write
159
--> 160 write(path_or_buf, self, file_format, **kwargs)
161
162 def get_cells_type(self, cell_type):
/usr/local/lib/python3.6/dist-packages/meshio/_helpers.py in write(filename, mesh, file_format, **kwargs)
144
145 # Write
--> 146 return writer(filename, mesh, **kwargs)
/usr/local/lib/python3.6/dist-packages/meshio/gmsh/main.py in (f, m, **kwargs)
109 {
110 "gmsh22": lambda f, m, **kwargs: write(f, m, "2.2", **kwargs),
--> 111 "gmsh": lambda f, m, **kwargs: write(f, m, "4.1", **kwargs),
112 },
113 )
/usr/local/lib/python3.6/dist-packages/meshio/gmsh/main.py in write(filename, mesh, fmt_version, binary, float_fmt)
100 )
101
--> 102 writer.write(filename, mesh, binary=binary, float_fmt=float_fmt)
103
104
/usr/local/lib/python3.6/dist-packages/meshio/gmsh/_gmsh41.py in write(filename, mesh, float_fmt, binary)
356
357 _write_entities(fh, cells, tag_data, mesh.cell_sets, mesh.point_data, binary)
--> 358 _write_nodes(fh, mesh.points, mesh.cells, mesh.point_data, float_fmt, binary)
359 _write_elements(fh, cells, tag_data, binary)
360 if mesh.gmsh_periodic is not None:
/usr/local/lib/python3.6/dist-packages/meshio/gmsh/_gmsh41.py in _write_nodes(fh, points, cells, point_data, float_fmt, binary)
609 if len(cells) != 1:
610 raise WriteError(
--> 611 "Specify entity information to deal with more than one cell type"
612 )
613
WriteError: Specify entity information to deal with more than one cell type
Preparing mesh and boundary files for Fenics
Falling back to gmsh format version 2.2, I could generate the mesh and boundary files like in the original post:
outfile_mesh = f"{prefix:s}_mesh.xdmf"
outfile_boundary = f"{prefix:s}_boundaries.xdmf"
# read input from infile
inmsh = meshio.read(f"{prefix:s}.msh")
# 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)
Just to figure out that while pygmsh allows you to assign a string label to physical groups it numbers them automatically (apparently starting at 0). This is fine, it just seems to be written nowhere to be found.
Modifying the original code at the definition of the Dirichlet Boundary Condition did the trick (full listing at the end):
# ...
bcs = [
do.DirichletBC(FS, do.Constant((0.0, 0.0)), outerwall, 0), # Choose physical group "index" zero here.
]
# ...
Conclusion
While this arguably still works and does the job it took me a few tries to figure out how. So, I hope this helps others at some point.
Full listing of the Fenics part
import dolfin as do
import numpy as np
import matplotlib.pyplot as plt
# 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)
do.plot(mesh); plt.show()
# Generate Function Space
FE = do.FiniteElement("RTE", mesh.ufl_cell(), 3)
FS = do.FunctionSpace(mesh, FE)
# Use markers for boundary conditions (watch for the change! ;-) )
bcs = [
do.DirichletBC(FS, do.Constant((0.0, 0.0)), outerwall, 0), # Choose physical group "index" zero here.
]
# 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 = 20
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 = -gamma^2, find gamma = sqrt(-r)
gamma = np.sqrt(r)
do.plot(f); plt.title(f"{gamma:f}"); plt.show()
computed_eigenvalues.append(r)
print(np.sort(np.array(computed_eigenvalues)))
