energies

triwarp.energies

Energies over a mesh: the scalar regularizers, and the operators the solvers minimize.

The line between this module and triwarp.laplacian is order: laplacian builds the first-order operators -- the cotangent stiffness matrix, its mass matrix and the intrinsic repairs that keep them finite -- and everything here is assembled out of those.

Four families, and the first is the only one that returns a number rather than a matrix:

• Scalar mesh regularizers. edge_length_loss, normal_consistency_loss and laplacian_smoothing_loss are the priors a mesh optimization adds to a data term -- one per edge length, one per dihedral angle and one per vertex Laplacian residual, each reduced to a single float. They are here rather than in triwarp.metrics, which hosts the data terms they pair with, because their machinery is this module's: laplacian_smoothing_loss's two cotangent variants consume cotmatrix and mass_matrix_entries, which is exactly the import set the quadratic forms below use.

• Smoothness energies over vertices. k_harmonic is the integrated k-harmonic form -- k = 1 is Dirichlet, k = 2 the biharmonic operator behind smooth interpolation -- and it distorts a field near the boundary, because clamping a biharmonic solve there is not a natural condition. hessian_energy and curved_hessian_energy are the alternatives that do not: both integrate a squared Hessian instead, so linear (respectively, locally linear) functions sit exactly in the null space, boundary or not.

• The edge-based Crouzeix-Raviart pair. crouzeix_raviart_cotmatrix and crouzeix_raviart_massmatrix put the degrees of freedom on edge midpoints rather than vertices, which is the nonconforming-FEM discretization curved_hessian_energy is built on. They are the siblings of cotmatrix and mass_matrix -- a reader arriving from libigl, where all four sit together, should start there.
• The LSCM operator. lscm_hessian is the (2n, 2n) form behind the least-squares conformal map, and vector_area_matrix is the boundary term that couples its two coordinate blocks. lscm is the solve; these are what it minimizes.

Every operator is assembled in float64 in a single bsr_from_triplets, because the conjugate-gradient solves they feed run in float64 for determinism and a float32 intermediate would be the accuracy floor. The three scalar regularizers are float32 throughout: they are read by a human or by an optimizer's stopping rule, not solved with.

edge_length_loss

edge_length_loss(
    vertices: array[vec3],
    faces: array[int32],
    target_length: float = 0.0,
) -> float

Mean squared deviation of the undirected edge lengths from a resting length.

The edge regularizer of the deformation losses: mean((||e|| - L0)^2) over the unique undirected edges, so an interior edge counts once rather than twice. At the default target_length = 0.0 it is the mean squared edge length, which is what a shrinking prior wants; a positive value pulls the mesh toward uniform edges of that size instead.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) positions.	required
faces	array[int32]	Length-3 * n_faces wp.int32 triangle index buffer.	required
target_length	float	Resting edge length L0.	0.0

Returns:

Type	Description
float	The mean, as a host scalar. 0.0 for a mesh with no edges.

Raises:

Type	Description
RuntimeError	If vertices and faces are not all on one device.

See Also

normal_consistency_loss laplacian_smoothing_loss edges_unique_length The per-edge lengths this reduces, if the distribution rather than the mean is wanted. triwarp.metrics The data terms these regularizers are added to, and the differentiable Chamfer family.

Notes

Matches pytorch3d.loss.mesh_edge_loss. Its per-mesh 1 / E weighting collapses to a plain mean for one mesh, which is triwarp's only case, so there is no batch weighting to port.

Source code in triwarp/energies.py

def edge_length_loss(
    vertices: wp.array[wp.vec3], faces: wp.array[wp.int32], target_length: float = 0.0
) -> float:
    """
    Mean squared deviation of the undirected edge lengths from a resting length.

    The edge regularizer of the deformation losses: ``mean((||e|| - L0)^2)`` over the **unique
    undirected** edges, so an interior edge counts once rather than twice. At the default
    ``target_length = 0.0`` it is the mean squared edge length, which is what a shrinking prior
    wants; a positive value pulls the mesh toward uniform edges of that size instead.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    target_length
        Resting edge length ``L0``.

    Returns
    -------
    float
        The mean, as a host scalar. ``0.0`` for a mesh with no edges.

    Raises
    ------
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`normal_consistency_loss`][triwarp.energies.normal_consistency_loss]
    [`laplacian_smoothing_loss`][triwarp.energies.laplacian_smoothing_loss]
    [`edges_unique_length`][triwarp.edges.edges_unique_length]
        The per-edge lengths this reduces, if the distribution rather than the mean is wanted.
    [`triwarp.metrics`][triwarp.metrics]
        The data terms these regularizers are added to, and the differentiable Chamfer family.

    Notes
    -----
    Matches ``pytorch3d.loss.mesh_edge_loss``. Its per-mesh ``1 / E`` weighting collapses to a
    plain mean for one mesh, which is triwarp's only case, so there is no batch weighting to port.
    """
    require_same_device(vertices=vertices, faces=faces)
    lengths = edges_unique_length(vertices, faces)
    if int(lengths.shape[0]) == 0:
        return 0.0
    deviations = wp.empty(int(lengths.shape[0]), dtype=wp.float32, device=lengths.device)
    wp.map(kernel_energies.squared_deviation, lengths, wp.float32(target_length), out=deviations)
    return float(tw.reduce.mean(deviations))

normal_consistency_loss

normal_consistency_loss(
    vertices: array[vec3], faces: array[int32]
) -> float

Mean 1 - cos(theta) over the pairs of faces sharing an edge.

The dihedral regularizer of the deformation losses, and the one that penalizes a fold: it is 0 for a flat pair, 1 at a right angle and 2 for a face doubled back on itself. Read it against edge_length_loss, which constrains the sizes; this one constrains the orientations and says nothing about scale.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) positions.	required
faces	array[int32]	Length-3 * n_faces wp.int32 triangle index buffer.	required

Returns:

Type	Description
float	The mean, as a host scalar. 0.0 for a mesh with no adjacent face pair.

Raises:

Type	Description
RuntimeError	If vertices and faces are not all on one device.

See Also

edge_length_loss laplacian_smoothing_loss face_adjacency_angles The per-pair dihedral angles this reduces.

Notes

Matches pytorch3d.loss.mesh_normal_consistency on edge-manifold input. The restriction is real and is not a tolerance -- the reference enumerates every pair of faces sharing an edge, so an edge with k incident faces contributes C(k, 2) terms where face_adjacency_angles reports one pair per adjacency. The two coincide exactly wherever every edge has at most two faces.

Source code in triwarp/energies.py

