KnotVector(
    n_basis_functions::Integer, 
    degree::Integer; 
    extent::Tuple{Number, Number} = (0,1), 
    distribution::Symbol = :equispaced)

Construct a clamped knot vector, i.e. the multiplicity of the first and last knot is degree + 1 and the other multiplicities are 1.

Arguments

• n_basis_functions: The number of basis functions that will be defined on this knot vector
• degree: The degree of the basis functions that will be defined on this knot vector

Keyword Arguments

• extent: A tuple (tmin, tmax) defining the extend of the knot vector
• distribution: The distribution of the internal knots. The options are :equispaced or :random
• backend: The KernelAbstractions backend of the arrays in the object. Defaults to CPU().
• float_type: The type of all floating point arrays. Defaults to Float32.
• int_type: The type of all integer arrays. Defaults to Int32.

SplineDimension(
    n_basis_functions::Integer,
    degree::Integer,
    n_sample_points::Integer;
    max_derivative_order::Integer = 0,
    knot_vector::Union{Nothing, KnotVector{Tv, Ti}} = nothing,
    backend::Backend = CPU(),
    float_type::Type{Tv} = Float32,
    int_type::Type{Ti} = Int32,
    kwargs...)::SplineDimension where {Tv <: AbstractFloat, Ti <: Integer}

Constructor for a SplineDimension. Optionally a knot_vector kwarg can be passed, otherwise a default knot vector is generated. For now by default the sample points are evenly spaced on the extent of the knot vector. Key word arguments are passed to the KnotVector constructor.

Inputs

• n_basis_functions: The number of basis functions of this spline dimension
• degree: The degree of the basis functions of this spline dimension
• n_sample_points: The number of points at which the domain of the basis functions will be sampled
• max_derivative_order: The maximum derivative order of the basis functions that will be computed in evaluate!. Defaults to 0. knot_vector: A knot vector on which the basis functions will be defined. Defaults to nothing, which means that a default clamped/open equally spaced knot vector will be defined.
• backend: The KernelAbstractions backend of the arrays in the object. Defaults to CPU(). NOTE: If a knot vector is supplied, its backend takes precedence.
• float_type: The type of all floating point arrays. Defaults to Float32.

SplineGrid(spline_dimensions::NTuple{Nin, <:SplineDimension{Tv, Ti}}, Nout::Integer)::SplineGrid{Nin, Tv, Ti} where {Nin, Tv, Ti}

Define a SplineGrid from an NTuple of spline dimensions and the number of output dimensions.

Inputs

• spline_dimensions: an NTuple of spline dimensions
• Nout: The number of output dimensions. I.e. the control points and thus the spline live in ℝ^Nout.

Create a SplineGrid but with preallocated weights to define a NURBS. See The SplineGrid constructors for more details.

rmeye(
    n::Integer;
    backend::Backend = CPU(),
    float_type::Type{Tv} = Float32,
    int_type::Type{Ti} = Int)::RefinementMatrix{Tv, Ti} where {Tv, Ti <: Integer}

Construct an identity refinement matrix.

Input

• n: The size of the identity matrix is n×n
• backend: The KernelAbstractions backend of the matrix data
• float_type: The value type of the matrix data
• int_type: The integer type of the matrix data

RefinementMatrix(
    spline_dimension::SplineDimension{Tv, Ti},
    knot_span_index::Integer,
    knot_new
    )::RefinementMatrix{Tv, Ti} where {Tv, Ti <: Integer}

Build a refinement matrix per row given a knot vector and a to be added new knot.

Inputs

• spline_dimension: The spline dimension whose knot vector the new knot is to be added to
• knot_span_index: The index such that knots_all[knot_span_index] ≤ knot_new < knots_all[knot_span_index + 1]
• knot_new: The value of the new knot

evaluate!(spline_dimension)

Per sample point, get the value of the spline_dimension.degree + 1 basis functions that have a non-zero value for that sample point. This is based on the Cox-de Boor recursion formula.

The l-th sample point t has sample index i, meaning that t ∈ [tᵢ, tᵢ₊₁). Therefore Bᵢ₀(t) = 1, Bⱼ₀(t) = 0 for j ≠ i. For degree k, t is in the domain of Bⱼₖ which is [tⱼ, tⱼ₊ₖ₊₁), for j = i - k, ..., i.

Arguments

• spline_dimension

