API reference ¶

Raises:

ValueError – If direction is out of range [0, dim).
ValueError – If the degree in the given direction is 0.

Example

>>> # First derivative of a 1D curve
>>> f_prime = f.derivative()
>>> # Partial derivative of a surface with respect to direction 1
>>> df_dv = surface.derivative(direction=1)
>>> # Second derivative (composable)
>>> f_double_prime = f.derivative().derivative()
>>> # Derivative preserving degree
>>> f_prime_same_deg = f.derivative(keep_degree=True)

elevate_degree(degree_increments)[source]¶

Elevate the polynomial degree of the B-spline.

Creates a new B-spline that represents the same mapping as the original one but with higher-order polynomial basis functions. This is achieved by increasing the degree in each parametric direction and adjusting the control points and knot vectors accordingly.

Parameters:

degree_increments (int | Sequence[int]) – Number of degrees to increase. If an integer, the same increment is applied to all parametric directions. If a sequence, must have length equal to the B-spline dimension.

Returns:

A new B-spline with elevated degrees.

Return type:

Raises:

ValueError – If any degree increment is negative.
ValueError – If all degree increments are zero.
ValueError – If the number of increments does not match the dimension.

References

Degree elevation of spline curves [Piegl and Tiller, 1997].

reduce_degree(degree_decrements)[source]¶

Reduce the polynomial degree of the B-spline via least-squares approximation.

Decomposes to Bézier segments, reduces each segment using bidiagonal least-squares (Givens QR), and coarsens the knot vector to restore the original continuity structure.

Unlike elevate_degree(), this operation is not exact in general: the result is an approximation of the original mapping.

Parameters:

degree_decrements (int | Sequence[int]) – Number of degrees to reduce. If an integer, the same decrement is applied to all parametric directions. If a sequence, must have length equal to the B-spline dimension.

Returns:

A new B-spline with reduced degrees.

Return type:

Raises:

ValueError – If any degree decrement is negative.
ValueError – If all degree decrements are zero.
ValueError – If the number of decrements does not match the dimension.
ValueError – If any decrement exceeds the current degree in that direction.

insert_knots(new_knots)[source]¶

Return a geometrically equivalent B-spline with additional knots inserted.

Parameters:

new_knots (npt.ArrayLike | Sequence[npt.ArrayLike | None]) – For a 1D B-spline, a flat non-empty 1D array-like of knot values to insert. For multi-dimensional B-splines, a sequence of length dim where each element is a 1D array-like of knots to insert in that direction, or None to skip that direction. At least one direction must have a non-empty array of knots to insert. Repeated values in an array insert the same knot multiple times.

Returns:

New B-spline with the same geometry and refined knot vectors.

Return type:

Raises:

ValueError – If the sequence length does not match dim (multi-dim case).
ValueError – If all directions have empty or None knot arrays.
ValueError – If any knot lies outside its direction’s domain.
ValueError – If any insertion would exceed maximum multiplicity.

References

Knot insertion and refinement [Piegl and Tiller, 1997].

remove_knots(knot_values, *, num=None, tol=None)[source]¶

Return a B-spline with specified interior knots removed.

Each listed knot value is removed up to num times (or as many as possible when num=None), provided the geometric deviation stays within tol.

Parameters:

knot_values (float | npt.ArrayLike | Sequence[npt.ArrayLike | None]) – For a 1D B-spline, a single float or a 1D array-like of distinct interior knot values to remove. For multi-dimensional B-splines, a sequence of length dim where each element is a 1D array-like of knot values to remove in that direction, or None to skip that direction. At least one direction must have a non-empty array of knot values.
num (int | None) – Maximum number of removals per distinct knot value. None (default) removes as many as possible (up to the current multiplicity, capped at the degree).
tol (float | None) – Maximum allowed geometric deviation for each removal step. None (default) uses 1e-10.

Returns:

New B-spline with the same geometry (within tolerance) and reduced knot vectors.

Return type:

Raises:

ValueError – If the B-spline is periodic in any direction.
ValueError – If the sequence length does not match dim (multi-dim case).
ValueError – If all directions have empty or None knot arrays.
ValueError – If any knot value is not found or is a boundary knot.

to_open_bspline()[source]¶

Return an open (clamped) non-periodic B-spline equivalent to this one.

Converts each parametric direction to an open representation by inserting knots at the domain boundaries until each has multiplicity degree + 1, then trimming any ghost knots outside the domain. Works for 1D and multi-dimensional B-splines, and handles periodic, unclamped non-periodic, and mixed cases.

For periodic splines the n_full = len(knots) - degree - 1 control points are reconstructed by modulo-wrapping the n_periodic stored control points (ctrl_full[i] = ctrl[i % n_periodic]), and the Oslo algorithm is applied to this full set. The resulting open B-spline represents the mathematical periodic function defined by the periodic knot vector and the wrapped control points.

Returns:: Open, non-periodic B-spline with clamped knot vectors.
Return type:: Bspline
Raises:: ValueError – If the B-spline is already open in every direction.

restrict(bounds)[source]¶

Return a B-spline restricted to a sub-region of the parametric domain.

Extracts the portion of the B-spline defined on the given sub-domain by inserting knots at the new boundaries until they reach multiplicity degree + 1, then extracting the relevant knot sub-vector and control points. Skips insertion when a bound already coincides with a clamped domain endpoint. Periodic directions are automatically converted to open form before restricting.

Parameters:

bounds (tuple[float, float] | Sequence[tuple[float, float] | None]) – For a 1D B-spline, a (lower, upper) tuple defining the new domain. For multi-dimensional B-splines, a sequence of length dim where each element is a (lower, upper) tuple for that direction, or None to keep the full domain in that direction. At least one direction must have non-None bounds that actually restrict the domain.

Returns:

New B-spline defined on the restricted domain.

Return type:

Raises:

ValueError – If the sequence length does not match dim (multi-dim).
ValueError – If all directions are None or match the full domain.
ValueError – If any bound lies outside its direction’s domain.
ValueError – If lower >= upper in any direction.

split(direction, value)[source]¶

Split the B-spline into two at a parameter value in one direction.

Inserts knots at value until the multiplicity reaches degree + 1, then extracts the left and right sub-splines. Periodic directions are automatically converted to open form before splitting.

Parameters:

direction (int) – Parametric direction along which to split. Must be in [0, dim).
value (float) – Parameter value at which to split. Must lie strictly inside the domain of the specified direction.

Returns:

A pair (left, right) of B-splines. The left sub-spline has domain [domain_start, value] and the right has domain [value, domain_end] in the split direction. Other directions are unchanged.

Return type:

tuple[Bspline, Bspline]

Raises:

ValueError – If direction is out of range [0, dim).
ValueError – If value is not strictly inside the domain.

Example

>>> left, right = spline.split(0, 0.5)
>>> left, right = surface.split(1, 0.3)

to_periodic(continuity=None)[source]¶

Return a periodic B-spline equivalent to this one.

Converts each non-periodic parametric direction to a periodic representation with the requested continuity at the seam (the domain boundary where the periodic function wraps around). Directions that are already periodic are left unchanged.

The conversion performs an exact change of basis: the Oslo algorithm defines the linear map from periodic control points (with modulo wrapping) to open control points, and inverting this map via QR factorization recovers the periodic representation.

Warning

This method assumes the B-spline represents an inherently periodic function — i.e., the function values and derivatives match at the domain boundaries up to the requested continuity order. If the function is not periodic, the result will not represent the same geometry. Gross violations raise ValueError (via a residual check), but small deviations may pass silently. Use to_open_bspline() on the result to verify round-trip fidelity if in doubt.

Parameters:

continuity (int | tuple[int | None, ...] | None) – Target continuity at the seam per direction. None (default) requests maximum regularity C^{p-1} in every non-periodic direction. An integer applies to all non-periodic directions; a tuple specifies per-direction values where None entries skip that direction (leave it unchanged). Integer values must satisfy 0 <= continuity <= degree - 1.

Returns:

B-spline with the requested directions made periodic.

Return type:

Raises:

ValueError – If no direction would be converted (all already periodic or all skipped via None in the tuple).
ValueError – If the function is not periodic (endpoint mismatch or residual exceeds tolerance).
ValueError – If continuity is out of range.

to_bezier(*, copy=True)[source]¶

Convert to an equivalent Bézier.

Extracts a Bezier from this B-spline when it represents a single polynomial patch (one element per direction). Periodic and non-open B-splines are automatically converted to open form before the check.

Parameters:: copy (bool) – If True (default), the control points are deep-copied into the new Bézier. If False, the Bézier shares the same underlying control point array when possible (direct extraction) or owns the freshly allocated array produced by the open-form conversion.
Returns:: Equivalent Bézier representation.
Return type:: Bezier
Raises:: ValueError – If the B-spline has more than one element per direction and cannot be represented as a single Bézier.

to_beziers()[source]¶

Decompose into Bézier patches.

Decomposes the B-spline into its constituent Bézier patches by applying Bézier extraction operators direction by direction. Periodic directions are automatically converted to open form first.

The result is cached: the first call computes and stores the decomposition; subsequent calls return the cached array. In-place mutations (reverse, permute_directions, transform with in_place=True) invalidate the cache.

Returns:: Array of Bezier objects with shape (*num_intervals) following the tensor-product interval structure. For a 1D B-spline with 3 intervals the shape is (3,). For a 2D surface with 3x2 intervals the shape is (3, 2). The Bézier control points are read-only.
Return type:: npt.NDArray[np.object_]

Example

>>> beziers = spline.to_beziers()
>>> beziers.shape
(3,)
>>> beziers[0].degree
(2,)

References

Bézier extraction of NURBS [Borden et al., 2011] (and its T-spline generalization [Scott et al., 2011]).

multiply(other)[source]¶

Return the exact pointwise product of this B-spline and another.

Given B-splines self and other over the same parametric domain, returns a new B-spline h such that h(t) = self(t) * other(t) for all t in the domain. The result lives in the product space of degree p_d + q_d per direction where p_d and q_d are the degrees of the two operands in direction d.

Both non-rational and rational (NURBS) operands are supported. A non-rational operand is promoted to rational (unit weights) when the other is rational.

Parameters:

other (Bspline) – The second B-spline operand. Must have the same dimension, dtype, rank, and parametric domain as self.

Returns:

A new B-spline representing self * other.

Return type:

Raises:

ValueError – If the operands have different dimensions.
ValueError – If the operands have different dtypes.
ValueError – If the operands have different ranks.
ValueError – If the operands have different parametric domains.

Note

The boundary structure of the operands is preserved per direction: both periodic → periodic, both non-open → non-open, either open → open.

Example

>>> h = f.multiply(g)
>>> h2 = f * g  # equivalent via __mul__

reverse(direction: int = 0, *, in_place: Literal[False] = False) → Bspline[source]¶

reverse(direction: int = 0, *, in_place: Literal[True]) → None

Reverse the orientation of one parametric direction.

Flips the control points along the given parametric axis and reflects the corresponding knot vector so that the mapping is reparametrised in the opposite sense along that direction.

Parameters:

direction (int) – Parametric direction to reverse. Must be in [0, dim). Defaults to 0.
in_place (bool) – If True, modify this B-spline in place and return None. If False (default), return a new B-spline.

Returns:

The reversed B-spline, or None when in_place=True.

Return type:

Bspline | None

Raises:

ValueError – If direction is out of range [0, dim).

Example

>>> rev = spline.reverse(direction=0)
>>> spline.reverse(direction=1, in_place=True)

permute_directions(permutation: Sequence[int], *, in_place: Literal[False] = False) → Bspline[source]¶

permute_directions(permutation: Sequence[int], *, in_place: Literal[True]) → None

Reorder the parametric directions according to a permutation.

Given a permutation [i_0, i_1, …], the new direction k is the old direction permutation[k]. For example, [1, 2, 0] on a 3D volume maps old direction 1 → new 0, old 2 → new 1, old 0 → new 2.

Parameters:

permutation (Sequence[int]) – A permutation of range(dim).
in_place (bool) – If True, modify this B-spline in place and return None. If False (default), return a new B-spline.

Returns:

The permuted B-spline, or None when in_place=True.

Return type:

Bspline | None

Raises:

ValueError – If permutation is not a valid permutation of range(dim).

Example

>>> surface.permute_directions([1, 0])  # swap u ↔ v

transform(affine: AffineTransform, *, in_place: Literal[False] = False) → Bspline[source]¶

transform(affine: AffineTransform, *, in_place: Literal[True]) → None

Apply an affine transformation to the control points.

For non-rational B-splines, every control point P is mapped to A @ P + b. For rational B-splines (NURBS) the weighted homogeneous coordinates are updated so that the projected geometry undergoes the same affine map while the weights are preserved.

Parameters:

affine (AffineTransform) – The affine transformation to apply.
in_place (bool) – If True, the control points are modified in place and None is returned. If False (default), a new Bspline is returned.

Returns:

The transformed B-spline, or None when in_place=True.

Return type:

Bspline | None

Raises:

ValueError – If the transform dimension does not match the geometric rank of the B-spline.

Example

>>> from pantr.transform import AffineTransform
>>> T = AffineTransform.translation([1.0, 2.0])
>>> shifted = spline.transform(T)

subdivide(n_subdivisions, regularity=None)[source]¶

Return a geometrically equivalent B-spline with uniformly refined knot vectors.

For every non-zero knot span in each active parametric direction, inserts n_subdivisions - 1 uniformly spaced knot values. Each value is repeated degree - regularity times so that the B-spline has C^regularity continuity at every inserted knot.

Parameters:

n_subdivisions (int | Sequence[int | None]) – Number of equal sub-spans per existing interval. A single int is applied to all directions; must be >= 2. A sequence of length dim provides per-direction counts; use None to skip a direction. At least one direction must have a count >= 2.
regularity (int | None) – Continuity order at every inserted knot. Applied uniformly across all active directions. Must be in [-1, degree - 1] for each active direction. None (default) uses degree - 1 per direction.

Returns:

New B-spline with refined knot vectors and same geometry.

Return type:

THBSplineSpaceRestriction

Raises:

ValueError – If the sequence length does not match dim.
ValueError – If any subdivision count is < 1.
ValueError – If no direction has a count >= 2.
ValueError – If regularity is outside the valid range for any active direction.

slice(axis, value)[source]¶

Slice the B-spline by fixing one parametric direction at a given value.

Reduces the parametric dimension by one using de Boor corner cutting on the control points. A volume becomes a surface, a surface becomes a curve, and a curve becomes a point (returned as a NumPy array).

When the parameter value coincides with a knot of multiplicity s, only p - s de Boor iterations are needed. At a C0 knot (s >= p) the result is obtained in O(1) by direct control point lookup.

Parameters:

axis (int) – Parametric direction to fix (0-indexed). Must be in [0, dim).
value (float) – Parameter value at which to slice. Must lie within the domain of the specified direction.

Returns:

A B-spline with dim - 1 dimensions when dim >= 2, or a NumPy array of shape (rank,) when dim == 1. Rational B-splines preserve the NURBS structure when dim >= 2; for dim == 1 the result is projected to physical coordinates.

Return type:

Bspline | npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If axis is out of range [0, dim).
ValueError – If value is outside the domain of the specified direction.

Example

>>> # Slice a surface at v=0.5 to get a curve
>>> curve = surface.slice(1, 0.5)
>>> # Slice a volume at w=0.3 to get a surface
>>> srf = volume.slice(2, 0.3)
>>> # Composable: volume -> surface -> curve -> point
>>> pt = volume.slice(2, 0.3).slice(1, 0.5).slice(0, 0.2)

boundary(axis, side)[source]¶

Extract the boundary of the B-spline along one parametric direction.

Returns the restriction of the B-spline to one end of the domain in the given direction.

Parameters:

axis (int) – Parametric direction (0-indexed). Must be in [0, dim).
side (int) – Which end of the domain: 0 for the start, 1 for the end.

Returns:

A B-spline with dim - 1 dimensions when dim >= 2, or a NumPy array of shape (rank,) when dim == 1.

Return type:

Bspline | npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If axis is out of range [0, dim).
ValueError – If side is not 0 or 1.

Example

>>> # Extract left boundary of a surface along direction 0
>>> left_curve = surface.boundary(0, 0)
>>> # Extract right boundary along direction 1
>>> right_curve = surface.boundary(1, 1)

plot(*, color=None, show_control_polygon=False, show_knot_lines=False, **plotter_kwargs)[source]¶

Quick interactive visualization of this B-spline (requires pyvista).

For finer control, use pantr.viz.Scene directly.

Parameters:

color (str | None) – Surface color.
show_control_polygon (bool) – Render control polygon (points and wireframe).
show_knot_lines (bool) – Render knot lines.
**plotter_kwargs (Any) – Additional keyword arguments for pv.Plotter().

Returns:

The pyvista Plotter after showing.

Return type:

object

Raises:

ImportError – If pyvista is not installed.

class pantr.bspline.BsplineSpace(spaces)[source]¶

Bases: object

A class representing a multi-dimensional B-spline space.

This space is defined by a set of B-spline spaces, one for each dimension.

This class provides methods to analyze B-spline properties, validate input parameters, compute various geometric characteristics of the spline, and access various properties of the B-spline.

Parameters:: spaces (Iterable[BsplineSpace1D])

_spaces¶

List of B-spline spaces, one for each dimension.

Type:: Iterable[BsplineSpace1D]

__init__(spaces)[source]¶

Initialize a B-spline space object.

Parameters:: spaces (Iterable[BsplineSpace1D]) – List of B-spline spaces, one for each dimension.
Raises:: ValueError – If the B-spline spaces have different data types.
Return type:: None

property dim: int¶

Get the dimension of the B-spline space.

Returns:: The dimension of the B-spline space.
Return type:: int

property spaces: tuple[BsplineSpace1D, ...]¶

Get the B-spline spaces.

Returns:: The B-spline spaces.
Return type:: tuple[BsplineSpace1D, …]

property degrees: tuple[int, ...]¶

Get the polynomial degree of the B-spline.

Returns:: The degree for each dimension.
Return type:: tuple[int, …]

property tolerance: float¶

Get the tolerance value used for numerical comparisons.

It is the maximum tolerance of the B-spline spaces.

Returns:: The tolerance value.
Return type:: float

property dtype: DTypeLike¶

Get the data type of the B-spline space.

Returns:: The numpy data type of the B-spline space.
Return type:: npt.DTypeLike

property num_basis: tuple[int, ...]¶

Get the number of basis functions for each dimension.

Returns:: The number of basis functions for each dimension.
Return type:: tuple[int, …]

property num_total_basis: int¶

Get the total number of basis functions.

Returns:: The total number of basis functions.
Return type:: int

property num_intervals: tuple[int, ...]¶

Get the number of intervals for each dimension.

Returns:: The number of intervals for each dimension.
Return type:: tuple[int, …]

property num_total_intervals: int¶

Get the total number of intervals.

Returns:: The total number of intervals.
Return type:: int

property domain: ndarray[tuple[Any, ...], dtype[float32 | float64]]¶

Get the domain of the B-spline space.

Returns:: The domain of the B-spline space. The shape is (dim, 2), where the last dimension contains the start and end values of the domain.
Return type:: npt.NDArray[np.float32 | np.float64]

has_Bezier_like_knots()[source]¶

Check if the knot vector represents a Bézier-like configuration.

A Bézier-like configuration has open ends and only one non-zero span for each dimension.

Returns:: True if knots have open ends and only one span.
Return type:: bool

Example

>>> bspline_1D = BsplineSpace1D([1, 1, 1, 3, 3, 3], 2)
>>> bspline_2D = BsplineSpace([bspline_1D, bspline_1D])
>>> bspline_2D.has_Bezier_like_knots()
True

tabulate_basis(pts, out_basis=None, out_first_basis=None)[source]¶

Tabulate the B-spline basis functions at the given points.

Parameters:

pts (npt.NDArray[np.float32 | np.float64] | PointsLattice) – The points at which to tabulate the basis functions. It can be a 2D array with shape (num_pts, dim) or a PointsLattice object.
out_basis (npt.NDArray[np.float32 | np.float64] | None) – Optional output array where the basis values will be stored. If None, a new array is allocated. Must have the correct shape and dtype if provided. This follows NumPy’s style for output arrays. Defaults to None.
out_first_basis (npt.NDArray[numpy.intp] | None) – Optional output array where the first basis indices will be stored. If None, a new array is allocated. Must have the correct shape and dtype numpy.intp if provided. This follows NumPy’s style for output arrays. Defaults to None.

Returns:

The basis function values and the first basis function indices.

In the case pts is a 2D array, the shape of the basis function values array is (num_pts, order[0], order[1], …, order[d-1]), where d is the dimension of the B-spline space and num_pts is the number of points. In the case pts is a PointsLattice object, the shape of the basis function values array is (num_pts_0, num_pts_1, …, num_pts_d, order[0], order[1], …, order[d-1]), where num_pts_i is the number of points in the i-th dimension.

The shape of the first basis function indices array is (num_pts, dim), if pts is a 2D array, or (num_pts_0, num_pts_1, …, num_pts_d, dim), if pts is a PointsLattice object.

If out_basis or out_first_basis was provided, the corresponding element of the tuple is the same array.

Return type:

tuple[npt.NDArray[np.float32 | np.float64], npt.NDArray[numpy.intp]]

Raises:

ValueError – If pts is not a 2D array or a PointsLattice object.
ValueError – If the pts dimension does not match the dimension of the B-spline space.
ValueError – If one or more points are outside the domain of the B-spline space, or if out_basis or out_first_basis is provided and has incorrect shape or dtype.

restrict(cell_ids)[source]¶

Return the bounding-box windowed sub-space spanning cell_ids.

The window is the per-axis multi-index bounding box of the requested knot-span cells (flat ids in C-order over num_intervals, the same convention as pantr.grid.tensor_product_grid() and SpanwiseElementExtraction). Each axis is windowed by slicing its knot vector (never re-clamped), so the windowed basis equals this space’s basis pointwise over the windowed cells.

Parameters:

cell_ids (npt.ArrayLike) – Flat knot-span cell ids to span; duplicates are ignored. Each must satisfy 0 <= cid < num_total_intervals.

Returns:

The windowed BsplineSpace and the read-only local_to_global_dof map of shape (windowed_space.num_total_basis,).

Return type:

BsplineSpaceRestriction

Raises:

ValueError – If cell_ids is empty or any axis is periodic.
IndexError – If any cell id is out of range [0, num_total_intervals).
TypeError – If cell_ids is not integer-valued.

class pantr.bspline.BsplineSpace1D(knots, degree, periodic=False, snap_knots=True)[source]¶

Bases: object

A class representing a 1D B-spline with configurable degree and knot vector.

This class provides methods to analyze B-spline properties, validate input parameters, compute various geometric characteristics of the spline, and access various properties of the B-spline.

Parameters:

knots (npt.ArrayLike)
degree (int)
periodic (bool)
snap_knots (bool | None)

_tol¶

Tolerance value for numerical comparisons.

Type:: np.float32 | np.float64

_knots¶

Knot vector defining the B-spline.

Type:: npt.NDArray[np.float32 | np.float64]

_degree¶

Polynomial degree of the B-spline.

Type:: int

_periodic¶

Whether the B-spline is periodic.

Type:: bool

__init__(knots, degree, periodic=False, snap_knots=True)[source]¶

Initialize a B-spline 1D object.

Parameters:

knots (npt.ArrayLike) – Knot vector defining the B-spline. Must be non-decreasing and have at least 2*degree+2 elements.
degree (int) – Polynomial degree of the B-spline. Must be non-negative.
periodic (bool) – Whether the B-spline is periodic. Defaults to False.
snap_knots (bool | None) – Whether to snap nearby knots to avoid numerical issues. Defaults to True.

Raises:

ValueError – If degree is negative, knots are insufficient, or knots are not non-decreasing.
TypeError – If knots cannot be converted to a numpy array.

Return type:

None

property degree: int¶

Get the polynomial degree of the B-spline.

Returns:: The degree.
Return type:: int

property knots: ndarray[tuple[Any, ...], dtype[float32 | float64]]¶

Get the knot vector.

Returns:: The knot vector.
Return type:: npt.NDArray[np.float32 | np.float64]

property periodic: bool¶

Check if the B-spline is periodic.

Returns:: True if periodic, False otherwise.
Return type:: bool

property tolerance: float¶

Get the tolerance value used for numerical comparisons.

Returns:: The tolerance value.
Return type:: float

property dtype: DTypeLike¶

Get the data type of the knot vector (and used in computations).

Returns:: The numpy data type of the knots.
Return type:: npt.DTypeLike

property num_basis: int¶

Get the number of basis functions.

This depends on the knot vector length and the degree, but also on whether the B-spline is periodic.

Returns:: Number of basis functions.
Return type:: int

get_unique_knots_and_multiplicity(in_domain=False)[source]¶

Get unique knots and their multiplicities.

Parameters:: in_domain (bool) – If True, only consider knots in the domain. Defaults to False.
Returns:: Tuple of (unique_knots, multiplicities) where unique_knots contains the distinct knot values and multiplicities contains the corresponding multiplicity counts.
Return type:: tuple[npt.NDArray[np.float32 | np.float64], npt.NDArray[numpy.intp]]

property num_intervals: int¶

Get the number of intervals in the domain.

Returns:: Number of intervals.
Return type:: int

Example

>>> space = BsplineSpace1D([0, 0, 0, 1, 2, 2, 2], 2)
>>> space.num_intervals
2

property domain: tuple[float32 | float64, float32 | float64]¶

Get the knot vector domain.

Returns:: Tuple of (start_value, end_value) defining the domain.
Return type:: tuple[np.float32 | np.float64, np.float64]

Example

>>> bspline = BsplineSpace1D([0, 0, 0, 1, 2, 2, 2], 2)
>>> bspline.domain
(0.0, 2.0)

has_left_end_open()[source]¶

Check if the left end of the B-spline is open.

A left end is open if the first degree+1 knots are equal.

Returns:: True if the left end is open, False otherwise.
Return type:: bool

has_right_end_open()[source]¶

Check if the right end of the B-spline is open.

A right end is open if the last degree+1 knots are equal.

Returns:: True if the right end is open, False otherwise.
Return type:: bool

has_open_knots()[source]¶

Check if the B-spline has open ends.

Returns:: True if both ends are open, False otherwise.
Return type:: bool

has_Bezier_like_knots()[source]¶

Check if the knot vector represents a Bézier-like configuration.

A Bézier-like configuration has open ends and only one non-zero span.

Returns:: True if knots have open ends and only one span.
Return type:: bool

Example

>>> bspline = BsplineSpace1D([1, 1, 1, 3, 3, 3], 2)
>>> bspline.has_Bezier_like_knots()
True

get_cardinal_intervals(out=None)[source]¶

Get boolean array indicating whether the intervals (non-zero spans) are cardinal or not.

An interval is cardinal if has the same length as the degree-1 previous and the degree-1 next intervals.

In the case of open knot vectors, this definition automatically discards the first degree-1 and the last degree-1 intervals.

Parameters:

out (npt.NDArray[bool] | None) – Optional output array where the result will be stored. If None, a new array is allocated. Must have the correct shape and dtype if provided. This follows NumPy’s style for output arrays. Defaults to None.

Returns:

Boolean array where True indicates cardinal intervals.: It has length equal to the number of intervals. If out was provided, returns the same array.

Return type:

npt.NDArray[bool]

Raises:

ValueError – If out is provided and has incorrect shape or dtype.

Example

>>> bspline = BsplineSpace1D([0, 0, 0, 1, 2, 3, 4, 5, 6, 6, 6], 2)
>>> bspline.get_cardinal_intervals()
array([False, False, True, True, False, False])

>>> bspline = BsplineSpace1D([0, 0, 0, 1, 2, 3, 4, 5, 5, 6, 6, 6], 2)
>>> bspline.get_cardinal_intervals()
array([False, False, True, False, False, False])

>>> bspline = BsplineSpace1D([0, 1, 2, 3, 4, 5, 6, 7, 10], 3)
>>> bspline.get_cardinal_intervals()
array([True, False])

tabulate_basis(pts, out_basis=None, out_first_basis=None, *, validate=True)[source]¶

Evaluate the B-spline basis functions at the given points.

Parameters:

pts (npt.ArrayLike) – Evaluation points.
out_basis (npt.NDArray[np.float32 | np.float64] | None) – Optional output array where the basis values will be stored. If None, a new array is allocated. Must have the correct shape and dtype if provided. This follows NumPy’s style for output arrays. Defaults to None.
out_first_basis (npt.NDArray[numpy.intp] | None) – Optional output array where the first basis indices will be stored. If None, a new array is allocated. Must have the correct shape and dtype numpy.intp if provided. This follows NumPy’s style for output arrays. Defaults to None.
validate (bool) – If True (default), check that every point lies inside the spline domain. Pass False only when the caller guarantees in-domain points; out-of-domain points are then undefined behavior. Defaults to True.

Returns:

tuple[

npt.NDArray[np.float32] | npt.NDArray[np.float64], npt.NDArray[numpy.intp]

]: Tuple containing:

basis_values: (npt.NDArray[np.float32] | npt.NDArray[np.float64]) Array of shape matching pts with the last dimension length (degree+1), containing the basis function values evaluated at each point. If out_basis was provided, returns the same array.
first_basis_indices: (npt.NDArray[numpy.intp]) 1D integer array indicating the index of the first nonzero basis function for each evaluation point. The length is the same as the number of evaluation points. If out_first_basis was provided, returns the same array.

Raises:

ValueError – If validate is True and any evaluation point is outside the B-spline domain, or if out_basis or out_first_basis is provided and has incorrect shape or dtype.

Return type:

tuple[ndarray[tuple[Any, …], dtype[float32 | float64]], ndarray[tuple[Any, …], dtype[int64]]]

Example

>>> bspline = BsplineSpace1D([0, 0, 0, 0.25, 0.7, 0.7, 1, 1, 1], 2)
>>> bspline.tabulate_basis([0.0, 0.5, 0.75, 1.0])
(array([[1.        , 0.        , 0.        ],
        [0.12698413, 0.5643739 , 0.30864198],
        [0.69444444, 0.27777778, 0.02777778],
        [0.        , 0.        , 1.        ]]),
 array([0, 1, 3, 3]))

References

Basis values are computed with the stable Cox-de Boor recurrence [de Boor, 2001].

tabulate_basis_derivatives(pts, n_deriv, out_deriv=None, out_first_basis=None, *, validate=True)[source]¶

Evaluate B-spline basis function derivatives at the given points.

Implements Algorithm A2.3 (DerBasisFuncs) from Piegl & Tiller. The 0th slice deriv_values[..., 0, :] is identical to the result of tabulate_basis(). For n_deriv > degree all rows beyond degree are identically zero.

Parameters:

pts (npt.ArrayLike) – Evaluation points.
n_deriv (int) – Maximum derivative order to compute (>= 0).
out_deriv (npt.NDArray[np.float32 | np.float64] | None) – Optional output array for derivative values. If None, a new array is allocated. Must have shape (*pts_shape, n_deriv+1, degree+1) and dtype matching pts if provided. Defaults to None.
out_first_basis (npt.NDArray[numpy.intp] | None) – Optional output array for first basis indices. If None, a new array is allocated. Must have shape pts_shape and dtype numpy.intp if provided. Defaults to None.
validate (bool) – If True (default), check that every point lies inside the spline domain. Pass False only when the caller guarantees in-domain points; out-of-domain points are then undefined behavior. Defaults to True.

Returns:

tuple[

npt.NDArray[np.float32] | npt.NDArray[np.float64], npt.NDArray[numpy.intp]

]: Tuple containing:

deriv_values: Array of shape (*pts_shape, n_deriv+1, degree+1). deriv_values[..., k, i] is the k-th derivative of the i-th local B-spline basis function at each point.
first_basis_indices: Integer array of shape pts_shape giving the global index of the first nonzero basis function for each point.

Raises:

ValueError – If n_deriv < 0, if validate is True and any evaluation point is outside the domain, or out_deriv / out_first_basis has incorrect shape or dtype.

Return type:

tuple[ndarray[tuple[Any, …], dtype[float32 | float64]], ndarray[tuple[Any, …], dtype[int64]]]

Example

>>> bspline = BsplineSpace1D([0, 0, 0, 1, 1, 1], 2)
>>> d, first = bspline.tabulate_basis_derivatives([0.5], n_deriv=1)
>>> d.shape
(1, 2, 3)
>>> d[0, 1, :]   # first derivatives at x=0.5: B0'=-1, B1'=0, B2'=1
array([-1.,  0.,  1.])

tabulate_Bezier_extraction_operators(out=None)[source]¶

Create Bézier extraction operators of the B-spline.

Parameters:

out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array where the result will be stored. If None, a new array is allocated. Must have the correct shape and dtype if provided. This follows NumPy’s style for output arrays. Defaults to None.

Returns:

Array of extraction matrices with shape

(n_intervals, degree+1, degree+1) where each matrix transforms Bernstein basis functions to B-spline basis functions for that interval.

Each matrix C[i, :, :] transforms Bernstein basis functions to B-spline basis functions for the i-th interval as

C[i, :, :] @ [Bernstein values] = [B-spline values in interval].

If out was provided, returns the same array.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If out is provided and has incorrect shape or dtype.

tabulate_Lagrange_extraction_operators(lagrange_variant=LagrangeVariant.EQUISPACES, out=None)[source]¶

Create Lagrange extraction operators of the B-spline.

Parameters:

lagrange_variant (LagrangeVariant) – Lagrange point distribution to use (e.g., equispaced, Gauss-Lobatto-Legendre, etc). Defaults to LagrangeVariant.EQUISPACES.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array where the result will be stored. If None, a new array is allocated. Must have the correct shape and dtype if provided. This follows NumPy’s style for output arrays. Defaults to None.

Returns:

Array of extraction matrices with shape

(n_intervals, degree+1, degree+1) where each matrix transforms Lagrange basis functions to B-spline basis functions for that interval.

Each matrix C[i, :, :] transforms Lagrange basis functions to B-spline basis functions for the i-th interval as

C[i, :, :] @ [Lagrange values] = [B-spline values in interval].

If out was provided, returns the same array.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If out is provided and has incorrect shape or dtype.

tabulate_cardinal_extraction_operators(out=None)[source]¶