def normal_consistency_loss(vertices: wp.array[wp.vec3], faces: wp.array[wp.int32]) -> float:
    """
    Mean ``1 - cos(theta)`` over the pairs of faces sharing an edge.

    The dihedral regularizer of the deformation losses, and the one that penalizes a fold: it is
    ``0`` for a flat pair, ``1`` at a right angle and ``2`` for a face doubled back on itself. Read
    it against [`edge_length_loss`][triwarp.energies.edge_length_loss], which constrains the
    *sizes*; this one constrains the *orientations* and says nothing about scale.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.

    Returns
    -------
    float
        The mean, as a host scalar. ``0.0`` for a mesh with no adjacent face pair.

    Raises
    ------
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`edge_length_loss`][triwarp.energies.edge_length_loss]
    [`laplacian_smoothing_loss`][triwarp.energies.laplacian_smoothing_loss]
    [`face_adjacency_angles`][triwarp.adjacency.face_adjacency_angles]
        The per-pair dihedral angles this reduces.

    Notes
    -----
    Matches ``pytorch3d.loss.mesh_normal_consistency`` on **edge-manifold** input. The restriction
    is real and is not a tolerance -- the reference enumerates *every* pair of faces sharing an
    edge, so an edge with ``k`` incident faces contributes ``C(k, 2)`` terms where
    [`face_adjacency_angles`][triwarp.adjacency.face_adjacency_angles] reports one pair per
    adjacency. The two coincide exactly wherever every edge has at most two faces.
    """
    require_same_device(vertices=vertices, faces=faces)
    angles = tw.adjacency.face_adjacency_angles(vertices, faces)
    if int(angles.shape[0]) == 0:
        return 0.0
    terms = wp.empty(int(angles.shape[0]), dtype=wp.float32, device=angles.device)
    wp.map(kernel_energies.one_minus_cosine, angles, out=terms)
    return float(tw.reduce.mean(terms))

laplacian_smoothing_loss

laplacian_smoothing_loss(
    vertices: array[vec3],
    faces: array[int32],
    method: Literal[
        "uniform", "cot", "cotcurv"
    ] = "uniform",
) -> float

Mean magnitude of a Laplacian residual per vertex, under one of three normalizations.

The smoothness regularizer of the deformation losses. The three methods are three different quantities, not one with a tuning knob, and they differ by roughly an order of magnitude:

• "uniform": || (A v)_i - v_i || with A the row-normalized 1-ring average (laplacian with equal_weight=True) -- the umbrella residual, which is a length and therefore scales with the mesh.
• "cot": the same residual against the cotangent-weighted neighbour average, || (L v)_i / s_i || with L the cotangent stiffness matrix and s_i its off-diagonal row sum. Also a length, and the geometry-aware version of the above.
• "cotcurv": || (L v)_i / (6 M_ii) || with M the barycentric lumped mass -- the mean curvature magnitude, so it carries units of one over length and is the largest of the three on a unit-scale mesh.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) positions.	required
faces	array[int32]	Length-3 * n_faces wp.int32 triangle index buffer.	required
method	Literal['uniform', 'cot', 'cotcurv']	Which normalization, as above.	'uniform'

Returns:

Type	Description
float	The mean, as a host scalar. 0.0 for an empty mesh.

Raises:

Type	Description
ValueError	If method is not one of "uniform", "cot" or "cotcurv".
RuntimeError	If vertices and faces are not all on one device.

See Also

edge_length_loss normal_consistency_loss cotmatrix The operator the two cotangent variants are built from. filter_laplacian The smoother that minimizes this, rather than measuring it.

Notes

Matches pytorch3d.loss.mesh_laplacian_smoothing on all three methods. Two conventions are inherited from it rather than chosen here, both because they are what makes the numbers comparable at all: the reference's cot and cotcurv read a cotangent Laplacian whose off-diagonal is twice triwarp's half-cotangent table and whose diagonal is identically zero, which cancels out of both ratios above -- and where a vertex's row sum is not positive its averaging is undefined, so "cot" falls back to || v_i || there, matching the reference's norm_w = 0 branch.

Source code in triwarp/energies.py

def laplacian_smoothing_loss(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    method: Literal["uniform", "cot", "cotcurv"] = "uniform",
) -> float:
    """
    Mean magnitude of a Laplacian residual per vertex, under one of three normalizations.

    The smoothness regularizer of the deformation losses. The three methods are **three different
    quantities**, not one with a tuning knob, and they differ by roughly an order of magnitude:

    - ``"uniform"``: ``|| (A v)_i - v_i ||`` with ``A`` the row-normalized 1-ring average
      ([`laplacian`][triwarp.laplacian.laplacian] with ``equal_weight=True``) -- the umbrella
      residual, which is a *length* and therefore scales with the mesh.
    - ``"cot"``: the same residual against the **cotangent-weighted** neighbour average,
      ``|| (L v)_i / s_i ||`` with ``L`` the cotangent stiffness matrix and ``s_i`` its
      off-diagonal row sum. Also a length, and the geometry-aware version of the above.
    - ``"cotcurv"``: ``|| (L v)_i / (6 M_ii) ||`` with ``M`` the barycentric lumped mass -- the mean
      curvature magnitude, so it carries units of one over length and is the largest of the three
      on a unit-scale mesh.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    method
        Which normalization, as above.

    Returns
    -------
    float
        The mean, as a host scalar. ``0.0`` for an empty mesh.

    Raises
    ------
    ValueError
        If ``method`` is not one of ``"uniform"``, ``"cot"`` or ``"cotcurv"``.
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`edge_length_loss`][triwarp.energies.edge_length_loss]
    [`normal_consistency_loss`][triwarp.energies.normal_consistency_loss]
    [`cotmatrix`][triwarp.laplacian.cotmatrix]
        The operator the two cotangent variants are built from.
    [`filter_laplacian`][triwarp.smoothing.filter_laplacian]
        The smoother that *minimizes* this, rather than measuring it.

    Notes
    -----
    Matches ``pytorch3d.loss.mesh_laplacian_smoothing`` on all three methods. Two conventions are
    inherited from it rather than chosen here, both because they are what makes the numbers
    comparable at all: the
    reference's ``cot`` and ``cotcurv`` read a cotangent Laplacian whose off-diagonal is **twice**
    triwarp's half-cotangent table and whose diagonal is identically zero, which cancels out of
    both ratios above -- and where a vertex's row sum is not positive its averaging is undefined,
    so ``"cot"`` falls back to ``|| v_i ||`` there, matching the reference's ``norm_w = 0`` branch.
    """
    require_same_device(vertices=vertices, faces=faces)
    if method not in ("uniform", "cot", "cotcurv"):
        raise ValueError(f'method must be "uniform", "cot" or "cotcurv", got {method!r}')
    n_vertices = int(vertices.shape[0])
    if n_vertices == 0 or int(faces.shape[0]) == 0:
        return 0.0
    device = vertices.device
    # The three methods initialize the two scale buffers three different ways, so each branch
    # allocates them holding what it needs -- ``wp.full`` where the value is a constant,
    # ``wp.empty`` where a launch writes every element, ``wp.zeros`` where the value is zero.
    if method == "uniform":
        operator = tw.laplacian.laplacian(vertices, faces, equal_weight=True)
        row_scale = wp.full(n_vertices, 1.0, dtype=wp.float32, device=device)
        self_scale = wp.full(n_vertices, -1.0, dtype=wp.float32, device=device)
    else:
        operator = cotmatrix(vertices, faces)
        row_scale = wp.empty(n_vertices, dtype=wp.float32, device=device)
        if method == "cot":
            self_scale = wp.empty(n_vertices, dtype=wp.float32, device=device)
            wp.launch(
                kernel_energies.cot_row_scales,
                dim=n_vertices,
                inputs=[operator.offsets, operator.columns, operator.values, row_scale, self_scale],
                device=device,
            )
        else:
            wp.map(
                kernel_energies.reciprocal_scaled_or_zero,
                mass_matrix_entries(vertices, faces),
                wp.float32(1.0 / 6.0),
                out=row_scale,
            )
            self_scale = wp.zeros(n_vertices, dtype=wp.float32, device=device)
    norms = wp.empty(n_vertices, dtype=wp.float32, device=device)
    wp.launch(
        kernel_energies.laplacian_residual_norms,
        dim=n_vertices,
        inputs=[
            operator.offsets,
            operator.columns,
            operator.values,
            vertices,
            row_scale,
            self_scale,
            norms,
        ],
        device=device,
    )
    return float(tw.reduce.mean(norms))