evaluate!(spline_grid::SplineGrid{Nin};
    derivative_order::NTuple{Nin, <:Integer} = ntuple(_ -> 0, Nin),
    control_points::AbstractArray = spline_grid.control_points,
    eval::AbstractArray = spline_grid.eval)

Evaluate the spline grid, that is: take the evaluated basis functions for each sample point for each SplineDimension, and compute the output grid on each sample point combination as a linear combination of control points with basis function products as coefficients.

If weights are supplied, compute the rational basis functions for NURBS as the control point coefficients.

Uses the control_points and eval arrays from the spline_grid by default, but different arrays can be specified as a convenience for optimization algorithms.

evaluate!(spline_grid::SplineGrid{Nin};
    derivative_order::NTuple{Nin, <:Integer} = ntuple(_ -> 0, Nin),
    control_points::AbstractArray = spline_grid.control_points,
    eval::AbstractArray = spline_grid.eval)

Evaluate the spline grid, that is: take the evaluated basis functions for each sample point for each SplineDimension, and compute the output grid on each sample point combination as a linear combination of control points with basis function products as coefficients.

If weights are supplied, compute the rational basis functions for NURBS as the control point coefficients.

Uses the control_points and eval arrays from the spline_grid by default, but different arrays can be specified as a convenience for optimization algorithms.

evaluate_adjoint!(spline_grid::AbstractSplineGrid{Nin, Nout, false, Tv};
    derivative_order::NTuple{Nin, <:Integer} = ntuple(_ -> 0, Nin),
    control_points::AbstractControlPointArray{Nin, Nout, Tv} = spline_grid.control_points,
    eval::AbstractArray = spline_grid.eval)::Nothing where {Nin, Nout, Tv}

evaluate the adjoint of the linear mapping control_points -> eval. This is a computation of the form eval -> control_points. If we write evaluate!(spline_grid) as a matrix vector multiplication eval = M * control_points, Then the adjoint is given by v -> M' * v. This mapping is used in fitting algorithms.

mult!(
    Y::V,
    As::NTuple{<:RefinementMatrix{Tv}, N},
    B::V,
    dims_refinement::NTuple{<:Integer, N})::Nothing where {Tv, V <: AbstractArray{Tv}, N}

'left-multiply' the array B by every refinement matrix along the specified dimension.

Inputs

• Y: The target array of the multiplication
• As: The refinement matrices Y will be multiplied by
• B: The array that will be multiplied by the refinement matrices
• dims_refinement: The dimension along which the refinement matrix multiplication will take place for each matrix

Evaluate the locally refined control points. For each local refinement, first apply the refinement matrix and then overwrite the desired control point values.

insert_knot(
knot_vector::KnotVector, 
knot_new::AbstractFloat)::Tuple{KnotVector, Integer}

Create a new knot vector with the new knot of multiplicity 1.

Inputs

• knot_vector: The knot vector object to which the knot will be added
• knot_new: The value of the new knot. Should not be part of the knot values already

Outputs

• knot_vector_new: The newly created knot vector with the added knot
• knot_span_index: The index such that knots_all[knot_span_index] ≤ knot_new < knots_all[knot_span_index + 1]

insert_knot(
    spline_dimension::SplineDimension{Tv, Ti},
    knot_new::AbstractFloat;
    recompute_sample_indices::Bool = true,
    evaluate::Bool = true)::Tuple{
        SplineDimension{Tv, Ti},
        RefinementMatrix{Tv, Ti}
    } where {Tv, Ti}

Create a new spline dimension whose knot vector has the new knot with multiplicity 1.

Inputs

• spline_dimension: The spline dimension the new knot will be added to
• knot_new: The value of the new knot
• recompute_sample_indices: Whether the indices of the sample points should be recomputed after the knot insertion. Defaults to true.
• evaluate: Whether the spline dimension should be evaluated after the knot insertion. Defaults to true.

Outputs

• spline_dimension_new: The newly created spline dimension with the same underlying memory except for the new knot vector.
• refinement_matrix: The sparse matrix which expresses the basis functions from before the knot insertion in terms of the basis functions after the knot insertion.

insert_knot(
    spline_grid::A,
    dim_refinement::Integer,
    knot_new::AbstractFloat;
    evaluate_spline_dimension::Bool = true)::Tuple{A, RefinementMatrix} where {A <: AbstractSplineGrid}

Create a new spline grid where a new knot is added to the knot vector underlying the indicated spline dimension.