Create cardinal B-spline extraction operators of the B-spline.

Parameters:

Returns:

Array of extraction matrices with shape

(n_intervals, degree+1, degree+1) where each matrix transforms cardinal spline basis functions to B-spline basis functions for that interval.

Each matrix C[i, :, :] transforms cardinal spline basis functions to B-spline basis functions for the i-th interval as

C[i, :, :] @ [cardinal values] = [B-spline values in interval].

If out was provided, returns the same array.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If out is provided and has incorrect shape or dtype.

insert_knots(new_knots)[source]¶

Return a new BsplineSpace1D with additional knots inserted.

Parameters:

new_knots (npt.ArrayLike) – 1D array-like of knot values to insert. Must be non-empty. Values must lie in [knots[degree], knots[-degree-1]]. Inserting a value already present increases its multiplicity; multiplicity cannot exceed degree + 1. Repeated values in new_knots insert that knot multiple times in one call.

Returns:

New space with the inserted knots.

Return type:

BsplineSpace1D

Raises:

ValueError – If new_knots is empty.
ValueError – If new_knots is not 1D.
ValueError – If any new knot lies outside the domain.
ValueError – If inserting a knot would exceed maximum multiplicity degree + 1.

subdivide(n_subdivisions, regularity=None)[source]¶

Return a new BsplineSpace1D with each knot span split into equal sub-spans.

For every non-zero knot interval [t_i, t_{i+1}), inserts n_subdivisions - 1 uniformly spaced knot values. Each value is repeated degree - regularity times so that the resulting B-spline has C^regularity continuity at every inserted knot.

Parameters:

n_subdivisions (int) – Number of equal sub-spans per existing knot interval. Must be >= 2.
regularity (int | None) – Continuity order at every inserted knot. Must be in [-1, degree - 1]. None (default) uses degree - 1, giving maximum continuity (multiplicity 1).

Returns:

New space with uniformly refined knot vector.

Return type:

BsplineSpace1D

Raises:

ValueError – If n_subdivisions < 2.
ValueError – If regularity is outside [-1, degree - 1].

restrict(interval_lo, interval_hi)[source]¶

Return the windowed sub-space over the intervals [interval_lo, interval_hi).

The windowed knot vector is a pure slice of this space’s knots (never re-clamped), chosen so the windowed basis functions are exactly the global basis functions supported on the windowed intervals; they therefore equal the global basis pointwise over those intervals.

Parameters:

interval_lo (int) – First interval (knot span) of the window, inclusive.
interval_hi (int) – One past the last interval of the window, exclusive. Must satisfy 0 <= interval_lo < interval_hi <= num_intervals.

Returns:

The windowed space and a read-only local_to_global_dof array of shape (windowed_space.num_basis,) mapping each windowed basis index to its index in this space.

Return type:

tuple[BsplineSpace1D, npt.NDArray[np.int64]]

Raises:

ValueError – If the space is periodic, or the interval range is invalid.

class pantr.bspline.BsplineSpaceRestriction(space, local_to_global_dof)[source]¶

Bases: NamedTuple

Result of BsplineSpace.restrict(): a windowed space and its DOF map.

space – the windowed BsplineSpace: per axis a pure knot-vector slice of the parent (never re-clamped), so its basis equals the parent’s pointwise over the windowed cells.
local_to_global_dof – read-only, shape (space.num_total_basis,), mapping each windowed DOF (flat, C-order over the windowed per-axis basis counts) to its flat index in the parent space.

Unlike pantr.grid.GridRestriction there is no in_subset mask: every windowed DOF is a genuine parent DOF (a windowed space spans a box of cells, so there are no fill DOFs).

Parameters:

space (BsplineSpace)
local_to_global_dof (ndarray[tuple[Any, ...], dtype[int64]])

space: BsplineSpace¶: Alias for field number 0

local_to_global_dof: ndarray[tuple[Any, ...], dtype[int64]]¶: Alias for field number 1

class pantr.bspline.CouplingGraph(num_vertices, xadj, adjncy, edge_weights, vertex_weights)[source]¶

Bases: NamedTuple

Cell-coupling graph of a space, in METIS / Scotch CSR adjacency format.

A typing.NamedTuple returned by coupling_graph(). The graph is undirected (symmetric adjacency) and has no self-loops:

num_vertices – number of cells (graph vertices).
xadj – CSR row pointers, shape (num_vertices + 1,); the neighbors of cell c are adjncy[xadj[c]:xadj[c + 1]]. (METIS xadj.) Read-only.
adjncy – concatenated neighbor cell ids, shape (xadj[-1],). (METIS adjncy.) Read-only.
edge_weights – per-adjacency-entry weight = number of basis functions the two cells share, aligned with adjncy. (METIS adjwgt.) Read-only.
vertex_weights – per-cell weight (assembly cost), shape (num_vertices,); uniform 1.0 unless cell_weights was given. (METIS vwgt.) Read-only.

Parameters:

num_vertices (int)
xadj (npt.NDArray[np.int64])
adjncy (npt.NDArray[np.int64])
edge_weights (npt.NDArray[np.int64])
vertex_weights (npt.NDArray[np.float64])

num_vertices: int¶: Alias for field number 0

xadj: npt.NDArray[np.int64]¶: Alias for field number 1

adjncy: npt.NDArray[np.int64]¶: Alias for field number 2

edge_weights: npt.NDArray[np.int64]¶: Alias for field number 3

vertex_weights: npt.NDArray[np.float64]¶: Alias for field number 4

class pantr.bspline.ExtractionStructView(compact_ops_1d, idx_maps_1d, is_identity_mask_1d, num_intervals, input_shape_per_dir, output_shape_per_dir, dim)[source]¶

Bases: NamedTuple

Immutable struct view of a SpanwiseElementExtraction for @njit callers.

A typing.NamedTuple bundling the compact per-direction operator storage, index maps, identity masks, and shape metadata into a single object that Numba can unbox. Each array field is homogeneous in dtype and array dimensionality across all directions (per-direction shapes may differ), so Numba represents those tuple fields as UniTuple inside an @njit function. This makes ExtractionStructView a drop-in replacement for the separate (ops_1d, idx_maps_1d, is_identity_mask_1d, …) bundle when calling the Layer-3 batch kernels in pantr.bspline._extraction_kernels from downstream Numba code.

Construct via make_struct_view(). Field semantics mirror the same-named members of SpanwiseElementExtraction:

compact_ops_1d — per-direction compact 3D operator arrays of shape (n_compact_k, n_out_k, n_in_k); only non-identity rows are stored. Always has at least one row (sentinel zeros) to ensure safe Numba indexing.
idx_maps_1d — per-direction compact index maps of shape (n_elements_k,).
is_identity_mask_1d — per-direction identity masks of shape (n_elements_k,).
num_intervals — per-direction number of elements.
input_shape_per_dir — per-direction input sizes (n_in_0, …, n_in_{d-1}).
output_shape_per_dir — per-direction output sizes (n_out_0, …, n_out_{d-1}).
dim — number of tensor-product directions d.

Parameters:

compact_ops_1d (tuple[ndarray[tuple[Any, ...], dtype[float32 | float64]], ...])
idx_maps_1d (tuple[ndarray[tuple[Any, ...], dtype[int64]], ...])
is_identity_mask_1d (tuple[ndarray[tuple[Any, ...], dtype[bool]], ...])
num_intervals (tuple[int, ...])
input_shape_per_dir (tuple[int, ...])
output_shape_per_dir (tuple[int, ...])
dim (int)

compact_ops_1d: tuple[ndarray[tuple[Any, ...], dtype[float32 | float64]], ...]¶: Alias for field number 0

idx_maps_1d: tuple[ndarray[tuple[Any, ...], dtype[int64]], ...]¶: Alias for field number 1

is_identity_mask_1d: tuple[ndarray[tuple[Any, ...], dtype[bool]], ...]¶: Alias for field number 2

num_intervals: tuple[int, ...]¶: Alias for field number 3

input_shape_per_dir: tuple[int, ...]¶: Alias for field number 4

output_shape_per_dir: tuple[int, ...]¶: Alias for field number 5

dim: int¶: Alias for field number 6

class pantr.bspline.LocalSpace(space, local_to_global_cell, local_to_global_dof, owned_cell_mask, owned_dof_mask, n_global_cells, n_global_dofs)[source]¶

Bases: NamedTuple

The rank-local windowed view of a distributed B-spline or THB space.

Produced by build_local(). A typing.NamedTuple bundling a windowed BsplineSpace or THBSplineSpace with the maps and masks relating it to the global space:

space – the windowed space over the bounding box of the rank’s owned cells and their support-closure halo; a real pantr object whose basis equals the global basis pointwise over the rank’s owned cells.
local_to_global_cell – read-only map (one entry per local cell) from each local cell to its global cell id.
local_to_global_dof – read-only (space.num_total_basis,) map from each local DOF to its global DOF id (-1 for a THB boundary DOF with no global counterpart).
owned_cell_mask – read-only boolean mask (one entry per local cell); True for the cells this rank owns (the rest are halo / bounding-box fill).
owned_dof_mask – read-only boolean (space.num_total_basis,) mask; True for the DOFs this rank owns.
n_global_cells – total cells in the global space’s grid.
n_global_dofs – total basis functions in the global space.

Parameters:

space (BsplineSpace | THBSplineSpace)
local_to_global_cell (npt.NDArray[np.int64])
local_to_global_dof (npt.NDArray[np.int64])
owned_cell_mask (npt.NDArray[np.bool_])
owned_dof_mask (npt.NDArray[np.bool_])
n_global_cells (int)
n_global_dofs (int)

space: BsplineSpace | THBSplineSpace¶: Alias for field number 0

local_to_global_cell: npt.NDArray[np.int64]¶: Alias for field number 1

local_to_global_dof: npt.NDArray[np.int64]¶: Alias for field number 2

owned_cell_mask: npt.NDArray[np.bool_]¶: Alias for field number 3

owned_dof_mask: npt.NDArray[np.bool_]¶: Alias for field number 4

n_global_cells: int¶: Alias for field number 5

n_global_dofs: int¶: Alias for field number 6

class pantr.bspline.MultiLevelExtraction(space, target='bezier')[source]¶

Bases: object

Per-element multi-level (Bézier) extraction for a THBSplineSpace.

Mirrors SpanwiseElementExtraction: it is constructed from a space and a target reference basis, caches the single-level per-level extractions, and exposes per-element operators via operator(). Because hierarchical refinement introduces a non-constant number of active functions per cell (and the hierarchical basis is not of tensor-product structure), the operators are ragged across cells; there is consequently no constant-shape tabulate / ops_1d.

For a cell with K = active_basis(cid).size active functions, degree p, and dimension d (so n = (p + 1) ** d single-level functions on the cell):

multilevel_operator() returns $M^\epsilon$ of shape (K, n) mapping the level-L tensor-product B-splines on the cell to the active hierarchical functions (independent of target).
operator() returns $C^\epsilon = M^\epsilon E^\epsilon$ of shape (K, n) mapping the target reference basis (Bernstein on $[0, 1]^d$ for "bezier") to the active hierarchical functions.

The operators’ rows are ordered as active_basis() (sorted global dof).

References