k_harmonic

k_harmonic(
    laplacian: BsrMatrix[float32],
    mass: Array1dFloat | None = None,
    k: int = 2,
) -> wps.BsrMatrix[wp.float32]

Integrated k-harmonic operator Q = (-L) (M^-1 (-L))^(k-1) from a Laplacian and a mass.

The quadratic form whose minimizers are k-harmonic functions: k == 1 gives the Dirichlet energy -L (positive semi-definite for a cotmatrix-sign Laplacian), k == 2 the biharmonic operator L M^-1 L behind harmonic's smooth interpolation, and so on (igl::harmonic_integrated_from_laplacian_and_mass). Like igl's, the composition is not numerically robust for k > 2 — the entries grow as the k-th power of the inverse mesh size — so high powers want a float64 laplacian.

Named for the operator, not the map. This used to be laplacian.harmonic_integrated, which put a second unrelated "harmonic" in the package beside harmonic — a map into the plane, not an operator. The two are related (harmonic minimizes this form with the boundary pinned) but they are not interchangeable, and neither docstring named the other.

Each power is assembled by one triplet pass over matching CSR rows — (A M^-1 B)_ij = sum_t A_ti M_t^-1 B_tj with both operands symmetric — followed by a single bsr_from_triplets, so the product is built without warp.sparse.bsr_mm.

Parameters:

Name	Type	Description	Default
laplacian	BsrMatrix[float32]	Square 1x1-block BSR Laplacian in igl's sign convention (negative diagonal, each row summing to zero), e.g. from cotmatrix or graph_laplacian.	required
mass	Array1dFloat | None	Length-n_vertices lumped mass diagonal, e.g. from mass_matrix_entries; cast to the Laplacian's scalar type if it differs. None (identity mass) composes plain powers of -L, the tutte convention. Zero entries are treated as killed degrees of freedom (their rows contribute nothing), matching igl::invert_diag.	None
k	int	Harmonic power (>= 1): 1 harmonic, 2 biharmonic, 3 triharmonic, ...	2

Returns:

Type	Description
BsrMatrix	Square (n_vertices, n_vertices) positive semi-definite operator in the Laplacian's scalar type on its device.

Raises:

Type	Description
ValueError	If k < 1.
RuntimeError	If laplacian and mass are not on the same device.

See Also

cotmatrix mass_matrix_entries hessian_energy harmonic

Source code in triwarp/energies.py

def k_harmonic(
    laplacian: wps.BsrMatrix[wp.float32], mass: twt.Array1dFloat | None = None, k: int = 2
) -> wps.BsrMatrix[wp.float32]:
    """
    Integrated k-harmonic operator ``Q = (-L) (M^-1 (-L))^(k-1)`` from a Laplacian and a mass.

    The quadratic form whose minimizers are k-harmonic functions: ``k == 1`` gives the Dirichlet
    energy ``-L`` (positive semi-definite for a [`cotmatrix`][triwarp.laplacian.cotmatrix]-sign
    Laplacian), ``k == 2`` the biharmonic operator ``L M^-1 L`` behind
    [`harmonic`][triwarp.parametrization.harmonic]'s smooth interpolation, and so on
    (``igl::harmonic_integrated_from_laplacian_and_mass``). Like igl's, the composition is not
    numerically robust for ``k > 2`` — the entries grow as the k-th power of the inverse mesh
    size — so high powers want a float64 ``laplacian``.

    **Named for the operator, not the map.** This used to be ``laplacian.harmonic_integrated``,
    which put a second unrelated "harmonic" in the package beside
    [`harmonic`][triwarp.parametrization.harmonic] — a *map* into the plane, not an operator. The
    two are related (``harmonic`` minimizes this form with the boundary pinned) but they are not
    interchangeable, and neither docstring named the other.

    Each power is assembled by one triplet pass over matching CSR rows —
    ``(A M^-1 B)_ij = sum_t A_ti M_t^-1 B_tj`` with both operands symmetric — followed by a single
    ``bsr_from_triplets``, so the product is built without ``warp.sparse.bsr_mm``.

    Parameters
    ----------
    laplacian
        Square 1x1-block BSR Laplacian in igl's sign convention (negative diagonal, each row
        summing to zero), e.g. from [`cotmatrix`][triwarp.laplacian.cotmatrix] or
        [`graph_laplacian`][triwarp.laplacian.graph_laplacian].
    mass
        Length-``n_vertices`` lumped mass diagonal, e.g. from
        [`mass_matrix_entries`][triwarp.laplacian.mass_matrix_entries]; cast to the Laplacian's
        scalar type if it differs. ``None`` (identity mass) composes plain powers of ``-L``, the
        [`tutte`][triwarp.parametrization.tutte] convention. Zero entries are treated as killed
        degrees of freedom (their rows contribute nothing), matching ``igl::invert_diag``.
    k
        Harmonic power (``>= 1``): 1 harmonic, 2 biharmonic, 3 triharmonic, ...

    Returns
    -------
    warp.sparse.BsrMatrix
        Square ``(n_vertices, n_vertices)`` positive semi-definite operator in the Laplacian's
        scalar type on its device.

    Raises
    ------
    ValueError
        If ``k < 1``.
    RuntimeError
        If ``laplacian`` and ``mass`` are not on the same device.

    See Also
    --------
    [`cotmatrix`][triwarp.laplacian.cotmatrix]
    [`mass_matrix_entries`][triwarp.laplacian.mass_matrix_entries]
    [`hessian_energy`][triwarp.energies.hessian_energy]
    [`harmonic`][triwarp.parametrization.harmonic]
    """
    require_same_device(laplacian=laplacian, mass=mass)
    if k < 1:
        raise ValueError(f"harmonic power k must be >= 1, got {k}.")
    negated = wps.bsr_axpy(x=laplacian, alpha=-1.0)
    if k == 1:
        return negated

    dtype = laplacian.values.dtype
    n_rows = int(laplacian.nrow)
    device = laplacian.values.device
    if mass is None:
        inverse_mass = wp.ones(n_rows, dtype=dtype, device=device)
    else:
        if mass.dtype != dtype:
            mass = tw.array.astype(mass, dtype)
        inverse_mass = wp.empty(n_rows, dtype=dtype, device=device)
        wp.map(kernel_energies.reciprocal_or_zero, mass, out=inverse_mass)

    operator = negated
    for _ in range(k - 1):
        operator = _diagonal_sandwich(operator, inverse_mass, negated)
    return operator