Inputs

• spline_grid: The spline grid to which the new knot will be added
• knot_new: The value of the knot to be added
• dim_refinement: The index of the spline dimension to which the knot will be added
• evaluate_spline_dimension: Whether the spline dimension to which the knot is added should be evaluated. Defaults to true.

Outputs

• spline_grid_new: The newly created spline grid with all the same underlying memory except for the updated knot vector and the control points.
• refinement_matrix: The sparse matrix which expresses the basis functions from before the knot insertion in terms of the basis functions after the knot insertion for the knot insertion dimension.

refine(
    spline_dimension::SplineDimension{Tv, Ti};
    knots_new::Union{Vector{<:AbstractFloat}, Nothing} = nothing)::Tuple{SplineDimension{Tv, Ti}, RefinementMatrix{Tv, Ti}} where {Tv, Ti}

Create a new spline dimension with multiple knots added to the underlying knot vector.

Inputs

• spline_dimension: The spline dimension the new knots will be added to
• knots_new: The vector of knots that will be added. Defaults to the midpoints of the knot spans of the vector underlying the spline dimension.

Outputs

• spline_dimension_new: The newly created spline dimension with the same underlying memory except for the new knot vector.
• refinement_matrix: The sparse matrix which expresses the basis functions from before the refinement in terms of the basis functions after the refinement.

refine(
    spline_grid::A,
    dim_refinement::Integer;
    knots_new::Union{AbstractVector{<:AbstractFloat}, Nothing} = nothing)::Tuple{A, RefinementMatrix} where {A <: AbstractSplineGrid}

Create a new spline grid where multiple knots are added to the knot vector underlying the indicated spline dimension.

Inputs

• spline_grid: The spline grid from which one of the knot vectors will be refined
• dim_refinement: The index of the spline dimension whose knot vector will be refined
• knots_new: The knots that will be added. Defaults to nothing, which internally is translated to all midpoints of the non-trivial knot spans of the knot vector that will be refined.

Outputs

• spline_grid_new: The newly created spline grid with all the same underlying memory except for the updated knot vector and the control points.
• refinement_matrix: The sparse matrix which expresses the basis functions from before the knot insertion in terms of the basis functions after the knot insertion for the knot insertion dimension.

