core

core.mesh

Generic shared-point unstructured mesh – the gmsh/meshio-style data model.

class nekmeshpy.core.mesh.Mesh

Bases: object

Shared-point mesh model: points (P,3) plus a {cell_type: connectivity} dict, named point_sets and cell_sets, and gmsh field_data. Serves as the meshio bridge.

__init__(points, cells=None, point_sets=None, cell_sets=None, field_data=None)

Parameters:

• points (NDArray[float64])

• cells (dict[str, NDArray[int64]] | None)

• point_sets (dict[str, NDArray[int64]] | None)

• cell_sets (dict[str, dict[str, NDArray[int64]]] | None)

• field_data (dict[str, Any] | None)

Return type:

None

property n_points: int

Number of points.

property cell_types: list[str]

The cell-type keys present.

n_cells(ctype=None)

Number of cells of type ctype, or the total when ctype is None.

Parameters:

ctype (str | None)

Return type:

int

to_meshio()

Return an equivalent meshio.Mesh (requires meshio).

Return type:

Any

classmethod from_meshio(m)

Build a Mesh from a meshio.Mesh.

Parameters:

m (Any)

Return type:

Mesh

write(path, file_format=None)

Write via meshio (.vtu, .msh, .xdmf, …).

Parameters:

• path (str)

• file_format (str | None)

Return type:

str

classmethod read(path, file_format=None)

Read any meshio-supported file into a Mesh.

Parameters:

• path (str)

• file_format (str | None)

Return type:

Mesh

core.boundary_condition

Named boundary conditions: a face tag name and its Nek5000 BC code.

class nekmeshpy.core.boundary_condition.BoundaryConditions

Bases: object

Boundary names to Nek BC codes, one table per solved field.

velocity is the momentum field and is required. temperature is the second field Nek’s .re2 carries once a mesh has a solid region (nelgv < nelgt); it covers every element, and a name left out of it stays conformal ('E  '). Each table maps a name to a BCSpec.

Each entry carries a cbc and a boundaryID (see Side). Within one table each is all-or-none: state it for every boundary or for none, and a table that states none takes the defaults. A table’s boundaryID set is its own – velocity’s and temperature’s are separate – and runs 1..k without a gap, though several boundaries may share one. .vtu paints it as bc_id; a table that states none is painted by tag_of() instead, so its boundaries still differ.

There are deliberately no presets here. A name-to-code table is a statement about one piece of geometry – which opening is the inlet, which surface is a measurement plane – so it belongs in the mesher that knows, next to the tags it names, where a reader meets it at the same time as the mesh. A built-in nek_default() put the carotid’s own vocabulary in the toolkit and let two other examples inherit a mapping neither of them stated.

Periodicity is not stated here: a 'P  ' code needs a correspondence between two named face groups, which periodic= supplies to the writer.

__init__(velocity, temperature=None)

Parameters:

• velocity (Mapping[str, str | int | tuple[str | None, int | None] | Mapping[str, str | int | tuple[str | None, int | None] | None]])

• temperature (Mapping[str, str | int | tuple[str | None, int | None] | Mapping[str, str | int | tuple[str | None, int | None] | None]] | None)

Return type:

None

property fields: tuple[str, ...]

The fields stated, velocity first.

tag_of(field, name)

The 1-based position of name in field’s own table: the id a viewer paints when the table states no boundaryID.

Parameters:

• field (str)

• name (str)

Return type:

int

numbered(field)

Whether field’s table states boundaryID (for every boundary, since it is all-or-none).

Parameters:

field (str)

Return type:

bool

has(field, name)

Whether field’s table states name.

Parameters:

• field (str)

• name (str)

Return type:

bool

side(field, name, region)

(cbc, boundaryID) field writes for a row of boundary name owned by an element in region; None means “write no row from this side”.

A name with one Side reads the same from everywhere. One with per-region entries must name every region it actually borders, since a missing key is far more likely a typo than an intent to drop the face.

Parameters:

• field (str)

• name (str)

• region (str)

Return type:

tuple[str, int] | None

names_with(field, cbc)

Names field gives cbc from at least one side.

Parameters:

• field (str)

• cbc (str)

Return type:

set[str]

core.topology

Mesh topology / validity checks, decoupled from the mesh containers.

class nekmeshpy.core.topology.TopologyReport

Bases: NamedTuple

Facet inventory and validity verdicts for an all-hex mesh.

n_elements: int

Number of hex elements.

n_faces: int

Distinct quad facets, counted orientation-free.

n_boundary_faces: int

Facets borne by exactly one hex – the domain boundary.

n_internal_faces: int

Facets shared by exactly two hexes.

n_nonmanifold_faces: int

Facets shared by three or more hexes – always a defect.

n_open_edges: int

Boundary-surface edges not bounded by exactly two boundary faces (leaks).

n_hanging_points: int

Boundary points lying strictly inside another boundary edge (T-junctions).

n_components: int

Connected components of the hex-adjacency graph.

watertight: bool

True when the boundary is a closed, non-empty 2-manifold. Says nothing about conformity: a T-junction interface still reads watertight.

conformal: bool

True when there are no hanging points.

static __new__(_cls, n_elements, n_faces, n_boundary_faces, n_internal_faces, n_nonmanifold_faces, n_open_edges, n_hanging_points, n_components, watertight, conformal)

Create new instance of TopologyReport(n_elements, n_faces, n_boundary_faces, n_internal_faces, n_nonmanifold_faces, n_open_edges, n_hanging_points, n_components, watertight, conformal)

Parameters:

• n_elements (int)

• n_faces (int)

• n_boundary_faces (int)

• n_internal_faces (int)

• n_nonmanifold_faces (int)

• n_open_edges (int)

• n_hanging_points (int)

• n_components (int)

• watertight (bool)

• conformal (bool)

nekmeshpy.core.topology.hex_report(points, hexes)

Topology / watertightness report for an all-hex mesh.

Parameters:

• points (NDArray[float64])

• hexes (NDArray[int64])

Return type:

TopologyReport

nekmeshpy.core.topology.is_watertight(points, hexes)

True if the all-hex mesh’s boundary is a closed 2-manifold.

Parameters:

• points (NDArray[float64])

• hexes (NDArray[int64])

Return type:

bool

nekmeshpy.core.topology.count_overlapping_pairs(points, hexes)

Number of hex pairs whose volumes geometrically overlap (see _count_overlapping_pairs) – independent of, and not part of, hex_report(): unlike watertight/conformal (a fixed-size facet-incidence scan), this is a geometric broad-then-narrow-phase search whose candidate count can run into the hundreds of thousands on a large mesh, so it is not folded into the fast report every is_watertight()/hexmesh.is_conforming call already pays for – call it explicitly where the extra work is wanted (a summary, typically, alongside those two).

Parameters:

• points (NDArray[float64])

• hexes (NDArray[int64])

Return type:

int

nekmeshpy.core.topology.is_overlap_free(points, hexes)

True if no two hexes geometrically overlap (see count_overlapping_pairs()).

Parameters:

• points (NDArray[float64])

• hexes (NDArray[int64])

Return type:

bool

nekmeshpy.core.topology.quad_report(points, quads)

Topology / watertightness report for an all-quad mesh – hex_report() one rung down: the facets are edges, the boundary is a curve, and the boundary’s own “open edge” is an open point (a boundary point not shared by exactly two boundary edges), which is what TopologyReport.n_open_edges counts here.

Parameters:

• points (NDArray[float64])

• quads (NDArray[int64])

Return type:

TopologyReport

nekmeshpy.core.topology.overlapping_quads(points, quads)

Sorted ids of the quads that take part in at least one overlapping pair (see count_overlapping_quad_pairs()): the which behind its count.

Parameters:

• points (NDArray[float64])

• quads (NDArray[int64])

Return type:

NDArray[int64]

nekmeshpy.core.topology.count_overlapping_quad_pairs(points, quads)

Number of quad pairs whose areas geometrically overlap – independent of, and not part of, quad_report(), for the same cost reason as count_overlapping_pairs().

Parameters:

• points (NDArray[float64])

• quads (NDArray[int64])

Return type:

int

nekmeshpy.core.topology.is_overlap_free_quads(points, quads)

True if no two quads geometrically overlap (see count_overlapping_quad_pairs()).

Parameters:

• points (NDArray[float64])

• quads (NDArray[int64])

Return type:

bool

class nekmeshpy.core.topology.TagReport

Bases: NamedTuple

How a mesh’s side-tag table lines up with its topological boundary – the same report at the hex rung (faces) and the quad rung (edges).

The two are independent by design – boundary is what connectivity says, a side-tag table is a named subset that may also name interior entities – so they can disagree in exactly two ways, and both are worth seeing. An n_untagged_boundary above zero means the export has open sides no boundary condition covers; an n_tagged_interior above zero means it will write boundary conditions onto sides that are not on the boundary at all. Either can be deliberate (a flux-measurement plane is a tagged interior face; so is a conjugate fluid/solid interface that keeps the fluid’s wall condition), so these are counts to recognise, not assertions.

n_rows: int

Named sides. One tag per shared entity, so this counts entities, not rows – the two could differ only while a tag was addressed by (element, side).

n_untagged_boundary: int

Sides on the topological boundary that carry no name.

n_tagged_interior: int

Named sides that are not on the topological boundary.

static __new__(_cls, n_rows, n_untagged_boundary, n_tagged_interior)

Create new instance of TagReport(n_rows, n_untagged_boundary, n_tagged_interior)

Parameters:

• n_rows (int)

• n_untagged_boundary (int)

• n_tagged_interior (int)

nekmeshpy.core.topology.surface_report(points, tris)

Topology report for a triangle surface mesh.

Parameters:

• points (NDArray[float64])

• tris (NDArray[int64])

Return type:

dict[str, Any]

nekmeshpy.core.topology.is_closed(points, tris)

True if the surface is a closed 2-manifold.

Parameters:

• points (NDArray[float64])

• tris (NDArray[int64])

Return type:

bool

nekmeshpy.core.topology.format_report(report, kind='hex')

Human-readable multi-line summary of a hex, quad or surface report.

kind ("hex" or "quad") only changes the vocabulary: a quad mesh’s facets are edges, and its boundary’s open edge is an open point.

Parameters:

• report (TopologyReport | dict[str, Any])

• kind (str)

Return type:

str

core.quality

The shared schema of the per-rung scaled-Jacobian quality summaries.

nekmeshpy.core.quality.POOR_THRESHOLD = 0.2

Scaled Jacobian below which an element is counted as poor. Part of the public schema: it names both the n_poor field and the poor (<...) report line.

nekmeshpy.core.quality.SCAN_ORDER = 7

The order a report double-checks the mesh at. This is the polynomial order the solver will run it at, which is the only reading that predicts anything, so there is no universal right value – change it to match your solver.

nekmeshpy.core.quality.SCAN_BUDGET = 15000000

Sampled points a default order_scan may spend, summed over its orders. The work is n_elements * (order+1)**dim, and at SCAN_ORDER that is 512 points a hex – so a mesh too big for the cap is reported as not checked rather than quietly costing minutes inside a routine report. Raise it (or call order_scan directly) when you want the answer anyway.

class nekmeshpy.core.quality.OrderScan

Bases: NamedTuple

What the stored map looks like read on lattices finer than the mesh’s own.

An element’s Jacobian is sampled, exact at its (order+1)**dim GLL nodes and silent between them, so a clean QualitySummary at the mesh’s own order is not a certificate. This is the same map at more points – what a solver running at its own polynomial order computes.

orders: tuple[int, ...]

Sampling orders actually read, ascending.

min_sj: tuple[float, ...]

Smallest scaled Jacobian at each of orders.

n_inverted: tuple[int, ...]

Elements with a non-positive scaled Jacobian at each of orders.

skipped: tuple[int, ...]

Orders the budget refused, ascending – unchecked, not clean.

property clean: bool

True when at least one order was read and none found an inverted element.

A scan that read nothing is not clean – it is unchecked, which is the same reason skipped exists rather than the budget dropping orders in silence.

property worst: tuple[int, float]

(order, min_sj) of the harshest reading.

static __new__(_cls, orders, min_sj, n_inverted, skipped)

Create new instance of OrderScan(orders, min_sj, n_inverted, skipped)

Parameters:

• orders (tuple[int, ...])

• min_sj (tuple[float, ...])

• n_inverted (tuple[int, ...])

• skipped (tuple[int, ...])

class nekmeshpy.core.quality.SolverMetrics

Bases: NamedTuple

The three numbers Nek5000/nekRS print as mesh metrics:, of the trilinear (.re2) map read at one polynomial order.

:attr:`jac_ratio` is not a scaled Jacobian in :class:`QualitySummary`’s sense. nekRS’s line reads scaled Jacobian but computes min(J)/max(J) of the Jacobian determinant over an element’s GLL nodes – how uniform the map is, blind to skew, where the angle-based metric is blind to a size gradient. The two agree on a cube and on nothing else, so 0.08 here and 0.22 there are the same mesh.

order: int

Polynomial order the map was sampled at.

spacing_min: float

Shortest / longest distance between neighbouring GLL nodes, over the mesh.

spacing_max: float

Alias for field number 2

jac_ratio_min: float

min(J)/max(J) per element – smallest, largest, mean.

jac_ratio_max: float

Alias for field number 4

jac_ratio_mean: float

Alias for field number 5

aspect_min: float

Longest over shortest of an element’s 12 corner-to-corner edges – smallest, largest, mean.

aspect_max: float

Alias for field number 7

aspect_mean: float

Alias for field number 8

static __new__(_cls, order, spacing_min, spacing_max, jac_ratio_min, jac_ratio_max, jac_ratio_mean, aspect_min, aspect_max, aspect_mean)

Create new instance of SolverMetrics(order, spacing_min, spacing_max, jac_ratio_min, jac_ratio_max, jac_ratio_mean, aspect_min, aspect_max, aspect_mean)

Parameters:

• order (int)

• spacing_min (float)

• spacing_max (float)

• jac_ratio_min (float)

• jac_ratio_max (float)

• jac_ratio_mean (float)

• aspect_min (float)

• aspect_max (float)

• aspect_mean (float)

class nekmeshpy.core.quality.QualitySummary

Bases: NamedTuple

Aggregate scaled-Jacobian statistics over one mesh’s elements.

n_elements: int