hessian_energy

hessian_energy(
    vertices: array[vec3],
    faces: array[int32],
    dtype: type = wp.float64,
    *,
    vertex_faces: tuple[array[int32], array[int32]]
    | None = None,
) -> wps.BsrMatrix[wp.float64]

Hessian smoothness energy with natural boundary conditions.

The mixed-FEM quadratic form Q = H^T M^-1 H of Stein et al. 2018, Natural Boundary Conditions for Smoothing in Geometry Processing: x' Q x integrates the squared Hessian of the piecewise-linear field x, so minimizing it smooths without the boundary distortion the clamped biharmonic operator (k_harmonic at k == 2) produces — linear functions are exactly in its null space, boundary or not. M is the Voronoi lumped mass with boundary degrees of freedom killed, per the reference.

Rather than materializing the sparse (9 n_faces, n_vertices) stacked Hessian H, the product is contracted analytically over its nine component pairs and assembled in one triplet pass per vertex: Q_ij = sum_k M_k^-1 sum_{f,g ni k} A_f A_g (g_fk . g_gk)(g_fi . g_gj) with g_fc corner c's hat-function gradient in face f. The per-vertex triplet count is 9 * valence^2, so cost is quadratic in valence. A degenerate face contributes nothing (igl emits NaN there).

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) mesh vertex positions.	required
faces	array[int32]	Length-3 * n_faces wp.int32 triangle index buffer.	required
dtype	type	Scalar block type of the assembled matrix: wp.float64 (default) or wp.float32. Float64 is the default, unlike the first-order operators in this module, because the entries scale as the inverse fourth power of the mesh size and the operator exists to be solved against.	float64
vertex_faces	tuple[array[int32], array[int32]] | None	Optional precomputed vertex_face_adjacency as (vertex_faces, offsets). Depends on the connectivity alone, so one CSR serves every incidence walk over the same mesh -- Trimesh.vertex_face_adjacency has it cached.	None

Returns:

Type	Description
BsrMatrix	Square (n_vertices, n_vertices) positive semi-definite energy matrix on vertices.device.

Raises:

Type	Description
RuntimeError	If vertices, faces and vertex_faces are not all on one device.

See Also

curved_hessian_energy k_harmonic cotmatrix

Notes

Matches igl::hessian_energy except on degenerate faces, which contribute nothing here and NaN there.

Source code in triwarp/energies.py

def hessian_energy(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    dtype: type = wp.float64,
    *,
    vertex_faces: tuple[wp.array[wp.int32], wp.array[wp.int32]] | None = None,
) -> wps.BsrMatrix[wp.float64]:
    """
    Hessian smoothness energy with natural boundary conditions.

    The mixed-FEM quadratic form ``Q = H^T M^-1 H`` of Stein et al. 2018, *Natural Boundary
    Conditions for Smoothing in Geometry Processing*: ``x' Q x`` integrates the squared Hessian of
    the piecewise-linear field ``x``, so minimizing it smooths **without** the boundary distortion
    the clamped biharmonic operator
    ([`k_harmonic`][triwarp.energies.k_harmonic] at ``k == 2``) produces —
    linear functions are exactly in its null space, boundary or not. ``M`` is the Voronoi lumped
    mass with boundary degrees of freedom killed, per the reference.

    Rather than materializing the sparse ``(9 n_faces, n_vertices)`` stacked Hessian ``H``, the
    product is contracted analytically over its nine component pairs and assembled in one triplet
    pass per vertex: ``Q_ij = sum_k M_k^-1 sum_{f,g ni k} A_f A_g (g_fk . g_gk)(g_fi . g_gj)``
    with ``g_fc`` corner ``c``'s hat-function gradient in face ``f``. The per-vertex triplet count
    is ``9 * valence^2``, so cost is quadratic in valence. A degenerate face contributes nothing
    (igl emits NaN there).

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    dtype
        Scalar block type of the assembled matrix: ``wp.float64`` (default) or ``wp.float32``.
        Float64 is the default, unlike the first-order operators in this module, because the
        entries scale as the inverse fourth power of the mesh size and the operator exists to be
        solved against.
    vertex_faces
        Optional precomputed [`vertex_face_adjacency`][triwarp.adjacency.vertex_face_adjacency] as
        ``(vertex_faces, offsets)``. Depends on the connectivity alone, so one CSR serves every
        incidence walk over the same mesh --
        [`Trimesh.vertex_face_adjacency`][triwarp.mesh.Trimesh.vertex_face_adjacency] has it cached.

    Returns
    -------
    warp.sparse.BsrMatrix
        Square ``(n_vertices, n_vertices)`` positive semi-definite energy matrix on
        ``vertices.device``.

    Raises
    ------
    RuntimeError
        If ``vertices``, ``faces`` and ``vertex_faces`` are not all on one device.

    See Also
    --------
    [`curved_hessian_energy`][triwarp.energies.curved_hessian_energy]
    [`k_harmonic`][triwarp.energies.k_harmonic]
    [`cotmatrix`][triwarp.laplacian.cotmatrix]

    Notes
    -----
    Matches ``igl::hessian_energy`` except on degenerate faces, which contribute nothing here and
    ``NaN`` there.
    """
    require_same_device(vertices=vertices, faces=faces, vertex_faces=vertex_faces)
    n_vertices = int(vertices.shape[0])
    n_faces = int(faces.shape[0]) // 3
    device = vertices.device
    if n_faces == 0:
        return tw.array.empty_square_bsr(n_vertices, dtype, device)

    gradients = wp.empty(3 * n_faces, dtype=wp.vec3d, device=device)
    areas = wp.empty(n_faces, dtype=wp.float64, device=device)
    mass = wp.zeros(n_vertices, dtype=wp.float64, device=device)
    wp.launch(
        kernel_energies.hessian_face_terms,
        dim=n_faces,
        inputs=[vertices, faces, gradients, areas, mass],
        device=device,
    )
    inverse_mass = _interior_inverse(vertices, faces, mass)

    vf_indices, vf_offsets = (
        vertex_faces
        if vertex_faces is not None
        else tw.adjacency.vertex_face_adjacency(faces, n_vertices=n_vertices)
    )
    counts = wp.empty(n_vertices, dtype=wp.int32, device=device)
    wp.launch(
        kernel_energies.hessian_energy_counts,
        dim=n_vertices,
        inputs=[vf_offsets, inverse_mass, counts],
        device=device,
    )
    # Host readback: only the device knows the scan total, and it sizes the triplet buffers.
    segment_offsets, n_triplets = tw.array.counts_to_offsets(counts, include_total=True)

    rows, cols, vals = tw.array.triplet_buffers(n_triplets, dtype, device)
    if n_triplets > 0:
        wp.launch(
            kernel_energies.HESSIAN_ENERGY_TRIPLETS[dtype],
            dim=n_vertices,
            inputs=[
                faces,
                vf_offsets,
                vf_indices,
                gradients,
                areas,
                inverse_mass,
                segment_offsets,
                rows,
                cols,
                vals,
            ],
            device=device,
        )
    return wps.bsr_from_triplets(
        n_vertices, n_vertices, rows, cols, vals, prune_numerical_zeros=False
    )