Multi-level Bézier extraction for hierarchical local refinement [D'Angella et al., 2018].

Parameters:

space (THBSplineSpace)
target (Target)

_space¶

The hierarchical space being extracted.

Type:: THBSplineSpace

_target¶

The single-level reference basis tag.

Type:: Target

_oslo¶

Cached per-level, per-direction two-scale (Oslo) matrices; _oslo[m][k] maps level m to level m+1 in direction k.

Type:: tuple[tuple[npt.NDArray[np.float64], …], …]

_ext¶

Cache of per-level single-level extractions, built lazily.

Type:: dict[int, SpanwiseElementExtraction]

_coeffs_cache (dict[tuple[int, tuple[int, ...], int], tuple[tuple[int, ...],: npt.NDArray[np.float64]]]): Memoized _element_coeffs results keyed by (origin_level, multi, target_level). A hierarchical function’s coefficients in a finer level’s basis are independent of the cell, but a function’s support covers up to (p + 1) ** d cells per level, so the per-cell operators would otherwise recompute each entry many times. Cached coefficient arrays are frozen read-only.

__init__(space, target='bezier')[source]¶

Create a multi-level extraction for a hierarchical space.

Parameters:

space (THBSplineSpace) – The truncated (or non-truncated) hierarchical space.
target (Target) – Single-level reference basis, one of "bezier", "lagrange", "cardinal". Defaults to "bezier".

Raises:

TypeError – If space is not a THBSplineSpace.
ValueError – If target is not a recognized tag.

Return type:

None

property space: THBSplineSpace¶

Get the underlying hierarchical space.

Returns:: The space supplied at construction time.
Return type:: THBSplineSpace

property target: Literal['bezier', 'lagrange', 'cardinal']¶

Get the single-level reference basis tag.

Returns:: One of "bezier", "lagrange", "cardinal".
Return type:: Target

property dim: int¶

Get the parametric dimension.

Returns:: Number of parametric directions.
Return type:: int

property dtype: type[float64]¶

Get the floating-point dtype of the operators.

Returns:: Always numpy.float64; Oslo matrices and all operators are computed in double precision.
Return type:: type[np.float64]

property num_elements: int¶

Get the number of active cells (elements).

Returns:: space.grid.num_cells.
Return type:: int

active_basis(cid)[source]¶

Return the global dofs labelling the rows of the operators on cell cid.

Parameters:

cid (int) – Active cell flat id in [0, num_elements).

Returns:

Sorted global hierarchical-dof indices (the operator rows), as returned by THBSplineSpace.active_basis().

Return type:

npt.NDArray[np.int64]

Raises:

IndexError – If cid is out of range.
RuntimeError – If the grid has been modified since construction.

multilevel_operator(cid, *, out=None)[source]¶

Return the multi-level extraction operator $M^\epsilon$ on cell cid.

$M^\epsilon$ (shape (K, n)) maps the level-L tensor-product B-splines with support on the cell to the active hierarchical functions (H^\epsilon = M^\epsilon N^{\epsilon,L}). Rows follow active_basis(); columns are the (p + 1) ** d single-level functions on the cell in C-order.

Parameters:

cid (int) – Active cell flat id in [0, num_elements).
out (npt.NDArray[np.float64] | None) – Optional output array of shape (K, n) where K = active_basis(cid).size and n = (p + 1) ** d. Allocated when None.

Returns:

The operator $M^\epsilon$.

Return type:

npt.NDArray[np.float64]

Raises:

IndexError – If cid is out of range.
ValueError – If out has the wrong shape, dtype, or is not writeable.
RuntimeError – If the grid has been modified since construction.

operator(cid, *, out=None)[source]¶

Return the multi-level Bézier extraction $C^\epsilon$ on cell cid.

$C^\epsilon = M^\epsilon E^\epsilon$ (shape (K, n)) maps the target reference basis on the cell to the active hierarchical functions (H^\epsilon = C^\epsilon B). For target="bezier", B is the Bernstein basis on $[0, 1]^d$. Rows follow active_basis().

Parameters:

cid (int) – Active cell flat id in [0, num_elements).
out (npt.NDArray[np.float64] | None) – Optional output array of shape (K, n) where K = active_basis(cid).size and n = (p + 1) ** d. Allocated when None.

Returns:

The operator $C^\epsilon$.

Return type:

npt.NDArray[np.float64]

Raises:

IndexError – If cid is out of range.
ValueError – If out has the wrong shape, dtype, or is not writeable.
RuntimeError – If the grid has been modified since construction.

class pantr.bspline.SpanwiseElementExtraction(space, target, *, lagrange_variant=LagrangeVariant.EQUISPACES)[source]¶

Bases: object

Tensor-product change-of-basis operator across B-spline elements.

For a BsplineSpace of dimension d and a chosen target basis, this class eagerly builds per-direction compact operator storage: only the non-identity rows of each direction’s extraction operator array are retained, reducing memory for identity-heavy spaces (e.g. cardinal spaces on uniform meshes). Per-element d-dimensional operators are never materialized unless explicitly requested via operator() or tabulate(): instead the apply-style methods dispatch to the matrix-free Kronecker kernels in pantr.bspline._extraction_kernels.

With the current 1D builders all per-direction operators are square of size (degree_k + 1, degree_k + 1). The class also supports non-square per-direction operators, so new 1D builders can plug in without changes.

The per-direction data is exposed as two complementary representations:

Compact (compact_ops_1d, idx_maps_1d, is_identity_mask_1d): primary storage, suitable for downstream @njit code that calls the Layer-3 batch kernels directly.
Dense (ops_1d): the full (n_elements_k, n_out_k, n_in_k) layout, reconstructed lazily from compact storage on first access.

Parameters:

space (BsplineSpace)
target (Target)
lagrange_variant (LagrangeVariant)

_space¶

Underlying multi-dimensional B-spline space.

Type:: BsplineSpace

_target¶

Target basis tag.

Type:: Target

_lagrange_variant¶

Point distribution used when target == "lagrange"; ignored otherwise.

Type:: LagrangeVariant

_compact_ops_1d¶

Per-direction compact 3D operator arrays of shape (n_compact_k, n_out_k, n_in_k); only non-identity rows are stored. Always has at least one row to ensure safe Numba indexing.

Type:: tuple[npt.NDArray[np.float32 | np.float64], …]

_idx_maps_1d¶

Per-direction compact index maps of shape (n_elements_k,); _idx_maps_1d[k][e] is the row index into _compact_ops_1d[k] for element e (undefined for identity elements, stored as 0).

Type:: tuple[npt.NDArray[np.intp], …]

_is_identity_mask_1d¶

Per-direction identity masks of shape (n_elements_k,).

Type:: tuple[npt.NDArray[bool], …]

__init__(space, target, *, lagrange_variant=LagrangeVariant.EQUISPACES)[source]¶

Build the per-direction operators and identity masks.

Parameters:

space (BsplineSpace) – Multi-dimensional B-spline space.
target (Target) – One of "bezier", "lagrange", "cardinal".
lagrange_variant (LagrangeVariant) – Point distribution used when target == "lagrange". Defaults to pantr.basis.LagrangeVariant.EQUISPACES.

Raises:

ValueError – If target is not a recognized tag.
NotImplementedError – If any direction of space is periodic; periodic support is deferred to a later version.

Return type:

None

property space: BsplineSpace¶

Get the underlying B-spline space.

Returns:: The space supplied at construction time.
Return type:: BsplineSpace

property target: Literal['bezier', 'lagrange', 'cardinal']¶

Get the target basis tag.

Returns:: One of "bezier", "lagrange", "cardinal".
Return type:: Target

property lagrange_variant: LagrangeVariant¶

Get the Lagrange point distribution used for "lagrange" target.

Returns:: The point distribution. Meaningless for other targets.
Return type:: LagrangeVariant

property dim: int¶

Get the number of tensor-product directions.

Returns:: The dimension d of the space.
Return type:: int

property dtype: DTypeLike¶

Get the floating-point dtype shared by all operators.

Returns:: The dtype inherited from the space (float32 or float64).
Return type:: npt.DTypeLike

property num_intervals: tuple[int, ...]¶

Get the per-direction number of elements (intervals).

Returns:: Length-d tuple (n_elements_0, …, n_elements_{d-1}).
Return type:: tuple[int, …]

property num_total_intervals: int¶

Get the total number of elements across the tensor-product grid.

Returns:: prod(num_intervals).
Return type:: int

property ops_1d: tuple[ndarray[tuple[Any, ...], dtype[float32 | float64]], ...]¶

Get the per-direction 1D operator arrays (dense, reconstructed lazily).

Reconstructs the full (n_elements_k, n_out_k, n_in_k) array from compact storage on first access and caches the result. Identity elements are filled with numpy.eye(n_out, n_in) (rectangular identity for non-square operators); non-identity elements are read from compact_ops_1d.

Returns:: Length-d tuple of read-only 3D arrays; ops_1d[k] has shape (n_elements_k, n_out_k, n_in_k). Intended for consumption by downstream @njit code when the full dense layout is required. For compact-aware downstream code, prefer compact_ops_1d and idx_maps_1d.
Return type:: tuple[npt.NDArray[np.float32 | np.float64], …]

property compact_ops_1d: tuple[ndarray[tuple[Any, ...], dtype[float32 | float64]], ...]¶

Get the per-direction compact operator arrays (non-identity rows only).

Returns:: Length-d tuple of read-only 3D arrays; compact_ops_1d[k] has shape (n_compact_k, n_out_k, n_in_k) where n_compact_k is the number of non-identity elements in direction k (at least 1 to ensure safe Numba indexing). Intended for downstream @njit code alongside idx_maps_1d and is_identity_mask_1d.
Return type:: tuple[npt.NDArray[np.float32 | np.float64], …]

property idx_maps_1d: tuple[ndarray[tuple[Any, ...], dtype[int64]], ...]¶

Get the per-direction compact index maps.

Returns:: Length-d tuple of read-only 1D integer arrays; idx_maps_1d[k] has shape (n_elements_k,) and idx_maps_1d[k][e] is the row index into compact_ops_1d [k] for element e. For identity elements the stored value is 0 (unused; the kernel short-circuits on is_identity_mask_1d). Intended for downstream @njit code.
Return type:: tuple[npt.NDArray[np.intp], …]

property is_identity_mask_1d: tuple[ndarray[tuple[Any, ...], dtype[bool]], ...]¶

Get the per-direction identity masks.

All three targets use structural (multiplicity-based) identity predicates. For "bezier", an element is identity iff both its boundary unique knots have multiplicity >= degree + 1; multiplicities are computed using space.tolerance. For "lagrange", the mask delegates to the Bézier mask when the Lagrange-to-Bernstein matrix equals I (e.g. degree == 1 with equispaced or GLL nodes), returns all-True for degree == 0, and all-False otherwise. For "cardinal", the mask is the structural output of BsplineSpace1D.get_cardinal_intervals().

Returns:: Length-d tuple of read-only 1D boolean arrays; is_identity_mask_1d[k][i] is True iff the i-th element in direction k has an identity operator.
Return type:: tuple[npt.NDArray[bool], …]

property input_shape_per_dir: tuple[int, ...]¶

Get the per-direction input sizes of each element’s operator.

Returns:: (n_in_0, …, n_in_{d-1}).
Return type:: tuple[int, …]

property output_shape_per_dir: tuple[int, ...]¶

Get the per-direction output sizes of each element’s operator.

Returns:: (n_out_0, …, n_out_{d-1}).
Return type:: tuple[int, …]

is_identity_at(cell_idx)[source]¶

Check whether the per-element operator is identity along every direction.

Parameters:: cell_idx (CellIndex) – Element index (flat or per-direction).
Returns:: True iff the d-dimensional operator at cell_idx is the identity (all per-direction operators are identity).
Return type:: bool

property num_identity_elements: int¶

Count elements whose per-direction operators are all identity.

Returns:: The number of fully-identity elements on the tensor-product grid.
Return type:: int

property is_identity: bool¶

Check whether every element on the grid has an identity operator.

Returns:: True iff all per-direction identity masks are all-True, meaning every element’s operator is the identity.
Return type:: bool

per_direction_identity_flags(cell_idx)[source]¶

Return the per-direction identity flags for a single element.

Parameters:: cell_idx (CellIndex) – Element index (flat or per-direction).
Returns:: Length-d tuple of identity flags for the element.
Return type:: tuple[bool, …]

apply(v, cell_idx, *, out=None, scratch=None)[source]¶

Compute out = M @ v for the element at cell_idx.

Here M = kron(M_0, …, M_{d-1}) with M_k the 1D operator at cell_idx in direction k; identity directions short-circuit.

Parameters:

v (npt.NDArray[np.float32 | np.float64]) – Input vector of shape (prod(input_shape_per_dir),).
cell_idx (CellIndex) – Element index (flat or per-direction).
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array of shape (prod(output_shape_per_dir),). Allocated if None.
scratch (npt.NDArray[np.float32 | np.float64] | None) – Optional scratch buffer. Allocated if None.

Returns:

The result array (the same array as out when out was provided).

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

NotImplementedError – If the space has more than 3 directions; specialized kernels only exist for d in {1, 2, 3}.

apply_transpose(v, cell_idx, *, out=None, scratch=None)[source]¶

Compute out = M^T @ v for the element at cell_idx.

Parameters:

v (npt.NDArray[np.float32 | np.float64]) – Input vector of shape (prod(output_shape_per_dir),).
cell_idx (CellIndex) – Element index.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array of shape (prod(input_shape_per_dir),).
scratch (npt.NDArray[np.float32 | np.float64] | None) – Optional scratch buffer.

Returns:

The result array.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

NotImplementedError – If the space has more than 3 directions; specialized kernels only exist for d in {1, 2, 3}.

apply_MT_K_M(K, cell_idx, *, out=None, scratch=None)[source]¶

Compute out = M^T @ K @ M for the element at cell_idx.

Parameters:

K (npt.NDArray[np.float32 | np.float64]) – Input matrix of shape (N_out, N_out) with N_out = prod(output_shape_per_dir).
cell_idx (CellIndex) – Element index.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output matrix of shape (N_in, N_in) with N_in = prod(input_shape_per_dir). Must not alias K.
scratch (npt.NDArray[np.float32 | np.float64] | None) – Optional scratch buffer.

Returns:

The result matrix.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

NotImplementedError – If the space has more than 3 directions; specialized kernels only exist for d in {1, 2, 3}.

apply_M_K_MT(K, cell_idx, *, out=None, scratch=None)[source]¶

Compute out = M @ K @ M^T for the element at cell_idx.

Parameters:

K (npt.NDArray[np.float32 | np.float64]) – Input matrix of shape (N_in, N_in) with N_in = prod(input_shape_per_dir).
cell_idx (CellIndex) – Element index.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output matrix of shape (N_out, N_out) with N_out = prod(output_shape_per_dir). Must not alias K.
scratch (npt.NDArray[np.float32 | np.float64] | None) – Optional scratch buffer.

Returns:

The result matrix.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

NotImplementedError – If the space has more than 3 directions; specialized kernels only exist for d in {1, 2, 3}.

operator(cell_idx, *, out=None)[source]¶

Materialize the full (N_out, N_in) operator for one element.

Assembles the full Kronecker product in memory using numpy.kron(). Prefer the matrix-free apply methods in production code when the full matrix is not needed explicitly.

Parameters:

cell_idx (CellIndex) – Element index.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output matrix of shape (N_out, N_in).

Returns:

The full Kronecker operator.

Return type:

npt.NDArray[np.float32 | np.float64]

tabulate(*, out=None)[source]¶

Materialize per-element operators for every element on the grid.

Parameters:: out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array of shape (num_total_intervals, N_out, N_in). Cells are ordered row-major over num_intervals (so flat index f corresponds to multi-index np.unravel_index(f, num_intervals)).
Returns:: Stacked per-element operators.
Return type:: npt.NDArray[np.float32 | np.float64]

apply_many(v, cell_indices, *, out=None, scratch=None)[source]¶

Compute out[c] = M_c @ v[c] for all cells in the batch.

M_c = kron(M_0[c_0], …, M_{d-1}[c_{d-1}]) with M_k[c_k] the 1D operator at element c_k in direction k; identity directions short-circuit per cell.

Parameters:

v (npt.NDArray[np.float32 | np.float64]) – Batch input vectors, shape (n_cells, N_in) with N_in = prod(input_shape_per_dir).
cell_indices (CellIndicesBatch) – Cell indices — flat 1-D array of shape (n_cells,) or per-direction 2-D array of shape (n_cells, d).
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array of shape (n_cells, N_out). Allocated if None.
scratch (npt.NDArray[np.float32 | np.float64] | None) – Optional per-cell scratch array of shape (n_cells, s) with s >= scratch_size_per_cell. Allocated if None.

Returns:

Result array of shape (n_cells, N_out).

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

NotImplementedError – If the space has more than 3 directions.

apply_transpose_many(v, cell_indices, *, out=None, scratch=None)[source]¶

Compute out[c] = M_c^T @ v[c] for all cells in the batch.

Parameters:

v (npt.NDArray[np.float32 | np.float64]) – Batch input vectors, shape (n_cells, N_out) with N_out = prod(output_shape_per_dir).
cell_indices (CellIndicesBatch) – Cell indices — flat or per-direction.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array of shape (n_cells, N_in). Allocated if None.
scratch (npt.NDArray[np.float32 | np.float64] | None) – Optional per-cell scratch array. Allocated if None.

Returns:

Result array of shape (n_cells, N_in).

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

NotImplementedError – If the space has more than 3 directions.

apply_MT_K_M_many(K, cell_indices, *, out=None, scratch=None)[source]¶

Compute out[c] = M_c^T @ K[c] @ M_c for all cells in the batch.

Parameters:

K (npt.NDArray[np.float32 | np.float64]) – Batch input matrices, shape (n_cells, N_out, N_out) with N_out = prod(output_shape_per_dir).
cell_indices (CellIndicesBatch) – Cell indices — flat or per-direction.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array of shape (n_cells, N_in, N_in). Must not alias K. Allocated if None.
scratch (npt.NDArray[np.float32 | np.float64] | None) – Optional per-cell scratch array. Allocated if None.

Returns:

Result array of shape (n_cells, N_in, N_in).

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

NotImplementedError – If the space has more than 3 directions.

apply_M_K_MT_many(K, cell_indices, *, out=None, scratch=None)[source]¶

Compute out[c] = M_c @ K[c] @ M_c^T for all cells in the batch.

Parameters:

K (npt.NDArray[np.float32 | np.float64]) – Batch input matrices, shape (n_cells, N_in, N_in) with N_in = prod(input_shape_per_dir).
cell_indices (CellIndicesBatch) – Cell indices — flat or per-direction.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array of shape (n_cells, N_out, N_out). Must not alias K. Allocated if None.
scratch (npt.NDArray[np.float32 | np.float64] | None) – Optional per-cell scratch array. Allocated if None.

Returns:

Result array of shape (n_cells, N_out, N_out).

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

NotImplementedError – If the space has more than 3 directions.

class pantr.bspline.THBSpline(space, control_points)[source]¶

Bases: object

An evaluable THB spline function f = Σ_i c_i φ_i.

The hierarchical analogue of Bspline: it pairs a THBSplineSpace with one control point per active hierarchical dof. Because the hierarchical basis is not of tensor-product structure, the control points stay a flat (num_total_basis,) (scalar) or (num_total_basis, rank) (vector-valued) array, rather than the (*num_basis, rank) grid of Bspline.

Parameters:

space (THBSplineSpace)
control_points (npt.ArrayLike)

_space¶

The hierarchical space.

Type:: THBSplineSpace

_control_points¶

Control points reshaped to (num_total_basis, rank) (read-only).

Type:: npt.NDArray[np.float64]

_scalar¶

Whether the field is scalar (control_points was 1-D).

Type:: bool

__init__(space, control_points)[source]¶

Create a THB spline function.

Parameters:

space (THBSplineSpace) – The hierarchical space.
control_points (npt.ArrayLike) – Control points of shape (num_total_basis,) (scalar) or (num_total_basis, rank) (vector-valued), ordered by global dof.

Raises:

TypeError – If space is not a THBSplineSpace.
ValueError – If control_points is not 1-D/2-D or its leading dimension does not equal space.num_total_basis.

Return type:

None

property space: THBSplineSpace¶

Get the underlying hierarchical space.

Returns:: The space supplied at construction time.
Return type:: THBSplineSpace

property control_points: ndarray[tuple[Any, ...], dtype[float64]]¶

Get the control points (read-only view).

Returns:: Shape (num_total_basis,) for a scalar field, (num_total_basis, rank) otherwise. The array is read-only; copy it before modifying.
Return type:: npt.NDArray[np.float64]

property dim: int¶

Get the parametric dimension.

Returns:: Number of parametric directions.
Return type:: int

property degree: tuple[int, ...]¶

Get the per-direction polynomial degrees.

Returns:: Degree per parametric direction (mirrors degree).
Return type:: tuple[int, …]

property rank: int¶

Get the output rank (number of value components).

Returns:: 1 for a scalar field; the number of components otherwise.
Return type:: int

property dtype: dtype[float64]¶

Get the floating-point dtype of the control points.

Returns:: Dtype of the stored control points, always numpy.float64.
Return type:: np.dtype[np.float64]

evaluate(pts, out=None)[source]¶

Evaluate the THB spline at pts.

Each point is located in its active leaf cell; the cell’s active-basis values (THBSplineSpace.tabulate_basis()) are combined with the matching control points. Mirrors evaluate().

Parameters:

pts (npt.ArrayLike) – Parametric points of shape (..., dim).
out (npt.NDArray[np.float64] | None) – Optional output array of the result shape ((...) for a scalar field, (..., rank) otherwise). Allocated when None.

Returns:

Values of shape (...) for a scalar field or (..., rank) for a vector-valued field.

Return type:

npt.NDArray[np.float64]

Raises:

ValueError – If pts does not have trailing dimension dim, any point lies outside the grid domain, or out has the wrong shape/dtype or is not writeable.
RuntimeError – If the grid has been modified since construction, or a cell has no active basis functions (inconsistent space).

evaluate_derivatives(pts, orders, out=None)[source]¶

Evaluate a mixed partial derivative of the THB spline at pts.

Computes the single mixed partial $\partial^{orders} f$, where orders[k] is the derivative order in parametric direction k (with respect to the parametric coordinates). Mirrors evaluate_derivatives().

Parameters:

pts (npt.ArrayLike) – Parametric points of shape (..., dim).
orders (int | Sequence[int]) – Per-direction derivative orders. A scalar is broadcast to every direction. Each entry must be >= 0.
out (npt.NDArray[np.float64] | None) – Optional output array of the result shape. Allocated when None.

Returns:

Derivative values of shape (...) for a scalar field or (..., rank) for a vector-valued field.

Return type:

npt.NDArray[np.float64]

Raises:

ValueError – If orders has the wrong length or a negative entry, if pts does not have trailing dimension dim, any point lies outside the grid domain, or out has the wrong shape/dtype or is not writeable.
RuntimeError – If the grid has been modified since construction, or a cell has no active basis functions (inconsistent space).

refine(cell_ids, *, admissible_class=2)[source]¶

Return the same function on a space with the marked cells refined.

Refines the underlying THBSplineSpace (see THBSplineSpace.refine()) and prolongs the control points onto it, so the returned spline represents the same function – hierarchical refinement is exact (the spaces are nested). This does not mutate self.

Parameters:

cell_ids (npt.ArrayLike) – Flat ids of active cells to refine.
admissible_class (int | None) – Admissibility class m >= 2 to maintain (graded refinement), or None for ungraded refinement. Defaults to 2. See THBSplineSpace.refine().

Returns:

A new spline on the refined space, equal to self as a function.

Return type:

THBSpline

Raises:

IndexError – If any id is outside [0, grid.num_cells).
ValueError – If admissible_class is an integer < 2.
RuntimeError – If the grid has been modified since construction.

refine_region(level, lo, hi, *, admissible_class=2)[source]¶

Return the same function on a space with a rectangular region refined.

Refines the underlying THBSplineSpace (see THBSplineSpace.refine_region()) and prolongs the control points onto it, so the returned spline represents the same function. This does not mutate self.

Parameters:

level (int) – Level at which the box lives. Must satisfy 0 <= level <= grid.max_level.
lo (Sequence[int]) – Per-direction start index (inclusive), in level-level coordinates.
hi (Sequence[int]) – Per-direction end index (exclusive), in level-level coordinates.
admissible_class (int | None) – Admissibility class m >= 2 to maintain, or None for ungraded refinement. Defaults to 2. See THBSplineSpace.refine_region().

Returns:

A new spline on the refined space, equal to self as a function.

Return type:

THBSpline

Raises:

ValueError – If admissible_class is an integer < 2, level is out of range, lo/hi have the wrong length, any lo[k] >= hi[k], or [lo, hi) lies outside the level domain.
RuntimeError – If the grid has been modified since construction.

class pantr.bspline.THBSplineSpace(root_space, grid, *, truncate=True, regularity=None)[source]¶

Bases: object

Hierarchical B-spline space on a HierarchicalGrid.

Built from a root BsplineSpace (level 0) and a HierarchicalGrid carrying the active-cell hierarchy. The per-level tensor-product spaces are obtained by uniformly subdividing the root space according to the grid’s per-direction factor. The active hierarchical basis is the Kraft selection [Kraft, 1997, Vuong et al., 2011]: a level-l tensor-product B-spline is active iff its support lies in the level-l subdomain $\Omega_l$ but not entirely in the finer subdomain $\Omega_{l+1}$.

With truncate=True (the default) the truncated hierarchical basis (THB) is built: each active function that straddles a finer-level refinement boundary has its components on active finer functions removed (Giannelli-Jüttler-Speleers truncation [Giannelli et al., 2012]), restoring the partition of unity. Only truncated functions store a coefficient vector (in the finest tensor-product basis their support reaches); untruncated functions remain plain tensor-product B-splines. With truncate=False the non-truncated hierarchical basis (HB) is built.

This space is a snapshot of the grid at construction time. Calling refine() on the underlying grid after construction invalidates this space; subsequent calls to active_basis() or tabulate_basis() will raise RuntimeError. Create a new THBSplineSpace from the updated grid instead.

Note

active_basis() lists functions whose untruncated support covers a cell; under truncation a few of those may evaluate to exactly zero on the cell. tabulate_basis() always returns the correct (possibly zero) values.

References

Adaptive isogeometric algorithms for hierarchical splines [Garau and Vázquez, 2018]. Per-element multi-level Bézier extraction (used for element assembly and visualization) is provided by MultiLevelExtraction, following D'Angella et al. [2018].

Parameters:

root_space (BsplineSpace)
grid (HierarchicalGrid)
truncate (bool)
regularity (int | Sequence[int | None] | None)

_root_space¶

The level-0 tensor-product space.

Type:: BsplineSpace

_grid¶

The active-cell hierarchy (snapshot reference).

Type:: HierarchicalGrid

_truncate¶

Whether the truncated (THB) basis is used; False for the plain hierarchical (HB) basis.

Type:: bool

_regularity¶

Per-direction continuity used when subdividing to build finer levels.

Type:: tuple[int | None, …]

_level_spaces¶

Per-level tensor-product spaces; index l is the root subdivided to level l.

Type:: tuple[BsplineSpace, …]

_support¶

Per-level, per-direction function-to-cell support arrays; _support[level][k] is the (first_basis, first_cell, last_cell) int64 triple (a _Support1D) for direction k at level.

Type:: tuple

_active_funcs¶

Per-level sorted flat (C-order) indices of the active tensor-product functions.

Type:: tuple[npt.NDArray[np.int64], …]

_func_offset¶

Per-level global-dof base; length num_levels + 1 (cumulative active-function counts).

Type:: npt.NDArray[np.int64]

_num_active¶

Total number of active hierarchical functions.

Type:: int

_grid_snapshot¶

(max_level, num_cells, version) captured at construction; used to detect post-construction grid mutations (the grid’s version counter catches mutations the other two cannot distinguish).

Type:: tuple[int, int, int]

_trunc¶

Map from global dof (int) to _TruncCoeffs; only truncated functions appear (empty when truncate=False).

Type:: dict

__init__(root_space, grid, *, truncate=True, regularity=None)[source]¶

Create a hierarchical B-spline space.

Parameters:

root_space (BsplineSpace) – The level-0 tensor-product B-spline space.
grid (HierarchicalGrid) – Hierarchical grid whose root knot-span grid matches root_space.
truncate (bool) – If True (default), build the truncated (THB) basis; if False, build the non-truncated hierarchical (HB) basis.
regularity (int | Sequence[int | None] | None) – Per-direction continuity at the knots inserted when subdividing to finer levels. A scalar is broadcast to every axis; None (default) uses maximal smoothness. Each non-None entry must satisfy -1 <= regularity[k] < degree[k].

Raises:

TypeError – If root_space is not a BsplineSpace or grid is not a HierarchicalGrid.
ValueError – If grid and root_space disagree on dimension or on the root knot-span grid, if regularity has the wrong length, or if any per-direction regularity value is out of range.

Return type:

None

property grid: HierarchicalGrid¶

Get the underlying hierarchical grid.

Returns:: The active-cell hierarchy this space is built on.
Return type:: HierarchicalGrid

property root_space: BsplineSpace¶

Get the level-0 tensor-product space.

Returns:: The root B-spline space.
Return type:: BsplineSpace

property dim: int¶

Get the parametric dimension.

Returns:: Number of parametric directions.
Return type:: int

property degrees: tuple[int, ...]¶

Get the per-direction polynomial degrees.

Returns:: Degree per direction (the same at every level).
Return type:: tuple[int, …]

property num_levels: int¶

Get the number of hierarchy levels at construction time.

Returns:: Number of levels; stable even if the grid is later refined.
Return type:: int

property truncate: bool¶

Get whether the hierarchical basis is truncated.

Returns:: True for the truncated (THB) basis, False for the plain hierarchical (HB) basis.
Return type:: bool

property num_total_basis: int¶

Get the total number of active hierarchical basis functions.

Mirrors num_total_basis (the hierarchical basis is not tensor-product, so there is no per-direction num_basis).

Returns:: Total active-function count across all levels.
Return type:: int

property num_basis_per_level: tuple[int, ...]¶

Get the number of active basis functions at each level.

Returns:: Active-function count per level.
Return type:: tuple[int, …]

property domain: npt.NDArray[np.float32 | np.float64]¶

Get the parametric domain bounds.

Returns:: Shape (dim, 2) [lo, hi] per direction (from the root space).
Return type:: npt.NDArray[np.float32 | np.float64]

property dtype: npt.DTypeLike¶

Get the floating-point dtype of the space.

Returns:: Always numpy.float64 (THB evaluation is float64).
Return type:: npt.DTypeLike

property tolerance: float¶

Get the numerical tolerance.

Returns:: The root space’s tolerance.
Return type:: float

level_space(level)[source]¶

Return the tensor-product space at level.

Returns a construction-time snapshot; not affected by subsequent grid mutations (no stale check is performed).

Parameters:: level (int) – Hierarchy level in [0, num_levels).
Returns:: The root space subdivided to level.
Return type:: BsplineSpace
Raises:: ValueError – If level is out of range.

active_function_indices(level)[source]¶

Return the flat indices of the active functions at level.

Returns a construction-time snapshot; not affected by subsequent grid mutations (no stale check is performed).

Parameters:: level (int) – Hierarchy level in [0, num_levels).
Returns:: Sorted flat (C-order) level-level function indices selected by the Kraft rule. A fresh copy is returned.
Return type:: npt.NDArray[np.int64]
Raises:: ValueError – If level is out of range.

active_basis(cid)[source]¶

Return the global dofs of the active functions whose support intersects cell cid.

Parameters:

cid (int) – Active cell flat id in [0, grid.num_cells).

Returns:

Sorted global hierarchical-dof indices of the functions whose support intersects cell cid.

Return type:

npt.NDArray[np.int64]

Raises:

IndexError – If cid is out of range [0, grid.num_cells).
RuntimeError – If the grid has been modified since construction.

restrict(cell_ids)[source]¶

Return the windowed sub-space over a subset of active cells.

Windows this space to the root-cell-aligned bounding box of cell_ids: the hierarchical grid is restricted (pantr.grid.HierarchicalGrid.restrict()), the root space is windowed (pantr.bspline.BsplineSpace.restrict()), and a new THBSplineSpace is rebuilt on the sub-grid (re-running the Kraft active-function selection and truncation).

Unlike the tensor-product pantr.bspline.BsplineSpace.restrict(), the windowed THB basis equals the global one only over the interior cells – those whose entire (cross-level) function-support-closure lies inside the window – because Kraft selection and truncation depend on the subdomain near the window boundary. Callers make the cells they care about interior by padding cell_ids with a support-closure halo.

Parameters:

cell_ids (npt.ArrayLike) – Active cell flat ids to span; duplicates ignored.

Returns:

The windowed THBSplineSpace and a read-only local_to_global_dof map; entry d is the global hierarchical dof of local dof d when the local function matches a globally-active function of the same level and multi-index, else -1. Values are exact over interior cells; functions near the window boundary may map to -1.

Return type:

Raises:

ValueError – If cell_ids is empty.
TypeError – If cell_ids is not integer-valued.
IndexError – If any cell id is out of range [0, grid.num_cells).

tabulate_basis(cid, pts, out_basis=None, out_dofs=None)[source]¶

Evaluate the active hierarchical functions on cell cid at pts.

Untruncated functions are a single tensor-product B-spline (the product of their 1D B-spline values). Truncated functions are evaluated from their stored coefficients in the finest tensor-product basis their support reaches. The returned columns are ordered as dofs (the sorted global dofs, equal to active_basis()); a listed truncated function may evaluate to exactly zero on the cell. Mirrors the (basis, first_basis) two-return of tabulate_basis().

Parameters:

cid (int) – Active cell flat id in [0, grid.num_cells).
pts (npt.ArrayLike) – Parametric points of shape (..., dim) lying in cell cid. Points outside the cell’s bounds raise ValueError.
out_basis (npt.NDArray[np.float64] | None) – Optional output array of shape (..., K) with K = active_basis(cid).size. Allocated when None.
out_dofs (npt.NDArray[np.int64] | None) – Optional output array of shape (K,) for the dofs. Allocated when None.

Returns:

(values, dofs) of shapes (..., K) and (K,).

Return type:

tuple[npt.NDArray[np.float64], npt.NDArray[np.int64]]

Raises:

IndexError – If cid is out of range [0, grid.num_cells).
ValueError – If pts does not have trailing dimension dim, if any point lies outside the bounds of cell cid, or if out_basis / out_dofs has the wrong shape, dtype, or is not writeable.
RuntimeError – If the grid has been modified since construction.

tabulate_basis_derivatives(cid, pts, orders, out_basis=None, out_dofs=None)[source]¶

Evaluate a mixed partial derivative of the active functions on cell cid.

Computes the single mixed partial $\partial^{orders}$ of each active hierarchical function, where orders[k] is the derivative order in parametric direction k (derivatives are with respect to the parametric coordinates). Untruncated functions differentiate as a tensor product of 1D B-spline derivatives; truncated functions apply their stored coefficients to the B-spline derivatives at their representation level. The returned columns are ordered as dofs (the sorted global dofs, equal to active_basis()).

Parameters:

cid (int) – Active cell flat id in [0, grid.num_cells).
pts (npt.ArrayLike) – Parametric points of shape (..., dim) lying in cell cid. Points outside the cell’s bounds raise ValueError.
orders (int | Sequence[int]) – Per-direction derivative orders. A scalar is broadcast to every direction. Each entry must be >= 0; orders exceeding the degree yield zero.
out_basis (npt.NDArray[np.float64] | None) – Optional output array of shape (..., K) with K = active_basis(cid).size. Allocated when None.
out_dofs (npt.NDArray[np.int64] | None) – Optional output array of shape (K,) for the dofs. Allocated when None.

Returns:

(values, dofs) of shapes (..., K) and (K,).

Return type:

tuple[npt.NDArray[np.float64], npt.NDArray[np.int64]]

Raises:

IndexError – If cid is out of range [0, grid.num_cells).
ValueError – If orders has the wrong length or a negative entry, if pts does not have trailing dimension dim, if any point lies outside cell cid, or if out_basis / out_dofs has the wrong shape, dtype, or is not writeable.
RuntimeError – If the grid has been modified since construction.

refine(cell_ids, *, admissible_class=2)[source]¶

Return a new space with the marked cells refined.

This method does not mutate self or its grid: a fresh grid is refined and a new THBSplineSpace is built; self and its grid are unchanged.

With admissible_class=m (the default m=2) the refinement is graded so the resulting mesh is admissible of class m (the truncated functions acting on any cell span at most m successive levels), following the recursive refinement-neighborhood algorithm of Carraturo et al. (2019). This assumes the current mesh is already admissible of class m (true for the root and for any mesh built via graded refine()). With admissible_class=None exactly the marked cells are refined (no grading).

Parameters:

cell_ids (npt.ArrayLike) – Flat ids of active cells to refine.
admissible_class (int | None) – Admissibility class m >= 2 to maintain, or None for ungraded refinement. Defaults to 2.

Returns:

A new space on the refined grid (same root_space, truncate, and regularity).

Return type:

Raises:

IndexError – If any id is outside [0, grid.num_cells).
ValueError – If admissible_class is an integer < 2.
RuntimeError – If the grid has been modified since construction.

refine_region(level, lo, hi, *, admissible_class=2)[source]¶

Return a new space with the active cells in a rectangular region refined.

The region is the integer cell-index box [lo, hi) at level (in level-level coordinates), matching the convention of pantr.grid.HierarchicalGrid.refine(). Only the currently-active leaf cells inside the box are refined; the rest of the box (already refined, or not present at level) is ignored. If the box contains no active leaf cells, the call is a no-op and returns a space equivalent to self. This is the region-based counterpart of refine(), which marks individual cells by flat id.

Like refine(), this does not mutate self or its grid: a fresh grid is refined and a new THBSplineSpace is returned. Calls chain, so successive regions refine progressively (graded by default).

Parameters:

level (int) – Level at which the box lives. Must satisfy 0 <= level <= grid.max_level.
lo (Sequence[int]) – Per-direction start index (inclusive), in level-level coordinates.
hi (Sequence[int]) – Per-direction end index (exclusive), in level-level coordinates.
admissible_class (int | None) – Admissibility class m >= 2 to maintain (graded refinement), or None for ungraded refinement. Defaults to 2. See refine().

Returns:

A new space on the refined grid (same root_space, truncate, and regularity).

Return type:

Raises:

ValueError – If admissible_class is an integer < 2, level is out of range, lo/hi have the wrong length, any lo[k] >= hi[k], or any part of [lo, hi) lies outside the level domain.
RuntimeError – If the grid has been modified since construction.

coarsen(cell_ids, *, admissible_class=2)[source]¶

Return a new space with the marked cells coarsened away.

A parent cell is reactivated (its children removed) only when all of its children are marked active leaves, mirroring the coarsening algorithm of Carraturo et al. (2019, Alg. 5). With admissible_class=None this is the exact inverse of refine(): space.refine(cells).coarsen(children_of(cells)) recovers space. With admissible_class=m the guard may suppress some coarsenings, so the recovery holds only when the guard permits them all.

With admissible_class=m (the default m=2) a parent is reactivated only if its coarsening neighborhood (Def. 3.5) is empty, so the resulting mesh stays admissible of class m. With admissible_class=None that guard is skipped.

The space is immutable: a fresh grid is coarsened and a new THBSplineSpace is built; self and its grid are unchanged. An empty cell_ids is valid and returns an unchanged copy of the space.

Parameters:

cell_ids (npt.ArrayLike) – Flat ids of active leaf cells to coarsen away. An empty array is valid and produces an unchanged copy.
admissible_class (int | None) – Admissibility class m >= 2 to maintain, or None to skip the admissibility guard. Defaults to 2.

Returns:

A new space on the coarsened grid (same root_space, truncate, and regularity).

Return type:

Raises:

IndexError – If any id is outside [0, grid.num_cells).
ValueError – If admissible_class is an integer < 2.
RuntimeError – If the grid has been modified since construction.

prolongation_to(fine)[source]¶

Return the prolongation matrix from this space to a refinement fine.

The hierarchical spaces are nested (V_h ⊆ V_h'), so every function of this (coarse) space lies in fine. The returned matrix P maps a coefficient vector in this space’s basis to the coefficients of the same function in fine’s basis: if u are coarse coefficients, P @ u are the fine coefficients.

It is built column-by-column following the local two-scale construction used in practice (Garau & Vazquez 2018; D’Angella et al. 2018): each coarse function is matched against only the fine functions over its support, expressed in the deepest level present there (not the global finest level), and reproduced by a small local least-squares solve. Functions far from the refinement yield trivial (identity) columns, so cost and sparsity follow the refined region.

Parameters:

fine (THBSplineSpace) – A refinement of this space (same root_space, factor, regularity, and truncate; more levels / refined cells).

Returns:

Matrix P of shape (fine.num_total_basis, self.num_total_basis).

Return type:

npt.NDArray[np.float64]

Raises:

TypeError – If fine is not a THBSplineSpace.
ValueError – If fine is not a refinement of this space (mismatched root/factor/regularity/truncation, fewer levels, or the prolongation residual is non-negligible).

restriction_to(coarse)[source]¶

Return the restriction matrix from this space to a coarsening coarse.

self must be a refinement of coarse. The restriction is the algebraic pseudo-inverse of the prolongation, R = pinv(P) with P = coarse.prolongation_to(self). It is assembly-free (no mass matrix) and satisfies R @ P == I, so restricting a prolonged coarse field recovers it exactly; for a general fine field R @ u_fine is the least-squares (coefficient-space) projection onto the coarse space.

Parameters:

coarse (THBSplineSpace) – A coarsening of this space (self is a refinement of coarse).

Returns:

Matrix R of shape (coarse.num_total_basis, self.num_total_basis).

Return type:

npt.NDArray[np.float64]

Raises:

TypeError – If coarse is not a THBSplineSpace.
ValueError – If self is not a refinement of coarse.

class pantr.bspline.THBSplineSpaceRestriction(space, local_to_global_dof, local_to_global_cell)[source]¶

Bases: NamedTuple

Result of THBSplineSpace.restrict(): a windowed THB space with its maps.

space – the windowed THBSplineSpace rebuilt on the restricted grid; its basis equals the global basis pointwise over interior cells (those whose function-support-closure lies inside the window).
local_to_global_dof – read-only (space.num_total_basis,) map; entry d is the global hierarchical dof of local dof d when the local function matches a globally-active function of the same level and multi-index, else -1 (local function active in the sub-space but absent from the global active set – arises near the window boundary where Kraft selection may differ). Values are exact over interior cells.
local_to_global_cell – read-only (space.grid.num_cells,) map; entry c is the global flat cell id of local cell c. Same ordering as space.grid.

Parameters:

space (THBSplineSpace)
local_to_global_dof (npt.NDArray[np.int64])
local_to_global_cell (npt.NDArray[np.int64])

space: THBSplineSpace¶: Alias for field number 0

local_to_global_dof: npt.NDArray[np.int64]¶: Alias for field number 1

local_to_global_cell: npt.NDArray[np.int64]¶: Alias for field number 2

pantr.bspline.build_local(global_space, partition, rank)[source]¶

Build the rank-local windowed space of a distributed B-spline or THB space.

Windows global_space to the cells owned by rank together with their support-closure halo, so the local basis equals the global basis pointwise over the rank’s owned cells. Dispatches on the space type: a BsplineSpace uses the tensor-product path (compute_halo() / dof_owner()); a THBSplineSpace uses the hierarchical path (cross-level support closure and ownership).

Parameters:

global_space (BsplineSpace | THBSplineSpace) – The global space (non-periodic).
partition (Partition) – Owner of every cell of the space’s grid; cell_owner length must equal the global cell count.
rank (int) – The rank whose local space is built; must be in [0, partition.n_parts).

Returns:

The windowed space plus local-to-global cell/DOF maps and ownership masks.

Return type:

LocalSpace

Raises:

ValueError – If global_space is a periodic BsplineSpace (THB spaces have no periodicity); or if rank is out of range, partition does not match the space’s cell count, or rank owns no cells.

pantr.bspline.compute_halo(space, owned_cells)[source]¶

Return the support-closure halo of owned_cells.

The halo is the set of knot-span cells, excluding the owned cells, covered by the support of any B-spline function non-zero on an owned cell. A rank owning owned_cells needs exactly these extra cells so every function touching its owned cells is fully represented over the union of the owned and halo cells. For open/uniform knots this is the degree-wide halo; general or repeated knots are handled exactly via the per-direction support.

Parameters:

space (BsplineSpace) – Tensor-product B-spline space (non-periodic). Its knot-span grid has num_total_intervals cells.
owned_cells (npt.ArrayLike) – Flat cell ids (C-order over num_intervals) owned by the rank. Duplicates are ignored.

Returns:

Sorted, read-only flat ids of the halo cells – those in the support closure but not in owned_cells.

Return type:

npt.NDArray[np.int64]

Raises:

ValueError – If any axis is periodic.
TypeError – If owned_cells is not integer-valued.
IndexError – If any owned cell id is out of range [0, num_total_intervals).

pantr.bspline.coupling_graph(space, *, cell_weights=None)[source]¶

Build the cell-coupling graph of a B-spline or THB-spline space.

Two cells are joined by an edge when they share at least one basis function; the edge weight is the number of shared functions, and each vertex (cell) carries an optional assembly-cost weight. The result is the dual graph a graph partitioner uses to minimize cross-rank DOF coupling.

Parameters:

space (BsplineSpace | THBSplineSpace) – The space whose cells to couple. A BsplineSpace must be non-periodic.
cell_weights (npt.ArrayLike | None) – Optional per-cell assembly cost, shape (num_cells,), finite and non-negative. None means uniform.

Returns:

The coupling graph in METIS / Scotch CSR adjacency format.

Return type:

CouplingGraph

Raises:

TypeError – If space is neither a BsplineSpace nor a THBSplineSpace.
ValueError – If space is a periodic BsplineSpace, or if cell_weights has the wrong shape or invalid values.

Example

>>> from pantr.bspline import coupling_graph, create_uniform_space
>>> space = create_uniform_space(2, 4)
>>> graph = coupling_graph(space)
>>> graph.num_vertices
4
>>> int(graph.xadj[-1]) == graph.adjncy.size
True

pantr.bspline.create_cardinal_knots(num_intervals, degree, dtype=<class 'numpy.float64'>)[source]¶

Create a knot vector for cardinal B-spline basis functions.

Cardinal B-splines are B-splines defined on uniform knot vectors with maximum continuity, where the basis functions in the central region have the same shape and are translated versions of each other.

Parameters:

num_intervals (int) – Number of intervals in the domain. Must be at least 1.
degree (int) – B-spline degree. Must be non-negative.
dtype (npt.DTypeLike) – Data type for the knot vector. It must be either float32 or float64. Defaults to np.float64.

Returns:

Cardinal B-spline knot vector: with uniform spacing.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If num_intervals < 1, degree < 0, or dtype is not float32/float64.

Example

>>> create_cardinal_knots(2, 2)
array([-2., -1.,  0.,  1.,  2.,  3., 4.])

pantr.bspline.create_from_bezier(bezier, *, copy=True)[source]¶

Create a B-spline from a Bézier.

Builds a Bspline with Bézier-like knot vectors ([0]*(p+1) + [1]*(p+1) per direction) whose control points are taken from the given Bézier.

Parameters:

bezier (Bezier) – The source Bézier.
copy (bool) – If True (default), the control points are deep-copied into the new B-spline. If False, the B-spline shares the same underlying control point array.

Returns:

Equivalent B-spline with Bézier-like knots.

Return type:

pantr.bspline.create_greville_lattice(space)[source]¶

Compute the tensor-product Greville abscissae as a PointsLattice.

Returns a PointsLattice whose per-direction arrays are the Greville abscissae of each 1D sub-space.

Parameters:: space (BsplineSpace) – The multi-dimensional B-spline space.
Returns:: Tensor-product grid of Greville abscissae.
Return type:: PointsLattice

Example

>>> from pantr.bspline import BsplineSpace1D, BsplineSpace, create_uniform_open_knots
>>> knots = create_uniform_open_knots(2, 2)
>>> s1d = BsplineSpace1D(knots, 2)
>>> space = BsplineSpace([s1d, s1d])
>>> lattice = create_greville_lattice(space)
>>> lattice.pts_per_dir[0]
array([0. , 0.25, 0.75, 1.  ])

pantr.bspline.create_thb_space(root, factor=2, *, truncate=True, regularity=None)[source]¶

Create a trivial (unrefined) THB-spline space from a B-spline space.

Convenience factory that wraps root in a single-level HierarchicalGrid (its knot-span grid, ready to subdivide by factor) and builds the corresponding THBSplineSpace. The result has one level and its active basis coincides with that of root; refine it with THBSplineSpace.refine() or THBSplineSpace.refine_region().

It is the ergonomic entry point for lifting an existing BsplineSpace into a single-level hierarchy, leaving the two-argument THBSplineSpace constructor for callers that build the hierarchical grid explicitly.

Parameters:

root (BsplineSpace) – The level-0 tensor-product B-spline space.
factor (int | Sequence[int]) – Per-direction subdivision factor used when the space is later refined. A scalar is broadcast to every axis. Each entry must be >= 1. Defaults to 2 (dyadic refinement).
truncate (bool) – If True (default), build the truncated (THB) basis; if False, build the non-truncated hierarchical (HB) basis.
regularity (int | Sequence[int | None] | None) – Per-direction continuity at the knots inserted when subdividing to finer levels. See THBSplineSpace. Defaults to None (maximal smoothness).

Returns:

A single-level THB space over root.

Return type:

Raises:

ValueError – If any factor entry is < 1, factor has the wrong length, or regularity is out of range for root’s degrees.

Example

>>> from pantr.bspline import create_uniform_space, create_thb_space
>>> thb = create_thb_space(create_uniform_space([2, 2], [8, 8]))
>>> thb.num_levels
1
>>> thb = thb.refine_region(0, [0, 0], [4, 4])  # refine the lower-left quarter
>>> thb.num_levels
2

pantr.bspline.create_uniform_open_knots(num_intervals, degree, continuity=None, domain=None, dtype=None)[source]¶

Create a uniform open knot vector.

An open knot vector has the first and last knots repeated (degree+1) times, ensuring the B-spline interpolates the first and last control points.

Parameters:

num_intervals (int) – Number of intervals in the domain. Must be non-negative.
degree (int) – B-spline degree. Must be non-negative.
continuity (Optional[int]) – Continuity level at interior knots. Must be between -1 and degree-1. Defaults to degree-1 (maximum continuity).
domain (Optional[tuple[np.float32 | np.float64 | float, np.float32 | np.float64 | float]]) – Domain boundaries as (start, end). Defaults to (0.0, 1.0) if not provided.
dtype (Optional[np.dtype]) – Data type for the knot vector. If None, inferred from start/end or defaults to float64.

Returns:

Open knot vector with uniform spacing.

Return type:

npt.NDArray[np.floating]

Raises:

ValueError – If any parameter is invalid.

Example

>>> create_uniform_open_knots(2, 2, domain=(0.0, 1.0))
array([0., 0., 0., 0.5, 1., 1., 1.])

pantr.bspline.create_uniform_periodic_knots(num_intervals, degree, continuity=None, domain=None, dtype=<class 'numpy.float64'>)[source]¶

Create a uniform periodic knot vector.

A periodic knot vector extends beyond the domain boundaries to ensure periodicity of the B-spline basis functions.

Parameters:

num_intervals (int) – Number of intervals in the domain. Must be non-negative.
degree (int) – B-spline degree. Must be non-negative.
continuity (Optional[int]) – Continuity level at interior knots. Must be between -1 and degree-1. Defaults to degree-1 (maximum continuity).
domain (Optional[tuple[np.float32 | np.float64 | float, np.float32 | np.float64 | float]]) – Domain boundaries as (start, end). Defaults to (0.0, 1.0) if not provided.
dtype (Optional[np.dtype]) – Data type for the knot vector. If None, inferred from the domain endpoints. Defaults to np.float64.

Returns:

Periodic knot vector with uniform spacing.

Return type:

npt.NDArray[np.floating]

Raises:

ValueError – If any parameter is invalid.

Example

>>> create_uniform_periodic_knots(2, 2, domain=(0.0, 1.0))
array([-1. , -0.5,  0. ,  0.5,  1. ,  1.5,  2. ])

pantr.bspline.create_uniform_space(degree, num_intervals, *, continuity=None, periodic=False, domain=None, dtype=<class 'numpy.float64'>)[source]¶

Create a tensor-product B-spline space with uniform knot vectors.

Scalar arguments are broadcast to all parametric directions. The parametric dimension is inferred from whichever argument is given as a sequence (they must all agree in length when more than one is a sequence).

Uses create_uniform_open_knots() for non-periodic directions and create_uniform_periodic_knots() for periodic ones.

Parameters:

degree (int | Sequence[int]) – Polynomial degree per direction.
num_intervals (int | Sequence[int]) – Number of elements per direction.
continuity (int | Sequence[int] | None) – Interior knot continuity per direction. Defaults to degree - 1 (maximum continuity).
periodic (bool | Sequence[bool]) – Whether each direction is periodic. Defaults to False.
domain (tuple[float32 | float64 | float, float32 | float64 | float] | Sequence[tuple[float32 | float64 | float, float32 | float64 | float]] | None) – Domain boundaries per direction as (start, end) tuples. A single tuple is broadcast. Defaults to (0.0, 1.0).
dtype (npt.DTypeLike) – Data type for the knot vectors. Defaults to np.float64.

Returns:

A tensor-product B-spline space.

Return type:

BsplineSpace

Raises:

ValueError – If sequence lengths are inconsistent.

Example

>>> space = create_uniform_space(3, 4, periodic=True, domain=(0.0, 2.0))
>>> space.dim
1
>>> space.degrees
(3,)

pantr.bspline.dof_owner(space, partition)[source]¶

Return the owner rank of every global DOF (lex-first-active-cell rule).

Each global B-spline DOF is owned by the rank that owns the active cell with the smallest flat id in the DOF’s support. A DOF whose support contains no active cell (cell_owner == -1 throughout) is a dead DOF, assigned -1.

Parameters:

space (BsplineSpace) – Tensor-product B-spline space (non-periodic). DOFs are flat-indexed in C-order over num_basis.
partition (Partition) – Owner of every knot-span cell; cell_owner must have length space.num_total_intervals.

Returns:

Read-only (num_total_basis,) owner rank per DOF; -1 for dead DOFs.

Return type:

npt.NDArray[np.int32]

Raises:

ValueError – If any axis is periodic.
ValueError – If partition.cell_owner length does not match space.num_total_intervals.

pantr.bspline.fit_bspline(values, nodes, space, *, tol=None)[source]¶

Construct a B-spline from pre-evaluated sample values at known nodes.

Recover B-spline coefficients by solving the collocation system. For tensor-product nodes, uses per-direction solves (Kronecker structure). For scattered nodes, builds the full collocation matrix and solves via SVD.

The output dtype is inferred from values.

Supports two point layouts:

Tensor-product (a PointsLattice, a single 1D array, or a sequence of 1D arrays): values must have shape (*n_pts_per_dir) (scalar) or (*n_pts_per_dir, rank) (vector).
Scattered (a 2D ndarray of shape (n_pts, dim)): values must have shape (n_pts,) (scalar) or (n_pts, rank) (vector).

Parameters:

values (npt.ArrayLike) – Sample values at the nodes.
nodes (PointsLattice | npt.NDArray[np.floating[Any]] | Sequence[npt.NDArray[np.floating[Any]]]) –
Interpolation nodes.
- A PointsLattice: tensor-product grid.
- A 1D ndarray: 1D tensor-product (single direction).
- A sequence of 1D ndarray values: N-D tensor-product.
- A 2D ndarray of shape (n_pts, dim): scattered points.
space (BsplineSpace) – The target B-spline space.
tol (float | None) – SVD truncation tolerance. If None, defaults to 100 * machine_epsilon. Only affects overdetermined or near-singular systems; square well-conditioned systems use a direct solve.

Returns:

A non-rational B-spline.

Return type:

Raises:

TypeError – If space is not a BsplineSpace.
ValueError – If nodes are inconsistent with space, or the system is underdetermined for scattered nodes.

Example

>>> import numpy as np
>>> from pantr.bspline import create_uniform_space, fit_bspline
>>> space = create_uniform_space(3, 4)
>>> nodes = np.linspace(0, 1, 20)
>>> vals = np.sin(nodes)
>>> b = fit_bspline(vals, [nodes], space)

pantr.bspline.get_greville_abscissae(space)[source]¶

Compute the Greville abscissae (knot averages) of a 1D B-spline space.

Each Greville abscissa is the average of degree consecutive internal knots: g_i = (1/p) * sum(knots[i+1 : i+p+1]) for i = 0, ..., n-1, where n is the number of basis functions and p is the degree.

For periodic spaces, the Greville points are computed from the full knot vector and then wrapped into the domain [a, b).

Parameters:

space (BsplineSpace1D) – The 1D B-spline space.

Returns:

Array of shape (num_basis,): containing one Greville abscissa per basis function.

Return type:

npt.NDArray[np.float32 | np.float64]

Example

>>> from pantr.bspline import BsplineSpace1D, create_uniform_open_knots
>>> knots = create_uniform_open_knots(4, 3)
>>> space = BsplineSpace1D(knots, 3)
>>> get_greville_abscissae(space)
array([0.  , 0.08333333, 0.25, 0.5 , 0.75, 0.91666667, 1.  ])

pantr.bspline.interpolate_bspline(func, space, *, nodes=None, boundary_derivatives=None, tol=None)[source]¶

Interpolate a callable onto a B-spline space.

Evaluate func on a tensor-product grid of interpolation nodes and recover B-spline coefficients by solving per-direction collocation systems (Kronecker structure).

Parameters:

func (Callable[..., npt.ArrayLike]) – Function to interpolate. Called as func(lattice) where lattice is a PointsLattice representing the tensor-product sampling grid. Must return an array of shape (n_total,) for scalar or (n_total, rank) for vector-valued functions, where n_total = prod(n_pts).
space (BsplineSpace) – The target B-spline space.
nodes (Literal['greville'] | PointsLattice | npt.NDArray[np.floating[Any]] | Sequence[npt.NDArray[np.floating[Any]]] | None) –
Interpolation node selection.
- None or "greville" (default): Greville abscissae (one per basis function per direction).
- A PointsLattice: custom tensor-product grid.
- A 1D ndarray: custom nodes for a 1D space.
- A sequence of 1D ndarray values: per-direction custom nodes.
boundary_derivatives (Sequence[tuple[int, ...] | None] | None) – Per-direction boundary derivative constraints. Each entry is (n_left, n_right) specifying how many derivative orders to constrain at left/right boundaries. The corresponding derivatives are set to zero. None entries skip that direction. Ignored for periodic directions. Defaults to None (no derivative constraints).
tol (float | None) – SVD truncation tolerance for the collocation solve. Singular values below tol * sigma_max are treated as zero. If None, defaults to 100 * machine_epsilon. Only affects overdetermined or near-singular systems; square well-conditioned systems use a direct solve.

Returns:

A non-rational B-spline whose evaluation approximates func.

Return type:

Raises:

TypeError – If space is not a BsplineSpace.
ValueError – If nodes is inconsistent with space, or the callable returns an unexpected shape.

Example

>>> import numpy as np
>>> from pantr.bspline import create_uniform_space, interpolate_bspline
>>> space = create_uniform_space(3, 4)
>>> b = interpolate_bspline(
...     lambda lat: lat.pts_per_dir[0] ** 2, space
... )
>>> b.degree
(3,)

pantr.bspline.l2_project_bspline(func, space, *, n_quad=None, quadrature='gauss-legendre', boundary_interpolation=False, tol=None)[source]¶

L2-project a callable function onto a B-spline space.

Assemble per-direction mass matrices and load vectors using per-element quadrature, then solve the normal equations via sequential 1D solves (Kronecker structure).

The output dtype is inferred from the return value of func.

Parameters:

func (Callable[..., npt.ArrayLike]) – Function to project. Called as func(lattice) where lattice is a PointsLattice of quadrature points. Must return an array of shape (n_total,) for scalar or (n_total, rank) for vector-valued functions.
space (BsplineSpace) – The target B-spline space.
n_quad (int | Sequence[int] | None) – Quadrature points per element per direction. Defaults to degree + 1.
quadrature (Literal["gauss-legendre", "gauss-lobatto"]) – Quadrature rule type. Defaults to "gauss-legendre".
boundary_interpolation (bool | Sequence[tuple[bool, bool]]) –
Replace boundary rows with interpolation conditions.
- False (default): pure L2 projection.
- True: interpolate at all non-periodic boundaries.
- A sequence of (left, right) bool pairs: per-direction boundary flags.
tol (float | None) – SVD truncation tolerance. If None, defaults to 100 * machine_epsilon.

Returns:

A non-rational B-spline whose evaluation is the L2 best approximation of func in the given space.

Return type:

Raises:

TypeError – If space is not a BsplineSpace.
ValueError – If n_quad or boundary_interpolation is inconsistent with space.

Example

>>> import numpy as np
>>> from pantr.bspline import create_uniform_space, l2_project_bspline
>>> space = create_uniform_space(3, 4)
>>> b = l2_project_bspline(
...     lambda lat: lat.pts_per_dir[0] ** 2, space
... )

pantr.bspline.make_struct_view(extraction)[source]¶

Bundle a SpanwiseElementExtraction into a Numba-passable struct view.

Shares the underlying per-direction arrays (no copies). The arrays are already marked read-only by SpanwiseElementExtraction, so the returned view is safe to pass into @njit code without risk of accidental mutation.

Parameters:: extraction (SpanwiseElementExtraction) – Source extraction object.
Returns:: Named tuple wrapping the extraction’s compact storage and shape metadata. Suitable for direct use as a single argument to @njit functions that call the Layer-3 batch kernels in pantr.bspline._extraction_kernels.
Return type:: ExtractionStructView

Example

>>> from pantr.bspline import BsplineSpace1D, BsplineSpace
>>> from pantr.bspline import SpanwiseElementExtraction, make_struct_view
>>> sp = BsplineSpace1D([0, 0, 0, 1, 2, 2, 2], 2)
>>> space = BsplineSpace([sp, sp])
>>> ext = SpanwiseElementExtraction(space, "bezier")
>>> view = make_struct_view(ext)
>>> view.dim
2

pantr.bspline.partition_graph(coupling, n_parts, *, backend='spectral', cell_active=None)[source]¶

Partition a cell-coupling graph into n_parts rank subdomains.

Balances CouplingGraph.vertex_weights across parts while minimizing the weight of cut edges (shared DOFs). Cells excluded by cell_active get owner -1 and are dropped from the graph.

Parameters:

coupling (CouplingGraph) – The cell-coupling graph (see coupling_graph()); its vertex_weights drive the load balance and edge_weights the cut cost.
n_parts (int) – Number of parts (ranks); must be >= 1.
backend (str) – "spectral" (default; recursive Fiedler bisection, no extra dependency, never leaves a rank empty) or "metis" (METIS via the optional pymetis package; higher-quality min-cut, but may leave a rank empty).
cell_active (npt.ArrayLike | None) – Optional boolean mask, shape (coupling.num_vertices,); inactive cells get owner -1 and are excluded. None means all active.

Returns:

A per-cell owner assignment with n_parts parts; -1 for inactive cells, otherwise a rank in range(n_parts).

Return type:

Partition

Raises:

TypeError – If coupling is not a CouplingGraph.
ValueError – If backend is unknown; if n_parts < 1; if cell_active has the wrong shape or marks no cell active; or if n_parts exceeds the number of active cells.
ImportError – If backend="metis" but pymetis is not installed.

pantr.bspline.quasi_interpolate_bspline(func, space, *, kind='llm')[source]¶

Quasi-interpolate a callable onto a tensor-product B-spline space.

Builds the Lee-Lyche-Mørken local projector: a local, point-value functional per basis function whose weights come from inverting a small collocation block. The operator reproduces the spline space (Q B_β = B_β) and is local.

Parameters:

func (Callable) – Function to quasi-interpolate. Called once on an (M, dim) point array and must return (M,) (scalar) or (M, rank) (vector-valued). Unlike the lattice-based interpolate_bspline() / l2_project_bspline(), func receives a flat point array because the QI samples scattered local points.
space (BsplineSpace) – The target tensor-product B-spline space.
kind (QIKind) – Quasi-interpolant kind. Only "llm" (Lee-Lyche-Mørken) is currently supported. Defaults to "llm".

Returns:

A non-rational B-spline whose evaluation quasi-interpolates func.

Return type:

Raises:

TypeError – If space is not a BsplineSpace.
ValueError – If kind is not recognized, or if func returns an output with an invalid shape (0-D, more than 2-D, or wrong leading dimension).

Note

All internal computation uses float64. The returned control points are cast to space.dtype at the end, so a float32 space incurs a precision loss on the final cast.

Example

>>> import numpy as np
>>> from pantr.bspline import create_uniform_space, quasi_interpolate_bspline
>>> space = create_uniform_space([2], [4])
>>> qi = quasi_interpolate_bspline(lambda p: p[:, 0] ** 2, space)
>>> bool(np.isclose(qi.evaluate(np.array([[0.3]]))[0, 0], 0.09))
True

pantr.bspline.quasi_interpolate_thb_spline(func, space, *, kind='llm')[source]¶

Quasi-interpolate a callable onto a THB spline space.

Assembles the Speleers-Manni hierarchical quasi-interpolant: every active dof gets the coefficient of the per-level Lee-Lyche-Mørken local projector, sampled on a level-l active leaf cell inside the function’s support. For the truncated (THB) basis this reproduces any THB spline exactly; for the non-truncated (HB) basis it remains a valid local approximant but is not an exact projector.

Parameters:

func (Callable) – Function to quasi-interpolate. Called once on an (M, dim) point array and must return (M,) (scalar) or (M, rank) (vector-valued).
space (THBSplineSpace) – The target hierarchical space.
kind (QIKind) – Quasi-interpolant kind. Only "llm" (Lee-Lyche-Mørken) is currently supported. Defaults to "llm".

Returns:

A THB spline whose evaluation quasi-interpolates func.

Return type:

THBSpline

Raises:

TypeError – If space is not a THBSplineSpace.
ValueError – If kind is not recognized, or if func returns an output with an invalid shape (0-D, more than 2-D, or wrong leading dimension).
RuntimeError – If the grid has been modified since space was constructed, or an active dof has no leaf cell at its level (inconsistent space).

References

Hierarchical quasi-interpolation in truncated hierarchical spaces [Speleers and Manni, 2016].

Basis function evaluation for various polynomial bases.

This package provides high-level functions for tabulating basis functions (Bernstein, Lagrange, cardinal B-spline, Legendre) in 1D and multi-dimensional settings. It handles input normalization and dispatches to specialized Layer 2 and Layer 3 implementations.

LagrangeVariant: Lagrange polynomial node variants.
tabulate_bernstein_1d(), tabulate_cardinal_bspline_1d(), tabulate_lagrange_1d(), tabulate_legendre_1d(): 1D tabulation.
tabulate_bernstein(), tabulate_cardinal_bspline(), tabulate_lagrange(): multi-dimensional tabulation.

class pantr.basis.LagrangeVariant(value, names=<not given>, *values, module=None, qualname=None, type=None, start=1, boundary=None)[source]¶

Bases: StrEnum

Enumeration for Lagrange polynomial variants.

EQUISPACES = 'equispaces'¶: Equispaced points.

GAUSS_LEGENDRE = 'gauss_legendre'¶: Gauss-Legendre points (roots of Legendre polynomial).

GAUSS_LOBATTO_LEGENDRE = 'gauss_lobatto_legendre'¶: Gauss-Lobatto-Legendre points.

CHEBYSHEV_1ST = 'chebyshev_1st'¶: Chebyshev 1st kind points (x = [pi*(k + 0.5)/npts for k in range(npts)]).

CHEBYSHEV_2ND = 'chebyshev_2nd'¶: Chebyshev 2nd kind points (x = [pi*k/(npts - 1) for k in range(npts)]).

pantr.basis.tabulate_bernstein(degrees, pts, funcs_order='C', out=None)[source]¶

Evaluate the Bernstein basis functions at the given points.

Evaluates Bernstein basis functions by combining 1D basis values across each dimension, supporting both general points (e.g., a 2D array of shape (n_pts, dim) for scattered points) and point lattices (points in tensor-product grids). Fully supports C/F-ordering for functions and points.

Parameters:

degrees (Iterable[int]) – Iterable of degrees of the B-spline basis functions.
pts (npt.ArrayLike | PointsLattice) – Points at which to evaluate the basis. It can be a 2D array of shape (n_points, dim) for scattered points, or a PointsLattice object.
funcs_order (Literal["C", "F"]) – Ordering of the basis functions: ‘C’ for C-order (last index varies fastest) or ‘F’ for Fortran-order (first index varies fastest). Defaults to ‘C’.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array where the result will be stored. If None, a new array is allocated. Must have the correct shape and dtype if provided. This follows NumPy’s style for output arrays. Defaults to None.

Returns:

Array of shape (n_points, n_basis_functions) containing the combined basis function values. If out was provided, returns the same array.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If any degree is negative, or if out is provided and has incorrect shape or dtype.

pantr.basis.tabulate_bernstein_1d(degree, pts, out=None)[source]¶

Evaluate the Bernstein basis polynomials of the given degree at the given points.

Parameters:

degree (int) – Degree of the Bernstein polynomials. Must be non-negative.
pts (npt.ArrayLike) – Evaluation points. Can be a scalar, list, or numpy array. Types different from float32 or float64 are automatically converted to float64.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array where the result will be stored. If None, a new array is allocated. Must have the correct shape and dtype if provided. This follows NumPy’s style for output arrays. Defaults to None.

Returns:

Evaluated basis functions, with the same shape as the input points and the last dimension equal to (degree + 1). If out was provided, returns the same array.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If degree is negative, or if out is provided and has incorrect shape or dtype.

Example

>>> tabulate_bernstein_1d(2, [0.0, 0.5, 0.75, 1.0])
array([[1.    , 0.    , 0.    ],
       [0.25  , 0.5   , 0.25  ],
       [0.0625, 0.375 , 0.5625],
       [0.    , 0.    , 1.    ]])

pantr.basis.tabulate_cardinal_bspline(degrees, pts, funcs_order='C', out=None)[source]¶

Evaluate the cardinal B-spline basis functions at the given points.

Evaluates cardinal B-spline basis functions by combining 1D basis values across each dimension, supporting both general points (e.g., a 2D array of shape (n_pts, dim) for scattered points) and point lattices (points in tensor-product grids). Fully supports C/F-ordering for functions and points.

Parameters:

degrees (Iterable[int]) – Iterable of degrees of the cardinal B-spline basis functions.
pts (npt.ArrayLike | PointsLattice) – Points at which to evaluate the basis. It can be a 2D array of shape (n_points, dim) for scattered points, or a PointsLattice object.
funcs_order (Literal["C", "F"]) – Ordering of the basis functions: ‘C’ for C-order (last index varies fastest) or ‘F’ for Fortran-order (first index varies fastest). Defaults to ‘C’.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array where the result will be stored. If None, a new array is allocated. Must have the correct shape and dtype if provided. This follows NumPy’s style for output arrays. Defaults to None.

Returns:

Array of shape (n_points, n_basis_functions) containing the combined basis function values. If out was provided, returns the same array.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If any degree is negative, or if out is provided and has incorrect shape or dtype.

pantr.basis.tabulate_cardinal_bspline_1d(degree, pts, out=None)[source]¶

Evaluate the cardinal B-spline basis polynomials of given degree at given points.

The cardinal B-spline basis is the set of B-spline basis functions defined on an interval of maximum continuity that has degree-1 contiguous knot spans on each side with the same length as the interval itself. These basis functions appear in the central knot spans in the case of maximum regularity uniform knot vectors.

Explicit expression: [ B_{n,i}(t) = (1/n!) * sum_{j=0}^{n-i} binom(n+1, j) * (-1)^j * (t + n - i - j)^n ] where ( B_{n,i}(t) ) is the B-spline basis function of degree ( n ) and index ( i ) at point ( t ), and ( binom(a, b) ) is the binomial coefficient.

Its actual implementation is based on the Cox-de Boor recursion formula for the central cardinal B-spline basis.

Parameters:

degree (int) – Degree of the B-spline basis. Must be non-negative.
pts (npt.ArrayLike) – Evaluation points. Can be a scalar, list, or numpy array. Types different from float32 or float64 are automatically converted to float64.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array where the result will be stored. If None, a new array is allocated. Must have the correct shape and dtype if provided. This follows NumPy’s style for output arrays. Defaults to None.

Returns:

Evaluated basis functions, with the same shape as the input points and the last dimension equal to (degree + 1). If out was provided, returns the same array.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If provided degree is negative, or if out is provided and has incorrect shape or dtype.

Example

>>> tabulate_cardinal_bspline_1d(2, [0.0, 0.5, 1.0])
array([[0.5  , 0.5  , 0.   ],
       [0.125, 0.75 , 0.125],
       [0.   , 0.5  , 0.5  ]])

pantr.basis.tabulate_lagrange(degrees, variant, pts, funcs_order='C', out=None)[source]¶

Evaluate the Lagrange basis functions at the given points.

Evaluates Lagrange basis functions by combining 1D basis values across each dimension, supporting both general points (e.g., a 2D array of shape (n_pts, dim) for scattered points) and point lattices (points in tensor-product grids). Fully supports C/F-ordering for functions and points.

Parameters:

degrees (Iterable[int]) – Iterable of degrees of the Lagrange basis functions.
variant (LagrangeVariant) – Variant of the Lagrange basis.
pts (npt.ArrayLike | PointsLattice) – Points at which to evaluate the basis. It can be a 2D array of shape (n_points, dim) for scattered points, or a PointsLattice object.
funcs_order (Literal["C", "F"]) – Ordering of the basis functions: ‘C’ for C-order (last index varies fastest) or ‘F’ for Fortran-order (first index varies fastest). Defaults to ‘C’.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array where the result will be stored. If None, a new array is allocated. Must have the correct shape and dtype if provided. This follows NumPy’s style for output arrays. Defaults to None.

Returns:

Array of shape (n_points, n_basis_functions) containing the combined basis function values. If out was provided, returns the same array.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If any degree is negative, or if out is provided and has incorrect shape or dtype.

pantr.basis.tabulate_lagrange_1d(degree, variant, pts, out=None)[source]¶

Evaluate Lagrange basis polynomials at points using the specified variant.

The polynomials are defined in the interval [0, 1] and are given by the formula: [ L_{n,i}(t) = prod_{j=0,, j neq i}^{n} frac{t - x_j}{x_i - x_j} ] where ( x_i ) are the points at which the basis is evaluated.

The variant determines the points at which the basis is evaluated.

Parameters:

degree (int) – Degree of the Lagrange basis. Must be non-negative.
variant (LagrangeVariant) – Variant of the Lagrange basis.
pts (npt.ArrayLike) – Evaluation points. Can be a scalar, list, or numpy array. Types different from float32 or float64 are automatically converted to float64.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array where the result will be stored. If None, a new array is allocated. Must have the correct shape and dtype if provided. This follows NumPy’s style for output arrays. Defaults to None.

Returns:

Evaluated basis functions, with the same shape as the input points and the last dimension equal to (degree + 1). If out was provided, returns the same array.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If provided degree is negative, or if out is provided and has incorrect shape or dtype.

pantr.basis.tabulate_legendre_1d(degree, pts, out=None)[source]¶

Evaluate the normalized Shifted Legendre basis polynomials at the given points.

The polynomials are defined on the interval [0, 1] and are orthonormal: [ int_0^1 tilde{p}_n(x) tilde{p}_m(x) , dx = delta_{nm} ]

Parameters:

degree (int) – Degree of the Legendre basis. Must be non-negative.
pts (npt.ArrayLike) – Evaluation points. Can be a scalar, list, or numpy array. Types different from float32 or float64 are automatically converted to float64.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array where the result will be stored. If None, a new array is allocated. Must have the correct shape and dtype if provided. This follows NumPy’s style for output arrays. Defaults to None.

Returns:

Evaluated basis functions, with the same shape as the input points and the last dimension equal to (degree + 1). If out was provided, returns the same array.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If provided degree is negative, or if out is provided and has incorrect shape or dtype.

Change of basis operators for various polynomial bases in 1D.

This module provides functions to create transformation matrices between different polynomial bases including Lagrange, Bernstein, cardinal B-spline, and monomial bases.

Architecturally, this module serves as the bridge between different basis types, providing pure mathematical functions to compute the exact $(degree+1, degree+1)$ transformation matrices (e.g., $M$ such that $new\_basis(x) = M @ old\_basis(x)$) without tying the dense numerical quadrature logic directly into the core Spline space objects.

pantr.change_basis.compute_bernstein_to_cardinal_1d(degree, dtype=<class 'numpy.float64'>, out=None)[source]¶

Create transformation matrix from Bernstein to cardinal B-spline basis.

Parameters:

degree (int) – Polynomial degree. Must be non-negative.
dtype (npt.DTypeLike) – Floating point type for the output matrix. Defaults to np.float64.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array where the result will be stored. If None, a new array is allocated. Must have shape (degree+1, degree+1) and dtype matching the dtype parameter if provided. This follows NumPy’s style for output arrays. Defaults to None.

Returns:

(degree+1, degree+1) transformation matrix C such that: C @ [Bernstein values] = [cardinal values]. If out was provided, returns the same array.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If degree is negative, dtype is not float32 or float64, or if out is provided and has incorrect shape or dtype.

pantr.change_basis.compute_bernstein_to_lagrange_1d(degree, lagrange_variant=LagrangeVariant.EQUISPACES, dtype=<class 'numpy.float64'>, out=None)[source]¶

Construct the matrix mapping Bernstein basis evaluations to Lagrange basis evaluations.

Note

Both Bernstein and Lagrange bases follow the standard ordering (see https://en.wikipedia.org/wiki/Bernstein_polynomial).

Parameters:

degree (int) – Polynomial degree. Must be at least 1.
lagrange_variant (LagrangeVariant) – Lagrange point distribution (e.g., equispaced, gauss lobatto legendre, etc). Defaults to LagrangeVariant.EQUISPACES.
dtype (npt.DTypeLike) – Floating point type for the output matrix. Defaults to np.float64.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array where the result will be stored. If None, a new array is allocated. Must have shape (degree+1, degree+1) and dtype matching the dtype parameter if provided. This follows NumPy’s style for output arrays. Defaults to None.

Returns:

(degree+1, degree+1) transformation matrix C such that: C @ [Bernstein values] = [Lagrange values]. If out was provided, returns the same array.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If degree is lower than 1, dtype is not float32 or float64, or if out is provided and has incorrect shape or dtype.

pantr.change_basis.compute_cardinal_to_bernstein_1d(degree, dtype=<class 'numpy.float64'>, out=None)[source]¶

Create transformation matrix from cardinal B-spline to Bernstein basis.

Parameters:

degree (int) – Polynomial degree. Must be non-negative.
dtype (npt.DTypeLike) – Floating point type for the output matrix. Defaults to np.float64.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array where the result will be stored. If None, a new array is allocated. Must have shape (degree+1, degree+1) and dtype matching the dtype parameter if provided. This follows NumPy’s style for output arrays. Defaults to None.

Returns:

(degree+1, degree+1) transformation matrix C such that: C @ [cardinal values] = [Bernstein values]. If out was provided, returns the same array.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If degree is negative, dtype is not float32 or float64, or if out is provided and has incorrect shape or dtype.

pantr.change_basis.compute_lagrange_to_bernstein_1d(degree, lagrange_variant=LagrangeVariant.EQUISPACES, dtype=<class 'numpy.float64'>, out=None)[source]¶

Construct the matrix mapping Lagrange basis evaluations to Bernstein basis evaluations.

Note

Both Bernstein and Lagrange bases follow the standard ordering (see https://en.wikipedia.org/wiki/Bernstein_polynomial).

Parameters:

degree (int) – Polynomial degree. Must be at least 1.
lagrange_variant (LagrangeVariant) – Lagrange point distribution (e.g., equispaced, gauss lobatto legendre, etc). Defaults to LagrangeVariant.EQUISPACES.
dtype (npt.DTypeLike) – Floating point type for the output matrix. Defaults to np.float64.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array where the result will be stored. If None, a new array is allocated. Must have shape (degree+1, degree+1) and dtype matching the dtype parameter if provided. This follows NumPy’s style for output arrays. Defaults to None.

Returns:

(degree+1, degree+1) transformation matrix C such that: C @ [Lagrange values] = [Bernstein values]. If out was provided, returns the same array.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If degree is lower than 1, dtype is not float32 or float64, or if out is provided and has incorrect shape or dtype.

pantr.change_basis.compute_monomial_to_bernstein_1d(degree, dtype=<class 'numpy.float64'>, out=None)[source]¶

Create transformation matrix from monomial to Bernstein basis on [0, 1].

Given a polynomial of degree degree written in the monomial basis on [0, 1], the returned matrix M converts its coefficient vector to the Bernstein basis: bern_coeffs = M @ mono_coeffs. Equivalently, on basis evaluations, monomial(x) = M.T @ bernstein(x).

The entries are M[i, j] = C(i, j) / C(degree, j) for j <= i, else 0, where C(n, k) is the binomial coefficient.

Parameters:

degree (int) – Polynomial degree. Must be non-negative.
dtype (npt.DTypeLike) – Floating point type for the output matrix. Defaults to np.float64.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array where the result will be stored. If None, a new array is allocated. Must have shape (degree+1, degree+1) and dtype matching the dtype parameter if provided. This follows NumPy’s style for output arrays. Defaults to None.

Returns:

(degree+1, degree+1) lower-triangular: transformation matrix M such that M @ [monomial coefficients] = [Bernstein coefficients]. If out was provided, returns the same array.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If degree is negative, dtype is not float32 or float64, or if out is provided and has incorrect shape or dtype.

Bezier geometric objects, approximation functions, and Bernstein polynomial root finding.

This module provides Bezier, which stores control points to represent a parametric Bezier curve, surface, or volume. Degree is derived from the control point array shape. Evaluation and manipulation methods use direct Bernstein algorithms implemented in _bezier_core, _bezier_eval, _bezier_derivative, _bezier_degree, and _bezier_product.

interpolate_bezier(), fit_bezier(): approximation functions.
create_from_bspline(): create a Bézier from a B-spline.

Root-finding exports:

find_roots() – find all roots (single or batch, auto-dispatch).
find_monotone_root() – fast solver for monotone polynomials (single or batch).

class pantr.bezier.Bezier(control_points, is_rational=False)[source]¶

Bases: object

A parametric Bézier curve, surface, or volume defined by control points.

Stores only control points and an is_rational flag. The polynomial degree in each parametric direction is inferred from the control point array shape: degree[d] = control_points.shape[d] - 1.

Parameters:

control_points (npt.ArrayLike)
is_rational (bool)

_control_points¶

Control point array with shape (*degrees_plus_1, rank), where the last axis is the output rank (including the homogeneous weight for rational).

Type:: npt.NDArray[np.float32 | np.float64]

_is_rational¶

Whether the Bézier is rational (last control point coordinate is a homogeneous weight).

Type:: bool

__init__(control_points, is_rational=False)[source]¶

Initialize a Bézier from control points.

Parameters:

control_points (npt.ArrayLike) – Control points. Shape must be at least 2D: (*degrees_plus_1, rank). A 1D input of shape (n,) is reshaped to (n, 1) (scalar field). Integer arrays are cast to float64.
is_rational (bool) – Whether the Bézier is rational. Defaults to False.

Raises:

ValueError – If the control points have fewer than 1 entry in any parametric direction.
ValueError – If the Bézier has rank smaller than 1.

Return type:

None

property dim: int¶

Get the parametric dimension of the Bézier.

Returns:: Number of parametric dimensions.
Return type:: int

property degree: tuple[int, ...]¶

Get the polynomial degrees per parametric direction.

Returns:: Polynomial degree for each parametric dimension, computed as shape[d] - 1.
Return type:: tuple[int, …]

property control_points: ndarray[tuple[Any, ...], dtype[float32 | float64]]¶

Get the control points of the Bézier.

Returns:: Control point array with shape (*degrees_plus_1, rank).
Return type:: npt.NDArray[np.float32 | np.float64]

property is_rational: bool¶

Check whether the Bézier is rational.

Returns:: True if the Bézier is rational (i.e., the last control point coordinate is a homogeneous weight), False otherwise.
Return type:: bool

property rank: int¶

Get the output rank of the Bézier.

The rank is the number of value dimensions produced by the mapping. For a scalar field it is 1; for a 3D curve it is 3. For rational Béziers the weight coordinate is excluded.

Returns:: Output rank of the Bézier.
Return type:: int

property dtype: DTypeLike¶

Get the floating-point dtype of the Bézier.

Returns:: The numpy dtype (float32 or float64) of the control point array.
Return type:: npt.DTypeLike

evaluate(pts, out=None)[source]¶

Evaluate the Bézier at the given parametric points.

Parameters:

pts (npt.NDArray[np.float32 | np.float64] | PointsLattice) – The parametric points at which to evaluate. For 1D Bézier, must be a 1D array of shape (n_pts,). For multi-dimensional Bézier, must be a 2D array of shape (n_pts, dim) or a PointsLattice.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional output array. Defaults to None.

Returns:

Bézier values at the given points.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If the points dtype does not match the Bézier dtype, or if out has incorrect shape or dtype.

evaluate_derivatives(pts, orders, out=None)[source]¶

Evaluate a specific partial derivative of the Bézier.

Computes the single partial derivative specified by orders, where orders[d] is the derivative order in parametric direction d. For rational Bézier the generalised quotient rule is applied.

Parameters:

pts (npt.NDArray[np.float32 | np.float64] | PointsLattice) – The parametric points at which to evaluate.
orders (int | Sequence[int]) – Derivative order(s). A single int is broadcast to all self.dim directions. A sequence must contain one non-negative integer per parametric direction.
out (npt.NDArray[np.float32 | np.float64] | None) – Optional pre-allocated output array. Defaults to None.

Returns:

Mixed partial derivative values.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If len(orders) != self.dim, if any order is negative, or if the points dtype does not match the Bézier dtype.

Example

>>> result = bezier.evaluate_derivatives(pts, 2)
>>> result = bezier.evaluate_derivatives(pts, [1, 2])

derivative(direction=0, *, keep_degree=False)[source]¶

Return a Bézier representing the first derivative in the given direction.

Computes the hodograph: a new Bézier whose value at every parametric point equals the partial derivative of this Bézier with respect to parametric direction direction.

For non-rational Bézier of degree p in direction d, the result has degree p - 1 (or p when keep_degree=True).

For rational Bézier, the quotient rule is applied, producing a rational Bézier of degree 2p in direction d (or the original degree when keep_degree=True).

Parameters:

direction (int) – Parametric direction for differentiation. Must be in [0, dim). Defaults to 0.
keep_degree (bool) – If True, the result preserves the same degree as the original Bézier by fusing derivative and degree elevation. This is useful, for instance, when computing derivatives of rational polynomials (in the numerator). Defaults to False.

Returns:

A new Bézier representing the derivative.

Return type:

Raises:

ValueError – If direction is out of range [0, dim).
ValueError – If the degree in the given direction is 0.

Example

>>> f_prime = f.derivative()
>>> df_dv = surface.derivative(direction=1)
>>> f_prime_same_deg = f.derivative(keep_degree=True)

elevate_degree(degree_increments)[source]¶

Elevate the polynomial degree of the Bézier.

Creates a new Bézier that represents the same mapping as the original but with higher polynomial degree.

Parameters:

Returns:

A new Bézier with elevated degrees.

Return type:

Raises:

ValueError – If any degree increment is negative.
ValueError – If all degree increments are zero.
ValueError – If the number of increments does not match the dimension.

reduce_degree(degree_decrements)[source]¶

Reduce the polynomial degree of the Bézier via least-squares approximation.

Creates a new Bézier whose degree is lower by the requested amount in each parametric direction. The reduction minimises the squared error under the Bernstein degree-elevation matrix using QR factorisation with Givens rotations.

Unlike elevate_degree(), this operation is not exact in general: the result is an approximation of the original mapping.

Parameters:

Returns:

A new Bézier with reduced degrees.

Return type:

Raises:

ValueError – If any degree decrement is negative.
ValueError – If all degree decrements are zero.
ValueError – If the number of decrements does not match the dimension.
ValueError – If any decrement exceeds the current degree in that direction.

minimize_degree(tol=None)[source]¶

Find the lowest degree that preserves accuracy within tolerance.

Iterates over each parametric direction and repeatedly tries to reduce the degree by 1. A reduction is accepted when the round-trip (reduce then elevate) relative L2 error stays below tol. For vector-valued Bézier, all rank components are checked simultaneously.

Parameters:: tol (float | None) – Relative tolerance for accepting a degree reduction. If None, uses a default based on machine epsilon (1e3 * eps).
Returns:: A new Bézier with the lowest degree that preserves accuracy. If no reduction is possible, returns a copy.
Return type:: Bezier

Example

>>> import numpy as np
>>> b = Bezier(np.array([3.0, 3.0, 3.0, 1.0]).reshape(4, 1))
>>> b.degree
(3,)
>>> b_min = b.minimize_degree()

multiply(other)[source]¶

Return the exact pointwise product of this Bézier and another.

Given Bézier self and other over the same parametric domain [0, 1]^dim, returns a new Bézier h such that h(t) = self(t) * other(t). The result has degree p_d + q_d per direction.

Parameters:: other (Bezier) – The second Bézier operand. Must have the same dimension, dtype, and rank as self.
Returns:: A new Bézier representing self * other.
Return type:: Bezier
Raises:: ValueError – If the operands have different dimensions, dtypes, or ranks.

Example

>>> h = f.multiply(g)
>>> h2 = f * g

compose(inner)[source]¶

Compose this Bézier with another: result(t) = self(inner(t)).

Computes the exact composition of two non-rational Bézier objects. The result is a new Bézier with parametric dimension equal to inner.dim, rank equal to self.rank, and degree sum(self.degree) * inner.degree[s] in each direction s.

Parameters:

inner (Bezier) – The inner Bézier (reparametrization). Must be non-rational and satisfy inner.rank == self.dim.

Returns:

A new Bézier representing self(inner(t)).

Return type:

Raises:

TypeError – If either Bézier is rational.
ValueError – If inner.rank != self.dim.
ValueError – If the operands have different dtypes.

Example

>>> f = Bezier(np.array([[0.0, 0.0], [1.0, 1.0], [2.0, 0.0]]))
>>> g = Bezier(np.array([[0.2], [0.8]]))
>>> h = f.compose(g)

reverse(direction: int = 0, *, in_place: Literal[False] = False) → Bezier[source]¶

reverse(direction: int = 0, *, in_place: Literal[True]) → None

Reverse the orientation of one parametric direction.

Flips the control points along the given parametric axis so that the mapping is reparametrised in the opposite sense along that direction.

Parameters:

direction (int) – Parametric direction to reverse. Must be in [0, dim). Defaults to 0.
in_place (bool) – If True, modify this Bézier in place and return None. If False (default), return a new Bézier.

Returns:

The reversed Bézier, or None when in_place=True.

Return type:

Bezier | None

Raises:

ValueError – If direction is out of range [0, dim).

Example

>>> rev = bezier.reverse(direction=0)
>>> bezier.reverse(direction=1, in_place=True)

permute_directions(permutation: Sequence[int], *, in_place: Literal[False] = False) → Bezier[source]¶

permute_directions(permutation: Sequence[int], *, in_place: Literal[True]) → None

Reorder the parametric directions according to a permutation.

Parameters:

permutation (Sequence[int]) – A permutation of range(dim).
in_place (bool) – If True, modify this Bézier in place and return None. If False (default), return a new Bézier.

Returns:

The permuted Bézier, or None when in_place=True.

Return type:

Bezier | None

Raises:

ValueError – If permutation is not a valid permutation of range(dim).

Example

>>> surface.permute_directions([1, 0])  # swap u ↔ v

transform(affine: AffineTransform, *, in_place: Literal[False] = False) → Bezier[source]¶

transform(affine: AffineTransform, *, in_place: Literal[True]) → None

Apply an affine transformation to the control points.

For non-rational Bézier, every control point P is mapped to A @ P + b. For rational Bézier the weighted homogeneous coordinates are updated so that the projected geometry undergoes the same affine map while the weights are preserved.

Parameters:

affine (AffineTransform) – The affine transformation to apply.
in_place (bool) – If True, the control points are modified in place and None is returned. If False (default), a new Bezier is returned.

Returns:

The transformed Bézier, or None when in_place=True.

Return type:

Bezier | None

Raises:

ValueError – If the transform dimension does not match the geometric rank of the Bézier.

Example

>>> from pantr.transform import AffineTransform
>>> T = AffineTransform.translation([1.0, 2.0])
>>> shifted = bezier.transform(T)

restrict(bounds)[source]¶

Return a Bézier restricted to a sub-region of [0, 1]^dim.

Extracts the portion of the Bézier defined on the given sub-interval and reparametrizes the result back to [0, 1]^dim. The returned Bézier has the same degree but different control points that encode the restricted mapping.

Uses two de Casteljau passes per direction for direct Bernstein coefficient computation without B-spline conversion. The order of the passes is chosen for numerical stability.

Parameters:

bounds (tuple[float, float] | Sequence[tuple[float, float] | None]) – For a 1D Bézier, a (lower, upper) tuple within [0, 1]. For multi-dimensional Bézier, a sequence of length dim where each element is a (lower, upper) tuple for that direction, or None to keep the full [0, 1] range. At least one direction must have non-None bounds that restrict the domain.

Returns:

New Bézier on [0, 1]^dim representing the restriction.

Return type:

Raises:

ValueError – If the sequence length does not match dim.
ValueError – If all directions are None or match the full domain.
ValueError – If any bound lies outside [0, 1].
ValueError – If lower >= upper in any direction.

split(direction, value)[source]¶

Split the Bézier into two at a parameter value in one direction.

Uses the de Casteljau algorithm to subdivide the Bézier into a left half (representing the original on [0, value]) and a right half (representing the original on [value, 1]), both reparametrized to [0, 1].

Parameters:

direction (int) – Parametric direction along which to split. Must be in [0, dim).
value (float) – Parameter value at which to split. Must lie strictly inside (0, 1).

Returns:

A pair (left, right) of Béziers on [0, 1]^dim.

Return type:

tuple[Bezier, Bezier]

Raises:

ValueError – If direction is out of range [0, dim).
ValueError – If value is not strictly inside (0, 1).

Example

>>> left, right = curve.split(0, 0.5)
>>> left, right = surface.split(1, 0.3)

slice(axis, value)[source]¶

Slice the Bézier by fixing one parametric direction at a given value.

Reduces the parametric dimension by one using the de Casteljau algorithm on the control points. A surface becomes a curve, a curve becomes a point (returned as a NumPy array).

At the boundary values 0 and 1 the result is obtained in O(1) by direct control point lookup.

Parameters:

axis (int) – Parametric direction to fix (0-indexed). Must be in [0, dim).
value (float) – Parameter value at which to slice. Must lie within [0, 1].

Returns:

A Bézier with dim - 1 dimensions when dim >= 2, or a NumPy array of shape (rank,) when dim == 1. Rational Béziers preserve the rational structure when dim >= 2; for dim == 1 the result is projected to physical coordinates.

Return type:

Bezier | npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If axis is out of range [0, dim).
ValueError – If value is outside [0, 1].

Example

>>> # Slice a surface at v=0.5 to get a curve
>>> curve = surface.slice(1, 0.5)
>>> # Composable: surface -> curve -> point
>>> pt = surface.slice(1, 0.5).slice(0, 0.2)

boundary(axis, side)[source]¶

Extract the boundary of the Bézier along one parametric direction.

Returns the restriction of the Bézier to one end of the [0, 1] domain in the given direction.

Parameters:

axis (int) – Parametric direction (0-indexed). Must be in [0, dim).
side (int) – Which end of the domain: 0 for the start, 1 for the end.

Returns:

A Bézier with dim - 1 dimensions when dim >= 2, or a NumPy array of shape (rank,) when dim == 1.

Return type:

Bezier | npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If axis is out of range [0, dim).
ValueError – If side is not 0 or 1.

Example

>>> # Extract left boundary of a surface along direction 0
>>> left_curve = surface.boundary(0, 0)

collapse_along_axis(axis, values)[source]¶

Collapse to a univariate Bézier along one parametric direction.

Fixes all parametric directions except axis at the given parameter values, producing a 1D Bézier whose control points are the Bernstein coefficients along axis. This is a tensor contraction: for each collapsed direction, the Bernstein basis is evaluated at the given value and contracted with the control point array.

Parameters:

axis (int) – Parametric direction to keep (0-indexed). Must be in [0, dim).
values (npt.ArrayLike) – Parameter values for all directions except axis. Must have length dim - 1 with all values in [0, 1]. values[i] corresponds to direction i for i < axis, and direction i + 1 for i >= axis.

Returns:

A 1D Bézier with degree self.degree[axis] and the same rank and rationality as the input.

Return type:

Raises:

ValueError – If dim < 2 (nothing to collapse).
ValueError – If axis is out of range [0, dim).
ValueError – If values does not have length dim - 1.
ValueError – If any value is outside [0, 1].

Example

>>> # Collapse a 3D volume along axis 1 at (u=0.3, w=0.7)
>>> curve = volume.collapse_along_axis(1, [0.3, 0.7])

to_bspline(*, copy=True)[source]¶

Convert to an equivalent B-spline with Bézier knot vectors.

Creates a Bspline with open knot vectors [0]*(p+1) + [1]*(p+1) in each parametric direction.

Parameters:: copy (bool) – If True (default), the control points are deep-copied into the new B-spline. If False, the B-spline shares the same underlying control point array.
Returns:: Equivalent B-spline representation.
Return type:: Bspline

plot(*, color=None, show_control_polygon=False, **plotter_kwargs)[source]¶

Quick interactive visualization of this Bézier (requires pyvista).

For finer control, use pantr.viz.Scene directly.

Parameters:

color (str | None) – Surface color.
show_control_polygon (bool) – Render control polygon (points and wireframe).
**plotter_kwargs (Any) – Additional keyword arguments for pv.Plotter().

Returns:

The pyvista Plotter after showing.

Return type:

object

Raises:

ImportError – If pyvista is not installed.

pantr.bezier.create_from_bspline(bspline, *, copy=True)[source]¶

Create a Bézier from a B-spline with Bézier-like knot vectors.

Validates that the B-spline has Bézier-like knots (open knots with num_basis == degree + 1 in each direction) and extracts the control points.

Parameters:

bspline (Bspline) – A B-spline with Bézier-like knot structure.
copy (bool) – If True (default), the control points are deep-copied into the new Bézier. If False, the Bézier shares the same underlying control point array.

Returns:

The equivalent Bézier.

Return type:

Raises:

ValueError – If the B-spline does not have Bézier-like knots.

pantr.bezier.find_monotone_root(bezier, *, tol=None)[source]¶

Find the unique root of one or more monotone scalar Bezier curves in [0, 1].

Uses a Newton/bisection hybrid with false-position initialization. Each Bezier must be monotone on [0, 1] (i.e. its derivative does not change sign).

When called with a single Bezier, returns a float. When called with a sequence of Beziers (batch mode), all curves must have the same degree and the function returns an array.

Parameters:

bezier (Bezier | Sequence[Bezier]) – A single 1-D (dim == 1) scalar (rank == 1) Bezier curve, or a sequence of such curves (all with the same degree). For rational Beziers, roots are found on the numerator polynomial.
tol (float | None) – Parameter-space termination tolerance. Defaults to tolerance.get_strict(bezier.dtype).

Returns:

(single mode) Root parameter in [0, 1], or NaN if no root: exists (no sign change across the interval).
npt.NDArray[np.float64]: (batch mode) 1-D array of shape: (n_polys,) with root values. Contains NaN where no root exists. Always float64.

Return type:

float

Raises:

TypeError – If bezier is not a Bezier instance (or sequence thereof).
ValueError – If any Bezier has dim != 1 or rank != 1, or tol is not positive. In batch mode, also raised if degrees are not uniform.

pantr.bezier.find_roots(bezier, *, tol=None)[source]¶

Find all roots of one or more scalar Bezier curves in [0, 1].

Auto-selects between Yuksel’s monotone-decomposition algorithm and Bezier clipping based on polynomial degree and coefficient dynamic range.

When called with a single Bezier, returns a sorted array of root parameters. When called with a sequence of Beziers (batch mode), all curves must have the same degree and the function returns a (roots, counts) tuple.

Parameters:

bezier (Bezier | Sequence[Bezier]) – A single 1-D (dim == 1) scalar (rank == 1) Bezier curve, or a sequence of such curves (all with the same degree). For rational Beziers, roots are found on the numerator polynomial (first homogeneous component).
tol (float | None) – Root-finding tolerance (bracket-width termination). Defaults to tolerance.get_strict(bezier.dtype).

Returns:

(single mode) Sorted array of root

parameters in [0, 1]. Empty if no roots exist. Always float64 regardless of input dtype.

tuple[npt.NDArray[np.float64], npt.NDArray[np.intp]]: (batch mode)

roots: padded array of shape (n_polys, max(degree, 1)) where only the first counts[i] entries per row are valid. Always float64.
counts: 1-D array of shape (n_polys,) with the number of valid roots per polynomial.

Return type:

npt.NDArray[np.float64]

Raises:

TypeError – If bezier is not a Bezier instance (or sequence thereof).
ValueError – If any Bezier has dim != 1 or rank != 1, or tol is not positive. In batch mode, also raised if degrees are not uniform.

Note

In batch mode, a simpler fixed-threshold dedup (no derivative-aware merge radius) is used compared to single mode. In rare edge cases involving near-duplicate roots, the two modes may report different root counts.

Example

>>> import numpy as np
>>> from pantr.bezier import Bezier, find_roots
>>> find_roots(Bezier([1.0, -1.0]))
array([0.5])

References

Yuksel’s monotone-decomposition root finder [Yuksel, 2022] and the Bézier clipping method [Sederberg and Nishita, 1990].

pantr.bezier.fit_bezier(values, nodes, *, degree=None, tol=None)[source]¶

Construct a Bézier from pre-evaluated sample values at known nodes.

The output dtype is inferred from values.

Supports two point layouts:

Tensor-product (a PointsLattice, a single 1D array, or a sequence of 1D arrays): values must have shape (*n_pts_per_dir) (scalar) or (*n_pts_per_dir, rank) (vector). Per-direction SVD pseudo-inverse is used (efficient).
Scattered (a 2D ndarray of shape (n_pts, dim)): values must have shape (n_pts,) (scalar) or (n_pts, rank) (vector). degree is required. A full tensor-product Vandermonde is built and solved via SVD.

Parameters:

values (npt.ArrayLike) – Sample values.
nodes (PointsLattice | npt.NDArray[np.floating[Any]] | Sequence[npt.NDArray[np.floating[Any]]]) –
Interpolation nodes.
- A PointsLattice: tensor-product grid.
- A 1D ndarray: 1D tensor-product (single direction).
- A sequence of 1D ndarray values: N-D tensor-product.
- A 2D ndarray of shape (n_pts, dim): scattered points.
degree (int | Sequence[int] | None) – Polynomial degree per direction. Required for scattered nodes. If None (default) for tensor-product nodes, degree = n_pts - 1 (exact interpolation). Must satisfy prod(degree + 1) <= n_pts for scattered nodes and degree < n_pts per direction for tensor-product.
tol (float | None) – SVD truncation tolerance. If None, uses a default based on machine epsilon.

Returns:

A non-rational Bézier.

Return type:

Bezier

Raises:

ValueError – If nodes are inconsistent with values, degree is invalid, or degree is missing for scattered nodes.

pantr.bezier.interpolate_bezier(func, n_pts, *, degree=None, nodes=None, tol=None)[source]¶

Interpolate a callable function into a Bézier in Bernstein form.

Evaluates func on a tensor-product grid of interpolation nodes and recovers the Bernstein coefficients via truncated SVD.

The parametric dimension is determined by the length of n_pts (when given as a sequence) or defaults to 1 (when given as a scalar int).

The output dtype is inferred from the return value of func.

Parameters:

func (Callable[..., npt.ArrayLike]) – Function to interpolate. Called as func(lattice) where lattice is a PointsLattice representing the tensor-product sampling grid. The callable may also accept a plain ndarray (for compatibility with Bezier.evaluate() signatures). Must return an array of shape (n_total,) for a scalar-valued function or (n_total, rank) for a vector-valued function, where n_total = prod(n_pts).
n_pts (int | Sequence[int]) – Number of sample points per parametric direction. A single int gives a 1D Bézier.
degree (int | Sequence[int] | None) – Polynomial degree per direction. If None (default), degree = n_pts - 1 (exact interpolation). If provided, must satisfy degree < n_pts in each direction; the result is a least-squares approximation.
nodes (Literal['chebyshev', 'uniform'] | PointsLattice | npt.NDArray[np.floating[Any]] | Sequence[npt.NDArray[np.floating[Any]]] | None) –
Interpolation node selection.
- None or "chebyshev" (default): modified Chebyshev-Lobatto nodes on [0, 1].
- "uniform": equispaced nodes on [0, 1].
- A PointsLattice: custom tensor-product grid.
- A 1D ndarray: custom nodes broadcast to all directions.
- A sequence of 1D ndarray values: per-direction custom nodes.
tol (float | None) – SVD truncation tolerance. If None, uses a default based on machine epsilon.

Returns:

A non-rational Bézier whose evaluation approximates func.

Return type:

Bezier

Raises:

ValueError – If n_pts values are < 1, degree >= n_pts, nodes is inconsistent with n_pts, or the callable returns an unexpected shape.

Example

>>> from pantr.bezier import Bezier
>>> import numpy as np
>>> b = interpolate_bezier(lambda lattice: lattice.get_all_points()[:, 0] ** 2, [5])
>>> b.degree
(4,)

Geometry & CAD¶

Constructive geometry for B-spline curves, surfaces, and volumes.

This package provides CAD-style functions for creating and combining B-spline geometric objects:

Primitives: create_line(), create_circle(), create_bilinear(), create_trilinear().
Derived primitives: create_rectangle(), create_disk(), create_cylinder().
Compatibility: make_compat().
Operations: create_extrusion(), create_revolution(), create_ruled(), create_sweep().
Coons blending: create_coons_surface(), create_coons_volume().
Assembly: join().

pantr.cad.create_bilinear(corners=None)[source]¶

Construct a bilinear B-spline surface from four corner points.

Creates a degree (1, 1), non-rational B-spline surface with 2 x 2 control points.

The corner ordering follows the tensor-product convention:

corners[0, 1]       corners[1, 1]
o------------------o
|  v               |
|  ^               |
|  |               |
|  +-------> u     |
o------------------o
corners[0, 0]       corners[1, 0]

Parameters:: corners (ArrayLike | None) – Array of shape (2, 2, rank) with rank <= 3. Coordinates shorter than 3D are zero-padded. If None, defaults to the unit square [-0.5, 0.5]^2 x {0} in the xy-plane.
Returns:: A 2D, degree-(1, 1), rank-3, non-rational B-spline surface.
Return type:: Bspline
Raises:: ValueError – If corners does not have shape (2, 2, rank) with 1 <= rank <= 3.

Example

>>> srf = create_bilinear()
>>> srf.degree
(1, 1)
>>> srf.dim
2

pantr.cad.create_circle(radius=1.0, center=None, angle=None)[source]¶

Construct a NURBS circular arc or full circle.

Creates a degree-2 rational B-spline curve in the xy-plane. The arc is split into spans of at most 90 degrees each. Interior knots have multiplicity 2 (equal to the degree), giving C0 continuity at arc junctions. This is the standard exact representation of conics using rational quadratic B-splines.

The number of spans depends on the sweep angle:

|sweep| <= 90: 1 span, 3 control points
90 < |sweep| <= 180: 2 spans, 5 control points
180 < |sweep| <= 270: 3 spans, 7 control points
270 < |sweep| <= 360: 4 spans, 9 control points

Parameters:

radius (float) – Circle radius. Defaults to 1.
center (ArrayLike | None) – Center point (up to 3D, zero-padded). If None, the circle is centered at the origin.
angle (float | tuple[float, float] | None) –
Sweep specification.
- None – full circle (360 degrees).
- float – arc from angle 0 to the given value (radians).
- (start, end) – arc from start to end (radians).

Returns:

A 1D, degree-2, rank-3, rational B-spline curve.

Return type:

Example

>>> crv = create_circle()
>>> crv.degree
(2,)
>>> crv.is_rational
True

pantr.cad.create_coons_surface(curves)[source]¶

Construct a Coons patch from four boundary curves.

Builds a bilinearly blended surface from four boundary curves using the formula S = R0 + R1 - B, where R0 and R1 are ruled surfaces and B is the bilinear interpolant of the four corners.

The curve layout is:

       C_u1
o--------------o
|  v           |
|  ^           |
|  |     C_v1  |
|  +---> u     |
o--------------o
C_v0   C_u0

Parameters:

curves (tuple[tuple[Bspline, Bspline], tuple[Bspline, Bspline]]) – ((C_v0, C_v1), (C_u0, C_u1)) – two pairs of opposite boundary curves. C_v0 and C_v1 run in the v-direction (left/right). C_u0 and C_u1 run in the u-direction (bottom/top). All must be 1D B-splines.

Returns:

A 2D B-spline surface.

Return type:

Raises:

ValueError – If any curve is not 1D.
ValueError – If corner points are not geometrically consistent.

pantr.cad.create_coons_volume(faces)[source]¶

Construct a trivariate Coons blending volume from 6 boundary faces.

Uses the trilinear Coons formula:

V = (R_u + R_v + R_w) - (B_uv + B_uw + B_vw) + T

where R are ruled volumes from opposite face pairs, B are bilinear blend volumes from edge quadruples, and T is the trilinear corner interpolant.

Edges and corners are derived automatically from face boundaries.

Face labelling convention:

faces[0] = (face_u0, face_u1): faces at u=0 and u=1, each a surface parameterized by (v, w).
faces[1] = (face_v0, face_v1): faces at v=0 and v=1, each a surface parameterized by (u, w).
faces[2] = (face_w0, face_w1): faces at w=0 and w=1, each a surface parameterized by (u, v).

Parameters:: faces (tuple[tuple[Bspline, Bspline], tuple[Bspline, Bspline], tuple[Bspline, Bspline]]) – Three pairs of opposite boundary faces.
Returns:: A 3D (trivariate) B-spline volume.
Return type:: Bspline
Raises:: ValueError – If any face is not a 2D B-spline.

pantr.cad.create_cylinder(radius=1.0, height=1.0, center=None, angle=None)[source]¶

Construct a cylindrical NURBS surface.

Built by extruding a circle along the z-axis.

Parameters:

radius (float) – Cylinder radius.
height (float) – Cylinder height along the z-axis.
center (ArrayLike | None) – Center of the base circle (up to 3D, zero-padded).
angle (float | tuple[float, float] | None) – Sweep specification (same as create_circle()).

Returns:

A 2D rational B-spline surface.

Return type:

Example

>>> cyl = create_cylinder(radius=2, height=5)
>>> cyl.dim
2

pantr.cad.create_disk(radius_inner=0.0, radius_outer=1.0, center=None, angle=None)[source]¶

Construct a disk or annular sector as a NURBS surface.

When radius_inner > 0, produces an annular sector via create_ruled() between inner and outer circular arcs. When radius_inner == 0, the inner boundary degenerates to a point.

Parameters:

radius_inner (float) – Inner radius. Use 0 for a full disk.
radius_outer (float) – Outer radius.
center (ArrayLike | None) – Center point (up to 3D, zero-padded). If None, centered at the origin.
angle (float | tuple[float, float] | None) – Sweep specification (same as create_circle()).

Returns:

A 2D rational B-spline surface.

Return type:

Example

>>> d = create_disk(radius_outer=2.0)
>>> d.dim
2

pantr.cad.create_extrusion(bspline, displacement)[source]¶

Extrude a B-spline curve or surface along a displacement vector.

Creates a new B-spline with one additional parametric dimension by translating the input along the given vector. The new direction is appended as the last parametric axis with degree 1 and knots [0, 0, 1, 1].

Parameters:

bspline (Bspline) – Input curve (dim=1) or surface (dim=2).
displacement (ArrayLike) – Translation vector (up to 3D, zero-padded).

Returns:

A B-spline with dim + 1 parametric dimensions.

Return type:

Raises:

ValueError – If bspline.dim > 2.

Example

>>> from pantr.cad import create_circle, create_extrusion
>>> cyl = create_extrusion(create_circle(), [0, 0, 1])
>>> cyl.dim
2

pantr.cad.create_line(p0=(0.0, 0.0, 0.0), p1=(1.0, 0.0, 0.0))[source]¶

Construct a straight-line B-spline curve between two points.

Creates a degree-1 (linear), non-rational B-spline curve with two control points. Input points shorter than 3 elements are zero-padded to 3D.

Parameters:

p0 (ArrayLike) – Start point. Defaults to the origin.
p1 (ArrayLike) – End point. Defaults to (1, 0, 0).

Returns:

A 1D, degree-1, rank-3, non-rational B-spline curve.

Return type:

Example

>>> crv = create_line([0, 0], [1, 1])
>>> crv.degree
(1,)
>>> crv.rank
3

pantr.cad.create_rectangle(corner=(0.0, 0.0, 0.0), width=1.0, height=1.0)[source]¶

Construct a closed rectangular B-spline curve in the xy-plane.

Creates a degree-1 curve with 5 control points (the first point repeated to close the loop) and 4 knot spans.

Parameters:

corner (ArrayLike) – Bottom-left corner (up to 3D, zero-padded).
width (float) – Rectangle width along the x-axis.
height (float) – Rectangle height along the y-axis.

Returns:

A 1D, degree-1, rank-3, non-rational closed curve.

Return type:

Example

>>> rect = create_rectangle([0, 0], 2, 3)
>>> rect.dim
1

pantr.cad.create_revolution(bspline, point, axis=2, angle=None)[source]¶

Revolve a B-spline curve or surface around an axis.

Creates a new B-spline with one additional parametric dimension (the angular direction, appended last). The input is first promoted to rational if needed, then transformed to a coordinate system aligned with the rotation axis, revolved via the circular arc construction, and transformed back.

The angular direction inherits the same span/continuity structure as create_circle(): one span per 90 degrees, C0 at arc junctions.

Parameters:

bspline (Bspline) – Input curve (dim=1) or surface (dim=2) with rank 3.
point (ArrayLike) – A point on the rotation axis (up to 3D, zero-padded).
axis (int | ArrayLike) – Rotation axis. An int in {0, 1, 2} selects a coordinate axis. An array-like of length 3 specifies an arbitrary axis direction (normalised internally).
angle (float | tuple[float, float] | None) – Sweep specification (same as create_circle()).

Returns:

A rational B-spline with dim + 1 dimensions.

Return type:

Raises:

ValueError – If bspline.dim > 2.
ValueError – If bspline.rank != 3.

pantr.cad.create_ruled(bspline1, bspline2)[source]¶

Construct a ruled surface or volume between two B-splines.

Creates a new B-spline by linearly interpolating control points between bspline1 (at parameter 0) and bspline2 (at parameter 1) along a new last parametric axis with degree 1.

The two inputs are first made compatible via make_compat() so they share the same degree and knot vectors. If one is rational and the other is not, the non-rational one is promoted.

Parameters:

bspline1 (Bspline) – First boundary (curve or surface, dim <= 2).
bspline2 (Bspline) – Second boundary (same dim as bspline1).

Returns:

A B-spline with dim + 1 parametric dimensions.

Return type:

Raises:

ValueError – If the inputs have different parametric dimensions.
ValueError – If either input has dim > 2.

Example

>>> from pantr.cad import create_circle, create_ruled
>>> annulus = create_ruled(create_circle(radius=0.5), create_circle(radius=1.0))
>>> annulus.dim
2

pantr.cad.create_sweep(section, trajectory)[source]¶

Construct the translational sweep of a section along a trajectory.

Creates a new B-spline by summing the section and trajectory geometries: S(u, v) = section(u) + trajectory(v). The trajectory direction is appended as the last parametric axis.

For rational B-splines the product formula applies: the result weight is the product of section and trajectory weights, and the weighted coordinates combine as w_s * C_t + w_t * C_s.

Parameters:

section (Bspline) – Section curve (dim=1) or surface (dim=2).
trajectory (Bspline) – Trajectory curve (dim=1).

Returns:

A B-spline with section.dim + 1 dimensions.

Return type:

Raises:

ValueError – If section.dim > 2.
ValueError – If trajectory.dim != 1.

pantr.cad.create_trilinear(corners=None)[source]¶

Construct a trilinear B-spline volume from eight corner points.

Creates a degree (1, 1, 1), non-rational B-spline volume with 2 x 2 x 2 control points.

The corner ordering follows the tensor-product convention:

   corners[0,1,1]       corners[1,1,1]
   o--------------------o
  /|                   /|
 / |                  / |                 w
o--------------------o  |                 ^  v
|  | corners[0,0,1]  |  | corners[1,0,1]  | /
|  |                 |  |                 |/
|  o-----------------|--o                 +------> u
| / corners[0,1,0]   | / corners[1,1,0]
|/                   |/
o--------------------o
corners[0,0,0]       corners[1,0,0]

Parameters:: corners (ArrayLike | None) – Array of shape (2, 2, 2, rank) with rank <= 3. Coordinates shorter than 3D are zero-padded. If None, defaults to the unit cube [-0.5, 0.5]^3 centered at the origin.
Returns:: A 3D, degree-(1, 1, 1), rank-3, non-rational B-spline volume.
Return type:: Bspline
Raises:: ValueError – If corners does not have shape (2, 2, 2, rank) with 1 <= rank <= 3.

Example

>>> vol = create_trilinear()
>>> vol.degree
(1, 1, 1)
>>> vol.dim
3

pantr.cad.join(bspline1, bspline2, axis)[source]¶

Join two B-splines along a parametric axis with C0 continuity.

The two inputs are made compatible on all non-join axes (degree and knots), elevated to the same degree on the join axis, and then concatenated. The shared boundary control points are averaged.

Optionally, knots are removed at the junction to recover higher smoothness when the geometry permits it.

Parameters:

bspline1 (Bspline) – First B-spline (left side of the join).
bspline2 (Bspline) – Second B-spline (right side of the join).
axis (int) – Parametric axis along which to join (0-indexed).

Returns:

A single B-spline spanning both inputs.

Return type:

Raises:

ValueError – If the inputs have different parametric dimensions.
ValueError – If axis is out of range.

pantr.cad.make_compat(*bsplines, axes=None)[source]¶

Make B-splines compatible along specified parametric axes.

Returns new B-splines that share the same polynomial degree and knot vector along the given axes. The geometric mapping of each B-spline is preserved exactly.

The algorithm proceeds in three stages:

Same domain – remap knot vectors to the common envelope [min(starts), max(ends)] per axis.
Same degree – elevate each B-spline to the maximum degree per axis.
Merge knots – insert knots so all B-splines share the same interior breakpoints with the maximum multiplicity.

Periodic B-splines are converted to open form before processing.

Parameters:

*bsplines (Bspline) – One or more B-splines to make compatible. All must have the same parametric dimension. If fewer than two are provided, returns them unchanged.
axes (int | Sequence[int] | None) – Parametric axes along which to operate. None (default) means all axes. A single int or a sequence of int selects specific axes.

Returns:

New B-splines with identical knot structure along the specified axes.

Return type:

list[Bspline]

Raises:

ValueError – If B-splines have different parametric dimensions.
ValueError – If any axis index is out of range.

Affine transformations for geometric objects.

Provides the AffineTransform class, which represents an affine map T(x) = A @ x + b in n-dimensional space, together with factory methods for common transformations (translation, rotation, scaling, mirroring, shear) and composition operators.

Main exports:

AffineTransform — logically immutable affine-transformation object.

class pantr.transform.AffineTransform(matrix, translation=None)[source]¶

Bases: object

An affine transformation T(x) = A x + b in n-dimensional space.

The transformation is defined by a square matrix A (the linear part) and a translation vector b. Instances are treated as immutable: every factory method and operator returns a new AffineTransform; no method mutates an existing instance.

Parameters:

matrix (npt.ArrayLike)
translation (npt.ArrayLike | None)

_matrix¶

The (n, n) linear part.

Type:: npt.NDArray[np.float64]

_translation¶

The (n,) translation.

Type:: npt.NDArray[np.float64]

__init__(matrix, translation=None)[source]¶

Create an affine transformation from a matrix and translation.

Parameters:

matrix (npt.ArrayLike) – The (n, n) linear part of the transformation. Must be a square 2-D array.
translation (npt.ArrayLike | None) – The (n,) translation vector. If None, defaults to the zero vector.

Raises:

ValueError – If matrix is not 2-D or not square.
ValueError – If translation length does not match the matrix dimension.

Return type:

None

Note

Both matrix and translation are stored as C-contiguous, read-only float64 arrays.

property dim: int¶

Get the spatial dimension of the transformation.

Returns:: Dimension n of the transformation.
Return type:: int

property matrix: ndarray[tuple[Any, ...], dtype[float64]]¶

Get the linear part of the transformation.

Returns:: Read-only (n, n) matrix.
Return type:: npt.NDArray[np.float64]

property offset: ndarray[tuple[Any, ...], dtype[float64]]¶

Get the translation (offset) part of the transformation.

Returns:: Read-only (n,) vector.
Return type:: npt.NDArray[np.float64]

property inverse: AffineTransform¶

Get the inverse transformation.

Computed once and cached; subsequent accesses are free. Safe because neither the matrix nor the translation is ever modified after construction.

Returns:: The inverse such that T @ T.inverse is the identity.
Return type:: AffineTransform
Raises:: ValueError – If the matrix is singular.

static identity(n)[source]¶

Create the identity transformation in n dimensions.

Parameters:: n (int) – Spatial dimension.
Returns:: The identity map.
Return type:: AffineTransform

static translation(offset)[source]¶

Create a pure translation.

Parameters:: offset (npt.ArrayLike) – Translation vector of length n.
Returns:: A translation by offset.
Return type:: AffineTransform

static scaling(factors, *, center=None)[source]¶

Create a scaling transformation.

Parameters:

factors (float | npt.ArrayLike) – If a scalar, isotropic scaling is applied and center (or a separate call) determines the dimension. If an array, anisotropic scaling along each axis.
center (npt.ArrayLike | None) – Optional center point. If given, the scaling is performed about this point rather than the origin.

Returns:

The scaling transformation.

Return type:

Raises:

ValueError – If factors is a scalar and center is None (dimension cannot be inferred).
ValueError – If factors is a scalar and center is not a 1-D array-like.
ValueError – If any factor is non-finite or zero (singular transform).
ValueError – If factors is an array and center has the wrong shape.

static rotation_2d(angle, *, center=None)[source]¶

Create a 2-D counter-clockwise rotation.

Parameters:

angle (float) – Rotation angle in radians.
center (npt.ArrayLike | None) – Optional center of rotation.

Returns:

The 2-D rotation.

Return type:

Raises:

ValueError – If angle is non-finite.

static rotation_3d(angle, axis=2, *, center=None)[source]¶

Create a 3-D rotation via the Rodrigues formula.

Parameters:

angle (float) – Rotation angle in radians.
axis (int | npt.ArrayLike) – Rotation axis. An int in {0, 1, 2} selects the corresponding coordinate axis (x, y, z). An array-like of length 3 specifies an arbitrary axis (will be normalised internally).
center (npt.ArrayLike | None) – Optional center of rotation.

Returns:

The 3-D rotation.

Return type:

Raises:

ValueError – If angle is non-finite.
ValueError – If an integer axis is not in {0, 1, 2}.
ValueError – If a vector axis does not have shape (3,).
ValueError – If a vector axis is zero or non-finite.

static mirror(normal, *, center=None)[source]¶

Create a reflection (mirror) across a hyperplane.

The hyperplane passes through the origin (or center) and has the given normal vector. The Householder formula is used: A = I - 2 n nᵀ.

Parameters:

normal (npt.ArrayLike) – Normal vector of the mirror plane. Will be normalised internally.
center (npt.ArrayLike | None) – Optional point on the mirror plane.

Returns:

The reflection.

Return type:

Raises:

ValueError – If normal is zero or non-finite.

static shear(dim, component, direction, factor)[source]¶

Create a shear transformation.

The resulting map adds factor * x[direction] to x[component], leaving all other components unchanged.

Parameters:

dim (int) – Spatial dimension.
component (int) – The axis that is modified.
direction (int) – The axis whose value drives the shear.
factor (float) – Shear magnitude.

Returns:

The shear transformation.

Return type:

Raises:

ValueError – If component equals direction.
ValueError – If component or direction is out of range.
ValueError – If factor is non-finite.

compose(other)[source]¶

Compose this transformation with other.

Returns the transformation self(other(x)).

Parameters:: other (AffineTransform) – The inner transformation.
Returns:: The composed transformation.
Return type:: AffineTransform
Raises:: ValueError – If the dimensions do not match.

Axis-aligned bounding boxes and geometry primitives.

This module exposes AABB, a lightweight, immutable value type for an axis-aligned bounding box in any spatial dimension ndim >= 1. It is the shared box / domain primitive for PaNTr (for example, the parametric domain of a spline space or the bounds of a grid cell) and for libraries built on PaNTr.

An AABB stores two read-only float64 corner arrays AABB.lo and AABB.hi of shape (ndim,). Entries may be finite or +/- numpy.inf (for unbounded axes); numpy.nan is rejected. A point x is inside the box when lo[i] <= x[i] <= hi[i] on every axis; a box with lo[i] > hi[i] on some axis is empty and reported by AABB.is_empty().

Main exports:

AABB – immutable axis-aligned bounding box.

class pantr.geometry.AABB(lo, hi)[source]¶

Bases: object

Axis-aligned bounding box in any spatial dimension ndim >= 1.

An AABB stores two 1-D arrays lo and hi of equal shape (ndim,). Entries may be finite or +/- numpy.inf (for unbounded axes); numpy.nan is rejected. Instances are frozen: the stored arrays are C-contiguous, float64, and flagged read-only, and __setattr__ rejects attempts to replace the bound attributes.

The sign convention is the natural one: a point x is “inside” the box when lo[i] <= x[i] <= hi[i] on every axis. A box with lo[i] > hi[i] on some axis is empty (no point satisfies the inequality); is_empty() reports this.

Parameters:

lo (ndarray[tuple[Any, ...], dtype[float64]])
hi (ndarray[tuple[Any, ...], dtype[float64]])

lo: ndarray[tuple[Any, ...], dtype[float64]]¶: Lower corner, shape (ndim,), read-only.

hi: ndarray[tuple[Any, ...], dtype[float64]]¶: Upper corner, shape (ndim,), read-only.

__init__(lo, hi)[source]¶

Build and validate an AABB from array-like corners.

The arguments are coerced to C-contiguous float64 arrays and checked for NaN. The resulting buffers are cached as read-only attributes on self.

Parameters:

lo (npt.ArrayLike) – Lower corner; array-like of finite-or-infinite floats. Any rank is accepted and ravelled to 1-D of length ndim.
hi (npt.ArrayLike) – Upper corner, same shape and semantics as lo.

Raises:

ValueError – If the shapes mismatch, contain a NaN, or the implied ndim is less than 1.
TypeError – If lo or hi has a non-numeric dtype.

Return type:

None

property ndim: int¶

Get the spatial dimensionality of the box.

Returns:: The number of axes (>= 1).
Return type:: int

is_empty()[source]¶

Check whether the box contains no points.

A box is empty iff lo[i] > hi[i] on at least one axis. Useful when constructing boxes manually, accepting boxes from external sources, or after pad() with a negative radius.

Returns:: True when the box is empty.
Return type:: bool

contains_point(x)[source]¶

Check whether point x lies inside or on the boundary of the box.

A point is inside iff lo[i] <= x[i] <= hi[i] on every axis. An empty box contains no points.

Parameters:

x (npt.ArrayLike) – Point to test; array-like of floats, ravelled to 1-D of length ndim.

Returns:

True when x is inside or on the boundary.

Return type:

bool

Raises:

ValueError – If the ravelled x does not have length ndim, or if x contains NaN.
TypeError – If x has a non-numeric dtype.

overlaps(other)[source]¶

Check whether self and other share at least one point.

Empty boxes (on either side) overlap nothing.

Parameters:: other (AABB) – The box to test against; must share ndim.
Returns:: True when the two boxes intersect, False otherwise.
Return type:: bool
Raises:: ValueError – If other.ndim != self.ndim.

union(other)[source]¶

Return the smallest axis-aligned box that contains self and other.

Empty boxes act as neutral elements: union(empty, x) == x.

Parameters:: other (AABB) – Box to union with; must share ndim.
Returns:: The union bounding box.
Return type:: AABB
Raises:: ValueError – If other.ndim != self.ndim.

intersect(other)[source]¶

Return the axis-aligned intersection, or None if disjoint.

Parameters:: other (AABB) – Box to intersect with; must share ndim.
Returns:: The intersection, or None when the two boxes do not overlap (including the case where either operand is empty).
Return type:: AABB | None
Raises:: ValueError – If other.ndim != self.ndim.

pad(r)[source]¶

Inflate the box by r on every axis (symmetric on both sides).

A scalar r expands every axis by the same amount; a length-ndim array expands each axis by its own entry. Padding an unbounded axis leaves it unbounded. Padding by a negative r shrinks the box and can produce an empty AABB, which is allowed and reported by is_empty().

Parameters:: r (float | npt.ArrayLike) – Per-side padding amount. Scalar or length-ndim array-like of finite reals.
Returns:: The padded box.
Return type:: AABB
Raises:: ValueError – If r has the wrong shape or non-finite entries.

transform(affine)[source]¶

Return the axis-aligned wrap of the image of self under affine.

For an affine map T(x) = A x + b, the tight axis-aligned box containing T(self) is computed from the per-entry sign of A: each output axis i receives contributions A[i, j] * lo[j] or A[i, j] * hi[j] – whichever gives the minimum / maximum – summed over j. Zero entries of A contribute nothing even when lo[j] / hi[j] is infinite; this preserves the correct wrap for unbounded axes that the transform projects out.

Parameters:: affine (_AffineMap) – The affine map to apply – any object exposing dim, matrix (A), and offset (b). Its dimension must match ndim.
Returns:: The axis-aligned wrap of the transformed box.
Return type:: AABB
Raises:: ValueError – If affine.dim != self.ndim, if affine.matrix is not square (ndim, ndim), or if the wrap produces NaN (e.g. a matrix containing NaN, or inf - inf arithmetic from opposing infinite bounds in the same row).

static unbounded(ndim)[source]¶

Build the everywhere-true AABB (-inf, +inf)^ndim.

Parameters:: ndim (int) – Spatial dimension (>= 1).
Returns:: The unbounded AABB.
Return type:: AABB
Raises:: ValueError – If ndim < 1.

static empty(ndim)[source]¶

Build an empty (zero-volume) AABB with lo > hi.

An empty AABB contains no points (detected by is_empty()) and acts as the neutral element of union() (union(empty, x) == x). Symmetric counterpart to unbounded().

Parameters:: ndim (int) – Spatial dimension (>= 1).
Returns:: An empty AABB (lo = +inf, hi = -inf).
Return type:: AABB
Raises:: ValueError – If ndim < 1.

static from_bounds(bounds)[source]¶

Build an AABB from a (ndim, 2) [[lo, hi], ...] array.

The dual of as_bounds(); useful when interoperating with numpy.histogramdd-style bounds arguments.

Parameters:

bounds (npt.ArrayLike) – (ndim, 2) array-like of [lo, hi] rows.

Returns:

The constructed AABB.

Return type:

AABB

Raises:

ValueError – If bounds does not have shape (ndim, 2) with ndim >= 1.
TypeError – If bounds has a non-numeric dtype.

as_bounds()[source]¶

Return the AABB as a (ndim, 2) [[lo, hi], ...] array.

Returns:: A freshly-allocated, writeable, C-contiguous (ndim, 2) array.
Return type:: npt.NDArray[np.float64]

Grids & quadrature¶

Structured cell grids for the PaNTr stack.

This package provides a small, performance-conscious grid layer: a partition of a parametric domain into cells with implicit (computed, not stored) connectivity. It is the shared grid abstraction consumed by immersed / unfitted discretizations, B-spline knot-span grids, and hierarchical refinement grids.

The TensorProductGrid is deliberately low-footprint: it stores only the per-axis breakpoint arrays and a little metadata, computing cell bounds, neighbours, and ids on demand. The only O(num_cells) structure – a BVH spatial index – is built lazily on the first query_aabb(), so a grid used purely as a B-spline’s knot grid stays proportional to its breakpoints, not its cell count.

Main exports:

Grid: abstract base class defining the grid contract and supplying axis-aligned box defaults for facets, neighbours, reference maps, batch point location, and spatial queries.
GridRestriction: result of Grid.restrict() – a windowed sub-grid plus local-to-global cell index maps.
TensorProductGrid: concrete tensor-product grid of axis-aligned boxes with per-axis breakpoints and row-major (C-order) cell ids.
HierarchicalGrid: hierarchical grid with a fixed per-direction subdivision factor; active cells stored as rectangular blocks per level.
uniform_grid(): build a uniform grid on a bounding box.
tensor_product_grid(): build the knot-span grid of a pantr.bspline.BsplineSpace.
hierarchical_grid(): build a HierarchicalGrid from a root TensorProductGrid and a subdivision factor.
BVH: bounding-volume hierarchy over cell AABBs, backing Grid.query_aabb().
CellTags, FacetTags: sparse, dolfinx-style named tag registries for cells and facets.
cell_quadrature(): map a pantr.quad.QuadratureRule from the unit cube onto a grid’s cells (per-cell points and weights).
overlay(): the coarsest TensorProductGrid refining two input tensor-product grids (union of per-axis breakpoints on their domain overlap).
Partition: a per-cell owner assignment for distributing a grid’s cells.
partition_grid(): split a grid into n_parts rank subdomains (native, dependency-free): the "block" Cartesian backend, the "rcb" recursive coordinate bisection backend (weight- and activity-aware), and an "auto" dispatch between them.

class pantr.grid.BVH(node_lo, node_hi, node_left, node_right, node_cell, *, n_cells)[source]¶

Bases: object

Bounding-volume hierarchy indexing a fixed collection of AABBs.

Instances are immutable: once built, the internal arrays are flagged read-only. Queries allocate fresh int64 output arrays per call.

Build by passing per-cell AABBs to from_cell_bounds(); direct construction from the raw array representation is supported via the default constructor but is mostly intended for tests and round-trip serialization.

Parameters:

node_lo (npt.NDArray[np.float64])
node_hi (npt.NDArray[np.float64])
node_left (npt.NDArray[np.int64])
node_right (npt.NDArray[np.int64])
node_cell (npt.NDArray[np.int64])
n_cells (int)

__init__(node_lo, node_hi, node_left, node_right, node_cell, *, n_cells)[source]¶

Store the raw BVH arrays after validating their shapes.

Callers should prefer from_cell_bounds(); this constructor is useful for tests that need to poke specific tree shapes.

Parameters:

node_lo (npt.NDArray[np.float64]) – Per-node AABB lo corners, shape (n_nodes, ndim).
node_hi (npt.NDArray[np.float64]) – Per-node AABB hi corners, shape (n_nodes, ndim).
node_left (npt.NDArray[np.int64]) – Left-child indices; -1 on leaves. Shape (n_nodes,).
node_right (npt.NDArray[np.int64]) – Right-child indices; -1 on leaves.
node_cell (npt.NDArray[np.int64]) – Leaf cell identifiers; -1 on internal nodes.
n_cells (int) – Number of indexed cells (leaves).

Raises:

TypeError – If any array has the wrong dtype.
ValueError – If shapes are inconsistent, ndim is < 1, or n_nodes != 2 * n_cells - 1 (0 when n_cells == 0).

Return type:

None

ndim: int¶: Spatial dimension of the indexed AABBs (>= 1).

n_cells: int¶: Number of cells indexed (equal to the number of leaves).

n_nodes: int¶: Total number of nodes (2 * n_cells - 1 for n_cells > 0, else 0).

classmethod from_cell_bounds(cell_lo, cell_hi)[source]¶

Build a BVH over n_cells axis-aligned cell AABBs.

Uses a top-down recursive median-of-longest-axis split. Cells are sorted by centroid on the longest axis; the median splits the list into two halves of equal size (+/- 1). Each leaf indexes exactly one cell.

Parameters:

cell_lo (npt.ArrayLike) – Per-cell lo corners; shape (n_cells, ndim) with ndim >= 1. Validated, not mutated.
cell_hi (npt.ArrayLike) – Per-cell hi corners; same shape and conventions as cell_lo. Each entry must satisfy cell_hi >= cell_lo.

Returns:

The constructed hierarchy.

Return type:

Self

Raises:

TypeError – If inputs cannot be cast to float64.
ValueError – If shapes are inconsistent, ndim is < 1, any cell has hi < lo, or the implied tree exceeds the internal stack depth.

query_aabb(aabb)[source]¶

Return the ids of every leaf cell whose AABB overlaps aabb.

Parameters:: aabb (AABB) – Query box; must match ndim.
Returns:: Overlapping cell ids. Order matches the internal preorder traversal; callers that need a particular order should sort the result.
Return type:: npt.NDArray[np.int64]
Raises:: ValueError – If aabb.ndim != self.ndim.

property node_lo: npt.NDArray[np.float64]¶

Get the read-only view of per-node AABB lo corners.

Returns:: Shape (n_nodes, ndim).
Return type:: npt.NDArray[np.float64]

property node_hi: npt.NDArray[np.float64]¶

Get the read-only view of per-node AABB hi corners.

Returns:: Shape (n_nodes, ndim).
Return type:: npt.NDArray[np.float64]

property node_left: npt.NDArray[np.int64]¶

Get the read-only view of per-node left-child indices.

Returns:: Shape (n_nodes,); -1 on leaves.
Return type:: npt.NDArray[np.int64]

property node_right: npt.NDArray[np.int64]¶

Get the read-only view of per-node right-child indices.

Returns:: Shape (n_nodes,); -1 on leaves.
Return type:: npt.NDArray[np.int64]

property node_cell: npt.NDArray[np.int64]¶

Get the read-only view of per-leaf cell identifiers.

Returns:: Shape (n_nodes,); -1 on internal nodes, 0 <= id < n_cells on leaves.
Return type:: npt.NDArray[np.int64]

class pantr.grid.CellTags(num_cells)[source]¶

Bases: object

Sparse named integer tags over a grid’s cells.

Each tag named name is a pair of parallel int64 arrays (ids, values) sorted by ids; a cell not listed in ids is untagged under name. Distinct tag names are independent. The owning grid’s cell count is exposed through the num_cells property.

Parameters:: num_cells (int)

__init__(num_cells)[source]¶

Create an empty cell-tag registry.

Parameters:: num_cells (int) – Number of cells in the owning grid (>= 0).
Raises:: ValueError – If num_cells is negative.
Return type:: None

property num_cells: int¶

Get the number of cells in the owning grid.

Returns:: The cell count; valid cell ids are [0, num_cells).
Return type:: int

set(name, ids, values)[source]¶

Create or replace the tag name with the association ids -> values.

Parameters:

name (str) – Tag name. Replaces any existing tag with the same name.
ids (npt.ArrayLike) – 1-D integer array-like of cell ids; each must satisfy 0 <= id < num_cells and be unique.
values (npt.ArrayLike) – Scalar integer (broadcast to every id) or a 1-D integer array-like of the same length as ids.

Raises:

TypeError – If ids or values is not integer-typed.
ValueError – If ids is not 1-D, contains duplicates, or has an id out of range, or if values has a length other than len(ids).

Return type:

None

property names: tuple[str, ...]¶

Get the registered tag names.

Returns:: Tag names in insertion order.
Return type:: tuple[str, …]

remove(name)[source]¶

Delete the tag name.

Parameters:: name (str) – Tag name.
Raises:: KeyError – If no tag named name exists.
Return type:: None

to_dense(name, *, fill=0, dtype=<class 'numpy.int64'>)[source]¶

Scatter tag name into a dense (num_cells,) array.

Untagged cells receive fill. Useful when a downstream Numba kernel wants a per-cell label array rather than the sparse representation.

Parameters:

name (str) – Tag name.
fill (int) – Value for untagged cells. Defaults to 0.
dtype (npt.DTypeLike) – Output integer dtype. Defaults to numpy.int64. Values are stored as int64 internally; if dtype is narrower than int64 and any stored value falls outside the dtype’s representable range, an OverflowError is raised rather than silently truncating the value.

Returns:

Fresh, writeable (num_cells,) array with dtype as the scalar type.

Return type:

npt.NDArray[Any]

Raises:

KeyError – If no tag named name exists.
TypeError – If dtype is not an integer or unsigned-integer dtype.
OverflowError – If any stored value cannot be represented exactly in dtype (only raised when dtype is narrower than int64).

class pantr.grid.FacetTags(num_cells, facets_per_cell)[source]¶

Bases: object

Sparse named integer tags over a grid’s local facets.

Each facet is addressed by a (cell_id, local_facet_id) key, with local_facet_id in [0, facets_per_cell). Each tag named name is a pair (keys, values) where keys is an (M, 2) int64 array of (cell_id, local_facet_id) rows and values is a length-M int64 array, sorted lexicographically by key. The owning grid’s cell count and per-cell facet count are exposed through the num_cells and facets_per_cell properties.

Parameters:

num_cells (int)
facets_per_cell (int)

__init__(num_cells, facets_per_cell)[source]¶

Create an empty facet-tag registry.

Parameters:

num_cells (int) – Number of cells in the owning grid (>= 0).
facets_per_cell (int) – Number of local facets per cell (>= 1).

Raises:

ValueError – If num_cells is negative or facets_per_cell is < 1.

Return type:

None

property num_cells: int¶

Get the number of cells in the owning grid.

Returns:: The cell count.
Return type:: int

property facets_per_cell: int¶

Get the number of local facets per cell.

Returns:: 2 * ndim for an axis-aligned box grid.
Return type:: int

set(name, keys, values)[source]¶

Create or replace the tag name with the association keys -> values.

Parameters:

name (str) – Tag name. Replaces any existing tag with the same name.
keys (npt.ArrayLike) – (M, 2) integer array-like of (cell_id, local_facet_id) rows; each row must be unique with 0 <= cell_id < num_cells and 0 <= local_facet_id < facets_per_cell.
values (npt.ArrayLike) – Scalar integer (broadcast to every key) or a 1-D integer array-like of length M.

Raises:

TypeError – If keys or values is not integer-typed.
ValueError – If keys does not have shape (M, 2), contains a duplicate or out-of-range key, or values has a length other than M.

Return type:

None

property names: tuple[str, ...]¶

Get the registered tag names.

Returns:: Tag names in insertion order.
Return type:: tuple[str, …]

remove(name)[source]¶

Delete the tag name.

Parameters:: name (str) – Tag name.
Raises:: KeyError – If no tag named name exists.
Return type:: None

to_dense(name, *, fill=0, dtype=<class 'numpy.int64'>)[source]¶

Scatter tag name into a dense (num_cells, facets_per_cell) array.

Untagged facets receive fill. Useful when a downstream Numba kernel wants a per-facet label array rather than the sparse representation.

Parameters:

name (str) – Tag name.
fill (int) – Value for untagged facets. Defaults to 0.
dtype (npt.DTypeLike) – Output integer dtype. Defaults to numpy.int64. Values are stored as int64 internally; if dtype is narrower than int64 and any stored value falls outside the dtype’s representable range, an OverflowError is raised rather than silently truncating the value.

Returns:

Fresh, writeable (num_cells, facets_per_cell) array with dtype as the scalar type.

Return type:

npt.NDArray[Any]

Raises:

KeyError – If no tag named name exists.
TypeError – If dtype is not an integer or unsigned-integer dtype.
OverflowError – If any stored value cannot be represented exactly in dtype (only raised when dtype is narrower than int64).

class pantr.grid.Grid[source]¶

Bases: ABC

Abstract structured cell grid with implicit connectivity.

See the module docstring for the contract and the set of methods a subclass must implement versus those provided as overridable defaults. The size metadata is exposed through the ndim and num_cells properties.

__init__()[source]¶

Initialize the lazy spatial-index and tag-registry slots.

Subclasses must call super().__init__() before use so the lazy cell_tags, facet_tags, and query_aabb() caches are available.

Return type:: None

abstract property ndim: int¶

Get the spatial dimension of the grid.

Returns:: Number of axes (>= 1).
Return type:: int

abstract property num_cells: int¶

Get the number of cells in this (local) grid.

Returns:: Non-negative cell count.
Return type:: int

abstract cell_bounds(cid)[source]¶

Return the axis-aligned (lo, hi) corners of cell cid.

Parameters:: cid (int) – Cell identifier; must satisfy 0 <= cid < num_cells.
Returns:: Fresh, writeable length-ndim float64 arrays.
Return type:: tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]
Raises:: IndexError – If cid is out of range.

abstract locate(pt)[source]¶

Return the cell containing pt, or None if pt is outside.

Parameters:: pt (npt.ArrayLike) – Length-ndim point in parametric coordinates.
Returns:: Containing cell id, or None when pt lies outside every cell.
Return type:: int | None
Raises:: ValueError – If pt does not have length ndim.

abstract neighbor_across_facet(cid, lfid)[source]¶

Return the cell across local facet lfid of cid, or None.

The result is None iff the facet lies on the grid’s outer boundary.

Parameters:

cid (int) – Cell identifier.
lfid (int) – Local facet identifier in [0, num_local_facets(cid)).

Returns:

Neighbouring cell id, or None on a boundary facet.

Return type:

int | None

Raises:

IndexError – If cid or lfid is out of range.

iter_cells()[source]¶

Yield every cell identifier exactly once, in id order.

Returns:: iter(range(num_cells)).
Return type:: Iterator[int]

cell_aabb(cid)[source]¶

Return cell cid’s axis-aligned bounding box.

Parameters:: cid (int) – Cell identifier.
Returns:: Equivalent to AABB(*self.cell_bounds(cid)).
Return type:: AABB
Raises:: IndexError – If cid is out of range.

cell_level(cid)[source]¶

Return the refinement level of cell cid.

Non-hierarchical grids always return 0; hierarchical backends use level >= 1 for refined cells.

Parameters:: cid (int) – Cell identifier.
Returns:: Refinement level (0 for a flat grid).
Return type:: int
Raises:: IndexError – If cid is out of range.

child_cells(cid)[source]¶

Return the immediate refinement children of cell cid.

For a flat (non-hierarchical) grid this is always empty.

Parameters:: cid (int) – Cell identifier.
Returns:: Child cell ids; empty for a flat grid.
Return type:: tuple[int, …]
Raises:: IndexError – If cid is out of range.

reference_map(cid)[source]¶

Return the affine map [0, 1]^ndim -> cell for cell cid.

For an axis-aligned cell this is T(u) = diag(hi - lo) @ u + lo.

Parameters:: cid (int) – Cell identifier.
Returns:: Push-forward from the unit cube to the cell.
Return type:: AffineTransform
Raises:: IndexError – If cid is out of range.

neighbors(cid)[source]¶

Return the facet-neighbour cell ids of cid.

Built by collecting neighbor_across_facet() over every local facet and dropping boundary facets (None).

Parameters:: cid (int) – Cell identifier.
Returns:: Neighbouring cell ids; length between ndim (corner cell) and 2 * ndim (interior cell) for a box grid.
Return type:: list[int]
Raises:: IndexError – If cid is out of range.

restrict(cell_ids)[source]¶

Return the structured sub-grid spanning a subset of cells.

The sub-grid is the smallest grid of the same kind that contains every requested cell – for a tensor-product grid, the multi-index bounding box of cell_ids. Its breakpoints are pure slices of this grid’s (never re-clamped or re-based), so a B-spline or field built on the sub-grid agrees pointwise with the global one over the shared cells.

The bounding box may contain cells that were not requested (when cell_ids is non-convex); GridRestriction.in_subset flags, per local cell, whether it was requested (True) or is bounding-box fill (False).

Restriction is an optional grid capability: this base implementation raises NotImplementedError. Grids with structured windowing (for example pantr.grid.TensorProductGrid) override it.

Parameters:

cell_ids (npt.ArrayLike) – Flat cell identifiers to span; duplicates are ignored. Each must satisfy 0 <= cid < num_cells.

Returns:

The windowed sub-grid, the local_to_global_cell map, and the in_subset mask.

Return type:

GridRestriction

Raises:

NotImplementedError – If this grid kind does not support restriction.
ValueError – If cell_ids is empty (in overriding implementations).
IndexError – If any cell id is out of range (in overriding implementations).
TypeError – If cell_ids is not integer-valued (in overriding implementations).

num_local_facets(cid)[source]¶

Return the number of local facets of cell cid.

Defaults to 2 * ndim (an axis-aligned box).

Parameters:: cid (int) – Cell identifier.
Returns:: 2 * ndim.
Return type:: int
Raises:: IndexError – If cid is out of range.

local_facet_axis_side(cid, lfid)[source]¶

Return (axis, side) for local facet lfid of cell cid.

Uses the conventional lfid = 2 * axis + side encoding, with axis in [0, ndim) and side in {0, 1} (0 = low face, 1 = high face).

Parameters:

cid (int) – Cell identifier.
lfid (int) – Local facet identifier in [0, 2 * ndim).

Returns:

(axis, side).

Return type:

tuple[int, int]

Raises:

IndexError – If cid or lfid is out of range.

local_facet_bounds(cid, lfid)[source]¶

Return the degenerate (lo, hi) AABB of facet lfid of cell cid.

On the facet’s normal axis both corners coincide (the cell’s lo corner for side == 0, hi corner for side == 1); the other axes span the cell’s extent.

Parameters:

cid (int) – Cell identifier.
lfid (int) – Local facet identifier in [0, 2 * ndim).

Returns:

Fresh, writeable length-ndim float64 arrays.

Return type:

tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]

Raises:

IndexError – If cid or lfid is out of range.

is_mesh_boundary_facet(cid, lfid)[source]¶

Return whether facet lfid of cell cid is on the grid’s outer boundary.

Defaults to neighbor_across_facet(cid, lfid) is None.

Parameters:

cid (int) – Cell identifier.
lfid (int) – Local facet identifier in [0, num_local_facets(cid)).

Returns:

True iff no neighbouring cell shares the facet.

Return type:

bool

Raises:

IndexError – If cid or lfid is out of range.

hanging_neighbors(cid, lfid)[source]¶

Return all active cells sharing facet lfid of cid.

For conforming grids (flat TensorProductGrid) this is equivalent to neighbor_across_facet() wrapped in a tuple. For hierarchical grids, a single coarse face may abut multiple fine cells (hanging nodes); this method returns all of them.

Subclasses with hanging-node support should override this method.

Parameters:

cid (int) – Cell identifier.
lfid (int) – Local facet identifier in [0, num_local_facets(cid)).

Returns:

All neighbouring cell ids across the facet; empty when the facet lies on the grid’s outer boundary.

Return type:

tuple[int, …]

Raises:

IndexError – If cid or lfid is out of range.

locate_many(points)[source]¶

Locate a batch of points, returning one cell id per point.

Points outside the grid map to -1. The default loops over locate(); subclasses may override with a vectorized kernel.

Parameters:: points (npt.ArrayLike) – (npts, ndim) array-like of points, or a single length-ndim point.
Returns:: Shape (npts,) cell ids; -1 for points outside the grid.
Return type:: npt.NDArray[np.int64]
Raises:: ValueError – If the trailing axis of points is not ndim.

query_aabb(aabb)[source]¶

Return the ids of every cell whose AABB overlaps aabb.

Backed by a pantr.grid.BVH over the grid’s cell AABBs, built lazily on first call and cached for the grid’s lifetime. The overlap test is inclusive on every axis, so cells touching aabb on any face, edge, or corner are included.

Parameters:: aabb (AABB) – Query box; must match ndim.
Returns:: Overlapping cell ids (unordered).
Return type:: npt.NDArray[np.int64]
Raises:: ValueError – If aabb.ndim != self.ndim.

cell_bvh()[source]¶

Return the cached pantr.grid.BVH over the grid’s cell AABBs.

Built lazily on first call from collect_cell_bounds and cached. Building the BVH materializes O(num_cells) node arrays, so it is deferred until an query_aabb() (or direct) call needs it – an untagged, un-queried grid never pays this cost.

Returns:: The grid’s spatial index over its cell AABBs.
Return type:: BVH

Warning

Not fully thread-safe. Under CPython (with the GIL), concurrent first calls may each build a valid BVH and the second write silently wins — the only cost is redundant construction. Under free-threaded Python 3.13+ (--disable-gil) the assignment is not atomic and a concurrent caller could observe a partially-written reference. Call this method once on the main thread before sharing the grid across threads.

collect_cell_bounds()[source]¶

Materialize per-cell (lo, hi) as (num_cells, ndim) arrays.

The default iterates cell_bounds() over every cell in id order. Subclasses with structure (for example, tensor-product grids) should override this with a vectorized construction.

Returns:: (cell_lo, cell_hi) of shape (num_cells, ndim).
Return type:: tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]

property cell_tags: CellTags¶

Get the grid’s sparse cell-tag registry (created lazily).

Returns:: A registry of named (cell_ids, values) associations, empty until first use.
Return type:: CellTags

property facet_tags: FacetTags¶

Get the grid’s sparse facet-tag registry (created lazily).

Returns:: A registry of named ((cell_id, local_facet_id), value) associations, empty until first use. Sized for 2 * ndim facets per cell (an axis-aligned box grid).
Return type:: FacetTags

class pantr.grid.GridRestriction(grid, local_to_global_cell, in_subset)[source]¶

Bases: NamedTuple

Result of Grid.restrict(): a windowed sub-grid plus index maps.

A typing.NamedTuple returned by Grid.restrict():

grid – the windowed sub-grid: the smallest structured grid containing every requested cell. For structured implementations (e.g. TensorProductGrid) this is the same concrete type as the grid that produced it.
local_to_global_cell – shape (grid.num_cells,) map from each sub-grid cell id to its id in the original grid, in the sub-grid’s own cell-id order. Read-only.
in_subset – shape (grid.num_cells,) boolean mask: True for cells that were in the requested cell_ids, False for bounding-box fill cells (present only when the request was non-convex). Read-only.

Parameters:

grid (Grid)
local_to_global_cell (npt.NDArray[np.int64])
in_subset (npt.NDArray[np.bool_])

grid: Grid¶: Alias for field number 0

local_to_global_cell: npt.NDArray[np.int64]¶: Alias for field number 1

in_subset: npt.NDArray[np.bool_]¶: Alias for field number 2

class pantr.grid.HierarchicalGrid(root, factor)[source]¶

Bases: Grid

Hierarchical grid with a fixed per-direction uniform subdivision factor.

Built on a root TensorProductGrid. Active cells are stored as non-overlapping rectangular blocks at each level; no per-cell data is kept. A new level is created by refine(); the grid starts with all root cells active at level 0.

Flat cell ids are assigned level-by-level (level 0 first), block-by-block within a level (sorted by lo tuple), and C-order within each block.

Parameters:

root (TensorProductGrid)
factor (int | Sequence[int])

_root¶

The level-0 root grid.

Type:: TensorProductGrid

_factor¶

Per-direction subdivision factor (>= 1).

Type:: tuple[int, …]

_blocks¶

_blocks[l] is a sorted list of non-overlapping (lo, hi) blocks of active cells at level l.

Type:: list[list[tuple[tuple[int, …], tuple[int, …]]]]

_level_base¶

_level_base[l] is the flat-id base of level l; length max_level + 2 (includes sentinel).

Type:: list[int]

_num_cells¶

Cached total active cell count.

Type:: int

_version¶

Monotonic mutation counter, incremented on every structural change (see version).

Type:: int

_packed_block_lo¶

Packed block lower bounds for the Numba kernels, shape (n_blocks_total, ndim), concatenated level by level in flat-id order (see _hier_core).

Type:: npt.NDArray[np.int64]

_packed_block_hi¶

Packed block upper bounds, same shape.

Type:: npt.NDArray[np.int64]

_packed_block_base¶

Flat cell id of each block’s first cell, shape (n_blocks_total,).

Type:: npt.NDArray[np.int64]

_packed_level_start¶

Block index range of each level, shape (max_level + 2,).

Type:: npt.NDArray[np.int64]

_root_knots_flat¶

Root per-axis breakpoints concatenated end to end (kernel descriptor).

Type:: npt.NDArray[np.float64]

_root_knot_starts¶

Per-axis start offset into _root_knots_flat, shape (ndim,).

Type:: npt.NDArray[np.int64]

__init__(root, factor)[source]¶

Create a hierarchical grid from a root and a subdivision factor.

Parameters:

root (TensorProductGrid) – The level-0 grid.
factor (int | Sequence[int]) – Per-direction subdivision factor. A scalar is broadcast to every axis. Each entry must be >= 1; a factor of 1 on an axis prevents subdivision in that direction.

Raises:

TypeError – If root is not a TensorProductGrid.
ValueError – If factor has the wrong length or any entry is < 1.

Return type:

None

property ndim: int¶

Get the spatial dimension.

Returns:: Number of axes (>= 1).
Return type:: int

property num_cells: int¶

Get the total number of active cells across all levels.

Returns:: Total active cell count.
Return type:: int

property root: TensorProductGrid¶

Get the level-0 root grid.

Returns:: The root grid used to construct this hierarchy.
Return type:: TensorProductGrid

property factor: tuple[int, ...]¶

Get the per-direction subdivision factor.

Returns:: Length-ndim tuple; factor[k] >= 1.
Return type:: tuple[int, …]

property max_level: int¶

Get the index of the deepest non-empty level.

Returns:: 0 before any refinement; increases with each refine() call that adds a new level.
Return type:: int

property version: int¶

Get the monotonic mutation counter of this grid.

Incremented on every structural change (refine(), coarsen()). Snapshot consumers (e.g. THBSplineSpace) compare it to detect any post-construction mutation – max_level and num_cells alone cannot distinguish compensating refine/coarsen pairs.

Returns:: The current mutation count (>= 1).
Return type:: int

cell_bounds(cid)[source]¶

Return the axis-aligned (lo, hi) corners of cell cid.

Parameters:: cid (int) – Flat cell identifier.
Returns:: Fresh, writeable length-ndim float64 arrays.
Return type:: tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]
Raises:: IndexError – If cid is out of range.

