Skip to content

Adaptivity — h, r and p on one problem

Download .py All tutorials

jNO adapts a discretisation three ways, all through the same adapt= slot on fem.solve:

  • h — jno.solve.remesh() / refine(): add elements where the error is.
  • r — jno.solve.relocate(): move a fixed node set down a differentiable objective.
  • p — jno.solve.enrich(): raise the local order, leaving the mesh alone, by switching interpolation covers on at the marked nodes (space="cover").

They are being compared, so they all run on the same problem, from the same coarse mesh, against the same reference, measured by the same functional:

\[-\Delta u = f \quad\text{on the L-shape},\qquad u=0 \ \text{on } \partial\Omega\]

with \(f\) a compact \(C^3\) bump in the lower-right arm. That carries two features of different character at once — the re-entrant corner's \(r^{2/3}\) singularity, and a smooth localized bump — so one run exercises both regimes. The boundary condition is homogeneous on purpose: a cover field's trace is only the P1 interpolant of an inhomogeneous \(g\), so \(u=0\) keeps that documented scope limit out of a comparison it would otherwise distort.

The solution with its corner and bump; the h-refined mesh; the enriched node set; and all four runs on one error-versus-DOF axis.

The result

method error \(E_\text{ref}-E\) active DOFs
coarse start 5.725e-03 144
h 6.727e-04 561 8.5× lower
r 9.147e-03 144 (+0) worse than the start
p 2.147e-04 313 (50 % enriched) 26.7× lower
h then p 5.773e-05 739 (70 %) 99.2× lower

p reaches 3.1× lower error than h with 56 % of the DOFs, and composing the two beats either alone.

The r row is worse than the start, and the objective is why

Relocation descends a mesh functional, and the default one — arclength equidistribution — targets resolution, not this error. It keeps the mesh sane (smallest angle 40.8° → 31.8°) and does not help here.

objective="energy" descends the Ritz functional \(J(v)=\tfrac12 a(v,v)-(f,v)\) instead. Galerkin orthogonality gives \(J_h-J_\text{exact}=\tfrac12\|u-u_h\|_E^2\), so lowering \(J_h\) lowers the error itself: on this problem it cuts it to 0.70× at fixed DOFs (smallest angle 42.7° → 28.1°), against 1.43× for the default. That needs the bump integrated accurately, jno.fem(..., quad_degree=4). At the default degree \(J_h\) is computed with the quadrated load, and the descent lowered it by moving vertices so the quadrature over-counted the bump. It landed above the converged energy, an improvement that exists only in the quadrature.

