Repository: hyrodium/BasicBSpline.jl (julia)

Issue:
Add support for NURBS (`RationalBSplineManifold`)
Like this:
```julia
struct RationalBSplineManifold{Dim,Deg,T,S<:Tuple,Dim₊₁} <: AbstractRationalBSplineManifold{Dim,Deg}
    bsplinespaces::S
    controlpoints::Array{T,Dim₊₁}
    weights::Array{T,Dim}
end

struct CustomRationalBSplineManifold{Dim,Deg,C,S<:Tuple} <: AbstractRationalBSplineManifold{Dim,Deg}
    bsplinespaces::S
    controlpoints::Array{C,Dim}
    weights::Array{T,Dim}
end
```

Pull request:
This PR fixes #142.

- [x] `RationalBSplineManifold` constructor
- [x] Add `bsplinespaces` methods etc.
- [x] Add `mapping` methods
- [x] Add `refinement` methods

Patch:
diff --git a/src/BasicBSpline.jl b/src/BasicBSpline.jl
index c70e61364..c471a4360 100644
--- a/src/BasicBSpline.jl
+++ b/src/BasicBSpline.jl
@@ -9,12 +9,13 @@ export AbstractKnotVector, KnotVector, UniformKnotVector
 export BSplineSpace, UniformBSplineSpace
 export dim, ⊑, ⊒, ⋢, ⋣, ≃
 export BSplineDerivativeSpace
-export BSplineManifold
+export BSplineManifold, RationalBSplineManifold
 export bsplinebasis₊₀, bsplinebasis₋₀, bsplinebasis
 export bsplinebasis′₊₀, bsplinebasis′₋₀, bsplinebasis′
 export bsplinesupport, domain, changebasis, expandspace
 export bsplinebasisall, intervalindex
-export refinement, bsplinespaces, controlpoints
+export bsplinespaces, controlpoints, weights
+export refinement
 export isdegenerate, isnondegenerate, exactdim
 export degree, knotvector
 export 𝔫
@@ -32,6 +33,7 @@ include("_DerivativeSpace.jl")
 include("_DerivativeBasis.jl")
 include("_ChangeBasis.jl")
 include("_BSplineManifold.jl")
+include("_RationalBSplineManifold.jl")
 include("_Refinement.jl")
 include("_Fitting.jl")
 
diff --git a/src/_BSplineManifold.jl b/src/_BSplineManifold.jl
index 7295ff8bc..18b843db0 100644
--- a/src/_BSplineManifold.jl
+++ b/src/_BSplineManifold.jl
@@ -1,6 +1,7 @@
 # B-spline manifold
 
-abstract type AbstractBSplineManifold{Dim, Deg} end
+abstract type AbstractManifold{Dim} end
+abstract type AbstractBSplineManifold{Dim, Deg} <: AbstractManifold{Dim} end
 
 dim(::AbstractBSplineManifold{Dim}) where Dim = Dim
 
@@ -25,7 +26,7 @@ end
 bsplinespaces(M::BSplineManifold) = M.bsplinespaces
 controlpoints(M::BSplineManifold) = M.controlpoints
 
-@generated function unsafe_mapping(M::BSplineManifold{1,Deg,S,T},t1::Real) where {Deg,S,T}
+@generated function unsafe_mapping(M::BSplineManifold{1,Deg},t1::Real) where {Deg}
     p1, = Deg
     exs = Expr[]
     for j1 in 1:p1
@@ -42,7 +43,7 @@ controlpoints(M::BSplineManifold) = M.controlpoints
     )
 end
 
-@generated function unsafe_mapping(M::BSplineManifold{2,Deg,S,T},t1,t2) where {Deg,S,T}
+@generated function unsafe_mapping(M::BSplineManifold{2,Deg},t1::Real,t2::Real) where {Deg}
     p1, p2 = Deg
     exs = Expr[]
     for j2 in 1:p2+1, j1 in 1:p1+1
@@ -61,7 +62,7 @@ end
     )
 end
 
-@generated function unsafe_mapping(M::BSplineManifold{3,Deg,S,T},t1,t2,t3) where {Deg,S,T}
+@generated function unsafe_mapping(M::BSplineManifold{3,Deg},t1::Real,t2::Real,t3::Real) where {Deg}
     p1, p2, p3 = Deg
     exs = Expr[]
     for j3 in 1:p3+1, j2 in 1:p2+1, j1 in 1:p1+1
@@ -80,20 +81,20 @@ end
     )
 end
 
-@inline function (M::AbstractBSplineManifold{1})(t1)
+@inline function (M::AbstractManifold{1})(t1)
     Ps = bsplinespaces(M)
     t1 in domain(Ps[1]) || throw(DomainError(t1, "The input $(t1) is out of range."))
     unsafe_mapping(M,t1)
 end
 
-@inline function (M::AbstractBSplineManifold{2})(t1,t2)
+@inline function (M::AbstractManifold{2})(t1,t2)
     Ps = bsplinespaces(M)
     t1 in domain(Ps[1]) || throw(DomainError(t1, "The input $(t1) is out of range."))
     t2 in domain(Ps[2]) || throw(DomainError(t2, "The input $(t2) is out of range."))
     unsafe_mapping(M,t1,t2)
 end
 
-@inline function (M::AbstractBSplineManifold{3})(t1,t2,t3)
+@inline function (M::AbstractManifold{3})(t1,t2,t3)
     Ps = bsplinespaces(M)
     t1 in domain(Ps[1]) || throw(DomainError(t1, "The input $(t1) is out of range."))
     t2 in domain(Ps[2]) || throw(DomainError(t2, "The input $(t2) is out of range."))
diff --git a/src/_BSplineSpace.jl b/src/_BSplineSpace.jl
index a1f743a50..cf2a33d19 100644
--- a/src/_BSplineSpace.jl
+++ b/src/_BSplineSpace.jl
@@ -1,6 +1,7 @@
 # B-spline space
 
-abstract type AbstractBSplineSpace{p,T} end
+abstract type AbstractFunctionSpace{T} end
+abstract type AbstractBSplineSpace{p,T} <: AbstractFunctionSpace{T} end
 
 # Broadcast like a scalar
 Base.Broadcast.broadcastable(P::AbstractBSplineSpace) = Ref(P)
diff --git a/src/_ChangeBasis.jl b/src/_ChangeBasis.jl
index 76f45dfb6..d735b5392 100644
--- a/src/_ChangeBasis.jl
+++ b/src/_ChangeBasis.jl
@@ -151,7 +151,7 @@ function _changebasis_I(P::AbstractBSplineSpace{p,T}, P′::AbstractBSplineSpace
     return A
 end
 
-function changebasis(P, P′)
+function changebasis(P::AbstractFunctionSpace, P′::AbstractFunctionSpace)
     P ⊆ P′ && return _changebasis_R(P, P′)
     P ⊑ P′ && return _changebasis_I(P, P′)
     throw(DomainError((P, P′),"𝒫[p,k] ⊆ 𝒫[p′,k′] or 𝒫[p,k] ⊑ 𝒫[p′,k′] must hold."))
diff --git a/src/_DerivativeSpace.jl b/src/_DerivativeSpace.jl
index 6990f3c59..a066266fb 100644
--- a/src/_DerivativeSpace.jl
+++ b/src/_DerivativeSpace.jl
@@ -1,9 +1,9 @@
 # Space of derivative of B-spline basis function
 
-struct BSplineDerivativeSpace{r, T<:AbstractBSplineSpace}
-    bsplinespace::T
-    function BSplineDerivativeSpace{r}(P::T) where {r, T<:AbstractBSplineSpace}
-        new{r,T}(P)
+struct BSplineDerivativeSpace{r, S<:AbstractBSplineSpace, T} <: AbstractFunctionSpace{T}
+    bsplinespace::S
+    function BSplineDerivativeSpace{r}(P::S) where {r, S<:AbstractBSplineSpace{p,T}} where {p,T}
+        new{r,S,T}(P)
     end
 end
 
diff --git a/src/_RationalBSplineManifold.jl b/src/_RationalBSplineManifold.jl
new file mode 100644
index 000000000..636133233
--- /dev/null
+++ b/src/_RationalBSplineManifold.jl
@@ -0,0 +1,91 @@
+# Rational B-spline manifold (NURBS)
+abstract type AbstractRationalBSplineManifold{Dim, Deg} <: AbstractManifold{Dim} end
+
+struct RationalBSplineManifold{Dim,Deg,C,S<:Tuple,T} <: AbstractRationalBSplineManifold{Dim,Deg}
+    bsplinespaces::S
+    controlpoints::Array{C,Dim}
+    weights::Array{T,Dim}
+    function RationalBSplineManifold(a::Array{C,Dim},w::Array{T,Dim},Ps::S) where {S<:Tuple,C,Dim,T<:Real}
+        for P in Ps
+            if !(P isa AbstractBSplineSpace)
+                throw(TypeError(:BSplineManifold,AbstractBSplineSpace,P))
+            end
+        end
+        if size(a) != dim.(Ps)
+            msg = "The size of control points array $(size(a)) and dimensions of B-spline spaces $(dim.(Ps)) must be equal."
+            throw(DimensionMismatch(msg))
+        end
+        Deg = degree.(Ps)
+        new{Dim,Deg,C,S,T}(Ps,a,w)
+    end
+end
+
+controlpoints(M::RationalBSplineManifold) = M.controlpoints
+weights(M::RationalBSplineManifold) = M.weights
+bsplinespaces(M::RationalBSplineManifold) = M.bsplinespaces
+
+@generated function unsafe_mapping(M::RationalBSplineManifold{1,Deg},t1::Real) where {Deg}
+    p1, = Deg
+    exs = Expr[]
+    for j1 in 1:p1
+        push!(exs, :(u += b1[$(1+j1)]*getindex(w,i1+$(j1))))
+        push!(exs, :(v += b1[$(1+j1)]*getindex(w,i1+$(j1))*getindex(a,i1+$(j1))))
+    end
+    Expr(:block,
+        :((P1,) = bsplinespaces(M)),
+        :(a = controlpoints(M)),
+        :(w = weights(M)),
+        :(i1 = intervalindex(P1,t1)),
+        :(b1 = bsplinebasisall(P1,i1,t1)),
+        :(u = b1[1]*getindex(w,i1)),
+        :(v = b1[1]*getindex(w,i1)*getindex(a,i1)),
+        exs...,
+        :(return v/u)
+    )
+end
+
+@generated function unsafe_mapping(M::RationalBSplineManifold{2,Deg},t1::Real,t2::Real) where {Deg}
+    p1, p2 = Deg
+    exs = Expr[]
+    for j2 in 1:p2+1, j1 in 1:p1+1
+        push!(exs, :(u += b1[$(j1)]*b2[$(j2)]*getindex(w,i1+$(j1-1),i2+$(j2-1))))
+        push!(exs, :(v += b1[$(j1)]*b2[$(j2)]*getindex(w,i1+$(j1-1),i2+$(j2-1))*getindex(a,i1+$(j1-1),i2+$(j2-1))))
+    end
+    deleteat!(exs,1)
+    deleteat!(exs,1)
+    Expr(
+        :block,
+        :((P1, P2) = bsplinespaces(M)),
+        :(a = controlpoints(M)),
+        :(w = weights(M)),
+        :((i1, i2) = (intervalindex(P1,t1), intervalindex(P2,t2))),
+        :((b1, b2) = (bsplinebasisall(P1,i1,t1), bsplinebasisall(P2,i2,t2))),
+        :(u = b1[1]*b2[1]*getindex(w,i1,i2)),
+        :(v = b1[1]*b2[1]*getindex(w,i1,i2)*getindex(a,i1,i2)),
+        exs...,
+        :(return v/u)
+    )
+end
+
+@generated function unsafe_mapping(M::RationalBSplineManifold{3,Deg},t1::Real,t2::Real,t3::Real) where {Deg}
+    p1, p2, p3 = Deg
+    exs = Expr[]
+    for j3 in 1:p3+1, j2 in 1:p2+1, j1 in 1:p1+1
+        push!(exs, :(u += b1[$(j1)]*b2[$(j2)]*b3[$(j3)]*getindex(w,i1+$(j1-1),i2+$(j2-1),i3+$(j3-1))))
+        push!(exs, :(v += b1[$(j1)]*b2[$(j2)]*b3[$(j3)]*getindex(w,i1+$(j1-1),i2+$(j2-1),i3+$(j3-1))*getindex(a,i1+$(j1-1),i2+$(j2-1),i3+$(j3-1))))
+    end
+    deleteat!(exs,1)
+    deleteat!(exs,1)
+    Expr(
+        :block,
+        :((P1, P2, P3) = bsplinespaces(M)),
+        :(a = controlpoints(M)),
+        :(w = weights(M)),
+        :((i1, i2, i3) = (intervalindex(P1,t1), intervalindex(P2,t2), intervalindex(P3,t3))),
+        :((b1, b2, b3) = (bsplinebasisall(P1,i1,t1), bsplinebasisall(P2,i2,t2), bsplinebasisall(P3,i3,t3))),
+        :(u = b1[1]*b2[1]*b3[1]*getindex(w,i1,i2,i3)),
+        :(v = b1[1]*b2[1]*b3[1]*getindex(w,i1,i2,i3)*getindex(a,i1,i2,i3)),
+        exs...,
+        :(return v/u)
+    )
+end
diff --git a/src/_Refinement.jl b/src/_Refinement.jl
index b5720f8d3..af9b19506 100644
--- a/src/_Refinement.jl
+++ b/src/_Refinement.jl
@@ -39,11 +39,48 @@ function refinement(M::BSplineManifold{3}, Ps′::Tuple{<:AbstractBSplineSpace,<
     return BSplineManifold(a′, Ps′)
 end
 
+function refinement(M::RationalBSplineManifold{1}, Ps′::Tuple{<:AbstractBSplineSpace})
+    Ps = bsplinespaces(M)
+    a = controlpoints(M)
+    w = weights(M)
+    (n1,) = dim.(Ps)
+    (n1′,) = dim.(Ps′)
+    (A1,) = changebasis.(Ps, Ps′)
+
+    w′ = [sum(A1[I₁,J₁] * w[I₁] for I₁ in 1:n1) for J₁ in 1:n1′]
+    a′ = [sum(A1[I₁,J₁] * a[I₁] * w[I₁] for I₁ in 1:n1) for J₁ in 1:n1′] ./ w′
+    return RationalBSplineManifold(a′, w′, Ps′)
+end
+function refinement(M::RationalBSplineManifold{2}, Ps′::Tuple{<:AbstractBSplineSpace,<:AbstractBSplineSpace})
+    Ps = bsplinespaces(M)
+    a = controlpoints(M)
+    w = weights(M)
+    (n1,n2) = dim.(Ps)
+    (n1′,n2′) = dim.(Ps′)
+    (A1,A2) = changebasis.(Ps, Ps′)
+
+    w′ = [sum(A1[I₁,J₁] * A2[I₂,J₂] * w[I₁,I₂] for I₁ in 1:n1, I₂ in 1:n2) for J₁ in 1:n1′, J₂ in 1:n2′]
+    a′ = [sum(A1[I₁,J₁] * A2[I₂,J₂] * a[I₁,I₂] * w[I₁,I₂] for I₁ in 1:n1, I₂ in 1:n2) for J₁ in 1:n1′, J₂ in 1:n2′] ./ w′
+    return RationalBSplineManifold(a′, w′, Ps′)
+end
+function refinement(M::RationalBSplineManifold{3}, Ps′::Tuple{<:AbstractBSplineSpace,<:AbstractBSplineSpace,<:AbstractBSplineSpace})
+    Ps = bsplinespaces(M)
+    a = controlpoints(M)
+    w = weights(M)
+    (n1,n2,n3) = dim.(Ps)
+    (n1′,n2′,n3′) = dim.(Ps′)
+    (A1,A2,A3) = changebasis.(Ps, Ps′)
+
+    w′ = [sum(A1[I₁,J₁] * A2[I₂,J₂] * A3[I₃,J₃] * w[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′]
+    a′ = [sum(A1[I₁,J₁] * A2[I₂,J₂] * A3[I₃,J₃] * a[I₁,I₂,I₃] * w[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′] ./ w′
+    return RationalBSplineManifold(a′, w′, Ps′)
+end
+
 @doc raw"""
 Refinement of B-spline manifold with additional degree and knotvector.
 """
-function refinement(M::BSplineManifold; p₊::Union{Nothing,NTuple{Dim,Int}}=nothing, k₊::Union{Nothing,NTuple{Dim,KnotVector{T}}}=nothing) where {Dim, T}
-    Ps = collect(bsplinespaces(M))
+function refinement(M::AbstractManifold{Dim}; p₊::Union{Nothing,NTuple{Dim,Int}}=nothing, k₊::Union{Nothing,NTuple{Dim,KnotVector{T}}}=nothing) where {Dim, T}
+    Ps = bsplinespaces(M)
     d = length(Ps)
     if isnothing(p₊)
         p₊ = zeros(Int, d)


Test patch:
diff --git a/test/runtests.jl b/test/runtests.jl
index ac3c4a991..83e1bc05c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,6 +6,11 @@ using Random
 using StaticArrays
 using GeometryBasics
 
+function rand_interval(D::ClosedInterval)
+    a,b = endpoints(D)
+    return a+(b-a)*rand()
+end
+
 @testset "BasicBSpline.jl" begin
     include("test_KnotVector.jl")
     include("test_BSplineSpace.jl")
@@ -16,5 +21,7 @@ using GeometryBasics
     include("test_Derivative.jl")
     include("test_ChangeBasis.jl")
     include("test_BSplineManifold.jl")
+    include("test_RationalBSplineManifold.jl")
+    include("test_Refinement.jl")
     include("test_Fitting.jl")
 end
diff --git a/test/test_RationalBSplineManifold.jl b/test/test_RationalBSplineManifold.jl
new file mode 100644
index 000000000..46aaa650c
--- /dev/null
+++ b/test/test_RationalBSplineManifold.jl
@@ -0,0 +1,90 @@
+@testset "RationalBSplineManifold" begin
+    @testset "1dim-arc" begin
+        a = [SVector(1,0), SVector(1,1), SVector(0,1), SVector(-1,1), SVector(-1,0)]
+        w = [1, 1/√2, 1, 1/√2, 1]
+
+        k = KnotVector([0,0,0,1,1,2,2,2])
+        p = 2
+        P = BSplineSpace{p}(k)
+
+        M = BasicBSpline.RationalBSplineManifold(a,w,(P,))
+        for _ in 1:100
+            t1 = 2rand()
+            @test norm(M(t1)) ≈ 1
+        end
+        @test_throws DomainError M(-0.1)
+        @test_throws DomainError M(2.1)
+    end
+
+    @testset "compatibility with BSplineManifold (1dim)" begin
+        k1 = KnotVector(rand(30)) + 5KnotVector(0,1)
+        p1 = 4
+        P1 = BSplineSpace{p1}(k1)
+        n1 = dim(P1)
+        a = randn(n1)
+        w1 = ones(n1)
+        w2 = 2ones(n1)
+
+        M = BasicBSpline.BSplineManifold(a,(P1,))
+        M1 = BasicBSpline.RationalBSplineManifold(a,w1,(P1,))
+        M2 = BasicBSpline.RationalBSplineManifold(a,w2,(P1,))
+        for _ in 1:100
+            t1 = rand()
+            @test M(t1) ≈ M1(t1) atol=1e-14
+            @test M(t1) ≈ M2(t1) atol=1e-14
+        end
+    end
+
+    @testset "compatibility with BSplineManifold (2dim)" begin
+        k1 = KnotVector(rand(30)) + 5KnotVector(0,1)
+        k2 = KnotVector(rand(30)) + 5KnotVector(0,1)
+        p1 = 2
+        p2 = 3
+        P1 = BSplineSpace{p1}(k1)
+        P2 = BSplineSpace{p2}(k2)
+        n1 = dim(P1)
+        n2 = dim(P2)
+        a = randn(n1,n2)
+        w1 = ones(n1,n2)
+        w2 = 2ones(n1,n2)
+
+        M = BasicBSpline.BSplineManifold(a,(P1,P2))
+        M1 = BasicBSpline.RationalBSplineManifold(a,w1,(P1,P2))
+        M2 = BasicBSpline.RationalBSplineManifold(a,w2,(P1,P2))
+        for _ in 1:100
+            t1 = rand()
+            t2 = rand()
+            @test M(t1,t2) ≈ M1(t1,t2) atol=1e-14
+            @test M(t1,t2) ≈ M2(t1,t2) atol=1e-14
+        end
+    end
+
+    @testset "compatibility with BSplineManifold (3dim)" begin
+        k1 = KnotVector(rand(30)) + 5KnotVector(0,1)
+        k2 = KnotVector(rand(30)) + 5KnotVector(0,1)
+        k3 = KnotVector(rand(30)) + 5KnotVector(0,1)
+        p1 = 2
+        p2 = 3
+        p3 = 4
+        P1 = BSplineSpace{p1}(k1)
+        P2 = BSplineSpace{p2}(k2)
+        P3 = BSplineSpace{p3}(k3)
+        n1 = dim(P1)
+        n2 = dim(P2)
+        n3 = dim(P3)
+        a = randn(n1,n2,n3)
+        w1 = ones(n1,n2,n3)
+        w2 = 2ones(n1,n2,n3)
+
+        M = BasicBSpline.BSplineManifold(a,(P1,P2,P3))
+        M1 = BasicBSpline.RationalBSplineManifold(a,w1,(P1,P2,P3))
+        M2 = BasicBSpline.RationalBSplineManifold(a,w2,(P1,P2,P3))
+        for _ in 1:100
+            t1 = rand()
+            t2 = rand()
+            t3 = rand()
+            @test M(t1,t2,t3) ≈ M1(t1,t2,t3) atol=1e-14
+            @test M(t1,t2,t3) ≈ M2(t1,t2,t3) atol=1e-14
+        end
+    end
+end
diff --git a/test/test_Refinement.jl b/test/test_Refinement.jl
new file mode 100644
index 000000000..10caabc66
--- /dev/null
+++ b/test/test_Refinement.jl
@@ -0,0 +1,195 @@
+@testset "Refinement" begin
+    p1,p2,p3 = 3,2,4
+
+    k1 = KnotVector(1:8)
+    k2 = KnotVector(1,3) + p2*KnotVector(1,3)
+    k3 = KnotVector(0,2) + p3*KnotVector(1,2)
+
+    P1 = BSplineSpace{p1}(k1)
+    P2 = BSplineSpace{p2}(k2)
+    P3 = BSplineSpace{p3}(k3)
+
+    P1′ = BSplineSpace{3}(k1+KnotVector(1.2,5.5,7.2,8.5))
+    P2′ = expandspace(P2, p₊=1, k₊=KnotVector(1.8))
+    P3′ = BSplineSpace{4}(UniformKnotVector(-3:6))
+
+    n1,n2,n3 = dim(P1),dim(P2),dim(P3)
+    D1,D2,D3 = domain(P1),domain(P2),domain(P3)
+
+    @test P1 ⊆ P1′
+    @test P2 ⊆ P2′
+    @test !(P3 ⊆ P3′)
+
+    @test !(P1 ⊑ P1′)
+    @test P2 ⊑ P2′
+    @test P3 ⊑ P3′
+
+    @testset "1dim" begin
+        a = [Point(i, rand()) for i in 1:n1]
+        p₊ = (1,)
+        k₊ = (KnotVector(1.8,1.9),)
+
+        @testset "BSplineManifold" begin
+            w = ones(n1)
+            M = BSplineManifold(a, (P1,))
+            M0 = refinement(M)
+            M1 = refinement(M, k₊=k₊)
+            # M2 = refinement(M, p₊=p₊)
+            # M3 = refinement(M, p₊=p₊, k₊=k₊)
+            M4 = refinement(M, (P1′,))
+            R = RationalBSplineManifold(a, w, (P1,))
+            R0 = refinement(R)
+            R1 = refinement(R, k₊=k₊)
+            # R2 = refinement(R, p₊=p₊)
+            # R3 = refinement(R, p₊=p₊, k₊=k₊)
+            R4 = refinement(R, (P1′,))
+
+            for _ in 1:100
+                t1 = rand_interval(D1)
+                @test M(t1) ≈ R(t1)
+                @test M(t1) ≈ M1(t1)
+                # @test M(t1) ≈ M2(t1)
+                # @test M(t1) ≈ M3(t1)
+                @test M(t1) ≈ M4(t1)
+                @test R(t1) ≈ R1(t1)
+                # @test R(t1) ≈ R2(t1)
+                # @test R(t1) ≈ R3(t1)
+                @test R(t1) ≈ R4(t1)
+            end
+        end
+
+        @testset "RationalBSplineManifold" begin
+            w = rand(n1) .+ 1
+            R = RationalBSplineManifold(a, w, (P1,))
+            R0 = refinement(R)
+            R1 = refinement(R, k₊=k₊)
+            # R2 = refinement(R, p₊=p₊)
+            # R3 = refinement(R, p₊=p₊, k₊=k₊)
+            R4 = refinement(R, (P1′,))
+
+            for _ in 1:100
+                t1 = rand_interval(D1)
+                @test R(t1) ≈ R1(t1)
+                # @test R(t1) ≈ R2(t1)
+                # @test R(t1) ≈ R3(t1)
+                @test R(t1) ≈ R4(t1)
+            end
+        end
+    end
+
+    @testset "2dim" begin
+        a = [Point(i, rand()) for i in 1:n1, j in 1:n2]
+        p₊ = (1,2)
+        k₊ = (KnotVector(1.8,1.9),KnotVector())
+
+        @testset "BSplineManifold" begin
+            w = ones(n1,n2)
+
+            M = BSplineManifold(a, (P1,P2))
+            M0 = refinement(M)
+            M1 = refinement(M, k₊=k₊)
+            # M2 = refinement(M, p₊=p₊)
+            # M3 = refinement(M, p₊=p₊, k₊=k₊)
+            M4 = refinement(M, (P1′,P2′))
+            R = RationalBSplineManifold(a, w, (P1,P2))
+            R0 = refinement(R)
+            R1 = refinement(R, k₊=k₊)
+            # R2 = refinement(R, p₊=p₊)
+            # R3 = refinement(R, p₊=p₊, k₊=k₊)
+            R4 = refinement(R, (P1′,P2′))
+
+            for _ in 1:100
+                t1 = rand_interval(D1)
+                t2 = rand_interval(D2)
+                @test M(t1,t2) ≈ R(t1,t2)
+                @test M(t1,t2) ≈ M1(t1,t2)
+                # @test M(t1,t2) ≈ M2(t1,t2)
+                # @test M(t1,t2) ≈ M3(t1,t2)
+                @test M(t1,t2) ≈ M4(t1,t2)
+                @test R(t1,t2) ≈ R1(t1,t2)
+                # @test R(t1,t2) ≈ R2(t1,t2)
+                # @test R(t1,t2) ≈ R3(t1,t2)
+                @test R(t1,t2) ≈ R4(t1,t2)
+            end
+        end
+
+        @testset "RationalBSplineManifold" begin
+            w = rand(n1,n2) .+ 1
+
+            R = RationalBSplineManifold(a, w, (P1,P2))
+            R0 = refinement(R)
+            R1 = refinement(R, k₊=k₊)
+            # R2 = refinement(R, p₊=p₊)
+            # R3 = refinement(R, p₊=p₊, k₊=k₊)
+            R4 = refinement(R, (P1′,P2′))
+
+            for _ in 1:100
+                t1 = rand_interval(D1)
+                t2 = rand_interval(D2)
+                @test R(t1,t2) ≈ R1(t1,t2)
+                # @test R(t1,t2) ≈ R2(t1,t2)
+                # @test R(t1,t2) ≈ R3(t1,t2)
+                @test R(t1,t2) ≈ R4(t1,t2)
+            end
+        end
+    end
+
+    @testset "3dim" begin
+        a = [Point(i, rand()) for i in 1:n1, j in 1:n2, k in 1:n3]
+        p₊ = (1,2,0)
+        k₊ = (KnotVector(1.8,1.9),KnotVector(),KnotVector(42))
+
+        @testset "BSplineManifold" begin
+            w = ones(n1,n2,n3)
+
+            M = BSplineManifold(a, (P1,P2,P3))
+            M0 = refinement(M)
+            M1 = refinement(M, k₊=k₊)
+            # M2 = refinement(M, p₊=p₊)
+            # M3 = refinement(M, p₊=p₊, k₊=k₊)
+            M4 = refinement(M, (P1′,P2′,P3′))
+            R = RationalBSplineManifold(a, w, (P1,P2,P3))
+            R0 = refinement(R)
+            R1 = refinement(R, k₊=k₊)
+            # R2 = refinement(R, p₊=p₊)
+            # R3 = refinement(R, p₊=p₊, k₊=k₊)
+            R4 = refinement(R, (P1′,P2′,P3′))
+
+            for _ in 1:100
+                t1 = rand_interval(D1)
+                t2 = rand_interval(D2)
+                t3 = rand_interval(D3)
+                @test M(t1,t2,t3) ≈ R(t1,t2,t3)
+                @test M(t1,t2,t3) ≈ M1(t1,t2,t3)
+                # @test M(t1,t2,t3) ≈ M2(t1,t2,t3)
+                # @test M(t1,t2,t3) ≈ M3(t1,t2,t3)
+                @test M(t1,t2,t3) ≈ M4(t1,t2,t3)
+                @test R(t1,t2,t3) ≈ R1(t1,t2,t3)
+                # @test R(t1,t2,t3) ≈ R2(t1,t2,t3)
+                # @test R(t1,t2,t3) ≈ R3(t1,t2,t3)
+                @test R(t1,t2,t3) ≈ R4(t1,t2,t3)
+            end
+        end
+
+        @testset "RationalBSplineManifold" begin
+            w = rand(n1,n2,n3) .+ 1
+
+            R = RationalBSplineManifold(a, w, (P1,P2,P3))
+            R0 = refinement(R)
+            R1 = refinement(R, k₊=k₊)
+            # R2 = refinement(R, p₊=p₊)
+            # R3 = refinement(R, p₊=p₊, k₊=k₊)
+            R4 = refinement(R, (P1′,P2′,P3′))
+
+            for _ in 1:100
+                t1 = rand_interval(D1)
+                t2 = rand_interval(D2)
+                t3 = rand_interval(D3)
+                @test R(t1,t2,t3) ≈ R1(t1,t2,t3)
+                # @test R(t1,t2,t3) ≈ R2(t1,t2,t3)
+                # @test R(t1,t2,t3) ≈ R3(t1,t2,t3)
+                @test R(t1,t2,t3) ≈ R4(t1,t2,t3)
+            end
+        end
+    end
+end