Number of elements the statistics cover.

min: float

Smallest per-element scaled Jacobian; <= 0 means inverted / degenerate.

max: float

Largest per-element scaled Jacobian; 1 is a perfect corner.

mean: float

Mean over the elements.

median: float

Median over the elements.

n_inverted: int

Elements with a non-positive scaled Jacobian (inverted or degenerate).

n_poor: int

Elements below POOR_THRESHOLD – includes n_inverted.

static __new__(_cls, n_elements, min, max, mean, median, n_inverted, n_poor)

Create new instance of QualitySummary(n_elements, min, max, mean, median, n_inverted, n_poor)

Parameters:

• n_elements (int)

• min (float)

• max (float)

• mean (float)

• median (float)

• n_inverted (int)

• n_poor (int)

core.measure

Metric measures over an element’s node block: extent, size, centroid.

The rung-agnostic kernel behind linemesh.length / quadmesh.area / hexmesh.volume and the bounds / centroid each rung spells the same way. Everything here takes a node block – (E, (order+1)**dim, 3) in the lexicographic lattice order the element_blocks queries return – so the same quadrature serves the corner (linear) and the curved order-N reading of a mesh, and nothing in this module knows a container.

class nekmeshpy.core.measure.Bounds

Bases: NamedTuple

The axis-aligned bounding box of a point set – min / max corners, with the derived readings (size, center, diagonal) that callers actually ask for.

min: NDArray[float64]

Per-axis minimum (3,).

max: NDArray[float64]

Per-axis maximum (3,).

property size: NDArray[float64]

Per-axis extent max - min (3,).

property center: NDArray[float64]

The box centre – the midpoint of the extremes, not a centroid.

property diagonal: float

Length of the box diagonal – the scalar “how big is this” reading.

static __new__(_cls, min, max)

Create new instance of Bounds(min, max)

Parameters:

• min (Point)

• max (Point)

nekmeshpy.core.measure.bounds_of(points)

The bounding box of points (any leading shape, trailing axis the 3 components). Raises on an empty set: a box around nothing has no corners.

Parameters:

points (NDArray[float64])

Return type:

Bounds

nekmeshpy.core.measure.centroid_of(blocks, dim)

The measure-weighted centroid integral x d(Omega) / integral d(Omega) of the elements blocks describes – the mass-property centroid, not the mean of the nodes.

Parameters:

• blocks (NDArray[float64])

• dim (int)

Return type:

NDArray[float64]

nekmeshpy.core.measure.corner_blocks(points, conn, dim)

(E, 2**dim, 3) node blocks of the straight-sided elements conn names – the order-1 block of a mesh read at its corners only, whatever order it stores.

Parameters:

• points (NDArray[float64])

• conn (NDArray[int64])

• dim (int)

Return type:

NDArray[float64]

nekmeshpy.core.measure.integrate(blocks, dim, *, moments=False)

(measure (E,), moment (E,3)) – each element’s integral d(Omega) and, when moments is asked for, its integral x d(Omega) (zeros otherwise).

At dim=3 the measure is the signed volume integral of the isoparametric map, so an inverted hex contributes negatively and a wholly mis-oriented mesh reads negative – the same sign convention as the scaled Jacobian. A length (dim=1) or a surface area (dim=2) embedded in 3-D has no such sign and comes back non-negative.

Parameters:

• blocks (NDArray[float64])

• dim (int)

• moments (bool)

Return type:

tuple[NDArray[float64], NDArray[float64]]

core.stations

Where a sweep’s stations are, and what lays out across them.

Rung-agnostic by nature: a station is a parameter along the sweep, and a level is the copy of the slice sitting at one. Validating a sweep’s fractions, refining them onto the node lattice, sampling a path at them, and laying a slice’s entity table out across them are all that one subject, which belongs to no rung in particular – every loft / sweep / extrude on the ladder reads from here.

class nekmeshpy.core.stations.Sampled

Bases: Protocol

What an evaluated sweep needs of a slice, at any rung. The three containers share no base class, so this is structural – it names the reads, nothing more.

property order: int

property n_points: int

property points: NDArray[float64]

__init__(*args, **kwargs)

nekmeshpy.core.stations.at_levels(table, level, per_level)

table laid out at every level of a sweep – level * per_level + table, level-major, keeping table’s trailing shape: (n_level,) x (k, ...) -> (n_level * k, ...).

Every id a sweep writes is one of these: a section point / edge / quad id, shifted by how many of that entity one level owns. The loft bodies above pick the table and the levels; this is the layout they share.

Parameters:

• table (NDArray[int64])

• level (NDArray[int64])

• per_level (int)

Return type:

NDArray[int64]

nekmeshpy.core.stations.check_fraction_count(fr, *, loop, name)

Raise unless fr carries enough fractions for at least one sweep layer (two, with loop=True – the last fraction is the wrap back onto the first profile rather than a level of its own).

Parameters:

• fr (NDArray[float64])

• loop (bool)

• name (str)

Return type:

None

nekmeshpy.core.stations.refined_lattice(fractions, order)

The n*order + 1 parameter positions of every node of the order-order chain graded by fractions (n = len(fractions) - 1 elements): element i’s node a sits at fr[i] + g[a]*(fr[i+1] - fr[i]) for the GLL nodes g on [0, 1], and the chain ends at fr[-1].

Parameters:

• fractions (NDArray[float64])

• order (int)

Return type:

NDArray[float64]

nekmeshpy.core.stations.spline_levels(blocks, t, *, loop)

Resample a per-level node block onto the sweep lattice through a cubic spline in the level index, rather than a straight blend between the two bounding levels.

blocks is any node array stacked over the levels – (n_lev,) + rest – and the levels are taken to sit at 0, 1, ... n_lev-1, which is exactly where loft’s own layers put them. So evaluating at t returns level k verbatim at integer t = k: the spline interpolates, it does not approximate, and a sweep resampled this way still passes through every profile it was given.

This is what a plain loft cannot do. Its sweep-direction interior nodes are a straight blend of the two profiles bounding their layer, so the surface is linear between profiles however high the order is – the trap the module docs call “high order in storage, linear in geometry”. A spline sees the neighbouring profiles too, so a station’s high-order nodes bend the way the stack bends, and a feature that turns sharply across a few profiles is resolved instead of being cut across.

loop closes the spline periodically, so the seam layer is no less smooth than any other – an open spline would leave the one place a closed sweep cannot show a joint. Fewer than three levels carry no curvature to fit, and fall back to the straight blend.

Parameters:

• blocks (NDArray[float64])

• t (NDArray[float64])

• loop (bool)

Return type:

NDArray[float64]

nekmeshpy.core.stations.sweep_lattice(fractions, order, *, loop, name)

Validate a sweep’s fractions and return (fr, node lattice).

Parameters:

• fractions (NDArray[float64])

• order (int)

• loop (bool)

• name (str)

Return type:

tuple[NDArray[float64], NDArray[float64]]

nekmeshpy.core.stations.split_evaluated(profs, order, *, loop, conn, noun, elems, name)

Validate a sweep whose slices were evaluated on the refined node lattice rather than handed in, and split them into (slices, sweep_nodes).

Every rung above the line does this identically – one slice per lattice value, all index-paired with the first, the wrap level dropped on a loop – so the only things that vary are what to compare (conn, which doubles as the element count) and what to call the pieces in the errors (noun / elems). Those nouns are not decoration: they are what the rungs’ own tests match on.

The lattice itself is not an argument: the slices are the sweep, one per lattice value by construction at every caller, so passing the parameters alongside them only offered a second thing to disagree with. What remains is the one property the list must have on its own – a whole number of layers of order levels each – and positions are reported as level indices, which every caller shares, rather than as parameter values, which only an evaluated sweep has.

Parameters:

• profs (Sequence[M])

• order (int)

• loop (bool)

• conn (Callable[[M], NDArray[int64]])

• noun (str)

• elems (str)

• name (str)

Return type:

tuple[list[M], list[list[M]]]

nekmeshpy.core.stations.sweep_path(path, tangent, tv)

Sample a sweep’s centreline (and, if given, its analytic derivative) on the station parameters tv, as (K,3) arrays.

Parameters:

• path (Callable[[NDArray[float64]], NDArray[float64]])

• tangent (Callable[[NDArray[float64]], NDArray[float64]] | None)

• tv (NDArray[float64])

Return type:

tuple[NDArray[float64], NDArray[float64] | None]

core.affine

Affine maps behind the rung-preserving translate / rotate / scale.

nekmeshpy.core.affine.ORIGIN: NDArray[float64] = array([0., 0., 0.])

The world origin – default centre for rotation() / scaling().

nekmeshpy.core.affine.Z_AXIS: NDArray[float64] = array([0., 0., 1.])

The +z axis – default rotation axis, matching the sweep factories’ axis=.

nekmeshpy.core.affine.apply(P, matrix, offset)

Map every coordinate of P (any leading shape, trailing axis 3) through the affine pair, returning a new array of the same shape. With matrix=None this is the exact translation P + offset.

Parameters:

• P (NDArray[float64])

• matrix (NDArray[float64] | None)

• offset (NDArray[float64])

Return type:

NDArray[float64]

nekmeshpy.core.affine.translation(vector)

The affine map that shifts by vector (3,).

Parameters:

vector (NDArray[float64] | Sequence[float])

Return type:

tuple[NDArray[float64] | None, NDArray[float64]]

nekmeshpy.core.affine.scaling(factor, center=array([0., 0., 0.]))

The affine map that scales about center by factor – a scalar (uniform) or a (3,) per-axis vector. A zero or negative component is rejected: it would collapse or invert every element.

Parameters:

• factor (float | NDArray[float64] | Sequence[float])

• center (NDArray[float64] | Sequence[float])

Return type:

tuple[NDArray[float64] | None, NDArray[float64]]

nekmeshpy.core.affine.rotation(angle, axis=array([0., 0., 1.]), center=array([0., 0., 0.]))

The affine map that rotates by angle radians about the line through center with direction axis (right-handed; axis need not be normalized).

Parameters:

• angle (float)

• axis (NDArray[float64] | Sequence[float])

• center (NDArray[float64] | Sequence[float])

Return type:

tuple[NDArray[float64] | None, NDArray[float64]]

nekmeshpy.core.affine.reflection(normal, point=array([0., 0., 0.]))

The affine map that mirrors through the plane with unit normal (need not be normalized) passing through point – the Householder matrix I - 2 n n^T.

Its determinant is -1, so it turns every element inside out: this is the map mirror applies to the coordinates, paired with a re-winding of the connectivity. Handing it to transform alone leaves an inverted mesh.

Parameters:

• normal (NDArray[float64] | Sequence[float])

• point (NDArray[float64] | Sequence[float])

Return type:

tuple[NDArray[float64] | None, NDArray[float64]]

core.frames

Moving orthonormal frames along a sampled curve – the placement half of a sweep.

nekmeshpy.core.frames.SEGMENT_TOL: float = 1e-12

Segments shorter than this fraction of the curve’s overall extent are treated as repeated points. Relative, so the check is scale-free: a curve modelled in metres and the same curve in millimetres are rejected (or accepted) identically.

nekmeshpy.core.frames.PARALLEL_TOL: float = 1e-12

A tangent whose component along up exceeds this in magnitude leaves too little of up in the cross-section plane for fixed_up() to normalize meaningfully; at 1 - 1e-12 the surviving perpendicular component is ~1.4e-6 of up, already at the edge of what double precision resolves in a direction.

nekmeshpy.core.frames.PLANAR_TOL: float = 1e-09

A profile point this far off its own best-fit plane, relative to the profile’s overall extent, is not a planar cross-section. Relative, so a section modelled in metres and the same section in millimetres are judged identically; loose enough to pass a section assembled from several factories, tight enough that a genuinely bowed profile is caught rather than sheared.

nekmeshpy.core.frames.CURVATURE_TOL: float = 1e-10

Below this the unit-tangent turning per sample is indistinguishable from round-off and frenet() has no principal normal to speak of. Absolute, because the quantity it bounds – |t[i+1] - t[i-1]| / 2 – is already dimensionless.

nekmeshpy.core.frames.tangents(points, *, loop=False)

Unit tangents (K,3) of the (K,3) sampled curve points.

Parameters:

• points (NDArray[float64])

• loop (bool)

Return type:

NDArray[float64]

nekmeshpy.core.frames.fixed_up(tangents, up)

(K,3,3) frames from Gram-Schmidt of a constant world up against each unit tangent: u = normalize(up - (up.w) w), v = w x u, w = the tangent.

Parameters:

• tangents (NDArray[float64])

• up (NDArray[float64] | Sequence[float])

Return type:

NDArray[float64]

nekmeshpy.core.frames.parallel_transport(points, tangents, up0, *, loop=False, distribute=True)

(K,3,3) rotation-minimizing frames along points, seeded by up0.

Parameters:

• points (NDArray[float64])

• tangents (NDArray[float64])

• up0 (NDArray[float64] | Sequence[float])

• loop (bool)

• distribute (bool)

Return type:

NDArray[float64]

nekmeshpy.core.frames.spin(u, tangents, angles)

Rotate each cross-section vector u[k] by angles[k] about tangents[k], right-handed about the tangent. u must already be perpendicular to its tangent; the result is renormalized, not reprojected.

Parameters:

• u (NDArray[float64])

• tangents (NDArray[float64])

• angles (NDArray[float64])

Return type:

NDArray[float64]

nekmeshpy.core.frames.holonomy(points, tangents, up0)

The residual twist, in radians, of transporting a frame once around the closed curve points (the closing segment points[K-1] -> points[0] included).

Parameters:

• points (NDArray[float64])

• tangents (NDArray[float64])

• up0 (NDArray[float64] | Sequence[float])

Return type:

float

nekmeshpy.core.frames.frenet(points, tangents)

(K,3,3) Frenet-Serret frames: u the principal normal, v the binormal, w the tangent.

Parameters:

• points (NDArray[float64])

• tangents (NDArray[float64])

Return type:

NDArray[float64]

nekmeshpy.core.frames.frame_transform(R_from, origin_from, R_to, origin_to)

The single rigid (matrix, offset) carrying the frame (R_from, origin_from) onto (R_to, origin_to), in the convention of nekmeshpy.core.affine.apply() – p @ matrix.T + offset.

Parameters:

• R_from (NDArray[float64])

