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,