Skip to content

Normalise Gamma_N for consistent penalty/Nitsche scaling #109

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
lmoresi merged 1 commit into from
Apr 7, 2026
Merged

Normalise Gamma_N for consistent penalty/Nitsche scaling #109

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 6 additions & 8 deletions src/underworld3/cython/petsc_generic_snes_solvers.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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:
Expand Down
22 changes: 15 additions & 7 deletions src/underworld3/discretisation/discretisation_mesh.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Comment on lines +1454 to +1467

Copilot AI Apr 7, 2026

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copilot uses AI. Check for mistakes.

@lmoresi lmoresi Apr 7, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.


@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
-------
Expand Down
6 changes: 4 additions & 2 deletions src/underworld3/utilities/_jitextension.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
Loading