Normalise Gamma_N for consistent penalty/Nitsche scaling - #109
Conversation
mesh.Gamma_N now returns Gamma / |Gamma| — a unit normal regardless of element size. mesh.Gamma remains the raw (un-normalised) PETSc face normal whose magnitude scales with edge length (2D) / face area (3D). This fixes a scaling issue where: - penalty * Gamma.dot(v) * Gamma had effective penalty ~ penalty * h² - Nitsche gamma/h * Gamma.dot(v) * Gamma scaled as h (should be 1/h) With normalised Gamma_N, both penalty and Nitsche terms are now mesh-independent at the symbolic level. Updated Nitsche BC code to use the normalised matrix form. JIT extension updated to register ccode on the raw Gamma base scalars. Underworld development team with AI support from Claude Code
There was a problem hiding this comment.
Pull request overview
This PR normalizes the mesh boundary facet normal used in boundary-integral formulations so penalty and Nitsche terms scale consistently across mesh resolution (removing implicit element-size dependence).
Changes:
- Update
mesh.Gamma_Nto return a unit normal (Gamma / |Gamma|) while keepingmesh.Gammaas the raw PETSc facet normal. - Update Nitsche BC implementations to use the normalized normal components from the new
Gamma_Nrepresentation. - Adjust JIT ccode registration to target raw
Gammabase scalars rather thanGamma_Ncomponents.
Reviewed changes
Copilot reviewed 3 out of 3 changed files in this pull request and generated 1 comment.
| File | Description |
|---|---|
src/underworld3/utilities/_jitextension.py |
Ensures ccode printing is registered for raw facet-normal base scalars used by normalized expressions. |
src/underworld3/discretisation/discretisation_mesh.py |
Implements normalized Gamma_N and clarifies Gamma as the un-normalized facet normal. |
src/underworld3/cython/petsc_generic_snes_solvers.pyx |
Updates Nitsche BC assembly to pull normal components via matrix indexing from mesh.Gamma_N. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| def Gamma_N(self) -> sympy.Matrix: | ||
| r"""Normalised boundary/surface normal as a row matrix. | ||
|
|
||
| Returns ``Gamma / |Gamma|`` so that the result is a unit normal | ||
| regardless of element size. Use this for penalty and Nitsche BCs | ||
| where mesh-independent scaling is required. | ||
|
|
||
| Returns | ||
| ------- | ||
| sympy.vector.CoordSys3D | ||
| The boundary coordinate system object. | ||
| sympy.Matrix | ||
| Row matrix of normalised boundary normal components. | ||
| """ | ||
| return self._Gamma | ||
| G = self.Gamma | ||
| return G / sympy.sqrt(G.dot(G)) |
There was a problem hiding this comment.
mesh.Gamma_N changes from returning the raw CoordSys3D (self._Gamma) to returning a normalized sympy.Matrix. This is an API-breaking change for any downstream code that previously used mesh.Gamma_N.x/.y/.z. If backward compatibility is required, consider adding a separate public property for the raw coordinate system (e.g., Gamma_CS/Gamma_raw_cs) or documenting a deprecation path and the new component access pattern explicitly.
There was a problem hiding this comment.
Acknowledged — this is intentional. The old CoordSys3D return type was misleading: Gamma_N implied a normalised normal but returned the raw un-normalised coordinate system. The name/contract now matches the implementation.
Internal call sites (Nitsche BC in both SNES_Vector and SNES_Stokes_SaddlePt) are updated. User code that accessed mesh.Gamma_N.x should use mesh.Gamma_N[0] (matrix indexing) or switch to mesh.Gamma for the raw normal.
The docstring documents the change and the relationship to mesh.Gamma.
|
Minor "bug-fix" or, more correctly, updating the interface so that Gamma always represents a unit vector. The petsc version is scaled for the surface area of the surface segment. Fine, but the symbolic versions just end up being unit vectors. This is the easiest way to make everything match. |
|
For copilot, this is a breaking change in that it breaks the bug that we had been seeing. |
Summary
mesh.Gamma_Nnow returnsGamma / |Gamma|— a unit normal regardless of element sizemesh.Gammaremains the raw un-normalised PETSc face normal (magnitude ~ edge length in 2D, face area in 3D)Why
The un-normalised
Gammain penalty BCs gives an effective penalty that scales with h²:penalty * Gamma.dot(v) * Gamma→ effective penalty ~ penalty * h²gamma * mu / hcombined with un-normalised Gamma scaled as h instead of 1/hWith normalised
Gamma_N, penalty and Nitsche terms are mesh-independent. Users should adjust penalty magnitude (no implicit h-scaling).Test plan
Gamma_Ngive mesh-independent constraint strengthUnderworld development team with AI support from Claude Code