Extend mesh.gmsh with ability to convert geo files

This patch generalizes the msh.gmsh function to include .geo input, in which
case the gmsh application is called to convert it to .msh on the fly.

.github/workflows/test.yaml

@@ -90,6 +90,10 @@ jobs:
       - name: Install Graphviz
         if: ${{ matrix.os == 'ubuntu-latest' }}
         run: sudo apt install -y graphviz
+      - name: Install Gmsh
+        # install gmsh via pip on windows and macos and via apt on linux, as
+        # the latter version is dynamically linked and requires libgl etc.
+        run: ${{ matrix.os == 'ubuntu-latest' && 'sudo apt install gmsh' || 'python -um pip install gmsh' }}
       - name: Install Nutils and dependencies
         id: install
         env:

nutils/mesh.py

@@ -12,6 +12,7 @@
 from .topology import Topology
 from math import comb
 from typing import Optional, Sequence, Tuple, Union
+from pathlib import Path
 import numpy
 import os
 import itertools
@@ -20,6 +21,9 @@
 import treelog as log
 import io
 import contextlib
+import tempfile
+import subprocess
+import shutil
 _ = numpy.newaxis
 
 # MESH GENERATORS
@@ -316,6 +320,11 @@ def parsegmsh(mshdata):
 
     msh = gmsh.main.read_buffer(mshdata)
 
+    return _simplex_args_from_meshio(msh)
+
+
+def _simplex_args_from_meshio(msh):
+
     if not msh.cell_sets:
         # Old versions of the gmsh file format repeat elements that have multiple
         # tags. To support this we edit the meshio data to bring it in the same
@@ -448,17 +457,30 @@ def parsegmsh(mshdata):
 
 
 @log.withcontext
-def gmsh(fname, *, space='X'):
+def gmsh(fname, *, dimension=None, order=1, numbers={}, space='X'):
     """Gmsh parser
 
-    Parser for Gmsh files in `.msh` format. Only files with physical groups are
-    supported. See the `Gmsh manual
+    Parser for Gmsh files in `.msh` or `.geo` format. Requires the meshio
+    module to be installed. Conversion of .geo files additionally requires the
+    gmsh application to be installed in the execution path.
+
+    Only files with physical groups are supported. See the `Gmsh manual
     <http://geuz.org/gmsh/doc/texinfo/gmsh.html>`_ for details.
 
     Parameters
     ----------
-    fname : :class:`str` or :class:`io.BufferedIOBase`
-        Path to mesh file or mesh file object.
+    fname : :class:`str` or :class:`pathlib.Path` or :class:`io.BufferedIOBase`
+        Path to .geo or .msh file, or a file object for .msh data.
+    dimension : :class:`int` (optional)
+        Spatial dimension. Specifying this signals to the gmsh function that
+        the input is a .geo file; otherwise the input is assumed to be in .msh
+        format. Valid values are 1, 2 and 3.
+    numbers : :class:`dict` (optional)
+        Only valid for .geo files: number definitions that are added to the
+        beginning of a .geo file via the -setnumber argument.
+    space : :class:`str` (optional)
+        Name of the Nutils topology to distinguish it from other topological
+        directions.
 
     Returns
     -------
@@ -468,8 +490,50 @@ def gmsh(fname, *, space='X'):
         Isoparametric map.
     """
 
