Repository: hyrodium/BasicBSpline.jl
Pull request #20 (2020-09-05 05:12:35)

Feature/specific bspline manifold
* Add `BSplineCurve`, `BSplineSurface` and `BSplineSolid`
* refactoring


Patch:
diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
index fd5df3e4..111704ad 100644
--- a/.JuliaFormatter.toml
+++ b/.JuliaFormatter.toml
@@ -1,3 +1,3 @@
 indent = 4
-margin = 100
+margin = 150
 always_for_in = true
diff --git a/src/BasicBSpline.jl b/src/BasicBSpline.jl
index 34306c0d..16a7c39a 100644
--- a/src/BasicBSpline.jl
+++ b/src/BasicBSpline.jl
@@ -15,6 +15,7 @@ export isproper, properdim
 export degree, knots
 export 𝔫
 export FastBSplineSpace, FastBSplineManifold
+export BSplineCurve, BSplineSurface, BSplineSolid
 export AbstractBSplineManifold, AbstractBSplineSpace
 export fittingcontrolpoints
 
diff --git a/src/_BSplineBasis.jl b/src/_BSplineBasis.jl
index aebc9946..774fc937 100644
--- a/src/_BSplineBasis.jl
+++ b/src/_BSplineBasis.jl
@@ -192,6 +192,8 @@ Modified version.
 =p\left(\frac{1}{k_{i+p}-k_{i}}B_{(i,p-1,k)}(t)-\frac{1}{k_{i+p+1}-k_{i+1}}B_{(i+1,p-1,k)}(t)\right)
 ```
 """
+bsplinebasis′
+
 function bsplinebasis′(P::BSplineSpace, t)::Array{Float64,1}
     p = P.degree
     k = P.knots
diff --git a/src/_BSplineManifold.jl b/src/_BSplineManifold.jl
index d8872792..d5326ff4 100644
--- a/src/_BSplineManifold.jl
+++ b/src/_BSplineManifold.jl
@@ -36,12 +36,6 @@ struct BSplineManifold <: AbstractBSplineManifold
     end
 end
 
-function BSplineManifold(Ps::AbstractArray{<:AbstractBSplineSpace,1}, a::Array{Array{T,1}} where T<:Real)
-    d̂ = length(a[1])
-    A = reshape(transpose(hcat(reshape(a,prod(size(a)))...)), size(a)..., d̂)
-    return BSplineManifold(Ps, A)
-end
-
 """
 convert AbstractBSplineManifold to BSplineManifold
 """
@@ -98,6 +92,14 @@ end
 @doc raw"""
 Calculate the dimension of B-spline manifold.
 """
-function dim(M::AbstractBSplineManifold)
-    length(M.bsplinespaces)
+dim
+
+dim(M::AbstractBSplineManifold) = length(M.bsplinespaces)
+
+function bsplinespaces(M::BSplineManifold)
+    return M.bsplinespaces
+end
+
+function controlpoints(M::BSplineManifold)
+    return M.controlpoints
 end
diff --git a/src/_BSplineSpace.jl b/src/_BSplineSpace.jl
index fb895aa9..4ab30887 100644
--- a/src/_BSplineSpace.jl
+++ b/src/_BSplineSpace.jl
@@ -80,7 +80,7 @@ Check inclusive relationship between B-spline spaces.
 \subseteq\mathcal{P}[p',k']|_{[\sharp k'_{p'+1},k'_{\sharp k'-p'}]}
 ```
 """
-function ⊑(P::AbstractBSplineSpace, P′::AbstractBSplineSpace)
+function issqsubset(P::AbstractBSplineSpace, P′::AbstractBSplineSpace)
     p = degree(P)
     k = knots(P)
     p′ = degree(P′)
@@ -99,18 +99,10 @@ function ⊑(P::AbstractBSplineSpace, P′::AbstractBSplineSpace)
     return inner_knots+p₊*unique(inner_knots) ⊆ inner_knots′
 end
 
-function ⊒(P::AbstractBSplineSpace, P′::AbstractBSplineSpace)
-    return ⊑(P′, P)
-end
-
-function ⋢(P::AbstractBSplineSpace, P′::AbstractBSplineSpace)
-    return !⊑(P, P′)
-end
-
-function ⋣(P::AbstractBSplineSpace, P′::AbstractBSplineSpace)
-    return !⊒(P, P′)
-end
-
+const ⊑ = issqsubset
+⊒(l, r) = r ⊑ l
+⋢(l, r) = !⊑(l, r)
+⋣(l, r) = r ⋢ l
 
 function iszeros(P::AbstractBSplineSpace)
     p = degree(P)
diff --git a/src/_FastBSplineBasis.jl b/src/_FastBSplineBasis.jl
index 21c6676c..82f4ab84 100644
--- a/src/_FastBSplineBasis.jl
+++ b/src/_FastBSplineBasis.jl
@@ -1,100 +1,81 @@
-# TODO: use macro for bsplienbasis(::FastBSplineSpace)
-function bsplinebasis(i::Integer,P::FastBSplineSpace{0},t::Real)
-    k_1, k_2, k_end = P.vector[i], P.vector[i+1], P.vector[end]
-    B_1 = Float64((k_1 ≤ t < k_2) | (k_1 ≤ t ≤ k_2 == k_end))
-
-    return B_1
+function _s(c::Char, i::Int)
+    string(c)*"_"*string(i)
 end
 
-function bsplinebasis(i::Integer,P::FastBSplineSpace{1},t::Real)
-    ÷(a,b) = ifelse(b == 0.0, 0.0, a/b)
-    k_1, k_2, k_3, k_end = P.vector[i], P.vector[i+1], P.vector[i+2], P.vector[end]
-    B_1, B_2 = Float64(k_1 ≤ t < k_2), Float64((k_2 ≤ t < k_3) | (k_2 ≤ t ≤ k_3 == k_end))
-
-    K_1, K_2 = (t-k_1)÷(k_2-k_1), (t-k_2)÷(k_3-k_2)
-    B_1 = K_1*B_1+(1-K_2)*B_2
-
-    return B_1
+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 bsplinebasis(i::Integer,P::FastBSplineSpace{2},t::Real)
-    ÷(a,b) = ifelse(b == 0.0, 0.0, a/b)
-    k_1, k_2, k_3, k_4, k_end = P.vector[i], P.vector[i+1], P.vector[i+2], P.vector[i+3], P.vector[end]
-    B_1, B_2, B_3 = Float64(k_1 ≤ t < k_2), Float64(k_2 ≤ t < k_3), Float64((k_3 ≤ t < k_4) | (k_3 ≤ t ≤ k_4 == k_end))
-
-    K_1, K_2, K_3 = (t-k_1)÷(k_2-k_1), (t-k_2)÷(k_3-k_2), (t-k_3)÷(k_4-k_3)
-    B_1, B_2 = K_1*B_1+(1-K_2)*B_2, K_2*B_2+(1-K_3)*B_3
-
-    K_1, K_2 = (t-k_1)÷(k_3-k_1), (t-k_2)÷(k_4-k_2)
-    B_1 = K_1*B_1+(1-K_2)*B_2
-
-    return B_1
+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 bsplinebasis(i::Integer,P::FastBSplineSpace{3},t::Real)
-    ÷(a,b) = ifelse(b == 0.0, 0.0, a/b)
-    k_1, k_2, k_3, k_4, k_5, k_end = P.vector[i], P.vector[i+1], P.vector[i+2], P.vector[i+3], P.vector[i+4], P.vector[end]
-    B_1, B_2, B_3, B_4 = Float64(k_1 ≤ t < k_2), Float64(k_2 ≤ t < k_3), Float64(k_3 ≤ t < k_4), Float64((k_4 ≤ t < k_5) | (k_4 ≤ t ≤ k_5 == k_end))
-
-    K_1, K_2, K_3, K_4 = (t-k_1)÷(k_2-k_1), (t-k_2)÷(k_3-k_2), (t-k_3)÷(k_4-k_3), (t-k_4)÷(k_5-k_4)
-    B_1, B_2, B_3 = K_1*B_1+(1-K_2)*B_2, K_2*B_2+(1-K_3)*B_3, K_3*B_3+(1-K_4)*B_4
-
-    K_1, K_2, K_3 = (t-k_1)÷(k_3-k_1), (t-k_2)÷(k_4-k_2), (t-k_3)÷(k_5-k_3)
-    B_1, B_2 = K_1*B_1+(1-K_2)*B_2, K_2*B_2+(1-K_3)*B_3
-
-    K_1, K_2 = (t-k_1)÷(k_4-k_1), (t-k_2)÷(k_5-k_2)
-    B_1 = K_1*B_1+(1-K_2)*B_2
-
-    return B_1
+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
 
-function bsplinebasis′₊₀(i::Integer,P::FastBSplineSpace{0},t::Real)
-    return 0.0
+for p in 0:MAX_DEGREE
+    @eval function bsplinebasis₊₀(i::Integer, P::FastBSplineSpace{$p}, 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.knotvector[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₋₀(i::Integer, P::FastBSplineSpace{$p}, 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.knotvector[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(i::Integer, P::FastBSplineSpace{$p}, t::Real)
+        ÷(a,b) = ifelse(b == 0.0, 0.0, a/b)
+        k_end = P.knotvector[end]
+        $(Meta.parse(join([_s('k',j) for j in 1:p+2],',')*" = "*join(["P.knotvector[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
 
-function bsplinebasis′₊₀(i::Integer,P::FastBSplineSpace{1},t::Real)
-    ÷(a,b) = ifelse(b == 0.0, 0.0, a/b)
-    k_1, k_2, k_3, k_end = P.vector[i], P.vector[i+1], P.vector[i+2], P.vector[end]
-    B_1, B_2 = Float64(k_1 ≤ t < k_2), Float64(k_2 ≤ t < k_3)
-
-    K_1, K_2 = 1÷(k_2-k_1), 1÷(k_3-k_2)
-    B_1 = K_1*B_1-K_2*B_2
-
-    return B_1
-end
-
-function bsplinebasis′₊₀(i::Integer,P::FastBSplineSpace{2},t::Real)
-    ÷(a,b) = ifelse(b == 0.0, 0.0, a/b)
-    k_1, k_2, k_3, k_4, k_end = P.vector[i], P.vector[i+1], P.vector[i+2], P.vector[i+3], P.vector[end]
-    B_1, B_2, B_3 = Float64(k_1 ≤ t < k_2), Float64(k_2 ≤ t < k_3), Float64(k_3 ≤ t < k_4)
-
-    K_1, K_2, K_3 = (t-k_1)÷(k_2-k_1), (t-k_2)÷(k_3-k_2), (t-k_3)÷(k_4-k_3)
-    B_1, B_2 = K_1*B_1+(1-K_2)*B_2, K_2*B_2+(1-K_3)*B_3
-
-    K_1, K_2 = 1÷(k_3-k_1), 1÷(k_4-k_2)
-    B_1 = 2*(K_1*B_1-K_2*B_2)
-
-    return B_1
+bsplinebasis′₊₀(i::Integer,P::FastBSplineSpace{0},t::Real) = 0.0
+bsplinebasis′₋₀(i::Integer,P::FastBSplineSpace{0},t::Real) = 0.0
+bsplinebasis′(i::Integer,P::FastBSplineSpace{0},t::Real) = 0.0
+
+for p in 1:MAX_DEGREE
+    @eval function bsplinebasis′₊₀(i::Integer, P::FastBSplineSpace{$p}, 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.knotvector[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′₋₀(i::Integer, P::FastBSplineSpace{$p}, 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.knotvector[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′(i::Integer, P::FastBSplineSpace{$p}, t::Real)
+    #     ÷(a,b) = ifelse(b == 0.0, 0.0, a/b)
+    #     k_end = P.knotvector[end]
+    #     $(Meta.parse(join([_s('k',j) for j in 1:p+2],',')*" = "*join(["P.knotvector[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′₊₀(i::Integer,P::FastBSplineSpace{3},t::Real)
-    ÷(a,b) = ifelse(b == 0.0, 0.0, a/b)
-    k_1, k_2, k_3, k_4, k_5, k_end = P.vector[i], P.vector[i+1], P.vector[i+2], P.vector[i+3], P.vector[i+4], P.vector[end]
-    B_1, B_2, B_3, B_4 = Float64(k_1 ≤ t < k_2), Float64(k_2 ≤ t < k_3), Float64(k_3 ≤ t < k_4), Float64(k_4 ≤ t < k_5)
-
-    K_1, K_2, K_3, K_4 = (t-k_1)÷(k_2-k_1), (t-k_2)÷(k_3-k_2), (t-k_3)÷(k_4-k_3), (t-k_4)÷(k_5-k_4)
-    B_1, B_2, B_3 = K_1*B_1+(1-K_2)*B_2, K_2*B_2+(1-K_3)*B_3, K_3*B_3+(1-K_4)*B_4
-
-    K_1, K_2, K_3 = (t-k_1)÷(k_3-k_1), (t-k_2)÷(k_4-k_2), (t-k_3)÷(k_5-k_3)
-    B_1, B_2 = K_1*B_1+(1-K_2)*B_2, K_2*B_2+(1-K_3)*B_3
-
-    K_1, K_2 = 1÷(k_4-k_1), 1÷(k_5-k_2)
-    B_1 = 3*(K_1*B_1-K_2*B_2)
-
-    return B_1
-end
-
-
 """
 TODO: faster.....
 """
diff --git a/src/_FastBSplineManifold.jl b/src/_FastBSplineManifold.jl
index 92c82539..09d1f765 100644
--- a/src/_FastBSplineManifold.jl
+++ b/src/_FastBSplineManifold.jl
@@ -5,7 +5,6 @@ B-spline manifold for lower polynomial degree
 TODO: make the field `bsplinespaces` to be conposite type, not abstract type, for performance
 """
 struct FastBSplineManifold <: AbstractBSplineManifold
