Demo: a bent pipe at order 3

A pipe bent through 90 degrees at order 3: an O-grid section of 20 quads swept through 8 layers, 160 hexes, every wall node on the exact torus. The interior nodes, vertices and curved nodes alike, are perturbed; the wall is held and the caps slide in their planes (boundary="slide"). Three repairs: TMOP on the straight-sided shadow (the order-1 route), the native optimizer on the curved mesh, and TMOP on the curved mesh, both TMOP runs with target="input". The stage and its knobs: TMOP via the MFEM backend.

Minimum scaled Jacobian read three ways: at the stored nodes, on the order-7 lattice, at the corners

Key points

• TMOP smooths the curved mesh as curved. The order stays 3 and every node moves: the minimum scaled Jacobian rises from 0.24 to 0.57 on the order-7 lattice and from 0.68 to 0.85 at the corners, with no hex inverted, in 9 Newton iterations.

• The wall is held, the caps slide. No wall node moves, and the cap nodes stay in their planes to 2e-16. Between its nodes the wall stays as close to the torus as the input’s order-3 wall: 0.020 % of R.

• The order-1 route loses the wall. Its result is straight-sided: off the torus by up to 9.3 % of R, and 10.5 % short in volume. Its scaled Jacobian reads highest, but it is the reading of a different, faceted geometry.

• Moving vertices cannot repair curved nodes. The native optimizer raises the corners from 0.68 to 0.81, but the perturbation of the curved nodes rides along: 0.47 on the order-7 lattice against 0.57 for TMOP, and a mean of 0.79 at the nodes against 0.90. Its wall vertices slide on the faceted input wall and the carried nodes are not projected back, so the wall drifts to 0.022 % of R.

• Every result is conformal and watertight.

How it is set up

The code below is the demo script’s own; every name is the landed API.

1. Build a curved mesh. order= is a factory argument: the wall circle takes it, the O-grid carries the curvature inward, and sweep evaluates the bend’s path at every node, so every wall node lies on the torus.

def path(t):
    return np.stack([R_BEND * (1 - np.cos(np.pi * t / 2)), 0 * t, R_BEND * np.sin(np.pi * t / 2)], axis=-1)


def tangent(t):
    return np.stack([np.sin(np.pi * t / 2), 0 * t, np.cos(np.pi * t / 2)], axis=-1)


wall = linemesh.circle(R, 4 * N_SIDE, element_tag="wall", order=ORDER)    # order=: the wall nodes on the circle
section = quadmesh.ogrid(wall, N_SIDE, RADIAL)                            # the fill carries the curvature inward
exact = hexmesh.sweep(section, path, np.linspace(0.0, 1.0, LAYERS + 1), origin=(0.0, 0.0, 0.0), tangent=tangent,
                      first_tag="inlet", last_tag="outlet", element_tags="fluid")   # the path on every node

A mesh built from explicit points (linemesh.loft, from_grid, a plain loft between profiles) is high order in storage and straight in geometry: smoothing it as curved gains nothing.

2. Constrain the boundary. boundary="slide" puts the caps on their axis-aligned planes and the wall on a curved target. Under TMOP each edge and face node takes the hold of its vertices: cap nodes slide in their planes, wall nodes are held, interior nodes are free. Pass Frozen constraints to fix more.

r = hs.resolve(mesh, [], boundary="slide")      # caps slide in their planes, the curved wall is held

3. Run TMOP on the curved mesh, as a schedule stage (TMOP = {"target": "input"} in the script). The knobs that matter at order N are quad_order (default 2N + 2) and sample_order (default max(N, 7)), both in the stage’s dict.

out, logs = hs.run_schedule(mesh, r, [("tmop", TMOP)])   # every node moves; out.order == mesh.order
log = logs[-1]     # log.reason: iterations, energy, what was held; log.after: min det, scaled J at order 7

Or as a reviewed round. The session keeps a curved mesh curved; keep_curved=False smooths its straight-sided shadow and returns a linear mesh.

s = hs.SmoothSession(mesh, [], boundary="slide")          # keep_curved=True by default
p = s.propose([("tmop", {"target": "input"})])
print(p.report)               # adds the curved map at sample_order and the .re2 corners, before -> after
s.accept(p)

4. Read the result at the solver’s order, and at the corners. A reading at the stored nodes is silent between them, and a .re2 file keeps the 8 corners of each hex only.

q7 = hexmesh.quality_summary(mesh, order=SAMPLE_ORDER)    # the curved map, sampled at the solver's order
qc = hexmesh.corner_summary(mesh)                          # the 8 corners alone: what .re2 writes
print("min scaled J: %.3f at order %d, %.3f at the corners" % (q7.min, SAMPLE_ORDER, qc.min))

5. Export. .vtu keeps the curved geometry for viewing; .re2 writes the corners.

for k, m in meshes.items():
    writer.to_vtu(m, "02_bent_pipe_%s.vtu" % k)        # VTK Lagrange cells: the curved geometry, for viewing
writer.to_re2(mesh, "02_bent_pipe.re2", bc=VEL_BC)     # Nek5000's mesh file: the 8 corners of each hex only

The numbers

One scaled Jacobian, read three ways: at the stored nodes, on the order-7 lattice (the order the solver runs at), and through the 8 corners alone (the straight-sided hex a .re2 file keeps). For the order-1 route the stored nodes are the corners. The wall is read between its nodes too, on a finer lattice, against the exact torus.

variant	order	min SJ, nodes	mean SJ, nodes	min SJ, order 7	inverted, order 7	min SJ, corners	wall off the torus, max (% of R)	time (s)
input	3	0.3273	0.6866	0.2429	0	0.6757	0.020	—
order-1 route	1	0.8507	0.8767	0.8507	0	0.8507	9.324	0.19
native optimizer	3	0.4779	0.7932	0.4699	0	0.8064	0.022	3.30
TMOP, curved	3	0.5736	0.8957	0.5736	0	0.8488	0.020	25.10

Bold: the best value in each column among the results (highest scaled Jacobian; lowest count, objective, time, change or distance); the input row is for reference.

Times: a laptop, running serial.

Run PYTHONPATH=. python examples/smoother/mfem/02_bent_pipe.py (needs pip install nekmeshpy[smooth]); it prints these numbers.