• origin_from (NDArray[float64])

• R_to (NDArray[float64])

• origin_to (NDArray[float64])

Return type:

tuple[NDArray[float64], NDArray[float64]]

nekmeshpy.core.frames.plane_frame(points, *, normal=None, origin=None, hint=None)

The local frame of a planar cross-section: the (R, origin) pair naming the plane the profile was authored in, ready to hand to frame_transform() as the _from half of a sweep.

Parameters:

• points (NDArray[float64])

• normal (NDArray[float64] | Sequence[float] | None)

• origin (NDArray[float64] | Sequence[float] | None)

• hint (NDArray[float64] | Sequence[float] | None)

Return type:

tuple[NDArray[float64], NDArray[float64]]

nekmeshpy.core.frames.sweep_placements(profile_points, path_points, *, orientation='transport', up=None, twist=0.0, close_twist=True, loop=False, origin=None, normal=None, path_tangents=None)

One rigid (matrix, offset) per path station, carrying a planar profile from its own plane onto the moving frame of the sampled curve path_points (K,3).

Parameters:

• profile_points (NDArray[float64])

• path_points (NDArray[float64])

• orientation (Literal['transport', 'fixed', 'frenet'])

• up (NDArray[float64] | Sequence[float] | None)

• twist (float)

• close_twist (bool)

• loop (bool)

• origin (NDArray[float64] | Sequence[float] | None)

• normal (NDArray[float64] | Sequence[float] | None)

• path_tangents (NDArray[float64] | None)

Return type:

list[tuple[NDArray[float64], NDArray[float64]]]

core.paths

A declarative 3-D turtle-walk path builder.

A path is a move table – line(), arc(), helix() – walked by walk() into a Path. The turtle carries a full orthonormal frame, not just a heading, which is what makes the vocabulary three-dimensional: tilt rolls a bend’s own plane about the direction of travel, so one arc verb reaches every direction out of the current frame, and roll spins the frame about that direction as the move travels. Every move starts at the previous one’s end point and keeps its frame, so a walk is C1 by construction – no fillet fitting, no corner rounding.

Each move is a constant screw, closed-form in arc length: exact total length (no quadrature), an exact analytic tangent (differencing a sampled path is only O(h**2), and worst exactly at a junction where the curvature jumps), and an s in [0, 1] that is true arc length, so an element length is a station count.

The frame the walk carried comes back on Path.up, and hexmesh.sweep_path holds it per station unless told otherwise – which is the only way an out-of-plane bend or a distributed twist can reach the mesh, since a sweep’s own generators fix the frame only up to how the path itself was rolled.

class nekmeshpy.core.paths.Arc

Bases: NamedTuple

A circular bend of radius through angle degrees, in the plane tilt degrees round from the frame’s own up.

radius: float

Radius of the bend.

angle: float

Signed turn in degrees – positive turns toward the frame’s left.

tilt: float

0 turns left, 180 right, 90 pitches up, anything between is oblique.

Type:

Degrees the bend plane is rolled about the heading

roll: float

Degrees of frame spin about the heading, spread over the whole bend – on top of the rotation the bend itself carries.

static __new__(_cls, radius, angle, tilt=0.0, roll=0.0)

Create new instance of Arc(radius, angle, tilt, roll)

Parameters:

• radius (float)

• angle (float)

• tilt (float)

• roll (float)

class nekmeshpy.core.paths.Helix

Bases: NamedTuple

A constant screw: a bend of radius through angle degrees while advancing rise per full turn along the bend axis.

radius: float

Radius measured from the helix axis.

angle: float

Signed turn in degrees about the helix axis.

rise: float

Advance along the helix axis per full 360-degree turn.

tilt: float

Degrees the helix axis is rolled about the heading, as Arc.tilt.

roll: float

Degrees of frame spin about the heading, spread over the whole run.

static __new__(_cls, radius, angle, rise=0.0, tilt=0.0, roll=0.0)

Create new instance of Helix(radius, angle, rise, tilt, roll)

Parameters:

• radius (float)

• angle (float)

• rise (float)

• tilt (float)

• roll (float)

nekmeshpy.core.paths.LENGTH_TOL = 1e-12

A move shorter than this carries no length to distribute a roll over, and would put two junctions on one s. Absolute rather than relative to the walk: the walk’s own extent is not known until every move has been read, and the quantity is a length that has to be divided by, so what matters is that it is above the reciprocal’s noise.

class nekmeshpy.core.paths.Line

Bases: NamedTuple

A straight run of length, spinning the frame roll degrees as it goes.

length: float

Arc length of the run.

roll: float

Degrees of frame spin about the heading, spread over the whole run.

static __new__(_cls, length, roll=0.0)

Create new instance of Line(length, roll)

Parameters:

• length (float)

• roll (float)

class nekmeshpy.core.paths.Path

Bases: NamedTuple

A curve in space, exposed as continuous callables of normalized arc length s in [0, 1] rather than as its own segment table – what hexmesh.lift.sweep_path consumes. walk() builds one; so can a caller with a parametrization of their own, which is why every field is a callable and only up may be left out.

centerline: Callable[[NDArray[float64]], NDArray[float64]]

(K,) in [0, 1] -> (K, 3) points.

tangent: Callable[[NDArray[float64]], NDArray[float64]]

(K,) in [0, 1] -> (K, 3) unit tangents, the analytic derivative of centerline. A direction: no origin enters it.

total_length: float

The exact total arc length.

break_fractions: NDArray[float64]

Normalized s of every junction between moves, strictly increasing.

up: Callable[[NDArray[float64]], NDArray[float64]] | None

(K,) in [0, 1] -> (K, 3) unit vectors perpendicular to tangent: the cross-section reference axis the path itself carries, so a sweep along it needs no orientation/up of its own. None on a path that names no frame, which leaves the sweep to build one.

static __new__(_cls, centerline, tangent, total_length, break_fractions, up=None)

Create new instance of Path(centerline, tangent, total_length, break_fractions, up)

Parameters:

• centerline (Callable[[FloatArray], PointArray])

• tangent (Callable[[FloatArray], PointArray])

• total_length (float)

• break_fractions (FloatArray)

• up (Callable[[FloatArray], PointArray] | None)

nekmeshpy.core.paths.arc(radius, angle, *, tilt=0.0, roll=0.0)

A bend of radius through angle degrees – tilt picks which way it bends (0 left, 180 right, 90 up), roll twists the section as it travels.

Parameters:

• radius (float)

• angle (float)

• tilt (float)

• roll (float)

Return type:

Arc

nekmeshpy.core.paths.helix(radius, angle, *, rise=0.0, tilt=0.0, roll=0.0)

A screw of radius through angle degrees, climbing rise per full turn along the axis tilt names.

Parameters:

• radius (float)

• angle (float)

• rise (float)

• tilt (float)

• roll (float)

Return type:

Helix

nekmeshpy.core.paths.line(length, *, roll=0.0)

A straight run of length, twisting roll degrees over its own length.

Parameters:

• length (float)

• roll (float)

Return type:

Line

nekmeshpy.core.paths.resolve_frame(path, orientation, up)

The (orientation, up) a sweep along path should run with.

A path built by walk() carries its own frame, and with no orientation asked for that frame is the answer – held per station, which is what "fixed" means. An explicit orientation (or an explicit up) overrides it, because a caller who names a generator is asking for that generator’s frame, not the path’s.

Parameters:

• path (Path)

• orientation (Literal['transport', 'fixed', 'frenet'] | None)

• up (NDArray[float64] | Sequence[float] | Callable[[NDArray[float64]], NDArray[float64]] | None)

Return type:

tuple[Literal[‘transport’, ‘fixed’, ‘frenet’], ~numpy._typing._array_like.NDArray[~numpy.float64] | ~collections.abc.Sequence[float] | ~typing.Callable[[~numpy._typing._array_like.NDArray[~numpy.float64]], ~numpy._typing._array_like.NDArray[~numpy.float64]] | None]

nekmeshpy.core.paths.sample_up(up, t)

Evaluate a callable up on the station parameters t, and pass anything else through untouched – the one place a path’s own frame becomes the (K,3) field frames.sweep_placements consumes.

Parameters:

• up (NDArray[float64] | Sequence[float] | Callable[[NDArray[float64]], NDArray[float64]] | None)

• t (NDArray[float64])

Return type:

NDArray[float64] | Sequence[float] | Callable[[NDArray[float64]], NDArray[float64]] | None

nekmeshpy.core.paths.walk(moves, *, start=(0.0, 0.0, 0.0), heading=(1.0, 0.0, 0.0), up=(0.0, 0.0, 1.0))

Walk moves in space from start, setting out along heading with the cross-section’s own up along up, into a Path that carries its own moving frame.

up is normalized and projected off heading – it names the frame’s phase, not a length, so only its component across the direction of travel is used, and only a seed parallel to heading is refused. From there the moves steer: a positive angle turns toward the walk’s left (up x heading), tilt rolls that bend plane about the heading (90 pitches up), and roll spins the frame right-handed about the heading over the move’s own length.

Parameters:

• moves (Sequence[Line | Arc | Helix])

• start (NDArray[float64] | Sequence[float])

• heading (NDArray[float64] | Sequence[float])

• up (NDArray[float64] | Sequence[float])

Return type:

Path

core.surfaces

Curves carried as their parametrization on a surface, rather than as points.

class nekmeshpy.core.surfaces.SurfaceCurve

Bases: NamedTuple

A curve living in a surface’s parameter domain.

g: Callable[[NDArray[float64]], NDArray[float64]]

(K,) curve parameters -> (K, 2) surface parameters. Vectorized.

fr: NDArray[float64]

The curve-parameter values of this curve’s own nodes, one per node. May run descending – that is how reverse() traverses a curve the other way, and what loft_fn accepts a descending sequence for.

static __new__(_cls, g, fr)

Create new instance of SurfaceCurve(g, fr)

Parameters:

• g (Callable[[FloatArray], FloatArray])

• fr (FloatArray)

nekmeshpy.core.surfaces.blend(a, b, lam)

The two curves interpolated in parameter space at lam – so the result is still exactly on the surface, which a point-space lerp would not be.

Parameters:

• a (SurfaceCurve)

• b (SurfaceCurve)

• lam (float)

Return type:

SurfaceCurve

nekmeshpy.core.surfaces.curve(g, fr)

A SurfaceCurve from a parametrization and its node parameters.

Parameters:

• g (Callable[[NDArray[float64]], NDArray[float64]])

• fr (NDArray[float64] | Sequence[float])

Return type:

SurfaceCurve

nekmeshpy.core.surfaces.node(c, i)

The (2,) surface parameters of node i.

Parameters:

• c (SurfaceCurve)

• i (int)

Return type:

NDArray[float64]

nekmeshpy.core.surfaces.reparam(c, pa, pb)

The straight pa -> pb parameter segment, but expressed on c’s own parameter domain and carrying c’s own nodes – so the two can be blend()-ed station by station.

Parameters:

• c (SurfaceCurve)

• pa (NDArray[float64] | Sequence[float])

• pb (NDArray[float64] | Sequence[float])

Return type:

SurfaceCurve

nekmeshpy.core.surfaces.reverse(c)

The same curve traversed the other way, by reversing its node parameters.

Parameters:

c (SurfaceCurve)

Return type:

SurfaceCurve

nekmeshpy.core.surfaces.ruled(pa, pb, n)

The straight segment from pa to pb in parameter space, with n elements evenly spaced in the parameter, on the domain [0, 1].

Parameters:

• pa (NDArray[float64] | Sequence[float])

• pb (NDArray[float64] | Sequence[float])

• n (int)

Return type:

SurfaceCurve

nekmeshpy.core.surfaces.segment(c, i0, i1)

The piece of c between nodes i0 and i1, as a parametrization on [0, 1] – ready for a Coons patch, which wants plain callables.

Parameters:

• c (SurfaceCurve)

• i0 (int)

• i1 (int)

Return type:

Callable[[NDArray[float64]], NDArray[float64]]

nekmeshpy.core.surfaces.shift(c, delta)

The same curve with a constant offset added in parameter space.

Parameters:

• c (SurfaceCurve)

• delta (NDArray[float64] | Sequence[float])

Return type:

SurfaceCurve

nekmeshpy.core.surfaces.spoke(start, end)

The straight parameter-space segment start -> end on [0, 1], as a plain callable – ruled() without a node table, for a Coons boundary.

Parameters:

• start (NDArray[float64] | Sequence[float])

• end (NDArray[float64] | Sequence[float])

Return type:

Callable[[NDArray[float64]], NDArray[float64]]

core.fields

Mesh-sizing fields (gmsh-style) and 1-D point distributions.

class nekmeshpy.core.fields.Field

Bases: object

Base sizing field: field(points) -> sizes mapping (k,3) points to k target edge sizes.

sample(point)

Target size at a single (3,) point.

Parameters:

point (NDArray[float64])

Return type:

float

class nekmeshpy.core.fields.ConstantField

Bases: Field

__init__(size)

Parameters:

size (float)

Return type:

None

class nekmeshpy.core.fields.AxisLinearField

Bases: Field

Size varies linearly along one axis from size0 at c0 to size1 at c1 (clamped outside).

__init__(axis, c0, size0, c1, size1)

Parameters:

• axis (int)

• c0 (float)

• size0 (float)

• c1 (float)

• size1 (float)

Return type:

None

class nekmeshpy.core.fields.DistanceField

Bases: Field

Size grows from size_near at the given points to size_far beyond dist_far (linear ramp on distance to the nearest source point).

__init__(points, size_near, dist_far, size_far)

Parameters:

• points (NDArray[float64])

• size_near (float)

• dist_far (float)

• size_far (float)

Return type:

None

class nekmeshpy.core.fields.MinField

Bases: Field

Pointwise minimum of several fields (the finest constraint wins).

__init__(*fields)

Parameters:

fields (Field)

Return type:

None

nekmeshpy.core.fields.geometric_spacing(n, ratio=1.0)

n+1 normalized positions in [0,1] with a geometric size ratio between consecutive cells (ratio==1 -> uniform).

Parameters:

• n (int)

• ratio (float)

Return type:

NDArray[float64]

nekmeshpy.core.fields.uniform_spacing(n)

n+1 uniformly spaced positions in [0, 1]; shorthand for geometric_spacing(n, 1.0).