locate(pt)[source]¶

Return the active leaf cell containing pt, or None.

Performs a top-down traversal from the root level.

Parameters:: pt (npt.ArrayLike) – Length-ndim point in parametric coordinates.
Returns:: Active cell flat id, or None when pt is outside the grid domain.
Return type:: int | None
Raises:: ValueError – If pt does not have length ndim.

locate_many(points)[source]¶

Locate a batch of points via the Numba top-down descent kernel.

Parameters:: points (npt.ArrayLike) – (npts, ndim) array-like of points, or a single length-ndim point.
Returns:: Shape (npts,) cell ids; -1 for points outside the grid (including points with NaN or infinite coordinates).
Return type:: npt.NDArray[np.int64]
Raises:: ValueError – If the trailing axis of points is not ndim.

neighbor_across_facet(cid, lfid)[source]¶

Return the cell across local facet lfid of cid, or None.

Handles hanging-node interfaces across any level difference (no 2:1 balance is enforced by refine()): when the neighbour is coarser, the active leaf covering the position – however many levels up – is returned. When the neighbour side is finer, the first active leaf descendant touching the face (lowest C-order along the face) is returned. Use hanging_neighbors() to retrieve all fine neighbours.

Parameters:

cid (int) – Cell identifier.
lfid (int) – Local facet identifier in [0, 2 * ndim).

