Smoothing a quad mesh¶

SmoothSession (nekmeshpy.quadmesh.smoother) smooths a quad mesh, flat or a surface in 3-D, one reviewed round at a time. It moves vertices, never connectivity (a curved mesh stays curved: High order), and every round reports what it did and can be undone. The order of work, step by step, is the workflow chart.

Quick start¶

import numpy as np
from nekmeshpy import linemesh, quadmesh
from nekmeshpy.io.writer import quad_to_vtu
from nekmeshpy.quadmesh.smoother import SmoothSession, OnCurve, circle, with_points

# a disk with a perturbed interior (any linear or curved QuadMesh works)
mesh = quadmesh.ogrid(linemesh.circle(1.0, 24, element_tag="wall"), 6, 4, wall_tag="wall")
jitter = np.random.default_rng(3).normal(size=mesh.points.shape) * [0.03, 0.03, 0]
mesh = with_points(mesh, mesh.points + jitter)

s = SmoothSession(mesh, [OnCurve("wall", circle([0, 0, 0], 1.0))], out_dir="smooth_rounds")
proposal = s.propose()        # runs s.suggest() on a copy; nothing is accepted yet
print(proposal.report)        # every metric before -> after, PASS / FAIL
s.accept(proposal)            # s.undo() restores the previous state bit for bit
quad_to_vtu(s.current, "smoothed.vtu")

call	does
`s.propose(schedule=None)`	a round on a copy (default: `s.suggest()`)
`s.accept(p)` / `s.undo(n=1)`	keep a round (refused if it failed) / pop accepted rounds
`s.run(n, review="print")`	`n` rounds unattended, stopping at a failure; `review="auto"` also stops when nothing more is gained
`s.history()`	the trend over rounds

The session is shared with the hex smoother (nekmeshpy.core.smoother.session).

Constraints¶

constraint	what it does	typical use
`Frozen(points)`	never moves	a feature placed by hand
`Tangent(points)`	slides along the boundary, never inward; nothing projects it back	a boundary with no geometry, under Laplace or Winslow
`OnSurface(points, surface=None, slide=True)`	slides tangentially, snapped back every step; no `surface`: fitted to the input faces; `slide=False`: projected once, then held	a curved wall, a planar face
`OnCurve(points, curve=None, slide=True)`	the same along a curve; no `curve`: a spline through the input boundary	a rim, a box crease
`KeepLayer(seeds, rings=1)`	freezes `rings` quad-rings off the seeds	a graded boundary layer
`SlideColumns(columns, direction)`	rigid columns, one shared degree of freedom each; cells between rows stay rectangles; frozen whole if any member is	rows across a wall; found by `constrain.grid_columns(mesh, seeds, direction, within=mask, outside_depth=1)`

points: ids, a boolean mask, or an edge-tag name. Targets: plane, sphere, line, circle, or any Surface(project, normal) / Curve(project, tangent); gradients use the target’s normal, never the mesh’s.
Precedence: Frozen > OnCurve > OnSurface > Tangent > Free; a vertex claimed by two surfaces or curves is frozen (and listed).
weights=w scales each quad’s shape term (0: shape irrelevant); the inversion penalty still applies.

boundary= (in resolve(mesh, constraints, boundary=...) and SmoothSession) is a shorthand for every boundary vertex no constraint names:

`boundary`	states	an unclaimed boundary vertex	use for
`"frozen"` (default)	`Frozen`	never moves	a boundary placed as meshed
`"free"`	nothing	moves as proposed: Laplace pulls a convex boundary in, the optimizer pushes it out	a boundary whose position does not matter
`"tangent"`	`Tangent`	slides along the boundary, never inward	Laplace, Winslow; no known geometry
`"slide"`	`OnCurve`, fitted	slides on a spline through the input boundary (surface quads too)	no known geometry, any method
—	`OnCurve` / `OnSurface` with a target	slides on the exact target	known geometry: always preferred

Corners (a turn above corner_deg=30 that is also a jump against its neighbours’) are frozen in every mode but "free". "tangent" keeps the boundary’s area under every method.

Warning

The shape metric seeks squares, so an unconstrained graded boundary layer gets flattened. Freeze it with KeepLayer.

How each piece works: Constraints and projection.

Schedule¶

s.suggest() proposes 10 x [laplace(3) + optimize(20)] plus two local lifts (user_schedule()); when any quad starts inverted, optimize stages with a wide-to-narrow delta come first. A schedule is a list of (stage, knobs).

stage	what it does	main knobs
`laplace`	moves each vertex toward its neighbours	`passes` 3, shrink `w` 0.97
`optimize`	descends the shape metric plus an inversion penalty; with a `delta` continuation (`untangle_schedule()`, a preset of `optimize`, not a separate stage) it pulls folded vertices through	`iters` 20, `kappa`, `eps`, `delta` 0.01 (the preset: 1.0 → 0.01, `iters` 60 each)
`laplace_local`	the failing quads only, lifting the smallest edges	`passes` 5, `w`
`winslow` / `winslow_local`	elliptic smoothing (all quads / the failing ones): an alternative to `laplace` for concave outlines and thin curved blocks; not in the suggested schedules, add it by hand: `("winslow", {"sweeps": 100})`	`sweeps` 200, `omega` 1, `tol` 1e-6
`tmop`	MFEM’s target-matrix optimizer, planar meshes, any order; optional, needs PyMFEM: TMOP via the MFEM backend	`target` `"input"`, `designed`