-    with util.binaryfile(fname) as f:
-        return simplex(**parsegmsh(f), space=space)
+    try:
+        import meshio
+    except ImportError as e:
+        raise Exception('function.gmsh requires the meshio module to be installed') from e
+
+    if dimension is None:
+        cached_read = cache.function(meshio.gmsh.main.read_buffer)
+        with util.binaryfile(fname) as f:
+            mesh = cached_read(f)
+        return simplex(**_simplex_args_from_meshio(mesh), space=space)
+
+    gmsh = shutil.which('gmsh')
+    if not gmsh:
+        raise RuntimeError('gmsh application does not appear to be installed')
+
+    if not (isinstance(fname, str) and fname.endswith('.geo') or isinstance(fname, Path) and fname.suffix == '.geo'):
+        raise ValueError(f'fname parameter should be a file name with extension .geo, got {fname!r}')
+
+    if not isinstance(dimension, int) or not 1 <= dimension <= 3:
+        raise ValueError(f'dimension parameter should be 1, 2 or 3, got {dimension!r}')
+
+    if not isinstance(order, int) or not order >= 1:
+        raise ValueError(f'order be a strictly positive integer, got {order!r}')
+
+    fid, mshpath = tempfile.mkstemp(suffix='.msh')
+    try:
+        os.close(fid) # release file for writing by gmsh (windows)
+        args = [gmsh, fname, '-o', mshpath, '-order', str(order), f'-{dimension}', '-bin']
+        for name, number in numbers.items():
+            args.extend(['-setnumber', name, str(number)])
+        status = subprocess.run(args, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True)
+        for line in status.stdout.splitlines():
+            try:
+                level, text = line.split(': ', 1)
+                getattr(log, level.strip().lower())(text)
+            except:
+                log.info(line)
+        if status.returncode != 0:
+            raise RuntimeError(f'gmsh failed with error code {status.returncode}')
+        mesh = meshio.read(mshpath)
+    finally:
+        os.unlink(mshpath)
+
+    return simplex(**_simplex_args_from_meshio(mesh), space=space)
 
 
 def simplex(nodes, cnodes, coords, tags, btags, ptags, *, space='X'):

nutils/testing.py

@@ -247,6 +247,15 @@ def assertAlmostEqual64(self, actual, desired, *, atol=2e-15, rtol=2e-3, dtype='
         status.extend(s[i:i+80] for i in range(0, len(s), 80))
         self.fail('\n'.join(status))
 
+    @contextlib.contextmanager
+    def skipIfRaises(self, E=Exception, message=None):
+        try:
+            yield
+        except E as e:
+            if message is None or str(e) == message:
+                self.skipTest(f'caught permissible error: {e}')
+            raise
+
 
 ContextTestCase = TestCase
 

tests/test_mesh.py

@@ -9,8 +9,12 @@ class gmsh(TestCase):
 
     def setUp(self):
         super().setUp()
-        path = pathlib.Path(__file__).parent/'test_mesh'/'mesh{0.ndims}d_p{0.degree}_v{0.version}.msh'.format(self)
-        self.domain, self.geom = mesh.gmsh(path)
+        path = pathlib.Path(__file__).parent/'test_mesh'
+        if self.version == 'geo':
+            with self.skipIfRaises(RuntimeError, 'gmsh application does not appear to be installed'):
+                self.domain, self.geom = mesh.gmsh(path/f'mesh{self.ndims}d.geo', order=self.degree, dimension=self.ndims)
+        else:
+            self.domain, self.geom = mesh.gmsh(path/f'mesh{self.ndims}d_p{self.degree}_v{self.version}.msh')
 
     @requires('meshio')
     def test_volume(self):
@@ -83,7 +87,7 @@ def test_refinesubset(self):
 
 
 for ndims in 2, 3:
-    for version in 2, 4:
+    for version in 2, 4, 'geo':
         for degree in range(1, 5 if ndims == 2 else 3):
             gmsh(ndims=ndims, version=version, degree=degree)
 
@@ -93,8 +97,12 @@ class gmshmanifold(TestCase):
 
     def setUp(self):
         super().setUp()
-        path = pathlib.Path(__file__).parent/'test_mesh'/'mesh3dmani_p{0.degree}_v{0.version}.msh'.format(self)
-        self.domain, self.geom = mesh.gmsh(path)
+        path = pathlib.Path(__file__).parent/'test_mesh'
+        if self.version == 'geo':
+            with self.skipIfRaises(RuntimeError, 'gmsh application does not appear to be installed'):
+                self.domain, self.geom = mesh.gmsh(path/'mesh3dmani.geo', order=self.degree, dimension=3)
+        else:
+            self.domain, self.geom = mesh.gmsh(path/f'mesh3dmani_p{self.degree}_v{self.version}.msh')
 
     @requires('meshio')
     def test_volume(self):
@@ -107,7 +115,7 @@ def test_length(self):
         self.assertAllAlmostEqual(length, 2*numpy.pi, places=1 if self.degree == 1 else 3)
 
 
-for version in 2, 4:
+for version in 2, 4, 'geo':
     for degree in 1, 2:
         gmshmanifold(version=version, degree=degree)