fix(mesh): boundary normals and domain membership must track a deformed mesh#264
Open
lmoresi wants to merge 2 commits into
Open
fix(mesh): boundary normals and domain membership must track a deformed mesh#264lmoresi wants to merge 2 commits into
lmoresi wants to merge 2 commits into
Conversation
…ed mesh Two related bugs where analytic/cached UNDEFORMED geometry survived a mesh.deform(), corrupting Nitsche/penalty free-slip BCs and semi-Lagrangian advection on free surfaces. (1) Projected boundary normals stayed radial on a deformed surface. mesh.Gamma's base scalars carry the C-code petsc_n[] (PETSc facet normal), defined only inside surface-integral kernels; _update_projected_normals point-evaluated it, so it fell back to the coordinate -> (x,y)/r = radial. Gamma_P1 therefore never tracked deformation and add_nitsche_bc constrained v.r_hat instead of v.n_hat_true, leaking throughflow proportional to surface tilt. Fix: new Mesh.boundary_normal(label) assembles the exact PETSc facet normals (computeCellGeometryFVM, outward, area-weighted) from ONLY that boundary's faces -> per-boundary (a corner shared with another boundary keeps its own one-sided normal, never an average across the discontinuity), deformation-tracking, smooth on a smooth boundary. add_nitsche_bc (Stokes and Scalar) now uses it; global Gamma_P1 is unchanged for back-compat. The field is refreshed on every deform so BCs that captured its .sym at setup read the new geometry. Parallel-safe (guards getStratumIS against dead label values on ranks holding no faces of a boundary). (2) points_in_domain()/function.evaluate() used the original boundary. The boundary-skeleton kd-tree was built from mesh._nav_coords, captured as a reference to the ORIGINAL coords in __init__ and never refreshed in nuke_coords_and_rebuild (only the per-cell control-point arrays were, on adapt). On a volume mesh _nav_coords stayed at the undeformed boundary, so a region that bulged past the original surface read as EXTERIOR -- stranding SL trace-back feet there and mis-locating evaluations, injecting the cold boundary value at topographic highs. Fix: invalidate the boundary kd-tree and refresh _nav_coords from the current DM coordinates on every deform. Tests: test_0056 (boundary_normal tracks a deformed annulus + a Cartesian free-surface corner) and test_0057 (deformed-domain membership + evaluation), both tier_a. test_1060 Nitsche free-slip passes serial and np2. Follow-up: ADD-reduce the unnormalised facet contributions at partition-seam boundary vertices before normalising (those few vertices currently use a this-rank-only stencil; parallel-safe and correct elsewhere). Underworld development team with AI support from Claude Code
Contributor
There was a problem hiding this comment.
Pull request overview
Fixes two deformation-related geometry caching bugs in Mesh: (1) free-slip/Nitsche boundary normals now track mesh deformation per-boundary (and remain corner-correct), and (2) domain membership / point-location now uses the deformed boundary geometry by invalidating and rebuilding navigation caches on deform().
Changes:
- Add
Mesh.boundary_normal(boundary)that assembles deformation-tracking PETSc facet normals per boundary and update Nitsche BCs to use it by default. - Fix deformed-domain membership / evaluation by invalidating the boundary kd-tree and refreshing
_nav_coordsafter coordinate rebuilds. - Add regression tests for deformed normals and deformed domain membership/evaluation.
Reviewed changes
Copilot reviewed 4 out of 4 changed files in this pull request and generated 6 comments.
| File | Description |
|---|---|
src/underworld3/discretisation/discretisation_mesh.py |
Adds per-boundary normal assembly, refreshes navigation coords and invalidates boundary control-point kd-tree on deform/rebuild. |
src/underworld3/cython/petsc_generic_snes_solvers.pyx |
Switches default Nitsche BC normal from Gamma_P1 to deformation-tracking boundary_normal. |
tests/test_0056_projected_normals_deform.py |
Regression test ensuring per-boundary normals track deformation and are corner-correct. |
tests/test_0057_deformed_domain_membership.py |
Regression test ensuring points_in_domain() and function.evaluate() track deformed geometry. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Comment on lines
+2509
to
+2511
| import underworld3 as uw | ||
| from scipy.spatial import cKDTree | ||
|
|
Comment on lines
+2560
to
+2580
| tree = cKDTree(coords) | ||
| nverts = cdim # P1 vertices per facet (2D edge=2, 3D tri=3) | ||
| for f in face_pts: | ||
| if dm.getSupportSize(f) != 1: | ||
| continue | ||
| area, cent, nrm = dm.computeCellGeometryFVM(f) | ||
| nrm = numpy.asarray(nrm)[:cdim] | ||
| cell = dm.getSupport(f)[0] | ||
| _, ccent, _ = dm.computeCellGeometryFVM(cell) | ||
| if numpy.dot(nrm, numpy.asarray(cent)[:cdim] | ||
| - numpy.asarray(ccent)[:cdim]) < 0: | ||
| nrm = -nrm | ||
| # Accumulate to the facet's P1 DOFs (its vertices) — found as the | ||
| # nearest DOFs to the facet centroid. This avoids indexing the local | ||
| # coordinate array by (vertex_point - vStart), which is only valid | ||
| # in serial (the parallel coordinate layout differs → out-of-range). | ||
| # Normalisation at the end makes the per-DOF weight (full vs share) | ||
| # irrelevant to the resulting direction. | ||
| _, idxs = tree.query(numpy.asarray(cent)[:cdim], k=nverts) | ||
| for idx in numpy.atleast_1d(idxs): | ||
| accum[idx] += area * nrm |
Comment on lines
+150
to
+154
| # normal, NOT a 45° average with the vertical walls. ∂h/∂x at x=0,1 is | ||
| # ±0.1π, so the true top normal there is tilted ~17° from vertical — and | ||
| # crucially has a SMALL... no: it must equal the one-sided top facet normal | ||
| # (the reference), i.e. err already bounds it. Assert the corner normal is | ||
| # not the 45° wall-average: |
Comment on lines
+4775
to
+4789
| # Surface normal components. By default use this boundary's OWN | ||
| # deformation-tracking facet normal (assembled from only this | ||
| # boundary's faces, so a corner shared with another boundary is not | ||
| # averaged across the discontinuity, and a deformed surface is | ||
| # followed correctly). The legacy global mesh.Gamma_P1 point-evaluates | ||
| # petsc_n and stays radial on a deformed surface — see | ||
| # Mesh.boundary_normal. | ||
| if normal is not None: | ||
| if isinstance(normal, sympy.MatrixBase): | ||
| n = [normal[i] for i in range(dim)] | ||
| else: | ||
| n = list(normal) | ||
| else: | ||
| Gamma_P1 = mesh.Gamma_P1 | ||
| n = [Gamma_P1[i] for i in range(dim)] | ||
| bnorm = mesh.boundary_normal(boundary) | ||
| n = [bnorm[i] for i in range(dim)] |
Member
Author
|
This fix is a critical component of the adaptive-meshing / free-surface work, not just the free-slip BC path. While validating mmpde adaptive convection we traced the long-standing "holes / giant empty cells in the adapted mesh + wrecked adaptive convection" failure to Bug 2 here (deformed-mesh point location). The chain:
With this PR, the metric evaluates correctly/positive on the deformed mesh, so the mover is fed a valid metric and adaptation is stable. Validation (with this fix): forced adaptation every step (39/39) at R=5 under vigorous stagnant-lid convection (Ra1e6, Δη1e3) → mesh clean every frame (folded=0, cell-area-ratio ~14 flat), vrms→18, Nu→1.86. The same configuration without the fix produced area-ratio spikes 189/96/875 (holes). So this unblocks adaptive convection (and, as noted, the free-surface dynamic-topography work). A small complementary hardening (mover bakes the metric from positive nodal values via monotone RBF/Shepard, so it can never go negative regardless of eval precision) will follow as a separate bugfix PR. |
- boundary_normal(): drop unused cKDTree import (tree is built in _assemble_boundary_normal). - _assemble_boundary_normal(): drop unused `import underworld3 as uw`; count P1 vertices PER FACET from its closure (vStart..vEnd) instead of assuming nverts == cdim, so non-simplex facets (3D quad/hex boundaries) are handled correctly. - add_nitsche_bc docstring: `normal` default is now the per-boundary, deformation-tracking mesh.boundary_normal(boundary), not mesh.Gamma_N. - tests: drop unused surf_now; tidy an editorial comment fragment. test_0056 + test_0057 pass (3 passed). Underworld development team with AI support from Claude Code
Member
Author
|
Addressed the Copilot review in e086e72:
|
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.
Two related bugs where analytic/cached undeformed geometry survived a
mesh.deform(), corrupting Nitsche/penalty free-slip BCs and semi-Lagrangian advection on free (deforming) surfaces. Surfaced while building a free-slip + stress-driven dynamic-topography free-surface convection scheme.Bug 1 — projected boundary normals stayed radial on a deformed surface
mesh.Gamma's base scalars carry the C-codepetsc_n[](the PETSc facet normal), which is only defined inside surface-integral kernels._update_projected_normalspoint-evaluated it, so it fell back to the coordinate →(x,y)/r= radial.Gamma_P1thus never tracked deformation andadd_nitsche_bcconstrainedv·r̂instead ofv·n̂_true, leaking throughflow ∝ surface tilt (verified ~9° off at 15% topography).Fix: new
Mesh.boundary_normal(label)assembles the exact PETSc facet normals (computeCellGeometryFVM, outward, area-weighted) from only that boundary's faces → per-boundary (a corner shared with another boundary keeps its own one-sided normal, never an average across the discontinuity — handles the box free-slip/no-slip corner and a Cartesian free-surface top-meets-wall corner), deformation-tracking, smooth on a smooth boundary.add_nitsche_bc(Stokes + Scalar) now uses it; the globalGamma_P1is unchanged for back-compat. The field is refreshed on everydeformso BCs that captured its.symat setup read the new geometry. Parallel-safe (guardsgetStratumISagainst dead label values on ranks holding no faces of a boundary).Bug 2 —
points_in_domain()/function.evaluate()used the original boundaryThe boundary-skeleton kd-tree is built from
mesh._nav_coords, which was captured as a reference to the original coords in__init__and never refreshed innuke_coords_and_rebuild(only the per-cell control-point arrays were, onadapt). On a volume mesh_nav_coordsstayed at the undeformed boundary, so a region that bulged past the original surface (r > r_oon a free surface) read as exterior — stranding SL trace-back feet there and mis-locating evaluations to the old boundary, injecting the cold boundary value at topographic highs.Fix: invalidate the boundary kd-tree and refresh
_nav_coordsfrom the current DM coordinates on everydeform.Tests
test_0056—boundary_normaltracks a deformed annulus + a Cartesian free-surface corner (tier_a).test_0057— deformed-domain membership + evaluation (tier_a).test_1060Nitsche free-slip passes serial and np2; broad serial regression (meshes, pointwise, evaluate, slip, Nitsche) green.Known follow-up
ADD-reduce the unnormalised facet contributions at partition-seam boundary vertices before normalising — those few vertices currently use a this-rank-only stencil (parallel-safe and correct everywhere else; flagged with a
TODOin_assemble_boundary_normal).Underworld development team with AI support from Claude Code