Smoothing a hex mesh (experimental)

nekmeshpy.hexmesh.smoother is the quad smoother one rung up: the same constraint vocabulary, the same element-shrink Laplace and shape optimizer, the same schedules, on the corners of a HexMesh. It moves vertices only (on a curved mesh the curved nodes follow); connectivity, element tags and face tags never change.

Warning

Experimental. It has been tested on small blocks only (a few hundred hexes) and is pure NumPy: expect it to be slow on a large mesh.

Quick start

import numpy as np
from nekmeshpy import hexmesh, linemesh, quadmesh, writer
from nekmeshpy.hexmesh import smoother as hs

mesh = hexmesh.extrude(quadmesh.ogrid(linemesh.circle(1.0, 16), 4, 2), 1.0, 3)
inner = np.setdiff1d(np.arange(mesh.n_points), hexmesh.boundary_points(mesh))
P = mesh.points.copy()
P[inner] += np.random.default_rng(0).normal(scale=0.03, size=(inner.size, 3))
mesh = hs.with_points(mesh, P)                # a perturbed interior

# walls slide on themselves, rims along the rims, corners stay
s = hs.SmoothSession(mesh, [], boundary="slide", out_dir="smooth_rounds")
p = s.propose()               # runs s.suggest() on a copy; nothing is accepted yet
print(p.report)               # every metric before -> after, PASS / FAIL
s.accept(p)                   # s.undo() restores the input bit for bit
writer.to_vtu(s.current, "smoothed.vtu")

On this 144-hex cylinder the round passes: the minimum scaled Jacobian goes from 0.14 to 0.86, no hex inverts, and every boundary vertex stays on its cap plane, its wall or its rim (about 6 s on a laptop, running serial).

The session is the quad smoother’s own (nekmeshpy.core.smoother.session, shared by both rungs): propose, accept, undo, run and history behave as in the quad guide.

What exists

piece	names
constraints	Frozen, Tangent, OnSurface, OnCurve, KeepLayer(seeds, rings), SlideColumns, frozen_region(mesh, tags); targets plane, sphere, surface_from_trimesh, line, circle, spline_curve, or any Surface / Curve
resolution	resolve(mesh, constraints, boundary="frozen" | "free" | "tangent" | "slide", corner_deg=30), snap, constrain_direction, on_target_report
Laplace	laplace(mesh, resolved, passes=3, w=0.97); laplace_local in a schedule
optimizer	optimize(mesh, resolved, iters=20, delta=0.01, kappa=1000, eps=0.05); optimize_local in a schedule
schedules	run_schedule, user_schedule, untangle_schedule, suggest_schedule; register_method(name, fn) adds a stage
TMOP (optional)	the tmop stage, MFEM’s target-matrix optimizer: keeps a designed layer’s size and aspect (TMOP via the MFEM backend)
session	SmoothSession(mesh, constraints, boundary=..., sample_order=, out_dir=): propose / accept / undo / run / history, shared with the quad rung
readings	quality, shape_stats, round_report

The constraint objects, Resolved, snap and constrain_direction are the quad smoother’s own. At this rung a string selection names a face tag, and KeepLayer counts hex rings.

The shape metric. At a hex corner with edge vectors a, b, c (right-handed), the optimizer descends phi = ((|a|² + |b|² + |c|²)/3)³ / det(a, b, c)²: Knupp’s mean-ratio shape measure cubed. It is 1 only at a cube corner, does not change under scaling or rotation, and is exactly the quad metric with the dimension raised from 2 to 3. Inversion has its own penalty, kappa * max(0, eps - s), on the signed scaled Jacobian s = det / (|a||b||c|). The barrier-free delta continuation (untangle_schedule) pulls folded hexes through, which the exact metric cannot. The gradient is analytic and checked against finite differences.

The boundary, in three layers

boundary="slide" holds every unclaimed boundary vertex on the boundary as it was meshed:

