{"id":87,"date":"2020-12-24T15:06:27","date_gmt":"2020-12-24T14:06:27","guid":{"rendered":"https:\/\/bowfinger.de\/blog\/?p=87"},"modified":"2020-12-24T15:07:17","modified_gmt":"2020-12-24T14:07:17","slug":"breaking-changes-2-pygmsh-fenics-and-meshio4-0-0","status":"publish","type":"post","link":"https:\/\/bowfinger.de\/blog\/2020\/12\/breaking-changes-2-pygmsh-fenics-and-meshio4-0-0\/","title":{"rendered":"Breaking Changes 2: Pygmsh, Fenics and Meshio>=4.0.0"},"content":{"rendered":"\n
Package versions this was tested with (2020-12-22):\ngmsh 4.7.1\nfenics 2019.1.0:latest\nmeshio 4.3.7\npygmsh 7.1.5<\/pre>\n\n\n\n
So, I thought, let’s do the article on Meshio>=4.0.0 and Fenics<\/a> and show how interchangeable the gmsh itself, gmsh Python API, and pygmsh<\/a> are. Well, I seem to not get it, but at least I did not manage to show that.<\/p>\n\n\n\n
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.<\/p>\n\n\n\n
Generate the mesh using pygmsh<\/h2>\n\n\n\n
import pygmsh\n# OpenCascade in pygmsh seems not to support extraction of lines from a rectangle (... to use with physical labels).\n# So, let's use the geo kernel:\nwith pygmsh.geo.Geometry() as geom:\n    r1 = geom.add_rectangle(0., 5e-3, 0., 2.5e-3, z=0.)\n    geom.add_physical(r1.lines, label=\"1\")\n    geom.add_physical(r1.surface, label=\"2\")\n\nmesh = geom.generate_mesh(dim=2)<\/code>\n# We'll use gmsh format version 2.2 here, as there's a problem\n# with writing nodes in the format version 4.1 here, that I cannot figure out\nmesh.write(\"test.msh\", file_format=\"gmsh22\")<\/code><\/pre>\n\n\n\n
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).<\/p>\n\n\n\n
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:<\/p>\n\n\n\n
>>> mesh.write(\"test.msh\", file_format=\"gmsh\") # That's gmsh41\n\n---------------------------------------------------------------------------\nWriteError                                Traceback (most recent call last)\n  in \n      13 # with writing nodes in the format version 4.1 here, that I cannot\n      14 # figure out\n ---> 15 mesh.write(\"test.msh\", file_format=\"gmsh\")\n \/usr\/local\/lib\/python3.6\/dist-packages\/meshio\/_mesh.py in write(self, path_or_buf, file_format, **kwargs)\n     158         from ._helpers import write\n     159 \n --> 160         write(path_or_buf, self, file_format, **kwargs)\n     161 \n     162     def get_cells_type(self, cell_type):\n \/usr\/local\/lib\/python3.6\/dist-packages\/meshio\/_helpers.py in write(filename, mesh, file_format, **kwargs)\n     144 \n     145     # Write\n --> 146     return writer(filename, mesh, **kwargs)\n \/usr\/local\/lib\/python3.6\/dist-packages\/meshio\/gmsh\/main.py in (f, m, **kwargs)\n     109     {\n     110         \"gmsh22\": lambda f, m, **kwargs: write(f, m, \"2.2\", **kwargs),\n --> 111         \"gmsh\": lambda f, m, **kwargs: write(f, m, \"4.1\", **kwargs),\n     112     },\n     113 )\n \/usr\/local\/lib\/python3.6\/dist-packages\/meshio\/gmsh\/main.py in write(filename, mesh, fmt_version, binary, float_fmt)\n     100             )\n     101 \n --> 102     writer.write(filename, mesh, binary=binary, float_fmt=float_fmt)\n     103 \n     104 \n \/usr\/local\/lib\/python3.6\/dist-packages\/meshio\/gmsh\/_gmsh41.py in write(filename, mesh, float_fmt, binary)\n     356 \n     357         _write_entities(fh, cells, tag_data, mesh.cell_sets, mesh.point_data, binary)\n --> 358         _write_nodes(fh, mesh.points, mesh.cells, mesh.point_data, float_fmt, binary)\n     359         _write_elements(fh, cells, tag_data, binary)\n     360         if mesh.gmsh_periodic is not None:\n \/usr\/local\/lib\/python3.6\/dist-packages\/meshio\/gmsh\/_gmsh41.py in _write_nodes(fh, points, cells, point_data, float_fmt, binary)\n     609         if len(cells) != 1:\n     610             raise WriteError(\n --> 611                 \"Specify entity information to deal with more than one cell type\"\n     612             )\n     613 \n WriteError: Specify entity information to deal with more than one cell type<\/code><\/pre>\n\n\n\n
Preparing mesh and boundary files for Fenics<\/h2>\n\n\n\n
Falling back to gmsh format version 2.2, I could generate the mesh and boundary files like in the original post:<\/p>\n\n\n\n
outfile_mesh = f\"{prefix:s}_mesh.xdmf\"\noutfile_boundary = f\"{prefix:s}_boundaries.xdmf\"\n\n# read input from infile\ninmsh = meshio.read(f\"{prefix:s}.msh\")\n# delete third (obj=2) column (axis=1), this strips the z-component\noutpoints = np.delete(arr=inmsh.points, obj=2, axis=1)\n\n# create (two dimensional!) triangle mesh file\noutmsh = meshio.Mesh(points=outpoints,\n                      cells=[('triangle', inmsh.get_cells_type(\"triangle\"))],\n                      cell_data={'Subdomain': [inmsh.cell_data_dict['gmsh:physical']['triangle']]},\n                      field_data=inmsh.field_data)\n\n# write mesh to file\nmeshio.write(outfile_mesh, outmsh)\n# create (two dimensional!) boundary data file\noutboundary = meshio.Mesh(points=outpoints,\n                           cells=[('line', inmsh.get_cells_type('line') )],\n                           cell_data={'Boundary': [inmsh.cell_data_dict['gmsh:physical']['line']]},\n                           field_data=inmsh.field_data)\n# write boundary data to file\nmeshio.write(filename=outfile_boundary, mesh=outboundary)<\/code><\/pre>\n\n\n\n
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.<\/p>\n\n\n\n
Modifying the original code at the definition of the Dirichlet Boundary Condition did the trick (full listing at the end):<\/p>\n\n\n\n
# ...\n\nbcs = [\n     do.DirichletBC(FS, do.Constant((0.0, 0.0)), outerwall, 0), # Choose physical group \"index\" zero here.\n     ]\n\n# ...<\/code><\/pre>\n\n\n\n
Conclusion<\/h2>\n\n\n\n
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.<\/p>\n\n\n\n
Full listing of the Fenics part<\/h2>\n\n\n\n
import dolfin as do\nimport numpy as np\nimport matplotlib.pyplot as plt\n\n# Import Mesh\nmesh = do.Mesh()\nwith do.XDMFFile(outfile_mesh) as meshfile, \\\n        do.XDMFFile(outfile_boundary) as boundaryfile:\n    meshfile.read(mesh)\n    mvc = do.MeshValueCollection(\"size_t\", mesh, 2)\n    boundaryfile.read(mvc, \"Boundary\")\n    outerwall = do.MeshFunction(\"size_t\", mesh, mvc)\n\ndo.plot(mesh); plt.show()\n\n# Generate Function Space\nFE = do.FiniteElement(\"RTE\", mesh.ufl_cell(), 3)\nFS = do.FunctionSpace(mesh, FE)\n\n# Use markers for boundary conditions (watch for the change! ;-) )\nbcs = [\n    do.DirichletBC(FS, do.Constant((0.0, 0.0)), outerwall, 0), # Choose physical group \"index\" zero here.\n    ]\n\n# Trial and Test functions\nE = do.TrialFunction(FS)\nEE = do.TestFunction(FS)\n\n# Helmholtz EVP (Ae = -k_co^2<\/em>B*e)\na = do.inner(do.curl(E), do.curl(EE))*do.dx\nb = do.inner(E, EE)*do.dx\n\n# For EVP make use of PETSc and SLEPc\ndummy = E[0]*do.dx\nA = do.PETScMatrix()\nB = do.PETScMatrix()\n\n# Assemble System\ndo.assemble_system(a, dummy, bcs, A_tensor=A)\ndo.assemble_system(b, dummy, bcs, A_tensor=B)\n\n# Apply Boundaries\n[bc.apply(B) for bc in bcs]\n[bc.apply(A) for bc in bcs]\n\n# Let SLEPc solve that generalised EVP\nsolver = do.SLEPcEigenSolver(A, B)\nsolver.parameters[\"solver\"] = \"krylov-schur\"\nsolver.parameters[\"tolerance\"] = 1e-16\nsolver.parameters[\"problem_type\"] = \"gen_hermitian\"\nsolver.parameters[\"spectrum\"] = \"target magnitude\"\nsolver.parameters[\"spectral_transform\"] = \"shift-and-invert\"\nsolver.parameters[\"spectral_shift\"] = -(2.np.pi\/10e-3)<\/em>*2\n\nneigs = 20\nsolver.solve(neigs)\nprint(f\"Found {solver.get_number_converged():d} solutions.\")\n\n# Return the computed eigenvalues in a sorted array\ncomputed_eigenvalues = []\n\nfor i in range(min(neigs, solver.get_number_converged())):\n    r, _, fieldRe, fieldIm = solver.get_eigenpair(i) # ignore the imaginary part\n    f = do.Function(FS)\n    f.vector()[:] = fieldRe\n    if np.abs(r) > 1.1:<\/code>\n        # With r = -gamma^2, find gamma = sqrt(-r)<\/code>\n        gamma = np.sqrt(r)<\/code>\n        do.plot(f); plt.title(f\"{gamma:f}\"); plt.show()<\/code>\n    computed_eigenvalues.append(r)<\/code>\n\nprint(np.sort(np.array(computed_eigenvalues)))<\/code><\/pre>\n\n\n\n
<\/p>\n","protected":false},"excerpt":{"rendered":"
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. However, I adapted my code a bit. So here are the code blocks…<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"templates\/template-full-width.php","format":"standard","meta":{"inline_featured_image":false,"footnotes":""},"categories":[5],"tags":[8,7,6,11,12],"class_list":["post-87","post","type-post","status-publish","format-standard","hentry","category-fem","tag-dolfin","tag-fenics","tag-finite-elements","tag-pygmsh","tag-python"],"_links":{"self":[{"href":"https:\/\/bowfinger.de\/blog\/wp-json\/wp\/v2\/posts\/87","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/bowfinger.de\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/bowfinger.de\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/bowfinger.de\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/bowfinger.de\/blog\/wp-json\/wp\/v2\/comments?post=87"}],"version-history":[{"count":7,"href":"https:\/\/bowfinger.de\/blog\/wp-json\/wp\/v2\/posts\/87\/revisions"}],"predecessor-version":[{"id":97,"href":"https:\/\/bowfinger.de\/blog\/wp-json\/wp\/v2\/posts\/87\/revisions\/97"}],"wp:attachment":[{"href":"https:\/\/bowfinger.de\/blog\/wp-json\/wp\/v2\/media?parent=87"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/bowfinger.de\/blog\/wp-json\/wp\/v2\/categories?post=87"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/bowfinger.de\/blog\/wp-json\/wp\/v2\/tags?post=87"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}