Transformations and Projections#

Geometric Transformations#

Linear and affine transformations on mesh geometry. Each function returns a new Mesh with transformed point coordinates and appropriately invalidated caches. Cached quantities such as normals and areas are automatically recomputed on next access.

All transformations are also available as methods on Mesh.

import numpy as np
from physicsnemo.mesh.primitives.surfaces import sphere_icosahedral

mesh = sphere_icosahedral.load(subdivisions=3)

# Via Mesh methods
translated = mesh.translate([1.0, 0.0, 0.0])
rotated = mesh.rotate(axis=[0, 0, 1], angle=np.pi / 4)
scaled = mesh.scale(2.0)
scaled_aniso = mesh.scale([2.0, 1.0, 0.5])

# Arbitrary linear transform
import torch
matrix = torch.eye(3) * 2
transformed = mesh.transform(matrix)

Linear and affine transformations for simplicial meshes.

This module implements geometric point transformations with intelligent cache handling. Topology caches survive coordinate-only changes; geometry caches are invalidated unless a transformation explicitly preserves or updates them.

Cached fields handled: - areas: point and cell cache categories - normals: point and cell cache categories - centroids: cell cache category only

physicsnemo.mesh.transformations.geometric.rotate( mesh: Mesh, angle: float, axis: Float[Tensor, 'n_spatial_dims'] | Sequence[float] | Literal['x', 'y', 'z'] | None = None, center: Float[Tensor, 'n_spatial_dims'] | Sequence[float] | None = None, transform_point_data: bool | TensorDict = False, transform_cell_data: bool | TensorDict = False, transform_global_data: bool | TensorDict = False, ) → Mesh[source]#

Rotate the mesh about an axis by a specified angle.

Call it as rotate(mesh, ...) or as mesh.rotate(...). The bound method supplies mesh automatically.

Parameters:

mesh (Mesh) – Input mesh to rotate.
angle (float) – Rotation angle in radians (counterclockwise, right-hand rule).
axis (Float[torch.Tensor, " n_spatial_dims"] or Sequence[float] or {"x", "y", "z"} or None) – Rotation axis vector. None for 2D, shape \((3,)\) for 3D. String literals "x", "y", "z" are converted to unit vectors (1,0,0), (0,1,0), (0,0,1) respectively.
center (Float[torch.Tensor, " n_spatial_dims"] or Sequence[float] or None) – Center point for rotation. If None, rotates about the origin.
transform_point_data (bool or TensorDict) – Controls transformation of point_data fields. See transform() for full semantics.
transform_cell_data (bool or TensorDict) – Same semantics as transform_point_data, for cell_data.
transform_global_data (bool or TensorDict) – Same semantics as transform_point_data, for global_data.

Returns:

New Mesh with rotated geometry.

Return type:

Mesh

Notes

Cache Handling:

areas: Unchanged (rotation preserves volumes)

centroids: Rotated

normals: Rotated

physicsnemo.mesh.transformations.geometric.rotation_matrix( angle: float, axis: Float[Tensor, 'n_spatial_dims'] | Sequence[float] | Literal['x', 'y', 'z'] | None, n_spatial_dims: int, device: device, dtype: dtype, ) → Float[Tensor, 'n_spatial_dims n_spatial_dims'][source]#

Build a rotation matrix from angle and axis.

Parameters:

angle (float) – Rotation angle in radians (counterclockwise, right-hand rule).
axis (Float[torch.Tensor, " n_spatial_dims"] or Sequence[float] or {"x", "y", "z"} or None) – Rotation axis. None for 2D, tensor/sequence/string for 3D.
n_spatial_dims (int) – Number of spatial dimensions.
device (torch.device) – Target device for the output matrix.
dtype (torch.dtype) – Target dtype for the output matrix.

Returns:

Rotation matrix, shape \((S, S)\).

Return type:

Float[torch.Tensor, “n_spatial_dims n_spatial_dims”]

physicsnemo.mesh.transformations.geometric.scale( mesh: Mesh, factor: float | Float[Tensor, 'n_spatial_dims'] | Sequence[float], center: Float[Tensor, 'n_spatial_dims'] | Sequence[float] | None = None, transform_point_data: bool | TensorDict = False, transform_cell_data: bool | TensorDict = False, transform_global_data: bool | TensorDict = False, assume_invertible: bool | None = None, ) → Mesh[source]#

Scale the mesh by specified factor(s).

Call it as scale(mesh, ...) or as mesh.scale(...). The bound method supplies mesh automatically.

Parameters:

mesh (Mesh) – Input mesh to scale.
factor (float or Float[torch.Tensor, " n_spatial_dims"] or Sequence[float]) – Scale factor(s). Scalar for uniform, vector for non-uniform.
center (Float[torch.Tensor, " n_spatial_dims"] or Sequence[float] or None) – Center point for scaling. If None, scales about the origin.
transform_point_data (bool or TensorDict) – Controls transformation of point_data fields. See transform() for full semantics.
transform_cell_data (bool or TensorDict) – Same semantics as transform_point_data, for cell_data.
transform_global_data (bool or TensorDict) – Same semantics as transform_point_data, for global_data.
assume_invertible (bool or None) –
Controls cache propagation:
- True: Assume all factors are non-zero, propagate caches (compile-safe)
- False: Assume some factor is zero, skip cache propagation (compile-safe)
- None: Check determinant at runtime (may cause graph breaks under torch.compile)

Returns:

New Mesh with scaled geometry.

Return type:

Mesh

Notes

Cache Handling:

areas: Scaled correctly. For non-isotropic transforms of codimension-1
embedded manifolds, per-element scaling is computed using normals.

centroids: Scaled

normals: Transformed by inverse-transpose (direction adjusted, magnitude normalized)

physicsnemo.mesh.transformations.geometric.scale_matrix( factor: float | Float[Tensor, 'n_spatial_dims'] | Sequence[float], n_spatial_dims: int, device: device, dtype: dtype, ) → Float[Tensor, 'n_spatial_dims n_spatial_dims'][source]#