Returns:

Neighbouring cell id, or None on a boundary facet.

Return type:

int | None

Raises:

IndexError – If cid or lfid is out of range.

hanging_neighbors(cid, lfid)[source]¶

Return all active neighbours across facet lfid of cid.

Equivalent to neighbor_across_facet() for conforming and coarser interfaces (a single neighbour, however many levels up). For a hanging (fine-side) interface, returns all active leaves touching the face, descending as many levels as the interface requires, ordered depth-first along the face (C-order over the non-axis directions at each level).

Parameters:

cid (int) – Cell identifier.
lfid (int) – Local facet identifier in [0, 2 * ndim).

Returns:

Neighbouring cell ids; empty on a boundary facet.

Return type:

tuple[int, …]

Raises:

IndexError – If cid or lfid is out of range.

restrict(cell_ids)[source]¶

Return the root-cell-aligned bounding-box sub-grid spanning cell_ids.

The window is the multi-index bounding box, in root-cell coordinates, of the root cells containing the requested leaves (a leaf at (level, midx) lives in root cell midx[k] // factor[k] ** level). The sub-grid’s root is the matching slice of this grid’s root breakpoints (never re-clamped) and it keeps the same factor; its active leaves are the per-level intersections of this grid’s blocks with the window.

Because the window is root-cell-aligned, restricting a single deep leaf returns the whole leaf-tiling of its root cell, with only the requested leaf flagged in GridRestriction.in_subset.

Parameters:

cell_ids (npt.ArrayLike) – Flat cell identifiers to span; duplicates are ignored. Each must satisfy 0 <= cid < num_cells.

Returns:

The windowed HierarchicalGrid, its local_to_global_cell map of shape (sub.num_cells,), and the boolean in_subset mask flagging requested versus bounding-box-fill cells.

Return type:

GridRestriction

Raises:

ValueError – If cell_ids is empty.
IndexError – If any cell id is out of range [0, num_cells).
TypeError – If cell_ids is not integer-valued.
RuntimeError – If an internal invariant is violated (should be unreachable).

cell_level(cid)[source]¶

Return the refinement level of cell cid.

Parameters:: cid (int) – Cell identifier.
Returns:: Refinement level (0 for unrefined root cells).
Return type:: int
Raises:: IndexError – If cid is out of range.

cell_multi_index(cid)[source]¶

Return the per-axis index of cell cid in its level’s coordinates.

Parameters:: cid (int) – Cell identifier.
Returns:: Length-ndim per-axis index tuple at the cell’s refinement level.
Return type:: tuple[int, …]
Raises:: IndexError – If cid is out of range.

is_active_leaf(level, midx)[source]¶

Return whether (level, midx) is an active (leaf) cell.

Parameters:

level (int) – Hierarchy level.
midx (Sequence[int]) – Per-axis index in level-level coordinates.

Returns:

True iff a cell with this level and multi-index is currently active (a leaf); False if it is out of range, not yet created, or has been refined away.

Return type:

bool

level_cells_per_axis(level)[source]¶

Return the per-axis cell count of the level-level grid.

This is a pure formula — level need not be an existing hierarchy level. Values above max_level return the count for the hypothetical finer grid that would result from additional uniform subdivision. This differs from active_blocks(), active_leaf_mask(), and subdomain_mask(), which all require level <= max_level.

Parameters:: level (int) – Hierarchy level. Must be >= 0; values above max_level are accepted and return the geometrically valid count.
Returns:: root.cells_per_axis[k] * factor[k] ** level for every axis k.
Return type:: tuple[int, …]
Raises:: ValueError – If level < 0.

active_blocks(level)[source]¶

Return the active-leaf blocks at level.

Parameters:: level (int) – Hierarchy level. Must be in [0, max_level].
Returns:: The sorted, non-overlapping (lo, hi) integer rectangles of active (leaf) cells at level, in level-level coordinates.
Return type:: tuple[tuple[tuple[int, …], tuple[int, …]], …]
Raises:: ValueError – If level is outside [0, max_level].

active_leaf_mask(level)[source]¶

Return a boolean mask of the active-leaf cells at level.

Parameters:: level (int) – Hierarchy level. Must be in [0, max_level].
Returns:: Fresh array of shape level_cells_per_axis(level); True where the level-level cell (level, midx) is an active (leaf) cell.
Return type:: npt.NDArray[np.bool_]
Raises:: ValueError – If level is outside [0, max_level].

subdomain_mask(level)[source]¶

Return a boolean mask of the level-level refined subdomain.

A level-level cell lies in the subdomain $\Omega_{level}$ (the region refined to at least level) iff it is not covered by an active leaf of a coarser level. The mask is computed at the level-level resolution by projecting every coarser active-leaf block up to level and clearing those cells.

Parameters:: level (int) – Hierarchy level. Must be in [0, max_level].
Returns:: Fresh array of shape level_cells_per_axis(level); True where the level-level cell lies in $\Omega_{level}$.
Return type:: npt.NDArray[np.bool_]
Raises:: ValueError – If level is outside [0, max_level].

Note

The mask is sized to the level-level cell grid and computed on demand; it is never stored.

refine(level, lo, hi)[source]¶

Refine the rectangular region [lo, hi) at level to level+1.

Union semantics: only the currently-active portion of [lo, hi) is refined. If the intersection with active blocks at level is empty, the call is a silent no-op. Overlapping calls therefore safely extend the refined region.

After the call all flat cell ids are reassigned (the BVH, cell tags, and facet tags are also invalidated).

Parameters:

level (int) – Refinement level at which the region lives. Must satisfy 0 <= level <= max_level.
lo (Sequence[int]) – Per-direction start index (inclusive), in level-level coordinates.
hi (Sequence[int]) – Per-direction end index (exclusive), in level-level coordinates.

Raises:

ValueError – If level is out of range, lo/hi have the wrong length, any lo[k] >= hi[k], or [lo, hi) falls entirely outside the level-level domain.

Return type:

None

References

Refinement and coarsening algorithms for adaptive hierarchical-spline meshes [Garau and Vázquez, 2018].

refine_cells(cell_ids)[source]¶

Refine a set of active cells using per-level bounding-box aggregation.

Groups cell_ids by their level, computes the bounding box of all cells at each level (the smallest axis-aligned rectangle containing them), and calls refine() once per level.

Parameters:: cell_ids (Sequence[int]) – Flat cell ids to refine. Cells from multiple levels are handled; repeated ids are silently ignored.
Raises:: IndexError – If any id in cell_ids is out of range.
Return type:: None

coarsen(level, lo, hi)[source]¶

Coarsen the rectangular region [lo, hi) at level (inverse of refine).

Reactivates the level-level cells in [lo, hi) and removes their level-(level+1) children. The region must be fully refined to exactly level ``level+1``: every child cell in [lo*factor, hi*factor) must be an active leaf at level+1 (none further refined, none still a leaf at level). Calling coarsen() with the same arguments as a preceding refine() exactly restores the grid.

After the call all flat cell ids are reassigned (the BVH, cell tags, and facet tags are invalidated).

Parameters:

level (int) – Level whose cells are reactivated. Must satisfy 0 <= level < max_level.
lo (Sequence[int]) – Per-direction start index (inclusive), in level-level coordinates.
hi (Sequence[int]) – Per-direction end index (exclusive), in level-level coordinates.

Raises:

ValueError – If level is out of range, lo/hi have the wrong length, any lo[k] >= hi[k], [lo, hi) is out of bounds, or the region is not fully refined to exactly level level+1.

Return type:

None

References

Coarsening algorithms for adaptive hierarchical-spline meshes [Garau and Vázquez, 2018].

collect_cell_bounds()[source]¶

Materialize (cell_lo, cell_hi) in flat-id order for the BVH.

Backed by a Numba kernel parallelizing over the active blocks; results are identical to calling cell_bounds() per cell.

Returns:: (cell_lo, cell_hi) of shape (num_cells, ndim).
Return type:: tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]

