Quality metrics + smoothing

Per-element quality computations and the pre-quadrangulation smoother.

mesh_quality

admesh.mesh_quality

mesh_quality(p: ArrayLike, t: ArrayLike, element: str = 'triangle') -> tuple[float, float, NDArray[np.float64]]

Element-quality statistics.

Parameters:

Name	Type	Description	Default
p	(array_like, shape(n, 2))	Nodal coordinates.	required
t	(array_like, shape(m, 3) or (m, 4))	Connectivity list (0-based, per the port's indexing convention).	required
element	('triangle', 'quad')	Element type.	"triangle"

Returns:

Name	Type	Description
min_q	float
mean_q	float
q	(ndarray, shape(m))	Per-element quality in [0, 1].

Source code in src/admesh/_stages/quality.py

def mesh_quality(
    p: ArrayLike,
    t: ArrayLike,
    element: str = "triangle",
) -> tuple[float, float, NDArray[np.float64]]:
    """Element-quality statistics.

    Parameters
    ----------
    p : array_like, shape (n, 2)
        Nodal coordinates.
    t : array_like, shape (m, 3) or (m, 4)
        Connectivity list (0-based, per the port's indexing convention).
    element : {"triangle", "quad"}
        Element type.

    Returns
    -------
    min_q : float
    mean_q : float
    q : ndarray, shape (m,)
        Per-element quality in [0, 1].
    """
    p = np.asarray(p, dtype=float)
    t = np.asarray(t, dtype=int)

    if p.ndim != 2 or p.shape[1] != 2:
        raise ValueError("p must be (n, 2)")

    x = p[:, 0]
    y = p[:, 1]

    if element == "triangle":
        if t.shape[1] != 3:
            raise ValueError("triangle connectivity must have 3 columns")
        a = np.hypot(x[t[:, 1]] - x[t[:, 0]], y[t[:, 1]] - y[t[:, 0]])
        b = np.hypot(x[t[:, 2]] - x[t[:, 1]], y[t[:, 2]] - y[t[:, 1]])
        c = np.hypot(x[t[:, 0]] - x[t[:, 2]], y[t[:, 0]] - y[t[:, 2]])
        denom = a * b * c
        with np.errstate(divide="ignore", invalid="ignore"):
            q = np.where(
                denom > 0,
                ((b + c - a) * (c + a - b) * (a + b - c)) / denom,
                0.0,
            )
    elif element == "quad":
        if t.shape[1] != 4:
            raise ValueError("quad connectivity must have 4 columns")
        fx = np.column_stack([x[t[:, (k + 1) % 4]] - x[t[:, k]] for k in range(4)])
        fy = np.column_stack([y[t[:, (k + 1) % 4]] - y[t[:, k]] for k in range(4)])
        ax = fx
        ay = fy
        bx = np.roll(fx, -1, axis=1)
        by = np.roll(fy, -1, axis=1)
        dot = ax * bx + ay * by
        na = np.hypot(ax, ay)
        nb = np.hypot(bx, by)
        with np.errstate(divide="ignore", invalid="ignore"):
            cos_theta = np.clip(dot / (na * nb), -1.0, 1.0)
        theta = np.arccos(cos_theta)
        per_angle = 1.0 - np.abs((np.pi / 2 - theta) / (np.pi / 2))
        q = np.prod(per_angle, axis=1)
    else:
        raise ValueError(f"unknown element type: {element!r}")

    q = np.clip(q, 0.0, 1.0)
    return float(q.min()), float(q.mean()), q

right_iso_quality

admesh.right_iso_quality

right_iso_quality(p: ArrayLike, t: ArrayLike) -> float

Mesh-wide right-isoceles quality score in [0, 1].

Companion to :func:mesh_quality (which scores deviation from equilateral). Reported side-by-side as a delta after running :func:admesh.quad_prep.smooth_for_quadrangulation (spec-004 SC-007).

Per-element score is the product of three terms in [0, 1]:

1. Leg-equality: 1 - |L1 - L2| / max(L1, L2)
2. Right-angle: 1 - |angle_apex - π/2| / (π/2)
3. Hypotenuse-fit: 1 - |L_hyp - sqrt(2) * (L1+L2)/2| / L_hyp

where L1, L2 are the two shortest sides (legs), L_hyp is the longest (hypotenuse), and angle_apex is the angle between the two legs. The mesh score is the unweighted mean over elements.

The existing :func:mesh_quality is NOT modified — this is purely additive (spec-004 FR-006).

Parameters:

Name	Type	Description	Default
p	(array_like, shape(N, 2))	Nodal coordinates.	required
t	(array_like, shape(M, 3))	Triangle connectivity (0-based).	required

Returns:

Type	Description
float	Mesh-wide right-isoceles quality, in [0, 1]. 1.0 means every element is exactly right-isoceles. Empty mesh returns 1.0 (vacuously perfect).

Source code in src/admesh/_stages/quality.py

def right_iso_quality(
    p: ArrayLike,
    t: ArrayLike,
) -> float:
    """Mesh-wide right-isoceles quality score in ``[0, 1]``.

    Companion to :func:`mesh_quality` (which scores deviation from
    equilateral). Reported side-by-side as a delta after running
    :func:`admesh.quad_prep.smooth_for_quadrangulation` (spec-004 SC-007).

    Per-element score is the product of three terms in ``[0, 1]``:

    1. Leg-equality:    ``1 - |L1 - L2| / max(L1, L2)``
    2. Right-angle:     ``1 - |angle_apex - π/2| / (π/2)``
    3. Hypotenuse-fit:  ``1 - |L_hyp - sqrt(2) * (L1+L2)/2| / L_hyp``

    where ``L1, L2`` are the two shortest sides (legs), ``L_hyp`` is the
    longest (hypotenuse), and ``angle_apex`` is the angle between the
    two legs. The mesh score is the unweighted mean over elements.

    The existing :func:`mesh_quality` is NOT modified — this is purely
    additive (spec-004 FR-006).

    Parameters
    ----------
    p : array_like, shape (N, 2)
        Nodal coordinates.
    t : array_like, shape (M, 3)
        Triangle connectivity (0-based).

    Returns
    -------
    float
        Mesh-wide right-isoceles quality, in ``[0, 1]``. ``1.0`` means
        every element is exactly right-isoceles. Empty mesh returns
        ``1.0`` (vacuously perfect).
    """
    p = np.asarray(p, dtype=float)
    t = np.asarray(t, dtype=int)

    if p.ndim != 2 or p.shape[1] != 2:
        raise ValueError("p must be (N, 2)")
    if t.ndim != 2 or t.shape[1] != 3:
        raise ValueError("t must be (M, 3)")

    if len(t) == 0:
        return 1.0

    x = p[:, 0]
    y = p[:, 1]

    a = np.hypot(x[t[:, 1]] - x[t[:, 0]], y[t[:, 1]] - y[t[:, 0]])
    b = np.hypot(x[t[:, 2]] - x[t[:, 1]], y[t[:, 2]] - y[t[:, 1]])
    c = np.hypot(x[t[:, 0]] - x[t[:, 2]], y[t[:, 0]] - y[t[:, 2]])

    sides = np.column_stack([a, b, c])
    sides_sorted = np.sort(sides, axis=1)
    L1 = sides_sorted[:, 0]
    L2 = sides_sorted[:, 1]
    L_hyp = sides_sorted[:, 2]

    with np.errstate(divide="ignore", invalid="ignore"):
        leg_eq = np.where(
            L2 > 0, 1.0 - np.abs(L1 - L2) / np.maximum(L2, 1e-300), 0.0
        )

    cos_apex = (L1 * L1 + L2 * L2 - L_hyp * L_hyp) / (2.0 * L1 * L2)
    cos_apex = np.where(np.isfinite(cos_apex), cos_apex, 1.0)
    cos_apex = np.clip(cos_apex, -1.0, 1.0)
    angle_apex = np.arccos(cos_apex)
    right_angle = 1.0 - np.abs(angle_apex - np.pi / 2.0) / (np.pi / 2.0)

    target_hyp = np.sqrt(2.0) * (L1 + L2) / 2.0
    with np.errstate(divide="ignore", invalid="ignore"):
        hyp_fit = np.where(
            L_hyp > 0, 1.0 - np.abs(L_hyp - target_hyp) / L_hyp, 0.0
        )

    leg_eq = np.clip(leg_eq, 0.0, 1.0)
    right_angle = np.clip(right_angle, 0.0, 1.0)
    hyp_fit = np.clip(hyp_fit, 0.0, 1.0)

    q = leg_eq * right_angle * hyp_fit
    q = np.where(np.isfinite(q), q, 0.0)
    q = np.clip(q, 0.0, 1.0)

    return float(q.mean())

smooth_for_quadrangulation

admesh.smooth_for_quadrangulation

smooth_for_quadrangulation(p: NDArray[float64], t: NDArray[int64], fd: Callable[[NDArray[float64]], NDArray[float64]], h: Callable[[NDArray[float64]], NDArray[float64]] | None = None, pair_hint: bool = True, n_outer: int = 2) -> tuple[NDArray[np.float64], NDArray[np.int64]]

Nudge a triangle mesh toward right-isoceles for downstream quad fusion.

Implements the SVD-invariant FEM target-Jacobian formulation (Formulation 1 in specs/004-quad-prep-smoother/research.md). Connectivity is preserved (FR-002): t_out is the same array object as t.

Parameters:

Name	Type	Description	Default
p	ndarray, shape (N, 2), dtype float64	Node coordinates of the input triangulation.	required
t	ndarray, shape (M, 3), dtype int64	Triangle index triples; CCW winding.	required
fd	callable	Signed-distance function fd(q) -> distances for q of shape (K, 2). Required (FR-013).	required
h	callable	Size field h(q) -> edge_lengths. When provided, the per- element target leg length tracks h(centroid) (FR-004). Default: uniform target.	None
pair_hint	bool	When True, run a greedy longest-edge pairing pre-pass and bias geometry toward mutual longest-edge alignment via a soft per-element stiffness penalty (FR-005).	True
n_outer	int	Number of outer iterations. Each outer pass does one solve plus one boundary projection. Must be >= 1.	2

Returns:

Name	Type	Description
p_out	ndarray, shape (N, 2), dtype float64
t_out	ndarray, shape (M, 3), dtype int64	Same array object as input t.

Raises:

Type	Description
ValueError	If fd is None (FR-013), or input shapes are invalid, or n_outer < 1.

Source code in src/admesh/quad_prep.py

def smooth_for_quadrangulation(
    p: NDArray[np.float64],
    t: NDArray[np.int64],
    fd: Callable[[NDArray[np.float64]], NDArray[np.float64]],
    h: Callable[[NDArray[np.float64]], NDArray[np.float64]] | None = None,
    pair_hint: bool = True,
    n_outer: int = 2,
) -> tuple[NDArray[np.float64], NDArray[np.int64]]:
    """Nudge a triangle mesh toward right-isoceles for downstream quad fusion.

    Implements the SVD-invariant FEM target-Jacobian formulation
    (Formulation 1 in ``specs/004-quad-prep-smoother/research.md``).
    Connectivity is preserved (FR-002): ``t_out`` is the same array
    object as ``t``.

    Parameters
    ----------
    p : ndarray, shape (N, 2), dtype float64
        Node coordinates of the input triangulation.
    t : ndarray, shape (M, 3), dtype int64
        Triangle index triples; CCW winding.
    fd : callable
        Signed-distance function ``fd(q) -> distances`` for ``q`` of
        shape ``(K, 2)``. **Required** (FR-013).
    h : callable, optional
        Size field ``h(q) -> edge_lengths``. When provided, the per-
        element target leg length tracks ``h(centroid)`` (FR-004).
        Default: uniform target.
    pair_hint : bool, default True
        When ``True``, run a greedy longest-edge pairing pre-pass and
        bias geometry toward mutual longest-edge alignment via a soft
        per-element stiffness penalty (FR-005).
    n_outer : int, default 2
        Number of outer iterations. Each outer pass does one solve plus
        one boundary projection. Must be ``>= 1``.

    Returns
    -------
    p_out : ndarray, shape (N, 2), dtype float64
    t_out : ndarray, shape (M, 3), dtype int64
        Same array object as input ``t``.

    Raises
    ------
    ValueError
        If ``fd is None`` (FR-013), or input shapes are invalid, or
        ``n_outer < 1``.
    """
    if fd is None:
        raise ValueError("fd is required; pass an SDF callable")

    p = np.asarray(p, dtype=np.float64)
    if p.ndim != 2 or p.shape[1] != 2:
        raise ValueError(f"p must be (N, 2); got shape {p.shape}")
    if t.ndim != 2 or t.shape[1] != 3:
        raise ValueError(f"t must be (M, 3); got shape {t.shape}")
    if int(n_outer) < 1:
        raise ValueError(f"n_outer must be >= 1; got {n_outer}")

    # Empty / single-element fast path (Edge Case in spec.md).
    if len(t) == 0 or len(t) == 1:
        return p.copy(), t

    # geps follows the distmesh convention: 1e-3 * median edge length.
    # Computed from the input mesh once and held fixed across outer
    # iterations.
    geps = _auto_geps(p, t)

    # Boundary node mask is determined by the *input* mesh (FR-003) and
    # held fixed across outer iterations — topology and the boundary
    # ring don't change.
    boundary_mask = _boundary_node_mask(p, fd, geps=geps)

    # Build pairing map once (topology is invariant).
    pairs = _build_pairing_map(p, t) if pair_hint else None

    p_cur = p.copy()
    for _ in range(int(n_outer)):
        p_cur = _smoother_step(
            p_cur, t, fd=fd, h=h, pairs=pairs, boundary_mask=boundary_mask
        )
        p_cur = _project_boundary_nodes(p_cur, fd, geps=geps, mask=boundary_mask)

    return p_cur, t