curved_hessian_energy

curved_hessian_energy(
    vertices: array[vec3],
    faces: array[int32],
    dtype: type = wp.float64,
) -> wps.BsrMatrix[wp.float64]

Curved Hessian smoothness energy on the Crouzeix-Raviart discretization.

igl::curved_hessian_energy, from Stein et al. 2020, A Smoothness Energy without Boundary Distortion for Curved Surfaces: where hessian_energy treats the surface as locally flat, this one carries the Gaussian curvature into the operator through a per-vertex angle-defect correction, so the energy is intrinsic to the curved surface rather than to its triangles' planes. Constant functions are exactly in its null space.

Assembled as Q = D^T M^-1 (L + K) M^-1 D over edge-based Crouzeix-Raviart vector elements — D the scalar-to-CR-vector gradient, M the CR vector mass, L the CR vector Laplacian and K the curvature correction — but contracted per face in one pass: L + K couples edges within a face only, so each face emits its own 6x6 block sandwiched between its edges' gradient rows (a fixed 144 triplets per face), and no intermediate (2 n_edges, ...) matrix or sparse product exists. Requires an edge-manifold mesh, like the igl original (which asserts it) — raises ValueError otherwise; a degenerate face contributes nothing.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) mesh vertex positions.	required
faces	array[int32]	Length-3 * n_faces wp.int32 triangle index buffer.	required
dtype	type	Scalar block type of the assembled matrix: wp.float64 (default) or wp.float32, with float64 the default for the same conditioning reason as hessian_energy.	float64

Returns:

Type	Description
BsrMatrix	Square (n_vertices, n_vertices) positive semi-definite energy matrix on vertices.device.

Raises:

Type	Description
ValueError	If the mesh is not edge-manifold (some edge is shared by more than two faces).
RuntimeError	If vertices and faces are not all on one device.

See Also

hessian_energy crouzeix_raviart_cotmatrix vertex_defects

Source code in triwarp/energies.py

def curved_hessian_energy(
    vertices: wp.array[wp.vec3], faces: wp.array[wp.int32], dtype: type = wp.float64
) -> wps.BsrMatrix[wp.float64]:
    """
    Curved Hessian smoothness energy on the Crouzeix-Raviart discretization.

    ``igl::curved_hessian_energy``, from Stein et al. 2020, *A Smoothness Energy without Boundary
    Distortion for Curved Surfaces*: where [`hessian_energy`][triwarp.energies.hessian_energy]
    treats the surface as locally flat, this one carries the Gaussian curvature into the operator
    through a per-vertex angle-defect correction, so the energy is intrinsic to the curved surface
    rather than to its triangles' planes. Constant functions are exactly in its null space.

    Assembled as ``Q = D^T M^-1 (L + K) M^-1 D`` over edge-based Crouzeix-Raviart vector elements
    — ``D`` the scalar-to-CR-vector gradient, ``M`` the CR vector mass, ``L`` the CR vector
    Laplacian and ``K`` the curvature correction — but contracted per face in one pass: ``L + K``
    couples edges within a face only, so each face emits its own 6x6 block sandwiched between its
    edges' gradient rows (a fixed 144 triplets per face), and no intermediate ``(2 n_edges, ...)``
    matrix or sparse product exists. Requires an edge-manifold mesh, like the igl original (which
    asserts it) — raises ``ValueError`` otherwise; a degenerate face contributes nothing.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    dtype
        Scalar block type of the assembled matrix: ``wp.float64`` (default) or ``wp.float32``,
        with float64 the default for the same conditioning reason as
        [`hessian_energy`][triwarp.energies.hessian_energy].

    Returns
    -------
    warp.sparse.BsrMatrix
        Square ``(n_vertices, n_vertices)`` positive semi-definite energy matrix on
        ``vertices.device``.

    Raises
    ------
    ValueError
        If the mesh is not edge-manifold (some edge is shared by more than two faces).
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`hessian_energy`][triwarp.energies.hessian_energy]
    [`crouzeix_raviart_cotmatrix`][triwarp.energies.crouzeix_raviart_cotmatrix]
    [`vertex_defects`][triwarp.vertices.vertex_defects]
    """
    require_same_device(vertices=vertices, faces=faces)
    n_vertices = int(vertices.shape[0])
    n_faces = int(faces.shape[0]) // 3
    device = vertices.device
    if n_faces == 0:
        return tw.array.empty_square_bsr(n_vertices, dtype, device)

    # Reused below for ``edges_unique`` too, so the manifold check costs no extra sort.
    edges_sorted = tw.edges.faces_to_edges(faces, sorted=True)
    if not tw.validation.is_edge_manifold(faces, edges_sorted=edges_sorted, n_vertices=n_vertices):
        raise ValueError(
            "mesh must be edge-manifold (every edge shared by at most two faces); the "
            "Crouzeix-Raviart discretization curved_hessian_energy is built on is undefined "
            "otherwise, like igl::curved_hessian_energy, which asserts it"
        )

    unique_edges, inverse = edges_unique(faces, edges_sorted, n_vertices=n_vertices)
    n_edges = int(unique_edges.shape[0])

    angles = wp.empty((n_faces, 3), dtype=wp.float64, device=device)
    angle_sums = wp.zeros(n_vertices, dtype=wp.float64, device=device)
    wp.launch(
        kernel_energies.internal_angles_and_sums,
        dim=n_faces,
        inputs=[vertices, faces, angles, angle_sums],
        device=device,
    )
    # Angle defect, zeroed on the boundary (curvature is only corrected at interior vertices),
    # weighted by the actual angle sum -- igl::cr_vector_curvature_correction's kappa scaling.
    kappa = wp.empty(n_vertices, dtype=wp.float64, device=device)
    wp.map(kernel_predicates.angle_defect, angle_sums, out=kappa)
    _zero_at_boundary(vertices, faces, kappa)
    scaled_kappa = wp.empty(n_vertices, dtype=wp.float64, device=device)
    # ``kappa / angle_sums``, zeroed instead of dividing by a non-positive angle sum --
    # ``reciprocal_scaled_or_zero(value, numerator)`` takes the denominator first.
    wp.map(kernel_energies.reciprocal_scaled_or_zero, angle_sums, kappa, out=scaled_kappa)

    mass = _cr_mass_diagonal(vertices, faces, inverse, n_edges, wp.float64)
    inverse_mass = wp.empty(n_edges, dtype=wp.float64, device=device)
    wp.map(kernel_energies.reciprocal_or_zero, mass, out=inverse_mass)

    edge_halfedges = wp.full((n_edges, 2), -1, dtype=wp.int32, device=device)
    cursor = wp.zeros(n_edges, dtype=wp.int32, device=device)
    wp.launch(
        kernel_energies.scatter_edge_halfedges,
        dim=3 * n_faces,
        inputs=[inverse, cursor, edge_halfedges],
        device=device,
    )
    vertex_slots = wp.full((n_edges, 4), -1, dtype=wp.int32, device=device)
    par = wp.zeros((n_edges, 4), dtype=wp.float64, device=device)
    perp = wp.zeros((n_edges, 4), dtype=wp.float64, device=device)
    wp.launch(
        kernel_energies.cr_gradient_rows,
        dim=n_edges,
        inputs=[vertices, faces, unique_edges, edge_halfedges, vertex_slots, par, perp],
        device=device,
    )

    n_triplets = 144 * n_faces
    rows, cols, vals = tw.array.triplet_buffers(n_triplets, dtype, device)
    wp.launch(
        kernel_energies.CURVED_HESSIAN_TRIPLETS[dtype],
        dim=n_faces,
        inputs=[
            vertices,
            faces,
            inverse,
            angles,
            scaled_kappa,
            inverse_mass,
            vertex_slots,
            par,
            perp,
            rows,
            cols,
            vals,
        ],
        device=device,
    )
    return wps.bsr_from_triplets(
        n_vertices, n_vertices, rows, cols, vals, prune_numerical_zeros=False
    )