class pantr.grid.Partition(cell_owner, n_parts)[source]¶

Bases: object

A per-cell owner assignment over a grid’s cells.

Records, for every cell, the rank that owns it – or -1 for an inactive cell excluded from the partition (e.g. an exterior / trimmed cell). The owner array is coerced to a read-only int32 array on construction and the object is otherwise immutable. Owners and counts are exposed through the cell_owner, n_parts, n_cells, and active_mask properties.

Parameters:

cell_owner (npt.ArrayLike)
n_parts (int)

__init__(cell_owner, n_parts)[source]¶

Build a partition from a per-cell owner array.

Parameters:

cell_owner (npt.ArrayLike) – Per-cell owner ranks (-1 for inactive cells); coerced to a read-only 1D int32 array.
n_parts (int) – Number of parts (ranks); must be >= 1.

Raises:

ValueError – If n_parts < 1, cell_owner is not 1D integer, or any owner is outside [-1, n_parts).

Return type:

None

property cell_owner: npt.NDArray[np.int32]¶

Get the read-only per-cell owner array.

Returns:: (n_cells,) owners; -1 for inactive cells.
Return type:: npt.NDArray[np.int32]

property n_parts: int¶

Get the number of parts (ranks).

Returns:: The part count (>= 1).
Return type:: int