-    # TODO: avoid abstract type FastBSplineSpace
     bsplinespaces::Array{T,1} where T <: FastBSplineSpace
     controlpoints::Array{Float64}
     function FastBSplineManifold(Ps::AbstractArray{<:AbstractBSplineSpace,1}, a::AbstractArray{<:Real})
@@ -20,17 +19,109 @@ struct FastBSplineManifold <: AbstractBSplineManifold
     end
 end
 
-function FastBSplineManifold(Ps::AbstractArray{<:AbstractBSplineSpace,1}, a::Array{Array{T,1}} where T<:Real)
-    d̂ = length(a[1])
-    A = reshape(transpose(hcat(reshape(a,prod(size(a)))...)), size(a)..., d̂)
-    return FastBSplineManifold(Ps, A)
+struct BSplineCurve{p1} <: AbstractBSplineManifold
+    knotvector1::Array{Float64,1}
+    controlpoints::Array{Float64,2}
+    function BSplineCurve(Ps::AbstractArray{<:AbstractBSplineSpace,1}, a::AbstractArray{<:Real})
+        p1, k1 = degree(Ps[1]), knots(Ps[1]).vector
+        if collect(size(a)[1:end-1]) ≠ dim.(Ps)
+            error("dimension does not match")
+        elseif p1 > MAX_DEGREE
+            return BSplineManifold(Ps, a)
+        else
+            a = convert(Array{Float64}, a)
+            new{p1}(k1,a)
+        end
+    end
+    function BSplineCurve{q1}(Ps::AbstractArray{<:AbstractBSplineSpace,1}, a::AbstractArray{<:Real}) where q1
+        return BSplineCurve(Ps,a)
+    end
+end
+
+struct BSplineSurface{p1,p2} <: AbstractBSplineManifold
+    knotvector1::Array{Float64,1}
+    knotvector2::Array{Float64,1}
+    controlpoints::Array{Float64,3}
+    function BSplineSurface(Ps::AbstractArray{<:AbstractBSplineSpace,1}, a::AbstractArray{<:Real})
+        p1, k1 = degree(Ps[1]), knots(Ps[1]).vector
+        p2, k2 = degree(Ps[2]), knots(Ps[2]).vector
+        if collect(size(a)[1:end-1]) ≠ dim.(Ps)
+            error("dimension does not match")
+        elseif p1 > MAX_DEGREE || p2 > MAX_DEGREE
+            return BSplineManifold(Ps, a)
+        else
+            a = convert(Array{Float64}, a)
+            new{p1,p2}(k1,k2,a)
+        end
+    end
+    function BSplineSurface{q1,q2}(Ps::AbstractArray{<:AbstractBSplineSpace,1}, a::AbstractArray{<:Real}) where q1 where q2
+        return BSplineSurface(Ps,a)
+    end
+end
+
+struct BSplineSolid{p1,p2,p3} <: AbstractBSplineManifold
+    knotvector1::Array{Float64,1}
+    knotvector2::Array{Float64,1}
+    knotvector3::Array{Float64,1}
+    controlpoints::Array{Float64,4}
+    function BSplineSolid(Ps::AbstractArray{<:AbstractBSplineSpace,1}, a::AbstractArray{<:Real})
+        p1, k1 = degree(Ps[1]), knots(Ps[1]).vector
+        p2, k2 = degree(Ps[2]), knots(Ps[2]).vector
+        p3, k3 = degree(Ps[3]), knots(Ps[3]).vector
+        if collect(size(a)[1:end-1]) ≠ dim.(Ps)
+            error("dimension does not match")
+        elseif p1 > MAX_DEGREE || p2 > MAX_DEGREE || p3 > MAX_DEGREE
+            return BSplineManifold(Ps, a)
+        else
+            a = convert(Array{Float64}, a)
+            new{p1,p2,p3}(k1,k2,k3,a)
+        end
+    end
+    function BSplineSolid{q1,q2,q3}(Ps::AbstractArray{<:AbstractBSplineSpace,1}, a::AbstractArray{<:Real}) where q1 where q2 where q3
+        return BSplineSolid(Ps,a)
+    end
+end
+
+for fname in (:BSplineCurve, :BSplineSurface, :BSplineSolid, :FastBSplineManifold, :BSplineManifold)
+    @eval function $fname(Ps::AbstractArray{<:AbstractBSplineSpace,1}, a::Array{Array{T,1}} where T<:Real)
+        d̂ = length(a[1])
+        A = reshape(transpose(hcat(reshape(a,prod(size(a)))...)), size(a)..., d̂)
+        return $fname(Ps, A)
+    end
+end
+
+function BSplineManifold(M::BSplineCurve{p1}) where p1
+    k1 = Knots(M.knotvector1)
+    P1 = BSplineSpace(p1, k1)
+    a = M.controlpoints
+    return BSplineManifold([P1], a)
+end
+
+function BSplineManifold(M::BSplineSurface{p1, p2}) where p1 where p2
+    k1 = Knots(M.knotvector1)
+    P1 = BSplineSpace(p1, k1)
+    k2 = Knots(M.knotvector2)
+    P2 = BSplineSpace(p2, k2)
+    a = M.controlpoints
+    return BSplineManifold([P1, P2], a)
+end
+
+function BSplineManifold(M::BSplineSolid{p1, p2, p3}) where p1 where p2 where p3
+    k1 = Knots(M.knotvector1)
+    P1 = BSplineSpace(p1, k1)
+    k2 = Knots(M.knotvector2)
+    P2 = BSplineSpace(p2, k2)
+    k3 = Knots(M.knotvector3)
+    P3 = BSplineSpace(p3, k3)
+    a = M.controlpoints
+    return BSplineManifold([P1, P2, P3], a)
 end
 
 """
 convert AbstractBSplineManifold to FastBSplineManifold
 """
 function FastBSplineManifold(M::AbstractBSplineManifold)
