From f752f626a47533bf1151d1b6f66f34cc594c2d26 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 15 Jun 2026 07:46:28 +1000 Subject: [PATCH] fix(meshing): stop dm_plex_gmsh_spacedim leaking across gmsh imports (+ cache guard) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A (root): the dm_plex_gmsh_* options are global, import-time scratch on the PETSc options DB. A manifold-surface generator sets dm_plex_gmsh_spacedim=2 and, if it raises before cleanup (e.g. SegmentedSphericalSurface2D errors in set_periodicity AFTER the import), the option leaks — and the NEXT gmsh import (say a SphericalShell) is read with spacedim=2, producing a 2-D mesh in a 3-D cdim. This surfaced as cross-test pollution: after test_0760_mesh_ot_adapt, a fresh SphericalShell loaded as 2-D ('cannot reshape array of size 856 into shape (3)') and the corrupt 2-D mesh was even written to the shell's cache. Fix centrally in _from_gmsh: clear the whole dm_plex_gmsh_* namespace in a finally after every import (they are meaningful only for the single import that follows), so a value set for one import cannot leak into the next — on success or failure. B (defence-in-depth): in the Mesh constructor, validate the DM coordinate array is divisible by cdim before reshaping; if not, raise a clear, actionable error ('stale/corrupt cached mesh ... delete .meshes/*.msh(.h5)') instead of an opaque numpy reshape failure. Pre-existing core-meshing bug (independent of #202); exposed by the surface tests building a SphericalShell right after a 2-D manifold generator. Full tier-A green; verified: failed manifold construct → next SphericalShell is correct 3-D; no false positives on valid meshes. Underworld development team with AI support from Claude Code --- .../discretisation/discretisation_mesh.py | 86 ++++++++++++++----- 1 file changed, 65 insertions(+), 21 deletions(-) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index bf8c21c1..2ba2c1c2 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -45,6 +45,35 @@ def _temporary_petsc_option(key, value): opts[key] = old_value +# The ``dm_plex_gmsh_*`` options configure how PETSc imports a Gmsh ``.msh`` +# file (coordinate space dimension, label handling, …). They are GLOBAL on the +# PETSc options database and are *import-time scratch* — meaningful only for the +# single import that follows. A mesh constructor that sets one (notably the +# manifold-surface generators set ``dm_plex_gmsh_spacedim = 2``) and then raises +# before cleanup leaks it, and the stale option silently corrupts the NEXT gmsh +# import (e.g. a subsequent SphericalShell is read with spacedim=2 → a 2-D mesh +# in a 3-D cdim → a "cannot reshape" crash, or a corrupt cache). ``_from_gmsh`` +# therefore clears the whole namespace after every import (see below). +_GMSH_IMPORT_OPTION_KEYS = ( + "dm_plex_gmsh_spacedim", + "dm_plex_gmsh_multiple_tags", + "dm_plex_gmsh_use_regions", + "dm_plex_gmsh_mark_vertices", +) + + +def _clear_gmsh_import_options(): + """Remove all ``dm_plex_gmsh_*`` import-scratch options from the global + PETSc options database (idempotent).""" + opts = PETSc.Options() + for key in _GMSH_IMPORT_OPTION_KEYS: + try: + if opts.hasName(key): + opts.delValue(key) + except Exception: + pass + + ## Add the ability to inherit an Enum, so we can add standard boundary ## types to ones that are supplied by the users / the meshing module ## https://stackoverflow.com/questions/46073413/python-enum-combination @@ -103,29 +132,24 @@ def _from_gmsh(filename, comm=None, markVertices=False, useRegions=True, useMult # this process is more efficient done on the root process and then distributed # we do this by saving the mesh as h5 which is more flexible to re-use later - if uw.mpi.rank == 0: - plex_0 = PETSc.DMPlex().createFromFile(filename, interpolate=True, comm=PETSc.COMM_SELF) - - plex_0.setName("uw_mesh") - plex_0.markBoundaryFaces("All_Boundaries", 1001) - - viewer = PETSc.ViewerHDF5().create(filename + ".h5", "w", comm=PETSc.COMM_SELF) - viewer(plex_0) - viewer.destroy() - - # ## Now add some metadata to the mesh (not sure how to do this with the Viewer) - - # import h5py, json - - # f = h5py.File('filename + ".h5",'r+') - - # boundaries_dict = {i.name: i.value for i in cs_mesh.boundaries} - # string_repr = json.dumps(boundaries_dict) + try: + if uw.mpi.rank == 0: + plex_0 = PETSc.DMPlex().createFromFile(filename, interpolate=True, comm=PETSc.COMM_SELF) - # g = f.create_group("metadata") - # g.attrs["boundaries"] = string_repr + plex_0.setName("uw_mesh") + plex_0.markBoundaryFaces("All_Boundaries", 1001) - # f.close() + viewer = PETSc.ViewerHDF5().create(filename + ".h5", "w", comm=PETSc.COMM_SELF) + viewer(plex_0) + viewer.destroy() + finally: + # The gmsh import options are import-time scratch — meaningful only for + # the createFromFile above. Clear the whole namespace so a value set by + # the caller for THIS import (notably ``dm_plex_gmsh_spacedim`` for a + # manifold-surface generator) cannot leak into the next gmsh import and + # silently produce a wrong-dimension mesh (e.g. a later SphericalShell + # read as 2-D). Runs on success or failure. + _clear_gmsh_import_options() # Now we have an h5 file and we can hand this to _from_plexh5 @@ -743,6 +767,26 @@ class replacement_boundaries(Enum): flush=True, ) + # Validate that the DM's coordinate array is consistent with the + # coordinate dimension before reshaping. A mismatch here almost always + # means a STALE or CORRUPT cached mesh file — e.g. a cached .h5 that was + # written as a lower-dimension mesh (the classic symptom of a leaked + # ``dm_plex_gmsh_spacedim`` during generation; see _from_gmsh). Raise a + # clear, actionable error instead of an opaque numpy "cannot reshape" + # failure. + _coord_size = self.dm.getCoordinatesLocal().array.size + if self.cdim and _coord_size % self.cdim != 0: + _src = f" ('{self.filename}')" if getattr(self, "filename", None) else "" + raise RuntimeError( + f"Mesh coordinate array (size {_coord_size}) is not divisible by " + f"the coordinate dimension cdim={self.cdim}, so it cannot be a " + f"valid set of {self.cdim}-D node coordinates. This usually " + f"indicates a stale or corrupt cached mesh file{_src} (e.g. a " + f"cache written at a different space dimension). Delete the " + f"cached '.meshes/*.msh' and '.msh.h5' for this mesh and " + f"regenerate." + ) + # Expose mesh points through special numpy array class with a callback # on all setter operations