crouzeix_raviart_cotmatrix

crouzeix_raviart_cotmatrix(
    vertices: array[vec3],
    faces: array[int32],
    cot_entries: Array2dFloat | None = None,
    dtype: type = wp.float32,
    *,
    unique_edges: Array2dInt32 | None = None,
    edge_map: array[int32] | None = None,
) -> wps.BsrMatrix[wp.float32]

Edge-based Crouzeix-Raviart cotangent stiffness matrix.

The nonconforming-FEM sibling of cotmatrix: degrees of freedom live on edge midpoints, so the matrix is (n_edges, n_edges) and each face couples its three edges pairwise with minus four times the half-cotangent at their shared corner (positive diagonal — the igl sign convention for this operator, opposite to cotmatrix's). Rows follow edges_unique's edge numbering. Requires an edge-manifold mesh, like the igl original (which asserts it) — raises ValueError otherwise.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) mesh vertex positions.	required
faces	array[int32]	Length-3 * n_faces wp.int32 triangle index buffer.	required
cot_entries	Array2dFloat | None	Optional precomputed (n_faces, 3) weights from cotmatrix_entries; computed in dtype when None.	None
dtype	type	Scalar block type of the assembled matrix: wp.float32 (default) or wp.float64.	float32
unique_edges	Array2dInt32 | None	Optional precomputed edge numbering from edges_unique — pass both or neither. Sharing it with crouzeix_raviart_massmatrix keeps the two operators on identical rows without recomputing the sort.	None
edge_map	Array2dInt32 | None	Optional precomputed edge numbering from edges_unique — pass both or neither. Sharing it with crouzeix_raviart_massmatrix keeps the two operators on identical rows without recomputing the sort.	None

Returns:

Type	Description
BsrMatrix	Square (n_edges, n_edges) stiffness matrix in 1x1-block BSR form on vertices.device.

Raises:

Type	Description
ValueError	If exactly one of unique_edges / edge_map is provided, or if the mesh is not edge-manifold (some edge is shared by more than two faces).
RuntimeError	If vertices, faces, cot_entries, unique_edges and edge_map are not all on one device.

See Also

crouzeix_raviart_massmatrix cotmatrix edges_unique

Notes

Matches igl::crouzeix_raviart_cotmatrix up to the edge numbering, which follows edges_unique rather than igl::unique_edge_map.

Source code in triwarp/energies.py

def crouzeix_raviart_cotmatrix(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    cot_entries: twt.Array2dFloat | None = None,
    dtype: type = wp.float32,
    *,
    unique_edges: twt.Array2dInt32 | None = None,
    edge_map: wp.array[wp.int32] | None = None,
) -> wps.BsrMatrix[wp.float32]:
    """
    Edge-based Crouzeix-Raviart cotangent stiffness matrix.

    The nonconforming-FEM sibling of [`cotmatrix`][triwarp.laplacian.cotmatrix]: degrees of
    freedom live on edge midpoints, so the matrix is ``(n_edges, n_edges)`` and each face couples
    its three edges pairwise with minus four times the half-cotangent at their shared corner
    (positive diagonal — the igl sign convention for this operator, opposite to ``cotmatrix``'s).
    Rows follow [`edges_unique`][triwarp.edges.edges_unique]'s edge numbering. Requires an
    edge-manifold mesh, like the igl original (which asserts it) — raises ``ValueError``
    otherwise.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    cot_entries
        Optional precomputed ``(n_faces, 3)`` weights from
        [`cotmatrix_entries`][triwarp.laplacian.cotmatrix_entries]; computed in ``dtype`` when
        ``None``.
    dtype
        Scalar block type of the assembled matrix: ``wp.float32`` (default) or ``wp.float64``.
    unique_edges, edge_map
        Optional precomputed edge numbering from
        [`edges_unique`][triwarp.edges.edges_unique] — pass both or neither. Sharing it with
        [`crouzeix_raviart_massmatrix`][triwarp.energies.crouzeix_raviart_massmatrix] keeps the
        two operators on identical rows without recomputing the sort.

    Returns
    -------
    warp.sparse.BsrMatrix
        Square ``(n_edges, n_edges)`` stiffness matrix in 1x1-block BSR form on
        ``vertices.device``.

    Raises
    ------
    ValueError
        If exactly one of ``unique_edges`` / ``edge_map`` is provided, or if the mesh is not
        edge-manifold (some edge is shared by more than two faces).
    RuntimeError
        If ``vertices``, ``faces``, ``cot_entries``, ``unique_edges`` and ``edge_map`` are not all
        on one device.

    See Also
    --------
    [`crouzeix_raviart_massmatrix`][triwarp.energies.crouzeix_raviart_massmatrix]
    [`cotmatrix`][triwarp.laplacian.cotmatrix]
    [`edges_unique`][triwarp.edges.edges_unique]

    Notes
    -----
    Matches ``igl::crouzeix_raviart_cotmatrix`` up to the edge numbering, which follows
    [`edges_unique`][triwarp.edges.edges_unique] rather than ``igl::unique_edge_map``.
    """
    require_same_device(
        vertices=vertices,
        faces=faces,
        cot_entries=cot_entries,
        unique_edges=unique_edges,
        edge_map=edge_map,
    )
    unique_edges, edge_map = _edge_numbering(vertices, faces, unique_edges, edge_map)
    n_edges = int(unique_edges.shape[0])
    n_faces = int(faces.shape[0]) // 3
    device = faces.device
    if n_faces == 0:
        return tw.array.empty_square_bsr(n_edges, dtype, device)

    if not tw.validation.is_edge_manifold(faces, n_vertices=int(vertices.shape[0])):
        raise ValueError(
            "mesh must be edge-manifold (every edge shared by at most two faces); the "
            "Crouzeix-Raviart discretization is undefined otherwise, like "
            "igl::crouzeix_raviart_cotmatrix, which asserts it"
        )

    if cot_entries is None:
        cot_entries = cotmatrix_entries(vertices, faces, dtype=dtype)

    n_triplets = 12 * n_faces
    rows, cols, vals = tw.array.triplet_buffers(n_triplets, dtype, device)
    wp.launch(
        kernel_energies.CROUZEIX_RAVIART_COTMATRIX_TRIPLETS[cot_entries.dtype, dtype],
        dim=n_faces,
        inputs=[edge_map, cot_entries, rows, cols, vals],
        device=device,
    )
    return wps.bsr_from_triplets(n_edges, n_edges, rows, cols, vals, prune_numerical_zeros=False)