-    FastBSplineManifold(BSplineSpace.(M.bsplinespaces), M.controlpoints)
+    return FastBSplineManifold(BSplineSpace.(M.bsplinespaces), M.controlpoints)
 end
 
 @doc raw"""
@@ -45,9 +136,54 @@ function mapping(M::FastBSplineManifold, t::Array{<:Real,1})
     return mapping(BSplineManifold(M),t)
 end
 
-@doc raw"""
-Calculate the dimension of B-spline manifold.
-"""
-function dim(M::FastBSplineManifold)
-    length(M.bsplinespaces)
+function mapping(M::BSplineCurve, t::Array{<:Real,1})
+    # TODO: faster
+    return mapping(BSplineManifold(M),t)
+end
+
+function mapping(M::BSplineSurface, t::Array{<:Real,1})
+    # TODO: faster
+    return mapping(BSplineManifold(M),t)
+end
+
+function mapping(M::BSplineSolid, t::Array{<:Real,1})
+    # TODO: faster
+    return mapping(BSplineManifold(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 [FastBSplineSpace(p1,Knots(M.knotvector1))]
+end
+
+function bsplinespaces(M::BSplineSurface{p1,p2}) where p1 where p2
+    return [FastBSplineSpace(p1,Knots(M.knotvector1)), FastBSplineSpace(p2,Knots(M.knotvector2))]
+end
+
+function bsplinespaces(M::BSplineSolid{p1,p2,p3}) where p1 where p2 where p3
+    return [FastBSplineSpace(p1,Knots(M.knotvector1)), FastBSplineSpace(p2,Knots(M.knotvector2)), FastBSplineSpace(p3,Knots(M.knotvector3))]
+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
index a44a2f24..287d9b2e 100644
--- a/src/_FastBSplineSpace.jl
+++ b/src/_FastBSplineSpace.jl
@@ -4,7 +4,7 @@
 B-spline space for lower polynomial degree
 """
 struct FastBSplineSpace{p} <: AbstractBSplineSpace