Keep graded layers with KeepLayer under every stage, Winslow most of all (Winslow smoothing).
kappa (1000) / eps (0.05) set the inversion penalty: ("optimize", {"iters": 20, "kappa": 1e4}) can unstall an inverted corner for about 10 % more time. A knob to try, not a better default.
adjust_schedule(schedule, "more_aggressive" | "more_conservative" | "same") changes shrink, passes and iterations, never a constraint.
Which stage when, measured side by side: Which method when.

The round report¶

#	section	shows
1	constraints	vertex counts per kind, targets, tolerances
2	what ran	each stage and why it stopped
3	quality	every metric before → after, marked better / WORSE / same: inverted and overlapping quads first, scaled J (min, mean, median, worst 1/5/10 %), aspect ratio, corner angles, size ratio, node spacing
4	geometry	are the vertices still on their targets (`on_target_report`)
5	propagate to hex	“quad-only” for now
6	undo	what `undo(1)` returns to

Warnings do not block acceptance; failures do (off target, more inverted or more overlapping quads). With out_dir, round_NN/before.vtu / after.vtu carry sj_min, phi_max, constraint_kind, dist_to_target, displacement.

Checks before smoothing¶

check	catches	opt out
gap check (`quadmesh.check_connectivity`: gencon’s rule, relative `connectivity_tol=0.2`)	points that should be one vertex, an edge in more than two quads, a quad naming a vertex twice, a collapsed edge, a pinch vertex; refuses to start	`allow_gaps=True` (stated in every report)
normal check	quads disagreeing with the surface normal (folded, reversed, wrong surface); `normal_report.n_disagree`	`check_normals=False`
overlap check (`quadmesh.overlapping_quads`)	overlapping input quads; recounted every round	`check_overlaps=False`

Curved surfaces¶

Claim the vertices with OnSurface: they slide and are projected back every step, never leaving the surface; on_target_report gives the largest distance.
Normals come from the surface: normals=, else surface=, else the OnSurface targets, and only then the quads’ winding. s.normal_report says which.
A tessellation becomes a surface through surface_from_trimesh(points, tris): an approximation (chord error remains, sharp edges blur), with one sign per surface voted against the mesh winding.
Judge against the surface, not only squareness:
- angle defect (quadmesh.angle_defect): on the perturbed test surfaces its error falls from about 0.15 to 0.01-0.02 rms, against a floor of about 0.01;
- sagitta: the chord error between vertices;
- overlapping quads: a branched star can have every Jacobian positive;
- Gauss-Bonnet residual (quadmesh.gauss_bonnet_residual): a topology identity; catches bad edge counts and pinch vertices, not poor representation.

from nekmeshpy.quadmesh.smoother import OnSurface, sphere, surface_from_trimesh

surf = sphere([0, 0, 0], 1.0)                   # analytic: plane, sphere, Surface(...)
surf = surface_from_trimesh(tri_points, tris)   # or a tessellation
s = SmoothSession(mesh, [OnSurface(np.arange(mesh.n_points), surf)], surface=surf)
print(s.normal_report)                          # normals' source; quads that disagree

Warning

Always give the surface. Normals from the quads’ winding are poisoned by one folded input quad: on a perturbed ellipsoid, smoothing then took the angle-defect error from 0.157 to 0.231 instead of 0.017. This, not the metric, was behind an earlier finding that smoothing worsens curvature.

High order¶

A curved input stays curved: the native stages move the vertices of its straight-sided shadow and the curved nodes follow the vertex displacement; a tmop stage moves every node of a planar mesh (TMOP via the MFEM backend). SmoothSession(..., keep_curved=False) smooths the shadow alone (highorder.linearize).
s.rebuild_high_order(order=None) rebuilds the curved nodes of the accepted vertices from the input’s curved nodes: each keeps its offset from the straight element and is projected back onto its target, which a carried node is not; order=N lifts a linear mesh. It discards a tmop round’s node moves.
Its HighOrderReport samples the scaled Jacobian at order 8 and the target distance before/after projection.
Limits: radial projection is not exactly Gauss-Lobatto arc spacing; curved validity is reported, not guaranteed; quads only.

Limitations¶

A smoothed surface is not carried into a hex mesh (no propagation, no block-state certificates). Hexes have their own, experimental smoother: Smoothing a hex mesh (experimental).
The shape metric seeks squares and undoes designed grading; protect layers with KeepLayer. A target-matrix metric that keeps the grading exists for planar meshes, at any order, through the optional TMOP stage.
Overlapping or folded input is reported, not repaired. A fold, mostly from projecting onto a curved surface, can leave a persistent overlap that, on perturbed capsules, longer schedules, a second pass or the delta continuation first did not remove (mechanism partly understood; a volume-preserving Laplacian is future work).
Size does not follow curvature (the metric aims at uniform size; a size target is being evaluated). With few irregular vertices part of the curvature error belongs to the connectivity: angle defects sum to 2 pi times the Euler characteristic.
Repairing an inverted quad is not guaranteed; the round says when one remains.

Note

Before/after viewers for flat, curved and real meshes, and what carrying a smoothed surface into hexes does, are in the Deep dives. Tests: tests/smoother/quad/.