Parameters:

n (int)

Return type:

NDArray[float64]

nekmeshpy.core.fields.symmetric_spacing(n, ratio=1.0)

n+1 positions in [0,1] clustered toward both ends: geometric spacing per half, mirrored (n must be even; ratio==1 -> uniform).

Parameters:

• n (int)

• ratio (float)

Return type:

NDArray[float64]

nekmeshpy.core.fields.validate_layers(positions)

Normalize a layer specification to the flattened array of normalized layer positions every sweep / fill factory consumes.

Parameters:

positions (int | NDArray[float64])

Return type:

NDArray[float64]

nekmeshpy.core.fields.gll_nodes(order)

The order+1 Gauss-Lobatto-Legendre nodes mapped to [0, 1], ascending.

Parameters:

order (int)

Return type:

NDArray[float64]

nekmeshpy.core.fields.gll_weights(order)

The order+1 GLL quadrature weights on [0, 1] (sum to 1), matching gll_nodes(). Used by the (deferred) order-N quality metric.

Parameters:

order (int)

Return type:

NDArray[float64]

nekmeshpy.core.fields.lagrange_matrix(nodes, eval_pts)

The (len(eval_pts), len(nodes)) matrix of Lagrange basis values: row i holds L_k(eval_pts[i]) for each node k. M @ f interpolates nodal values f to eval_pts; each row sums to 1 (partition of unity).

Parameters:

• nodes (NDArray[float64])

• eval_pts (NDArray[float64])

Return type:

NDArray[float64]

nekmeshpy.core.fields.lagrange_derivative_matrix(nodes, eval_pts)

The (len(eval_pts), len(nodes)) matrix of Lagrange basis derivatives: entry (i, k) is L'_k(eval_pts[i]).

Parameters:

• nodes (NDArray[float64])

• eval_pts (NDArray[float64])

Return type:

NDArray[float64]

nekmeshpy.core.fields.distribution_from_field(field, p0, p1, max_cells=200)

Graded positions along p0 -> p1 where each cell length approximates the field’s target size. Greedy walk stepping by the local size, then rescaled to land on p1. Returns normalized positions in [0,1].

Parameters:

• field (Field)

• p0 (NDArray[float64])

• p1 (NDArray[float64])

• max_cells (int)

Return type:

NDArray[float64]

core.selection

Resolving a select / remove argument to a mask over a mesh’s elements.

nekmeshpy.core.selection.Selection

one tag name, several tag names, a boolean mask over the elements, or element ids.

Type:

A select / remove argument

alias of str | Sequence[str] | NDArray[str_] | NDArray[bool] | NDArray[int64] | Sequence[int]

nekmeshpy.core.selection.mask_for_selection(which, tags)

The (n_elements,) boolean mask a select / remove argument names.

which is a tag name or a list / array of them (resolved by Tags.flatnonzero; a name absent from the mesh’s vocabulary is an error, since a silent empty selection is almost always a typo), a ready (n_elements,) boolean mask, or an array of element ids.

Parameters:

• which (str | Sequence[str] | NDArray[str_] | NDArray[bool] | NDArray[int64] | Sequence[int])

• tags (Tags)

Return type:

NDArray[bool]

core.smoother

The rung-agnostic smoother engine behind quadmesh.smoother and hexmesh.smoother.

core.smoother.constraints

Smoother constraints on points: what each vertex may do, and how the target geometry holds it.

Rung-agnostic: nothing here knows a quad from a hex. A smoother never guesses. Every vertex gets one mobility – FREE, ON_SURFACE, ON_CURVE, TANGENT or FROZEN – resolved once from the constraints the caller states (resolve, in each rung, which owns the connectivity that needs), and the two operations the smoothing methods share act on it:

• snap() projects constrained vertices back onto their target (x += w (P(x) - x); w = 1 is the full projection) and restores frozen ones;

• constrain_direction() removes from a gradient or a step whatever the target forbids – the component along the surface normal, everything but the tangent on a curve, all of it on a frozen vertex – using the target’s own normal, never the mesh’s.

Strict geometry mode (the default) is both, every step. Precedence is FROZEN > ON_CURVE > ON_SURFACE > TANGENT > FREE; a vertex claimed by two surfaces, or by two curves (a corner), is frozen and listed in the report – it is never given an order-dependent projection implicitly. Unclaimed boundary vertices are frozen by default, so a block’s seam cannot drift (boundary="free" opens them, "tangent" lets them move along the boundary only).

nekmeshpy.core.smoother.constraints.CORNER_DEG = 30.0

The turning angle, in degrees, above which a boundary vertex is a corner. It is frozen by resolve(..., boundary="tangent"), and frozen, ending its chains, by boundary="slide".

nekmeshpy.core.smoother.constraints.Select

A constraint names its vertices as point ids, a boolean mask, or an edge-tag name (the points of every edge carrying it – “this tagged loop lies on this curve”).

alias of str | NDArray[int64] | NDArray[bool]

class nekmeshpy.core.smoother.constraints.Surface

Bases: object

A target surface: project(P) the closest point, normal(P) the unit normal at that foot point, both (N, 3) -> (N, 3).

project: Callable[[NDArray[float64]], NDArray[float64]]

normal: Callable[[NDArray[float64]], NDArray[float64]]

name: str = 'surface'

__init__(project, normal, name='surface')

Parameters:

• project (Callable[[NDArray[float64]], NDArray[float64]])

• normal (Callable[[NDArray[float64]], NDArray[float64]])

• name (str)

Return type:

None

class nekmeshpy.core.smoother.constraints.Curve

Bases: object

A target curve: project(P) the closest point, tangent(P) the unit tangent at that foot point.

project: Callable[[NDArray[float64]], NDArray[float64]]

tangent: Callable[[NDArray[float64]], NDArray[float64]]

name: str = 'curve'

__init__(project, tangent, name='curve')

Parameters:

• project (Callable[[NDArray[float64]], NDArray[float64]])

• tangent (Callable[[NDArray[float64]], NDArray[float64]])

• name (str)

Return type:

None

nekmeshpy.core.smoother.constraints.plane(point, normal)

The plane through point with normal.

Parameters:

• point (NDArray[float64])

• normal (NDArray[float64])

Return type:

Surface

nekmeshpy.core.smoother.constraints.sphere(center, radius)

The sphere of radius about center; the normal points outward.

Parameters:

• center (NDArray[float64])

• radius (float)

Return type:

Surface

nekmeshpy.core.smoother.constraints.line(point, direction)

The straight line through point along direction.

Parameters:

• point (NDArray[float64])

• direction (NDArray[float64])

Return type:

Curve

nekmeshpy.core.smoother.constraints.circle(center, radius, axis=(0.0, 0.0, 1.0))

The circle of radius about center in the plane normal to axis; the tangent runs counter-clockwise about axis.

Parameters:

• center (NDArray[float64])

• radius (float)

• axis (NDArray[float64])

Return type:

Curve

class nekmeshpy.core.smoother.constraints.Frozen

Bases: object

Never moves.

points: str | NDArray[int64] | NDArray[bool]

__init__(points)

Parameters:

points (str | NDArray[int64] | NDArray[bool])

Return type:

None

class nekmeshpy.core.smoother.constraints.Tangent

Bases: object

Moves only along the boundary: a boundary vertex whose move is cut down to the tangent of the circle through it and its two boundary neighbours (circle_tangent(), from the current positions). Nothing pulls it back onto a curve – that is OnCurve. A corner (see boundary_corners()) or a vertex without two boundary neighbours is frozen instead; a vertex that is not on the boundary is an error.

points: str | NDArray[int64] | NDArray[bool]

__init__(points)

Parameters:

points (str | NDArray[int64] | NDArray[bool])

Return type:

None

class nekmeshpy.core.smoother.constraints.OnSurface

Bases: object

Stays on surface. slide=True moves tangentially then re-projected; slide=False projects once onto the surface and then holds still. Without a surface it is fitted: the input quads touching these vertices, as a triangle surface (surface_from_trimesh()).

points: str | NDArray[int64] | NDArray[bool]

surface: Surface | None = None

slide: bool = True

__init__(points, surface=None, slide=True)

Parameters:

• points (str | NDArray[int64] | NDArray[bool])

• surface (Surface | None)

• slide (bool)

Return type:

None

class nekmeshpy.core.smoother.constraints.OnCurve

Bases: object

Stays on curve, moving along it only when slide=True. Without a curve it is fitted: the boundary chains through these vertices, as meshed, are cut at corners (boundary_corners()) and at boundary vertices not named here, and each piece gets a spline through it (spline_curve(), closed for a loop); a corner is frozen.

points: str | NDArray[int64] | NDArray[bool]

curve: Curve | None = None

slide: bool = True

__init__(points, curve=None, slide=True)

Parameters:

• points (str | NDArray[int64] | NDArray[bool])

• curve (Curve | None)

• slide (bool)

Return type:

None

class nekmeshpy.core.smoother.constraints.KeepLayer

Bases: object

Keep a boundary layer exactly as meshed: every vertex within rings quad-rings of seeds (the interface or wall vertices, seeds included) is frozen, so the layer’s thickness, its growth ratio and its tangential spacing cannot change. The layers keep their geometry exactly at the price of not sliding; one ring protects the first row of cells off the seeds, two the first two. A square-seeking metric will otherwise inflate a thin layer (measured on a blanket-module section: edges next to an interface grew up to 34x).

seeds: str | NDArray[int64] | NDArray[bool]

rings: int = 1

__init__(seeds, rings=1)

Parameters:

• seeds (str | NDArray[int64] | NDArray[bool])

• rings (int)

Return type:

None

class nekmeshpy.core.smoother.constraints.SlideColumns

Bases: object

Rigid columns that slide together along direction: every column is a set of vertices (e.g. the rows across a wall’s thickness) with one shared degree of freedom, so they all move by the same amount along direction – the cells between the rows stay rectangles while the spacing along the wall adapts. Each column also stays on its own straight line. If any member of a column ends up frozen (a corner, a vertex that must not leave the boundary), the whole column is frozen: a rigid unit cannot move with one member pinned.

columns: Sequence[NDArray[int64]]

direction: Sequence[float]

__init__(columns, direction)

Parameters:

• columns (Sequence[NDArray[int64]])

• direction (Sequence[float])

Return type:

None

class nekmeshpy.core.smoother.constraints.Resolved

Bases: NamedTuple

Per-vertex mobility, resolved once. kind[i] is FREE/ON_SURFACE/ON_CURVE/FROZEN/TANGENT; target[i] indexes surfaces or curves for the first two; anchor is where a frozen vertex lives (the input position, or its one-time projection for slide=False).

kind: NDArray[int64]

Alias for field number 0

target: NDArray[int64]

Alias for field number 1

anchor: NDArray[float64]

Alias for field number 2

surfaces: tuple[Surface, ...]

Alias for field number 3

curves: tuple[Curve, ...]

Alias for field number 4

notes: tuple[str, ...]

Human-readable notes on every conflict and how it was resolved, never silently.

tie: NDArray[int64]

tie[i] is the rigid column vertex i belongs to (-1 none); tie_dir[i] its slide direction.

tie_dir: NDArray[float64]

Alias for field number 7

nbr: NDArray[int64]

nbr[i] the two boundary neighbours of a TANGENT vertex (-1 for every other vertex).

bfaces: NDArray[int64]

(F, 4) the outward boundary faces a hex TANGENT vertex with no nbr slides in the tangent plane of (empty on a quad mesh, whose TANGENT vertices all have two boundary neighbours).

counts()

Return type:

dict[str, int]

static __new__(_cls, kind, target, anchor, surfaces, curves, notes, tie, tie_dir, nbr, bfaces)

Create new instance of Resolved(kind, target, anchor, surfaces, curves, notes, tie, tie_dir, nbr, bfaces)

Parameters:

• kind (IntArray)

• target (IntArray)

• anchor (PointArray)

• surfaces (tuple[Surface, ...])

• curves (tuple[Curve, ...])

• notes (tuple[str, ...])

• tie (IntArray)

• tie_dir (PointArray)

• nbr (IntArray)

• bfaces (IntArray)

nekmeshpy.core.smoother.constraints.circle_tangent(xp, x, xn)

Unit tangent at x of the circle through xp, x, xn (rows), oriented xp -> xn.

With a = xp - x and c = xn - x the tangent at the middle point of the circumscribed circle is |a|^2 c - |c|^2 a: perpendicular to the radius, in the plane of the three points, exact on a circle and second-order on uneven spacing. Collinear points give the edge direction; the degenerate case xp == xn (no direction at all) returns a zero row.

Parameters:

• xp (NDArray[float64])

• x (NDArray[float64])

• xn (NDArray[float64])

Return type:

NDArray[float64]

nekmeshpy.core.smoother.constraints.turning_deg(P, ids, nb, normals=None)

Turning angle (degrees) of the boundary at ids, whose neighbours are nb (n, 2). With normals ((len(P), 3), e.g. the area-weighted vertex normals of the mesh) both edges are first projected onto the tangent plane at the vertex, so the turn is the one within the surface and its own bending does not count.

Parameters:

• P (NDArray[float64])

• ids (NDArray[int64])

• nb (NDArray[int64])

• normals (NDArray[float64] | None)

Return type:

NDArray[float64]

nekmeshpy.core.smoother.constraints.boundary_corners(P, adj, corner_deg=30.0, normals=None)

The corners of the boundary graph adj (vertex -> neighbours along it): a vertex without exactly two neighbours, or one whose turn is a jump – above corner_deg and above three times the mean turn of its two neighbours. The second test is what separates a corner from a coarsely sampled smooth curve: on a circle of N points every vertex turns 360/N, which can be far above corner_deg, but its neighbours turn about the same (a 3:1 ellipse on 12 points peaks at 2.4), so the ratio is near 1; at a real corner the neighbours (edge interiors) turn little, so the ratio is large. (A closed loop of only four points is such a circle.) normals makes the turns geodesic, see turning_deg().

Parameters:

• P (NDArray[float64])

• adj (dict[int, list[int]])

• corner_deg (float)

• normals (NDArray[float64] | None)

Return type:

set[int]

nekmeshpy.core.smoother.constraints.spline_curve(X, closed, name='boundary')