-    vector::Array{Float64,1}
+    knotvector::Array{Float64,1}
     function FastBSplineSpace(p::Integer, knots::Knots)
         if p < 0
             error("degree of polynominal must be non-negative")
@@ -35,7 +35,7 @@ function degree(P::FastBSplineSpace{p}) where p
 end
 
 function knots(P::FastBSplineSpace)
-    return Knots(P.vector)
+    return Knots(P.knotvector)
 end
 
 @doc raw"""
diff --git a/src/_Fitting.jl b/src/_Fitting.jl
index 14ae44fa..6683a59b 100644
--- a/src/_Fitting.jl
+++ b/src/_Fitting.jl
@@ -89,6 +89,7 @@ function fittingcontrolpoints_2dim(func::Function, Ps::Array{T,1} where T<:Abstr
         else
             S = 0.0
             for i1 in rng[1], i2 in rng[2]
+                # TODO: faster
                 S += GaussianQuadrature_2dim(t->(bsplinebasis(I[1],Ps[1],t[1]) * bsplinebasis(I[2],Ps[2],t[2]) * bsplinebasis(J[1],Ps[1],t[1]) * bsplinebasis(J[2],Ps[2],t[2])), [k[1][i1]..k[1][i1+1], k[2][i2]..k[2][i2+1]], nodes, weights)
             end
             return S
diff --git a/src/_Knots.jl b/src/_Knots.jl
index d126da55..cb6fd477 100644
--- a/src/_Knots.jl
+++ b/src/_Knots.jl
@@ -43,6 +43,8 @@ Base.length(k::Knots) = length(k.vector)
 Base.firstindex(k) = 1
 Base.lastindex(k) = length(k)
 Base.unique(k::Knots) = Knots(unique(k.vector))
+Base.iterate(k::Knots) = iterate(k.vector)
+Base.iterate(k::Knots, i::Integer) = iterate(k.vector, i)
 
 function Base.:⊆(k::Knots, k′::Knots)
     K′ = copy(k′.vector)
diff --git a/src/_Refinement.jl b/src/_Refinement.jl
index c5af8024..5910ae20 100644
--- a/src/_Refinement.jl
+++ b/src/_Refinement.jl
@@ -1,5 +1,5 @@
 # Refinement
-function changebasis(P::BSplineSpace, P′::BSplineSpace)::Array{Float64,2}
+function changebasis(P::AbstractBSplineSpace, P′::AbstractBSplineSpace)::Array{Float64,2}
     p = degree(P)
     k = knots(P)
     p′ = degree(P′)
@@ -34,37 +34,39 @@ function changebasis(P::BSplineSpace, P′::BSplineSpace)::Array{Float64,2}
     else
         Q = [1:W[1]-1,[W[i]:W[i+1]-1 for i ∈ 1:length(W)-1]...,W[end]:n′]
     end
-    l = length(Q)
-    L = length.(Q)
+    λ = length(Q)
+    Λ = length.(Q)
     Ãᵖ = [Aᵖ[:,q] for q ∈ Q]
 
-    for ȷ ∈ 2:l-1
-        if L[ȷ] == 1
+    for ȷ ∈ 2:λ-1
+        if Λ[ȷ] == 1
             Ãᵖ[ȷ] .= NaN
         end
     end
-    for ȷ ∈ 1:l-1
-        if L[ȷ] ≥ 2
+    for ȷ ∈ 1:λ-1
+        if Λ[ȷ] ≥ 2
             t = k′[W[ȷ]]
-            # TODO: define for FastBSplineSpace
-            Ãᵖ[ȷ][:,end] = bsplinebasis₋₀(BSplineSpace(p,k),t)
+            for i in 1:n
+                Ãᵖ[ȷ][i,end] = bsplinebasis₋₀(i,P,t)
+            end
         end
     end
-    for ȷ ∈ 2:l
-        if L[ȷ] ≥ 2
+    for ȷ ∈ 2:λ
+        if Λ[ȷ] ≥ 2
             t = k′[W[ȷ-1]+p]
-            # TODO: define for FastBSplineSpace
-            Ãᵖ[ȷ][:,1] = bsplinebasis₊₀(BSplineSpace(p,k),t)
+            for i in 1:n
+                Ãᵖ[ȷ][i,1] = bsplinebasis₊₀(i,P,t)
+            end
         end
     end
-    for ȷ ∈ 1:l
-        if L[ȷ] ≥ 3
+    for ȷ ∈ 1:λ
+        if Λ[ȷ] ≥ 3
             r = Q[ȷ]
             A₊ = copy(Ãᵖ[ȷ])
             A₋ = copy(Ãᵖ[ȷ])
-            for j ∈ 1:L[ȷ]-2
+            for j ∈ 1:Λ[ȷ]-2
                 A₊[:,j+1] = A₊[:,j]+Δ[:,j+r[1]]
-                A₋[:,L[ȷ]-j] = A₋[:,L[ȷ]-j+1]-Δ[:,L[ȷ]-j+r[1]]
+                A₋[:,Λ[ȷ]-j] = A₋[:,Λ[ȷ]-j+1]-Δ[:,Λ[ȷ]-j+r[1]]
             end
             Ãᵖ[ȷ] = (A₊+A₋)/2
         end
@@ -73,17 +75,12 @@ function changebasis(P::BSplineSpace, P′::BSplineSpace)::Array{Float64,2}
     return Aᵖ .* Float64[bsplinesupport(j,P′) ⊆ bsplinesupport(i,P) for i ∈ 1:n, j ∈ 1:n′]
 end
 
-function changebasis(P::AbstractBSplineSpace, P′::AbstractBSplineSpace)
-    return changebasis(BSplineSpace(P), BSplineSpace(P′))
-end
-
-
 @doc raw"""
 Refinement of B-spline manifold with given B-spline spaces.
 """
 function refinement(M::AbstractBSplineManifold, Ps′::Array{T,1} where T <: AbstractBSplineSpace)
-    Ps = M.bsplinespaces
-    𝒂 = M.controlpoints
+    Ps = bsplinespaces(M)
+    𝒂 = controlpoints(M)
     d̂ = size(𝒂)[end]
     d = length(Ps)
     n = dim.(Ps)
@@ -106,8 +103,8 @@ end
 Refinement of B-spline manifold with additional degree and knots.
 """
 function refinement(M::AbstractBSplineManifold; p₊::Union{Nothing,AbstractArray{<:Integer,1}}=nothing, k₊::Union{Nothing,Array{Knots,1}}=nothing)