property n_cells: int¶

Get the total number of cells (active and inactive).

Returns:: Length of cell_owner.
Return type:: int

property active_mask: npt.NDArray[np.bool_]¶

Get a boolean mask of the active cells (owned by some rank).

Returns:: Fresh (n_cells,) mask; True where the cell owner is not -1.
Return type:: npt.NDArray[np.bool_]

owned_cells(rank)[source]¶

Return the flat ids of the cells owned by rank, ascending.

Parameters:: rank (int) – Owner rank in [0, n_parts).
Returns:: Sorted cell ids with cell_owner == rank.
Return type:: npt.NDArray[np.int64]
Raises:: ValueError – If rank is outside [0, n_parts).

class pantr.grid.TensorProductGrid(breakpoints)[source]¶

Bases: Grid

Tensor-product grid of axis-aligned boxes with per-axis breakpoints.

Cells are numbered in row-major (C) order over cells_per_axis (last axis varies fastest). See the module docstring for the footprint and construction notes. Size and geometry metadata are exposed through the ndim, num_cells, cells_per_axis, breakpoints, and bounds properties.

Parameters:: breakpoints (Sequence[npt.ArrayLike])

__init__(breakpoints)[source]¶

Build a tensor-product grid from per-axis breakpoint vectors.

Parameters:

breakpoints (Sequence[npt.ArrayLike]) – One strictly increasing float64 array-like per axis, each of length cells_per_axis[d] + 1 >= 2.

Raises:

ValueError – If breakpoints is empty, any axis has fewer than two entries, or any axis is non-finite or not strictly increasing.
TypeError – If a breakpoint array cannot be cast to float64.

Return type:

None

property ndim: int¶

Get the spatial dimension of the grid.

Returns:: Number of axes (>= 1).
Return type:: int

property num_cells: int¶

Get the total number of cells.

Returns:: Product of the per-axis cell counts.
Return type:: int

property cells_per_axis: tuple[int, ...]¶

Get the per-axis cell counts.

Returns:: Length-ndim tuple of per-axis counts.
Return type:: tuple[int, …]

property breakpoints: tuple[npt.NDArray[np.float64], ...]¶

Get the per-axis strictly increasing breakpoint arrays.

Returns:: Read-only float64 arrays; breakpoints[d] has length cells_per_axis[d] + 1.
Return type:: tuple[npt.NDArray[np.float64], …]

property bounds: npt.NDArray[np.float64]¶

Get the per-axis [lo, hi] extremes.

Returns:: Read-only (ndim, 2) array.
Return type:: npt.NDArray[np.float64]

property is_uniform: bool¶

Get whether every axis has uniform breakpoint spacing.

Returns:: True iff each axis’s spacing is constant to within an absolute tolerance (_UNIFORM_SPACING_ATOL).
Return type:: bool

cell_multi_index(cid)[source]¶

Return the per-axis indices (i_0, ..., i_{ndim-1}) of cell cid.

Parameters:: cid (int) – Flat cell identifier.
Returns:: Length-ndim per-axis index tuple (C-order).
Return type:: tuple[int, …]
Raises:: IndexError – If cid is out of range.

flat_cell_index(multi)[source]¶

Map per-axis cell indices to the flat cell identifier (C-order).

Parameters:

multi (Sequence[int]) – Length-ndim per-axis indices; each entry must satisfy 0 <= i_k < cells_per_axis[k].

Returns:

Flat cell identifier.

Return type:

int

Raises:

ValueError – If len(multi) != ndim.
IndexError – If any per-axis index is out of range.

cell_bounds(cid)[source]¶

Return cell cid’s axis-aligned (lo, hi) corners.

Parameters:: cid (int) – Cell identifier.
Returns:: Fresh, writeable length-ndim float64 arrays.
Return type:: tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]
Raises:: IndexError – If cid is out of range.

locate(pt)[source]¶

Return the cell containing pt, or None if pt is outside.

A point exactly on an interior breakpoint is assigned to the lower-indexed cell sharing that face; a point on the outer boundary is assigned to the adjacent boundary cell. Non-finite coordinates (NaN or infinity) are outside every cell.

Parameters:: pt (npt.ArrayLike) – Length-ndim point.
Returns:: Containing cell id, or None when pt lies outside the grid domain.
Return type:: int | None
Raises:: ValueError – If pt does not have length ndim.

locate_many(points)[source]¶

Locate a batch of points via the Numba per-axis search kernel.

Parameters:: points (npt.ArrayLike) – (npts, ndim) array-like of points, or a single length-ndim point.
Returns:: Shape (npts,) cell ids; -1 for points outside the grid (including points with NaN or infinite coordinates).
Return type:: npt.NDArray[np.int64]
Raises:: ValueError – If the trailing axis of points is not ndim.

neighbor_across_facet(cid, lfid)[source]¶

Return the cell across local facet lfid of cid, or None.

Uses the lfid = 2 * axis + side encoding and per-axis arithmetic.

Parameters:

cid (int) – Cell identifier.
lfid (int) – Local facet identifier in [0, 2 * ndim).

Returns:

Neighbouring cell id, or None on a boundary facet.

Return type:

int | None

Raises:

IndexError – If cid or lfid is out of range.

collect_cell_bounds()[source]¶

Materialize per-cell (lo, hi) in C-order via per-axis broadcasting.

Returns:: (cell_lo, cell_hi) of shape (num_cells, ndim) in cell-id order.
Return type:: tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]

restrict(cell_ids)[source]¶

Return the multi-index bounding-box sub-grid spanning cell_ids.

The window is the axis-aligned bounding box of the requested cells in multi-index space. Each axis of the sub-grid uses a pure slice of this grid’s breakpoints, breakpoints[d][imin_d : imax_d + 2] – never re-based or re-clamped – so the sub-grid’s cells coincide exactly with the matching cells of this grid.

Parameters:

cell_ids (npt.ArrayLike) – Flat cell identifiers to span; duplicates are ignored. Each must satisfy 0 <= cid < num_cells.

Returns:

The windowed TensorProductGrid, its local_to_global_cell map of shape (sub.num_cells,), and the in_subset mask flagging requested versus bounding-box-fill cells.

Return type:

GridRestriction

Raises:

ValueError – If cell_ids is empty.
IndexError – If any cell id is out of range [0, num_cells).
TypeError – If cell_ids is not integer-valued.

pantr.grid.cell_quadrature(grid, rule, cells=None)[source]¶

Map a reference QuadratureRule onto grid cells.

For each selected cell, the rule’s unit-cube points are affinely mapped into the cell box and its weights are scaled by the cell volume.

Parameters:

grid (Grid) – The grid whose cells are integrated. grid.ndim must equal rule.ndim.
rule (QuadratureRule) – Reference rule on [0, 1]^ndim.
cells (npt.ArrayLike | None) – Cell ids to map the rule onto. None (the default) selects every cell in id order. Otherwise a 1D integer array-like of ids, each in [0, num_cells) (a scalar is treated as a single id); order and duplicates are preserved.

Returns:

A pair (points, weights) where points has shape (num_selected, num_points, ndim) and weights has shape (num_selected, num_points), ordered to match cells.

Return type:

tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]

Raises:

ValueError – If rule.ndim != grid.ndim, cells is not a 1D integer array-like, or any id is outside [0, num_cells).

Note

Cell boxes for the whole grid are materialized once (an O(num_cells * ndim) temporary), then the requested subset is selected; this keeps the common all-cells path vectorized.

pantr.grid.hierarchical_grid(root, factor)[source]¶

Build a HierarchicalGrid from a root grid and a subdivision factor.

Parameters:

root (TensorProductGrid) – The level-0 grid.
factor (int | Sequence[int]) – Per-direction subdivision factor. A scalar is broadcast to all axes. Each entry must be >= 1.

Returns:

A new hierarchical grid starting with all root cells active at level 0.

Return type:

HierarchicalGrid

Raises:

TypeError – If root is not a TensorProductGrid.
ValueError – If factor has the wrong length or any entry is < 1.

pantr.grid.overlay(grid_a, grid_b)[source]¶

Return the tensor-product overlay of grid_a and grid_b.

The overlay’s per-axis breakpoints are the sorted union of both inputs’ breakpoint arrays restricted to the intersection of their domains. Breakpoints closer than the default float64 tolerance (pantr.tolerance.get_default()) are merged into one. The result is the coarsest TensorProductGrid that refines both inputs.

Parameters:

grid_a (TensorProductGrid) – First input grid.
grid_b (TensorProductGrid) – Second input grid; must share ndim with grid_a and have a non-empty domain intersection on every axis.

Returns:

The overlay grid.

Return type:

TensorProductGrid

Raises:

TypeError – If either argument is not a TensorProductGrid.
ValueError – If the grids have different ndim, or if their domains do not overlap on some axis (the per-axis intersection is empty or degenerate).

Example

>>> import numpy.testing as npt
>>> from pantr.grid import uniform_grid, overlay
>>> a = uniform_grid([[0.0, 1.0]], 2)
>>> b = uniform_grid([[0.0, 1.0]], 3)
>>> npt.assert_allclose(
...     overlay(a, b).breakpoints[0],
...     [0.0, 1/3, 0.5, 2/3, 1.0],
... )

pantr.grid.partition_grid(grid, n_parts, *, backend='auto', cell_weights=None, cell_active=None)[source]¶

Partition a grid’s cells into n_parts rank subdomains.

Returns a Partition (a plain per-cell owner array) for the serial, communication-free distribution of a grid – no MPI is involved. See the module docstring for the "block", "rcb", and "auto" backends.

The cell_weights and cell_active hooks let a consumer (e.g. an immersed/unfitted code) drive the partition without pantr storing any classification: per-cell assembly cost via cell_weights (the "rcb" backend balances total weight), and an active subset via cell_active (inactive cells get owner -1). The "block" backend supports neither.

Parameters:

grid (Grid) – The grid to partition. The "block" backend requires a TensorProductGrid; "rcb" accepts any grid.
n_parts (int) – Number of parts (ranks); must be >= 1.
backend (str) – "auto" (default), "block", or "rcb".
cell_weights (npt.ArrayLike | None) – Optional per-cell cost, shape (num_cells,), finite and non-negative. None means uniform. Used by "rcb"; rejected by "block".
cell_active (npt.ArrayLike | None) – Optional boolean mask, shape (num_cells,); inactive cells get owner -1 and are excluded from partitioning. None means all active. Used by "rcb"; rejected by "block".

Returns:

A per-cell owner assignment with n_parts parts. Cells excluded by cell_active have owner -1; every active cell is assigned a rank in range(n_parts) and no rank is left empty.

Return type:

Partition

Raises:

ValueError – If n_parts < 1; if backend is unknown; if "block" is used on a non-TensorProductGrid or with weights/activity; if n_parts cannot be factored onto the axes ("block"); if cell_weights / cell_active have the wrong shape or invalid values; or if n_parts exceeds the number of active cells ("rcb").

Example

>>> from pantr.grid import partition_grid, uniform_grid
>>> grid = uniform_grid([[0.0, 1.0], [0.0, 1.0]], [4, 4])
>>> part = partition_grid(grid, 4)
>>> part.n_parts
4
>>> int(part.cell_owner.min()), int(part.cell_owner.max())
(0, 3)

pantr.grid.tensor_product_grid(space)[source]¶

Build the knot-span grid of a pantr.bspline.BsplineSpace.

The grid’s per-axis breakpoints are the unique in-domain knots of each 1-D sub-space, so its cells are exactly the space’s knot spans and its cell ids match pantr.bspline.SpanwiseElementExtraction on the same space.

Parameters:: space (BsplineSpace) – A tensor-product B-spline space. Periodic directions are rejected: a periodic knot vector does not map cleanly to a bounded grid.
Returns:: A grid whose cells are the knot spans of space.
Return type:: TensorProductGrid
Raises:: ValueError – If any direction of space is periodic.

pantr.grid.uniform_grid(bounds, cells)[source]¶

Build a uniform TensorProductGrid on a bounding box.

Parameters:

bounds (npt.ArrayLike) – (ndim, 2) array-like of per-axis [lo, hi] pairs with ndim >= 1.
cells (int | Sequence[int]) – Number of cells per axis. A scalar is broadcast to every axis; a length-ndim sequence gives per-axis counts. Each count must be >= 1.

Returns:

A uniform tensor-product grid.

Return type:

TensorProductGrid

Raises:

ValueError – If bounds does not have shape (ndim, 2), any axis has lo >= hi, cells has the wrong length, or any count is < 1.
TypeError – If bounds cannot be cast to float64.

Quadrature rules and evaluation grid helpers for 1D integration.

This module provides:

1D quadrature rules on [0, 1]: trapezoidal (equispaced), Gauss-Legendre, Gauss-Lobatto-Legendre, Chebyshev-Gauss (1st and 2nd kind).
PointsLattice: a multi-dimensional tensor-product evaluation grid.
create_lagrange_points_lattice(): factory for Lagrange-node lattices.
QuadratureRule: a d-dimensional quadrature rule on the unit cube, with tensor_product_quadrature() and gauss_legendre_quadrature() factories; the reference rule consumed by pantr.grid.cell_quadrature().

class pantr.quad.PointsLattice(pts_per_dir)[source]¶

Bases: object

A tensor-product grid of evaluation points in multiple dimensions.

Stores one 1D array of coordinates per spatial direction and provides helpers for constructing the full set of grid points or querying grid metadata.

Parameters:: pts_per_dir (Iterable[npt.NDArray[np.float32 | np.float64]])

_pts_per_dir¶

One 1D coordinate array per spatial dimension. All arrays share the same dtype.

Type:: tuple[npt.NDArray[np.float32 | np.float64], …]

__init__(pts_per_dir)[source]¶

Initialize the points lattice.

Parameters:: pts_per_dir (Iterable[npt.NDArray[np.float32 | np.float64]]) – The points per dimension. All points must have the same dtype.
Raises:: ValueError – If the dimension is less than 1 or the points have different dtypes.
Return type:: None

property dim: int¶

Get the dimension of the points lattice.

Returns:: Number of spatial dimensions.
Return type:: int

property dtype: DTypeLike¶

Get the dtype of the points lattice.

Returns:: The numpy floating-point dtype shared by all coordinate arrays.
Return type:: npt.DTypeLike

property pts_per_dir: tuple[ndarray[tuple[Any, ...], dtype[float32 | float64]], ...]¶

Get the points per dimension.

Returns:: One 1D coordinate array for each spatial dimension.
Return type:: tuple[npt.NDArray[np.float32 | np.float64], …]

get_all_points(order='C')[source]¶

Get all points in the points lattice.

Parameters:

order (Literal["C", "F"]) – The order of the points. Defaults to “C”. “C” means the last index varies fastest, “F” means the first index varies fastest.

Returns:

The dim-dimensional points: in the lattice. It has shape: (n_pts, dim).

Return type:

npt.NDArray[np.float32 | np.float64]

class pantr.quad.QuadratureRule(points, weights)[source]¶

Bases: object

Immutable quadrature rule on the unit cube [0, 1]^ndim.

Bundles quadrature points and weights as the reference rule that pantr.grid.cell_quadrature() affinely maps onto each cell of a grid. Points lie in the closed unit cube; the factory-built rules (tensor_product_quadrature(), gauss_legendre_quadrature()) have weights summing to 1 (the measure of the unit cube), so the rule integrates the constant 1 exactly. The stored arrays are read-only.

Parameters:

points (npt.ArrayLike)
weights (npt.ArrayLike)

__init__(points, weights)[source]¶

Build and validate a quadrature rule on the unit cube.

Parameters:

points (npt.ArrayLike) – (num_points, ndim) array-like; every coordinate must lie in [0, 1].
weights (npt.ArrayLike) – (num_points,) array-like of weights.

Raises:

ValueError – If points is not 2D, weights is not 1D, their lengths disagree, either is empty, any value is non-finite, or any point lies outside [0, 1].

Return type:

None

property ndim: int¶

Get the spatial dimension of the rule.

Returns:: Number of axes (>= 1).
Return type:: int

property num_points: int¶

Get the number of quadrature points.

Returns:: Point count (>= 1).
Return type:: int

property points: ndarray[tuple[Any, ...], dtype[float64]]¶

Get the quadrature points on the unit cube.

Returns:: Read-only (num_points, ndim) array in [0, 1]^ndim.
Return type:: npt.NDArray[np.float64]

property weights: ndarray[tuple[Any, ...], dtype[float64]]¶

Get the quadrature weights.

Returns:: Read-only (num_points,) array.
Return type:: npt.NDArray[np.float64]

pantr.quad.create_lagrange_points_lattice(lagrange_variant, n_pts_per_dir, dtype=<class 'numpy.float64'>)[source]¶

Create a Lagrange points lattice for tensor-product evaluation.

Builds a PointsLattice whose per-direction coordinate arrays are the Lagrange nodes of the specified variant on [0, 1].

Parameters:

lagrange_variant (LagrangeVariant) – The variant of the Lagrange basis (e.g., equispaced, Gauss-Legendre, Gauss-Lobatto-Legendre, etc.).
n_pts_per_dir (Iterable[int]) – Number of points per spatial dimension. Each value must be at least 1.
dtype (npt.DTypeLike) – Floating-point dtype for the coordinates. Must be float32 or float64. Defaults to np.float64.

Returns:

A lattice whose per-direction coordinate arrays are the Lagrange nodes for the given variant and point counts.

Return type:

PointsLattice

Raises:

ValueError – If any value in n_pts_per_dir is less than 1.

pantr.quad.gauss_legendre_quadrature(ndim, npts)[source]¶

Build a tensor-product Gauss-Legendre QuadratureRule on the unit cube.

Parameters:

ndim (int) – Number of axes (>= 1).
npts (int | Sequence[int]) – Points per axis. A scalar is broadcast to every axis; a length-ndim sequence gives per-axis counts. Each count must be >= 1.

Returns:

The tensor product of per-axis npts-point Gauss-Legendre rules. Exact for tensor-product polynomials of per-axis degree 2 * npts - 1; weights sum to 1.

Return type:

QuadratureRule

Raises:

ValueError – If ndim < 1, npts is a sequence of the wrong length, or any count is < 1.

References

Nodes and weights follow the classical Gauss-Legendre construction [Golub and Welsch, 1969].

pantr.quad.get_chebyshev_gauss_1st_kind_1d(n_pts, dtype=<class 'numpy.float64'>)[source]¶

Get Chebyshev-Gauss quadrature of the first kind on [0, 1] for the given number of points.

If n_pts == 1, the nodes are [0.5] and the weights are [1.0].

Parameters:

n_pts (int) – Number of quadrature points. Must be at least 1.
dtype (npt.DTypeLike) – Floating dtype for nodes/weights; float32 or float64. Defaults to float64.

Returns:

The nodes and weights.

Return type:

tuple[npt.NDArray[np.float32 | np.float64], npt.NDArray[np.float32 | np.float64]]

Raises:

ValueError – If n_pts is less than 1 or dtype is not float32 or float64.

pantr.quad.get_chebyshev_gauss_2nd_kind_1d(n_pts, dtype=<class 'numpy.float64'>)[source]¶

Get Chebyshev-Gauss quadrature of the second kind on [0, 1] for the given number of points.

The rule integrates against the Chebyshev second-kind weight function mapped to [0, 1]: $\int_0^1 f(x) \sqrt{1 - (2x - 1)^2}\, dx \approx \sum_k w_k f(x_k)$, exactly for polynomials of degree up to 2 * n_pts - 1. The nodes are the mapped roots of the Chebyshev polynomial of the second kind $U_{n_pts}$ – all interior (no endpoints). For endpoint-including Chebyshev-Lobatto interpolation nodes, use get_modified_chebyshev_nodes_1d() instead.

Parameters:

n_pts (int) – Number of quadrature points. Must be at least 2.
dtype (npt.DTypeLike) – Floating dtype for nodes/weights; float32 or float64. Defaults to float64.

Returns:

The nodes (ascending, interior) and weights.

Return type:

tuple[npt.NDArray[np.float32 | np.float64], npt.NDArray[np.float32 | np.float64]]

Raises:

ValueError – If n_pts is less than 2 or dtype is not float32 or float64.

pantr.quad.get_gauss_legendre_1d(n_pts, dtype=<class 'numpy.float64'>)[source]¶

Get Gauss-Legendre quadrature nodes on [0, 1] for the given number of points.

Parameters:

n_pts (int) – The number of points. Must be at least 1.
dtype (npt.DTypeLike) – The dtype of the nodes. If must be float32 or float64. Defaults to float64.

Returns:

The nodes and weights.

Return type:

tuple[npt.NDArray[np.float32 | np.float64], npt.NDArray[np.float32 | np.float64]]

Raises:

ValueError – If n_pts is less than 1 or dtype is not float32 or float64.

pantr.quad.get_gauss_lobatto_legendre_1d(n_pts, dtype=<class 'numpy.float64'>)[source]¶

Get Gauss-Lobatto-Legendre quadrature nodes on [0, 1] for the given number of points.

Parameters:

n_pts (int) – The number of points. Must be at least 2.
dtype (npt.DTypeLike) – The dtype of the nodes. If must be float32 or float64. Defaults to float64.

Returns:

The nodes and weights.

Return type:

tuple[npt.NDArray[np.float32 | np.float64], npt.NDArray[np.float32 | np.float64]]

Raises:

ValueError – If n_pts is less than 2 or dtype is not float32 or float64.

pantr.quad.get_modified_chebyshev_nodes_1d(n_pts, dtype=<class 'numpy.float64'>)[source]¶

Get modified Chebyshev nodes on [0, 1] for the given number of points.

Returns Chebyshev nodes of the second kind (Chebyshev-Lobatto points) mapped to [0, 1]. These include both endpoints and are suitable for polynomial interpolation into the Bernstein basis.

Unlike the quadrature functions in this module, this returns only nodes (no weights) since it is intended for interpolation, not integration.

Parameters:

n_pts (int) – Number of nodes. Must be at least 2.
dtype (npt.DTypeLike) – Floating dtype; float32 or float64. Defaults to float64.

Returns:

Array of shape (n_pts,) with nodes in [0, 1], starting at 0 and ending at 1.

Return type:

npt.NDArray[np.float32 | np.float64]

Raises:

ValueError – If n_pts < 2 or dtype is not float32 or float64.

pantr.quad.get_tanh_sinh_1d(n_pts, dtype=<class 'numpy.float64'>)[source]¶

Get tanh-sinh quadrature nodes on [0, 1] for the given number of points.

Tanh-sinh (double-exponential) quadrature clusters nodes near the endpoints of the interval, making it well suited for integrands with endpoint singularities or steep boundary layers. The scheme is symmetric and nodes near the endpoints that are indistinguishable from 0 or 1 in floating-point arithmetic are snapped to the boundary, so the effective number of returned nodes may be less than n_pts.

Parameters:

n_pts (int) – Requested number of quadrature points. Must be at least 1.
dtype (npt.DTypeLike) – Floating-point dtype for the output arrays. Must be float32 or float64. Defaults to float64.

Returns:

(nodes, weights) on [0, 1]. Both arrays have the same length, which may be less than n_pts due to endpoint snapping. Weights sum to 1.

Return type:

tuple[npt.NDArray[np.float32 | np.float64], npt.NDArray[np.float32 | np.float64]]

Raises:

ValueError – If n_pts < 1 or dtype is not float32/float64.

Example

>>> nodes, weights = get_tanh_sinh_1d(5)
>>> nodes.shape[0] <= 5
True
>>> abs(weights.sum() - 1.0) < 1e-14
True

pantr.quad.get_trapezoidal_1d(n_pts, dtype=<class 'numpy.float64'>)[source]¶

Get trapezoidal quadrature nodes on [0, 1] for the given number of points.

If n_pts == 1, the nodes are [0.5] and the weights are [1.0].

Parameters:

n_pts (int) – The number of points. Must be at least 1.
dtype (npt.DTypeLike) – The dtype of the nodes. If must be float32 or float64. Defaults to float64.

Returns:

The nodes and weights.

Return type:

tuple[npt.NDArray[np.float32 | np.float64], npt.NDArray[np.float32 | np.float64]]

Raises:

ValueError – If n_pts is less than 1 or dtype is not float32 or float64.

pantr.quad.tensor_product_quadrature(rules)[source]¶

Build a tensor-product QuadratureRule from per-axis 1D rules.