1. Patches. Boundary faces are grouped into patches: two faces sharing an edge are in one patch when their normals differ by at most corner_deg. A vertex inside a patch slides on it. A planar patch is a plane; a curved one is its own input faces, read as a triangle surface.

2. Edges. A vertex where two patches meet slides along their intersection: a line when straight, else a spline through the input vertices (a cylinder rim is a closed one).

3. Corners. A vertex where three or more patches meet, or where an edge turns by a jump (more than corner_deg and three times its neighbours’ turn), is frozen.

The boundary shorthand, for Laplace and the optimizer alike (Constraints and projection):

boundary	an unclaimed boundary vertex
"frozen" (default)	holds still
"free"	moves as proposed: Laplace shrinks a convex boundary
"tangent"	a face vertex moves in its patch’s tangent plane, an edge vertex along the edge’s tangent; corners frozen. No target: the volume is kept, the wall can drift (A free curved wall)
"slide"	the three layers above: fitted OnSurface per patch, fitted OnCurve per patch edge, corners frozen

In the Laplace pass a boundary face is never shrunk across its own normal, the 3-D form of the quad rule. Each pass is guarded: a pass that adds an inverted hex, or lowers the minimum scaled Jacobian by more than max_sj_drop, is rolled back.

Session

hs.SmoothSession is the shared session (nekmeshpy.core.smoother.session.SmoothSession) bound to the hex rung: it holds the mesh, its constraints and a stack of accepted states. propose(schedule) runs on a copy and returns a Proposal (candidate mesh, report, schedule); accept pushes it (refused if the round failed or is stale), undo(n) pops, and run(n, review="print" | "auto") is the no-wait path. With out_dir= each round writes round_NN/before.vtu and after.vtu with per-hex sj_min and per-point displacement (for a curved mesh, per-hex sj_min, sj_corner_min and node_move_max). It takes a linear or a curved mesh and keeps its order: on a curved one the native stages carry the curved nodes and a tmop stage moves them all (keep_curved=False smooths the straight-sided shadow instead).

sess = hs.SmoothSession(mesh, [], boundary="frozen", sample_order=2, out_dir="rounds")
prop = sess.propose(hs.user_schedule(outer=1, n_opt=10))   # runs on a copy
print(prop.report)                                         # failures would block accept
if prop.report.passed:
    sess.accept(prop)
print(sess.history())
sess.undo()                                                # back to the input, bit for bit

Any stage registered with register_method (a third-party backend included) runs through the session by name.

The round report

round_report(before, after, resolved, stages=logs) (or prop.report of a session round) has the quad report’s sections except 5, “propagate”, which does not arise at the top rung. It gives:

• the constraints and the stages that ran;

• every corner metric before → after;

• one row per region (element_tags), with counts below THRESHOLDS = (0.1, 0.2, 0.3);

• the sampled reading: the .re2 trilinear map read on the GLL lattice of sample_order (default SCAN_ORDER, the solver’s); for a curved mesh also its curved map at sample_order and the corners alone, before → after;

• the distance of every target’s vertices (planes that held show 0);

• whether every frozen vertex (kept layers, frozen regions, corners) is bit-identical.

A moved frozen vertex, a vertex off its target or more inverted hexes fail the round.

Limits

• Curved meshes: only TMOP smooths them as curved. A tmop stage moves every node and keeps the order; the native stages move the vertices and carry the curved nodes along, without projecting them back onto a curved wall (TMOP via the MFEM backend, demo). There is no native curved optimization and no hex rebuild_high_order yet.

• Slow and lightly tested. Small blocks only so far; larger meshes will be tested step by step.

• A curved wall is approximated by its input faces: a vertex sliding on it lands on the faceted surface, so pass your own Surface for an exact wall.

• The metric seeks cubes and flattens designed grading: freeze layers with KeepLayer or frozen_region, or run the optional TMOP stage with target="input", which keeps them (demo).

• The session has no input gate: a hex mesh is conforming by construction.