Draft UWTN 2026-016: Boundary conditions on non-planar boundaries

.gitignore

@@ -54,3 +54,8 @@ articles/*/*.pdf
 # macOS
 .DS_Store
 
+
+# Underworld writes a mesh cache beside the script that asked for a mesh. It is
+# regenerated by running the example, and one note's examples committed 650 kB
+# of it before this line existed.
+.meshes/

WRITING-PLAN.md

@@ -220,9 +220,9 @@ beginning of a solvers paper, but neither is being written to fit one.
 
 ### R1. Boundary conditions on non-planar boundaries
 
-Status: not started. **One post**, not a series. Rescoped 2026-08-17: rotated
-boundary conditions are the answer, but the *question* is the better frame, and
-it is the one a reader arrives with.
+Status: drafted as UWTN 2026-016, in review on PR #26. **One post**, not a
+series. Rescoped 2026-08-17: rotated boundary conditions are the answer, but the
+*question* is the better frame, and it is the one a reader arrives with.
 
 On a box, "no flow through this wall" is a component of the velocity and you
 constrain it. On an annulus, a sphere, a boundary with topography, or any mesh
@@ -451,6 +451,6 @@ Listed as candidates, not commitments.
 | 6 | F5 The comparison | Last, with Thyagarajulu's benchmarks as its evidence |
 | 7 | S1, S2 free surface | After the discussion about how to split it |
 | — | C1 Launching from any repository | Standalone; the capability is live and undocumented |
-| — | R1 BCs on non-planar boundaries | Standalone; write whenever it suits. Much of the evidence exists — see the rescoped entry |
+| — | R1 BCs on non-planar boundaries | Drafted as UWTN 2026-016, in review on PR #26 |
 | ✓ | G1 Setting up FMG | Published 2026-08-17 (UWTN 2026-014) |
 | 8 | #1 Release announcement | Written last; links to everything |

articles/boundary-conditions-on-non-planar-boundaries/examples/README-stress-test.md

@@ -0,0 +1,49 @@
+# The examples, and what each one measures
+
+Every table and figure in the note is produced by one of these. Run them from
+this directory with an Underworld3 environment on the path.
+
+| script | what it produces |
+|---|---|
+| `leak.py` | the constraint tables — how much flow each treatment lets through, under refinement (`sweep`) and against its own parameter (`params`) |
+| `stress.py` | the surface-stress tables against the exact `uw.analytic.CylindricalStokes` answer: `sweep`, `params`, `locking`, `control` |
+| `solcx.py` | the lateral-viscosity half, against `uw.analytic.SolCx`: `sweep`, `contrast`, `params`, `control` |
+| `generate-locking-figure.py` | `figures/locking.png` and `figures/banner.png` |
+| `generate-rotated-basis.py` | the data behind `rotated-basis.typ`, which draws `figures/rotated-basis.svg` |
+
+## The traps, all paid for once
+
+- **`v.array` is `(N, 1, dim)`.** It broadcasts silently against `(N, dim)`
+  normals and returns ~1e-16 projections for a 1e-2 velocity, with no error.
+  `np.squeeze` it. Every number in the first leak run was this artefact.
+- **Give every metric a negative control and run it.** Two free-slip circles
+  leave the rigid rotation unconstrained, and that nullspace is purely
+  tangential — so a radial leak metric read 2e-14 on a solve that had diverged
+  with `|u| = 2.7e5`. The inner boundary is now held.
+- **Check `snes.getConvergedReason() > 0` before tabulating anything.** Diverged
+  runs leave plausible numbers in the array; two nearly reached the note.
+- **Vertex against edge midpoint.** On a curved boundary, vertex values of
+  `sigma_nn` carry the O(h) facet error and midpoints are superconvergent
+  (underworld3#414). `stress.py` splits them and the note says which it uses.
+- **Select a boundary trace by the mesh LABEL, not by a radius band.** A band
+  that narrows with the mesh admits a different node set at each resolution, and
+  the earlier version of this comparison drifted by four nodes between methods
+  because of it.
+- **Enclosed domain**: the multiplier, the pressure and the traction are each
+  determined only up to a constant. Compare deviations, which is what topography
+  is anyway.
+- **A `Piecewise` viscosity inside a boundary penalty term does not solve** on
+  SolCx, at any magnitude tried.
+
+## What is still open
+
+- The spherical case. `uw.analytic.Zhong2008` publishes
+  `.response().surface_topography = 0.4191904156575601` for the default
+  degree-2 case (load r = 0.775, r_inner = 0.55, isoviscous), which is the right
+  oracle for dynamic topography proper. It is a 3-D shell, and that is the cost
+  step this note stopped short of.
+- Behr (2004) reports non-physical recirculation at curved walls even with the
+  consistent normal. We have not looked for it here.
+- underworld3 issues #607 (the multiplier misses the augmented-Lagrangian share)
+  and #608 (the rotated constraint at a corner shared with a component
+  condition) both came out of these runs and are open.

articles/boundary-conditions-on-non-planar-boundaries/examples/generate-locking-figure.py

@@ -0,0 +1,106 @@
+"""The figure for "A constraint that is satisfied, and wrong".
+
+Three solves of the same annulus problem, rendered with the same colour scale:
+the exact solution, a direct penalty written against the facet normal, and the
+same penalty written against the measure-weighted node normal. The coefficient
+is 1e6 in both penalty panels, so the only difference between them is the
+normal.
+
+Writes:
+    figures/locking.png   the three-panel figure
+    figures/banner.png    a wide two-panel crop for the article banner
+
+    python3 generate-locking-figure.py
+
+Run against underworld3 `bugfix/multiplier-traction` (PR #617); the
+constrained solver's `traction()` is the fix this note prompted.
+"""
+import pathlib
+
+import numpy as np
+import pyvista as pv
+
+pv.OFF_SCREEN = True
+
+import underworld3 as uw
+import underworld3.visualisation as vis
+
+import stress as S
+
+CELL = 0.075
+PENALTY = 1.0e6
+HERE = pathlib.Path(__file__).resolve().parent
+FIGURES = HERE.parent / "figures"
+
+
+def speed_mesh(mode):
+    """(pyvista mesh carrying |u|, mesh edges, max speed) for one treatment."""
+    mesh, stokes, v, exact = S.build(mode, cell=CELL, penalty=PENALTY)
+    stokes.solve()
+    assert S.converged(stokes), "%s did not converge" % mode
+    pv_v = vis.meshVariable_to_pv_mesh_object(v)
+    u = np.squeeze(np.asarray(v.array))
+    pv_v.point_data["speed"] = np.linalg.norm(u, axis=1)
+    edges = vis.mesh_to_pv_mesh(mesh).extract_all_edges()
+    return pv_v, edges, exact, v
+
+
+def exact_speed_mesh(reference):
+    """The same object built from the exact velocity, for the first panel."""
+    mesh, stokes, v, exact = reference
+    field = uw.discretisation.MeshVariable("Uplot", v.mesh, v.mesh.dim, degree=2)
+    field.data[:] = exact.evaluate("velocity", field.coords)
+    pv_v = vis.meshVariable_to_pv_mesh_object(field)
+    pv_v.point_data["speed"] = np.linalg.norm(
+        np.squeeze(np.asarray(field.array)), axis=1)
+    return pv_v
+
+
+def panel(plotter, index, pv_mesh, edges, title, clim, zoom=1.3):
+    plotter.subplot(0, index)
+    plotter.set_background("white")
+    plotter.add_mesh(pv_mesh, scalars="speed", cmap="RdBu_r", clim=clim,
+                     show_edges=False, lighting=False, show_scalar_bar=False)
+    plotter.add_mesh(edges, color="black", line_width=0.4, lighting=False)
+    plotter.add_text(title, position="upper_left", font_size=12, color="black")
+    plotter.view_xy()
+    plotter.camera.zoom(zoom)
+
+
+def main():
+    facet_pv, facet_edges, exact, facet_v = speed_mesh("penalty")
+    node_pv, node_edges, _exact, node_v = speed_mesh("penalty_node")
+
+    # The exact field on the node-normal run's mesh, which is the same mesh.
+    truth = uw.discretisation.MeshVariable("Utruth", node_v.mesh, node_v.mesh.dim,
+                                           degree=2)
+    truth.data[:] = _exact.evaluate("velocity", truth.coords)
+    truth_pv = vis.meshVariable_to_pv_mesh_object(truth)
+    truth_pv.point_data["speed"] = np.linalg.norm(
+        np.squeeze(np.asarray(truth.array)), axis=1)
+
+    top = float(truth_pv.point_data["speed"].max())
+    clim = (0.0, top)
+    print("colour scale 0 to %.4e" % top)
+    for name, pv_mesh in (("exact", truth_pv), ("facet", facet_pv), ("node", node_pv)):
+        print("%-6s max speed %.4e" % (name, float(pv_mesh.point_data["speed"].max())))
+
+    FIGURES.mkdir(exist_ok=True)
+
+    plotter = pv.Plotter(off_screen=True, shape=(1, 3), window_size=(1650, 620))
+    panel(plotter, 0, truth_pv, node_edges, "exact", clim)
+    panel(plotter, 1, facet_pv, facet_edges, "penalty, facet normal", clim)
+    panel(plotter, 2, node_pv, node_edges, "penalty, node normal", clim)
+    plotter.screenshot(str(FIGURES / "locking.png"))
+    plotter.close()
+
+    banner = pv.Plotter(off_screen=True, shape=(1, 2), window_size=(1600, 560))
+    panel(banner, 0, facet_pv, facet_edges, "", clim, zoom=1.9)
+    panel(banner, 1, truth_pv, node_edges, "", clim, zoom=1.9)
+    banner.screenshot(str(FIGURES / "banner.png"))
+    banner.close()
+    print("wrote", FIGURES / "locking.png", "and", FIGURES / "banner.png")
+
+
+if __name__ == "__main__":
+    main()

articles/boundary-conditions-on-non-planar-boundaries/examples/generate-rotated-basis.py

@@ -0,0 +1,102 @@
+"""Geometry for the rotated-boundary-conditions figure.
+
+The cetz skill's rule: geometry computation happens in Python and Typst just
+draws. Doing the triangulation here rather than in the figure is what stops
+nodes being left out of the mesh -- an earlier version connected nodes by a
+distance threshold and silently missed several.
+
+Emits the schema in the skill's `underworld-bridge.md`:
+
+    {"vertices": [[x, y], ...],
+     "triangles": [[i, j, k], ...],
+     "surface": [i, ...],              indices of the constrained nodes
+     "frames": [{"p": [x, y], "n": [nx, ny], "t": [tx, ty]}, ...],
+     "curve": [[x, y], ...]}           the surface, finely sampled
+
+Run:  python3 generate-rotated-basis.py
+"""
+import json
+import pathlib
+
+import numpy as np
+from scipy.spatial import Delaunay
+
+OUT = pathlib.Path(__file__).with_name("rotated-basis-data.json")
+
+X0, X1 = -3.5, 3.5
+BASE = -2.9
+NX = 9          # columns of nodes
+NY = 4          # rows, surface included
+
+
+def surface_y(x):
+    """A deformed surface: rises on the left, falls on the right, with an
+    inflection between, so the normal swings through a wide range and the
+    curvature changes sign. No global rotation straightens this out, which is
+    the reason the figure exists."""
+    x = np.asarray(x, dtype=float)
+    return (0.95 * np.exp(-(((x + 1.75) / 1.30) ** 2))
+            - 0.80 * np.exp(-(((x - 1.70) / 1.15) ** 2))
+            + 0.75)
+
+
+def surface_slope(x, eps=1.0e-4):
+    return (surface_y(x + eps) - surface_y(x - eps)) / (2 * eps)
+
+
+# Nodes on a grid warped to sit under the surface. Columns are staggered on
+# alternate rows so the Delaunay triangulation comes out as triangles rather
+# than as near-degenerate right angles on a perfect lattice.
+xs = np.linspace(X0, X1, NX)
+pts = []
+surface_idx = []
+for row in range(NY):
+    frac = row / (NY - 1)                      # 0 at the base, 1 at the surface
+    offset = 0.0 if row % 2 == 0 else 0.5 * (xs[1] - xs[0])
+    cols = xs + offset
+    if row == NY - 1:
+        cols = xs                              # surface row unstaggered
+    for x in cols:
+        if x < X0 - 1e-9 or x > X1 + 1e-9:
+            continue
+        top = float(surface_y(x))
+        y = BASE + (top - BASE) * frac
+        if row == NY - 1:
+            surface_idx.append(len(pts))
+        pts.append([float(x), float(y)])
+
+pts = np.array(pts)
+tri = Delaunay(pts)
+
+# Drop the slivers Delaunay leaves along a non-convex top edge: any triangle
+# whose centroid sits above the surface is outside the domain.
+keep = []
+for simplex in tri.simplices:
+    c = pts[simplex].mean(axis=0)
+    if c[1] <= float(surface_y(c[0])) + 1.0e-9:
+        keep.append([int(i) for i in simplex])
+
+frames = []
+for i in surface_idx:
+    x = pts[i][0]
+    m = float(surface_slope(x))
+    n = np.array([-m, 1.0])
+    n /= np.linalg.norm(n)
+    t = np.array([1.0, m])
+    t /= np.linalg.norm(t)
+    frames.append({"p": [float(pts[i][0]), float(pts[i][1])],
+                   "n": [float(n[0]), float(n[1])],
+                   "t": [float(t[0]), float(t[1])]})
+
+curve_x = np.linspace(X0, X1, 121)
+data = {
+    "vertices": [[float(a), float(b)] for a, b in pts],
+    "triangles": keep,
+    "surface": [int(i) for i in surface_idx],
+    "frames": frames,
+    "curve": [[float(a), float(b)] for a, b in zip(curve_x, surface_y(curve_x))],
+}
+OUT.write_text(json.dumps(data, indent=1))
+print("wrote %s: %d vertices, %d triangles, %d surface nodes"
+      % (OUT.name, len(data["vertices"]), len(data["triangles"]),
+         len(data["surface"])))