Build a diagonal scale matrix from a factor specification.

Parameters:

factor (float or Float[torch.Tensor, " n_spatial_dims"] or Sequence[float]) – Scale factor(s). Scalar for uniform, vector for non-uniform.
n_spatial_dims (int) – Number of spatial dimensions.
device (torch.device) – Target device for the output matrix.
dtype (torch.dtype) – Target dtype for the output matrix.

Returns:

Diagonal scale matrix, shape \((S, S)\).

Return type:

Float[torch.Tensor, “n_spatial_dims n_spatial_dims”]

Raises:

ValueError – If factor is a vector whose length does not match n_spatial_dims.

physicsnemo.mesh.transformations.geometric.transform( mesh: Mesh, matrix: Float[Tensor, 'new_n_spatial_dims n_spatial_dims'], transform_point_data: bool | TensorDict = False, transform_cell_data: bool | TensorDict = False, transform_global_data: bool | TensorDict = False, assume_invertible: bool | None = None, ) → Mesh[source]#

Apply a linear transformation to the mesh.

Call it as transform(mesh, ...) or as mesh.transform(...). The bound method supplies mesh automatically.

Parameters:

mesh (Mesh) – Input mesh to transform.
matrix (Float[torch.Tensor, "new_n_spatial_dims n_spatial_dims"]) – Transformation matrix, shape \((S', S)\).
transform_point_data (bool or TensorDict) – Controls transformation of point_data fields. True transforms all compatible fields; False transforms none; a TensorDict (or dict) with scalar bool leaves selectively transforms only the named fields.
transform_cell_data (bool or TensorDict) – Same semantics as transform_point_data, for cell_data.
transform_global_data (bool or TensorDict) – Same semantics as transform_point_data, for global_data.
assume_invertible (bool or None) –
Controls cache propagation for square matrices:
- True: assume matrix is invertible and propagate caches (compile-safe). This is a promise, not a check. If matrix is in fact singular, the inverse-transpose step silently yields non-finite values instead of raising, so the propagated normals and areas caches – and anything derived from them, such as the sum of cell_areas – come back as NaN. Use False or None unless you know the matrix is non-singular.
- False: assume matrix is singular and skip cache propagation (compile-safe). Caches are dropped and recomputed lazily on demand, which is always correct, just slower.
- None (default): test abs(det(matrix)) > 1e-10 at runtime and take one of the branches above. Safe for singular input, but the test reads a device scalar back to the host, which synchronizes on CUDA and may cause graph breaks under torch.compile.

Returns:

New Mesh with transformed geometry and appropriately updated caches.

Return type:

Mesh

Notes

Cache Handling:

areas: For square invertible matrices:

Full-dimensional meshes: scaled by |det|

Codimension-1 manifolds: per-element scaling using |det| * ||M^{-T} n||

Higher codimension: invalidated

centroids: Always transformed

normals: For square invertible matrices, transformed by inverse-transpose

physicsnemo.mesh.transformations.geometric.translate( mesh: Mesh, offset: Float[Tensor, 'n_spatial_dims'] | Sequence[float], ) → Mesh[source]#

Apply a translation to the mesh.

Translation only affects point positions and centroids. Vector/tensor fields are unchanged by translation (they represent directions, not positions).

Call it as translate(mesh, ...) or as mesh.translate(...). The bound method supplies mesh automatically.

Parameters:

mesh (Mesh) – Input mesh to translate.
offset (Float[torch.Tensor, " n_spatial_dims"] or Sequence[float]) – Translation vector, shape \((S,)\).

Returns:

New Mesh with translated geometry.

Return type:

Mesh

Notes

Cache Handling:

areas: Unchanged

centroids: Translated

normals: Unchanged

Deformations#

The deform namespace provides six deformation families:

Dense displacement through displace(), backed by displace_points().
Sobolev-filtered dense displacement through sobolev_deform(), backed by sobolev_deform_points().
Compact sparse-control morphing through morph(), backed by morph_points().
Global radial-basis deformation through radial_basis_function_deform(), backed by radial_basis_function_deform_points().
Lattice free-form deformation through free_form_deform(), backed by free_form_deform_points().
Nearest-surface conformance through shrinkwrap(), backed by shrinkwrap_points().

Each operation is also available as a method on Mesh.

Dense displacement accepts a tensor or a point-data key (including a nested tuple key). The operation returns a new mesh without changing mesh.points. Assigning a point-data key, as in the second example below, is a separate, explicit mutation of the source mesh’s attached data.

displacement = torch.zeros_like(mesh.points)
displacement[:, 2] = 0.05
displaced = mesh.displace(0.5 * displacement)

# Point-data fields can drive the same operation.
mesh.point_data["design_displacement"] = displacement
displaced_from_data = mesh.displace("design_displacement")

Sobolev-Filtered Deformation#

sobolev_deform() turns a raw displacement \(d\) into a spatially smooth displacement \(u\) by solving

\[(M + \ell^2 K)u = M d.\]

Here \(M=\bar m I\) is a uniform vertex mass matrix scaled by the mean positive lumped P1 mass, and \(K\) is the P1 stiffness matrix. The uniform mass makes the filter self-adjoint in standard Euclidean vertex coordinates, so the same operator smooths the displacement and its adjoint. length_scale is \(\ell\) in the same physical units as mesh.points. Larger values suppress variation over longer distances. A zero length scale applies the raw displacement exactly at unfixed vertices. This discrete Helmholtz filter is motivated by the implicit formulation for node-based shape optimization studied by Najian Asl and Bletzinger.

reference_points = mesh.points.detach()
candidate_vertices = reference_points.clone().requires_grad_()
raw_displacement = candidate_vertices - reference_points
smooth = mesh.sobolev_deform(
    raw_displacement,
    length_scale=0.2,
)

objective = smooth.points.square().mean()
objective.backward()
smooth_vertex_adjoint = candidate_vertices.grad

The forward deformation is the uniform-mass discrete Helmholtz solve itself. Its first-order implicit backward solves the matching adjoint system. Optimize candidate vertex coordinates by subtracting fixed reference points and passing that offset as raw_displacement. The adjoint with respect to those candidate vertices receives the same Sobolev filter. Gradients with respect to the reference coordinates also include the geometric dependence of the P1 stiffness and mass scale. Higher-order derivatives are not supported for a positive length scale. Each ambient displacement component is filtered independently.

fixed_points accepts a boolean point mask or a point-data key. True entries receive zero displacement, which imposes a homogeneous Dirichlet condition. Other mesh boundaries use the natural homogeneous Neumann condition. With no fixed points, constant displacements pass through to solver precision.

The Torch and Warp implementations use matrix-free preconditioned conjugate gradients. The Warp implementation runs on CUDA and supplies an explicit implicit-adjoint backward with an analytic geometry vector-Jacobian product. By default, CUDA segments, triangles, and tetrahedra select Warp when it is available. CPU meshes and higher-dimensional simplices select Torch. max_iterations and tolerance control the iterative solve. The operation raises an error if either the forward or adjoint solve does not reach the requested tolerance. At positive length scales, cells must be finite, nondegenerate simplices. Isolated points receive their raw displacement. CUDA Graph capture is not supported because P1 operator assembly and solver diagnostics are not capture-safe. Both backends support torch.compile. Warp CUDA results and point gradients may vary at roundoff between runs.

Before and After

The panels show the initial objective adjoint with respect to candidate vertex coordinates. The raw field contains checkerboard-scale oscillation. Applying the Sobolev deformation filters that adjoint over the mesh while retaining the fixed boundary. Both panels use the same color and arrow scale.

Raw noisy vertex adjoint and smoother Sobolev-filtered adjoint

Three-Dimensional Surface Example

The 3D figure applies the same workflow to a triangulated sheet with points shaped (N, 3). Its objective pulls the center upward while the boundary remains fixed. The reproducible figure source is docs/img/mesh/sobolev_adjoint_field_3d.py.

The arrows show the normalized negative adjoint, which is the gradient-descent update direction. Each panel applies that field to the sheet geometry for visibility. The raw panel is corrugated by vertex-scale oscillations. The Sobolev panel gives the same broad upward pull with a smooth surface. Orange points mark the fixed boundary, and both panels use one shared scale.

Raw and Sobolev-filtered upward updates on a fixed-boundary sheet

Sparse controls are useful when only a small set of design handles is known. A control point is a location in world coordinates, and its control displacement is a vector rather than a destination coordinate. Control points do not need to be mesh vertices, although selecting vertices makes their prescribed movement directly visible in the result.

Single-Control Morphing#

Indexing one vertex produces a coordinate vector with shape (3,). The morph API instead expects (n_controls, n_spatial_dims), so unsqueeze(0) adds the control dimension and gives shape (1, 3). It is not a batch dimension.

top_index = mesh.points[:, 2].argmax()
control_points = mesh.points[top_index].unsqueeze(0)  # (1, 3)
control_displacements = mesh.points.new_tensor(
    [[0.0, 0.0, 0.5]], requires_grad=True
)
single_morph = mesh.morph(
    control_points,
    control_displacements,
    radius=1.0,
)

# Autograd continues through the returned point coordinates.
objective = single_morph.points.square().mean()
objective.backward()

control_displacements is differentiable. An optimizer can learn it from any differentiable loss computed from single_morph.points. A model can also predict the displacements.

Without point weights, a mesh vertex exactly at a unique control moves by its prescribed displacement. point_weights can scale or mask the final movement. Duplicate controls at the same coordinate contribute their mean displacement.

Multiple-Control Morphing#

Advanced indexing retains the control dimension when several vertices are selected. Each row of control_points pairs with the same row of control_displacements and, when supplied, one entry of radius.

bottom_index = mesh.points[:, 2].argmin()
control_indices = torch.stack((top_index, bottom_index))
control_points = mesh.points[control_indices]  # (2, 3)
control_displacements = mesh.points.new_tensor(
    [[0.0, 0.0, 0.5], [0.0, 0.0, -0.5]]
)
radii = mesh.points.new_tensor([1.0, 1.0])

multiple_morph = mesh.morph(
    control_points,
    control_displacements,
    radius=radii,
)

The radius is a Euclidean support distance in the same coordinate units as the mesh. A control’s influence vanishes smoothly at its support boundary. Where supports overlap, all active controls are evaluated together using a stationary zero-displacement background. The result is not a simple sum or average. Points outside every support remain unchanged. Put simultaneous controls in one call, because applying several morphs sequentially evaluates later fields on already modified coordinates and is therefore order-dependent.

The kernel keyword names the compact radial kernel used by the field. "wendland_c2" is currently the supported value and the default.

Every tensor-valued radius must remain finite and strictly positive. Its values are not validated at runtime. When a model learns the radius, use a positive parameterization such as torch.nn.functional.softplus(raw_radius) + radius_epsilon rather than optimizing an unconstrained radius directly. Floating point_weights are used as supplied and may be signed or greater than one.

Visualization

The panels compare the original sphere with the single-control and multiple-control examples above. Green markers identify the displaced handle locations, while arrows and labels show the prescribed displacement directions and magnitudes.

Original sphere and single-control and multiple-control sphere morphing

Global Radial-Basis Deformation#

Radial-basis deformation fits one global displacement field through sparse handles. With zero smoothing and a nonsingular control layout, the fitted field interpolates each prescribed control displacement up to solver precision. Optional point weights are applied after interpolation. Unlike compact Shepard morphing, every control generally influences every point. Fixed controls are therefore useful as anchors.

The standard affine polynomial tail reproduces affine displacement fields. The controls must affinely span the coordinate space, and the augmented system must be nonsingular. This formulation follows the thin-plate-spline interpolant described by Bookstein [1].

[1] F. L. Bookstein, “Principal warps: thin-plate splines and the decomposition of deformations,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 11, no. 6, pp. 567-585, 1989. https://doi.org/10.1109/34.24792

Two-Dimensional Example

A mesh with n_spatial_dims=2 requires at least three non-collinear controls when the affine tail is enabled. This example fixes the four corners of a triangulated square and moves a fifth handle at the midpoint of its upper edge.

import torch
from physicsnemo.mesh.primitives.planar import unit_square

mesh_2d = unit_square.load(subdivisions=4)
controls_2d = mesh_2d.points.new_tensor(
    [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [0.5, 1.0]]
)
displacements_2d = torch.zeros_like(controls_2d)
displacements_2d[-1] = controls_2d.new_tensor([0.15, 0.35])

deformed_2d = mesh_2d.radial_basis_function_deform(
    controls_2d,
    displacements_2d,
    kernel="thin_plate_spline",
    polynomial=True,
    smoothing=0.0,
)

handle_index = torch.linalg.vector_norm(
    mesh_2d.points - controls_2d[-1], dim=1
).argmin()
torch.testing.assert_close(
    deformed_2d.points[handle_index],
    controls_2d[-1] + displacements_2d[-1],
    atol=2.0e-5,
    rtol=2.0e-5,
)

The corner controls remain fixed while the upper edge and interior deform smoothly. The hollow green marker in the output panel shows the handle’s source position.

Two-dimensional global RBF deformation of a triangulated square

Three-Dimensional Example

In three dimensions, the affine tail requires at least four controls spanning the coordinate space. Practical layouts usually use more. The sphere example uses its six axis-extreme vertices as controls.

control_indices = torch.stack(
    (
        mesh.points[:, 2].argmax(),
        mesh.points[:, 2].argmin(),
        mesh.points[:, 0].argmax(),
        mesh.points[:, 0].argmin(),
        mesh.points[:, 1].argmax(),
        mesh.points[:, 1].argmin(),
    )
)
rbf_controls = mesh.points[control_indices]
rbf_displacements = torch.zeros_like(rbf_controls)
rbf_displacements[0, 2] = 0.55

exact_rbf = mesh.radial_basis_function_deform(
    rbf_controls,
    rbf_displacements,
    kernel="thin_plate_spline",
    polynomial=True,
    smoothing=0.0,
)

torch.testing.assert_close(
    exact_rbf.points[control_indices],
    rbf_controls + rbf_displacements,
    atol=2.0e-5,
    rtol=2.0e-5,
)

smoothing=0.0 interpolates the controls up to solver precision. A positive smoothing value adds diagonal regularization. This deliberately relaxes interpolation accuracy.

Both evaluation backends use PyTorch for the dense coefficient solve. The Warp backend fuses evaluation over the mesh points without materializing its full point/control kernel matrix.

Orange controls are fixed anchors. Green controls mark moved handles, and the arrows show their displacement directions. The labels give the prescribed magnitudes. With zero smoothing, the deformed surfaces interpolate all six control displacements up to solver precision.

Original sphere with one-handle and two-handle RBF deformations

Lattice Free-Form Deformation#

free_form_deform() defines a regular array of control displacements over an axis-aligned evaluation box and deforms every point inside the box by tensor-product basis interpolation. Compared with sparse morphing, the design parameters form a structured grid of fixed size, which suits parametric shape optimization. A lattice of zeros is exactly the identity, and the same lattice deforms any geometry embedded in the box.

When origin and extent are omitted, the box spans the mesh bounds. Each coordinate axis must have positive range. Planar or linear geometry embedded in a higher-dimensional space therefore needs an explicit positive extent. Validating an automatically derived extent synchronizes with the device and is not CUDA Graph capture-safe. For capture, pass both origin and extent as device tensors.

basis="bernstein" provides classic global-support free-form deformation for coarse lattices. basis="bspline" provides local four-node-per-axis support and scales to fine lattices for local sculpting. Its first and last coefficient planes lie one knot spacing outside the evaluation box. basis="linear", "cubic_hermite", and "quintic_hermite" instead use the two neighboring nodes per axis and reproduce every control displacement at its lattice node. For a local cell coordinate \(t\), their upper-node weights are \(t\), \(3t^2-2t^3\) (cubic Hermite), and \(6t^5-15t^4+10t^3\) (quintic Hermite), respectively. The lower-node weight is one minus the upper-node weight. The resulting fields are C0, C1, and C2 across cell boundaries, respectively. Perlin introduced the quintic blend in Improving Noise to eliminate the cubic blend’s second-derivative discontinuities.

# A 4x4x4 Bernstein lattice spans the mesh bounds.
# Zero displacements start at the identity.
control_displacements = torch.zeros(4, 4, 4, 3, requires_grad=True)
deformed = mesh.free_form_deform(control_displacements)

# Autograd continues through the returned point coordinates.
objective = deformed.points.square().mean()
objective.backward()

Points outside the lattice box are unchanged. The deformation is generally not continuous across the box boundary. A sufficient condition for a fixed exterior is to zero the outermost coefficient plane on every Bernstein or node-interpolating face. For cubic B-splines, zero the first and last three coefficient planes on every axis. origin and extent are non-differentiable lattice parameters. Optimize control_displacements instead.

Visualization

The panels compare the original sphere with two lattice deformations. A coarse Bernstein lattice tapers the whole sphere because every node acts globally, while three independently displaced regions in a finer cubic B-spline lattice produce two local bulges and a local indentation.

Sphere with control lattice, Bernstein taper, and local B-spline sculpt

Nearest-Surface Shrinkwrap#

shrinkwrap() projects source vertices to their closest points on a triangle target. An optional signed offset keeps the result a specified distance from the target along its oriented face normals. point_weights can fix selected vertices or apply a partial projection. max_distance leaves vertices unchanged when no target surface is close enough.

The following call is the central operation in the curved-panel example:

conformed = source.shrinkwrap(
    panel_target,
    point_weights=movable,
    max_distance=0.34,
)

The example constructs a swept triangle panel and adds a smooth springback error to form the source sheet. The boolean movable mask keeps its green root attachment strip fixed. The movable vertices lie on the target surface. Color shows how far each source vertex moved.

Curved triangle panel, lifted source sheet, and shrinkwrapped result

Nearest-face selection and closest-feature changes are discrete. With those choices fixed, gradients propagate through source points, selected target vertices, floating point weights, and a tensor-valued offset. At a shared edge or vertex, adjacent faces can provide different normals. Use consistently oriented target faces and avoid placing offset-sensitive samples exactly on shared features.

Torch provides the reference search. Warp accelerates nearest-face search on CPU and CUDA. Both backends replay the selected point-to-triangle projection with PyTorch in the input dtype. Shrinkwrap is available on Mesh, not DomainMesh, because one source mesh is projected onto one target surface.

Float64 targets use the Torch search because Warp searches in float32. Safe float32 coordinates are searched unchanged. Warp falls back to Torch for unsafe coordinate magnitudes or face geometry.

Shrinkwrap performs data-dependent validation and nearest-face search setup. CUDA executions with either backend are not supported inside CUDA Graph capture.

Shape-Optimization Constraint Example

In shape optimization, an otherwise acceptable design iterate can locally cross an admissible clearance surface. This example constructs a closed triangulated low-profile enclosure with an inset lid and chamfered walls. A boolean mask selects only the vertices above a horizontal clearance plane.

feasible = candidate.shrinkwrap(
    admissible_envelope,
    point_weights=design_region,
)

Before repair, a smooth orange dome protrudes through the blue clearance plane while the broad gold lid deformation remains admissible. Shrinkwrap projects only the violating vertices back to the limit. After repair, the cap is green, the gold optimized surface is retained, and an orange wire outline marks the former excess.

All connectivity is triangular. The example checks the repaired-cap residual, exact preservation of admissible vertices, closed two-manifold source and result connectivity, consistent winding, positive enclosed volume, adjoints, and Torch versus Warp agreement. Shrinkwrap restores this geometric constraint. It does not run the optimizer or guarantee valid cells for arbitrary inputs.

Selective repair of an optimized enclosure crossing a clearance plane

Domain Meshes#

morph() evaluates one world-coordinate control field on the interior and every named boundary. With point_weights=None, coincident component points receive identical motion. Domain point weights must instead be a point-data key (or nested tuple key) present in every component. Raw weight tensors are rejected because component point counts can differ. Every resolved field must use one common dtype across the domain: bool for a hard mask, or the same floating dtype as the mesh points. Coincident points remain coincident under a point-weight key only when their resolved values also match.

radial_basis_function_deform() follows the same component and point-weight rules while fitting one global RBF field. The coefficient system is solved once, and the combined interior and boundary points are evaluated together before the component meshes are rebuilt.

import torch
from physicsnemo.mesh import DomainMesh, Mesh

interior = Mesh(
    points=torch.tensor([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]]),
    cells=torch.tensor([[0, 1, 2]]),
    point_data={"design_weight": torch.tensor([1.0, 0.8, 0.5])},
)
wall = Mesh(
    points=interior.points[:2],
    cells=torch.tensor([[0, 1]]),
    point_data={"design_weight": torch.tensor([1.0, 0.8])},
)
domain = DomainMesh(interior=interior, boundaries={"wall": wall})

domain_controls = interior.points[[0]]  # (1, 2)
domain_displacements = interior.points.new_tensor([[0.0, 0.25]])
morphed_domain = domain.morph(
    domain_controls,
    domain_displacements,
    radius=1.25,
    point_weights="design_weight",
    implementation="torch",
)

# Equal point weights keep the shared wall vertices coincident.
assert torch.allclose(
    morphed_domain.interior.points[:2],
    morphed_domain.boundaries["wall"].points,
)

free_form_deform() follows the same pattern for lattice free-form deformation:

The operation evaluates one lattice field over the combined interior and boundary points.
The default box spans the combined component bounds.

The combined bounds must have positive range on every coordinate axis unless an explicit extent is supplied.

Every deformation preserves connectivity and attached point, cell, global, and domain data. These operations treat attached vector and tensor fields as Lagrangian data and do not push them forward. They discard geometry-dependent caches and recompute them lazily. They retain topology caches.

Warning

Deformations do not detect or repair inverted, degenerate, or self-intersecting output cells. Call validate() or validate() explicitly when a deformation could compromise validity.

For optimization-time geometric penalties on a fixed topology, see Differentiable Deformation Energies.

Displace every mesh point by a dense vector field.

Computes points + displacement without changing connectivity, optionally multiplying the displacement by point_weights[..., None]. displacement and point_weights may be raw tensors or keys (including nested tuple keys) in point_data.

Call it as displace(mesh, ...) or as mesh.displace(...). The bound method supplies mesh automatically.

Parameters:

mesh (Mesh) – Mesh whose points are displaced. The source mesh is not modified.
displacement (str, tuple[str, ...], or torch.Tensor) – Dense displacement vectors with shape (mesh.n_points, mesh.n_spatial_dims), or a point-data key resolving to such a tensor. The tensor and mesh.points must have the same float32 or float64 dtype and device.
point_weights (str, tuple[str, ...], torch.Tensor, or None, optional) – Optional bool or floating-point weights with shape (mesh.n_points,), or a point-data key resolving to those point weights. Floating-point weights may be signed or greater than one. Default is None.
implementation ({"torch"} or None, optional) – Backend override. None selects Torch for dense displacement.

Returns:

New mesh with displaced points and unchanged connectivity and attached fields.

Return type:

Mesh

Notes

Attached fields are treated as Lagrangian data and are not pushed forward. Geometry-dependent caches are invalidated and topology caches are retained. The operation does not detect or repair inverted, degenerate, or self-intersecting cells; call validate() explicitly when needed.

physicsnemo.mesh.transformations.deform.free_form_deform( mesh: Mesh, control_displacements: Float[Tensor, '*lattice_resolution n_spatial_dims'], *, origin: Float[Tensor, 'n_spatial_dims'] | Sequence[float] | None = None, extent: Float[Tensor, 'n_spatial_dims'] | Sequence[float] | None = None, basis: Literal['bernstein', 'bspline', 'linear', 'cubic_hermite', 'quintic_hermite'] = 'bernstein', point_weights: str | tuple[str, ...] | Bool[Tensor, 'n_points'] | Float[Tensor, 'n_points'] | None = None, implementation: Literal['torch', 'warp'] | None = None, ) → Mesh[source]#

Deform a mesh with a control-point lattice by free-form deformation.

An n_1 x ... x n_D array of control displacements defines a field over the axis-aligned box [origin, origin + extent]. Mesh points inside the box move with the tensor-product basis interpolation of those values. Points outside the box are unchanged. A lattice of zero displacements is exactly the identity, and a constant lattice translates every point inside the box.

Call it as free_form_deform(mesh, ...) or as mesh.free_form_deform(...). The bound method supplies mesh automatically.

Parameters:

mesh (Mesh) – Mesh whose points are deformed. The source mesh is not modified.
control_displacements (torch.Tensor) – Displacement vectors, not destination coordinates, for every lattice node, with shape (n_1, ..., n_D, mesh.n_spatial_dims) and the same float32 or float64 dtype and device as mesh.points. Each axis needs at least two nodes for "bernstein" and the node-interpolating bases, and four for "bspline".
origin (torch.Tensor, sequence of float, or None, optional) – Minimum corner of the lattice box with shape (mesh.n_spatial_dims,). None uses the minimum corner of the mesh bounds. For repeated GPU calls with an explicit box, create origin and extent once as device tensors. Reuse them to avoid recreating and transferring sequence values. Default is None.
extent (torch.Tensor, sequence of float, or None, optional) – Edge lengths of the lattice box with the same accepted shapes as origin. Every value must be finite and strictly positive. The operation does not validate tensor values at runtime. None sizes the box from origin to the maximum corner of the mesh bounds. Validating a derived extent synchronizes with the device and is not CUDA Graph capture-safe. For capture, pass both origin and extent as device tensors. Every point-coordinate axis must have positive range when the extent is derived. Supply an explicit extent for lower-dimensional geometry embedded in a higher-dimensional space. Default is None.
basis ({"bernstein", "bspline", "linear", "cubic_hermite", "quintic_hermite"}, optional) –
Per-axis basis family:
- "bernstein" provides classic global-support FFD. Every lattice node influences every point inside the box.
- "bspline" uses a uniform cubic B-spline with local four-node-per-axis support and C2 continuity between knot spans. Coefficient index i corresponds to local coordinate (i - 1) / (n - 3). The first and last coefficient planes lie outside the evaluation box.
- "linear" uses upper-node weight \(s(t)=t\) within each lattice cell. It interpolates every node. It is continuous (C0) across cell boundaries, but its slope can jump.
- "cubic_hermite" uses the cubic Hermite blend \(s(t)=3t^2-2t^3\). Its first derivative vanishes at both cell endpoints. This gives C1 continuity across cell boundaries.
- "quintic_hermite" uses the quintic Hermite blend \(s(t)=6t^5-15t^4+10t^3\). Its first and second derivatives vanish at both endpoints. This gives C2 continuity across cell boundaries. Perlin introduced this improved interpolant in “Improving Noise” [1].
The node-interpolating bases use only the two neighboring nodes per axis. Here, t is the local cell coordinate in [0, 1]. The upper-node weight is \(s(t)\), and the lower-node weight is \(1-s(t)\). Default is "bernstein".
point_weights (str, tuple[str, ...], torch.Tensor, or None, optional) – Optional bool or floating mesh-point weights with shape (mesh.n_points,), or a point_data key resolving to those point weights. All weights must match the point device. Floating weights must also match the point dtype. Default is None.
implementation ({"torch", "warp"} or None, optional) – Backend override. None selects Torch on CPU. On CUDA, it selects Warp when available and otherwise Torch.

Returns:

New mesh with deformed points and unchanged connectivity and fields.

Return type:

Mesh

Raises:

TypeError – If tensors, lattice values, or point weights have unsupported types or dtypes.
ValueError – If shapes, devices, lattice parameters, point weights, or basis are invalid.
KeyError – If a point-data key or implementation name is not found.
ImportError – If an explicitly requested backend is unavailable.

Notes

The operation treats attached fields as Lagrangian data and does not push them forward. It invalidates geometry-dependent caches and retains topology caches. The deformation is generally not continuous across the lattice box boundary. To keep the exterior fixed, zero the outermost coefficient plane on every Bernstein or node-interpolating face. For cubic B-splines, zero the first and last three coefficient planes on every axis. origin and extent are non-differentiable lattice parameters. Optimize control_displacements instead. The operation does not detect or repair inverted, degenerate, or self-intersecting cells. Call validate() explicitly when needed.

References

[1] Perlin, K. (2002). “Improving Noise.” ACM Transactions on Graphics, 21(3), 681-682. https://doi.org/10.1145/566654.566636

physicsnemo.mesh.transformations.deform.morph( mesh: Mesh, control_points: Tensor, control_displacements: Tensor, *, radius: float | Tensor, point_weights: str | tuple[str, ...] | Tensor | None = None, kernel: Literal['wendland_c2'] = 'wendland_c2', implementation: Literal['torch', 'warp'] | None = None, ) → Mesh[source]#

Morph a mesh from sparse, compactly supported control displacements.

Control influence uses a Wendland-C2 compact Shepard field with a stationary zero-displacement background. Each control’s influence vanishes smoothly at its support boundary. Where supports overlap, all active controls and the background are blended together; the result is not a simple sum or average. The field is zero outside the union of all supports.

With no point weights, the field at a unique control coordinate is exactly that control’s displacement. Duplicate controls at one coordinate contribute their mean displacement. A control may be anywhere in world coordinates and need not coincide with a mesh point.

Call it as morph(mesh, ...) or as mesh.morph(...). The bound method supplies mesh automatically.

Parameters:

mesh (Mesh) – Mesh whose points are morphed. The source mesh is not modified.
control_points (torch.Tensor) – World-coordinate controls with shape (n_controls, mesh.n_spatial_dims) and the same dtype and device as mesh.points.
control_displacements (torch.Tensor) – Displacement vectors, not destination coordinates, with exactly the same shape, dtype, and device as control_points.
radius (float or torch.Tensor) – Support distance in mesh coordinate units. Supply one scalar for every control or a tensor with shape (n_controls,) that matches the control dtype and device. Every tensor value must remain positive and finite; values are not validated at runtime.
point_weights (str, tuple[str, ...], torch.Tensor, or None, optional) – Optional bool or floating mesh-point weights with shape (mesh.n_points,), or a point_data key resolving to those point weights. These are query-point weights, not per-control values. Default is None.
kernel ({"wendland_c2"}, optional) – Compact radial kernel used to blend control displacements. Default is "wendland_c2".
implementation ({"torch", "warp"} or None, optional) – Backend override. None selects Torch on CPU and Warp on CUDA when Warp is available, otherwise Torch.

Returns:

New mesh with morphed points and unchanged connectivity and attached fields.

Return type:

Mesh

Notes

Attached fields are treated as Lagrangian data and are not pushed forward. Geometry-dependent caches are invalidated and topology caches are retained. Parameterize learned radii to remain positive, for example as torch.nn.functional.softplus(raw_radius) + eps. The operation does not detect or repair inverted, degenerate, or self-intersecting cells; call validate() explicitly when needed.

physicsnemo.mesh.transformations.deform.radial_basis_function_deform( mesh: Mesh, control_points: Float[Tensor, 'n_controls n_spatial_dims'], control_displacements: Float[Tensor, 'n_controls n_spatial_dims'], *, kernel: Literal['thin_plate_spline'] = 'thin_plate_spline', polynomial: bool = True, smoothing: float = 0.0, point_weights: str | tuple[str, ...] | Bool[Tensor, 'n_points'] | Float[Tensor, 'n_points'] | None = None, implementation: Literal['torch', 'warp'] | None = None, ) → Mesh[source]#

Deform a mesh with a global thin-plate-spline RBF field.

A thin-plate-spline radial field is fitted to the prescribed sparse control displacements and evaluated at every mesh point. With the default affine polynomial tail, zero smoothing, and a nonsingular control layout, the unweighted field interpolates every control displacement up to solver precision.

Call it as radial_basis_function_deform(mesh, ...) or as mesh.radial_basis_function_deform(...). The bound method supplies mesh automatically.

Parameters:

mesh (Mesh) – Mesh whose points are deformed. The source mesh is not modified.
control_points (torch.Tensor) – World-coordinate controls with shape (n_controls, mesh.n_spatial_dims) and the same dtype and device as mesh.points.
control_displacements (torch.Tensor) – Displacement vectors, not destination coordinates, with exactly the same shape, dtype, and device as control_points.
kernel ({"thin_plate_spline"}, optional) – Radial kernel used by the interpolant. Default is "thin_plate_spline".
polynomial (bool, optional) – Add the standard affine polynomial tail and side constraints. This reproduces affine displacement fields. The controls must affinely span the coordinate space, and the augmented system must be nonsingular. Default is True.
smoothing (float, optional) – Nonnegative diagonal regularization added to the radial system. Zero gives exact interpolation for a nonsingular control layout up to solver precision. Positive values relax interpolation accuracy. Default is 0.0.
point_weights (str, tuple[str, ...], torch.Tensor, or None, optional) – Optional bool or floating mesh-point weights with shape (mesh.n_points,), or a point_data key resolving to those weights. Values scale or mask the fitted field after interpolation. Bool weights must be on the same device as the mesh points. Floating weights must have the same dtype and device as the mesh points.
implementation ({"torch", "warp"} or None, optional) – Evaluation-backend override. Both backends use PyTorch for the dense coefficient solve. None selects Torch on CPU and Warp on CUDA when Warp is available, otherwise Torch.

Returns:

New mesh with deformed points and unchanged connectivity and attached fields.

Return type:

Mesh

Raises:

TypeError – If control tensors or Python arguments have unsupported types, or if tensor dtypes are unsupported or mismatched.
ValueError – If tensor shapes, devices, control layout, point weights, or RBF options are invalid.
KeyError – If a point-data key is missing or implementation does not name a registered backend.
ImportError – If an explicitly requested backend is unavailable.
RuntimeError – If runtime validation or coefficient fitting fails, including for a singular system or during CUDA Graph capture.

Notes

The field has global support. Unlike compact Shepard morphing, every control generally influences every mesh point. Attached fields are treated as Lagrangian data and are not pushed forward. Geometry-dependent caches are invalidated and topology caches are retained. The operation does not detect or repair inverted, degenerate, or self-intersecting cells. Call validate() explicitly when needed. Coefficient fitting is not supported inside CUDA Graph capture because the singular-system check requires host interaction.

physicsnemo.mesh.transformations.deform.shrinkwrap( mesh: Mesh, target: Mesh, *, offset: float | Float[torch.Tensor, ''] = 0.0, max_distance: float | None = None, point_weights: str | tuple[str, ...] | Bool[torch.Tensor, ' n_points'] | Float[torch.Tensor, ' n_points'] | None = None, implementation: Literal['torch', 'warp'] | None = None, ) → Mesh[source]#

Project mesh vertices onto the nearest locations of a target surface.

For each source vertex \(x_i\), Shrinkwrap selects the closest point \(p_i\) on target and returns

\[x'_i = x_i + w_i\left(p_i + \delta n_i - x_i\right),\]

where \(w_i\) is an optional point weight, \(\delta\) is offset, and \(n_i\) is the oriented unit normal of the selected target face. Source connectivity and attached data are retained.

Call it as shrinkwrap(mesh, target, ...) or as mesh.shrinkwrap(target, ...). Neither source nor target is modified.

Parameters:

mesh (Mesh) – Source mesh whose 3D point coordinates are projected. Its connectivity may have any supported topology.
target (Mesh) – Nonempty triangle surface embedded in 3D. Source and target points must share their float32 or float64 dtype and device. Target connectivity may use any non-bool integer dtype supported by Mesh and is normalized to int64 for this operation.
offset (float or torch.Tensor, optional) – Signed scalar distance from the target along the selected face normal. A scalar tensor must match the mesh point dtype and device and may require gradients. The value must be finite. Positive values follow target face winding. Default is 0.0.
max_distance (float or None, optional) – Positive finite nearest-surface search radius in mesh coordinate units. Vertices without a target strictly closer than this distance remain unchanged. Values that round to zero in the mesh point dtype are rejected. None performs an unbounded search. Default is None.
point_weights (str, tuple[str, ...], torch.Tensor, or None, optional) – Optional bool or floating per-source-point weights, or a key/path in point_data resolving to one. Zero leaves a vertex unchanged and one applies the full projection. Floating values are not clamped. Default is None.
implementation ({"torch", "warp"} or None, optional) – Nearest-face backend. None selects Torch on CPU and Warp on CUDA when available.

Returns:

New mesh with projected points, preserved connectivity and data, invalidated geometry caches, and retained topology caches.

Return type:

Mesh

Notes

Nearest-face selection is discrete. Away from target face, edge, vertex, and distance-cutoff transitions, gradients propagate through source points, selected target vertices, floating point weights, and tensor-valued offset. Degenerate target faces are ignored. Target orientation must be consistent when using a nonzero offset.

Shrinkwrap does not detect or prevent inverted or self-intersecting source cells. Partial point weights interpolate toward the projected surface and generally do not place the resulting vertex exactly on it. Shrinkwrap is not supported inside CUDA Graph capture with either backend. Float64 targets use the Torch search because Warp searches in float32. Safe float32 coordinates are searched unchanged. Warp falls back to Torch for unsafe coordinate magnitudes or face geometry.

physicsnemo.mesh.transformations.deform.sobolev_deform( mesh: Mesh, displacement: str | tuple[str, ...] | Float[torch.Tensor, 'n_points n_spatial_dims'], *, length_scale: float, fixed_points: str | tuple[str, ...] | Bool[torch.Tensor, ' n_points'] | None = None, max_iterations: int = 128, tolerance: float | None = None, implementation: Literal['torch', 'warp'] | None = None, ) → Mesh[source]#

Deform a mesh with a uniform-mass P1 Sobolev displacement.

Filters a dense per-vertex displacement by solving

\[(M + \ell^2 K)u = M d\]

and returns a new mesh with points mesh.points + u. Here \(M=\bar m I\) is a uniform vertex mass matrix scaled by the mean positive lumped P1 mass. \(K\) is the P1 stiffness matrix, and \(\ell\) is length_scale in mesh coordinate units. The uniform mass makes the filter self-adjoint for PyTorch vertex tensors. Connectivity and attached fields are unchanged.

Call it as sobolev_deform(mesh, ...) or as mesh.sobolev_deform(...). The bound method supplies mesh automatically.

Parameters:

mesh (Mesh) – Simplicial mesh whose points are deformed. The source is not modified.
displacement (str, tuple[str, ...], or torch.Tensor) – Raw displacement with shape (mesh.n_points, mesh.n_spatial_dims), or a point_data key resolving to one. It must match the point dtype and device.
length_scale (float) – Nonnegative physical smoothing length. Zero applies the raw displacement directly at unfixed points.
fixed_points (str, tuple[str, ...], torch.Tensor, or None, optional) – Optional bool mask with shape (mesh.n_points,), or a point-data key resolving to one. True entries receive zero displacement. Default is None.
max_iterations (int, optional) – Maximum PCG iterations. Default is 128.
tolerance (float or None, optional) – Positive relative residual tolerance. None selects a dtype-dependent default. Default is None.
implementation ({"torch", "warp"} or None, optional) – Backend override. None selects Torch on CPU. On CUDA, it selects Warp for segments, triangles, and tetrahedra when available. It otherwise selects Torch, with a one-time RuntimeWarning when Warp is unavailable. The Warp backend requires CUDA tensors.

Returns:

New mesh with Sobolev-filtered points, unchanged connectivity, and unchanged attached fields.

Return type:

Mesh

Notes

Unfixed mesh boundaries use the natural homogeneous Neumann condition. Constant displacements are retained when no points are fixed. Isolated points receive their raw displacement.

Both backends participate in autograd through the source points and the raw displacement. Their reverse-mode derivatives solve the adjoint Helmholtz system, which makes the operation suitable for smooth vertex-based optimization. The Warp backend evaluates the geometry vector-Jacobian product analytically. A forward or adjoint solve that does not reach tolerance within max_iterations raises a RuntimeError. Warp supports segment, triangle, and tetrahedron cells. Higher-dimensional simplices use Torch by default. Warp CUDA results and point gradients may vary at roundoff between runs. CUDA Graph capture is not supported because P1 operator assembly and solver diagnostics are not capture-safe.

Geometry caches are invalidated and topology caches are retained. At positive length scales, cells must be finite, nondegenerate simplices. The operation does not detect inverted or self-intersecting output cells.

Projections#

Spatial dimension manipulation – changing the embedding dimension of a mesh without altering its manifold dimension.

embed() – add spatial dimensions (non-destructive; for example, 2D mesh to 3D by appending zero coordinates)
extrude() – sweep a manifold to create a mesh one dimension higher (for example, a triangle mesh extruded to a prism mesh)
project() – reduce spatial dimensions (lossy; drops coordinate axes)

Projection operations for mesh extrusion, embedding, and spatial dimension manipulation.

This module provides functionality for: - Embedding meshes in higher-dimensional spaces (non-destructive) - Projecting meshes to lower-dimensional spaces (lossy) - Extruding manifolds to higher dimensions