crouzeix_raviart_massmatrix

crouzeix_raviart_massmatrix(
    vertices: array[vec3],
    faces: array[int32],
    dtype: type = wp.float32,
    *,
    unique_edges: Array2dInt32 | None = None,
    edge_map: array[int32] | None = None,
) -> wps.BsrMatrix[wp.float32]

Edge-based Crouzeix-Raviart mass matrix.

Diagonal (n_edges, n_edges): each face donates a third of its area to each of its three edges, so an interior edge's entry is a third of its two incident faces' summed area. Rows follow edges_unique's edge numbering, the same numbering crouzeix_raviart_cotmatrix uses.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) mesh vertex positions.	required
faces	array[int32]	Length-3 * n_faces wp.int32 triangle index buffer.	required
dtype	type	Scalar block type of the assembled matrix: wp.float32 (default) or wp.float64.	float32
unique_edges	Array2dInt32 | None	Optional precomputed edge numbering from edges_unique — pass both or neither.	None
edge_map	Array2dInt32 | None	Optional precomputed edge numbering from edges_unique — pass both or neither.	None

Returns:

Type	Description
BsrMatrix	Square (n_edges, n_edges) diagonal mass matrix in 1x1-block BSR form on vertices.device.

Raises:

Type	Description
ValueError	If exactly one of unique_edges / edge_map is provided.
RuntimeError	If vertices, faces, unique_edges and edge_map are not all on one device.

See Also

crouzeix_raviart_cotmatrix mass_matrix edges_unique

Notes

Matches igl::crouzeix_raviart_massmatrix up to the edge numbering, as above.

Source code in triwarp/energies.py

def crouzeix_raviart_massmatrix(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    dtype: type = wp.float32,
    *,
    unique_edges: twt.Array2dInt32 | None = None,
    edge_map: wp.array[wp.int32] | None = None,
) -> wps.BsrMatrix[wp.float32]:
    """
    Edge-based Crouzeix-Raviart mass matrix.

    Diagonal ``(n_edges, n_edges)``: each face donates a third of its area to each of its three
    edges, so an interior edge's entry is a third of its two incident faces' summed area. Rows
    follow [`edges_unique`][triwarp.edges.edges_unique]'s edge numbering, the same numbering
    [`crouzeix_raviart_cotmatrix`][triwarp.energies.crouzeix_raviart_cotmatrix] uses.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.
    dtype
        Scalar block type of the assembled matrix: ``wp.float32`` (default) or ``wp.float64``.
    unique_edges, edge_map
        Optional precomputed edge numbering from
        [`edges_unique`][triwarp.edges.edges_unique] — pass both or neither.

    Returns
    -------
    warp.sparse.BsrMatrix
        Square ``(n_edges, n_edges)`` diagonal mass matrix in 1x1-block BSR form on
        ``vertices.device``.

    Raises
    ------
    ValueError
        If exactly one of ``unique_edges`` / ``edge_map`` is provided.
    RuntimeError
        If ``vertices``, ``faces``, ``unique_edges`` and ``edge_map`` are not all on one device.

    See Also
    --------
    [`crouzeix_raviart_cotmatrix`][triwarp.energies.crouzeix_raviart_cotmatrix]
    [`mass_matrix`][triwarp.laplacian.mass_matrix]
    [`edges_unique`][triwarp.edges.edges_unique]

    Notes
    -----
    Matches ``igl::crouzeix_raviart_massmatrix`` up to the edge numbering, as above.
    """
    require_same_device(
        vertices=vertices, faces=faces, unique_edges=unique_edges, edge_map=edge_map
    )
    unique_edges, edge_map = _edge_numbering(vertices, faces, unique_edges, edge_map)
    n_edges = int(unique_edges.shape[0])

    return wps.bsr_diag(diag=_cr_mass_diagonal(vertices, faces, edge_map, n_edges, dtype))

lscm_hessian

lscm_hessian(
    vertices: array[vec3], faces: array[int32]
) -> wps.BsrMatrix[wp.float64]

LSCM Hessian Q = -repdiag(L, 2) - 2 A.

Assembles the (2n, 2n) symmetric operator behind the least-squares conformal map, where L is the cotangent Laplacian cotmatrix (negative-diagonal convention), repdiag(L, 2) is the block-diagonal [[L, 0], [0, L]], and A is the boundary vector_area_matrix. Built natively in float64 in a single bsr_from_triplets (the within-quadrant repdiag triplets and the cross-quadrant -2 A triplets never collide), so it feeds the float64 conjugate-gradient solve directly. Matches the Q returned by igl.lscm exactly.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) mesh vertex positions.	required
faces	array[int32]	Length-3 * n_faces wp.int32 triangle index buffer.	required

Returns:

Type	Description
BsrMatrix	Square (2n, 2n) float64 matrix in 1x1-block BSR form on vertices.device.

Raises:

Type	Description
RuntimeError	If vertices and faces are not all on one device.

See Also

lscm vector_area_matrix cotmatrix

Notes

