Repository: hyrodium/BasicBSpline.jl
Pull request #107 (2021-12-22 17:05:47)

cleanup


Patch:
diff --git a/Makefile b/Makefile
deleted file mode 100644
index 5201dd92..00000000
--- a/Makefile
+++ /dev/null
@@ -1,3 +0,0 @@
-preview:
-	julia --project=docs -e 'using Pkg;Pkg.develop(PackageSpec(path=pwd()));Pkg.instantiate();include("docs/make.jl");'
-	julia -e 'using LiveServer; serve(dir="docs/build")'
diff --git a/src/BasicBSpline.jl b/src/BasicBSpline.jl
index 245e1175..2f9e040e 100644
--- a/src/BasicBSpline.jl
+++ b/src/BasicBSpline.jl
@@ -7,7 +7,7 @@ using GeometryBasics
 
 export Point
 export Knots, ♯, BSplineSpace, dim, ⊑, ⊒, ⋢, ⋣, ≃
-export BSplineManifold
+export BSplineManifold, CustomBSplineManifold
 export bsplinebasis₊₀, bsplinebasis₋₀, bsplinebasis
 export bsplinebasis′₊₀, bsplinebasis′₋₀, bsplinebasis′
 export bsplinesupport, bsplineunity, changebasis, lower
@@ -16,7 +16,6 @@ export refinement, bsplinespaces, controlpoints
 export isproper, properdim
 export degree, knots
 export 𝔫
-export FastBSplineSpace, FastBSplineManifold
 export BSplineCurve, BSplineSurface, BSplineSolid
 export AbstractBSplineManifold, AbstractBSplineSpace
 export fittingcontrolpoints
@@ -29,9 +28,6 @@ include("_BSplineBasis.jl")
 include("_ChangeBasis.jl")
 include("_BSplineManifold.jl")
 include("_CustomBSplineManifold.jl")
-include("_FastBSplineSpace.jl")  # TODO: remove this
-include("_FastBSplineBasis.jl")  # TODO: remove this
-include("_FastBSplineManifold.jl")  # TODO: remove this
 include("_Refinement.jl")
 include("_Integral.jl")
 include("_Fitting.jl")
diff --git a/src/_BSplineSpace.jl b/src/_BSplineSpace.jl
index da4ee59b..d20c2c0a 100644
--- a/src/_BSplineSpace.jl
+++ b/src/_BSplineSpace.jl
@@ -11,12 +11,12 @@ This type `BSplineSpace` is slower than `FastBSplineSpace`, but this type is not
 """
 struct BSplineSpace{p, T<:Real} <: AbstractBSplineSpace{p,T}
     knots::Knots{T}
-    global unsafe_bsplinespace(::Type{Val{p}}, k::Knots{T}) where {p,T} = new{p,T}(k)
+    global unsafe_bsplinespace(::Val{p}, k::Knots{T}) where {p,T} = new{p,T}(k)
 end
 function BSplineSpace{p}(k::Knots) where p
     # TOOD: add error handling like this:
     # throw(DomainError(p, "degree of polynominal must be non-negative"))
-    unsafe_bsplinespace(Val{p}, k)
+    unsafe_bsplinespace(Val{p}(), k)
 end
 
 """
