From 90123353c575dfcfd0b591aca7ae8b61300d4e21 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 7 Apr 2026 13:54:49 -0700 Subject: [PATCH] Normalise Gamma_N: unit boundary normal for consistent penalty scaling MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 --- .../cython/petsc_generic_snes_solvers.pyx | 14 +++++------- .../discretisation/discretisation_mesh.py | 22 +++++++++++++------ src/underworld3/utilities/_jitextension.py | 6 +++-- 3 files changed, 25 insertions(+), 17 deletions(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 62c411630..88d09c7f7 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -2193,10 +2193,9 @@ class SNES_Vector(SolverBaseClass): mesh = self.mesh dim = mesh.dim - # Surface normal components - n = [mesh.Gamma_N.x, mesh.Gamma_N.y] - if dim == 3: - n.append(mesh.Gamma_N.z) + # Surface normal components (normalised) + Gamma_N = mesh.Gamma_N + n = [Gamma_N[i] for i in range(dim)] # Constraint direction: defaults to surface normal if direction is not None: @@ -3233,16 +3232,15 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): mesh = self.mesh dim = mesh.dim - # Surface normal components. By default use PETSc's facet normal. + # Surface normal components. By default use normalised PETSc facet normal. if normal is not None: if isinstance(normal, sympy.MatrixBase): n = [normal[i] for i in range(dim)] else: n = list(normal) else: - n = [mesh.Gamma_N.x, mesh.Gamma_N.y] - if dim == 3: - n.append(mesh.Gamma_N.z) + Gamma_N = mesh.Gamma_N + n = [Gamma_N[i] for i in range(dim)] # Constraint direction: defaults to surface normal if direction is not None: diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 4b11dbd6e..7802dd938 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -1451,19 +1451,27 @@ def N(self) -> sympy.vector.CoordSys3D: return self._N @property - def Gamma_N(self) -> sympy.vector.CoordSys3D: - r"""SymPy coordinate system for boundary/surface coordinates. + 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)) @property - def Gamma(self) -> sympy.vector.CoordSys3D: - r"""Boundary coordinate scalars as a row matrix. + def Gamma(self) -> sympy.Matrix: + r"""Raw (un-normalised) boundary coordinate scalars as a row matrix. + + The magnitude scales with face edge length (2D) or face area (3D). + For a unit normal, use :attr:`Gamma_N` instead. Returns ------- diff --git a/src/underworld3/utilities/_jitextension.py b/src/underworld3/utilities/_jitextension.py index d9eb50559..a15e76778 100644 --- a/src/underworld3/utilities/_jitextension.py +++ b/src/underworld3/utilities/_jitextension.py @@ -775,8 +775,10 @@ def _basescalar_ccode(self, printer): return f"petsc_x[{idx}]" type(mesh.N.x)._ccode = _basescalar_ccode - if type(mesh.Gamma_N.x) is not type(mesh.N.x): - type(mesh.Gamma_N.x)._ccode = _basescalar_ccode + # Gamma base scalars (un-normalised face normal) — ensure ccode is registered + Gamma_scalars = mesh._Gamma.base_scalars() + if type(Gamma_scalars[0]) is not type(mesh.N.x): + type(Gamma_scalars[0])._ccode = _basescalar_ccode # Create a custom functions replacement dictionary. # Note that this dictionary is really just to appease Sympy,