Submesh infrastructure, region labels, and projected normals - #119
Merged
Conversation
Contributor
There was a problem hiding this comment.
Pull request overview
This PR adds core infrastructure to support subdomain (“submesh”) workflows in Underworld3, including extracting region-based submeshes from a parent DMPlex, transferring data between parent/submesh variables, adding region labels to internal-boundary meshes, introducing projected P1 boundary normals, and guarding against mixed-mesh expressions reaching the JIT/solver layer.
Changes:
- Add
Mesh.extract_region()plus parent↔submesh coordination hooks (coordinate sync, re-extraction on adaptation) and parent/submesh data transfer helpers (restrict/prolongate,copy_into/add_into). - Add region-label creation/round-tripping for internal-boundary mesh types (annulus/spherical shell/box internal boundary) and persist regions metadata in HDF5.
- Introduce mixed-mesh expression detection in solvers and
uw.function.evaluate(), and switch Nitsche defaults to use projected normals (Gamma_P1).
Reviewed changes
Copilot reviewed 28 out of 28 changed files in this pull request and generated 4 comments.
Show a summary per file
| File | Description |
|---|---|
| tests/parallel/test_0770_submesh_extract_mpi.py | New MPI regression tests for extract/restrict/prolongate/copy_into and mesh-mixing error. |
| src/underworld3/meshing/spherical.py | Adds regions enum and gmsh physical groups for Inner/Outer volumes on internal-boundary spherical shells; deprecates boundary_normals. |
| src/underworld3/meshing/segmented.py | Deprecates boundary_normals assignment in segmented spherical meshes. |
| src/underworld3/meshing/geographic.py | Deprecates boundary_normals assignment in geographic meshes. |
| src/underworld3/meshing/cartesian.py | Adds region enums/physical groups for internal-boundary box meshes and centroid-based region labels for BoxInternalBoundary. |
| src/underworld3/meshing/annulus.py | Adds regions for internal-boundary annulus and physical groups; deprecates boundary_normals. |
| src/underworld3/function/functions_unit_system.py | Adds mixed-mesh expression checking to uw.function.evaluate(). |
| src/underworld3/function/expressions.py | Adds extract_meshes() utility to detect mesh usage in expressions. |
| src/underworld3/discretisation/enhanced_variables.py | Adds copy_into() / add_into() convenience wrappers for parent↔submesh transfers. |
| src/underworld3/discretisation/discretisation_mesh_variables.py | Clarifies MeshVariable.write() as a low-level collective API. |
| src/underworld3/discretisation/discretisation_mesh.py | Implements extract_region, parent/submesh tracking, adaptation re-extraction, projected normals (Gamma_P1), region persistence, and variable transfer on adapt. |
| src/underworld3/cython/petsc_generic_snes_solvers.pyx | Adds solver-side mixed-mesh checks; switches Nitsche defaults to Gamma_P1; adds optional mask; introduces Region-DS plumbing. |
| src/underworld3/cython/petsc_extras.pxi | Declares PETSc APIs for DMPlexFilter and Region DS functions. |
| src/underworld3/cython/petsc_discretisation.pyx | Adds petsc_dm_filter_by_label() wrapper around DMPlexFilter; improves submesh-from-label wrapper docstring. |
| docs/examples/submesh_investigation/*.py | Adds investigation scripts/demos for submesh, Region DS experiments, and comparisons. |
| docs/developer/design/submesh-solver-architecture.md | New design document describing the intended submesh solver architecture and explored PETSc alternatives. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Comment on lines
+4647
to
+4662
| def set_active_region(self, region_label_name, region_label_value): | ||
| """Configure the solver to assemble only on cells in the given region. | ||
|
|
||
| Cells NOT in the specified region get a trivial DS (no volume | ||
| contributions) and their DOFs should be pinned via Dirichlet BCs. | ||
|
|
||
| Parameters | ||
| ---------- | ||
| region_label_name : str | ||
| DM label name for the INACTIVE region (e.g., "Outer"). | ||
| region_label_value : int | ||
| Label stratum value for the inactive region (e.g., 102). | ||
| """ | ||
| self._inactive_region_label = region_label_name | ||
| self._inactive_region_value = region_label_value | ||
| self.is_setup = False |
There was a problem hiding this comment.
set_active_region(region_label_name, region_label_value) stores region_label_value but _setup_region_ds never uses it (only the label name is consulted). This makes the API misleading and prevents selecting a specific stratum when a label has multiple values. Either remove the unused parameter or use it to select the intended stratum/value when configuring the region DS.
Comment on lines
260
to
268
| class boundary_normals(Enum): | ||
| Lower = new_mesh.CoordinateSystem.unit_e_0 | ||
| Upper = new_mesh.CoordinateSystem.unit_e_0 | ||
| Left = new_mesh.CoordinateSystem.unit_e_1 | ||
| Right = new_mesh.CoordinateSystem.unit_e_1 | ||
| Centre = None | ||
|
|
||
| new_mesh.boundary_normals = boundary_normals | ||
| # boundary_normals deprecated — use mesh.Gamma_P1 for boundary normals | ||
|
|
There was a problem hiding this comment.
The code defines a boundary_normals Enum but no longer assigns it to new_mesh.boundary_normals (only a deprecation comment remains). Since Mesh.init still exposes boundary_normals and existing user scripts may rely on mesh.boundary_normals., this is a breaking change rather than a deprecation. Consider keeping the attribute for at least one release cycle (e.g., assign it and emit a DeprecationWarning on access) while guiding users toward mesh.Gamma_P1.
Comment on lines
+1412
to
+1444
| def _build_dof_map(self, parent_var, sub_var): | ||
| """Build a DOF-level index mapping between parent and submesh variables. | ||
|
|
||
| Uses coordinate matching on DOF coordinates (exact match from | ||
| DMPlexFilter shared nodes). Cached per variable pair. | ||
|
|
||
| Returns (sub_rows, parent_rows) — numpy arrays of matching DOF indices. | ||
| """ | ||
| import numpy as np | ||
| from scipy.spatial import cKDTree | ||
|
|
||
| key = (id(parent_var), id(sub_var)) | ||
| if key in self._dof_maps: | ||
| return self._dof_maps[key] | ||
|
|
||
| tree = cKDTree(sub_var.coords) | ||
| dists, indices = tree.query(parent_var.coords) | ||
| matched = dists < 1.0e-10 | ||
|
|
||
| # indices[matched] maps parent row → sub row | ||
| parent_rows = np.where(matched)[0] | ||
| sub_rows = indices[matched] | ||
|
|
||
| if len(sub_rows) != sub_var.data.shape[0]: | ||
| import warnings | ||
| warnings.warn( | ||
| f"DOF mapping: matched {len(sub_rows)} of " | ||
| f"{sub_var.data.shape[0]} submesh DOFs" | ||
| ) | ||
|
|
||
| result = (sub_rows, parent_rows) | ||
| self._dof_maps[key] = result | ||
| return result |
There was a problem hiding this comment.
restrict/prolongate currently build DOF maps via per-rank KDTree coordinate matching (querying parent_var.coords against sub_var.coords). This is not MPI-safe when the submesh repartitions differently from the parent: a submesh-owned DOF may not exist in the same rank’s parent local(+ghost) vector, so it will never be matched/transferred. Since extract_region already stores subpoint_is from DMPlexFilter, prefer constructing a PETSc IS / VecScatter based on subpoint_is + Section offsets so data transfer is topology-based (and partition-independent) rather than coordinate-based.
Comment on lines
+1653
to
+1659
| for i in range(self.cdim): | ||
| n.data[:, i] = uw.function.evaluate(Gamma[i], n.coords).flatten() | ||
|
|
||
| mag = numpy.sqrt(numpy.sum(n.data ** 2, axis=1)) | ||
| nonzero = mag > 1.0e-30 | ||
| n.data[nonzero] /= mag[nonzero, numpy.newaxis] | ||
|
|
There was a problem hiding this comment.
_update_projected_normals updates MeshVariable data via multiple slice assignments (per component + normalisation). Each assignment triggers the NDArray_With_Callback PETSc sync machinery, which can cause multiple redundant LocalToGlobal/GlobalToLocal scatters and MPI syncs. Consider wrapping the updates in uw.synchronised_array_update() (or computing into a temporary array and calling pack_raw_data_to_petsc once) to avoid unnecessary PETSc communication and improve performance.
Suggested change
| for i in range(self.cdim): | |
| n.data[:, i] = uw.function.evaluate(Gamma[i], n.coords).flatten() | |
| mag = numpy.sqrt(numpy.sum(n.data ** 2, axis=1)) | |
| nonzero = mag > 1.0e-30 | |
| n.data[nonzero] /= mag[nonzero, numpy.newaxis] | |
| projected_normals = numpy.empty_like(n.data) | |
| for i in range(self.cdim): | |
| projected_normals[:, i] = uw.function.evaluate(Gamma[i], n.coords).flatten() | |
| mag = numpy.sqrt(numpy.sum(projected_normals ** 2, axis=1)) | |
| nonzero = mag > 1.0e-30 | |
| projected_normals[nonzero] /= mag[nonzero, numpy.newaxis] | |
| with uw.synchronised_array_update(n): | |
| n.data[...] = projected_normals |
lmoresi
force-pushed
the
feature/region-ds-cell-labels
branch
from
April 22, 2026 06:17
52bf5f0 to
51ad4c5
Compare
Replace single-surface + embed approach with two separate gmsh surfaces sharing the internal boundary curve loop. Each surface gets a Physical Group (Inner=101, Outer=102) that PETSc imports as DM labels. Changes: - annulus.py: Create s_inner and s_outer surfaces, add regions enum, attach mesh.regions attribute - discretisation_mesh.py: Serialize/restore regions enum in HDF5 metadata (same pattern as boundaries) - New test_region_ds_reference.py: Rock-only annulus Stokes reference solution for verifying future Region DS subdomain solving All 21 boundary integral tests pass unchanged. Region cell counts match expected geometry (inner/total ~ 0.375 for default radii). Underworld development team with AI support from Claude Code
Solves Stokes on full AnnulusInternalBoundary mesh with low-viscosity air (eta=1e-3) and compares inner-region norms against rock-only reference. Uses P0-like DG1 MeshVariable for element-wise viscosity. Results confirm the known limitation: penalty free-slip on the internal boundary does not reproduce the rock-only solution because both sides contribute to the stress integral. This motivates the Region DS approach which eliminates air-side assembly entirely. Key findings: - Saddle preconditioner must reflect viscosity field (1/eta), not constant - With proper preconditioner: 1 SNES iteration, seconds to solve - With constant preconditioner: ~800s on same problem - Inner-region velocity ~2x reference (expected with bilateral penalty) Underworld development team with AI support from Claude Code
- viz_region_ds_comparison.py: Interactive jupytext notebook comparing rock-only reference vs air-layer solve with pyvista visualization - Add write_timestep checkpointing to both solve scripts - Add docstring note to MeshVariable.write indicating it is a low-level method; prefer mesh.write_timestep() for normal usage Underworld development team with AI support from Claude Code
- test_region_ds_pinned_air.py: Dirichlet v=0 on all air DOFs via "Outer" region label. Velocity error drops from 157% to 11%. - test_region_ds_nitsche.py: Nitsche BC on internal boundary with viscosity contrast — gives same 157% error as simple penalty, confirming bilateral assembly is the structural problem. The pinned-air result validates the Region DS concept: eliminating air-side contributions (here via Dirichlet pinning) dramatically improves rock-region accuracy. Remaining 11% error is from mesh differences and penalty BC strength, not the approach. Underworld development team with AI support from Claude Code
DMPlexFilter extracts the inner region from AnnulusInternalBoundary as a full-dimension submesh with exact node positions. Solving Stokes on the extracted submesh reproduces the rock-only reference to machine precision (relative error ~1e-12). New infrastructure: - petsc_extras.pxi: Declare DMPlexFilter, DMSetRegionDS, DMGetRegionDS - petsc_discretisation.pyx: petsc_dm_filter_by_label() wrapper - petsc_generic_snes_solvers.pyx: set_active_region() method (Region DS approach — segfaults during assembly, needs further investigation) Test scripts: - test_region_ds_submesh.py: Submesh solve with machine-precision match - test_region_ds_phase3.py: Region DS attempt (incomplete — PETSc segfault) - test_region_ds_pinned_interior.py: Air-interior Dirichlet variant The submesh approach is the viable path for subdomain solving. Underworld development team with AI support from Claude Code
Key findings from comparing submesh vs full-mesh penalty solutions: 1. Air incompressibility constrains radial flow 77x more than penalty alone — the divergence-free air layer acts as a near-rigid boundary 2. Different penalty forms (Gamma vs unit_rvec) are secondary 3. Null space is negligible (~3e-6 relative) With matched penalty (1e6), submesh gives machine-precision match. The air layer provides physics-based free-slip enforcement beyond what the penalty alone achieves. Underworld development team with AI support from Claude Code
With continuous P1 pressure, the 1000x viscosity jump at the internal boundary smears pressure across elements, corrupting the velocity field well into the rock interior. Discontinuous pressure handles each side independently — much better velocity pattern. Updated comparison scripts and notebook to show all three cases: - Blue: rock-only submesh (reference) - Red: air-layer with continuous pressure - Green: air-layer with discontinuous pressure Underworld development team with AI support from Claude Code
Tested DMComposite with rock/air sub-DMs from DMPlexFilter: - Scatter/gather works, interface overlap confirmed (204 shared points) - Composite Vec concatenates sub-DM DOFs — interface DOFs are duplicated - Interface synchronisation still needed after each solve - Conclusion: DMComposite is for combining separate problems, not subdividing one mesh. Direct subpoint IS approach is simpler. Also: updated comparison scripts with normalised Gamma_N and dP1. Underworld development team with AI support from Claude Code
Design doc for multi-domain equation systems in UW3. Documents: - Use cases (air/rock, gravity, surface evolution, multi-physics) - PETSc alternatives investigated (DMComposite, PCFIELDSPLIT, DomainDecomposition) — none directly fit - Chosen approach: DMPlexFilter + subpoint IS + UW3-level restrict/prolongate - Implementation plan: extract_region, restrict/prolongate, solver integration, user API - Open questions: DM lifecycle, point-to-DOF IS, mesh adaptation Also adds bootstrap viscosity test (running in background) and DMComposite probe script. Underworld development team with AI support from Claude Code
Design document for multi-domain equation systems using DMPlexFilter submesh extraction. Covers: - Use cases (air/rock, gravity, surface evolution, multi-physics) - Design principles (one field multiple solvers, mesh lineage, restrict/prolongate, automatic boundary mapping) - PETSc investigation results (DMComposite, PCFIELDSPLIT, DomainDecomposition — none fit exactly, DMPlexFilter + subpoint IS is the right approach) - Implementation plan and open questions - Findings: dP1 required for viscosity contrasts, normalised Gamma_N Also adds bootstrap viscosity test (stepping through eta contrasts using previous solution as initial guess). Underworld development team with AI support from Claude Code
The immediate implementation is just Mesh.extract_region() wrapping DMPlexFilter. Mesh-to-mesh data transfer uses the existing uw.function.evaluate() path — no new infrastructure needed. IS-based restrict/prolongate and solver auto-detection are future optimisations. The two-mesh pattern with explicit evaluate is clearer for users and works with existing code. Underworld development team with AI support from Claude Code
Underworld development team with AI support from Claude Code
The subpoint IS from DMPlexFilter gives exact point correspondence — direct index mapping, no kd-tree, no interpolation. This should be the primary mechanism for parent-submesh data transfer, not a future optimisation. evaluate() remains for unrelated mesh pairs. Underworld development team with AI support from Claude Code
Drop the auto-managed globals concept. Each mesh owns its own MeshVariables. Data moves between meshes via explicit restrict/ prolongate calls. No hidden magic — the user controls data flow. DMComposite can't be used for solving (only block coupling), so the copy must happen. Make it easy and correct. Underworld development team with AI support from Claude Code
Two cases for parent→submesh synchronisation: - Coordinate deformation (ALE): subpoint IS still valid, restrict parent coords to submesh, rebuild geometry. Can auto-detect via mesh version counter. - Topology change (adaptation): subpoint IS invalidated, must re-extract submesh. Parent notifies registered submeshes via weak references (same pattern as _registered_swarms). Underworld development team with AI support from Claude Code
Expressions passed to a solver must only contain symbols from that solver's mesh. The JIT evaluates against one DM's auxiliary vector. Users must restrict cross-mesh data before building expressions. Detect and raise error if meshes are mixed. Underworld development team with AI support from Claude Code
One mesh per expression is a core constraint, not an implementation detail. Moved to design principles section. All MeshVariable symbols in a solver expression must share the solver's mesh. Restrict cross-mesh data first. Detect and error on mesh mismatch. Underworld development team with AI support from Claude Code
New method on the Mesh class wrapping DMPlexFilter. Extracts cells
matching a region label and returns a new Mesh with:
- parent reference to the source mesh
- subpoint_is (PETSc IS mapping submesh points -> parent points)
- boundaries inherited from parent labels
- coordinate system inherited from parent
Usage:
rock_mesh = full_mesh.extract_region("Inner")
rock_mesh.parent # full_mesh
rock_mesh.subpoint_is # IS for restrict/prolongate
Tested: Stokes solve on extracted submesh works end-to-end.
All 21 boundary integral tests pass.
Underworld development team with AI support from Claude Code
Mesh methods for copying MeshVariable data between parent and submesh: - restrict(parent_var, sub_var): parent → submesh DOFs - prolongate(sub_var, parent_var): submesh → parent DOFs - Both support mode="replace" (INSERT) and mode="add" (ADD_VALUES) Uses coordinate matching (cKDTree) on DOF coordinates from DMPlexFilter shared nodes. Mapping is cached per variable pair. Zero error on P1 and P2 variables tested. All 21 boundary integral tests pass. Underworld development team with AI support from Claude Code
The NDArray_With_Callback.copy() preserves the callback subclass, causing spurious callback errors when modifying the copy. Using numpy.array() instead produces a plain ndarray. Data is written back through pack_raw_data_to_petsc() which properly syncs the PETSc Vec. No callback warnings, zero error on P1 and P2 variables. Underworld development team with AI support from Claude Code
User-facing methods on MeshVariable for pushing data between
related meshes:
v_full.copy_into(v_rock) # parent → submesh (restrict)
v_rock.copy_into(v_full) # submesh → parent (prolongate)
v_rock.add_into(v_full) # prolongate with ADD_VALUES
Detects parent/submesh relationship automatically. Raises clear
error if meshes are unrelated. Zero error, no callback warnings.
Underworld development team with AI support from Claude Code
New extract_meshes() function in expressions.py finds all meshes referenced by MeshVariable symbols in a sympy expression. Uses the dynamically-created function class 'meshvar' weakref to trace back to the source mesh. Solver._check_expression_meshes() runs at build time and raises a clear ValueError if any expression contains variables from a foreign mesh, with guidance to use copy_into() for data transfer. Previously this produced a cryptic PrintMethodNotImplementedError deep in the JIT compiler. All 21 boundary integral tests pass. Underworld development team with AI support from Claude Code
Raises ValueError if an expression contains MeshVariable symbols from multiple meshes, before reaching the JIT compiler. Same extract_meshes() utility used by the solver check. All 21 boundary integral tests pass. Underworld development team with AI support from Claude Code
When _deform_mesh is called on a parent mesh, all registered submeshes automatically update their coordinates via a cached vertex index map built at extract_region time. - _build_vertex_map: coordinate matching at extraction (topology-based, survives subsequent deformations) - sync_coordinates_from_parent: copies parent coords at mapped indices, calls _deform_mesh on submesh to rebuild geometry - Parent tracks submeshes via _registered_submeshes (WeakSet) Tested: 1.1x scale deformation propagates with machine precision. All 21 boundary integral tests pass. Underworld development team with AI support from Claude Code
Notebook now loads rock-only, eta=1e-3, and eta=1e-5 (from bootstrap) checkpoints. Pyvista overlay shows all three at matched nodes. Underworld development team with AI support from Claude Code
All 6 tests pass on 2 and 4 MPI ranks: - extract_region produces valid submesh in parallel - restrict (parent→submesh) correct across ranks - prolongate (submesh→parent) correct, air DOFs untouched - copy_into works both directions - Stokes solve on extracted submesh converges - Mixed-mesh expression error detected correctly Underworld development team with AI support from Claude Code
Pass cell region labels to MMG via the rgLabel parameter of adaptMetric. MMG treats interfaces between different cell regions as required boundaries and preserves them during remeshing. - Create _CellRegions_ label from mesh.regions before adaptation - Pass as rgLabel to adaptMetric alongside the existing bdLabel - Reconstruct named region labels (Inner/Outer) on the adapted mesh Tested: internal boundary at r=1.0 preserved to 4e-6 deviation, zero cells crossing the interface, all labels survive, and extract_region works on the adapted mesh. All 21 boundary integral tests pass. Underworld development team with AI support from Claude Code
When mesh.adapt() runs, registered submeshes are automatically re-extracted from the adapted parent via _re_extract_from_parent(): - New DM from DMPlexFilter on adapted parent - Vertex map and DOF maps rebuilt - MeshVariables reinitialised on new DM (reset to zero) - Solvers marked for rebuild - Extraction label stored at extract_region time for re-extraction The submesh Python object is updated in-place — external references remain valid. Variables need reinitialisation after adaptation, same as the parent mesh pattern. Tested: adapt + re-extract + restrict all work with zero error. All 21 boundary integral tests pass. Underworld development team with AI support from Claude Code
Previously, mesh.adapt() reset all variables to zero. Now variables are interpolated from the old mesh to the adapted mesh: - Parent variables: evaluated via uw.function.evaluate at new coords before the DM swap, then restored after reinitialisation - Submesh variables: backed up as (coords, data) arrays before re-extraction, then interpolated via IDW to new submesh coords Transfer error is small (interpolation, not zero) — the data is preserved rather than lost. Users no longer need to manually reinitialise variables after adaptation. All 21 boundary integral tests pass. Underworld development team with AI support from Claude Code
End-to-end test of the multi-mesh data flow: 1. Extract rock submesh from full mesh 2. Set temperature on rock submesh 3. Prolongate T to full mesh (zero in air) 4. Solve Poisson gravity on full mesh 5. Restrict gravity to rock submesh 6. Solve Stokes on rock submesh with buoyancy 7. Prolongate velocity back to full mesh Exercises: extract_region, restrict, prolongate, Poisson on full mesh, Stokes on submesh, checkpointing both meshes. Underworld development team with AI support from Claude Code
- add_nitsche_bc gains optional mask parameter for one-sided application (multiplies all Nitsche terms by an element-wise DG variable) - Note: masked Nitsche/penalty on internal boundaries does NOT work reliably due to PETSc support[0] ordering — the mask evaluates from whichever cell owns the face, not a chosen side - Coupled submesh demonstration (thermal-Stokes with gravity) works correctly using the extract_region / restrict / prolongate path Conclusion: submesh approach is the correct path for internal boundary problems. Nitsche on internal faces cannot be made one-sided without PETSc-level changes to face ownership. Underworld development team with AI support from Claude Code
New mesh property that projects PETSc face normals (Gamma) onto a
continuous P1 field and normalises. Gives smooth, consistently-oriented
unit normals that work for any geometry without analytical formulas.
Key advantages over raw Gamma_N:
- Converges with mesh refinement (Gamma_N penalty diverges in 3D)
- 8x better alignment with true normals at boundary quadrature points
- Consistent orientation (no sign flips on inner boundaries)
- Fast: Nitsche with Gamma_P1 is 2-3x faster than with analytical normals
- Automatic: rebuilt lazily on first access, invalidated on deformation
Recommended for penalty and Nitsche BCs on curved boundaries:
stokes.add_nitsche_bc("Upper", direction=mesh.Gamma_P1,
normal=mesh.Gamma_P1, gamma=10, theta=1)
3D spherical convergence verified at three resolutions.
All 21 boundary integral tests pass.
Underworld development team with AI support from Claude Code
Both SNES_Vector and SNES_Stokes_SaddlePt add_nitsche_bc methods now default to mesh.Gamma_P1 (projected, normalised P1 normals) instead of mesh.Gamma_N (raw normalised PETSc face normals). This gives correct convergence on 3D spherical shells where the old Gamma_N default diverged with mesh refinement. Users can still override with direction= and normal= parameters. All 21 boundary integral tests pass. Underworld development team with AI support from Claude Code
…ry_normals Region labels (Inner/Outer Physical Groups): - SphericalShellInternalBoundary: OCC fragment creates two shell volumes sharing the internal surface. Region labels + extract_region verified. - BoxInternalBoundary: region Physical Groups added to gmsh but useRegions=False (needs careful integration with _dm_unstack_bcs path). TODO: enable useRegions for Box meshes. - AnnulusInternalBoundary: already had regions (unchanged). boundary_normals deprecated: - All mesh factories: replaced `new_mesh.boundary_normals = boundary_normals` with deprecation comment. Use mesh.Gamma_P1 instead. - boundary_normals enum definitions left in place for backward compat but no longer assigned to the mesh. Investigation scripts moved from tests/ to docs/examples/submesh_investigation/ (these are exploration scripts, not pytest tests). All 21 boundary integral tests pass. Underworld development team with AI support from Claude Code
BoxInternalBoundary uses a different label import path (Face Sets +
_dm_unstack_bcs) that conflicts with useRegions=True. Instead of
changing the import mechanism, classify cells by centroid position
after mesh construction: cells below zintCoord = Inner, above = Outer.
Works for both 2D and 3D, simplex and structured meshes.
extract_region("Inner") verified on 2D box.
All 21 boundary integral tests pass.
Underworld development team with AI support from Claude Code
PETSc 3.25 added an MPI_Comm argument to DMPlexFilter. Added UW_DMPlexFilter wrapper with version guard. Also resolved rebase conflict: _check_expression_meshes now runs before the fast-path early return in _build(). Underworld development team with AI support from Claude Code (https://claude.com/claude-code)
The flux property triggers tensor contraction (sympy.tensorcontraction) which fails for scalar diffusion models where the C-tensor shape is incompatible before full solver setup. Check constitutive model parameter expressions instead — these are always safe to inspect. Underworld development team with AI support from Claude Code (https://claude.com/claude-code)
lmoresi
force-pushed
the
feature/region-ds-cell-labels
branch
from
April 22, 2026 12:07
7619437 to
7dc50d4
Compare
Member
Author
|
Note - most of the commits above are a result of rebasing to updated development. |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Infrastructure for subdomain solving via PETSc
DMPlexFiltersubmeshes, region labels on all internal-boundary mesh types, projected P1 boundary normals, and mixed-mesh safety checks.Submesh extraction and data transfer
Mesh.extract_region(label, value)— creates a submesh that shares nodes with the parent via a subpoint IS.restrict/prolongate— IS-based parent↔submesh data transfer (machine-precision round-trip).copy_into/add_intoonEnhancedMeshVariable— user-friendly wrappers over the IS transfer._check_expression_meshesin solvers andextract_meshes()utility;uw.function.evaluate()gains the same check. Protects against silently using variables from the wrong mesh in an expression.Mesh adaptation
_re_extract_from_parent).evaluate.sync_coordinates_from_parentfor deformation.Region labels on internal-boundary meshes
AnnulusInternalBoundary: two surfaces sharing the internal curve loop, region Physical Groups,regionsenum, HDF5 save/restore.SphericalShellInternalBoundary: OCC fragment with bbox-radius classification (Inner/Outer).BoxInternalBoundary: centroid-based classification (cells belowzintCoord→Inner, above →Outer).Projected P1 boundary normals
mesh.Gamma_P1: projected P1 normals fromGamma_N. UnlikeGamma_N(which flips sign on inner boundaries and scales with mesh size),Gamma_P1is consistent and converges under refinement.Gamma_P1— resolution study shows Nitsche+Gamma_P1converges (7.39 → 7.16 → 6.92) where penalty+Gamma_Ndiverges (1.46 → 1.89 → 2.16).add_nitsche_bcgains amaskparameter for selective application.Tests and demos
Design docs
docs/developer/design/submesh-solver-architecture.md.Test plan
test_0502_boundary_integrals.py) passNotes for reviewers
boundary_normalsonSphericalShellInternalBoundaryis deprecated in favour ofGamma_P1.Underworld development team with AI support from Claude Code