diff --git a/src/_FastBSplineBasis.jl b/src/_FastBSplineBasis.jl
deleted file mode 100644
index e1949302..00000000
--- a/src/_FastBSplineBasis.jl
+++ /dev/null
@@ -1,227 +0,0 @@
-# Faster B-spline basis function
-
-function _s(c::Char, i::Int)
-    string(c) * "_" * string(i)
-end
-
-function _code_K(p, q)
-    join([_s('K', j) for j in 1:p-q+2], ',') *
-    " = " *
-    join(["(t-" * _s('k', j) * ")÷(" * _s('k', j + q) * "-" * _s('k', j) * ")" for j in 1:p-q+2], ',')
-end
-
-function _code_K′(p)
-    join([_s('K', j) for j in 1:2], ',') *
-    " = " *
-    join(["1÷(" * _s('k', j + p) * "-" * _s('k', j) * ")" for j in 1:2], ',')
-end
-
-function _code_B(p, q)
-    join([_s('B', j) for j in 1:p-q+1], ',') *
-    " = " *
-    join([_s('K', j) * "*" * _s('B', j) * "+(1-" * _s('K', j + 1) * ")*" * _s('B', j + 1) for j in 1:p-q+1], ',')
-end
-
-for p in 0:MAX_DEGREE
-    @eval function bsplinebasis₊₀(P::FastBSplineSpace{$p}, i::Integer, t::Real)
-        ÷(a, b) = ifelse(b == 0.0, 0.0, a / b)
-        $(Meta.parse(join([_s('k', j) for j in 1:p+2], ',') * " = " * join(["P.knots[i+" * string(j) * "]" for j in 0:p+1], ',')))
-        $(Meta.parse(
-            join([_s('B', j) for j in 1:p+1], ',') * " = " * join(["Float64(" * _s('k', j) * " ≤ t < " * _s('k', j + 1) * ")" for j in 1:p+1], ','),
-        ))
-        $(Meta.parse("begin " * join([_code_K(p, q) * "\n" * _code_B(p, q) for q in 1:p], '\n') * " end"))
-        return B_1
-    end
-    @eval function bsplinebasis₋₀(P::FastBSplineSpace{$p}, i::Integer, t::Real)
-        ÷(a, b) = ifelse(b == 0.0, 0.0, a / b)
-        $(Meta.parse(join([_s('k', j) for j in 1:p+2], ',') * " = " * join(["P.knots[i+" * string(j) * "]" for j in 0:p+1], ',')))
-        $(Meta.parse(
-            join([_s('B', j) for j in 1:p+1], ',') * " = " * join(["Float64(" * _s('k', j) * " < t ≤ " * _s('k', j + 1) * ")" for j in 1:p+1], ','),
-        ))
-        $(Meta.parse("begin " * join([_code_K(p, q) * "\n" * _code_B(p, q) for q in 1:p], '\n') * " end"))
-        return B_1
-    end
-    @eval function bsplinebasis(P::FastBSplineSpace{$p}, i::Integer, t::Real)
-        ÷(a, b) = ifelse(b == 0.0, 0.0, a / b)
-        k_end = P.knots[end]
-        $(Meta.parse(join([_s('k', j) for j in 1:p+2], ',') * " = " * join(["P.knots[i+" * string(j) * "]" for j in 0:p+1], ',')))
-        $(Meta.parse(
-            join([_s('B', j) for j in 1:p+1], ',') * " = " * join(["Float64(" * _s('k', j) * " ≤ t < " * _s('k', j + 1) * ")" for j in 1:p+1], ','),
-        ))
-        $(Meta.parse(_s('B', p + 1) * " += Float64(t == " * _s('k', p + 2) * " == k_end)"))
-        $(Meta.parse("begin " * join([_code_K(p, q) * "\n" * _code_B(p, q) for q in 1:p], '\n') * " end"))
-        return B_1
-    end
-end
-
-bsplinebasis′₊₀(P::FastBSplineSpace{0}, i::Integer, t::Real) = 0.0
-bsplinebasis′₋₀(P::FastBSplineSpace{0}, i::Integer, t::Real) = 0.0
-bsplinebasis′(P::FastBSplineSpace{0}, i::Integer, t::Real) = 0.0
-
-for p in 1:MAX_DEGREE
-    @eval function bsplinebasis′₊₀(P::FastBSplineSpace{$p}, i::Integer, t::Real)
-        ÷(a, b) = ifelse(b == 0.0, 0.0, a / b)
-        $(Meta.parse(join([_s('k', j) for j in 1:p+2], ',') * " = " * join(["P.knots[i+" * string(j) * "]" for j in 0:p+1], ',')))
-        $(Meta.parse(
-            join([_s('B', j) for j in 1:p+1], ',') * " = " * join(["Float64(" * _s('k', j) * " ≤ t < " * _s('k', j + 1) * ")" for j in 1:p+1], ','),
-        ))
-        $(Meta.parse("begin " * join([_code_K(p, q) * "\n" * _code_B(p, q) for q in 1:p-1], '\n') * " end"))
-        $(Meta.parse(_code_K′(p)))
-        B_1 = $p * (K_1 * B_1 - K_2 * B_2)
-        return B_1
-    end
-    @eval function bsplinebasis′₋₀(P::FastBSplineSpace{$p}, i::Integer, t::Real)
-        ÷(a, b) = ifelse(b == 0.0, 0.0, a / b)
-        $(Meta.parse(join([_s('k', j) for j in 1:p+2], ',') * " = " * join(["P.knots[i+" * string(j) * "]" for j in 0:p+1], ',')))
-        $(Meta.parse(
-            join([_s('B', j) for j in 1:p+1], ',') * " = " * join(["Float64(" * _s('k', j) * " < t ≤ " * _s('k', j + 1) * ")" for j in 1:p+1], ','),
-        ))
-        $(Meta.parse("begin " * join([_code_K(p, q) * "\n" * _code_B(p, q) for q in 1:p-1], '\n') * " end"))
-        $(Meta.parse(_code_K′(p)))
-        B_1 = $p * (K_1 * B_1 - K_2 * B_2)
-        return B_1
-    end
-    # @eval function bsplinebasis′(P::FastBSplineSpace{$p}, i::Integer, t::Real)
-    #     ÷(a,b) = ifelse(b == 0.0, 0.0, a/b)
-    #     k_end = P.knots[end]
-    #     $(Meta.parse(join([_s('k',j) for j in 1:p+2],',')*" = "*join(["P.knots[i+"*string(j)*"]" for j in 0:p+1],',')))
-    #     $(Meta.parse(join([_s('B',j) for j in 1:p+1],',')*" = "*join(["Float64("*_s('k',j)*" ≤ t < "*_s('k',j+1)*")" for j in 1:p+1],',')))
-    #     $(Meta.parse(_s('B',p+1)*" += Float64(t == "*_s('k',p+2)*" == k_end)"))
-    #     $(Meta.parse("begin "*join([_code_K(p,q)*"\n"*_code_B(p,q) for q in 1:p-1],'\n')*" end"))
-    #     $(Meta.parse(_code_K′(p)))
-    #     B_1 = $p*(K_1*B_1-K_2*B_2)
-    #     return B_1
-    # end
-end
-
-function bsplinebasis′(P::FastBSplineSpace, i::Integer, t::Real)::Float64
-    # TODO: faster
-    ÷(a, b) = ifelse(b == 0.0, 0.0, a / b)
-
-    p = degree(P)
-    k = knots(P)
-
-    B_1, B_2 = bsplinebasis(lower(P), i, t), bsplinebasis(lower(P), i+1, t)
-
-    K_1, K_2 = 1 ÷ (k[i+p] - k[i]), 1 ÷ (k[i+p+1] - k[i+1])
-    B_1 = p * (K_1 * B_1 - K_2 * B_2)
-
-    return B_1
-end
-
-# TODO: use @eval macro to produce code
-"""
-Returns the value of ``B_{(i,0,k)}(t)``.
-Assumption:
-* ```k_{i} ≤ t < k_{i+1}``
-"""
-function _bsb0(::Union{Vector{<:Real},Knots}, ::Real, ::Integer)
-    return 1.0
-end
-
-"""
-Returns the values of ``B_{(i,1,k)}(t), B_{(i+1,1,k)}(t)``.
-Assumption:
-* ``k_{i} ≤ t < k_{i+1}``
-"""
-function _bsb1(k::Union{Vector{<:Real},Knots}, t::Real, i::Integer)
-    B1 = (k[i+1]-t)/(k[i+1]-k[i])
-    B2 = (t-k[i])/(k[i+1]-k[i])
-    return B1, B2
-end
-
-"""
-Returns the values of ``B_{(i,2,k)}(t), B_{(i+1,2,k)}(t), B_{(i+2,2,k)}(t)``.
-Assumption:
-* ``k_{i} ≤ t < k_{i+1}``
-"""
-@inline function _bsb2(k::Union{Vector{<:Real},Knots}, t::Real, i::Integer)
-    B = _bsb1(k, t, i)
-
-    B1 = (k[i+1]-t)/(k[i+1]-k[i-1]) * B[1]
-    B2 = (t-k[i-1])/(k[i+1]-k[i-1]) * B[1] +
-         (k[i+2]-t)/(k[i+2]-k[i]) * B[2]
-    B3 = (t-k[i])/(k[i+2]-k[i]) * B[2]
-    return B1, B2, B3
-end
-
-"""
-Returns the values of ``B_{(i,3,k)}(t), B_{(i+1,3,k)}(t), B_{(i+2,3,k)}(t), B_{(i+3,3,k)}(t)``.
-Assumption:
-* ``k_{i} ≤ t < k_{i+1}``
-"""
-@inline function _bsb3(k::Union{Vector{<:Real},Knots}, t::Real, i::Integer)
-    B = _bsb2(k, t, i)
-
-    B1 = (k[i+1]-t)/(k[i+1]-k[i-2]) * B[1]
-    B2 = (t-k[i-2])/(k[i+1]-k[i-2]) * B[1] +
-         (k[i+2]-t)/(k[i+2]-k[i-1]) * B[2]
-    B3 = (t-k[i-1])/(k[i+2]-k[i-1]) * B[2] +
-         (k[i+3]-t)/(k[i+3]-k[i]) * B[3]
-    B4 = (t-k[i])/(k[i+3]-k[i]) * B[3]
-    return B1, B2, B3, B4
-end
-
-"""
-Returns the values of ``B_{(i,4,k)}(t), B_{(i+1,4,k)}(t), B_{(i+2,4,k)}(t), B_{(i+3,4,k)}(t), B_{(i+4,4,k)}(t)``.
-Assumption:
-* ``k_{i} ≤ t < k_{i+1}``
-"""
-@inline function _bsb4(k::Union{Vector{<:Real},Knots}, t::Real, i::Integer)
-    B = _bsb3(k, t, i)
-
-    B1 = (k[i+1]-t)/(k[i+1]-k[i-3]) * B[1]
-    B2 = (t-k[i-3])/(k[i+1]-k[i-3]) * B[1] +
-         (k[i+2]-t)/(k[i+2]-k[i-2]) * B[2]
-    B3 = (t-k[i-2])/(k[i+2]-k[i-2]) * B[2] +
-         (k[i+3]-t)/(k[i+3]-k[i-1]) * B[3]
-    B4 = (t-k[i-1])/(k[i+3]-k[i-1]) * B[3] +
-         (k[i+4]-t)/(k[i+4]-k[i]) * B[4]
-    B5 = (t-k[i])/(k[i+4]-k[i]) * B[4]
-    return B1, B2, B3, B4, B5
-end
-
-"""
-Returns the values of ``B_{(i,5,k)}(t), B_{(i+1,5,k)}(t), B_{(i+2,5,k)}(t), B_{(i+3,5,k)}(t), B_{(i+4,5,k)}(t), B_{(i+5,5,k)}(t)``.
-Assumption:
-* ``k_{i} ≤ t < k_{i+1}``
-"""
-@inline function _bsb5(k::Union{Vector{<:Real},Knots}, t::Real, i::Integer)
-    B = _bsb4(k, t, i)
-
-    B1 = (k[i+1]-t)/(k[i+1]-k[i-4]) * B[1]
-    B2 = (t-k[i-4])/(k[i+1]-k[i-4]) * B[1] +
-         (k[i+2]-t)/(k[i+2]-k[i-3]) * B[2]
-    B3 = (t-k[i-3])/(k[i+2]-k[i-3]) * B[2] +
-         (k[i+3]-t)/(k[i+3]-k[i-2]) * B[3]
-    B4 = (t-k[i-2])/(k[i+3]-k[i-2]) * B[3] +
-         (k[i+4]-t)/(k[i+4]-k[i-1]) * B[4]
-    B5 = (t-k[i-1])/(k[i+4]-k[i-1]) * B[4] +
-         (k[i+5]-t)/(k[i+5]-k[i]) * B[5]
-    B6 = (t-k[i])/(k[i+5]-k[i]) * B[5]
-    return B1, B2, B3, B4, B5, B6
-end
-
-@inline function _bsplinebasis(P::FastBSplineSpace{0}, t::Real, i::Integer)
-    return _bsb0(P.knots,t,i)
-end
-
-@inline function _bsplinebasis(P::FastBSplineSpace{1}, t::Real, i::Integer)
-    return _bsb1(P.knots,t,i)
-end
-
-@inline function _bsplinebasis(P::FastBSplineSpace{2}, t::Real, i::Integer)
-    return _bsb2(P.knots,t,i)
-end
-
-@inline function _bsplinebasis(P::FastBSplineSpace{3}, t::Real, i::Integer)
-    return _bsb3(P.knots,t,i)
-end
-
-@inline function _bsplinebasis(P::FastBSplineSpace{4}, t::Real, i::Integer)
-    return _bsb4(P.knots,t,i)
-end
-
-@inline function _bsplinebasis(P::FastBSplineSpace{5}, t::Real, i::Integer)
-    return _bsb5(P.knots,t,i)
-end
diff --git a/src/_FastBSplineManifold.jl b/src/_FastBSplineManifold.jl
deleted file mode 100644
index 5ee2c99c..00000000
--- a/src/_FastBSplineManifold.jl
+++ /dev/null
@@ -1,256 +0,0 @@
-# Faster B-spline manifold
-
-"""
-B-spline manifold for lower polynomial degree
-TODO: make the field `bsplinespaces` to be conposite type, not abstract type, for performance
-"""
-struct FastBSplineManifold{Dim,T} <: AbstractBSplineManifold{Dim}
-    bsplinespaces::Vector{<:FastBSplineSpace}
-    controlpoints::Array{T,Dim} where T<:Point
-    function FastBSplineManifold(Ps::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{T}) where T
-        Ps = FastBSplineSpace.(Ps)
-        Dim = length(Ps)
-        if collect(size(a)) ≠ dim.(Ps)
-            throw(DimensionMismatch())
-        else
-            P = convert(Vector{FastBSplineSpace}, Ps)
-            new{Dim,T}(P, float(a))
-        end
-    end
-end
-
-struct BSplineCurve{p1,T<:Point} <: AbstractBSplineManifold{1}
-    bsplinespace1::FastBSplineSpace{p1}
-    controlpoints::Array{T,1}
-    function BSplineCurve(P1::AbstractBSplineSpace, a::AbstractArray{T,1}) where T
-        p1 = degree(P1)
-        if size(a) ≠ (dim(P1),)
-            throw(DimensionMismatch())
-        else
-            new{p1,T}(P1, float(a))
-        end
-    end
-end
-function BSplineCurve(P::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{<:Point,1})
-    if length(P) ≠ 1
-        throw(DimensionMismatch())
-    else
-        BSplineCurve(P[1], a)
-    end
-end
-
-struct BSplineSurface{p1,p2,T<:Point} <: AbstractBSplineManifold{2}
-    bsplinespace1::FastBSplineSpace{p1}
-    bsplinespace2::FastBSplineSpace{p2}
-    controlpoints::Array{T,2}
-    function BSplineSurface(P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, a::AbstractArray{T,2}) where T
-        p1, p2 = degree(P1), degree(P2)
-        if size(a) ≠ (dim(P1), dim(P2))
-            throw(DimensionMismatch())
-        else
-            new{p1,p2,T}(P1, P2, float(a))
-        end
-    end
-end
-function BSplineSurface(P::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{<:Point,2})
-    if length(P) ≠ 2
-        throw(DimensionMismatch())
-    else
-        BSplineSurface(P[1], P[2], a)
-    end
-end
-
-struct BSplineSolid{p1,p2,p3,T<:Point} <: AbstractBSplineManifold{3}
-    bsplinespace1::FastBSplineSpace{p1}
-    bsplinespace2::FastBSplineSpace{p2}
-    bsplinespace3::FastBSplineSpace{p3}
-    controlpoints::Array{T,3}
-    function BSplineSolid(P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, P3::AbstractBSplineSpace, a::AbstractArray{T,3}) where T
-        p1, p2, p3 = degree(P1), degree(P2), degree(P3)
-        if size(a) ≠ (dim(P1), dim(P2), dim(P3))
-            throw(DimensionMismatch())
-        else
-            new{p1,p2,p3,T}(P1, P2, P3, float(a))
-        end
-    end
-end
-function BSplineSolid(P::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{<:Point,3})
-    if length(P) ≠ 3
-        throw(DimensionMismatch())
-    else
-        BSplineSolid(P[1], P[2], P[3], a)
-    end
-end
-
-"""
-convert AbstractBSplineManifold to FastBSplineManifold
-"""
-function FastBSplineManifold(M::AbstractBSplineManifold)
-    return FastBSplineManifold(BSplineSpace.(M.bsplinespaces), M.controlpoints)
-end
-
-function (M::FastBSplineManifold)(t1::Real)
-    P = bsplinespaces(M)
-    a = controlpoints(M)
-    P1 = P[1]
-    n1 = dim(P1)
-    sum([bsplinebasis(P1,i1,t1)*a[i1] for i1 in 1:n1])
-end
-
-function (M::FastBSplineManifold)(t1::Real, t2::Real)
-    P = bsplinespaces(M)
-    a = controlpoints(M)
-    P1, P2 = P[1], P[2]
-    n1, n2 = dim.(P)
-    sum([bsplinebasis(P1,i1,t1)*bsplinebasis(P2,i2,t2)*a[i1] for i1 in 1:n1, i2 in 1:n2])
-end
-
-function (M::FastBSplineManifold)(t1::Real, t2::Real, t3::Real)
-    P = bsplinespaces(M)
-    a = controlpoints(M)
-    P1, P2, P3 = P
-    n1, n2, n2 = dim.(P)
-    sum([bsplinebasis(P1,i1,t1)*bsplinebasis(P2,i2,t2)*bsplinebasis(P3,i3,t3)*a[i1] for i1 in 1:n1, i2 in 1:n2, i3 in 1:n3])
-end
-
-function (M::BSplineCurve{p1})(t1::Real) where {p1}
-    P1, = bsplinespaces(M)
-    a = controlpoints(M)
-    n1 = dim(P1)
-    k1 = P1.knots
-    j1 = _knotindex(P1, t1)
-    b1 = _bsplinebasis(P1,t1,j1)
-
-    f1 = j1-p1
-    l1 = j1
-    S1 = b1[1]*a[f1]
-    for i1 in f1+1:l1
-        S1 += b1[i1-f1+1]*a[i1]
-    end
-
-    return S1
-end
-
-function (M::BSplineSurface{p1,p2})(t1::Real,t2::Real) where {p1} where {p2}
-    P1, P2 = bsplinespaces(M)
-    a = controlpoints(M)
-    n1, n2 = dim(P1), dim(P2)
-    k1, k2 = P1.knots, P2.knots
-    j1, j2 = _knotindex(P1, t1), _knotindex(P2, t2)
-    b1 = _bsplinebasis(P1,t1,j1)
-    b2 = _bsplinebasis(P2,t2,j2)
-
-    f1, f2 = j1-p1, j2-p2
-    l1, l2 = j1, j2
-    S2 = b2[1]*a[f1,f2]
-    for i2 in f2+1:l2
-        S2 += b2[i2-f2+1]*a[f1,i2]
-    end
-    S1 = b1[1]*S2
-    for i1 in f1+1:l1
-        S2 = b2[1]*a[i1,f2]
-        for i2 in f2+1:l2
-            S2 += b2[i2-f2+1]*a[i1,i2]
-        end
-        S1 += b1[i1-f1+1]*S2
-    end
-
-    return S1
-end
-
-function (M::BSplineSolid{p1,p2,p3})(t1::Real,t2::Real,t3::Real) where {p1} where {p2} where {p3}
-    P1, P2, P3 = M.bsplinespace1, M.bsplinespace2, M.bsplinespace3
-    a = controlpoints(M)
-    n1, n2, n3 = dim(P1), dim(P2), dim(P3)
-    k1, k2, k3 = P1.knots, P2.knots, P3.knots
-    j1, j2, j3 = _knotindex(P1, t1), _knotindex(P2, t2), _knotindex(P3, t3)
-    b1 = _bsplinebasis(P1,t1,j1)
-    b2 = _bsplinebasis(P2,t2,j2)
-    b3 = _bsplinebasis(P3,t3,j3)
-
-    f1, f2, f3 = j1-p1, j2-p2, j3-p3
-    l1, l2, l3 = j1, j2, j3
-    S3 = b3[1]*a[f1,f2,f3]
-    for i3 in f3+1:l3
-        S3 += b3[i3-f3+1]*a[f1,f2,i3]
-    end
-    S2 = b1[1]*b2[1]*S3
-    for i2 in f2+1:l2
-        S3 = b3[1]*a[f1,i2,f3]
-        for i3 in f3+1:l3
-            S3 += b3[i3-f3+1]*a[f1,i2,i3]
-        end
-        S2 += b1[1]*b2[i2-f2+1]*S3
-    end
-    S1 = S2
-    for i1 in f1+1:l1
-        S3 = b3[1]*a[i1,f2,f3]
-        for i3 in f3+1:l3
-            S3 += b3[i3-f3+1]*a[i1,f2,i3]
-        end
-        S2 = b1[i1-f1+1]*b2[1]*S3
-        for i2 in f2+1:l2
-            S3 = b3[1]*a[i1,i2,f3]
-            for i3 in f3+1:l3
-                S3 += b3[i3-f3+1]*a[i1,i2,i3]
-            end
-            S2 += b1[i1-f1+1]*b2[i2-f2+1]*S3
-        end
-        S1 += S2
-    end
-
-    return S1
-end
-
-# function (M::Union{BSplineCurve,BSplineSurface,BSplineSolid})(t::AbstractVector{<:Real})
-#     return M(t...)
-# end
-
-function (M::BSplineCurve)(t::AbstractVector{<:Real})
-    return M(t...)
-end
-
-function (M::BSplineSurface)(t::AbstractVector{<:Real})
-    return M(t...)
-end
-
-function (M::BSplineSolid)(t::AbstractVector{<:Real})
-    return M(t...)
-end
-
-dim(M::FastBSplineManifold) = length(M.bsplinespaces)
-dim(M::BSplineCurve) = 1
-dim(M::BSplineSurface) = 2
-dim(M::BSplineSolid) = 3
-
-function bsplinespaces(M::FastBSplineManifold)
-    return M.bsplinespaces
-end
-
-function bsplinespaces(M::BSplineCurve{p1}) where {p1}
-    return (M.bsplinespace1, )
-end
-
-function bsplinespaces(M::BSplineSurface{p1,p2}) where {p1} where {p2}
-    return (M.bsplinespace1, M.bsplinespace2)
-end
-
-function bsplinespaces(M::BSplineSolid{p1,p2,p3}) where {p1} where {p2} where {p3}
-    return (M.bsplinespace1, M.bsplinespace2, M.bsplinespace3)
-end
-
-function controlpoints(M::FastBSplineManifold)
-    return M.controlpoints
-end
-
-function controlpoints(M::BSplineCurve{p1}) where {p1}
-    return M.controlpoints
-end
-
-function controlpoints(M::BSplineSurface{p1,p2}) where {p1} where {p2}
-    return M.controlpoints
-end
-
-function controlpoints(M::BSplineSolid{p1,p2,p3}) where {p1} where {p2} where {p3}
-    return M.controlpoints
-end
diff --git a/src/_FastBSplineSpace.jl b/src/_FastBSplineSpace.jl
deleted file mode 100644
index c25e0088..00000000
--- a/src/_FastBSplineSpace.jl
+++ /dev/null
@@ -1,51 +0,0 @@
-# Faster B-spline space
-
-@doc raw"""
-B-spline space for lower polynomial degree
-```math
-\mathcal{P}[p,k]
-```
-This type `FastBSplineSpace` is faster than `BSplineSpace`, but the degree must be equal or less than `MAX_DEGREE`.
-"""
-struct FastBSplineSpace{p} <: AbstractBSplineSpace{p,Float64}
-    knots::Knots
-    function FastBSplineSpace(p::Integer, knots::Knots)
-        if p < 0
-            throw(DomainError(p, "degree of polynominal must be non-negative"))
-        elseif p > MAX_DEGREE
-            throw(DomainError(p, "FastBSpline supports only degree 0,...,$(MAX_DEGREE)"))
-        end
-        new{p}(knots)
-    end
-end
-function FastBSplineSpace{q}(p::Integer, knots::Knots) where {q}
-    FastBSplineSpace(p,knots)
-end
-function FastBSplineSpace{p}(knots::Knots) where {p}
-    FastBSplineSpace(p,knots)
-end
-
-"""
-convert AbstractBSplineSpace to FastBSplineSpace
-"""
-function FastBSplineSpace(P::AbstractBSplineSpace)
-    return FastBSplineSpace(degree(P), knots(P))
-end
-
-# """
-# convert AbstractBSplineSpace to FastBSplineSpace
-# """
-# function Base.convert(::Type{<:FastBSplineSpace}, P::AbstractBSplineSpace)
-#     return FastBSplineSpace(degree(P), knots(P))
-# end
-
-function degree(::FastBSplineSpace{p}) where {p}
-    return p
-end
-
-function knots(P::FastBSplineSpace)
-    return P.knots
-end
-
-lower(::Type{FastBSplineSpace{p}}) where p = FastBSplineSpace{p-1}
-lower(P::FastBSplineSpace{p}) where {p} = FastBSplineSpace{p-1}(knots(P))
diff --git a/src/_Fitting.jl b/src/_Fitting.jl
index e0216826..fc683506 100644
--- a/src/_Fitting.jl
+++ b/src/_Fitting.jl
@@ -9,8 +9,7 @@ Assumption:
 * ``i ≤ j``
 * ``1 ≤ n ≤ p-j+i+1``
 """
-function _b_b_int(P::AbstractBSplineSpace, i, j, n, gl)
-    p = degree(P)
+function _b_b_int(P::BSplineSpace{p}, i, j, n, gl) where p
     I = knots(P)[j+n-1]..knots(P)[j+n]
 
     f(t) = bsplinebasis(P, i, t) * bsplinebasis(P, j, t)
@@ -30,9 +29,7 @@ Calculate
 \end{aligned}
 ```
 """
-function _b_b_int_R(P::AbstractBSplineSpace, i, j, gl)
-    p = degree(P)
-    k = knots(P)
+function _b_b_int_R(P::BSplineSpace{p}, i, j, gl) where p
     Δ = j - i
     if Δ < -p
         return 0.0
@@ -61,9 +58,7 @@ Calculate
 \end{aligned}
 ```
 """
-function _b_b_int_I(P::AbstractBSplineSpace, i, j, gl)
-    p = degree(P)
-    k = knots(P)
+function _b_b_int_I(P::BSplineSpace{p}, i, j, gl) where p
     n = dim(P)
     Δ = j - i
     if Δ < -p
@@ -85,25 +80,22 @@ function _b_b_int_I(P::AbstractBSplineSpace, i, j, gl)
     end
 end
 
-function innerproduct_R(P::AbstractBSplineSpace)
-    p = degree(P)
+function innerproduct_R(P::BSplineSpace{p}) where p
     n = dim(P)
     nip = p + 1
     gl = GaussLegendre(nip)
     return [_b_b_int_R(P, i, j, gl) for i in 1:n, j in 1:n]
 end
 
-function innerproduct_I(P::AbstractBSplineSpace)
-    p = degree(P)
+function innerproduct_I(P::BSplineSpace{p}) where p
     n = dim(P)
     nip = p + 1
     gl = GaussLegendre(nip)
     return [_b_b_int_I(P, i, j, gl) for i in 1:n, j in 1:n]
 end
 
-function _f_b_int_R(func, i1, P1::AbstractBSplineSpace, gl1)
+function _f_b_int_R(func, i1, P1::BSplineSpace{p1}, gl1) where {p1}
     k1 = knots(P1)
-    p1 = degree(P1)
     F(t1) = bsplinebasis(P1, i1, t1) * func(t1)
 
     f1,l1 = i1, i1+p1
@@ -115,10 +107,9 @@ function _f_b_int_R(func, i1, P1::AbstractBSplineSpace, gl1)
     return S1
 end
 
-function _f_b_int_I(func, i1, P1::AbstractBSplineSpace, gl1)
+function _f_b_int_I(func, i1, P1::BSplineSpace{p1}, gl1) where {p1}
     k1 = knots(P1)
     m1 = length(k1)
-    p1 = degree(P1)
     F(t1) = bsplinebasis(P1, i1, t1) * func(t1)
 
     f1,l1 = max(i1, p1+1), min(i1+p1, m1-p1-1)
@@ -130,9 +121,9 @@ function _f_b_int_I(func, i1, P1::AbstractBSplineSpace, gl1)
     return S1
 end
 
-function _f_b_int_R(func, i1, i2, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, gl1, gl2)
+function _f_b_int_R(func, i1, i2, P::Tuple{<:AbstractBSplineSpace{p1}, <:AbstractBSplineSpace{p2}}, gl1, gl2) where {p1,p2}
+    P1, P2 = P
     k1, k2 = knots(P1), knots(P2)
-    p1, p2 = degree(P1), degree(P2)
     F(t1, t2) = bsplinebasis(P1, i1, t1) * bsplinebasis(P2, i2, t2) * func(t1, t2)
 
     f1, l1 = i1, i1+p1
@@ -153,10 +144,10 @@ function _f_b_int_R(func, i1, i2, P1::AbstractBSplineSpace, P2::AbstractBSplineS
     return S1
 end
 
-function _f_b_int_I(func, i1, i2, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, gl1, gl2)
+function _f_b_int_I(func, i1, i2, P::Tuple{<:AbstractBSplineSpace{p1}, <:AbstractBSplineSpace{p2}}, gl1, gl2) where {p1,p2}
+    P1, P2 = P
     k1, k2 = knots(P1), knots(P2)
     m1, m2 = length(k1), length(k2)
-    p1, p2 = degree(P1), degree(P2)
     F(t1, t2) = bsplinebasis(P1, i1, t1) * bsplinebasis(P2, i2, t2) * func(t1, t2)
 
     f1, l1 = max(i1, p1+1), min(i1+p1, m1-p1-1)
@@ -177,9 +168,9 @@ function _f_b_int_I(func, i1, i2, P1::AbstractBSplineSpace, P2::AbstractBSplineS
     return S1
 end
 
-function _f_b_int_R(func, i1, i2, i3, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, P3::AbstractBSplineSpace, gl1, gl2, gl3)
+function _f_b_int_R(func, i1, i2, i3, P::Tuple{<:AbstractBSplineSpace{p1}, <:AbstractBSplineSpace{p2}, <:AbstractBSplineSpace{p3}}, gl1, gl2, gl3) where {p1,p2,p3}
+    P1, P2, P3 = P
     k1, k2, k3 = knots(P1), knots(P2), knots(P3)
-    p1, p2, p3 = degree(P1), degree(P2), degree(P3)
     F(t1, t2, t3) = bsplinebasis(P1, i1, t1) * bsplinebasis(P2, i2, t2) * bsplinebasis(P3, i3, t3) * func(t1, t2, t3)
 
     f1, l1 = i1, i1+p1
@@ -217,10 +208,10 @@ function _f_b_int_R(func, i1, i2, i3, P1::AbstractBSplineSpace, P2::AbstractBSpl
     return S1
 end
 
-function _f_b_int_I(func, i1, i2, i3, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, P3::AbstractBSplineSpace, gl1, gl2, gl3)
+function _f_b_int_I(func, i1, i2, i3, P::Tuple{<:AbstractBSplineSpace{p1}, <:AbstractBSplineSpace{p2}, <:AbstractBSplineSpace{p3}}, gl1, gl2, gl3) where {p1,p2,p3}
+    P1, P2, P3 = P
     k1, k2, k3 = knots(P1), knots(P2), knots(P3)
     m1, m2, m3 = length(k1), length(k2), length(k3)
-    p1, p2, p3 = degree(P1), degree(P2), degree(P3)
     F(t1, t2, t3) = bsplinebasis(P1, i1, t1) * bsplinebasis(P2, i2, t2) * bsplinebasis(P3, i3, t3) * func(t1, t2, t3)
 
     f1, l1 = max(i1, p1+1), min(i1+p1, m1-p1-1)
@@ -258,8 +249,7 @@ function _f_b_int_I(func, i1, i2, i3, P1::AbstractBSplineSpace, P2::AbstractBSpl
     return S1
 end
 
-function innerproduct_R(func, P1::AbstractBSplineSpace)
-    p1 = degree(P1)
+function innerproduct_R(func, P1::BSplineSpace{p1}) where {p1}
     n1 = dim(P1)
     nip1 = p1 + 1
     gl1 = GaussLegendre(nip1)
@@ -267,8 +257,7 @@ function innerproduct_R(func, P1::AbstractBSplineSpace)
     return b
 end
 
-function innerproduct_I(func, P1::AbstractBSplineSpace)
-    p1 = degree(P1)
+function innerproduct_I(func, P1::BSplineSpace{p1}) where {p1}
     n1 = dim(P1)
     nip1 = p1 + 1
     gl1 = GaussLegendre(nip1)
@@ -276,43 +265,40 @@ function innerproduct_I(func, P1::AbstractBSplineSpace)
     return b
 end
 
-function innerproduct_R(func, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace)
-    p1, p2 = degree(P1), degree(P2)
-    n1, n2 = dim(P1), dim(P2)
+function innerproduct_R(func, P::Tuple{<:AbstractBSplineSpace{p1}, <:AbstractBSplineSpace{p2}}) where {p1,p2}
+    n1, n2 = dim.(P)
     nip1 = p1 + 1
     nip2 = p2 + 1
     gl1 = GaussLegendre(nip1)
     gl2 = GaussLegendre(nip2)
-    b = [_f_b_int_R(func, i1, i2, P1, P2, gl1, gl2) for i1 in 1:n1, i2 in 1:n2]
+    b = [_f_b_int_R(func, i1, i2, P, gl1, gl2) for i1 in 1:n1, i2 in 1:n2]
     return b
 end
 
-function innerproduct_I(func, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace)
-    p1, p2 = degree(P1), degree(P2)
-    n1, n2 = dim(P1), dim(P2)
+function innerproduct_I(func, P::Tuple{<:AbstractBSplineSpace{p1}, <:AbstractBSplineSpace{p2}}) where {p1,p2}
+    n1, n2 = dim.(P)
     nip1 = p1 + 1
     nip2 = p2 + 1
     gl1 = GaussLegendre(nip1)
     gl2 = GaussLegendre(nip2)
-    b = [_f_b_int_I(func, i1, i2, P1, P2, gl1, gl2) for i1 in 1:n1, i2 in 1:n2]
+    b = [_f_b_int_I(func, i1, i2, P, gl1, gl2) for i1 in 1:n1, i2 in 1:n2]
     return b
 end
 
-function innerproduct_R(func, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, P3::AbstractBSplineSpace)
-    p1, p2, p3 = degree(P1), degree(P2), degree(P3)
-    n1, n2, n3 = dim(P1), dim(P2), dim(P3)
+function innerproduct_R(func, P::Tuple{<:AbstractBSplineSpace{p1}, <:AbstractBSplineSpace{p2}, <:AbstractBSplineSpace{p3}}) where {p1,p2,p3}
+    n1, n2, n3 = dim.(P)
     nip1 = p1 + 1
     nip2 = p2 + 1
     nip3 = p3 + 1
     gl1 = GaussLegendre(nip1)
     gl2 = GaussLegendre(nip2)
     gl3 = GaussLegendre(nip3)
-    b = [_f_b_int_R(func, i1, i2, i3, P1, P2, P3, gl1, gl2, gl3) for i1 in 1:n1, i2 in 1:n2, i3 in 1:n3]
+    b = [_f_b_int_R(func, i1, i2, i3, P, gl1, gl2, gl3) for i1 in 1:n1, i2 in 1:n2, i3 in 1:n3]
     return b
 end
 
-function innerproduct_I(func, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, P3::AbstractBSplineSpace)
-    p1, p2, p3 = degree(P1), degree(P2), degree(P3)
+function innerproduct_I(func, P::Tuple{<:AbstractBSplineSpace{p1}, <:AbstractBSplineSpace{p2}, <:AbstractBSplineSpace{p3}}) where {p1,p2,p3}
+    P1, P2, P3 = P
     n1, n2, n3 = dim(P1), dim(P2), dim(P3)
     nip1 = p1 + 1
     nip2 = p2 + 1
@@ -320,14 +306,15 @@ function innerproduct_I(func, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace
     gl1 = GaussLegendre(nip1)
     gl2 = GaussLegendre(nip2)
     gl3 = GaussLegendre(nip3)
-    b = [_f_b_int_I(func, i1, i2, i3, P1, P2, P3, gl1, gl2, gl3) for i1 in 1:n1, i2 in 1:n2, i3 in 1:n3]
+    b = [_f_b_int_I(func, i1, i2, i3, P, gl1, gl2, gl3) for i1 in 1:n1, i2 in 1:n2, i3 in 1:n3]
     return b
 end
 
 """
 * func: Real -> ℝ-vector space
 """
-function fittingcontrolpoints(func, P1::AbstractBSplineSpace; domain=:I)
+function fittingcontrolpoints(func, P::Tuple{<:AbstractBSplineSpace{p1}}; domain=:I) where {p1}
+    P1, = P
     if domain == :I
         b = innerproduct_I(func, P1)
         A = innerproduct_I(P1)
@@ -341,14 +328,15 @@ end
 """
 * func: (Real,Real) -> ℝ-vector space
 """
-function fittingcontrolpoints(func, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace; domain=:I)
-    n1, n2 = dim(P1), dim(P2)
+function fittingcontrolpoints(func, P::Tuple{<:AbstractBSplineSpace{p1}, <:AbstractBSplineSpace{p2}}; domain=:I) where {p1,p2}
+    P1, P2 = P
+    n1, n2 = dim.(P)
     if domain == :I
         A1, A2 = innerproduct_I(P1), innerproduct_I(P2)
-        b = innerproduct_I(func, P1, P2)
+        b = innerproduct_I(func, P)
     elseif domain == :R
         A1, A2 = innerproduct_R(P1), innerproduct_R(P2)
-        b = innerproduct_R(func, P1, P2)
+        b = innerproduct_R(func, P)
     end
     A = [A1[i1, j1] * A2[i2, j2] for i1 in 1:n1, i2 in 1:n2, j1 in 1:n1, j2 in 1:n2]
     _A = reshape(A, n1 * n2, n1 * n2)
@@ -359,29 +347,18 @@ end
 """
 * func: (Real,Real,Real) -> ℝ-vector space
 """
-function fittingcontrolpoints(func, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, P3::AbstractBSplineSpace; domain=:I)
+function fittingcontrolpoints(func, P::Tuple{<:AbstractBSplineSpace{p1}, <:AbstractBSplineSpace{p2}, <:AbstractBSplineSpace{p3}}; domain=:I) where {p1,p2,p3}
+    P1, P2, P3 = P
     n1, n2, n3 = dim(P1), dim(P2), dim(P3)
     if domain == :I
         A1, A2, A3 = innerproduct_I(P1), innerproduct_I(P2), innerproduct_I(P3)
-        b = innerproduct_I(func, P1, P2, P3)
+        b = innerproduct_I(func, P)
     elseif domain == :R
         A1, A2, A3 = innerproduct_R(P1), innerproduct_R(P2), innerproduct_R(P3)
-        b = innerproduct_R(func, P1, P2, P3)
+        b = innerproduct_R(func, P)
     end
     A = [A1[i1, j1] * A2[i2, j2] * A3[i3, j3] for i1 in 1:n1, i2 in 1:n2, i3 in 1:n3, j1 in 1:n1, j2 in 1:n2, j3 in 1:n3]
     _A = reshape(A, n1 * n2 * n3, n1 * n2 * n3)
     _b = reshape(b, n1 * n2 * n3)
     return reshape(inv(_A) * _b, n1, n2, n3)
 end
-
-
-"""
-Approximate given function by linear combination of B-spline functions.
-This function returns its control points.
-* func: Vector{<:Real} -> ℝ-vector space
-"""
-function fittingcontrolpoints(func, Ps::AbstractVector{<:AbstractBSplineSpace}; domain=:I)
-    # TODO: currently, this function only supports for 1-dim, 2-dim and 3-dim B-spline manifold.
-    _func(t...) = func([t...])
-    return fittingcontrolpoints(_func, Ps..., domain=domain)
-end
diff --git a/src/_Refinement.jl b/src/_Refinement.jl
index 1391977d..ee1cc865 100644
--- a/src/_Refinement.jl
+++ b/src/_Refinement.jl
@@ -1,64 +1,76 @@
 # Refinement
 
 # TODO: general dimension
-
-
-refinement(M::BSplineManifold, Ps′::Tuple) = refinement(M, Ps′...)
+# TODO: Update docstrings
 
 @doc raw"""
 Refinement of B-spline manifold with given B-spline spaces.
 """
-function refinement(M::BSplineManifold{1}, Ps′::AbstractBSplineSpace...)
+refinement
+
+function refinement(M::BSplineManifold{1}, Ps′::Tuple{<:AbstractBSplineSpace})
     Ps = bsplinespaces(M)
     a = controlpoints(M)
-    n = dim.(Ps)
-    n′ = dim.(Ps′)
     d = size(a,2)  # 2 == Dim+1
+    (n1,) = dim.(Ps)
+    (n1′,) = dim.(Ps′)
+    (A1,) = changebasis.(Ps, Ps′)
 
-    A = changebasis.(Ps, Ps′)
-    a′ = [sum(A[1][I₁, J₁] * a[I₁,i] for I₁ in 1:n[1]) for J₁ in 1:n′[1], i in 1:d]
+    a′ = [sum(A1[I₁,J₁] * a[I₁,j] for I₁ in 1:n1) for J₁ in 1:n1′, j in 1:d]
     return BSplineManifold(Ps′, a′)
 end
-
-@doc raw"""
-Refinement of B-spline manifold with given B-spline spaces.
-"""
-function refinement(M::BSplineManifold{2}, Ps′::AbstractBSplineSpace...)
+function refinement(M::BSplineManifold{2}, Ps′::Tuple{<:AbstractBSplineSpace,<:AbstractBSplineSpace})
     Ps = bsplinespaces(M)
     a = controlpoints(M)
-    n = dim.(Ps)
-    n′ = dim.(Ps′)
     d = size(a,3)  # 3 == Dim+1
+    (n1,n2) = dim.(Ps)
+    (n1′,n2′) = dim.(Ps′)
+    (A1,A2) = changebasis.(Ps, Ps′)
 
-    A = changebasis.(Ps, Ps′)
-    a′ = [sum(A[1][I₁, J₁] * A[2][I₂, J₂] * a[I₁, I₂,i] for I₁ in 1:n[1], I₂ in 1:n[2]) for J₁ in 1:n′[1], J₂ in 1:n′[2], i in 1:d]
+    a′ = [sum(A1[I₁,J₁] * A2[I₂,J₂] * a[I₁,I₂,j] for I₁ in 1:n1, I₂ in 1:n2) for J₁ in 1:n1′, J₂ in 1:n2′, j in 1:d]
     return BSplineManifold(Ps′, a′)
 end
+function refinement(M::BSplineManifold{3}, Ps′::Tuple{<:AbstractBSplineSpace,<:AbstractBSplineSpace,<:AbstractBSplineSpace})
+    Ps = bsplinespaces(M)
+    a = controlpoints(M)
+    d = size(a,4)  # 4 == Dim+1
+    (n1,n2,n3) = dim.(Ps)
+    (n1′,n2′,n3′) = dim.(Ps′)
+    (A1,A2,A3) = changebasis.(Ps, Ps′)
 
+    a′ = [sum(A1[I₁,J₁] * A2[I₂,J₂]* A3[I₃,J₃] * a[I₁,I₂,I₃,j] for I₁ in 1:n1, I₂ in 1:n2, I₃ in 1:n3) for J₁ in 1:n1′, J₂ in 1:n2′, J₃ in 1:n3′, j in 1:d]
+    return BSplineManifold(Ps′, a′)
+end
 
-@doc raw"""
-Refinement of B-spline manifold with given B-spline spaces.
-"""
-function refinement(M::AbstractBSplineManifold, Ps′::Vector{<:AbstractBSplineSpace})
-    Ps = collect(bsplinespaces(M))
+function refinement(M::CustomBSplineManifold{1}, Ps′::Tuple{<:AbstractBSplineSpace})
+    Ps = bsplinespaces(M)
     a = controlpoints(M)
-    d = length(Ps)
-    n = dim.(Ps)
-    n′ = dim.(Ps′)
-
-    A = changebasis.(Ps, Ps′)
-    # TODO: general dimension
-    if d == 1
-        a′ = [sum(A[1][I₁, J₁] * a[I₁] for I₁ in 1:n[1]) for J₁ in 1:n′[1]]
-    elseif d == 2
-        a′ = [sum(A[1][I₁, J₁] * A[2][I₂, J₂] * a[I₁, I₂] for I₁ in 1:n[1], I₂ in 1:n[2]) for J₁ in 1:n′[1], J₂ in 1:n′[2]]
-    elseif d == 3
-        a′ = [
-            sum(A[1][I₁, J₁] * A[2][I₂, J₂] * A[3][I₃, J₃] * a[I₁, I₂, I₃] for I₁ in 1:n[1], I₂ in 1:n[2], I₃ in 1:n[3])
-            for J₁ in 1:n′[1], J₂ in 1:n′[2], J₃ in 1:n′[3]
-        ]
-    end
-    return (typeof(M).name.wrapper)(Ps′, a′)
+    (n1,) = dim.(Ps)
+    (n1′,) = dim.(Ps′)
+    (A1,) = changebasis.(Ps, Ps′)
+
+    a′ = [sum(A1[I₁,J₁] * a[I₁] for I₁ in 1:n1) for J₁ in 1:n1′]
+    return CustomBSplineManifold(Ps′, a′)
+end
+function refinement(M::CustomBSplineManifold{2}, Ps′::Tuple{<:AbstractBSplineSpace,<:AbstractBSplineSpace})
+    Ps = bsplinespaces(M)
+    a = controlpoints(M)
+    (n1,n2) = dim.(Ps)
+    (n1′,n2′) = dim.(Ps′)
+    (A1,A2) = changebasis.(Ps, Ps′)
+
+    a′ = [sum(A1[I₁,J₁] * A2[I₂,J₂] * a[I₁,I₂] for I₁ in 1:n1, I₂ in 1:n2) for J₁ in 1:n1′, J₂ in 1:n2′]
+    return CustomBSplineManifold(Ps′, a′)
+end
+function refinement(M::CustomBSplineManifold{3}, Ps′::Tuple{<:AbstractBSplineSpace,<:AbstractBSplineSpace,<:AbstractBSplineSpace})
+    Ps = bsplinespaces(M)
+    a = controlpoints(M)
+    (n1,n2,n3) = dim.(Ps)
+    (n1′,n2′,n3′) = dim.(Ps′)
+    (A1,A2,A3) = changebasis.(Ps, Ps′)
+
+    a′ = [sum(A1[I₁,J₁] * A2[I₂,J₂]* A3[I₃,J₃] * a[I₁,I₂,I₃] for I₁ in 1:n1, I₂ in 1:n2, I₃ in 1:n3) for J₁ in 1:n1′, J₂ in 1:n2′, J₃ in 1:n3′]
+    return CustomBSplineManifold(Ps′, a′)
 end
 
 
@@ -82,7 +94,7 @@ function refinement(M::BSplineManifold{Dim}; p₊::Union{Nothing,NTuple{Dim,Int}
     p = degree(P);
     k = knots(P);
     k_unique = unique(k[1+p:end-p]);
-    BSplineSpace{p + p₊[i]}(k + p₊[i] * k_unique + k₊[i])) for i in eachindex(Ps)]
+    BSplineSpace{p + p₊[i]}(k + p₊[i] * k_unique + k₊[i])) for i in 1:Dim]
 
     return refinement(M, tuple(Ps′...))
 end
@@ -90,7 +102,7 @@ end
 @doc raw"""
 Refinement of B-spline manifold with additional degree and knots.
 """
-function refinement(M::AbstractBSplineManifold; p₊::Union{Nothing,AbstractVector{<:Integer}} = nothing, k₊::Union{Nothing,Vector{<:Knots}} = nothing)
+function refinement(M::CustomBSplineManifold; p₊::Union{Nothing,NTuple{Dim,Int}}=nothing, k₊::Union{Nothing,NTuple{Dim,Knots{T}}}=nothing) where {Dim, T}
     Ps = collect(bsplinespaces(M))
     𝒂 = controlpoints(M)
     d = length(Ps)
@@ -110,7 +122,7 @@ function refinement(M::AbstractBSplineManifold; p₊::Union{Nothing,AbstractVect
     p = degree(P);
     k = knots(P);
     k_unique = unique(k[1+p:end-p]);
-    typeof(P).name.wrapper{p + p₊[i]}(k + p₊[i] * k_unique + k₊[i])) for i in eachindex(Ps)]
+    BSplineSpace{p + p₊[i]}(k + p₊[i] * k_unique + k₊[i])) for i in 1:Dim]
 
-    return refinement(M, Ps′)
+    return refinement(M, tuple(Ps′...))
 end


Test patch:
diff --git a/test/runtests.jl b/test/runtests.jl
index 41441c17..2cd37da5 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,7 +9,6 @@ using Random
     include("test_BSplineSpace.jl")
     include("test_BSplineBasis.jl")
     include("test_ChangeBasis.jl")
-    include("test_FastBSplineSpace.jl")
     include("test_BSplineManifold.jl")
     include("test_Integral.jl")
     include("test_Fitting.jl")
diff --git a/test/test_BSplineBasis.jl b/test/test_BSplineBasis.jl
index ffe5a8bc..3217cda5 100644
--- a/test/test_BSplineBasis.jl
+++ b/test/test_BSplineBasis.jl
@@ -21,27 +21,22 @@ end
         v = rand(10)
         k = Knots(v) + (p+1)*Knots(0, 1) + Knots(v[2], v[3])
         P = BSplineSpace{p}(k)
-        fP = FastBSplineSpace(p, k)
         n = dim(P)
 
-        @test degree(P) == degree(fP) == 0
+        @test degree(P) == 0
 
         # Check the values of B-spline basis funcitons
         for i in 1:n, t in rand(5)
             for t in rand(5)
                 @test (bsplinebasis₊₀(P, i, t)
-                     ≈ bsplinebasis₊₀(fP, i, t)
                      ≈ bsplinebasis₋₀(P, i, t)
-                     ≈ bsplinebasis₋₀(fP, i, t)
-                     ≈ bsplinebasis(P, i, t)
-                     ≈ bsplinebasis(fP, i, t))
+                     ≈ bsplinebasis(P, i, t))
             end
             for t in k
                 t₊₀ = nextfloat(t)
                 t₋₀ = prevfloat(t)
-                @test bsplinebasis₊₀(P, i, t) ≈ bsplinebasis₊₀(fP, i, t) ≒ bsplinebasis₊₀(P, i, t₊₀) ≒ bsplinebasis₊₀(fP, i, t₊₀)
-                @test bsplinebasis₋₀(P, i, t) ≈ bsplinebasis₋₀(fP, i, t) ≒ bsplinebasis₋₀(P, i, t₋₀) ≒ bsplinebasis₋₀(fP, i, t₋₀)
-                @test bsplinebasis(P, i, t) ≈ bsplinebasis(fP, i, t)
+                @test bsplinebasis₊₀(P, i, t) ≒ bsplinebasis₊₀(P, i, t₊₀)
+                @test bsplinebasis₋₀(P, i, t) ≒ bsplinebasis₋₀(P, i, t₋₀)
             end
         end
 
@@ -49,18 +44,14 @@ end
         for i in 1:n
             for t in rand(5)
                 @test (bsplinebasis′₊₀(P, i, t)
-                     ≈ bsplinebasis′₊₀(fP, i, t)
                      ≈ bsplinebasis′₋₀(P, i, t)
-                     ≈ bsplinebasis′₋₀(fP, i, t)
-                     ≈ bsplinebasis′(P, i, t)
-                     ≈ bsplinebasis′(fP, i, t))
+                     ≈ bsplinebasis′(P, i, t))
             end
             for t in k
                 t₊₀ = nextfloat(t)
                 t₋₀ = prevfloat(t)
-                @test bsplinebasis′₊₀(P, i, t) ≈ bsplinebasis′₊₀(fP, i, t) ≒ bsplinebasis′₊₀(P, i, t₊₀) ≒ bsplinebasis′₊₀(fP, i, t₊₀)
-                @test bsplinebasis′₋₀(P, i, t) ≈ bsplinebasis′₋₀(fP, i, t) ≒ bsplinebasis′₋₀(P, i, t₋₀) ≒ bsplinebasis′₋₀(fP, i, t₋₀)
-                @test bsplinebasis′(P, i, t) ≈ bsplinebasis′(fP, i, t)
+                @test bsplinebasis′₊₀(P, i, t) ≒ bsplinebasis′₊₀(P, i, t₊₀)
+                @test bsplinebasis′₋₀(P, i, t) ≒ bsplinebasis′₋₀(P, i, t₋₀)
             end
         end
 
@@ -68,12 +59,12 @@ end
         @testset "derivative" begin
             Random.seed!(57)
             for i in 1:n, t in rand(5)
-                b′₊₀ = bsplinebasis′₊₀(fP, i, t)
-                b′₋₀ = bsplinebasis′₋₀(fP, i, t)
-                b′   = bsplinebasis′(fP, i, t)
-                a′₊₀ = (bsplinebasis₊₀(fP, i, t+Δt) - bsplinebasis₊₀(fP, i, t-Δt)) / 2Δt
-                a′₋₀ = (bsplinebasis₋₀(fP, i, t+Δt) - bsplinebasis₋₀(fP, i, t-Δt)) / 2Δt
-                a′   = (bsplinebasis(fP, i, t+Δt) - bsplinebasis(fP, i, t-Δt)) / 2Δt
+                b′₊₀ = bsplinebasis′₊₀(P, i, t)
+                b′₋₀ = bsplinebasis′₋₀(P, i, t)
+                b′   = bsplinebasis′(P, i, t)
+                a′₊₀ = (bsplinebasis₊₀(P, i, t+Δt) - bsplinebasis₊₀(P, i, t-Δt)) / 2Δt
+                a′₋₀ = (bsplinebasis₋₀(P, i, t+Δt) - bsplinebasis₋₀(P, i, t-Δt)) / 2Δt
+                a′   = (bsplinebasis(P, i, t+Δt) - bsplinebasis(P, i, t-Δt)) / 2Δt
                 @test b′₊₀ ≈ b′₋₀ ≈ b′ ≈ a′₊₀ ≈ a′₋₀ ≈ a′
             end
         end
@@ -86,27 +77,22 @@ end
         v = rand(10)
         k = Knots(v) + (p+1)*Knots(0, 1) + Knots(v[2], v[3])
         P = BSplineSpace{p}(k)
-        fP = FastBSplineSpace(p, k)
         n = dim(P)
 
-        @test degree(P) == degree(fP) == 1
+        @test degree(P) == 1
 
         # Check the values of B-spline basis funcitons
         for i in 1:n, t in rand(5)
             for t in rand(5)
                 @test (bsplinebasis₊₀(P, i, t)
-                     ≈ bsplinebasis₊₀(fP, i, t)
                      ≈ bsplinebasis₋₀(P, i, t)
-                     ≈ bsplinebasis₋₀(fP, i, t)
-                     ≈ bsplinebasis(P, i, t)
-                     ≈ bsplinebasis(fP, i, t))
+                     ≈ bsplinebasis(P, i, t))
             end
             for t in k
                 t₊₀ = nextfloat(t)
                 t₋₀ = prevfloat(t)
-                @test bsplinebasis₊₀(P, i, t) ≈ bsplinebasis₊₀(fP, i, t) ≒ bsplinebasis₊₀(P, i, t₊₀) ≒ bsplinebasis₊₀(fP, i, t₊₀)
-                @test bsplinebasis₋₀(P, i, t) ≈ bsplinebasis₋₀(fP, i, t) ≒ bsplinebasis₋₀(P, i, t₋₀) ≒ bsplinebasis₋₀(fP, i, t₋₀)
-                @test bsplinebasis(P, i, t) ≈ bsplinebasis(fP, i, t)
+                @test bsplinebasis₊₀(P, i, t) ≒ bsplinebasis₊₀(P, i, t₊₀)
+                @test bsplinebasis₋₀(P, i, t)≒ bsplinebasis₋₀(P, i, t₋₀)
             end
         end
 
@@ -114,18 +100,14 @@ end
         for i in 1:n
             for t in rand(5)
                 @test (bsplinebasis′₊₀(P, i, t)
-                     ≈ bsplinebasis′₊₀(fP, i, t)
                      ≈ bsplinebasis′₋₀(P, i, t)
-                     ≈ bsplinebasis′₋₀(fP, i, t)
-                     ≈ bsplinebasis′(P, i, t)
-                     ≈ bsplinebasis′(fP, i, t))
+                     ≈ bsplinebasis′(P, i, t))
             end
             for t in k
                 t₊₀ = nextfloat(t)
                 t₋₀ = prevfloat(t)
-                @test bsplinebasis′₊₀(P, i, t) ≈ bsplinebasis′₊₀(fP, i, t) ≒ bsplinebasis′₊₀(P, i, t₊₀) ≒ bsplinebasis′₊₀(fP, i, t₊₀)
-                @test bsplinebasis′₋₀(P, i, t) ≈ bsplinebasis′₋₀(fP, i, t) ≒ bsplinebasis′₋₀(P, i, t₋₀) ≒ bsplinebasis′₋₀(fP, i, t₋₀)
-                @test bsplinebasis′(P, i, t) ≈ bsplinebasis′(fP, i, t)
+                @test bsplinebasis′₊₀(P, i, t) ≒ bsplinebasis′₊₀(P, i, t₊₀)
+                @test bsplinebasis′₋₀(P, i, t) ≒ bsplinebasis′₋₀(P, i, t₋₀)
             end
         end
 
@@ -133,12 +115,12 @@ end
         @testset "derivative" begin
             Random.seed!(57)
             for i in 1:n, t in rand(5)
-                b′₊₀ = bsplinebasis′₊₀(fP, i, t)
-                b′₋₀ = bsplinebasis′₋₀(fP, i, t)
-                b′   = bsplinebasis′(fP, i, t)
-                a′₊₀ = (bsplinebasis₊₀(fP, i, t+Δt) - bsplinebasis₊₀(fP, i, t-Δt)) / 2Δt
-                a′₋₀ = (bsplinebasis₋₀(fP, i, t+Δt) - bsplinebasis₋₀(fP, i, t-Δt)) / 2Δt
-                a′   = (bsplinebasis(fP, i, t+Δt) - bsplinebasis(fP, i, t-Δt)) / 2Δt
+                b′₊₀ = bsplinebasis′₊₀(P, i, t)
+                b′₋₀ = bsplinebasis′₋₀(P, i, t)
+                b′   = bsplinebasis′(P, i, t)
+                a′₊₀ = (bsplinebasis₊₀(P, i, t+Δt) - bsplinebasis₊₀(P, i, t-Δt)) / 2Δt
+                a′₋₀ = (bsplinebasis₋₀(P, i, t+Δt) - bsplinebasis₋₀(P, i, t-Δt)) / 2Δt
+                a′   = (bsplinebasis(P, i, t+Δt) - bsplinebasis(P, i, t-Δt)) / 2Δt
                 @test b′₊₀ ≈ b′₋₀ ≈ b′ ≈ a′₊₀ ≈ a′₋₀ ≈ a′
             end
         end
@@ -151,27 +133,22 @@ end
         v = rand(10)
         k = Knots(v) + (p+1)*Knots(0, 1) + Knots(v[2], v[3])
         P = BSplineSpace{p}(k)
-        fP = FastBSplineSpace(p, k)
         n = dim(P)
 
-        @test degree(P) == degree(fP) == 2
+        @test degree(P) == 2
 
         # Check the values of B-spline basis funcitons
         for i in 1:n, t in rand(5)
             for t in rand(5)
                 @test (bsplinebasis₊₀(P, i, t)
-                     ≈ bsplinebasis₊₀(fP, i, t)
                      ≈ bsplinebasis₋₀(P, i, t)
-                     ≈ bsplinebasis₋₀(fP, i, t)
-                     ≈ bsplinebasis(P, i, t)
-                     ≈ bsplinebasis(fP, i, t))
+                     ≈ bsplinebasis(P, i, t))
             end
             for t in k
                 t₊₀ = nextfloat(t)
                 t₋₀ = prevfloat(t)
-                @test bsplinebasis₊₀(P, i, t) ≈ bsplinebasis₊₀(fP, i, t) ≒ bsplinebasis₊₀(P, i, t₊₀) ≒ bsplinebasis₊₀(fP, i, t₊₀)
-                @test bsplinebasis₋₀(P, i, t) ≈ bsplinebasis₋₀(fP, i, t) ≒ bsplinebasis₋₀(P, i, t₋₀) ≒ bsplinebasis₋₀(fP, i, t₋₀)
-                @test bsplinebasis(P, i, t) ≈ bsplinebasis(fP, i, t)
+                @test bsplinebasis₊₀(P, i, t) ≒ bsplinebasis₊₀(P, i, t₊₀)
+                @test bsplinebasis₋₀(P, i, t)≒ bsplinebasis₋₀(P, i, t₋₀)
             end
         end
 
@@ -179,18 +156,14 @@ end
         for i in 1:n
             for t in rand(5)
                 @test (bsplinebasis′₊₀(P, i, t)
-                     ≈ bsplinebasis′₊₀(fP, i, t)
                      ≈ bsplinebasis′₋₀(P, i, t)
-                     ≈ bsplinebasis′₋₀(fP, i, t)
-                     ≈ bsplinebasis′(P, i, t)
-                     ≈ bsplinebasis′(fP, i, t))
+                     ≈ bsplinebasis′(P, i, t))
             end
             for t in k
                 t₊₀ = nextfloat(t)
                 t₋₀ = prevfloat(t)
-                @test bsplinebasis′₊₀(P, i, t) ≈ bsplinebasis′₊₀(fP, i, t) ≒ bsplinebasis′₊₀(P, i, t₊₀) ≒ bsplinebasis′₊₀(fP, i, t₊₀)
-                @test bsplinebasis′₋₀(P, i, t) ≈ bsplinebasis′₋₀(fP, i, t) ≒ bsplinebasis′₋₀(P, i, t₋₀) ≒ bsplinebasis′₋₀(fP, i, t₋₀)
-                @test bsplinebasis′(P, i, t) ≈ bsplinebasis′(fP, i, t)
+                @test bsplinebasis′₊₀(P, i, t) ≒ bsplinebasis′₊₀(P, i, t₊₀)
+                @test bsplinebasis′₋₀(P, i, t) ≒ bsplinebasis′₋₀(P, i, t₋₀)
             end
         end
 
@@ -198,12 +171,12 @@ end
         @testset "derivative" begin
             Random.seed!(57)
             for i in 1:n, t in rand(5)
-                b′₊₀ = bsplinebasis′₊₀(fP, i, t)
-                b′₋₀ = bsplinebasis′₋₀(fP, i, t)
-                b′   = bsplinebasis′(fP, i, t)
-                a′₊₀ = (bsplinebasis₊₀(fP, i, t+Δt) - bsplinebasis₊₀(fP, i, t-Δt)) / 2Δt
-                a′₋₀ = (bsplinebasis₋₀(fP, i, t+Δt) - bsplinebasis₋₀(fP, i, t-Δt)) / 2Δt
-                a′   = (bsplinebasis(fP, i, t+Δt) - bsplinebasis(fP, i, t-Δt)) / 2Δt
+                b′₊₀ = bsplinebasis′₊₀(P, i, t)
+                b′₋₀ = bsplinebasis′₋₀(P, i, t)
+                b′   = bsplinebasis′(P, i, t)
+                a′₊₀ = (bsplinebasis₊₀(P, i, t+Δt) - bsplinebasis₊₀(P, i, t-Δt)) / 2Δt
+                a′₋₀ = (bsplinebasis₋₀(P, i, t+Δt) - bsplinebasis₋₀(P, i, t-Δt)) / 2Δt
+                a′   = (bsplinebasis(P, i, t+Δt) - bsplinebasis(P, i, t-Δt)) / 2Δt
                 @test b′₊₀ ≈ b′₋₀ ≈ b′ ≈ a′₊₀ ≈ a′₋₀ ≈ a′
             end
         end
@@ -216,27 +189,22 @@ end
         v = rand(10)
         k = Knots(v) + (p+1)*Knots(0, 1) + Knots(v[2], v[3])
         P = BSplineSpace{p}(k)
-        fP = FastBSplineSpace(p, k)
         n = dim(P)
 
-        @test degree(P) == degree(fP) == 3
+        @test degree(P) == 3
 
         # Check the values of B-spline basis funcitons
         for i in 1:n, t in rand(5)
             for t in rand(5)
                 @test (bsplinebasis₊₀(P, i, t)
-                     ≈ bsplinebasis₊₀(fP, i, t)
                      ≈ bsplinebasis₋₀(P, i, t)
-                     ≈ bsplinebasis₋₀(fP, i, t)
-                     ≈ bsplinebasis(P, i, t)
-                     ≈ bsplinebasis(fP, i, t))
+                     ≈ bsplinebasis(P, i, t))
             end
             for t in k
                 t₊₀ = nextfloat(t)
                 t₋₀ = prevfloat(t)
-                @test bsplinebasis₊₀(P, i, t) ≈ bsplinebasis₊₀(fP, i, t) ≒ bsplinebasis₊₀(P, i, t₊₀) ≒ bsplinebasis₊₀(fP, i, t₊₀)
-                @test bsplinebasis₋₀(P, i, t) ≈ bsplinebasis₋₀(fP, i, t) ≒ bsplinebasis₋₀(P, i, t₋₀) ≒ bsplinebasis₋₀(fP, i, t₋₀)
-                @test bsplinebasis(P, i, t) ≈ bsplinebasis(fP, i, t)
+                @test bsplinebasis₊₀(P, i, t) ≒ bsplinebasis₊₀(P, i, t₊₀)
+                @test bsplinebasis₋₀(P, i, t)≒ bsplinebasis₋₀(P, i, t₋₀)
             end
         end
 
@@ -244,18 +212,14 @@ end
         for i in 1:n
             for t in rand(5)
                 @test (bsplinebasis′₊₀(P, i, t)
-                     ≈ bsplinebasis′₊₀(fP, i, t)
                      ≈ bsplinebasis′₋₀(P, i, t)
-                     ≈ bsplinebasis′₋₀(fP, i, t)
-                     ≈ bsplinebasis′(P, i, t)
-                     ≈ bsplinebasis′(fP, i, t))
+                     ≈ bsplinebasis′(P, i, t))
             end
             for t in k
                 t₊₀ = nextfloat(t)
                 t₋₀ = prevfloat(t)
-                @test bsplinebasis′₊₀(P, i, t) ≈ bsplinebasis′₊₀(fP, i, t) ≒ bsplinebasis′₊₀(P, i, t₊₀) ≒ bsplinebasis′₊₀(fP, i, t₊₀)
-                @test bsplinebasis′₋₀(P, i, t) ≈ bsplinebasis′₋₀(fP, i, t) ≒ bsplinebasis′₋₀(P, i, t₋₀) ≒ bsplinebasis′₋₀(fP, i, t₋₀)
-                @test bsplinebasis′(P, i, t) ≈ bsplinebasis′(fP, i, t)
+                @test bsplinebasis′₊₀(P, i, t) ≒ bsplinebasis′₊₀(P, i, t₊₀)
+                @test bsplinebasis′₋₀(P, i, t) ≒ bsplinebasis′₋₀(P, i, t₋₀)
             end
         end
 
@@ -263,12 +227,12 @@ end
         @testset "derivative" begin
             Random.seed!(57)
             for i in 1:n, t in rand(5)
-                b′₊₀ = bsplinebasis′₊₀(fP, i, t)
-                b′₋₀ = bsplinebasis′₋₀(fP, i, t)
-                b′   = bsplinebasis′(fP, i, t)
-                a′₊₀ = (bsplinebasis₊₀(fP, i, t+Δt) - bsplinebasis₊₀(fP, i, t-Δt)) / 2Δt
-                a′₋₀ = (bsplinebasis₋₀(fP, i, t+Δt) - bsplinebasis₋₀(fP, i, t-Δt)) / 2Δt
-                a′   = (bsplinebasis(fP, i, t+Δt) - bsplinebasis(fP, i, t-Δt)) / 2Δt
+                b′₊₀ = bsplinebasis′₊₀(P, i, t)
+                b′₋₀ = bsplinebasis′₋₀(P, i, t)
+                b′   = bsplinebasis′(P, i, t)
+                a′₊₀ = (bsplinebasis₊₀(P, i, t+Δt) - bsplinebasis₊₀(P, i, t-Δt)) / 2Δt
+                a′₋₀ = (bsplinebasis₋₀(P, i, t+Δt) - bsplinebasis₋₀(P, i, t-Δt)) / 2Δt
+                a′   = (bsplinebasis(P, i, t+Δt) - bsplinebasis(P, i, t-Δt)) / 2Δt
                 @test b′₊₀ ≈ b′₋₀ ≈ b′ ≈ a′₊₀ ≈ a′₋₀ ≈ a′
             end
         end
@@ -281,27 +245,22 @@ end
         v = rand(10)
         k = Knots(v) + (p+1)*Knots(0, 1) + Knots(v[2], v[3])
         P = BSplineSpace{p}(k)
-        fP = FastBSplineSpace(p, k)
         n = dim(P)
 
-        @test degree(P) == degree(fP) == 4
+        @test degree(P) == 4
 
         # Check the values of B-spline basis funcitons
         for i in 1:n, t in rand(5)
             for t in rand(5)
                 @test (bsplinebasis₊₀(P, i, t)
-                     ≈ bsplinebasis₊₀(fP, i, t)
                      ≈ bsplinebasis₋₀(P, i, t)
-                     ≈ bsplinebasis₋₀(fP, i, t)
-                     ≈ bsplinebasis(P, i, t)
-                     ≈ bsplinebasis(fP, i, t))
+                     ≈ bsplinebasis(P, i, t))
             end
             for t in k
                 t₊₀ = nextfloat(t)
                 t₋₀ = prevfloat(t)
-                @test bsplinebasis₊₀(P, i, t) ≈ bsplinebasis₊₀(fP, i, t) ≒ bsplinebasis₊₀(P, i, t₊₀) ≒ bsplinebasis₊₀(fP, i, t₊₀)
-                @test bsplinebasis₋₀(P, i, t) ≈ bsplinebasis₋₀(fP, i, t) ≒ bsplinebasis₋₀(P, i, t₋₀) ≒ bsplinebasis₋₀(fP, i, t₋₀)
-                @test bsplinebasis(P, i, t) ≈ bsplinebasis(fP, i, t)
+                @test bsplinebasis₊₀(P, i, t) ≒ bsplinebasis₊₀(P, i, t₊₀)
+                @test bsplinebasis₋₀(P, i, t)≒ bsplinebasis₋₀(P, i, t₋₀)
             end
         end
 
@@ -309,18 +268,14 @@ end
         for i in 1:n
             for t in rand(5)
                 @test (bsplinebasis′₊₀(P, i, t)
-                     ≈ bsplinebasis′₊₀(fP, i, t)
                      ≈ bsplinebasis′₋₀(P, i, t)
-                     ≈ bsplinebasis′₋₀(fP, i, t)
-                     ≈ bsplinebasis′(P, i, t)
-                     ≈ bsplinebasis′(fP, i, t))
+                     ≈ bsplinebasis′(P, i, t))
             end
             for t in k
                 t₊₀ = nextfloat(t)
                 t₋₀ = prevfloat(t)
-                @test bsplinebasis′₊₀(P, i, t) ≈ bsplinebasis′₊₀(fP, i, t) ≒ bsplinebasis′₊₀(P, i, t₊₀) ≒ bsplinebasis′₊₀(fP, i, t₊₀)
-                @test bsplinebasis′₋₀(P, i, t) ≈ bsplinebasis′₋₀(fP, i, t) ≒ bsplinebasis′₋₀(P, i, t₋₀) ≒ bsplinebasis′₋₀(fP, i, t₋₀)
-                @test bsplinebasis′(P, i, t) ≈ bsplinebasis′(fP, i, t)
+                @test bsplinebasis′₊₀(P, i, t) ≒ bsplinebasis′₊₀(P, i, t₊₀)
+                @test bsplinebasis′₋₀(P, i, t) ≒ bsplinebasis′₋₀(P, i, t₋₀)
             end
         end
 
@@ -328,12 +283,12 @@ end
         @testset "derivative" begin
             Random.seed!(57)
             for i in 1:n, t in rand(5)
-                b′₊₀ = bsplinebasis′₊₀(fP, i, t)
-                b′₋₀ = bsplinebasis′₋₀(fP, i, t)
-                b′   = bsplinebasis′(fP, i, t)
-                a′₊₀ = (bsplinebasis₊₀(fP, i, t+Δt) - bsplinebasis₊₀(fP, i, t-Δt)) / 2Δt
-                a′₋₀ = (bsplinebasis₋₀(fP, i, t+Δt) - bsplinebasis₋₀(fP, i, t-Δt)) / 2Δt
-                a′   = (bsplinebasis(fP, i, t+Δt) - bsplinebasis(fP, i, t-Δt)) / 2Δt
+                b′₊₀ = bsplinebasis′₊₀(P, i, t)
+                b′₋₀ = bsplinebasis′₋₀(P, i, t)
+                b′   = bsplinebasis′(P, i, t)
+                a′₊₀ = (bsplinebasis₊₀(P, i, t+Δt) - bsplinebasis₊₀(P, i, t-Δt)) / 2Δt
+                a′₋₀ = (bsplinebasis₋₀(P, i, t+Δt) - bsplinebasis₋₀(P, i, t-Δt)) / 2Δt
+                a′   = (bsplinebasis(P, i, t+Δt) - bsplinebasis(P, i, t-Δt)) / 2Δt
                 @test b′₊₀ ≈ b′₋₀ ≈ b′ ≈ a′₊₀ ≈ a′₋₀ ≈ a′
             end
         end
@@ -346,34 +301,30 @@ end
         ts = rand(10)
 
         for p in 0:BasicBSpline.MAX_DEGREE
-            fP = FastBSplineSpace(p,k)
             P = BSplineSpace{p}(k)
             for t in ts
-                j = BasicBSpline._knotindex(fP,t)
-                J = intervalindex(P,t)
-                _B = BasicBSpline._bsplinebasis(fP,t,j)
-                _b = bsplinebasisall(P,J,t)
-                # _B′ = BasicBSpline._bsplinebasis′(fP,t,j)
-                B = [bsplinebasis(fP,i,t) for i in j-p:j]
-                B′ = [bsplinebasis′(fP,i,t) for i in j-p:j]
+                j = intervalindex(P,t)
+                _B = bsplinebasisall(P,j,t)
+                B = [bsplinebasis(P,i,t) for i in j:j+p]
+                B′ = [bsplinebasis′(P,i,t) for i in j:j+p]
                 @test norm(collect(_B) - B) < ε
 
-                _B = [bsplinebasis(P,i,t) for i in j-p:j]
+                _B = [bsplinebasis(P,i,t) for i in j:j+p]
                 @test norm(_B - B) < ε
 
-                _B = [bsplinebasis₊₀(P,i,t) for i in j-p:j]
+                _B = [bsplinebasis₊₀(P,i,t) for i in j:j+p]
                 @test norm(_B - B) < ε
 
-                _B = [bsplinebasis₋₀(P,i,t) for i in j-p:j]
+                _B = [bsplinebasis₋₀(P,i,t) for i in j:j+p]
                 @test norm(_B - B) < ε
 
-                _B′ =[ bsplinebasis′(P,i,t) for i in j-p:j]
+                _B′ =[ bsplinebasis′(P,i,t) for i in j:j+p]
                 @test norm(_B′ - B′) < ε
 
-                _B′ =[ bsplinebasis′₊₀(P,i,t) for i in j-p:j]
+                _B′ =[ bsplinebasis′₊₀(P,i,t) for i in j:j+p]
                 @test norm(_B′ - B′) < ε
 
-                _B′ =[ bsplinebasis′₋₀(P,i,t) for i in j-p:j]
+                _B′ =[ bsplinebasis′₋₀(P,i,t) for i in j:j+p]
                 @test norm(_B′ - B′) < ε
             end
         end
diff --git a/test/test_BSplineManifold.jl b/test/test_BSplineManifold.jl
index 8c8ac98e..0a3d247e 100644
--- a/test/test_BSplineManifold.jl
+++ b/test/test_BSplineManifold.jl
@@ -45,22 +45,22 @@
             end
         end
 
-        @testset "FastBSplineManifold-1dim" begin
+        @testset "CustomBSplineManifold-1dim" begin
             Random.seed!(42)
 
-            P1 = FastBSplineSpace(1, Knots([0, 0, 1, 1]))
+            P1 = BSplineSpace{1}(Knots([0, 0, 1, 1]))
             n1 = dim(P1) # 2
             a = [Point(i, rand()) for i in 1:n1]  # n1 × n2 array of d̂-dim vector.
-            M = FastBSplineManifold([P1], a)
+            M = CustomBSplineManifold((P1,), a)
             @test dim(M) == 1
 
-            P1′ = FastBSplineSpace(2, Knots([-2, 0, 0, 1, 1, 2]))
-            p₊ = [1]
-            k₊ = [Knots()]
+            P1′ = BSplineSpace{2}(Knots([-2, 0, 0, 1, 1, 2]))
+            p₊ = (1,)
+            k₊ = (Knots(),)
 
             @test P1 ⊑ P1′
 
-            M′ = refinement(M, [P1′])
+            M′ = refinement(M, (P1′,))
             M′′ = refinement(M, p₊=p₊, k₊=k₊)
             ts = [[rand()] for _ in 1:10]
             for t in ts
@@ -68,30 +68,6 @@
                 @test M(t...) ≈ M′′(t...)
             end
         end
-
-        @testset "BSplineCurve" begin
-            Random.seed!(42)
-
-            P1 = FastBSplineSpace(1, Knots([0, 0, 1, 1]))
-            n1 = dim(P1) # 2
-            a = [Point(i, rand()) for i in 1:n1]  # n1 × n2 array of d̂-dim vector.
-            M = BSplineCurve([P1], a)
-            @test dim(M) == 1
-
-            P1′ = FastBSplineSpace(2, Knots([-2, 0, 0, 1, 1, 2]))
-            p₊ = [1]
-            k₊ = [Knots()]
-
-            @test P1 ⊑ P1′
-
-            M′ = refinement(M, [P1′])
-            M′′ = refinement(M, p₊=p₊, k₊=k₊)
-            ts = [[rand()] for _ in 1:10]
-            for t in ts
-                @test M(t) ≈ M′(t)
-                @test M(t) ≈ M′′(t)
-            end
-        end
     end
 
     @testset "2dim" begin
@@ -124,26 +100,26 @@
             end
         end
 
-        @testset "FastBSplineManifold-2dim" begin
+        @testset "CustomBSplineManifold-2dim" begin
             Random.seed!(42)
 
-            P1 = FastBSplineSpace(1, Knots([0, 0, 1, 1]))
-            P2 = FastBSplineSpace(1, Knots([1, 1, 2, 3, 3]))
+            P1 = BSplineSpace{1}(Knots([0, 0, 1, 1]))
+            P2 = BSplineSpace{1}(Knots([1, 1, 2, 3, 3]))
             n1 = dim(P1) # 2
             n2 = dim(P2) # 3
             a = [Point(i, j) for i in 1:n1, j in 1:n2]  # n1 × n2 array of d̂-dim vector.
-            M = FastBSplineManifold([P1, P2], a)
+            M = CustomBSplineManifold((P1, P2), a)
             @test dim(M) == 2
 
-            P1′ = FastBSplineSpace(2, Knots([0, 0, 0, 1, 1, 1]))
-            P2′ = FastBSplineSpace(1, Knots([1, 1, 2, 1.45, 3, 3]))
-            p₊ = [1, 0]
-            k₊ = [Knots(), Knots(1.45)]
+            P1′ = BSplineSpace{2}(Knots([0, 0, 0, 1, 1, 1]))
+            P2′ = BSplineSpace{1}(Knots([1, 1, 2, 1.45, 3, 3]))
+            p₊ = (1, 0)
+            k₊ = (Knots(), Knots(1.45))
 
             @test P1 ⊆ P1′
             @test P2 ⊆ P2′
 
-            M′ = refinement(M, [P1′, P2′])
+            M′ = refinement(M, (P1′, P2′))
             M′′ = refinement(M, p₊=p₊, k₊=k₊)
             ts = [[rand(), 1 + 2 * rand()] for _ in 1:10]
             for t in ts
@@ -151,33 +127,5 @@
                 @test M(t...) ≈ M′′(t...)
             end
         end
-
-        @testset "BSplineSurface" begin
-            Random.seed!(42)
-
-            P1 = FastBSplineSpace(1, Knots([0, 0, 1, 1]))
-            P2 = FastBSplineSpace(1, Knots([1, 1, 2, 3, 3]))
-            n1 = dim(P1) # 2
-            n2 = dim(P2) # 3
-            a = [Point(i, j) for i in 1:n1, j in 1:n2]  # n1 × n2 array of d̂-dim vector.
-            M = BSplineSurface([P1, P2], a)
-            @test dim(M) == 2
-
-            P1′ = FastBSplineSpace(2, Knots([0, 0, 0, 1, 1, 1]))
-            P2′ = FastBSplineSpace(1, Knots([1, 1, 2, 1.45, 3, 3]))
-            p₊ = [1, 0]
-            k₊ = [Knots(), Knots(1.45)]
-
-            @test P1 ⊆ P1′
-            @test P2 ⊆ P2′
-
-            M′ = refinement(M, [P1′, P2′])
-            M′′ = refinement(M, p₊=p₊, k₊=k₊)
-            ts = [[rand(), 1 + 2 * rand()] for _ in 1:10]
-            for t in ts
-                @test M(t) ≈ M′(t)
-                @test M(t) ≈ M′′(t)
-            end
-        end
     end
 end
diff --git a/test/test_FastBSplineSpace.jl b/test/test_FastBSplineSpace.jl
deleted file mode 100644
index e294aea4..00000000
--- a/test/test_FastBSplineSpace.jl
+++ /dev/null
@@ -1,31 +0,0 @@
-# TODO: remove these tests
-@testset "FastBSplineSpace" begin
-    P1 = FastBSplineSpace(2, Knots([1, 3, 5, 6, 8, 9]))
-    @test bsplinesupport(P1, 2) == 3..8
-    @test dim(P1) == 3
-    @test properdim(P1) == 3
-    @test isproper(P1) == true
-
-    P2 = FastBSplineSpace(1, Knots([1, 3, 3, 3, 8, 9]))
-    @test isproper(P2) == false
-    @test properdim(P2) == 3
-
-    i = 2
-    k = Knots([5, 12, 13, 13, 14])
-    p = 2
-    P = FastBSplineSpace(p, k)
-    @test bsplinesupport(P) == [5..13, 12..14]
-    @test bsplinesupport(P, i) == 12..14
-
-    @test isproper(FastBSplineSpace(2, Knots([1, 3, 5, 6, 8, 9])))
-    @test !isproper(FastBSplineSpace(1, Knots([1, 3, 3, 3, 8, 9])))
-
-    @test dim(FastBSplineSpace(2, Knots([1, 3, 5, 6, 8, 9]))) == 3
-
-    P1 = FastBSplineSpace(1, Knots([1, 3, 5, 8]))
-    P2 = FastBSplineSpace(1, Knots([1, 3, 5, 6, 8, 9]))
-    P3 = FastBSplineSpace(2, Knots([1, 1, 3, 3, 5, 5, 8, 8]))
-    @test P1 ⊆ P2
-    @test P1 ⊆ P3
-    @test P2 ⊈ P3
-end
diff --git a/test/test_Fitting.jl b/test/test_Fitting.jl
index 53a91095..b535aa39 100644
--- a/test/test_Fitting.jl
+++ b/test/test_Fitting.jl
@@ -6,22 +6,22 @@
 
         p1 = 2
         k1 = Knots(rand(3)) + (p1 + 1) * Knots([0, 1])
-        P1 = FastBSplineSpace(p1, k1)
+        P1 = BSplineSpace{p1}(k1)
         n1 = dim(P1)
         a_org = [Point(i1, rand()) for i1 in 1:n1]
-        M = BSplineCurve([P1], a_org)
+        M = CustomBSplineManifold((P1,), a_org)
 
         p1′ = p1 + 1
         k1′ = k1 + unique(k1) + Knots(rand(2))
-        P1′ = FastBSplineSpace(p1′, k1′)
+        P1′ = BSplineSpace{p1′}(k1′)
 
-        M′ = refinement(M, [P1′])
+        M′ = refinement(M, (P1′,))
         a_ref = M′.controlpoints
 
-        a_fit = fittingcontrolpoints(M, [P1′])
+        a_fit = fittingcontrolpoints(M, (P1′,))
         @test norm(a_fit - a_ref) < ε
 
-        a_fit = fittingcontrolpoints(M, [P1′], domain=:R)
+        a_fit = fittingcontrolpoints(M, (P1′,), domain=:R)
         @test norm(a_fit - a_ref) < ε
     end
 
@@ -30,19 +30,19 @@
 
         p1 = 2
         k1 = Knots(rand(3)) + Knots([0, 1]) + p1 * Knots([-rand(), 1 + rand()])
-        P1 = FastBSplineSpace(p1, k1)
+        P1 = BSplineSpace{p1}(k1)
         n1 = dim(P1)
         a_org = [Point(i1, rand()) for i1 in 1:n1]
-        M = BSplineCurve([P1], a_org)
+        M = CustomBSplineManifold((P1,), a_org)
 
         p1′ = p1 + 1
         k1′ = k1 + unique(k1[1+p1:end-p1]) + Knots(rand(2))
-        P1′ = FastBSplineSpace(p1′, k1′)
+        P1′ = BSplineSpace{p1′}(k1′)
 
-        M′ = refinement(M, [P1′])
+        M′ = refinement(M, (P1′,))
         a_ref = M′.controlpoints
 
-        a_fit = fittingcontrolpoints(M, [P1′])
+        a_fit = fittingcontrolpoints(M, (P1′,))
 
         @test norm(a_fit - a_ref) < ε
     end
@@ -52,29 +52,29 @@
 
         p1 = 2
         k1 = Knots(rand(3)) + (p1 + 1) * Knots([0, 1])
-        P1 = FastBSplineSpace(p1, k1)
+        P1 = BSplineSpace{p1}(k1)
         n1 = dim(P1)
         p2 = 1
         k2 = Knots(rand(4)) + (p2 + 1) * Knots([0, 1])
-        P2 = FastBSplineSpace(p2, k2)
+        P2 = BSplineSpace{p2}(k2)
         n2 = dim(P2)
         a_org = [Point(i1, i2, rand()) for i1 in 1:n1, i2 in 1:n2]
-        M = BSplineSurface([P1, P2], a_org)
+        M = CustomBSplineManifold((P1, P2), a_org)
 
         p1′ = p1 + 1
         k1′ = k1 + unique(k1) + Knots(rand(2))
-        P1′ = FastBSplineSpace(p1′, k1′)
+        P1′ = BSplineSpace{p1′}(k1′)
         p2′ = p2 + 1
         k2′ = k2 + unique(k2) + Knots(rand(2))
-        P2′ = FastBSplineSpace(p2′, k2′)
+        P2′ = BSplineSpace{p2′}(k2′)
 
-        M′ = refinement(M, [P1′, P2′])
+        M′ = refinement(M, (P1′, P2′))
         a_ref = M′.controlpoints
 
-        a_fit = fittingcontrolpoints(M, [P1′, P2′])
+        a_fit = fittingcontrolpoints(M, (P1′, P2′))
         @test norm(a_fit - a_ref) < ε
 
-        a_fit = fittingcontrolpoints(M, [P1′, P2′], domain=:R)
+        a_fit = fittingcontrolpoints(M, (P1′, P2′), domain=:R)
         @test norm(a_fit - a_ref) < ε
     end
 
@@ -83,26 +83,26 @@
 
         p1 = 2
         k1 = Knots(rand(3)) + Knots([0, 1]) + p1 * Knots([-rand(), 1 + rand()])
-        P1 = FastBSplineSpace(p1, k1)
+        P1 = BSplineSpace{p1}(k1)
         n1 = dim(P1)
         p2 = 1
         k2 = Knots(rand(4)) + Knots([0, 1]) + p2 * Knots([-rand(), 1 + rand()])
-        P2 = FastBSplineSpace(p2, k2)
+        P2 = BSplineSpace{p2}(k2)
         n2 = dim(P2)
         a_org = [Point(i1, i2, rand()) for i1 in 1:n1, i2 in 1:n2]
-        M = BSplineSurface([P1, P2], a_org)
+        M = CustomBSplineManifold((P1, P2), a_org)
 
         p1′ = p1 + 1
         k1′ = k1 + unique(k1[1+p1:end-p1]) + Knots(rand(2))
-        P1′ = FastBSplineSpace(p1′, k1′)
+        P1′ = BSplineSpace{p1′}(k1′)
         p2′ = p2 + 1
         k2′ = k2 + unique(k2[1+p2:end-p2]) + Knots(rand(2))
-        P2′ = FastBSplineSpace(p2′, k2′)
+        P2′ = BSplineSpace{p2′}(k2′)
 
-        M′ = refinement(M, [P1′, P2′])
+        M′ = refinement(M, (P1′, P2′))
         a_ref = M′.controlpoints
 
-        a_fit = fittingcontrolpoints(M, [P1′, P2′])
+        a_fit = fittingcontrolpoints(M, (P1′, P2′))
         @test norm(a_fit - a_ref) < ε
     end
 
@@ -111,36 +111,36 @@
 
         p1 = 2
         k1 = Knots(rand(3)) + (p1 + 1) * Knots([0, 1])
-        P1 = FastBSplineSpace(p1, k1)
+        P1 = BSplineSpace{p1}(k1)
         n1 = dim(P1)
         p2 = 1
         k2 = Knots(rand(4)) + (p2 + 1) * Knots([0, 1])
-        P2 = FastBSplineSpace(p2, k2)
+        P2 = BSplineSpace{p2}(k2)
         n2 = dim(P2)
         p3 = 2
         k3 = Knots(rand(5)) + (p3 + 1) * Knots([0, 1])
-        P3 = FastBSplineSpace(p3, k3)
+        P3 = BSplineSpace{p3}(k3)
         n3 = dim(P3)
         a_org = [Point(i1, i2, i3, rand()) for i1 in 1:n1, i2 in 1:n2, i3 in 1:n3]
-        M = BSplineSolid([P1, P2, P3], a_org)
+        M = CustomBSplineManifold((P1, P2, P3), a_org)
 
         p1′ = p1 + 1
         k1′ = k1 + unique(k1) + Knots(rand(2))
-        P1′ = FastBSplineSpace(p1′, k1′)
+        P1′ = BSplineSpace{p1′}(k1′)
         p2′ = p2 + 1
         k2′ = k2 + unique(k2) + Knots(rand(2))
-        P2′ = FastBSplineSpace(p2′, k2′)
+        P2′ = BSplineSpace{p2′}(k2′)
         p3′ = p3 + 1
         k3′ = k3 + unique(k3) + Knots(rand(2))
-        P3′ = FastBSplineSpace(p3′, k3′)
+        P3′ = BSplineSpace{p3′}(k3′)
 
-        M′ = refinement(M, [P1′, P2′, P3′])
+        M′ = refinement(M, (P1′, P2′, P3′))
         a_ref = M′.controlpoints
 
-        a_fit = fittingcontrolpoints(M, [P1′, P2′, P3′])
+        a_fit = fittingcontrolpoints(M, (P1′, P2′, P3′))
         @test norm(a_fit - a_ref) < ε
 
-        a_fit = fittingcontrolpoints(M, [P1′, P2′, P3′], domain=:R)
+        a_fit = fittingcontrolpoints(M, (P1′, P2′, P3′), domain=:R)
         @test norm(a_fit - a_ref) < ε
     end
 
@@ -149,33 +149,33 @@
 
         p1 = 2
         k1 = Knots(rand(3)) + Knots([0, 1]) + p1 * Knots([-rand(), 1 + rand()])
-        P1 = FastBSplineSpace(p1, k1)
+        P1 = BSplineSpace{p1}(k1)
         n1 = dim(P1)
         p2 = 1
         k2 = Knots(rand(4)) + Knots([0, 1]) + p1 * Knots([-rand(), 1 + rand()])
-        P2 = FastBSplineSpace(p2, k2)
+        P2 = BSplineSpace{p2}(k2)
         n2 = dim(P2)
         p3 = 2
         k3 = Knots(rand(5)) + Knots([0, 1]) + p1 * Knots([-rand(), 1 + rand()])
-        P3 = FastBSplineSpace(p3, k3)
+        P3 = BSplineSpace{p3}(k3)
         n3 = dim(P3)
         a_org = [Point(i1, i2, i3, rand()) for i1 in 1:n1, i2 in 1:n2, i3 in 1:n3]
-        M = BSplineSolid([P1, P2, P3], a_org)
+        M = CustomBSplineManifold((P1, P2, P3), a_org)
 
         p1′ = p1 + 1
         k1′ = k1 + unique(k1[1+p1:end-p1]) + Knots(rand(2))
-        P1′ = FastBSplineSpace(p1′, k1′)
+        P1′ = BSplineSpace{p1′}(k1′)
         p2′ = p2 + 1
         k2′ = k2 + unique(k2[1+p2:end-p2]) + Knots(rand(2))
-        P2′ = FastBSplineSpace(p2′, k2′)
+        P2′ = BSplineSpace{p2′}(k2′)
         p3′ = p3 + 1
         k3′ = k3 + unique(k3[1+p3:end-p3]) + Knots(rand(2))
-        P3′ = FastBSplineSpace(p3′, k3′)
+        P3′ = BSplineSpace{p3′}(k3′)
 
-        M′ = refinement(M, [P1′, P2′, P3′])
+        M′ = refinement(M, (P1′, P2′, P3′))
         a_ref = M′.controlpoints
 
-        a_fit = fittingcontrolpoints(M, [P1′, P2′, P3′])
+        a_fit = fittingcontrolpoints(M, (P1′, P2′, P3′))
         @test norm(a_fit - a_ref) < 1.0e-6
     end
 end