A smooth Curve through the ordered points X (2-D or 3-D): a cubic spline in the chord-length parameter, periodic when closed and not-a-knot otherwise (a straight line for two points). The foot point is the nearest of a dense sampling refined by a few Newton steps on the spline itself, kept only where it lowers the distance, so a projection never ends farther than the sample it started from.

Parameters:

• X (NDArray[float64])

• closed (bool)

• name (str)

Return type:

Curve

nekmeshpy.core.smoother.constraints.snap(points, r, weight=1.0)

points with every constrained vertex moved weight of the way to its target (1 = onto it) and every frozen vertex put back at its anchor. A new array.

Parameters:

• points (NDArray[float64])

• r (Resolved)

• weight (float)

Return type:

NDArray[float64]

nekmeshpy.core.smoother.constraints.constrain_direction(d, points, r)

d (P, 3) (a gradient or a step) with the forbidden part removed, using the target’s own normal/tangent at each vertex’s foot point. A new array.

Parameters:

• d (NDArray[float64])

• points (NDArray[float64])

• r (Resolved)

Return type:

NDArray[float64]

nekmeshpy.core.smoother.constraints.vertex_length(points, corners)

(P,) the shortest edge at each vertex – the length unit every step, finite difference and tolerance of the optimizer is measured in.

Parameters:

• points (NDArray[float64])

• corners (NDArray[int64])

Return type:

NDArray[float64]

class nekmeshpy.core.smoother.constraints.OnTargetReport

Bases: NamedTuple

The geometry self-check of one state: how far constrained vertices sit from their targets and how far frozen ones moved, in absolute units and in local edge lengths.

n_checked: int

Alias for field number 0

max_distance: float

Alias for field number 1

mean_distance: float

Alias for field number 2

max_distance_in_edges: float

Alias for field number 3

max_frozen_move: float

Alias for field number 4

tol: float

Alias for field number 5

passed: bool

Alias for field number 6

static __new__(_cls, n_checked, max_distance, mean_distance, max_distance_in_edges, max_frozen_move, tol, passed)

Create new instance of OnTargetReport(n_checked, max_distance, mean_distance, max_distance_in_edges, max_frozen_move, tol, passed)

Parameters:

• n_checked (int)

• max_distance (float)

• mean_distance (float)

• max_distance_in_edges (float)

• max_frozen_move (float)

• tol (float)

• passed (bool)

nekmeshpy.core.smoother.constraints.on_target_report(points, corners, r, *, tol=1e-06)

Distance of every ON_SURFACE/ON_CURVE vertex to its target, and the largest displacement of a frozen vertex from its anchor; passes when the former is within tol local edge lengths and the latter is exactly 0 (to tol of the edge).

Parameters:

• points (NDArray[float64])

• corners (NDArray[int64])

• r (Resolved)

• tol (float)

Return type:

OnTargetReport

nekmeshpy.core.smoother.constraints.surface_from_trimesh(points, tris, *, k=8, name='trimesh')

A Surface read off a triangle mesh – a CAD tessellation or a scanned wall.

project returns the closest point on the k triangles whose centroids are nearest the query (exact point-triangle distance within that candidate set); normal blends the area-weighted vertex normals with the barycentric weights of that foot point, so it is continuous across triangles. It is an approximation of the true surface: the triangles’ chord error remains, sharp edges are smoothed over by the blended normal, and a k that is too small on a very graded tessellation can miss the nearest triangle. Cost: one KD-tree query and k triangle distances per point.

Parameters:

• points (NDArray[float64])

• tris (NDArray[int64])

• k (int)

• name (str)

Return type:

Surface

nekmeshpy.core.smoother.constraints.claimed_surfaces(r)

(vertex ids, surface) for every surface that holds vertices in r – the groups surface_normals reads its normals from.

Parameters:

r (Resolved)

Return type:

list[tuple[NDArray[int64], Surface]]

nekmeshpy.core.smoother.constraints.targets(r)

(kind, targets) for the two kinds that project onto a geometric object.

Parameters:

r (Resolved)

Return type:

list[tuple[int, Sequence[Surface | Curve]]]

nekmeshpy.core.smoother.constraints.unit_rows(v)

Parameters:

v (NDArray[float64])

Return type:

NDArray[float64]

core.smoother.metric

The smoother’s shape metric, its objective and analytic gradient, for any element dimension.

Every reading is taken off the corner Jacobian J = [e_1 ... e_d] of an element corner – the d edge vectors leaving it, right-handed – so one function serves quads (d = 2) and hexes (d = 3) and the rung only builds the edge vectors from its connectivity and scatters the gradient back onto its vertices.

Signed scaled Jacobian (s = 1 only for d orthogonal edges, s <= 0 inverted or degenerate):

s = det J / (|e_1| ... |e_d|)

Shape metric (Knupp’s mean-ratio measure to the power d; phi >= 1, 1 only for a square or cube corner):

phi = (|J|_F^2 / d)^d / det J^2 = K / s^2,        K = ((sum |e_i|^2) / d)^d / prod |e_i|^2

phi is blind to orientation, so inversion is a separate signed penalty kappa max(0, eps - s). The barrier-free phi_delta = K c / s_delta^2 (s_delta = (s + sqrt(s^2 + 4 delta^2)) / 2, c = s_delta(1)^2) removes the wall at a zero Jacobian that a line search cannot cross; delta = 0 is the exact metric. The objective is the sum over elements of the corner mean (1 / 2^d) of w (phi - 1) + kappa max(0, eps - s); a vertex’s gradient involves only its incident elements, so it does not depend on the element count.

A quad is a surface in 3-D with 2-D topology: det J is read through a reference normal n per element, det J = (e_1 x e_2) . n (normal=), so a mesh on a curved surface is never assumed planar. For d = 2 the exact metric is evaluated as (A + B)^2 / (4 |e_1 x e_2|^2) and for d = 3 as K / s^2; the two are the same function, written the way each is best conditioned, and the tests check them against the single determinant formula above.

nekmeshpy.core.smoother.metric.KAPPA = 1000.0

Inversion-penalty weight per corner, on the sum objective (see module docstring).

nekmeshpy.core.smoother.metric.EPS = 0.05

Scaled-Jacobian threshold below which a corner is penalised.

nekmeshpy.core.smoother.metric.DELTA = 0.01

Barrier-removal width (in units of the signed sine s) the optimizer uses to untangle: phi blows up as a corner’s Jacobian goes through zero, so a plain phi is a wall between an inverted corner and a valid one that no line search crosses. Knupp’s remedy replaces s by s_d = (s + sqrt(s^2 + 4 d^2)) / 2, positive and finite for every s and decreasing smoothly through zero. delta = 0 is the exact metric.

class nekmeshpy.core.smoother.metric.CornerMetrics

Bases: NamedTuple

Per-corner readings, each (n_elements, corners_per_element) in the element’s own corner order.

phi: NDArray[float64]

phi >= 1, 1 at a square / cube corner; blind to orientation.

sj: NDArray[float64]

Signed scaled Jacobian in [-1, 1]; <= 0 is inverted or degenerate.

static __new__(_cls, phi, sj)

Create new instance of CornerMetrics(phi, sj)

Parameters:

• phi (FloatArray)

• sj (FloatArray)

nekmeshpy.core.smoother.metric.regularised(s, delta)

(s_d, c, ds_d/ds): the barrier-free sine, the constant s_d(1)^2 that makes a perfect corner score exactly 1, and the derivative of s_d with respect to s.

Parameters:

• s (NDArray[float64])

• delta (float)

Return type:

tuple[NDArray[float64], float, NDArray[float64]]

nekmeshpy.core.smoother.metric.scatter(n, index, values)

Sum values (n_elements, k, 3) onto n vertices by index (n_elements, k) – the vertex assembly a shared-id mesh needs no gather-scatter library for.

Parameters:

• n (int)

• index (NDArray[int64])

• values (NDArray[float64])

Return type:

NDArray[float64]

nekmeshpy.core.smoother.metric.corner_metrics(edges, normals=None, *, delta=0.0)

phi and the signed scaled Jacobian at every corner. edges is the d edge vectors leaving each corner, each (n_elements, corners, 3); normals (n_elements, 3) is required for d = 2. delta > 0 returns the barrier-free phi_d the optimizer descends on; the default is the exact metric the report prints.

Parameters:

• edges (Sequence[NDArray[float64]])

• normals (NDArray[float64] | None)

• delta (float)

Return type:

CornerMetrics

nekmeshpy.core.smoother.metric.corner_gradient(edges, normals=None, *, kappa=1000.0, eps=0.05, delta=0.0, weights=None)

dF/d e_j for every edge vector, each (n_elements, corners, 3) (reference normals held fixed): the gradient with respect to the d edge vectors leaving each corner. The corner’s own vertex receives minus their sum and each edge’s far end its own term – assemble_gradient() does that with the rung’s connectivity.

Parameters:

• edges (Sequence[NDArray[float64]])

• normals (NDArray[float64] | None)

• kappa (float)

• eps (float)

• delta (float)

• weights (NDArray[float64] | None)

Return type:

list[NDArray[float64]]

nekmeshpy.core.smoother.metric.objective_value(m, dim, *, kappa=1000.0, eps=0.05, weights=None)

F = sum_e (1/2^d) sum_c [ w_e (phi - 1) + kappa max(0, eps - s) ] for the readings m – 0 only for an all-square (cube), all-positive mesh. weights (n_elements,) scale the shape term per element (default 1): 0 makes an element’s shape irrelevant while the inversion penalty still protects it. Divide by the element count for the mean the report shows.

Parameters:

• m (CornerMetrics)

• dim (int)

• kappa (float)

• eps (float)

• weights (NDArray[float64] | None)

Return type:

float

nekmeshpy.core.smoother.metric.evaluate(m, dim, *, kappa=1000.0, eps=0.05, weights=None)

(F, n_inverted) from one set of corner readings: the objective and the number of elements with a corner whose signed scaled Jacobian is <= 0.

Parameters:

• m (CornerMetrics)

• dim (int)

• kappa (float)

• eps (float)

• weights (NDArray[float64] | None)

Return type:

tuple[float, int]

nekmeshpy.core.smoother.metric.assemble_gradient(n_points, corner_index, far_index, grads)

The vertex gradient (n_points, 3) from per-corner edge gradients: corner_index (n_elements, k) the vertex at each corner, far_index[j] the vertex at the far end of edge j of that corner, grads[j] from corner_gradient(). The far end receives grads[j], the corner itself -sum_j grads[j].

Parameters:

• n_points (int)

• corner_index (NDArray[int64])

• far_index (Sequence[NDArray[int64]])

• grads (Sequence[NDArray[float64]])

Return type:

NDArray[float64]

core.smoother.quality

Quality readings of a smoothed mesh, from per-corner and per-element numbers any rung can supply.

The rung builds the corner metrics (nekmeshpy.core.smoother.metric), its edge lengths and its neighbour table; the reductions here turn them into the numbers a report prints. The rung keeps its own Quality tuple (its element count is named n_quads / n_hexes) and fills it from reduce_corners(); ShapeStats has the same fields at every rung and is shared.

class nekmeshpy.core.smoother.quality.ShapeStats

Bases: NamedTuple

Edge- and size-based readings that complement a rung’s Quality: aspect ratio (longest over shortest edge of an element), node spacing (the shortest and longest edge anywhere), corner angles, and the size jump between neighbouring elements (size_ratio_*: the larger over the smaller of the two characteristic sizes of every pair of elements sharing an edge or face, so 1 is a perfectly graded mesh). The scaled Jacobian does not see that last one at all.

aspect_max: float

Alias for field number 0

aspect_p95: float

Alias for field number 1

spacing_min: float

Alias for field number 2

spacing_max: float

Alias for field number 3

angle_min_deg: float

Alias for field number 4

angle_max_deg: float

Alias for field number 5

size_ratio_max: float

Alias for field number 6

size_ratio_p95: float

Alias for field number 7

size_ratio_mean: float

Alias for field number 8

static __new__(_cls, aspect_max, aspect_p95, spacing_min, spacing_max, angle_min_deg, angle_max_deg, size_ratio_max, size_ratio_p95, size_ratio_mean)

Create new instance of ShapeStats(aspect_max, aspect_p95, spacing_min, spacing_max, angle_min_deg, angle_max_deg, size_ratio_max, size_ratio_p95, size_ratio_mean)

Parameters:

• aspect_max (float)

• aspect_p95 (float)

• spacing_min (float)

• spacing_max (float)

• angle_min_deg (float)

• angle_max_deg (float)

• size_ratio_max (float)

• size_ratio_p95 (float)

• size_ratio_mean (float)

nekmeshpy.core.smoother.quality.history_row(q)

(n_inverted, sj_min, objective_mean) of a rung’s Quality – what a stage’s history logs.

Parameters:

q (Any)

Return type:

tuple[int, float, float]

nekmeshpy.core.smoother.quality.with_scaled_jacobian(q, sj)

A rung’s Quality q (a corner reading) with its scaled-Jacobian fields – inverted count, min, mean, median and the 1% / 5% / 10% quantiles – replaced by those of the per-element values sj (n,): how a curved mesh is read on the nodes it stores. The shape metric phi and the objective stay the corner reading (the metric is defined on corners).

Parameters:

• q (Any)

• sj (NDArray[float64])

Return type:

Any

nekmeshpy.core.smoother.quality.reduce_corners(m, objective)

The fields of a rung’s Quality from corner readings m and the objective F, in order: element count, inverted count, then the worst-corner scaled Jacobian’s min / mean / median / 1% / 5% / 10% quantiles, the worst-corner phi’s mean and max, and F per element. Per-element values are the worst corner (sj its minimum, phi its maximum).

Parameters:

• m (CornerMetrics)

• objective (float)

Return type:

tuple[int, int, float, float, float, float, float, float, float, float, float]

nekmeshpy.core.smoother.quality.failing_mask(m, edge, *, sj_target=0.05, short_edge=0.1)

(n_elements,) the elements a local stage should work on: a corner with signed scaled Jacobian below sj_target, or an edge (edge (n_elements, k) lengths) shorter than short_edge times the median.

Parameters:

• m (CornerMetrics)

• edge (NDArray[float64])

• sj_target (float)

• short_edge (float)

Return type:

NDArray[bool]

nekmeshpy.core.smoother.quality.corner_angles(u, v)

Angle in degrees between the edge vectors u and v ((n, k, 3) each) at every corner.