refine(
    spline_grid::AbstractSplineGrid{Nin, Nout, false, Tv, Ti},
    spline_dimension_new::SplineDimension{Tv, Ti},
    dim_refinement::Integer,
    refinement_matrix::RefinementMatrix{Tv, Ti} where {Nin, Nout, Tv, Ti}

Update the spline grid with a refined spline dimension in the specified dimension with the associated refinement matrix.

Given a refinement step (by default the last refinement), add new control points that overwrite the result from the refinement matrix. These new values are chosen such that the spline geometry does not change.

Inputs

• control_points: The locally refined control points where new control points will be activated
• refinement_indices: The indices in the control point grid at the refinement_index level which will be overwritten
• refinement_index: The index of the refinement after which new control points will be activated

Convenience wrapper of activate_local_refinement! for a spline grid

Instead of supplying explicit indices of the control points for local refinement activation, supply a range of indices per dimension.

error_informed_local_refinement!(
    spline_grid::AbstractSplineGrid{Nin, Nout, HasWeights, Tv, Ti},
    error::AbstractArray;
    threshold_factor::Number = 1.0
    )::Nothing where {Nin, Nout, HasWeights, Tv, Ti}

Refine the last level of the locally refined spline grid informed by the error array which has the same shape as spline_grid.eval. This is done by:

• mapping the error back onto the control points by using the adjoint of the refinement matrices multiplication
• summing over the output dimensions to obtain a single number per control point stored in control_grid_error
• activating each control point whose value is bigger than threshold = threshold_factor * mean(control_grid_error)

Perform deactivate_control_points for every refinement level except the last one in reverse order. For more details see deactivate_overwritten_control_points!(::LocallyRefinedControlPoints, ::Integer).

deactivate_overwritten_control_points!(
    control_points::LocallyRefinedControlPoints,
    local_refinement_level::Integer)::Nothing

Deactivate control points whose effect is completely overwritten. The procedure works as follows:

The 'forward' computation to process local refinement is as follows:

B = (O₂ ∘ L ∘ O₁)(A),

where:

• A is the control point grid at local_refinement_level
• B is the control point grid at local_refinement_level + 1
• O₁ is the overwriting operation of the active control points at local_refinement_level
• L is the linear operation consisting of multiplying by refinement matrices along specified dimensions
• O₂ is the overwriting operation of the active control points at local_refinement_level + 1

To figure out whether for any overwriting element of O₁ its effect is still present in the final array, we perform the following computation:

A = (L* ∘ O₂)(B),

where:

• B is initialized with a special number Flag, which is a simple wrapper of a Boolean effectively saying whether this number is relevant or not. B is initialized with all Flag(true)
• The overwriting values of O₂ are initialized with all Flag(false)
• L* is the adjoint version of L

After this computation we can read of in A at the overwriting locations of O₁ whether each number is still having an effect on B or not.

KnotVector(knots, multiplicities)

Defines a knot vector.

Arguments

• knot_values: The values in the knot vector. Must be strictly increasing.
• multiplicities: The multiplicity of each knot in knots.

SplineDimension(
    degree, 
    max_derivative_order, 
    knot_vector, 
    sample_points, 
    sample_indices, 
    eval, 
    eval_prev)

Defines the set of basis functions for a single dimension, and how it is sampled.

Arguments

• degree: The degree of the piecewise polynomial basis functions.
• max_derivative_order: The maximum derivative order of the basis functions that will be computed.
• knot_vector: The knot vector on which the basis functions are defined.
• sample_points: The points in the domain of the basis functions where they are sampled. Must
• lie within the boundaries of the knot vector.
• sample_indices: The indices i of the sample points t in the knot vector such that knot_vector.knots[i] ≤ t < knot_vector.knots[i + 1]`
• eval: An array of shape (length(sample_points), degree + 1, max_derivative + 1), with per sample point the values of those basis functions
• eval_prev: Helper array for intermediate results in the basis function computations whose support the sample point is in, and the derivatives if requested.

RefinementMatrix(m, n, row_pointer, column_start, nzval)

The refinement matrix is a sparse matrix encoding for matrices which:

1. Have consecutive non-zeros in all rows and columns
2. Have at least 1 nonzero in every row and column

The non-zeros are stored in a dense vector per row.

Fields

• m: The number of rows of the matrix
• n: The number of columns of the matrix
• row_pointer: row_pointer[i] indicates where in nzval the data for the i-th row starts
• column_start: column_start[i] indicates at which column the first nonzero for the i-th row is
• nzval: The nonzero values in the matrix

DefaultControlPoints(control_points)

A thin wrapper around an array of control point values.

LocalRefinement(
    dims_refinement,
    refinement_matrices,
    refinement_indices,
    refinement_values)

The data for performing a local refinement; one or more refinement matrices and the to be overwritten control point values after multiplication with the refinement matrices along the specified dimensions.

Fields

• dims_refinement: The dimensions along which will be refined
• refinement_matrices: The matrices that perform the refinement
• refinement_indices: The indices of the control points that will be overwritten after the refinement
• refinement_values: The new values of the control points that will be rewritten after refinement

LocallyRefinedControlPoints(
    control_points_base,
    control_points_refined,
    local_refinements)

All data required to perform multiple refinement steps and overwriting control point values along the way, yielding a Truncated Hierarchical Basis (THB) spline.

Fields

• control_points_base: The densely defined control points at the basis of the hierarchy
• control_points_refined: The intermediate and final control point arrays after applying the local refinements
• local_refinements: Data for local refinement for each step in the hierarchy

SplineGrid(
    spline_dimensions,
    control_points,
    weights,
    eval)

The SplineGrid is the central object of the SplineGrids.jl package, containing all information to evaluate the defined spline on the defined grid.

Fields

• spline_dimensions: A SplineDimension per dimension of the spline, containing data to evaluate basis functions.
• control points: The points that define the shape of the spline, and in how many dimensions it is embedded.
• weights: Control point weights to define NURBS.
• eval: The array where the evaluated spline grid is stored.

decompress(
    spline_dimension::SplineDimension{Tv}; derivative_order::Integer = 0) where {Tv}

Transform spline_dimension.eval into a matrix of shape (n_sample_points, n_points - degree - 1) which explicitly gives the value for each basis function at each sample point.

Get the number of control points within the control point object.

add_default_local_refinement(spline_grid)

Refine in the default way in every dimension, i.e. bisecting every non-trivial knot span. Yields a spline grid with a LocallyRefinedControlPoints object for the control_points field.