-    Ps = M.bsplinespaces
-    𝒂 = M.controlpoints
+    Ps = bsplinespaces(M)
+    𝒂 = controlpoints(M)
     d = length(Ps)
     d̂ = size(𝒂)[end]
     n = dim.(Ps)


Test patch:
diff --git a/test/runtests.jl b/test/runtests.jl
index 6ffa694d..d21f0bc7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -97,7 +97,92 @@ using Test
         @test P2 ⊈ P3
     end
 
-    @testset "BSplineManifold" begin
+    @testset "0th degree basis" begin
+        p = 0
+        k = Knots(rand(10))+Knots(0,1)
+        P = BSplineSpace(p,k)
+        fP = FastBSplineSpace(p,k)
+        n = dim(P)
+
+        @test prod([bsplinebasis₊₀(i,P,t) ≈ bsplinebasis₊₀(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis₊₀(i,P,t) ≈ bsplinebasis₊₀(i,fP,t) for i in 1:n, t in k])
+        @test prod([bsplinebasis₋₀(i,P,t) ≈ bsplinebasis₋₀(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis₋₀(i,P,t) ≈ bsplinebasis₋₀(i,fP,t) for i in 1:n, t in k])
+        @test prod([bsplinebasis(i,P,t) ≈ bsplinebasis(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis(i,P,t) ≈ bsplinebasis(i,fP,t) for i in 1:n, t in k])
+
+        @test prod([bsplinebasis′₊₀(i,P,t) ≈ bsplinebasis′₊₀(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis′₊₀(i,P,t) ≈ bsplinebasis′₊₀(i,fP,t) for i in 1:n, t in k])
+        @test prod([bsplinebasis′₋₀(i,P,t) ≈ bsplinebasis′₋₀(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis′₋₀(i,P,t) ≈ bsplinebasis′₋₀(i,fP,t) for i in 1:n, t in k])
+        @test prod([bsplinebasis′(i,P,t) ≈ bsplinebasis′(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis′(i,P,t) ≈ bsplinebasis′(i,fP,t) for i in 1:n, t in k])
+    end
+
+    @testset "1st degree basis" begin
+        p = 1
+        k = Knots(rand(10))+2*Knots(0,1)
+        P = BSplineSpace(p,k)
+        fP = FastBSplineSpace(p,k)
+        n = dim(P)
+        @test prod([bsplinebasis₊₀(i,P,t) ≈ bsplinebasis₊₀(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis₊₀(i,P,t) ≈ bsplinebasis₊₀(i,fP,t) for i in 1:n, t in k])
+        @test prod([bsplinebasis₋₀(i,P,t) ≈ bsplinebasis₋₀(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis₋₀(i,P,t) ≈ bsplinebasis₋₀(i,fP,t) for i in 1:n, t in k])
+        @test prod([bsplinebasis(i,P,t) ≈ bsplinebasis(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis(i,P,t) ≈ bsplinebasis(i,fP,t) for i in 1:n, t in k])
+
+        @test prod([bsplinebasis′₊₀(i,P,t) ≈ bsplinebasis′₊₀(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis′₊₀(i,P,t) ≈ bsplinebasis′₊₀(i,fP,t) for i in 1:n, t in k])
+        @test prod([bsplinebasis′₋₀(i,P,t) ≈ bsplinebasis′₋₀(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis′₋₀(i,P,t) ≈ bsplinebasis′₋₀(i,fP,t) for i in 1:n, t in k])
+        @test prod([bsplinebasis′(i,P,t) ≈ bsplinebasis′(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis′(i,P,t) ≈ bsplinebasis′(i,fP,t) for i in 1:n, t in k])
+    end
+
+    @testset "2nd degree basis" begin
+        p = 2
+        k = Knots(rand(10))+3*Knots(0,1)
+        P = BSplineSpace(p,k)
+        fP = FastBSplineSpace(p,k)
+        n = dim(P)
+        @test prod([bsplinebasis₊₀(i,P,t) ≈ bsplinebasis₊₀(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis₊₀(i,P,t) ≈ bsplinebasis₊₀(i,fP,t) for i in 1:n, t in k])
+        @test prod([bsplinebasis₋₀(i,P,t) ≈ bsplinebasis₋₀(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis₋₀(i,P,t) ≈ bsplinebasis₋₀(i,fP,t) for i in 1:n, t in k])
+        @test prod([bsplinebasis(i,P,t) ≈ bsplinebasis(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis(i,P,t) ≈ bsplinebasis(i,fP,t) for i in 1:n, t in k])
+
+        @test prod([bsplinebasis′₊₀(i,P,t) ≈ bsplinebasis′₊₀(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis′₊₀(i,P,t) ≈ bsplinebasis′₊₀(i,fP,t) for i in 1:n, t in k])
+        @test prod([bsplinebasis′₋₀(i,P,t) ≈ bsplinebasis′₋₀(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis′₋₀(i,P,t) ≈ bsplinebasis′₋₀(i,fP,t) for i in 1:n, t in k])
+        @test prod([bsplinebasis′(i,P,t) ≈ bsplinebasis′(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis′(i,P,t) ≈ bsplinebasis′(i,fP,t) for i in 1:n, t in k])
+    end
+
+    @testset "3rd degree basis" begin
+        p = 3
+        k = Knots(rand(10))+4*Knots(0,1)
+        P = BSplineSpace(p,k)
+        fP = FastBSplineSpace(p,k)
+        n = dim(P)
+        @test prod([bsplinebasis₊₀(i,P,t) ≈ bsplinebasis₊₀(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis₊₀(i,P,t) ≈ bsplinebasis₊₀(i,fP,t) for i in 1:n, t in k])
+        @test prod([bsplinebasis₋₀(i,P,t) ≈ bsplinebasis₋₀(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis₋₀(i,P,t) ≈ bsplinebasis₋₀(i,fP,t) for i in 1:n, t in k])
+        @test prod([bsplinebasis(i,P,t) ≈ bsplinebasis(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis(i,P,t) ≈ bsplinebasis(i,fP,t) for i in 1:n, t in k])
+
+        @test prod([bsplinebasis′₊₀(i,P,t) ≈ bsplinebasis′₊₀(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis′₊₀(i,P,t) ≈ bsplinebasis′₊₀(i,fP,t) for i in 1:n, t in k])
+        @test prod([bsplinebasis′₋₀(i,P,t) ≈ bsplinebasis′₋₀(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis′₋₀(i,P,t) ≈ bsplinebasis′₋₀(i,fP,t) for i in 1:n, t in k])
+        @test prod([bsplinebasis′(i,P,t) ≈ bsplinebasis′(i,fP,t) for i in 1:n, t in rand(5)])
+        @test prod([bsplinebasis′(i,P,t) ≈ bsplinebasis′(i,fP,t) for i in 1:n, t in k])
+    end
+
+    @testset "BSplineManifold-2dim" begin
         P1 = BSplineSpace(1, Knots([0, 0, 1, 1]))
         P2 = BSplineSpace(1, Knots([1, 1, 2, 3, 3]))
         n1 = dim(P1) # 2
@@ -117,7 +202,7 @@ using Test
         @test mapping(M, t) ≈ mapping(M′, t)
     end
 
-    @testset "FastBSplineManifold" begin
+    @testset "FastBSplineManifold-2dim" begin
         P1 = FastBSplineSpace(1, Knots([0, 0, 1, 1]))
         P2 = FastBSplineSpace(1, Knots([1, 1, 2, 3, 3]))
         n1 = dim(P1) # 2
@@ -137,6 +222,26 @@ using Test
         @test mapping(M, t) ≈ mapping(M′, t)
     end
 
+    @testset "BSplineSurface" begin
+        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
+        𝒂 = [[i, j] for i in 1:n1, j in 1:n2]  # n1 × n2 array of d̂ array.
+        M = BSplineSurface([P1, P2], 𝒂)
+        @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]))
+
+        @test P1 ⊆ P1′
+        @test P2 ⊆ P2′
+
+        M′ = refinement(M, [P1′, P2′])
+        t = [0.82, 1.8]
+        @test mapping(M, t) ≈ mapping(M′, t)
+    end
+
     @testset "Fitting" begin
         p1 = 2
         k1 = p1 * Knots([0, 1]) + Knots(rand(5))