Parameters:

• u (NDArray[float64])

• v (NDArray[float64])

Return type:

NDArray[float64]

nekmeshpy.core.smoother.quality.pair_ratios(size, pairs)

(M,) larger-over-smaller ratio of the element sizes size for every neighbour pair in pairs (M, 2).

Parameters:

• size (NDArray[float64])

• pairs (NDArray[int64])

Return type:

NDArray[float64]

nekmeshpy.core.smoother.quality.shape_stats(edge, angles, ratios)

ShapeStats from edge (n, k) lengths, corner angles (degrees, any shape) and the neighbour size ratios (pair_ratios()).

Parameters:

• edge (NDArray[float64])

• angles (NDArray[float64])

• ratios (NDArray[float64])

Return type:

ShapeStats

core.smoother.schedule

Smoothing stages, schedules and the method registry, shared by every rung.

A stage is one smoothing method run once and described by a StageLog. A schedule is an ordered list of (method, knobs). Each rung owns one Registry – the names its run_schedule knows (laplace, optimize, …) – and a third party plugs another backend in with register under a new name. The three preset schedules are the same at every rung.

class nekmeshpy.core.smoother.schedule.StageLog

Bases: NamedTuple

What one smoothing stage did – the part of the round report this stage owns.

method: str

Alias for field number 0

passes_requested: int

Alias for field number 1

passes_done: int

Alias for field number 2

stopped: str

Alias for field number 3

reason: str

Alias for field number 4

before: Any

The rung’s Quality before and after the stage.

after: Any

Alias for field number 6

history: tuple[tuple[int, int, float, float], ...]

Per pass (pass, n_inverted, sj_min, objective_mean) of the accepted iterates.

points: NDArray[float64] | None

The vertex positions after this stage, only when run_schedule(..., record=True) asked for them.

static __new__(_cls, method, passes_requested, passes_done, stopped, reason, before, after, history, points=None)

Create new instance of StageLog(method, passes_requested, passes_done, stopped, reason, before, after, history, points)

Parameters:

• method (str)

• passes_requested (int)

• passes_done (int)

• stopped (str)

• reason (str)

• before (Any)

• after (Any)

• history (tuple[tuple[int, int, float, float], ...])

• points (PointArray | None)

nekmeshpy.core.smoother.schedule.Schedule

A schedule is an ordered list of (method, knobs); the methods are names in a registry.

alias of Sequence[tuple[str, dict[str, Any]]]

class nekmeshpy.core.smoother.schedule.Method

Bases: NamedTuple

One registered smoothing method: fn(mesh, resolved, **knobs) -> (mesh, StageLog). local: the schedule runner hands it elements= (the failing elements) when it runs; weights: it accepts weights= and receives the runner’s.

fn: Callable[[...], tuple[Any, StageLog]]

Alias for field number 0

local: bool

Alias for field number 1

weights: bool

Alias for field number 2

static __new__(_cls, fn, local, weights)

Create new instance of Method(fn, local, weights)

Parameters:

• fn (Callable[..., tuple[Any, StageLog]])

• local (bool)

• weights (bool)

class nekmeshpy.core.smoother.schedule.Registry

Bases: object

The stage names one rung’s schedule runner knows. entries is the live name -> Method table.

__init__(entries=None)

Parameters:

entries (dict[str, Method] | None)

Return type:

None

entries: dict[str, Method]

register(name, fn, *, local=False, weights=False, replace=False)

Make fn available to run() as the stage name – the hook another backend plugs in by. fn(mesh, resolved, **knobs) must return a new mesh and a StageLog, leave mesh untouched and honour resolved. An existing name is refused unless replace=True.

Parameters:

• name (str)

• fn (Callable[[...], tuple[Any, StageLog]])

• local (bool)

• weights (bool)

• replace (bool)

Return type:

None

names()

The stage names, sorted.

Return type:

list[str]

run(mesh, resolved, schedule, *, failing, weights=None, common=None, record=False)

Run the stages in order, each on the previous stage’s result. failing(mesh) picks the elements a local stage works on, at the moment it runs; common knobs (e.g. the reference normals) go to every stage unless the schedule overrides them. Returns the final mesh and one StageLog per stage; mesh is untouched. record=True keeps a copy of the vertex positions after every stage in StageLog.points.

Parameters:

• mesh (Any)

• resolved (Any)

• schedule (Sequence[tuple[str, dict[str, Any]]])

• failing (Callable[[Any], NDArray[bool]])

• weights (NDArray[float64] | None)

• common (dict[str, Any] | None)

• record (bool)

Return type:

tuple[Any, list[StageLog]]

nekmeshpy.core.smoother.schedule.run_carried(registry, mesh, resolved, schedule, *, curved, linearize, carry, failing, weights=None, common=None, record=False)

Registry.run() for a curved mesh. A stage named in curved (tmop) runs on the curved mesh and returns one; every other stage is linear: each maximal run of them is handed the straight-sided linearize(mesh) and the curved nodes follow the result through carry(curved, moved). Carrying is additive (a node keeps its offset from the straight-sided element), so one carry per run of linear stages equals one per stage. mesh is untouched.

Parameters:

• registry (Registry)

• mesh (Any)

• resolved (Any)

• schedule (Sequence[tuple[str, dict[str, Any]]])

• curved (Collection[str])

• linearize (Callable[[Any], Any])

• carry (Callable[[Any, Any], Any])

• failing (Callable[[Any], NDArray[bool]])

• weights (NDArray[float64] | None)

• common (dict[str, Any] | None)

• record (bool)

Return type:

tuple[Any, list[StageLog]]

nekmeshpy.core.smoother.schedule.user_schedule(outer=10, n_lap=3, n_opt=20, *, w=0.97, final_local=True)

The tuned shape – loops of loops: outer times laplace(n_lap) then optimize(n_opt), closed by two local Laplace passes that lift the smallest edges. Only a starting suggestion; the caller proposes and the user edits it.

Parameters:

• outer (int)

• n_lap (int)

• n_opt (int)

• w (float)

• final_local (bool)

Return type:

list[tuple[str, dict[str, Any]]]

nekmeshpy.core.smoother.schedule.untangle_schedule(deltas=(1.0, 0.3, 0.1, 0.03, 0.01), iters=60)

Optimizer stages with a shrinking barrier width: a wide delta first smooths the landscape enough for a folded vertex to be pulled through, then each narrower one sharpens it back toward the exact metric.

Parameters:

• deltas (Sequence[float])

• iters (int)

Return type:

list[tuple[str, dict[str, Any]]]

nekmeshpy.core.smoother.schedule.suggest_schedule(start, *, outer=10)

The schedule to propose for a mesh in the state start (a rung’s Quality): the tuned shape, preceded by the untangling continuation when any element starts inverted – a Laplace pass on a tangled mesh has nothing to repair it with. A suggestion, never silently run.

Parameters:

• start (Any)

• outer (int)

Return type:

list[tuple[str, dict[str, Any]]]

core.smoother.session

The interactive layer of the smoother, shared by every rung: one SmoothSession holds the mesh, its constraints and a stack of accepted states, and every round answers the same questions (RoundReport) before anything is accepted.

The session never edits in place and never accepts silently: SmoothSession.propose() runs a schedule on a copy and returns a Proposal (report + candidate mesh + before/after files); the caller – a person at the prompt, or an assistant reading the report to a person – decides; accept() pushes it, undo() pops it. Meshes are immutable by contract, so a state is a reference and undo costs nothing. review="print" / "auto" in run() are the fast paths that do not wait for an answer.

Everything that depends on the mesh type sits behind a RungAdapter – how a mesh is linearised and its curved nodes rebuilt (a curved mesh stays curved through a round: see the rung’s run_schedule), how constraints resolve, which input checks run, how quality is read, how a round’s files are written, and which stages exist (the rung’s own registry, so a stage registered there – including a third party’s – runs through the session with no change here). The session itself is only the state stack, the round counter and the verdict. core/ imports no rung, so the adapter is the only door to one.

class nekmeshpy.core.smoother.session.Count

Bases: NamedTuple

An integer a rung counts per mesh that must not rise in a round (quads in an overlap, …). note is the %d warning printed while after is above zero.

label: str

Alias for field number 0

before: int

Alias for field number 1

after: int

Alias for field number 2

note: str

Alias for field number 3

static __new__(_cls, label, before, after, note)

Create new instance of Count(label, before, after, note)

Parameters:

• label (str)

• before (int)

• after (int)

• note (str)

class nekmeshpy.core.smoother.session.Audit

Bases: NamedTuple

A rung’s own part of a verdict: extra warnings / failures, structured data the rung’s report type reads back, and text lines appended to section 3 (quality) and section 4 (geometry) of the printout.

warnings: tuple[str, ...]

Alias for field number 0

failures: tuple[str, ...]

Alias for field number 1

data: dict[str, Any]

Alias for field number 2

detail3: tuple[str, ...]

Alias for field number 3

detail4: tuple[str, ...]

Alias for field number 4

static __new__(_cls, warnings=(), failures=(), data={}, detail3=(), detail4=())

Create new instance of Audit(warnings, failures, data, detail3, detail4)

Parameters:

• warnings (tuple[str, ...])

• failures (tuple[str, ...])

• data (dict[str, Any])

• detail3 (tuple[str, ...])

• detail4 (tuple[str, ...])

class nekmeshpy.core.smoother.session.RoundFields

Bases: NamedTuple

The fields of a RoundReport, in order.

round: int

Alias for field number 0

constraints: str

Alias for field number 1

stages: tuple[StageLog, ...]

Alias for field number 2

before: Any

Alias for field number 3

after: Any

Alias for field number 4

shape_before: Any

Alias for field number 5

shape_after: Any

Alias for field number 6

geometry_before: Any

Alias for field number 7

geometry_after: Any

Alias for field number 8

propagate: str

Alias for field number 9

undo: str

Alias for field number 10

warnings: tuple[str, ...]

Alias for field number 11

failures: tuple[str, ...]

Alias for field number 12

files: dict[str, str]

Alias for field number 13

counts: tuple[Count, ...]

Alias for field number 14

data: dict[str, Any]

Alias for field number 15

detail3: tuple[str, ...]

Alias for field number 16

detail4: tuple[str, ...]

Alias for field number 17

static __new__(_cls, round, constraints, stages, before, after, shape_before, shape_after, geometry_before, geometry_after, propagate, undo, warnings, failures, files, counts=(), data={}, detail3=(), detail4=())

Create new instance of RoundFields(round, constraints, stages, before, after, shape_before, shape_after, geometry_before, geometry_after, propagate, undo, warnings, failures, files, counts, data, detail3, detail4)

Parameters:

• round (int)

• constraints (str)

• stages (tuple[StageLog, ...])

• before (Any)

• after (Any)

• shape_before (Any)

• shape_after (Any)

• geometry_before (Any)

• geometry_after (Any)

• propagate (str)

• undo (str)

• warnings (tuple[str, ...])

• failures (tuple[str, ...])

• files (dict[str, str])

• counts (tuple[Count, ...])

• data (dict[str, Any])

• detail3 (tuple[str, ...])

• detail4 (tuple[str, ...])

class nekmeshpy.core.smoother.session.RoundReport

Bases: RoundFields

What a round owes, plus the verdict. failures block acceptance; warnings do not. Printing it gives the full, honest before/after table. A rung subclasses this (no new fields) to set the words (noun / plural / angle_word) and to read its Audit data back as attributes.

noun = 'element'

plural = 'elements'

angle_word = 'corner'

propagate_label = 'propagate'

property passed: bool

metric_rows()

Return type:

list[tuple[str, float, float, str]]

class nekmeshpy.core.smoother.session.Proposal

Bases: NamedTuple

A round that has been run but not accepted: the candidate mesh, its report, the schedule.

mesh: Any

Alias for field number 0

report: RoundReport

Alias for field number 1

schedule: list[tuple[str, dict[str, Any]]]

Alias for field number 2

static __new__(_cls, mesh, report, schedule)

Create new instance of Proposal(mesh, report, schedule)

Parameters:

• mesh (Any)

• report (RoundReport)

• schedule (list[tuple[str, dict[str, Any]]])

class nekmeshpy.core.smoother.session.RungAdapter

Bases: object

Everything the session needs from one mesh type. Subclass it per rung; the defaults are the “nothing to do” answers, so a rung overrides only what it has. report_type is the RoundReport subclass it fills.

report_type

alias of RoundReport

propagate = ''

Section 5 of the report (“propagate to the blocks above?”); empty when the question does not arise.

no_rebuild = 'this rung cannot rebuild a high-order mesh yet'

The explanation printed in place of a rebuilt high-order mesh when the rung has none.

preflight(mesh)

Refuse or warn about the input mesh (a gap gate, a curved mesh the rung cannot take). Raises.

Parameters:

mesh (Any)

Return type:

None

linearize(mesh)

The straight-sided order-1 shadow of a curved mesh (what the linear stages of a schedule run on).

Parameters:

mesh (Any)

Return type:

Any

resolve(mesh, constraints, boundary)

The per-vertex mobility (Resolved) of constraints on mesh.

Parameters:

• mesh (Any)

• constraints (list[Any])

• boundary (str)

Return type:

Any

prepare(mesh, resolved)

Called once with the linear mesh and its first resolution (reference normals, input checks).

Parameters:

• mesh (Any)

• resolved (Any)

Return type:

None

rebuild(base, points, resolved, order)

The curved mesh of the smoothed vertices points carried on base (see SmoothSession.rebuild_high_order()). Returns (mesh, report).

Parameters:

• base (Any)

• points (Any)

• resolved (Any)

• order (int | None)

Return type:

tuple[Any, Any]

quality(mesh)

Parameters:

mesh (Any)

Return type:

Any

shape_stats(mesh)

Parameters:

mesh (Any)

Return type:

Any

on_target(mesh, resolved, tol)

Parameters:

• mesh (Any)

• resolved (Any)

• tol (float)

Return type:

Any

counts(mesh)

(label, value, note) per integer the verdict watches (it may not rise); empty by default.

Parameters:

mesh (Any)

Return type:

tuple[tuple[str, int, str], …]

audit(before, after, resolved)

The rung’s own warnings / failures and report sections for the round before -> after.

Parameters:

• before (Any)

• after (Any)

• resolved (Any)

Return type:

Audit