Matches igl::lscm_hessian.

Source code in triwarp/energies.py

def lscm_hessian(
    vertices: wp.array[wp.vec3], faces: wp.array[wp.int32]
) -> wps.BsrMatrix[wp.float64]:
    """
    LSCM Hessian ``Q = -repdiag(L, 2) - 2 A``.

    Assembles the ``(2n, 2n)`` symmetric operator behind the least-squares conformal map, where
    ``L`` is the cotangent Laplacian [`cotmatrix`][triwarp.laplacian.cotmatrix] (negative-diagonal
    convention), ``repdiag(L, 2)`` is the block-diagonal ``[[L, 0], [0, L]]``, and ``A`` is the
    boundary [`vector_area_matrix`][triwarp.energies.vector_area_matrix]. Built natively in
    float64 in a single ``bsr_from_triplets`` (the within-quadrant repdiag triplets and the
    cross-quadrant ``-2 A`` triplets never collide), so it feeds the float64 conjugate-gradient
    solve directly. Matches the ``Q`` returned by ``igl.lscm`` exactly.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.

    Returns
    -------
    warp.sparse.BsrMatrix
        Square ``(2n, 2n)`` float64 matrix in 1x1-block BSR form on ``vertices.device``.

    Raises
    ------
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`lscm`][triwarp.parametrization.lscm]
    [`vector_area_matrix`][triwarp.energies.vector_area_matrix]
    [`cotmatrix`][triwarp.laplacian.cotmatrix]

    Notes
    -----
    Matches ``igl::lscm_hessian``.
    """
    require_same_device(vertices=vertices, faces=faces)
    n = int(vertices.shape[0])
    device = vertices.device
    laplacian = cotmatrix(vertices, faces, dtype=wp.float64)
    # The real compressed-CSR entry count is offsets[-1], not laplacian.nnz: bsr_from_triplets
    # reports nnz as the (over-allocated) triplet capacity, so sizing by nnz would leave an
    # uninitialized gap in the wp.empty buffers that bsr_from_triplets reads back as garbage.
    n_entries = int(read_scalar(laplacian.offsets, n))
    boundary = tw.boundary.oriented_boundary_edges(vertices, faces)
    n_be = int(boundary.shape[0])

    # Combined triplet buffers: 2 per Laplacian entry (the two diagonal blocks) plus 4 per oriented
    # boundary edge (the vector-area cross-quadrant terms). Every slot is written, so wp.empty.
    total = 2 * n_entries + 4 * n_be
    rows, cols, vals = tw.array.triplet_buffers(total, wp.float64, device)
    wp.launch(
        kernel_energies.neg_repdiag2_triplets,
        dim=n,
        inputs=[
            laplacian.offsets,
            laplacian.columns,
            laplacian.values,
            wp.int32(n),
            rows,
            cols,
            vals,
        ],
        device=device,
    )
    if n_be > 0:
        # The ``-2 A`` term shares its triplet kernel with
        # [`vector_area_matrix`][triwarp.energies.vector_area_matrix] but writes into a slice
        # of the combined buffer: assembling ``A`` as its own matrix and adding it would need a
        # second build plus a ``bsr_axpy``, where this one pass over exact-size buffers does.
        _vector_area_triplets(
            boundary, n, -2.0, rows[2 * n_entries :], cols[2 * n_entries :], vals[2 * n_entries :]
        )
    return wps.bsr_from_triplets(2 * n, 2 * n, rows, cols, vals, prune_numerical_zeros=False)

vector_area_matrix

vector_area_matrix(
    vertices: array[vec3], faces: array[int32]
) -> wps.BsrMatrix[wp.float64]

Boundary vector-area matrix A: the signed area enclosed by the UV boundary curve.

Assembles the (2n, 2n) matrix that turns the [u; v] quadratic form into the signed area enclosed by the boundary UV curve: for each oriented boundary edge (i, j) (from the face winding, via oriented_boundary_edges) it adds the cross-quadrant entries (i+n, j, -1/4), (j, i+n, -1/4), (i, j+n, +1/4), (j+n, i, +1/4). On a closed mesh (no boundary) A is the zero matrix. Built natively in float64 in a single bsr_from_triplets.

Parameters:

Name	Type	Description	Default
vertices	array[vec3]	(n_vertices,) mesh vertex positions; only the count and device are used.	required
faces	array[int32]	Length-3 * n_faces wp.int32 triangle index buffer.	required

Returns:

Type	Description
BsrMatrix	Square (2n, 2n) float64 matrix in 1x1-block BSR form on vertices.device.

Raises:

Type	Description
RuntimeError	If vertices and faces are not all on one device.

See Also

lscm_hessian lscm oriented_boundary_edges

Notes

Matches igl::vector_area_matrix.

Source code in triwarp/energies.py

def vector_area_matrix(
    vertices: wp.array[wp.vec3], faces: wp.array[wp.int32]
) -> wps.BsrMatrix[wp.float64]:
    """
    Boundary vector-area matrix ``A``: the signed area enclosed by the UV boundary curve.

    Assembles the ``(2n, 2n)`` matrix that turns the ``[u; v]`` quadratic form into the signed area
    enclosed by the boundary UV curve: for each **oriented** boundary edge ``(i, j)`` (from the face
    winding, via [`oriented_boundary_edges`][triwarp.boundary.oriented_boundary_edges]) it adds the
    cross-quadrant entries ``(i+n, j, -1/4)``, ``(j, i+n, -1/4)``, ``(i, j+n, +1/4)``,
    ``(j+n, i, +1/4)``. On a closed mesh (no boundary) ``A`` is the zero matrix. Built natively in
    float64 in a single ``bsr_from_triplets``.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions; only the count and device are used.
    faces
        Length-``3 * n_faces`` ``wp.int32`` triangle index buffer.

    Returns
    -------
    warp.sparse.BsrMatrix
        Square ``(2n, 2n)`` float64 matrix in 1x1-block BSR form on ``vertices.device``.

    Raises
    ------
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`lscm_hessian`][triwarp.energies.lscm_hessian]
    [`lscm`][triwarp.parametrization.lscm]
    [`oriented_boundary_edges`][triwarp.boundary.oriented_boundary_edges]

    Notes
    -----
    Matches ``igl::vector_area_matrix``.
    """
    require_same_device(vertices=vertices, faces=faces)
    n = int(vertices.shape[0])
    device = vertices.device
    boundary = tw.boundary.oriented_boundary_edges(vertices, faces)
    n_be = int(boundary.shape[0])
    if n_be == 0:
        return tw.array.empty_square_bsr(2 * n, wp.float64, device)

    rows, cols, vals = tw.array.triplet_buffers(4 * n_be, wp.float64, device)
    _vector_area_triplets(boundary, n, 1.0, rows, cols, vals)
    return wps.bsr_from_triplets(2 * n, 2 * n, rows, cols, vals, prune_numerical_zeros=False)