Each axis contributes a 1D rule (nodes, weights) on [0, 1]; the d-dimensional rule is their tensor product. Points are enumerated in row-major (C) order – the last axis varies fastest, matching pantr.grid.TensorProductGrid cell ids – and each weight is the product of the corresponding per-axis weights.

Parameters:: rules (Sequence[tuple[npt.ArrayLike, npt.ArrayLike]]) – One (nodes, weights) pair per axis. Within a pair, nodes and weights must be 1D of equal, non-zero length, with nodes in [0, 1].
Returns:: The tensor-product rule on [0, 1]^ndim with ndim == len(rules) and num_points the product of the per-axis point counts.
Return type:: QuadratureRule
Raises:: ValueError – If rules is empty, or any axis pair is not a matching pair of non-empty 1D arrays, or (via QuadratureRule) any node lies outside [0, 1].

Tolerance utilities for floating-point comparisons in IGA applications.

class pantr.tolerance.ToleranceInfo[source]¶

Bases: TypedDict

A TypedDict holding comprehensive tolerance and precision information.

pantr.tolerance.get_conservative(dtype)[source]¶

Get a conservative tolerance for robust floating-point comparisons.

Parameters:

dtype (npt.DTypeLike) – NumPy floating-point data type.

Returns:

Conservative tolerance value for the given dtype. Used when: robustness is more important than precision.

Return type:

float

Raises:

ValueError – If dtype is not a supported floating-point type.

pantr.tolerance.get_default(dtype)[source]¶

Get a reasonable default tolerance for floating-point comparisons.

Parameters:: dtype (npt.DTypeLike) – NumPy floating-point data type or numpy scalar type.
Returns:: Recommended tolerance value for the given dtype.
Return type:: float
Raises:: ValueError – If dtype is not a supported floating-point type.

Example

>>> get_default(np.float32)
1e-06
>>> get_default("float64")
1e-12

pantr.tolerance.get_info(dtype)[source]¶

Get comprehensive tolerance information for a dtype.

Parameters:

dtype (npt.DTypeLike) – NumPy floating-point data type.

Returns:

Dictionary containing tolerance information including: machine epsilon, default/strict/conservative tolerances, precision bits, and min/max values for the dtype.

Return type:

ToleranceInfo

Raises:

ValueError – If dtype is not a supported floating-point type.

pantr.tolerance.get_machine_epsilon(dtype)[source]¶

Get machine epsilon for a given floating-point dtype.

Machine epsilon is the smallest positive number that, when added to 1.0, produces a result different from 1.0. It represents the relative error in floating-point arithmetic for the given precision.

Parameters:: dtype (npt.DTypeLike) – NumPy floating-point data type.
Returns:: Machine epsilon for the given dtype.
Return type:: float
Raises:: ValueError – If dtype is not a supported floating-point type.

pantr.tolerance.get_strict(dtype)[source]¶

Get a strict tolerance for high-precision floating-point comparisons.

Parameters:

dtype (npt.DTypeLike) – NumPy floating-point data type.

Returns:

Strict tolerance value for the given dtype. Typically used for: parametric coordinates requiring high precision.

Return type:

float

Raises:

ValueError – If dtype is not a supported floating-point type.

Parallelism & visualization¶

Thread-pool control for the Numba kernels (see Parallelism):

pantr.get_num_threads()[source]¶

Return the current number of threads used by parallel kernels.

Returns:: Active Numba thread-pool size.
Return type:: int

pantr.set_num_threads(n)[source]¶

Set the number of threads used by parallel kernels.

Calling this marks the thread count as explicitly configured: default policies (e.g. the per-rank MPI default in pantr.mpi) will never override it afterwards.

Parameters:: n (int) – Desired thread count. Must be >= 1 and at most numba.config.NUMBA_NUM_THREADS.
Raises:: ValueError – If n is less than 1 or exceeds the maximum.
Return type:: None

pantr.num_threads(n, *, limit_blas=False)[source]¶

Context manager that temporarily sets the parallel thread count.

On entry the Numba thread-pool size is changed to n; on exit the previous value is restored. When limit_blas is True, BLAS/LAPACK thread pools (OpenBLAS, MKL, …) are also limited to n threads for the duration of the block, via threadpoolctl.

Parameters:

n (int) – Desired thread count for the block.
limit_blas (bool) – If True, also limit BLAS/LAPACK threads via threadpoolctl. Defaults to False.

Yields:

None

Return type:

Generator[None, None, None]

Note

Using this context manager permanently marks the thread count as explicitly configured (via set_num_threads() on both entry and exit). Default policies – e.g. the per-rank MPI default in pantr.mpi – will not override the thread count for the rest of the process lifetime, even after the with block exits.

Example

>>> with pantr.num_threads(1):
...     # all pantr operations run serially here
...     pass

Optional MPI-parallel distribution layer for PaNTr.

This subpackage hosts the MPI-dependent and dolfinx-bridge code for distributing B-spline and THB-spline spaces across ranks. It ships with every install of PaNTr but is backend-gated: the serial core (pantr.grid, pantr.bspline, …) never imports it (enforced by an import-linter contract), and import pantr.mpi succeeds even when mpi4py is absent – only helpers that genuinely need MPI raise at call time, via require_mpi().

mpi4py is the optional backend, not extra PaNTr code: a plain pip install pantr is serial-only and MPI-free, while pip install "pantr[mpi]" simply installs mpi4py for you (equivalently pip install mpi4py; an MPI library is also required).

Main exports: - HAS_MPI: whether mpi4py was importable at module load time. - mpi_available(): live runtime check for mpi4py availability. - require_mpi(): lazily import and return the mpi4py.MPI module, or raise. - from_dolfinx(): build a pantr.grid.Partition from a dolfinx mesh. - DistributedSpace: the per-rank handle to an MPI-distributed space. - create_distributed_space(): build a DistributedSpace from a space. - DistributedFunction: the per-rank handle to an MPI-distributed function. - create_distributed_function(): build a DistributedFunction from a function. - configure_threads(): explicitly set this rank’s Numba thread count. - interpolate_bspline_distributed(): distributed B-spline collocation interpolation. - fit_bspline_distributed(): distributed B-spline fit from pre-evaluated values. - l2_project_bspline_distributed(): distributed B-spline L2 projection. - quasi_interpolate_bspline_distributed(): distributed B-spline quasi-interpolation. - quasi_interpolate_thb_spline_distributed(): distributed THB-spline quasi-interpolation.

A process-level thread policy coordinates MPI with PaNTr’s Numba parallelism: the first use of any MPI-engaging entry point limits this process to one Numba thread per rank, unless threads were explicitly configured (NUMBA_NUM_THREADS, pantr.set_num_threads, pantr.num_threads, or configure_threads()). Every new MPI entry point added to this package must call _ensure_default_thread_policy() from pantr.mpi._thread_policy before any other work – there is no structural chokepoint, since communicators are duck-typed.

pantr.mpi.HAS_MPI: Final[bool] = True¶

Whether mpi4py was importable at module load time.

This value is frozen at import and will not reflect environment changes made afterwards. Use mpi_available() for a live probe. In code paths that actually need MPI, always call require_mpi() — it checks freshly and raises a clear error if MPI is unavailable or broken.

Type:: bool

class pantr.mpi.DistributedFunction(global_function, distributed_space)[source]¶

Bases: object

Per-rank handle to an MPI-distributed spline function.

Pairs a DistributedSpace with a global control-point field and holds this rank’s local function – the windowed space’s Bspline / THBSpline whose control points are the global field sliced to the rank’s local DOFs. The local function equals the global one pointwise over the rank’s owned cells. Construction performs no MPI communication (it reuses the distributed space’s windowing).

Parameters:

global_function (Bspline | THBSpline)
distributed_space (DistributedSpace)

_global_function¶

The undistributed global function.

Type:: Bspline | THBSpline

_distributed_space¶

The distributed space it is defined on.

Type:: DistributedSpace

_local¶

This rank’s local function, or None.

Type:: Bspline | THBSpline | None

__init__(global_function, distributed_space)[source]¶

Build the distributed function by localizing the global control points.

Parameters:

global_function (Bspline | THBSpline) – The global function, identical on every rank. Its space must be the one distributed_space distributes.
distributed_space (DistributedSpace) – The distributed space whose global_space is exactly global_function.space.

Raises:

ValueError – If global_function.space is not the distributed space’s global_space.

Return type:

None

property global_function: Bspline | THBSpline¶

Get the undistributed global function.

Returns:: The global function passed at construction.
Return type:: Bspline | THBSpline

property distributed_space: DistributedSpace¶

Get the distributed space this function is defined on.

Returns:: The underlying distributed space.
Return type:: DistributedSpace

property local: Bspline | THBSpline | None¶

Get this rank’s local function, or None if it owns no cells.

Returns:: The rank-local function on the windowed space (control points sliced to local DOFs), or None when owns_cells is False. Scalar vs. vector kind is preserved.
Return type:: Bspline | THBSpline | None

property rank: int¶

Get this rank’s id within the communicator.

Returns:: The rank id (== distributed_space.rank).
Return type:: int

property n_parts: int¶

Get the number of ranks.

Returns:: The number of parts (== distributed_space.n_parts).
Return type:: int

property owns_cells: bool¶

Report whether this rank owns at least one cell.

Returns:: True iff local is not None.
Return type:: bool

class pantr.mpi.DistributedSpace(global_space, partition, comm)[source]¶

Bases: object

A B-spline or THB-spline space distributed across an MPI communicator.

The per-rank, SPMD handle: every rank constructs its own DistributedSpace from the same global space and partition, and holds only its own LocalSpace (local). A rank that owns no cells (an over-provisioned run, or a partition that excludes some ranks) has local set to None and owns_cells False instead of failing.

Construction performs no MPI communication: the partition is assumed identical on every rank (guaranteed by pantr’s deterministic partitioners and by from_dolfinx()), so each rank can window the global space locally.

The public surface exposes:

comm – the MPI communicator passed at construction.
rank – this rank’s id within the communicator.
n_parts – number of ranks the space is distributed over.
global_space – the undistributed global space.
partition – the cell-ownership partition.
local – this rank’s windowed local space (None if it owns no cells).
owns_cells – whether this rank owns at least one cell.
owned_cells – read-only sorted global cell ids owned by this rank.

Parameters:

global_space (BsplineSpace | THBSplineSpace)
partition (Partition)
comm (Any)

__init__(global_space, partition, comm)[source]¶

Build the distributed space, windowing the global space to this rank.

Parameters:

global_space (BsplineSpace | THBSplineSpace) – The global (non-periodic) space, identical on every rank.
partition (Partition) – Owner of every cell of the space’s grid, identical on every rank. Its n_parts must equal comm.size.
comm (Any) – An MPI communicator (e.g. mpi4py.MPI.COMM_WORLD). Only its rank and size attributes are read.

Raises:

ValueError – If global_space is a periodic BsplineSpace; if comm.size != partition.n_parts; if partition does not match the global space’s cell count; or if comm.rank is outside [0, comm.size) (delegated to owned_cells()).

Return type:

None

Note

Unlike build_local(), construction succeeds when this rank owns no cells – local is set to None and owns_cells to False.

Note

Construction engages the default MPI thread policy: unless threads were explicitly configured, this process’s Numba thread pool is limited to one thread per rank. See pantr.mpi.configure_threads().

property comm: Any¶

Get the MPI communicator this space is distributed over.

Returns:: The communicator passed at construction.
Return type:: Any

property rank: int¶

Get this rank’s id within the communicator.

Returns:: comm.rank.
Return type:: int

property n_parts: int¶

Get the number of ranks (parts) the space is distributed over.

Returns:: partition.n_parts (equal to comm.size at construction time).
Return type:: int

property global_space: BsplineSpace | THBSplineSpace¶

Get the undistributed global space.

Returns:: The global space passed at construction.
Return type:: BsplineSpace | THBSplineSpace

property partition: Partition¶

Get the cell-ownership partition.

Returns:: The partition passed at construction.
Return type:: Partition

property local: LocalSpace | None¶

Get this rank’s local windowed space, or None if it owns no cells.

Returns:: The rank-local space (windowed space plus local-to-global cell/DOF maps and ownership masks), or None when owns_cells is False. See LocalSpace for all fields.
Return type:: LocalSpace | None

Note

Callers must narrow the type before use: assert ds.local is not None or if ds.local is not None. owns_cells does not statically narrow this property in mypy.

property owns_cells: bool¶

Report whether this rank owns at least one cell.

Returns:: True iff local is not None.
Return type:: bool

property owned_cells: npt.NDArray[np.int64]¶

Get the global ids of the cells this rank owns.

Returns:: Read-only sorted global cell ids owned by this rank (empty if the rank owns none).
Return type:: npt.NDArray[np.int64]

localize(control_points)[source]¶

Restrict a global coefficient field to this rank’s local function.

Slices the global control points (identical on every rank, one per global DOF) down to this rank’s local DOFs via local’s local_to_global_dof map, and wraps them on the windowed local space. The returned function equals the global one pointwise over the rank’s owned cells, so per-element evaluation and assembly are local. Reuse a single distributed space to localize many fields (right-hand side, solution, residual, …).

Parameters:: control_points (npt.ArrayLike) – Global control points, shape (n_global_dofs,) for a scalar field or (n_global_dofs, rank) for a vector field, where n_global_dofs == global_space.num_total_basis.
Returns:: This rank’s local function (a Bspline for a tensor-product space, a THBSpline for a hierarchical one), or None if the rank owns no cells. Scalar vs. vector kind is preserved.
Return type:: Bspline | THBSpline | None
Raises:: ValueError – If control_points’s leading dimension is not global_space.num_total_basis.

pantr.mpi.configure_threads(threads_per_rank=1, *, limit_blas=False)[source]¶

Explicitly set the number of Numba threads this process (MPI rank) uses.

The explicit alternative to the default policy: call it on every rank – at any point, including before any pantr.mpi object is built – to choose how many threads each rank may use in a hybrid MPI + threads run. Marks the thread count as explicitly configured, so the default one-thread policy never overrides it. Also the right call when running serial PaNTr under your own mpiexec without any pantr.mpi machinery.

Parameters:

threads_per_rank (int) – Numba threads this process may use. Must be >= 1 and at most numba.config.NUMBA_NUM_THREADS. Defaults to 1.
limit_blas (bool) – If True, also persistently limit this process’s BLAS/LAPACK thread pools (OpenBLAS, MKL, …) to threads_per_rank via the optional threadpoolctl package; emits a UserWarning when it is not installed. A successive call replaces the previous limit. Defaults to False.

Raises:

ValueError – If threads_per_rank is less than 1 or exceeds numba.config.NUMBA_NUM_THREADS.

Return type:

None

Note

Numba’s thread mask is thread-local: this call governs kernels launched from the calling thread only. Call it on each thread individually if running PaNTr kernels from multiple threads; in SPMD practice, calling it once on the rank’s main thread is sufficient.

pantr.mpi.create_distributed_function(global_function, comm, *, method='grid', backend=None, cell_weights=None, cell_active=None)[source]¶

Build an MPI-distributed function directly from a global function.

The function-level counterpart of create_distributed_space(): distributes the function’s space across comm and slices the control points to each rank, returning a DistributedFunction whose local is this rank’s Bspline / THBSpline.

Parameters:

global_function (Bspline | THBSpline) – The global function to distribute, identical on every rank.
comm (Any) – An MPI communicator (e.g. mpi4py.MPI.COMM_WORLD); only its rank and size are read.
method (str) – Partitioning strategy, "grid" (default) or "graph". See create_distributed_space().
backend (str | None) – Partitioner backend; None selects each method’s default. See create_distributed_space().
cell_weights (npt.ArrayLike | None) – Per-cell assembly-cost weights. Defaults to None.
cell_active (npt.ArrayLike | None) – Boolean per-cell activity mask. Defaults to None.

Returns:

The per-rank distributed-function handle.

Return type:

Raises:

ValueError – If method/backend are invalid, or the derived partition is incompatible with comm (as raised downstream).

Example

>>> from mpi4py import MPI  
>>> from pantr.bspline import Bspline, create_uniform_space  
>>> from pantr.mpi import create_distributed_function  
>>> space = create_uniform_space([2, 2], [8, 8])  
>>> f = Bspline(space, coeffs)  
>>> df = create_distributed_function(f, MPI.COMM_WORLD)  
>>> local_f = df.local  # this rank's windowed function  

pantr.mpi.create_distributed_space(global_space, comm, *, method='grid', backend=None, cell_weights=None, cell_active=None)[source]¶

Build an MPI-distributed space directly from a global space.

Convenience wrapper over the explicit flow (derive grid -> partition -> DistributedSpace). It derives the space’s cell grid (tensor_product_grid() for a BsplineSpace, or THBSplineSpace.grid for a hierarchical space), partitions it across comm.size ranks, and returns the per-rank handle. The partition is deterministic, so every rank computes the same one without communication.

The explicit three-step flow remains available for full control; this factory just removes the boilerplate for the common case.

Parameters:

global_space (BsplineSpace | THBSplineSpace) – The global space to distribute, identical on every rank. A BsplineSpace must be non-periodic.
comm (Any) – An MPI communicator (e.g. mpi4py.MPI.COMM_WORLD). Only its rank and size are read; it is duck-typed.
method (str) – Partitioning strategy. "grid" (default) splits the cell grid geometrically via partition_grid(); "graph" partitions the cell-coupling graph (coupling_graph() + partition_graph()) to minimize cross-rank DOF coupling.
backend (str | None) – Partitioner backend. None (default) selects each method’s default – "auto" for "grid" ("block" or "rcb"), "spectral" for "graph" (or "metis", needing the metis extra).
cell_weights (npt.ArrayLike | None) – Per-cell assembly-cost weights. For "grid" they balance the geometric split; for "graph" they become the coupling graph’s vertex weights. Defaults to None (uniform).
cell_active (npt.ArrayLike | None) – Boolean per-cell activity mask; inactive cells get owner -1 and drop out of the partition. Defaults to None (all active).

Returns:

The per-rank distributed-space handle.

Return type:

DistributedSpace

Raises:

TypeError – If method="grid" and global_space is neither a BsplineSpace nor a THBSplineSpace.
ValueError – If method is not "grid" or "graph"; if backend is invalid for the chosen method; or if the derived partition is incompatible with comm (e.g. comm.size mismatch), as raised downstream.

Example

>>> from mpi4py import MPI  
>>> from pantr.bspline import create_uniform_space  
>>> from pantr.mpi import create_distributed_space  
>>> space = create_uniform_space([2, 2], [8, 8])  
>>> ds = create_distributed_space(space, MPI.COMM_WORLD)  

pantr.mpi.fit_bspline_distributed(values, nodes, distributed_space, *, values_distributed=False, tol=None)[source]¶

Construct a distributed B-spline from pre-evaluated sample values at known nodes.

The MPI-parallel counterpart of fit_bspline(). The collocation solve is cheap (small 1D SVD solves); for tensor-product nodes the value field is the only large object, so this function supports values that arrive already distributed (one contiguous block of the flattened grid per rank), gathered with a single allgather before the replicated Kronecker solve.

Two value layouts are supported via values_distributed:

values_distributed=False (default, replicated): every rank passes the full value field (shape (*n_pts_per_dir) or (*n_pts_per_dir, rank)), exactly as serial fit_bspline(). No allgather is needed; each rank solves locally.
values_distributed=True (pre-distributed): each rank passes only its block-assigned chunk of the C-flattened grid values, shape (block_len,) or (block_len, rank). The blocks are concatenated (one allgather) into the full field before the solve. The per-rank split must match interpolate_bspline_distributed() (contiguous blocks of the C-order flattened grid, the first n_total % n_parts blocks one element longer).

Parameters:

values (npt.ArrayLike) – Sample values. Full tensor-product field when values_distributed is False; this rank’s flattened block when True.
nodes (PointsLattice | npt.NDArray[np.floating[Any]] | Sequence[npt.NDArray[np.floating[Any]]]) –
Tensor-product interpolation nodes (identical on every rank).
- A PointsLattice: tensor-product grid.
- A 1D ndarray: 1D tensor-product (single direction).
- A sequence of 1D ndarray values: N-D tensor-product.
distributed_space (DistributedSpace) – The distributed space to fit onto. Its global_space must be a BsplineSpace.
values_distributed (bool) – Whether values is this rank’s flattened block (True) or the full replicated field (False). Defaults to False.
tol (float | None) – SVD truncation tolerance. If None, defaults to 100 * machine_epsilon. Defaults to None.

Returns:

A distributed function whose global_function holds the full assembled global coefficient field (identical on every rank), and whose local reproduces the fit over the rank’s owned cells.

Return type:

Raises:

TypeError – If distributed_space.global_space is not a BsplineSpace.
ValueError – If nodes are scattered (a 2D ndarray) while values_distributed is True (scattered points have no tensor-product block layout), or if nodes / values are inconsistent with the space.

Note

Scattered nodes (a 2D ndarray of shape (n_pts, dim)) are only supported in the replicated path (values_distributed=False); they fall back to the serial fit_bspline() solve run identically on every rank (correct, but with no parallel speedup). The output dtype follows global_space.dtype.

Example

>>> from mpi4py import MPI  
>>> from pantr.bspline import create_greville_lattice  
>>> from pantr.bspline import create_uniform_space  
>>> from pantr.mpi import create_distributed_space  
>>> from pantr.mpi import fit_bspline_distributed  
>>> space = create_uniform_space([2, 2], [8, 8])  
>>> ds = create_distributed_space(space, MPI.COMM_WORLD)  
>>> lat = create_greville_lattice(space)  
>>> vals = (lat.get_all_points()[:, 0] ** 2)  # full field on every rank  
>>> dfn = fit_bspline_distributed(vals, lat, ds)  

pantr.mpi.from_dolfinx(mesh, n_cells, *, dolfinx_to_pantr=None)[source]¶

Build a cell Partition from a dolfinx mesh.

Reads the locally-owned cells on every rank and MPI-allgathers their original (input) global indices over mesh.comm to assemble one global, redundant per-cell owner array – the rank that owns each cell in dolfinx. The result is identical on every rank. Pantr cells absent from the mesh (e.g. exterior / trimmed cells excluded by an immersed consumer) get owner -1.

The correspondence between a dolfinx cell and a pantr grid cell is its original cell index (mesh.topology.original_cell_index): with dolfinx_to_pantr=None the mesh is assumed to have been built in pantr’s C-order cell numbering, so the original index is the pantr flat cell id. Pass dolfinx_to_pantr to map the original dolfinx global cell index to a pantr cell id explicitly.

Parameters:

mesh (Any) – A dolfinx.mesh.Mesh. The required interface is duck-typed (no dolfinx import): mesh.comm (an MPI communicator with size and allgather), mesh.topology.dim, mesh.topology.index_map(dim).size_local, and mesh.topology.original_cell_index (must support integer-slice indexing).
n_cells (int) – Total number of cells in the pantr grid (e.g. grid.num_cells); must be >= 1.
dolfinx_to_pantr (npt.ArrayLike | None) – Optional integer map from a dolfinx original global cell index to a pantr cell id. None means identity (the mesh is in pantr C-order).

Returns:

A per-cell owner assignment with mesh.comm.size parts; -1 for pantr cells absent from the mesh.

Return type:

Partition

Raises:

ValueError – If n_cells < 1; if topology.index_map(dim).size_local is negative; if dolfinx_to_pantr is not a 1D integer array, is not injective, or an original index falls outside it; if any (original or mapped) pantr cell id is outside [0, n_cells); if allgather returns an inconsistent number of entries or a non-1D result for any rank; if any gathered cell id is out of range; or if two ranks own the same cell (inconsistent mesh/map).

Note

Calling this engages the default MPI thread policy: unless threads were explicitly configured, this process’s Numba thread pool is limited to one thread per rank. See pantr.mpi.configure_threads().

pantr.mpi.interpolate_bspline_distributed(func, distributed_space, *, nodes=None, boundary_derivatives=None, tol=None)[source]¶

Interpolate a callable onto a distributed tensor-product B-spline space.

The MPI-parallel counterpart of interpolate_bspline(). The collocation solve is a cheap sequence of small 1D SVD solves; the cost that motivates parallelism is evaluating func on the tensor-product node grid (whose size is the global DOF count). Each rank evaluates func on its block-assigned chunk of the flattened grid, a single allgather assembles the full value field, and the Kronecker solve runs replicated on every rank. The returned DistributedFunction carries the global coefficient field (identical on every rank, equal to the serial result) and its local reproduces the serial interpolant over the rank’s owned cells.

Unlike serial interpolate_bspline(), whose func receives a PointsLattice, here func receives a flat (M, dim) point array – a contiguous chunk of the grid is not itself a tensor product, so a flat array is the natural distributed convention (and matches the distributed quasi-interpolants).

Construction requires one MPI collective (comm.allgather) after each rank’s local evaluation.

In (M, rank) shapes below, rank is func’s vector output dimension (the number of value components), not the MPI rank.

Parameters:

func (Callable) – Function to interpolate. Called on a flat (M, dim) point array; must return (M,) (scalar) or (M, rank) (vector-valued). On a rank with no assigned points it is called on an empty (0, dim) block and must return a correspondingly empty (0,) / (0, rank) result.
distributed_space (DistributedSpace) – The distributed space to interpolate onto. Its global_space must be a BsplineSpace.
nodes (Literal['greville'] | PointsLattice | npt.NDArray[np.floating[Any]] | Sequence[npt.NDArray[np.floating[Any]]] | None) –
Interpolation node selection, identical to interpolate_bspline().
- None or "greville" (default): Greville abscissae.
- A PointsLattice: custom tensor-product grid.
- A 1D ndarray: custom nodes for a 1D space.
- A sequence of 1D ndarray values: per-direction custom nodes.
boundary_derivatives (Sequence[tuple[int, ...] | None] | None) – Per-direction (n_left, n_right) boundary-derivative constraints (set to zero), or None. Defaults to None.
tol (float | None) – SVD truncation tolerance for the collocation solve. If None, defaults to 100 * machine_epsilon. Defaults to None.

Returns:

A distributed function whose global_function holds the full assembled global coefficient field (identical on every rank), and whose local interpolates func over this rank’s owned cells.

Return type:

Raises:

TypeError – If distributed_space.global_space is not a BsplineSpace.
ValueError – If nodes is inconsistent with the space, or func returns an output with an invalid shape.

Note

Tensor-product (lattice / Greville) nodes are the parallelized path. All internal computation uses float64; the global control points are cast to global_space.dtype before assembly, consistent with serial interpolate_bspline().

Example

>>> from mpi4py import MPI  
>>> import numpy as np  
>>> from pantr.bspline import create_uniform_space  
>>> from pantr.mpi import create_distributed_space  
>>> from pantr.mpi import interpolate_bspline_distributed  
>>> space = create_uniform_space([2, 2], [8, 8])  
>>> ds = create_distributed_space(space, MPI.COMM_WORLD)  
>>> dfn = interpolate_bspline_distributed(  
...     lambda p: np.sin(p[:, 0]) * np.cos(p[:, 1]), ds
... )

pantr.mpi.l2_project_bspline_distributed(func, distributed_space, *, n_quad=None, quadrature='gauss-legendre', boundary_interpolation=False, tol=None)[source]¶

L2-project a callable onto a distributed tensor-product B-spline space.

The MPI-parallel counterpart of l2_project_bspline(). L2 assembly is per-element and maps directly onto the cell partition: each rank evaluates func only on the quadrature points of its owned cells and contracts them into a per-component load tensor; a single allreduce sums the global load across ranks; and the replicated Kronecker solve recovers the global coefficients. The returned DistributedFunction agrees with the serial L2 projection pointwise over every owned cell.

The per-direction mass matrices are partition-independent and built identically on every rank (cheap n_dofs_i x n_dofs_i systems), so only the load is communicated. boundary_interpolation rows are handled after the reduce: each boundary trace is partition-independent, so every rank recomputes the same boundary row (the boundary trace spans the whole face, so it is not attributable to a single owning cell – every rank recomputes it identically from the global lattice). Construction requires one MPI collective (comm.allreduce) after the local assembly.

Parameters:

func (Callable[..., npt.ArrayLike]) – Function to project. Called as func(lattice) where lattice is a PointsLattice of quadrature points (the serial convention); must return an array of shape (n_total,) for scalar or (n_total, rank) for vector-valued functions.
distributed_space (DistributedSpace) – The distributed space to project onto. Its global_space must be a BsplineSpace.
n_quad (int | Sequence[int] | None) – Quadrature points per element per direction. Defaults to degree + 1.
quadrature (Literal["gauss-legendre", "gauss-lobatto"]) – Quadrature rule type. Defaults to "gauss-legendre".
boundary_interpolation (bool | Sequence[tuple[bool, bool]]) – Replace boundary rows with interpolation conditions. False (default) is a pure L2 projection, True interpolates at all non-periodic boundaries, and a sequence of (left, right) pairs sets per-direction flags.
tol (float | None) – SVD truncation tolerance for the per-direction solves. If None, defaults to 100 * machine_epsilon.

Returns:

A distributed function whose local L2-projects func over this rank’s owned cells, and whose global_function holds the full assembled global coefficient field (identical on every rank after the allreduce).

Return type:

Raises:

TypeError – If distributed_space.global_space is not a BsplineSpace.
ValueError – If n_quad or boundary_interpolation is inconsistent with the global space, if func returns an output with an invalid shape, or if the reduced load does not match the expected stacked shape (*num_basis, n_components) (a symptom of func returning inconsistent shapes across ranks).

Note

func MUST be rank-independent: for a given quadrature lattice it must return the same shape and dtype on every rank. The reduction (comm.allreduce) is a collective that every rank must reach with a matching contribution; if func raised on a subset of ranks (e.g. a shape error seen by some ranks only) those ranks would abort before the collective and deadlock the rest, and mixed dtypes would corrupt the reduction.

Note

The output dtype is inferred from the return value of func (as in the serial l2_project_bspline()). Unlike the serial path – where Bspline raises on a control-point/space dtype mismatch – the distributed path is more lenient: it coerces the assembled global control points to global_space.dtype rather than raising.

Example

>>> from mpi4py import MPI  
>>> import numpy as np  
>>> from pantr.bspline import create_uniform_space  
>>> from pantr.mpi import create_distributed_space  
>>> from pantr.mpi import l2_project_bspline_distributed  
>>> space = create_uniform_space([2, 2], [8, 8])  
>>> ds = create_distributed_space(space, MPI.COMM_WORLD)  
>>> dfn = l2_project_bspline_distributed(  
...     lambda lat: np.sin(lat.pts_per_dir[0]), ds
... )
>>> local = dfn.local  # rank-local Bspline on the windowed space  

pantr.mpi.mpi_available()[source]¶

Report whether mpi4py can be imported in this environment.

Uses importlib.util.find_spec(), so it never triggers an actual import or MPI initialization.

Returns:: True if mpi4py is installed and importable, False otherwise.
Return type:: bool

pantr.mpi.quasi_interpolate_bspline_distributed(func, distributed_space, *, kind='llm')[source]¶

Quasi-interpolate a callable onto a distributed tensor-product B-spline space.

The MPI-parallel counterpart of quasi_interpolate_bspline(). Each rank evaluates func only on the interior points required by its owned basis functions (Lee-Lyche-Mørken functionals); a single allgather at the end assembles the global coefficient field. The returned DistributedFunction agrees with the serial quasi-interpolant pointwise over every owned cell.

Construction requires one MPI collective (comm.allgather) after the local computation. Per-rank func evaluation is purely local: no rank ever evaluates func outside its windowed parametric sub-domain.

Parameters:

func (Callable) – Function to quasi-interpolate. Called on a flat (M, dim) point array; must return (M,) (scalar) or (M, rank) (vector-valued).
distributed_space (DistributedSpace) – The distributed space to interpolate onto. Its global_space must be a BsplineSpace.
kind (QIKind) – Quasi-interpolant kind. Only "llm" (Lee-Lyche-Mørken) is currently supported. Defaults to "llm".

Returns:

A distributed function whose local quasi-interpolates func over this rank’s owned cells, and whose global_function holds the full assembled global coefficient field (identical on every rank after the allgather).

Return type:

Raises:

TypeError – If distributed_space.global_space is not a BsplineSpace.
ValueError – If kind is not recognized, or if func returns an output with an invalid shape (0-D, more than 2-D, or wrong leading dimension).

Note

All internal computation uses float64. The global control points are cast to global_space.dtype before assembly, consistent with the serial quasi_interpolate_bspline().

Example

>>> from mpi4py import MPI  
>>> import numpy as np  
>>> from pantr.bspline import create_uniform_space  
>>> from pantr.mpi import create_distributed_space  
>>> from pantr.mpi import quasi_interpolate_bspline_distributed  
>>> space = create_uniform_space([2, 2], [8, 8])  
>>> ds = create_distributed_space(space, MPI.COMM_WORLD)  
>>> dfn = quasi_interpolate_bspline_distributed(  
...     lambda p: np.sin(p[:, 0]) * np.cos(p[:, 1]), ds
... )
>>> local = dfn.local  # rank-local Bspline on the windowed space  

pantr.mpi.quasi_interpolate_thb_spline_distributed(func, distributed_space, *, kind='llm')[source]¶

Quasi-interpolate a callable onto a distributed THB-spline space.

The MPI-parallel counterpart of quasi_interpolate_thb_spline(). Each rank runs the serial Speleers-Manni hierarchical quasi-interpolant on its windowed local space and keeps only its owned DOFs’ coefficients; a single allgather at the end assembles the global coefficient field. The returned DistributedFunction agrees with the serial quasi-interpolant pointwise over every owned cell.

The windowed local space covers the support closure of every owned DOF, so the per-level functional’s level-l active leaf cell – the one the serial routine samples func on – always lies inside the rank’s windowed parametric domain. No rank ever evaluates func outside that domain. Construction requires one MPI collective (comm.allgather) after the local computation.

Parameters:

func (Callable) – Function to quasi-interpolate. Called on a flat (M, dim) point array; must return (M,) (scalar) or (M, rank) (vector-valued).
distributed_space (DistributedSpace) – The distributed space to interpolate onto. Its global_space must be a THBSplineSpace.
kind (QIKind) – Quasi-interpolant kind. Only "llm" (Lee-Lyche-Mørken) is currently supported. Defaults to "llm".

Returns:

Return type: