{"id":45,"date":"2020-12-22T17:54:58","date_gmt":"2020-12-22T16:54:58","guid":{"rendered":"https:\/\/bowfinger.de\/blog\/?p=45"},"modified":"2020-12-24T15:07:03","modified_gmt":"2020-12-24T14:07:03","slug":"breaking-changes-fenics-and-meshio4-0-0","status":"publish","type":"post","link":"https:\/\/bowfinger.de\/blog\/2020\/12\/breaking-changes-fenics-and-meshio4-0-0\/","title":{"rendered":"Breaking Changes: 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\n\nCheck your versions with the code provided here<\/a>.<\/pre>\n\n\n\n
Meshio<\/a> is a great piece of software if you’re into converting meshes, or for that matter, if you’re using Fenics<\/a> 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.<\/p>\n\n\n\n
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<\/a> 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<\/a> broke with Meshio version 4.0.0 and newer. It’s a subtle change in Meshio, which makes it non-obvious.<\/p>\n\n\n\n
So, with no futher ado, here are code blocks that<\/p>\n\n\n\n
  1. Generate a simple mesh using Gmsh<\/a>
    1. A geo file<\/li>
    2. A code block using Gmsh’s python API<\/li><\/ol><\/li>
    3. Convert msh<\/code>-file to two XDMF files (mesh and boundary markers) for Fenics to digest<\/li>
    4. Import of the XDMF files into Fenics and a minimal use-case<\/li><\/ol>\n\n\n\n
      Generate a simple mesh using Gmsh<\/h2>\n\n\n\n
      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.<\/p>\n\n\n\n
      The goal is to build a rectangle of 5 mm \u00d7 2.5 mm with two physical groups (the boundary and the inner area).<\/p>\n\n\n\n
      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?):<\/p>\n\n\n\n
      SetFactory(\"OpenCASCADE\");\nRectangle(1) = {-1.3, 0.5, 0, 5e-3, 2.5e-3, 0};\nPhysical Surface(2) = {1};\nPhysical Curve(1) = {1,2,3,4};<\/code><\/pre>\n\n\n\n
      Save this file and generate a .msh<\/code> file using the following command (we’re forcing a 2D mesh with -2<\/code> and gmsh format version 4.1 using -format msh41<\/code> here):<\/p>\n\n\n\n
      gmsh -2 -format msh41 test.geo test.msh<\/code><\/pre>\n\n\n\n
      And last but not least using Gmsh’s Python API it looks like this:<\/p>\n\n\n\n
      import gmsh\nimport numpy as np\n\ngmsh.initialize()\ngmsh.model.add(\"test\")\ngmsh.logger.start()\n\nr1 = gmsh.model.occ.addRectangle(0, 0, 0, 5e-3, 2.5e-3)\ngmsh.model.occ.synchronize()\n\npg1 = gmsh.model.addPhysicalGroup(2, [1], tag=2) # inner surface\npg2 = gmsh.model.addPhysicalGroup(1, [1, 2, 3, 4], tag=1) # outer wall\n\ngmsh.model.mesh.generate(dim=2)\n\ngmsh.write(\"test.msh\")<\/code><\/pre>\n\n\n\n
      This example is closely based on an example in gmsh’s docs (t16.geo<\/a>, t16.py<\/a>). I would not say I’d know what gmsh.initialize<\/code>, gmsh.model.add<\/code>, gmsh.logger.start<\/code> and gmsh.model.occ.synchronize<\/code> actually do.<\/p>\n\n\n\n
      Convert msh<\/code>-file to two XDMF files (mesh and boundary markers) for Fenics to digest<\/h2>\n\n\n\n
      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<\/code>-files: test_mesh.xdmf<\/code> and test_boundaries.xdmf<\/code>.<\/p>\n\n\n\n
      The code below is shamelessly stolen mostly from dokken<\/a>‘s post<\/a> and Frankenstein-copy-pasted together with other contributions from forum users.<\/p>\n\n\n\n
      import meshio\n\n# Let's introduce some symbolic names:\ninfile_mesh = \"test.msh\"\noutfile_mesh = \"test_mesh.xdmf\"\noutfile_boundary = \"test_boundaries.xdmf\"\n\n# read input from infile\ninmsh = meshio.read(infile_mesh)\n\n# delete third (obj=2) column (axis=1), this strips the z-component\n# outpoints = np.delete(arr=inmsh.points, obj=2, axis=1), 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\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\n# write boundary data to file\nmeshio.write(filename=outfile_boundary, mesh=outboundary)<\/code><\/pre>\n\n\n\n
      So, after running this code snipped, you should see two files: test_mesh.xdmf<\/code> and test_boundaries.xdmf<\/code><\/p>\n\n\n\n
      Import of the XDMF files into Fenics and a minimal use-case<\/h2>\n\n\n\n
      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:<\/p>\n\n\n\n
      import dolfin as do\nimport numpy as np\nimport matplotlib.pyplot as plt # for plotting mesh\/results\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\n# Uncomment, if you want to see the mesh\n# do.plot(mesh); plt.show()\n\n# Generate Function Space\nFE = do.FiniteElement(\"RTE\", mesh.ufl_cell(), 1)\nFS = do.FunctionSpace(mesh, FE)\n\n# Use markers for boundary conditions\nbcs = [\n     do.DirichletBC(FS, do.Constant((0.0, 0.0)), outerwall, 1),\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# 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 = 10\nsolver.solve(neigs)\n\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\n    if np.abs(r) > 1.1:<\/code>\n    # With r = k_co^2, find k_co<\/code> = sqrt(r)<\/code>\n    k_co<\/code> = np.sqrt(r)<\/code>\n    do.plot(f); plt.title(f\"{gamma:f}\"); plt.show()<\/code>\n    computed_eigenvalues.append(r)<\/code>\n\nprint(\"All computed eigenvalues:\")\nprint(np.sort(np.array(computed_eigenvalues)))<\/code><\/pre>\n\n\n\n
      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<\/code> is 1 or smaller.<\/p>\n\n\n\n
      Conclusion<\/h2>\n\n\n\n
      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.<\/p>\n\n\n\n
      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<\/a>.<\/p>\n\n\n\n
      I am happy and thankful for any feedback!<\/p>\n\n\n\n
      How to check your package versions<\/h2>\n\n\n\n
      import fenics\nimport meshio\nimport gmsh\nimport pygmsh\n\nprint(\"fenics:\", fenics.__version__,\n      \"meshio:\", meshio.__version__,\n      \"gmsh:\", gmsh.__version__,\n      \"pygmsh:\", pygmsh.__version__)<\/code><\/pre>\n","protected":false},"excerpt":{"rendered":"
      Meshio and Fenics are great pieces of software. However, both packages are in active development and consequently breaking changes and sometimes a lack of up-to-date documentation are usual. Here’s my attempt to give newbies like me some code to start from…<\/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,10,9,11,12],"class_list":["post-45","post","type-post","status-publish","format-standard","hentry","category-fem","tag-dolfin","tag-fenics","tag-finite-elements","tag-gmsh","tag-meshio","tag-pygmsh","tag-python"],"_links":{"self":[{"href":"https:\/\/bowfinger.de\/blog\/wp-json\/wp\/v2\/posts\/45","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=45"}],"version-history":[{"count":39,"href":"https:\/\/bowfinger.de\/blog\/wp-json\/wp\/v2\/posts\/45\/revisions"}],"predecessor-version":[{"id":98,"href":"https:\/\/bowfinger.de\/blog\/wp-json\/wp\/v2\/posts\/45\/revisions\/98"}],"wp:attachment":[{"href":"https:\/\/bowfinger.de\/blog\/wp-json\/wp\/v2\/media?parent=45"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/bowfinger.de\/blog\/wp-json\/wp\/v2\/categories?post=45"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/bowfinger.de\/blog\/wp-json\/wp\/v2\/tags?post=45"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}