Add gmsh functionality

.github/workflows/test.yaml

@@ -93,6 +93,12 @@ jobs:
         run: |
           sudo apt install -y graphviz
           echo "NUTILS_TESTING_REQUIRES=$NUTILS_TESTING_REQUIRES app:dot" >> $GITHUB_ENV
+      - 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' }}
+          echo "NUTILS_TESTING_REQUIRES=$NUTILS_TESTING_REQUIRES app:gmsh" >> $GITHUB_ENV
       - name: Install Nutils and dependencies
         id: install
         env:

nutils/_graph.py

@@ -6,6 +6,7 @@
 import abc
 import html
 import re
+import shutil
 
 Metadata = TypeVar('Metadata')
 GraphvizColorCallback = Callable[['Node'], Optional[str]]
@@ -61,10 +62,13 @@
         return ''.join(itertools.chain(['digraph {bgcolor="darkgray";'], _generate_graphviz_subgraphs(subgraph_children, nodes, None, id_gen, 0), edges, ['}']))
 
     def export_graphviz(self, *, fill_color: Optional[GraphvizColorCallback] = None, dot_path: str = 'dot', image_type: str = 'svg') -> None:
+        dot = shutil.which(dot_path)
+        if not dot:
+            raise RuntimeError(f'cannot find graphviz application {dot_path!r}')
         src = self.generate_graphviz_source(fill_color=fill_color)
         with treelog.infofile('dot.'+image_type, 'wb') as img:
             src = src.replace(';', ';\n')
-            status = subprocess.run([dot_path, '-Gstart=1', '-T'+image_type], input=src.encode(), stdout=subprocess.PIPE)
+            status = subprocess.run([dot, '-Gstart=1', '-T'+image_type], input=src.encode(), stdout=subprocess.PIPE)
             if status.returncode:
                 for i, line in enumerate(src.split('\n'), 1):
                     print('{:04d}  {}'.format(i, line))

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 @@
 
     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,19 +457,30 @@
 
 
 @log.withcontext
-def gmsh(fname, name='gmsh', *, 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.
-    name : :class:`str` or :any:`None`
-        Name of parsed topology, defaults to 'gmsh'.
+    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
     -------
@@ -470,11 +490,53 @@
         Isoparametric map.
     """
 
-    with util.binaryfile(fname) as f:
-        return simplex(name=name, **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}')
 
-def simplex(nodes, cnodes, coords, tags, btags, ptags, name='simplex', *, space='X'):
+    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'):
     '''Simplex topology.
 
     Parameters
@@ -499,8 +561,6 @@
         preserving order.
     ptags : :class:`dict`
         Dictionary of name->node numbers referencing the ``nodes`` table.
-    name : :class:`str`
-        Name of simplex topology.
 
     Returns
     -------
@@ -510,6 +570,14 @@
         Geometry function.
     '''
 
+    used = numpy.zeros(len(coords), dtype=bool)
+    used[nodes] = True
+    if not used.all():
+        log.warning(f'{len(used) - used.sum()} vertices are not used')
+        for name, points in ptags.items():
+            if not used[points].all():
+                log.warning(f'point {name} is not connected')
+
     nverts = len(coords)
     nelems, ncnodes = cnodes.shape
     ndims = nodes.shape[1] - 1

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':
+            self.require_application('gmsh')
+            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':
+            self.require_application('gmsh')
+            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)