describe(resolved, *, sj_target, geom_tol)

Section 1 of the report: what the constraints were.

Parameters:

• resolved (Any)

• sj_target (float)

• geom_tol (float)

Return type:

str

run(mesh, resolved, schedule, *, sj_target, weights)

Run schedule on mesh through the rung’s registry (mesh is untouched).

Parameters:

• mesh (Any)

• resolved (Any)

• schedule (Sequence[tuple[str, dict[str, Any]]])

• sj_target (float)

• weights (Any)

Return type:

tuple[Any, list[StageLog]]

write(directory, before, after, resolved)

Write the round’s files into directory (exists); returns {"before": path, "after": path}.

Parameters:

• directory (str)

• before (Any)

• after (Any)

• resolved (Any)

Return type:

dict[str, str]

class nekmeshpy.core.smoother.session.SmoothSession

Bases: object

One smoothing job on one mesh, through a RungAdapter. A rung’s session subclasses this with its own constructor (the options its adapter takes) and binds the adapter; see the module docstring for the flow. A curved mesh stays curved (keep_curved); pass keep_curved=False to smooth its straight-sided shadow.

__init__(mesh, constraints, adapter, *, boundary='frozen', sj_target=0.05, geom_tol=1e-06, weights=None, out_dir=None, name='smooth', keep_curved=True)

Parameters:

• mesh (Any)

• constraints (list[Any])

• adapter (RungAdapter)

• boundary (str)

• sj_target (float)

• geom_tol (float)

• weights (Any)

• out_dir (str | None)

• name (str)

• keep_curved (bool)

Return type:

None

original: Any

property current: Any

property round: int

property resolved: Any

rebuild_high_order(order=None)

The curved mesh of the current accepted vertices: the original high-order nodes carried along and re-projected onto the targets. With a linear input and order=N the mesh is lifted to order N. The session state is not changed. Raises NotImplementedError at a rung without a rebuild.

Parameters:

order (int | None)

Return type:

tuple[Any, Any]

quality()

Return type:

Any

suggest()

The schedule to propose for the current state (never run silently).

Return type:

list[tuple[str, dict[str, Any]]]

set_constraints(constraints, *, boundary=None)

Change the constraints – a state of its own, so undo() restores the old ones.

Parameters:

• constraints (list[Any])

• boundary (str | None)

Return type:

None

propose(schedule=None)

Run schedule (default: suggest()) on a copy; the session is unchanged.

Parameters:

schedule (Sequence[tuple[str, dict[str, Any]]] | None)

Return type:

Proposal

accept(proposal)

Push an accepted proposal (refused if it failed or is stale).

Parameters:

proposal (Proposal)

Return type:

None

reject(proposal)

Drop a proposal. Nothing was pushed, so there is nothing to undo; the method exists so a reviewer’s answer has a name (and is checked for staleness like accept()).

Parameters:

proposal (Proposal)

Return type:

None

undo(n=1)

Pop n accepted states (constraint changes count as states); the input is the floor.

Parameters:

n (int)

Return type:

Any

history()

The trend over accepted states: inverted elements, min scaled J, objective per element.

Return type:

str

run(rounds, *, review='print', adjust='same', schedule=None, stall=0.001, ideal=1e-06)

The fast paths, which do not wait for an answer. "print" runs rounds rounds on the same trajectory (or adjust-ed each time), accepting every passing round, printing each report, and stopping at the first FAIL. "auto" additionally stops when the mean objective is already below ideal (a square lattice is 0; nothing left to gain), when a round improves it by less than stall (relative), or when it raises a regression warning – the interactive path is propose() + a human + accept().

Parameters:

• rounds (int)

• review (str)

• adjust (str)

• schedule (Sequence[tuple[str, dict[str, Any]]] | None)

• stall (float)

• ideal (float)

Return type:

list[RoundReport]

nekmeshpy.core.smoother.session.adjust_schedule(schedule, mode)

The “more aggressive / more conservative / same” step, as one table of knobs and never as a loosened constraint: Laplace shrink w (0.99 / keep / 0.92), Laplace passes (2 / keep / 6), optimizer iterations (10 / keep / 40) and optimizer step cap alpha_max (0.25 / keep / keep).

Parameters:

• schedule (Sequence[tuple[str, dict[str, Any]]])

• mode (str)

Return type:

list[tuple[str, dict[str, Any]]]

nekmeshpy.core.smoother.session.build_report(adapter, before, after, resolved, logs, *, rnd=0, undo='', propagate='', files=None, sj_target=0.05, geom_tol=1e-06)

The RoundReport of smoothing before into after under resolved, every reading taken through adapter and the verdict from verdict().

Parameters:

• adapter (RungAdapter)

• before (Any)

• after (Any)

• resolved (Any)

• logs (list[StageLog])

• rnd (int)

• undo (str)

• propagate (str)

• files (dict[str, str] | None)

• sj_target (float)

• geom_tol (float)

Return type:

RoundReport

nekmeshpy.core.smoother.session.verdict(report_type, logs, before, after, shape_before, shape_after, geometry, counts=(), audit=None)

(warnings, failures) of a round – the one place it is decided. Fails: off target (geometry is the after reading), more inverted elements, a watched count that rose, anything the rung’s audit fails. Warns: what remains, any metric that got worse, a stage rolled back, Laplace and the optimizer fighting.

Parameters:

• report_type (type[RoundReport])

• logs (list[StageLog])

• before (Any)

• after (Any)

• shape_before (Any)

• shape_after (Any)

• geometry (Any)

• counts (tuple[Count, ...])

• audit (Audit | None)

Return type:

tuple[list[str], list[str]]

core.smoother.descent

One nonlinear conjugate-gradient driver with a line search, for every rung’s shape optimizer.

The optimizer is the same loop at the quad and the hex rung: Fletcher-Reeves on the objective, the constraints enforced every step. What differs – how the objective and its gradient are evaluated on the elements in play, how a gradient or a step is cut down to what the targets allow, how a trial is snapped back onto the targets – is passed in as callables, so this module imports no mesh type.

Everything is measured in one unit, the local edge length l_v (shortest incident edge, frozen at stage start): the variables are u_v = dx_v / l_v, the search direction is normalised so its largest move is one l_v, and the line-search variable alpha is that move in edge lengths, capped at alpha_max. Nothing depends on the element count or the model’s units.

class nekmeshpy.core.smoother.descent.DescentResult

Bases: NamedTuple

What minimize() ended with: the last accepted positions P (not yet snapped a final time), the objective f and inverted count ninv there, one (iteration, n_inverted, sj_min, objective_mean) per accepted iterate, why it stopped and how many iterations were accepted.

P: NDArray[float64]

Alias for field number 0

f: float

Alias for field number 1

ninv: int

Alias for field number 2

history: list[tuple[int, int, float, float]]

Alias for field number 3

stopped: str

Alias for field number 4

reason: str

Alias for field number 5

done: int

Alias for field number 6

static __new__(_cls, P, f, ninv, history, stopped, reason, done)

Create new instance of DescentResult(P, f, ninv, history, stopped, reason, done)

Parameters:

• P (PointArray)

• f (float)

• ninv (int)

• history (list[tuple[int, int, float, float]])

• stopped (str)

• reason (str)

• done (int)

nekmeshpy.core.smoother.descent.minimize(P, *, evaluate, gradient, project, snap, summary, ell, movable, tangent, n_used, iters, alpha_max=1.0, restart='halve', tol_step=0.0001, tol_rel=1e-06)

Conjugate-gradient descent from P (already snapped onto the targets).

evaluate(Q) is (F, n_inverted) over the elements in play; gradient(P) the analytic dF/dx (P, 3) over the same; project(d, P) removes from a gradient or step d whatever the targets forbid at positions P; snap(Q) puts a trial back on the targets; summary(P) is the (n_inverted, sj_min, objective_mean) of the whole mesh logged per accepted iterate. ell (P, 1) is the local length unit, movable the vertices that may move, tangent those whose direction is re-cut after the CG mix, n_used the element count for the mean-excess test.

Each step: take the gradient, project it, form the CG direction, and line-search alpha on F(snap(x + alpha l s)) – the objective after the snap, so the search sees what the mesh will actually become. A trial that would add an inverted element is never accepted.

Stops on iters, a step below tol_step edge lengths, a relative drop of F below tol_rel, a mean excess below 1e-10 (already ideal), a zero projected gradient, or a line search that finds no step lowering F. restart="halve" halves beta on every second iteration, "none" is plain Fletcher-Reeves.

Parameters:

• P (NDArray[float64])

• evaluate (Callable[[NDArray[float64]], tuple[float, int]])

• gradient (Callable[[NDArray[float64]], NDArray[float64]])

• project (Callable[[NDArray[float64], NDArray[float64]], NDArray[float64]])

• snap (Callable[[NDArray[float64]], NDArray[float64]])

• summary (Callable[[NDArray[float64]], tuple[int, float, float]])

• ell (NDArray[float64])

• movable (NDArray[bool])

• tangent (NDArray[bool])

• n_used (int)

• iters (int)

• alpha_max (float)

• restart (str)

• tol_step (float)

• tol_rel (float)

Return type:

DescentResult

core.smoother.blocks

Per-element node blocks of a curved mesh: how a native stage carries the curved nodes, and how they reach TMOP.

Per-element node blocks, as arrays: the rung-agnostic half of smoothing a curved (order N) mesh.

A curved mesh stores its nodes once per entity (a shared vertex, a shared edge’s interior, a shared face’s interior, a private element interior); the smoother thinks in per-element blocks (E, (N+1)**dim, 3) in nekmeshpy’s lexicographic Gauss-Lobatto order (axis 0 fastest). A rung turns a mesh into blocks with its own element_blocks and back with its own with_blocks; this module holds what is the same at every rung:

• scatter_blocks() – blocks back onto the flat node table, through an id table that the rung gets by running its own element_blocks on node ids instead of coordinates (so every orientation convention of the storage is reused, none re-derived, and no coordinate search is made);

• carry_blocks() – the curved nodes follow a displacement of the element’s corners, interpolated multilinearly, so a shared edge or face (which sees only its own corners) stays conformal;

• node_fixed() – which nodes inherit a vertex constraint: a node on an entity is held in a component when every vertex of that entity holds it, and an element’s own interior is never held;

• min_det_blocks() – the Jacobian determinant of a block on a Gauss-Lobatto lattice (the quadrature points of the TMOP run).

nekmeshpy.core.smoother.blocks.lattice_weights(order, dim)

((order+1)**dim, 2**dim) multilinear weights of an element’s corners at its Gauss-Lobatto lattice (corners in the rung’s connectivity order).

Parameters:

• order (int)

• dim (int)

Return type:

NDArray[float64]

nekmeshpy.core.smoother.blocks.carry_blocks(blocks, displacement, order, dim)

blocks (E, M, 3) moved by the corner displacement (E, 2**dim, 3), interpolated multilinearly to every node: x_new = x_old + sum_c w_c(x) d_c. A node on an edge or face sees only that entity’s corners (the other weights are exactly zero), so two elements sharing it agree; an element whose corners did not move is returned bit for bit.

Parameters:

• blocks (NDArray[float64])

• displacement (NDArray[float64])

• order (int)

• dim (int)

Return type:

NDArray[float64]

nekmeshpy.core.smoother.blocks.scatter_blocks(ids, blocks, flat)

The flat node table flat (n, 3) with blocks (E, M, 3) written at ids (E, M). A node shared by several elements is written by each; the blocks of a conformal mesh agree (to rounding), and the last writer wins. Nodes no element references keep their value in flat. A copy.

Parameters:

• ids (NDArray[int64])

• blocks (NDArray[float64])

• flat (NDArray[float64])

Return type:

NDArray[float64]

nekmeshpy.core.smoother.blocks.node_fixed(vertex_fixed, cells, order, dim)

(E, M, k) the components held at every node of every block, from vertex_fixed (P, k) (the components held at each vertex) and the corner ids cells (E, 2**dim).

A node lies on one lowest-dimensional entity of its element – a vertex, an edge, a face, or the interior. It is held in a component when every vertex of that entity holds it; an element’s interior is held in a component only when all its corners do (a frozen element stays as it is), so it is free otherwise. The entity is read from the lattice index: an index at an end (0 or order) pins that axis, and the entity’s vertices are the corners agreeing on every pinned axis.

Parameters:

• vertex_fixed (NDArray[bool])

• cells (NDArray[int64])

• order (int)

• dim (int)

Return type:

NDArray[bool]

nekmeshpy.core.smoother.blocks.min_det_blocks(coords, order, dim, npts)

(E,) the smallest Jacobian determinant of each block over the npts-per-axis Gauss-Lobatto lattice.

coords (E, M, dim) are the node coordinates in a frame of the element’s own dimension; the map is the order-order tensor Lagrange interpolant, differentiated analytically, the determinant taken with respect to the reference cube [0, 1]**dim.

Parameters:

• coords (NDArray[float64])

• order (int)

• dim (int)

• npts (int)

Return type:

NDArray[float64]

core.smoother.mfem

The optional MFEM / TMOP bridge, shared by both rungs’ ("tmop", {...}) stage. It needs PyMFEM (pip install mfem); the module itself imports without it, and calling it without says so. The guide: TMOP via the MFEM backend.

The MFEM / TMOP bridge: one smoothing method for quads and hexes, written once.

TMOP (the Target-Matrix Optimization Paradigm of MFEM) minimises a sum over elements of a metric of the map from a target element to the actual one, with a Newton or L-BFGS solver. This module hands a mesh to it as plain arrays and takes the vertex positions back, so the quad rung and the hex rung each register a three-line wrapper (("tmop", {...}) in a schedule) and share everything else here.

Optional. MFEM is imported lazily, inside the functions that need it, so this module (and every smoother module that registers it) imports without PyMFEM; calling it without raises an ImportError naming the install (pip install mfem, the optional smooth extra). PyMFEM is a large binary wheel and wraps a C++ library that aborts the process on some invalid input, so the bridge checks what it can before handing over (tmop() refuses an inverted input for every metric but the untangling 313).

What it does, in order:

• Arrays in. points (N, 3) and cells (C, 4) quads or (C, 8) hexes in nekmeshpy’s corner order, which is MFEM’s SQUARE / CUBE vertex order (NEK_TO_MFEM states the permutation explicitly, and the tests check that a round trip keeps every Jacobian positive). TMOP needs the element dimension to equal the space dimension, so a quad mesh must be planar: its points are rewritten in an orthonormal frame of that plane (plane_frame()) and a quad mesh on a curved 3-D surface is refused – smooth it with the rung’s own optimizer.

• Constraints to essential DOFs (fixed_components()), one rule per mobility kind of Resolved:

  • FROZEN -> every component fixed;

  • ON_SURFACE on a plane whose normal is a coordinate axis -> that one component fixed, the vertex slides in the plane (the quad rung has no surface constraint to speak of: a vertex “on” the mesh’s own plane stays free);

  • ON_CURVE on a line along a coordinate axis -> the other components fixed, the vertex slides on the line;

  • anything else – a curved surface or curve, TANGENT, a rigid slide_column – is held fixed (for now: the optimizer cannot yet follow a curved target); the count is in the stage’s reason.

  Every constrained vertex is first projected onto its target (snap()), so a component is fixed at the target’s coordinate, not wherever the input happened to be.

• Targets – what “good” means per element, as the target Jacobian W of the metric:

  • "input": W is a box with the element’s own edge lengths and ideal (orthogonal) skew – the element is asked to become orthogonal at the size and aspect it already has, so a graded mesh stays graded (MFEM’s GIVEN_SHAPE_AND_SIZE from per-element size and aspect-ratio fields);

  • "input_smoothed": the same, with each element’s edge lengths first relaxed sweeps times toward its face neighbours’ (weight sweep_weight), so the grading itself is smoothed too;

  • "ideal": the ideal shape at the element’s own size (IDEAL_SHAPE_GIVEN_SIZE): wants cubes / squares, which inflates a thin layer – use with designed and a quality you can afford.

  designed elements (a boundary layer, a region to leave alone) are instead given their own input Jacobian as target, so only their shape is asked to stay and a layer keeps its thickness when its two faces are held. An element whose every vertex is fully fixed carries zero weight.

• Metric: metric is a TMOP metric id; the defaults are 2 (2-D shape) and 303 (3-D shape). 80 and 333 (shape + size, metric_weight) work; 313 untangles (shifted barrier, finite differences, one integrator: target="ideal").

• Solver: "newton" (MINRES, l1-Jacobi) or "lbfgs".

Orders. tmop() is the order-1 entry (the vertices are the nodes). tmop_blocks() is the curved one: a mesh of order N goes in as per-element node blocks, an MFEM mesh with an H1 Gauss-Lobatto nodal space of order N is built over the order-1 cells, every node is a degree of freedom, and the optimised blocks come back at the same order. Before optimising it reads the geometry back from MFEM’s own element map and compares it with the blocks (1e-12 of the extent), so a disagreement between the two libraries’ node orders is an error, not a bad mesh. A vertex is held as at order 1 (at its input coordinates), the nodes on an edge or face inherit the hold of the vertices of that entity, interior nodes are free. Targets are built from each element’s corner extents, the quadrature order defaults to 2N + 2 (Gauss-Lobatto, so it contains the nodes), and the log carries a BlockReading: minimum determinant at the quadrature points, scaled Jacobian at the nodes and on the finer sample_order lattice. Not available at order N: sliding along a curved target (held fixed, as at order 1).

nekmeshpy.core.smoother.mfem.available()

Whether PyMFEM can be imported (nothing is imported to find out).

Return type:

bool

nekmeshpy.core.smoother.mfem.NEK_TO_MFEM: dict[int, NDArray[int64]] = {2: array([0, 1, 2, 3]), 3: array([0, 1, 2, 3, 4, 5, 6, 7])}

The MFEM vertex order of a quad (SQUARE) and a hex (CUBE), as nekmeshpy’s own corner indices: mfem_vertex[k] = corner[NEK_TO_MFEM[dim][k]]. It is the identity today (a quad runs 0-1-2-3 cyclic, a hex 0-3 on the bottom face and 4-7 above them); it is a table so that the one place the two libraries’ conventions meet is stated, and tested, instead of assumed.

class nekmeshpy.core.smoother.mfem.Frame

Bases: NamedTuple

An orthonormal frame of the mesh’s space: coords = (points - origin) @ basis (N, dim), and points = origin + coords @ basis.T for a mesh that lies in the frame’s span.

origin: NDArray[float64]

Alias for field number 0

basis: NDArray[float64]

Alias for field number 1

static __new__(_cls, origin, basis)

Create new instance of Frame(origin, basis)

Parameters:

• origin (FloatArray)

• basis (FloatArray)

nekmeshpy.core.smoother.mfem.plane_frame(points, cells)

The frame the optimizer works in. A hex mesh keeps the coordinate axes. A quad mesh must lie in one plane (to 1e-9 of its extent; otherwise a ValueError, since TMOP needs the element dimension to equal the space dimension): when that plane is a coordinate plane the frame is its two axes, else the plane’s principal axes; the second one is flipped if needed so the quads keep positive area.

Parameters:

• points (NDArray[float64])

• cells (NDArray[int64])

Return type:

Frame

class nekmeshpy.core.smoother.mfem.FixedMap

Bases: NamedTuple

The essential components of a Resolved.

points: NDArray[float64]

(N, 3) the input with every constrained vertex projected onto its target.

fixed: NDArray[bool]

(N, dim) True where the component (in the working frame) is held.

notes: tuple[str, ...]

One line per rule that applied, e.g. "12 vertices slide in axis-aligned planes".

counts: dict[str, int]

frozen, plane, line, held_curved, held_other.

Type:

{rule: vertex count}

static __new__(_cls, points, fixed, notes, counts)

Create new instance of FixedMap(points, fixed, notes, counts)

Parameters:

• points (PointArray)

• fixed (BoolArray)

• notes (tuple[str, ...])

• counts (dict[str, int])

nekmeshpy.core.smoother.mfem.fixed_components(points, r, basis=None)

Map r to essential components in the frame basis (3, dim) (default the 3 coordinate axes).

The rules are in the module docstring. A target counts as an axis-aligned plane when its normal is the same unit coordinate axis at every vertex it holds (a sphere or a triangulated surface is not, even where one vertex happens to sit on an axis), and an axis-aligned line likewise for its tangent.

Parameters:

• points (NDArray[float64])

• r (Resolved)

• basis (NDArray[float64] | None)

Return type:

FixedMap

nekmeshpy.core.smoother.mfem.element_extents(coords, cells)

(C, dim) the mean edge length of each element along each of its reference axes.

Parameters:

• coords (NDArray[float64])

• cells (NDArray[int64])

Return type:

NDArray[float64]

nekmeshpy.core.smoother.mfem.regions(element_tags, designed=None)

(attributes, designed_mask) from a mesh’s dense element tags ((C,) strings, "" untagged): each distinct tag is an integer attribute (1-based, in sorted order), and designed – a tag name, a list of names, a boolean (C,) mask or None – becomes the mask tmop() takes. A name no element carries is an error.

Parameters:

• element_tags (Any)

• designed (Any)

Return type:

tuple[NDArray[int64], NDArray[bool] | None]

nekmeshpy.core.smoother.mfem.to_mfem_mesh(coords, cells, attributes=None)

An in-memory mfem.Mesh of coords (N, dim) and cells (quads or hexes in nekmeshpy corner order), vertex i of the array being vertex i of the mesh and element e being element e. attributes (C,) ints become the element attributes (default 1).

Parameters:

• coords (NDArray[float64])

• cells (NDArray[int64])

• attributes (NDArray[int64] | None)

Return type:

Any

nekmeshpy.core.smoother.mfem.from_mfem_mesh(mesh)

(coords (N, dim), cells in nekmeshpy corner order, attributes) read back from an mfem.Mesh – the inverse of to_mfem_mesh() (the vertex coordinates, not any nodal grid function).

Parameters:

mesh (Any)

Return type:

tuple[NDArray[float64], NDArray[int64], NDArray[int64]]

nekmeshpy.core.smoother.mfem.min_det(coords, cells, n=3)

The smallest Jacobian determinant of the multilinear elements over an n-per-axis lattice (Gauss-Lobatto-like: ends included) – the check made before MFEM is entered, since a negative determinant makes most TMOP metrics abort the process. n = 3 covers the corners and the midpoints.

Parameters:

• coords (NDArray[float64])

• cells (NDArray[int64])

• n (int)

Return type:

float

class nekmeshpy.core.smoother.mfem.TmopEnergy

Bases: NamedTuple

The reading a stage logs when no rung quality was supplied: the TMOP energy only.

energy: float

Alias for field number 0

n_inverted: int

Alias for field number 1

sj_min: float

Alias for field number 2

objective_mean: float

Alias for field number 3

static __new__(_cls, energy, n_inverted=-1, sj_min=nan, objective_mean=nan)

Create new instance of TmopEnergy(energy, n_inverted, sj_min, objective_mean)

Parameters:

• energy (float)

• n_inverted (int)

• sj_min (float)

• objective_mean (float)

nekmeshpy.core.smoother.mfem.tmop(points, cells, resolved, *, attributes=None, designed=None, target='input', metric=None, metric_weight=0.5, solver='newton', iters=50, rtol=1e-08, sweeps=5, sweep_weight=0.5, quad_order=4, order=1, quality=None, verbose=False)

One TMOP run on the order-1 mesh (points, cells) under resolved: the new vertex positions (N, 3) and a StageLog (its reason carries the energy before and after, the Newton / L-BFGS iterations and what was fixed). The rules are in the module docstring.

attributes (C,) ints label the elements (regions); designed selects the elements that keep their own input shape as target: a boolean (C,) mask, or integer attribute values (needs attributes). target is "input", "input_smoothed" (sweeps relaxations of weight sweep_weight on the target edge lengths) or "ideal". metric defaults to DEFAULT_METRIC; metric_weight is the weight of the combined metrics 80 and 333. solver is "newton" or "lbfgs", stopping at iters or a relative residual rtol. quad_order is the Gauss-Lobatto quadrature order (4 includes the corners, so the determinant check of the line search covers the geometry a .re2 file keeps). quality(points), when given, reads the rung’s own Quality for the log (and enables the rollback: a run that ends with more inverted elements than it began with is refused and the input returned).

Raises ImportError without PyMFEM, ValueError for a non-planar quad mesh, an inverted input (except metric 313) or a bad knob, or order != 1 (tmop_blocks() is the curved entry).

Parameters:

• points (NDArray[float64])

• cells (NDArray[int64])

• resolved (Resolved)

• attributes (NDArray[int64] | None)

• designed (NDArray[bool] | Sequence[int] | None)

• target (str)

• metric (int | None)

• metric_weight (float)

• solver (str)

• iters (int)

• rtol (float)

• sweeps (int)

• sweep_weight (float)

• quad_order (int)

• order (int)

• quality (Callable[[NDArray[float64]], Any] | None)

• verbose (bool)

Return type:

tuple[NDArray[float64], StageLog]

class nekmeshpy.core.smoother.mfem.BlockReading

Bases: NamedTuple

The reading a curved stage logs, taken on the node blocks. det_min / n_inverted are at the TMOP quadrature points (the Jacobian the line search guards), sj_min at the stored nodes and sj_sampled on the sample_order lattice – the order the solver will run at, which is the one that can see a fold between the nodes.

energy: float

Alias for field number 0

n_inverted: int

Alias for field number 1

sj_min: float

Alias for field number 2

objective_mean: float

Alias for field number 3

det_min: float

Alias for field number 4

sj_sampled: float

Alias for field number 5

sample_order: int

Alias for field number 6

static __new__(_cls, energy, n_inverted, sj_min, objective_mean, det_min, sj_sampled, sample_order)

Create new instance of BlockReading(energy, n_inverted, sj_min, objective_mean, det_min, sj_sampled, sample_order)

Parameters:

• energy (float)

• n_inverted (int)

• sj_min (float)

• objective_mean (float)

• det_min (float)

• sj_sampled (float)

• sample_order (int)

nekmeshpy.core.smoother.mfem.tmop_blocks(blocks, cells, resolved, *, order, attributes=None, designed=None, target='input', metric=None, metric_weight=0.5, solver='newton', iters=50, rtol=1e-08, sweeps=5, sweep_weight=0.5, quad_order=None, sample_order=None, verbose=False)

One TMOP run on a curved mesh given as node blocks: every Gauss-Lobatto node is a degree of freedom, and the new blocks come back at the same order.

blocks (E, (order+1)**dim, 3) are the elements’ nodes in nekmeshpy’s lexicographic order (a rung’s element_blocks); cells (E, 2**dim) the corner ids into resolved (the order-1 topology). Shared nodes must agree between the elements that share them – checked, by reading the blocks back from MFEM before anything is optimised (to 1e-12 of the mesh extent), which also catches any mismatch between the two libraries’ node orders.

What moves. Every node except the held ones. A vertex is held as at order 1 (fixed_components(): frozen in every component, an axis-aligned plane in its normal, an axis-aligned line in the two others, anything curved or tangent in all of them), at its input coordinates – vertices are not re-snapped onto their targets here. The nodes on an edge or face inherit the hold of their entity’s vertices (node_fixed(): held in a component when every vertex of the entity holds it), so a frozen wall keeps its curved edge and face nodes bit for bit and a plane keeps its nodes on the plane; interior nodes are free. Sliding along a curved target is not available, the vertices and entities on one are held.

Targets, metric, solver as tmop(); the targets are built from each element’s corner extents (one target per element, read at every quadrature point). quad_order is the Gauss-Lobatto quadrature order (default 2*order + 2, which includes the nodes and grows with the order); sample_order (default max(order, SCAN_ORDER)) is the lattice the logged scaled Jacobian is also read on. The log carries a BlockReading before and after, and a run that raises the number of elements with a non-positive determinant at the quadrature points, or does not lower the energy, is refused and the input returned.

Parameters:

• blocks (NDArray[float64])

• cells (NDArray[int64])

• resolved (Resolved)

• order (int)

• attributes (NDArray[int64] | None)

• designed (NDArray[bool] | Sequence[int] | None)

• target (str)

• metric (int | None)

• metric_weight (float)

• solver (str)

• iters (int)

• rtol (float)

• sweeps (int)

• sweep_weight (float)

• quad_order (int | None)

• sample_order (int | None)

• verbose (bool)

Return type:

tuple[NDArray[float64], StageLog]