(Until #114 objective="energy" descended \(a(u_h,u_h)\) whatever the load. With a source that is \(-2J_h\), so it walked away from the solution: the true error rose 3.6× and the smallest angle fell to 3.2°.)

Where each method spent its DOFs

This is the interesting half, and one problem with two features is what makes it visible:

p put its covers:       corner  75% | bump 100% | elsewhere  38%
h put its elements:     corner 0.0257 | bump 0.0135 | elsewhere 0.0289    (mean cell size)

Both concentrate on the bump, because at this resolution the bump dominates the error budget — the smooth feature is simply where the error is. What separates them is the price: p buys its accuracy at 1.5 extra DOFs per marked node, h by adding nodes and elements outright.

The corner is what neither fully resolves alone, and it is why h then p wins: refine the mesh where regularity is the limit, then raise the order where smoothness pays.

Choosing between them

the solution is… reach for why
singular (corner, crack, shock) h order cannot buy a rate the regularity does not support
smooth, resolved, the feature moves r no new DOFs; the mesh follows the feature
smooth, localized, under-resolved p order is far cheaper per DOF than nodes

The first row is measured, not asserted: on the pure singular-mode L-shape (all error at the corner, nothing smooth to enrich) h cuts the energy error by 72 % while p manages 4 % for 75 % more DOFs. \(r^{2/3}\notin H^2\), and a higher-order space buys its rate from smoothness the solution does not have.

What to notice

  • All three are one argument. remesh(), relocate() and enrich() all return an AdaptSpec for the same adapt= slot; the weak form — the physics — is untouched by the choice.
  • The functional has to be space-agnostic. \(E=\tfrac12\int|\nabla u_h|^2\) is read off the assembled form, 0.5 * u @ fem.eval(stiff, u), so P1, a relocated mesh and an enriched space are measured the same way. A geometric formula on the vertex values is blind to p: a cover's coefficients are the local gradient, so they move \(u_h\) between nodes and barely move the nodal values. Measured that way, enriching 5 % of the nodes and enriching all of them are indistinguishable.
  • The reference must out-resolve everything measured against it. A P1 reference is not good enough here — an enriched run beats a much larger P1 reference outright, which shows up as a negative error. The reference is a fully enriched (third-order) solve, and an assert fails if it ever fails to bound a run.
  • theta means the same thing everywhere — Dörfler bulk marking, over cells for h and over nodes for p. For p it runs over the unenriched nodes only: splitting a cell drops its indicator, so an h-loop can re-rank the whole field each round, but enriching a node does not move a geometric criterion at all, and ranking globally would re-mark the same nodes forever.
  • n_dofs in a p-run's history is the ACTIVE count. The padded layout gives every node its cover slots whether or not it is enriched, so the total never changes; an unenriched node has its slots pinned, and max_dofs budgets against the free count.
  • A fresh enrich starts from plain P1. Calling it on a field that is already uniformly space="cover" therefore removes enrichment before adding it back selectively — the active DOF count goes down at that call (912 → 739 in the h-then-p run above). A fresh cover field carries no mask, and a loop that began fully enriched would have nothing left to select. Calling enrich a second time is the other case: it resumes from the mask the first call left, so two budgeted runs compose instead of the second undoing the first.
  • Two API frictions worth knowing. adapt= does not take the linear= slot on a steady problem, so a non-default linear solver is passed positionally as solve_fn; and an adapt= driver returns the vertex view, which for a cover field is only the value slots — re-solve on the final space (the loop leaves fem bound to it) when you need the full coefficient vector.

Full script

r"""**h**, **r** and **p** adaptivity on ONE problem — same mesh, same reference, same functional.

jNO adapts a discretisation three ways, all through the same ``fem.solve(adapt=...)`` slot:

* **h** — ``jno.solve.remesh()`` / ``refine()``: *add* elements where the error is.
* **r** — ``jno.solve.relocate()``: *move* a fixed node set down a differentiable objective.
* **p** — ``jno.solve.enrich()``: *raise the local order*, leaving the mesh alone, by switching
  interpolation covers on at the marked nodes (``space="cover"``).

They are measured against each other here, so they all run on the SAME problem: an L-shape driven by a
compact smooth source,

    -lap u = f  on the L-shape,   u = 0 on the whole boundary,
    f = a C^3 bump supported in the lower-right arm,

which carries **two** features of different character at once -- the re-entrant corner's ``r^(2/3)``
singularity, and a smooth localized bump. The boundary condition is homogeneous on purpose: a cover
field's trace is only the P1 interpolant of an INHOMOGENEOUS ``g``, so ``u = 0`` keeps that documented
scope limit out of a comparison it would otherwise distort.

Everything is measured with one functional, ``E = 1/2 integral |grad u_h|^2``, read straight off the
assembled form with ``fem.eval``, so it means the same thing for P1, for a relocated mesh and for an
enriched space. The geometric alternative -- a formula on the vertex values -- is blind to p: a cover's
coefficients are the local gradient, so they move ``u_h`` BETWEEN nodes and barely move the nodal
values. Measured that way, enriching 5% of the nodes and enriching all of them are indistinguishable.

With a source and homogeneous data the Galerkin solution minimises ``J = 1/2 a(v,v) - (f,v)`` and
``J_h = -E_h``, so the energy RISES toward the truth and the error is ``E_ref - E_h >= 0``.

Run::

    JAX_PLATFORMS=cpu pixi run -e fem python docs/tutorial_examples/08_fem_and_varpinns/adaptive_l_shape.py
"""

from __future__ import annotations

import os

import numpy as np

os.environ.setdefault("JAX_PLATFORMS", "cpu")  # small FEM solves; keep off the GPU

import jax  # noqa: E402

jax.config.update("jax_enable_x64", True)

import jax.numpy as jnp  # noqa: E402

import jno  # noqa: E402
from jno.utils.solver.fem_cover import cover_block  # noqa: E402

L_SHAPE = [(0, 0), (1.0, 0), (1.0, 0.5), (0.5, 0.5), (0.5, 1.0), (0, 1.0)]  # reentrant corner at (0.5, 0.5)
BUMP = (0.74, 0.22, 0.17, 60.0)  # cx, cy, R, amplitude -- the source, compactly supported inside the arm
CORNER = (0.5, 0.5)
H = 0.075  # the common coarse mesh every method starts from
MARGIN = 0.04
BLK = cover_block(2)  # 3 in 2-D: a cover node carries its value plus two cover coefficients


def source(x, y, xp):
    """A C^3 bump, ``(1-s)^4`` with ``s = |x-c|^2/R^2``, supported strictly inside the lower-right arm."""
    cx, cy, R, A = BUMP
    s = ((x - cx) ** 2 + (y - cy) ** 2) / R**2
    return A * xp.where(s < 1.0, 1.0 - s, 0.0) ** 4


def _movable(x, y):
    """Interior L-shape vertices: inside the material, off the outer box and the two notch edges."""
    inside = ~((x > 0.5) & (y > 0.5))
    off_box = (x > MARGIN) & (x < 1 - MARGIN) & (y > MARGIN) & (y < 1 - MARGIN)
    off_notch = ~(((jnp.abs(x - 0.5) < MARGIN) & (y > 0.5 - MARGIN)) | ((jnp.abs(y - 0.5) < MARGIN) & (x > 0.5 - MARGIN)))
    return inside & off_box & off_notch


def build(size, space="Lagrange", movable=False):
    d = jno.shape.polygon(L_SHAPE, size=size).domain()
    if movable:
        # `.trainable()` on a spatial coordinate turns that region's mesh VERTICES into a design
        # variable: the assembler routes them into the element geometry, so `fem.solve()` becomes
        # differentiable in the node positions. Literal, per component. Must be called BEFORE jno.fem.
        xm, ym, _ = d.variable("mov", where=_movable, split=True)
        xm.trainable(name="ix")
        ym.trainable(name="iy")
    u, phi = d.fem_symbols(space=space)
    xi, yi, _ = d.variable("interior", split=True)
    xb, yb, _ = d.variable("boundary", split=True)
    ui, vi = u.bind(x=xi, y=yi), phi.bind(x=xi, y=yi)
    stiff = ui.x * vi.x + ui.y * vi.y
    fem = jno.fem([stiff - source(xi, yi, jno.np) * vi, u(xb, yb) - 0.0])
    return d, fem, stiff


def solve_fn(a, b):
    """Sparse-direct, handed to ``adapt=`` positionally.

    ``adapt=`` does not take the ``linear=`` slot on a steady problem, so a non-default linear solver
    has to arrive as ``solve_fn``. The reference reaches ~12k DOFs, where a dense factorisation is a
    gigabyte; an enriched system is also the least well-conditioned of the three."""
    from jno.utils.solver.linear import sparse_lu_solve

    return sparse_lu_solve(a, b)


def energy(fem, stiff, sol):
    """``E = 1/2 integral |grad u_h|^2``, from the assembled form: ``1/2 u . (A u)``.

    Space-agnostic ON PURPOSE -- ``fem.eval`` integrates the weak form with whatever basis the field
    actually carries, so P1, a relocated mesh and an enriched space are measured by the same rule."""
    return 0.5 * float(np.dot(np.asarray(sol).reshape(-1), np.asarray(fem.eval(stiff, sol)).reshape(-1)))


def active_dofs(d, fem):
    """DOFs the solve actually carries. An unenriched cover node has its slots PINNED, so the padded
    size ``(1+dim) * n_nodes`` is the same every round and would hide the whole effect of p."""
    pins = getattr(d, "_fem_native_dirichlet_pairs", None) or []
    return int(fem.dofs) - len({int(i) for i, _ in pins})


def _min_detj(pts, cells):
    v = pts[cells]
    a, b = v[:, 1] - v[:, 0], v[:, 2] - v[:, 0]
    return float(np.min(a[:, 0] * b[:, 1] - a[:, 1] * b[:, 0]))


def _near(pts, centre, radius):
    return np.hypot(pts[:, 0] - centre[0], pts[:, 1] - centre[1]) < radius


def _min_angle(pts, cells):
    """Smallest angle in the mesh, in degrees -- the quality number relocation can quietly destroy."""
    v = pts[cells]
    a = np.linalg.norm(v[:, 1] - v[:, 0], axis=1)
    b = np.linalg.norm(v[:, 2] - v[:, 1], axis=1)
    c = np.linalg.norm(v[:, 0] - v[:, 2], axis=1)
    cos = np.stack(
        [
            (b * b + c * c - a * a) / (2 * b * c),
            (a * a + c * c - b * b) / (2 * a * c),
            (a * a + b * b - c * c) / (2 * a * b),
        ],
        axis=1,
    )
    return float(np.degrees(np.arccos(np.clip(cos, -1.0, 1.0))).min())


# --- the reference, and the common coarse start ----------------------------------------------------------
# The reference must OUT-RESOLVE everything measured against it. A P1 reference is not good enough for
# this comparison: an enriched run beats a P1 reference of several times the DOF count outright, which
# shows up as a negative error. A fully enriched (cover) reference is third order instead. The drift from
# h=0.020 to h=0.015 is 4.8e-06 -- well under the smallest error below, which is what this is quoted to.
d_ref, fem_ref, stiff_ref = build(0.015, space="cover")
E_REF = energy(fem_ref, stiff_ref, solve_fn(fem_ref.A, fem_ref.b))

d0, fem0, stiff0 = build(H)
E0, n0 = energy(fem0, stiff0, solve_fn(fem0.A, fem0.b)), active_dofs(d0, fem0)
pts0, tris0 = np.asarray(d0.mesh.points)[:, :2], np.asarray(d0.mesh.cells_dict["triangle"])

# (1) h -- ADD elements. One call runs solve -> ZZ estimate -> Dorfler mark -> local remesh.
d_h, fem_h, stiff_h = build(H)
fem_h.solve(solve_fn, adapt=jno.solve.remesh(theta=0.6, max_iters=4, refine_factor=1.7))
E_h, n_h = energy(fem_h, stiff_h, solve_fn(fem_h.A, fem_h.b)), active_dofs(d_h, fem_h)
pts_h, tris_h = np.asarray(d_h.mesh.points)[:, :2], np.asarray(d_h.mesh.cells_dict["triangle"])

# (2) r -- RELOCATE a fixed node set, down a mesh functional, through the differentiable solve.
#     The default objective (arclength equidistribution) targets RESOLUTION, not this error, and here
#     it does not help (measured 1.43x the coarse error). `objective="energy"` descends the Ritz
#     functional J = 1/2 a(u,u) - (f,u) instead, whose discrete value is the energy-norm error up to a
#     constant, and it cuts that error to 0.70x at fixed DOFs -- but only with the bump integrated
#     accurately, `jno.fem(..., quad_degree=4)`. At the default degree the descent lowers J by inflating
#     the QUADRATURE of the bump rather than the accuracy of u, and lands above the converged energy.
d_r, fem_r, stiff_r = build(H, movable=True)
pts_r0 = np.asarray(d_r.mesh.points)[:, :2].copy()  # coarse start, for the animation
sol_r = np.asarray(fem_r.solve(solve_fn, adapt=jno.solve.relocate(max_iters=30, lr=5e-4))).reshape(-1)
E_r, n_r = energy(fem_r, stiff_r, sol_r), active_dofs(d_r, fem_r)
pts_r, tris_r = np.asarray(d_r.mesh.points)[:, :2], np.asarray(d_r.mesh.cells_dict["triangle"])

# (3) p -- RAISE THE ORDER on the same mesh. The loop returns the VERTEX VIEW, which for a cover field
#     is only the value slots, so re-solve on the final space for the full coefficient vector (the loop
#     leaves `fem` bound to it -- one solve, not a rebuild).
d_p, fem_p, stiff_p = build(H, space="cover")
fem_p.solve(solve_fn, adapt=jno.solve.enrich(theta=0.8, max_iters=5))
E_p, n_p = energy(fem_p, stiff_p, solve_fn(fem_p.A, fem_p.b)), active_dofs(d_p, fem_p)
mask_p = np.asarray(d_p._fem_enriched_nodes, dtype=bool)
pts_p = np.asarray(d_p.mesh.points)[:, :2]

# (4) hp -- the two composed across successive solves: refine the mesh, then enrich on it.
#     NOTE a FRESH `enrich` starts from plain P1, so calling it on a field that was already uniformly
#     `space="cover"` -- as this one is, since h ran on it first -- REMOVES enrichment before adding it
#     back selectively, and the active DOF count goes DOWN at that call. A fresh cover field carries no
#     mask, and a loop that began fully enriched would have nothing left to select; enrich called a
#     SECOND time does resume from the mask it left.
d_hp, fem_hp, stiff_hp = build(H, space="cover")
fem_hp.solve(solve_fn, adapt=jno.solve.remesh(theta=0.6, max_iters=3, refine_factor=1.7))
n_hp_after_h = active_dofs(d_hp, fem_hp)
fem_hp.solve(solve_fn, adapt=jno.solve.enrich(theta=0.8, max_iters=4))
E_hp, n_hp = energy(fem_hp, stiff_hp, solve_fn(fem_hp.A, fem_hp.b)), active_dofs(d_hp, fem_hp)
mask_hp = np.asarray(d_hp._fem_enriched_nodes, dtype=bool)

err = lambda E: E_REF - E  # noqa: E731  -- energy rises toward the truth; see the module docstring
print(f"ONE problem, one mesh, one reference.   E_ref = {E_REF:.6f}   (cover, {active_dofs(d_ref, fem_ref)} dofs)")
print(f"  coarse start   : {err(E0):.3e}   ({n0} dofs)")
print(f"  h-adaptivity   : {err(E_h):.3e}   ({n_h} dofs)   {err(E0) / err(E_h):5.1f}x lower")
print(f"  r-adaptivity   : {err(E_r):.3e}   ({n_r} dofs, +0)   {err(E0) / err(E_r):5.1f}x   <- not r's regime here")
print(f"  p-adaptivity   : {err(E_p):.3e}   ({n_p} dofs, {mask_p.mean():.0%} enriched)   {err(E0) / err(E_p):5.1f}x lower")
print(
    f"  h then p       : {err(E_hp):.3e}   ({n_hp} dofs, {mask_hp.mean():.0%} enriched)   {err(E0) / err(E_hp):5.1f}x lower"
)

# --- where did each method spend its DOFs? ---------------------------------------------------------------
# The two features want different treatment, so this is the interesting half of the comparison.
cx, cy, R, _A = BUMP
near_c, in_b = _near(pts_p, CORNER, 0.22), _near(pts_p, (cx, cy), R)
print("\np put its covers:")
print(
    f"  corner {mask_p[near_c].mean():5.0%} | bump {mask_p[in_b].mean():5.0%} | elsewhere {mask_p[~(near_c | in_b)].mean():5.0%}"
)

cent = pts_h[tris_h].mean(axis=1)
e1, e2 = pts_h[tris_h][:, 1] - pts_h[tris_h][:, 0], pts_h[tris_h][:, 2] - pts_h[tris_h][:, 0]
hsz = np.sqrt(0.5 * np.abs(e1[:, 0] * e2[:, 1] - e1[:, 1] * e2[:, 0]))
nc, ib = _near(cent, CORNER, 0.22), _near(cent, (cx, cy), R)
print("h put its elements (mean cell size, smaller = more refined):")
print(f"  corner {hsz[nc].mean():.4f} | bump {hsz[ib].mean():.4f} | elsewhere {hsz[~(nc | ib)].mean():.4f}")

assert err(E0) > 0 and err(E_p) > 0, "the reference must bound every run -- otherwise it is not a reference"
assert err(E_h) < 0.25 * err(E0), "h-adaptivity should sharply cut the energy error"
assert err(E_p) < err(E_h), "p should beat h per DOF on this smooth-source problem"
assert err(E_hp) < err(E_p), "h then p should beat either alone"
assert 0.0 < mask_p.mean() < 1.0, "the p run enriched everything or nothing -- nothing was chosen"
assert _min_detj(pts_r, tris_r) > 0, "r-adaptivity must not tangle the mesh"
# Not tangling is a low bar -- `objective="energy"` on this source-driven problem passed it while
# flattening the smallest angle from 40.8 degrees to 3.2 (it lowers the integral of |grad u|^2 by
# squashing elements). Assert the mesh stays USABLE, not merely non-inverted.
_ang0, _ang_r = _min_angle(pts0, tris0), _min_angle(pts_r, tris_r)
print(f"\nmesh quality (smallest angle):  start {_ang0:.1f} deg  ->  after relocation {_ang_r:.1f} deg")
assert _ang_r > 0.6 * _ang0, f"relocation wrecked the mesh: {_ang0:.1f} deg -> {_ang_r:.1f} deg"