From dbfeb9dd77b02ec8f7d22911d3ffc4a89d352f6b Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Mon, 20 Jan 2025 22:25:52 +0100 Subject: [PATCH] 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 | 2 + nutils/mesh.py | 78 +++++++++++++++++++++++++++++++++---- tests/test_mesh.py | 24 ++++++++---- 3 files changed, 90 insertions(+), 14 deletions(-) diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index 73c9a65fb..338450def 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -90,6 +90,8 @@ jobs: - name: Install Graphviz if: ${{ matrix.os == 'ubuntu-latest' }} run: sudo apt install -y graphviz + - name: Install Gmsh + run: python -um pip install gmsh - name: Install Nutils and dependencies id: install env: diff --git a/nutils/mesh.py b/nutils/mesh.py index 9d62b5b54..5029a064f 100644 --- a/nutils/mesh.py +++ b/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 `_ 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(f'gmsh 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'): diff --git a/tests/test_mesh.py b/tests/test_mesh.py index ab20041ba..12dbd4ea4 100644 --- a/tests/test_mesh.py +++ b/tests/test_mesh.py @@ -1,6 +1,6 @@ from nutils import mesh, function, element, transform, topology from nutils.testing import TestCase, parametrize, requires -import pathlib +import pathlib, shutil import numpy @@ -9,8 +9,13 @@ 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': + if not shutil.which('gmsh'): + self.skipTest('gmsh application is not 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 +88,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 +98,13 @@ 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': + if not shutil.which('gmsh'): + self.skipTest('gmsh application is not 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 +117,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)