Repository: hyrodium/BasicBSpline.jl
Pull request #42 (2020-12-12 14:46:23)

refactoring, update docs


Patch:
diff --git a/README.md b/README.md
index 5832d867..541a3d3a 100644
--- a/README.md
+++ b/README.md
@@ -18,9 +18,9 @@ This package provides basic (mathematical) operations for [B-spline](https://en.
     * >Currently this package's support is best for B-splines and also supports irregular grids.
     * But seems like no method for B-spline manifold.
 * [ApproXD.jl](https://github.com/floswald/ApproXD.jl)
-    * Wrapper for the dierckx Fortran library.
     * Its functions are similar to Interpolations.jl.
 * [Dierckx.jl](https://github.com/kbarbary/Dierckx.jl)
+    * Wrapper for the dierckx Fortran library.
     * Only 1-d or 2-d B-spline manifold are supported.
     * 5 or less degree of polynomial are supported.
 * **[BasicBSpline.jl](https://github.com/hyrodium/BasicBSpline.jl) (this package)**
@@ -104,9 +104,9 @@ k2 = Knots(-10:10)+p2*Knots(-10,10)
 P1 = FastBSplineSpace(p1, k1)
 P2 = FastBSplineSpace(p2, k2)
 
-f(u) = [2u[1]+sin(u[1])+cos(u[2])+u[2]/2, 3u[2]+sin(u[2])+sin(u[1])/2+u[1]^2/6]/5
+f(u1,u2) = [2u1 + sin(u1) + cos(u2) + u2 / 2, 3u2 + sin(u2) + sin(u1) / 2 + u1^2 / 6] / 5
 
-a = fittingcontrolpoints(f, [P1,P2])
+a = fittingcontrolpoints(f, P1, P2)
 M = BSplineManifold([P1,P2],a)
 save_png("fitting.png", M, unitlength=50, up=10, down=-10, left=-10, right=10)
 ```
@@ -136,9 +136,9 @@ p = 3
 k = Knots(range(-2π,2π,length=8))+p*Knots(-2π,2π)
 P = FastBSplineSpace(p, k)
 
-f(u) = [u[1],sin(u[1])]
+f(u) = [u,sin(u)]
 
-a = fittingcontrolpoints(f, [P])
+a = fittingcontrolpoints(f, P)
 M = BSplineManifold([P],a)
 save_svg("sine_curve.svg", M, unitlength=50, up=2, down=-2, left=-8, right=8)
 ```
diff --git a/docs/src/bsplineexport.jl b/docs/src/bsplineexport.jl
index 32f2c8d4..3b5f9eb7 100644
--- a/docs/src/bsplineexport.jl
+++ b/docs/src/bsplineexport.jl
@@ -36,9 +36,9 @@ k2 = Knots(-10:10) + p2 * Knots(-10, 10)
 P1 = FastBSplineSpace(p1, k1)
 P2 = FastBSplineSpace(p2, k2)
 
-f(u) = [2u[1] + sin(u[1]) + cos(u[2]) + u[2] / 2, 3u[2] + sin(u[2]) + sin(u[1]) / 2 + u[1]^2 / 6] / 5
+f(u1,u2) = [2u1 + sin(u1) + cos(u2) + u2 / 2, 3u2 + sin(u2) + sin(u1) / 2 + u1^2 / 6] / 5
 
-a = fittingcontrolpoints(f, [P1, P2])
+a = fittingcontrolpoints(f, P1, P2)
 M = BSplineManifold([P1, P2], a)
 save_png("docs/src/img/fitting.png", M, unitlength = 50, up = 10, down = -10, left = -10, right = 10)
 
@@ -61,9 +61,8 @@ p = 3
 k = Knots(range(-2π, 2π, length = 8)) + p * Knots(-2π, 2π)
 P = FastBSplineSpace(p, k)
 
-f(u) = [u[1], sin(u[1])]
-
-a = fittingcontrolpoints(f, [P])
+f(u) = [u, sin(u)]
+a = fittingcontrolpoints(f, P)
 M = BSplineManifold([P], a)
 save_svg("docs/src/img/sine_curve.svg", M, unitlength = 50, up = 2, down = -2, left = -8, right = 8)
 
diff --git a/docs/src/bsplineplots.jl b/docs/src/bsplineplots.jl
index a865ffef..64cb022f 100644
--- a/docs/src/bsplineplots.jl
+++ b/docs/src/bsplineplots.jl
@@ -58,8 +58,8 @@ plot(
 savefig("docs/src/img/subbsplineplot.png")
 
 ##
-A12 = BasicBSpline.changebasis(P1, P2)
-A13 = BasicBSpline.changebasis(P1, P3)
+A12 = changebasis(P1, P2)
+A13 = changebasis(P1, P3)
 
 plot(
     plot([t -> bsplinebasis₊₀(i, P1, t) for i in 1:dim(P1)], 1, 9, ylims = (0, 1), legend = false),
diff --git a/docs/src/img/sine_curve.svg b/docs/src/img/sine_curve.svg
index 5acd5b8f..b261835e 100644
--- a/docs/src/img/sine_curve.svg
+++ b/docs/src/img/sine_curve.svg
@@ -1,7 +1,9 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <svg xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" width="800pt" height="200pt" viewBox="0 0 800 200" version="1.1">
-<g id="surface4">
+<g id="surface11">
+<rect x="0" y="0" width="800" height="200" style="fill:rgb(100%,100%,100%);fill-opacity:1;stroke:none;"/>
 <path style="fill:none;stroke-width:2;stroke-linecap:butt;stroke-linejoin:miter;stroke:rgb(100%,0%,0%);stroke-opacity:1;stroke-miterlimit:10;" d="M 85.839844 99.535156 C 115.761719 70.269531 145.679688 43.554688 175.601562 50.191406 C 205.519531 56.832031 235.441406 96.828125 265.359375 122.097656 C 295.28125 147.371094 325.199219 157.917969 355.121094 140.371094 C 385.039062 122.820312 414.960938 77.179688 444.878906 59.628906 C 474.800781 42.082031 504.71875 52.628906 534.640625 77.902344 C 564.558594 103.171875 594.480469 143.167969 624.398438 149.808594 C 654.320312 156.445312 684.238281 129.730469 714.160156 100.464844 "/>
+<path style="fill:none;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke:rgb(30%,30%,30%);stroke-opacity:1;stroke-miterlimit:10;" d="M 85.839844 99.535156 L 115.761719 70.269531 L 175.601562 16.835938 L 265.359375 136.820312 L 355.121094 168.464844 L 444.878906 31.535156 L 534.640625 63.179688 L 624.398438 183.164062 L 684.238281 129.730469 L 714.160156 100.464844 "/>
 <path style=" stroke:none;fill-rule:nonzero;fill:rgb(10%,10%,10%);fill-opacity:1;" d="M 88.839844 99.535156 C 88.839844 103.535156 82.839844 103.535156 82.839844 99.535156 C 82.839844 95.535156 88.839844 95.535156 88.839844 99.535156 "/>
 <path style=" stroke:none;fill-rule:nonzero;fill:rgb(10%,10%,10%);fill-opacity:1;" d="M 118.761719 70.269531 C 118.761719 74.269531 112.761719 74.269531 112.761719 70.269531 C 112.761719 66.269531 118.761719 66.269531 118.761719 70.269531 "/>
 <path style=" stroke:none;fill-rule:nonzero;fill:rgb(10%,10%,10%);fill-opacity:1;" d="M 178.601562 16.835938 C 178.601562 20.835938 172.601562 20.835938 172.601562 16.835938 C 172.601562 12.835938 178.601562 12.835938 178.601562 16.835938 "/>
@@ -12,6 +14,5 @@
 <path style=" stroke:none;fill-rule:nonzero;fill:rgb(10%,10%,10%);fill-opacity:1;" d="M 627.398438 183.164062 C 627.398438 187.164062 621.398438 187.164062 621.398438 183.164062 C 621.398438 179.164062 627.398438 179.164062 627.398438 183.164062 "/>
 <path style=" stroke:none;fill-rule:nonzero;fill:rgb(10%,10%,10%);fill-opacity:1;" d="M 687.238281 129.730469 C 687.238281 133.730469 681.238281 133.730469 681.238281 129.730469 C 681.238281 125.730469 687.238281 125.730469 687.238281 129.730469 "/>
 <path style=" stroke:none;fill-rule:nonzero;fill:rgb(10%,10%,10%);fill-opacity:1;" d="M 717.160156 100.464844 C 717.160156 104.464844 711.160156 104.464844 711.160156 100.464844 C 711.160156 96.464844 717.160156 96.464844 717.160156 100.464844 "/>
-<path style="fill:none;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke:rgb(30%,30%,30%);stroke-opacity:1;stroke-miterlimit:10;" d="M 85.839844 99.535156 L 115.761719 70.269531 L 175.601562 16.835938 L 265.359375 136.820312 L 355.121094 168.464844 L 444.878906 31.535156 L 534.640625 63.179688 L 624.398438 183.164062 L 684.238281 129.730469 L 714.160156 100.464844 "/>
 </g>
 </svg>
diff --git a/docs/src/index.md b/docs/src/index.md
index ed37ab7f..33b3805c 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -14,9 +14,9 @@ This package provides basic (mathematical) operations for [B-spline](https://en.
     * >Currently this package's support is best for B-splines and also supports irregular grids.
     * But seems like no method for B-spline manifold.
 * [ApproXD.jl](https://github.com/floswald/ApproXD.jl)
-    * Wrapper for the dierckx Fortran library.
     * Its functions are similar to Interpolations.jl.
 * [Dierckx.jl](https://github.com/kbarbary/Dierckx.jl)
+    * Wrapper for the dierckx Fortran library.
     * Only 1-d or 2-d B-spline manifold are supported.
     * 5 or less degree of polynomial are supported.
 * **[BasicBSpline.jl](https://github.com/hyrodium/BasicBSpline.jl) (this package)**
@@ -118,9 +118,9 @@ k2 = Knots(-10:5:10)+p2*Knots(-10,10)
 P1 = FastBSplineSpace(p1, k1)
 P2 = FastBSplineSpace(p2, k2)
 
-f(u) = [2u[1]+sin(u[1])+cos(u[2])+u[2]/2, 3u[2]+sin(u[2])+sin(u[1])/2+u[1]^2/6]/5
+f(u1,u2) = [2u1 + sin(u1) + cos(u2) + u2 / 2, 3u2 + sin(u2) + sin(u1) / 2 + u1^2 / 6] / 5
 
-a = fittingcontrolpoints(f, [P1,P2])
+a = fittingcontrolpoints(f, P1, P2)
 M = BSplineManifold([P1,P2],a)
 save_png("docs/src/img/fitting_coarse.png", M, unitlength=50, up=10, down=-10, left=-10, right=10)
 ```
@@ -132,9 +132,9 @@ p = 3
 k = Knots(range(-2π,2π,length=8))+p*Knots(-2π,2π)
 P = FastBSplineSpace(p, k)
 
-f(u) = [u[1],sin(u[1])]
+f(u) = [u,sin(u)]
 
-a = fittingcontrolpoints(f, [P])
+a = fittingcontrolpoints(f, P)
 M = BSplineManifold([P],a)
 save_svg("docs/src/img/sine_curve.svg", M, unitlength=50, up=2, down=-2, left=-8, right=8)
 ```
diff --git a/docs/src/math.md b/docs/src/math.md
index efd797d8..e9ec080d 100644
--- a/docs/src/math.md
+++ b/docs/src/math.md
@@ -438,6 +438,18 @@ save_png("2dim.png", M) # save image
 ```
 ![](img/2dim.png)
 
+```julia
+points = [M([u,v]) for u in range(3.0,6.0,length=50), v in range(3.0,6.0,length=50)]
+X = [point[1] for point in points]
+Y = [point[2] for point in points]
+scene = Scene(resolution=(1000,1000))
+Makie.surface!(X,Y)
+Makie.xlims!(-5,5)
+Makie.ylims!(-5,5)
+save("2dim_makie.png", scene)
+```
+
+![](img/2dim_makie.png)
 
 ## Affine commutativity
 !!! info "Thm.  Affine commutativity"
diff --git a/src/_BSplineBasis.jl b/src/_BSplineBasis.jl
index 9df1285a..23ca4e92 100644
--- a/src/_BSplineBasis.jl
+++ b/src/_BSplineBasis.jl
@@ -18,7 +18,7 @@ Right-sided limit version.
 \end{aligned}
 ```
 """
-function bsplinebasis₊₀(P::BSplineSpace, t::Real)::Array{Float64,1}
+function bsplinebasis₊₀(P::BSplineSpace, t::Real)::Vector{Float64}
     p = P.degree
     k = P.knots
 
@@ -49,7 +49,7 @@ Left-sided limit version.
 \end{aligned}
 ```
 """
-function bsplinebasis₋₀(P::BSplineSpace, t::Real)::Array{Float64,1}
+function bsplinebasis₋₀(P::BSplineSpace, t::Real)::Vector{Float64}
     p = P.degree
     k = P.knots
 
@@ -81,7 +81,7 @@ Modified version.
 \end{aligned}
 ```
 """
-function bsplinebasis(P::BSplineSpace, t::Real)::Array{Float64,1}
+function bsplinebasis(P::BSplineSpace, t::Real)::Vector{Float64}
     p = P.degree
     k = P.knots
 
@@ -156,7 +156,7 @@ Right-sided limit 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)
 ```
 """
-function bsplinebasis′₊₀(P::BSplineSpace, t)::Array{Float64,1}
+function bsplinebasis′₊₀(P::BSplineSpace, t)::Vector{Float64}
     p = P.degree
     k = P.knots
 
@@ -177,7 +177,7 @@ Left-sided limit 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)
 ```
 """
-function bsplinebasis′₋₀(P::BSplineSpace, t)::Array{Float64,1}
+function bsplinebasis′₋₀(P::BSplineSpace, t)::Vector{Float64}
     p = P.degree
     k = P.knots
 
@@ -200,7 +200,7 @@ Modified version.
 """
 bsplinebasis′
 
-function bsplinebasis′(P::BSplineSpace, t)::Array{Float64,1}
+function bsplinebasis′(P::BSplineSpace, t)::Vector{Float64}
     p = P.degree
     k = P.knots
 
diff --git a/src/_BSplineManifold.jl b/src/_BSplineManifold.jl
index 01743274..aa5076cf 100644
--- a/src/_BSplineManifold.jl
+++ b/src/_BSplineManifold.jl
@@ -23,16 +23,16 @@ end
 B-spline manifold for general polynomial degree
 """
 struct BSplineManifold <: AbstractBSplineManifold
-    bsplinespaces::Array{BSplineSpace,1}
+    bsplinespaces::Vector{BSplineSpace}
     controlpoints::Array{Float64}
     function BSplineManifold(Ps::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{<:Real})
         Ps = BSplineSpace.(Ps)
         if collect(size(a)[1:end-1]) ≠ dim.(Ps)
-            error("dimension does not match")
+            throw(DimensionMismatch())
         else
-            P = convert(Array{BSplineSpace,1}, Ps)
-            a = convert(Array{Float64}, a)
-            new(P, a)
+            P = convert(Vector{BSplineSpace}, Ps)
+            a′ = convert(Array{Float64}, a)
+            new(P, a′)
         end
     end
 end
diff --git a/src/_BSplineSpace.jl b/src/_BSplineSpace.jl
index 445329b9..cf5eefd9 100644
--- a/src/_BSplineSpace.jl
+++ b/src/_BSplineSpace.jl
@@ -13,7 +13,7 @@ struct BSplineSpace <: AbstractBSplineSpace
     knots::Knots
     function BSplineSpace(degree::Integer, knots::Knots)
         if degree < 0
-            error("degree of polynominal must be non-negative")
+            throw(DomainError(p, "degree of polynominal must be non-negative"))
         end
         new(degree, knots)
     end
@@ -126,4 +126,4 @@ function _knotindex(P,t)
     p = degree(P)
     k = knots(P)
     return _knotindex(k[1+p:end-p], t)+p
-end
\ No newline at end of file
+end
diff --git a/src/_FastBSplineManifold.jl b/src/_FastBSplineManifold.jl
index d0e3f377..81ef1cd7 100644
--- a/src/_FastBSplineManifold.jl
+++ b/src/_FastBSplineManifold.jl
@@ -5,16 +5,16 @@ B-spline manifold for lower polynomial degree
 TODO: make the field `bsplinespaces` to be conposite type, not abstract type, for performance
 """
 struct FastBSplineManifold <: AbstractBSplineManifold
-    bsplinespaces::Array{<:FastBSplineSpace,1}
+    bsplinespaces::Vector{<:FastBSplineSpace}
     controlpoints::Array{Float64}
     function FastBSplineManifold(Ps::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{<:Real})
         Ps = FastBSplineSpace.(Ps)
         if collect(size(a)[1:end-1]) ≠ dim.(Ps)
-            error("dimension does not match")
+            throw(DimensionMismatch())
         else
-            P = convert(Array{FastBSplineSpace,1}, Ps)
-            a = convert(Array{Float64}, a)
-            new(P, a)
+            P = convert(Vector{FastBSplineSpace}, Ps)
+            a′ = convert(Array{Float64}, a)
+            new(P, a′)
         end
     end
 end
@@ -22,25 +22,28 @@ end
 struct BSplineCurve{p1} <: AbstractBSplineManifold
     bsplinespace1::FastBSplineSpace{p1}
     controlpoints::Array{Float64,2}
-    function BSplineCurve(P1::AbstractBSplineSpace, a::AbstractArray{<:Real})
+    function BSplineCurve(P1::AbstractBSplineSpace, a::AbstractArray{<:Real,2})
         p1 = degree(P1)
         if size(a)[1:end-1] ≠ (dim(P1),)
-            error("dimension does not match")
+            throw(DimensionMismatch())
         else
-            a = convert(Array{Float64}, a)
-            new{p1}(P1, a)
+            a′ = convert(Array{Float64}, a)
+            new{p1}(P1, a′)
         end
     end
 end
+function BSplineCurve(P1,::AbstractBSplineSpace, a::AbstractArray{<:AbstractVector{<:Real},1})
+    a′ = _arrayofvector2array(a)
+    return BSplineCurve(P1, a′)
+end
 function BSplineCurve(Ps::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{<:Real})
-    # p1, k1 = degree(Ps[1]), knots(Ps[1])
     if collect(size(a)[1:end-1]) ≠ dim.(Ps)
-        error("dimension does not match")
-    # elseif p1 > MAX_DEGREE
-    #     return BSplineManifold(Ps, a)
+        throw(DimensionMismatch())
+    elseif length(Ps) ≠ 1
+        throw(DimensionMismatch())
     else
-        a = convert(Array{Float64}, a)
-        return BSplineCurve(Ps[1],a)
+        a′ = convert(Array{Float64}, a)
+        return BSplineCurve(Ps[1],a′)
     end
 end
 function BSplineCurve{q1}(Ps::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{<:Real}) where {q1}
@@ -51,26 +54,28 @@ struct BSplineSurface{p1,p2} <: AbstractBSplineManifold
     bsplinespace1::FastBSplineSpace{p1}
     bsplinespace2::FastBSplineSpace{p2}
     controlpoints::Array{Float64,3}
-    function BSplineSurface(P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, a::AbstractArray{<:Real})
+    function BSplineSurface(P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, a::AbstractArray{<:Real,3})
         p1, p2 = degree(P1), degree(P2)
         if size(a)[1:end-1] ≠ (dim(P1), dim(P2))
-            error("dimension does not match")
+            throw(DimensionMismatch())
         else
-            a = convert(Array{Float64}, a)
-            new{p1,p2}(P1, P2, a)
+            a′ = convert(Array{Float64}, a)
+            new{p1,p2}(P1, P2, a′)
         end
     end
 end
+function BSplineSurface(P1,::AbstractBSplineSpace, P2,::AbstractBSplineSpace, a::AbstractArray{<:AbstractVector{<:Real},2})
+    a′ = _arrayofvector2array(a)
+    return BSplineSurface(P1, P2, a′)
+end
 function BSplineSurface(Ps::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{<:Real})
-    # p1, k1 = degree(Ps[1]), knots(Ps[1])
-    # p2, k2 = degree(Ps[2]), knots(Ps[2])
     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)
+        throw(DimensionMismatch())
+    elseif length(Ps) ≠ 2
+        throw(DimensionMismatch())
     else
-        a = convert(Array{Float64}, a)
-        return BSplineSurface(Ps[1], Ps[2], a)
+        a′ = convert(Array{Float64}, a)
+        return BSplineSurface(Ps[1], Ps[2], a′)
     end
 end
 function BSplineSurface{q1,q2}(Ps::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{<:Real}) where {q1} where {q2}
@@ -85,61 +90,64 @@ struct BSplineSolid{p1,p2,p3} <: AbstractBSplineManifold
     function BSplineSolid(P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, P3::AbstractBSplineSpace, a::AbstractArray{<:Real})
         p1, p2, p3 = degree(P1), degree(P2), degree(P3)
         if size(a)[1:end-1] ≠ (dim(P1), dim(P2), dim(P3))
-            error("dimension does not match")
+            throw(DimensionMismatch())
         else
-            a = convert(Array{Float64}, a)
-            new{p1,p2,p3}(P1, P2, P3, a)
+            a′ = convert(Array{Float64}, a)
+            new{p1,p2,p3}(P1, P2, P3, a′)
         end
     end
 end
+function BSplineSolid(P1,::AbstractBSplineSpace, P2,::AbstractBSplineSpace, P3,::AbstractBSplineSpace, a::AbstractArray{<:AbstractVector{<:Real},3})
+    a′ = _arrayofvector2array(a)
+    return BSplineSolid(P1, P2, P3, a′)
+end
 function BSplineSolid(Ps::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{<:Real})
-    # p1, k1 = degree(Ps[1]), knots(Ps[1])
-    # p2, k2 = degree(Ps[2]), knots(Ps[2])
-    # p3, k3 = degree(Ps[3]), knots(Ps[3])
     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)
+        throw(DimensionMismatch())
+    elseif length(Ps) ≠ 3
+        throw(DimensionMismatch())
     else
-        a = convert(Array{Float64}, a)
-        BSplineSolid(Ps[1],Ps[2],Ps[3],a)
+        a′ = convert(Array{Float64}, a)
+        BSplineSolid(Ps[1],Ps[2],Ps[3],a′)
     end
 end
 function BSplineSolid{q1,q2,q3}(Ps::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{<:Real}) where {q1} where {q2} where {q3}
     return BSplineSolid(Ps, a)
 end
 
+function _arrayofvector2array(a::AbstractArray{<:AbstractVector{R}, d})::Array{R,d+1} where d where R<:Real
+    d̂ = length(a[1])
+    s = size(a)
+    N = prod(s)
+    a_2dim = [a[i][j] for i in 1:N, j in 1:d̂]
+    a′ = reshape(a_2dim, s..., d̂)
+    return a′
+end
+
 for fname in (:BSplineCurve, :BSplineSurface, :BSplineSolid, :FastBSplineManifold, :BSplineManifold)
-    @eval function $fname(Ps::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{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)
+    @eval function $fname(Ps::AbstractVector{<:AbstractBSplineSpace}, a::AbstractArray{<:AbstractVector{<:Real}})
+        a′ = _arrayofvector2array(a)
+        return $fname(Ps, a′)
     end
 end
 
 function BSplineManifold(M::BSplineCurve{p1}) where {p1}
-    k1 = M.knots1
-    P1 = BSplineSpace(p1, k1)
+    P1 = M.bsplinespace1
     a = M.controlpoints
     return BSplineManifold([P1], a)
 end
 
 function BSplineManifold(M::BSplineSurface{p1,p2}) where {p1} where {p2}
-    k1 = M.knots1
-    P1 = BSplineSpace(p1, k1)
-    k2 = M.knots2
-    P2 = BSplineSpace(p2, k2)
+    k1 = M.bsplinespace1
+    P2 = M.bsplinespace2
     a = M.controlpoints
     return BSplineManifold([P1, P2], a)
 end
 
 function BSplineManifold(M::BSplineSolid{p1,p2,p3}) where {p1} where {p2} where {p3}
-    k1 = M.knots1
-    P1 = BSplineSpace(p1, k1)
-    k2 = M.knots2
-    P2 = BSplineSpace(p2, k2)
-    k3 = M.knots3
-    P3 = BSplineSpace(p3, k3)
+    P1 = M.bsplinespace1
+    P2 = M.bsplinespace2
+    P3 = M.bsplinespace3
     a = M.controlpoints
     return BSplineManifold([P1, P2, P3], a)
 end
@@ -156,11 +164,10 @@ function (M::FastBSplineManifold)(t::AbstractVector{<:Real})
     return BSplineManifold(M)(t)
 end
 
-function (M::BSplineCurve{p1})(t::AbstractVector{<:Real}) where {p1}
+function (M::BSplineCurve{p1})(t1::Real) where {p1}
     P1, = bsplinespaces(M)
     a = controlpoints(M)
     n1 = dim(P1)
-    t1, = t
     k1 = P1.knots
     j1 = _knotindex(P1, t1)
     b1 = _bsplinebasis(P1,t1,j1)
@@ -179,11 +186,10 @@ function (M::BSplineCurve{p1})(t::AbstractVector{<:Real}) where {p1}
     return S1
 end
 
-function (M::BSplineSurface{p1,p2})(t::AbstractVector{<:Real}) where {p1} where {p2}
+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)
-    t1, t2 = t
     k1, k2 = P1.knots, P2.knots
     j1, j2 = _knotindex(P1, t1), _knotindex(P2, t2)
     b1 = _bsplinebasis(P1,t1,j1)
@@ -212,11 +218,10 @@ function (M::BSplineSurface{p1,p2})(t::AbstractVector{<:Real}) where {p1} where
     return S1
 end
 
-function (M::BSplineSolid{p1,p2,p3})(t::AbstractVector{<:Real}) where {p1} where {p2} where {p3}
+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)
-    t1, t2, t3 = t
     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)
@@ -263,6 +268,22 @@ function (M::BSplineSolid{p1,p2,p3})(t::AbstractVector{<:Real}) where {p1} where
     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
diff --git a/src/_FastBSplineSpace.jl b/src/_FastBSplineSpace.jl
index 8275f8dd..4044ed15 100644
--- a/src/_FastBSplineSpace.jl
+++ b/src/_FastBSplineSpace.jl
@@ -7,9 +7,9 @@ struct FastBSplineSpace{p} <: AbstractBSplineSpace
     knots::Knots
     function FastBSplineSpace(p::Integer, knots::Knots)
         if p < 0
-            error("degree of polynominal must be non-negative")
+            throw(DomainError(p, "degree of polynominal must be non-negative"))
         elseif p > MAX_DEGREE
-            error("FastBSpline supports only degree 0 , ... , $(MAX_DEGREE)")
+            throw(DomainError(p, "FastBSpline supports only degree 0,...,$(MAX_DEGREE)"))
         end
         new{p}(knots)
     end
@@ -25,12 +25,12 @@ 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
+# """
+# convert AbstractBSplineSpace to FastBSplineSpace
+# """
+# function Base.convert(::Type{<:FastBSplineSpace}, P::AbstractBSplineSpace)
+#     return FastBSplineSpace(degree(P), knots(P))
+# end
 
 function degree(P::FastBSplineSpace{p}) where {p}
     return p
diff --git a/src/_Fitting.jl b/src/_Fitting.jl
index 71151557..4b1adc1a 100644
--- a/src/_Fitting.jl
+++ b/src/_Fitting.jl
@@ -104,7 +104,7 @@ end
 function _f_b_int_R(func, i1, P1::AbstractBSplineSpace, nip1, nodes1, weights1)
     k1 = knots(P1)
     p1 = degree(P1)
-    F(t1) = bsplinebasis(i1, P1, t1) * func([t1])
+    F(t1) = bsplinebasis(i1, P1, t1) * func(t1)
 
     f1,l1 = i1, i1+p1
 
@@ -119,7 +119,7 @@ function _f_b_int_I(func, i1, P1::AbstractBSplineSpace, nip1, nodes1, weights1)
     k1 = knots(P1)
     m1 = length(k1)
     p1 = degree(P1)
-    F(t1) = bsplinebasis(i1, P1, t1) * func([t1])
+    F(t1) = bsplinebasis(i1, P1, t1) * func(t1)
 
     f1,l1 = max(i1, p1+1), min(i1+p1, m1-p1-1)
 
@@ -133,7 +133,7 @@ end
 function _f_b_int_R(func, i1, i2, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, nip1, nip2, nodes1, nodes2, weights1, weights2)
     k1, k2 = knots(P1), knots(P2)
     p1, p2 = degree(P1), degree(P2)
-    F(t1, t2) = bsplinebasis(i1, P1, t1) * bsplinebasis(i2, P2, t2) * func([t1, t2])
+    F(t1, t2) = bsplinebasis(i1, P1, t1) * bsplinebasis(i2, P2, t2) * func(t1, t2)
 
     f1, l1 = i1, i1+p1
     f2, l2 = i2, i2+p2
@@ -157,7 +157,7 @@ function _f_b_int_I(func, i1, i2, P1::AbstractBSplineSpace, P2::AbstractBSplineS
     k1, k2 = knots(P1), knots(P2)
     m1, m2 = length(k1), length(k2)
     p1, p2 = degree(P1), degree(P2)
-    F(t1, t2) = bsplinebasis(i1, P1, t1) * bsplinebasis(i2, P2, t2) * func([t1, t2])
+    F(t1, t2) = bsplinebasis(i1, P1, t1) * bsplinebasis(i2, P2, t2) * func(t1, t2)
 
     f1, l1 = max(i1, p1+1), min(i1+p1, m1-p1-1)
     f2, l2 = max(i2, p2+1), min(i2+p2, m2-p2-1)
@@ -180,7 +180,7 @@ end
 function _f_b_int_R(func, i1, i2, i3, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, P3::AbstractBSplineSpace, nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3)
     k1, k2, k3 = knots(P1), knots(P2), knots(P3)
     p1, p2, p3 = degree(P1), degree(P2), degree(P3)
-    F(t1, t2, t3) = bsplinebasis(i1, P1, t1) * bsplinebasis(i2, P2, t2) * bsplinebasis(i3, P3, t3) * func([t1, t2, t3])
+    F(t1, t2, t3) = bsplinebasis(i1, P1, t1) * bsplinebasis(i2, P2, t2) * bsplinebasis(i3, P3, t3) * func(t1, t2, t3)
 
     f1, l1 = i1, i1+p1
     f2, l2 = i2, i2+p2
@@ -221,7 +221,7 @@ function _f_b_int_I(func, i1, i2, i3, P1::AbstractBSplineSpace, P2::AbstractBSpl
     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(i1, P1, t1) * bsplinebasis(i2, P2, t2) * bsplinebasis(i3, P3, t3) * func([t1, t2, t3])
+    F(t1, t2, t3) = bsplinebasis(i1, P1, t1) * bsplinebasis(i2, P2, t2) * bsplinebasis(i3, P3, t3) * func(t1, t2, t3)
 
     f1, l1 = max(i1, p1+1), min(i1+p1, m1-p1-1)
     f2, l2 = max(i2, p2+1), min(i2+p2, m2-p2-1)
@@ -325,15 +325,18 @@ function innerproduct_I(func, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace
 end
 
 """
-* func: Array{Real,1} -> ℝ-linear space
+* func: Real -> ℝ-vector space
 """
-function fittingcontrolpoints_1dim(func, P1::AbstractBSplineSpace)
+function fittingcontrolpoints(func, P1::AbstractBSplineSpace)
     b = innerproduct_I(func, P1)
     A = innerproduct_I(P1)
     return inv(A) * b
 end
 
-function fittingcontrolpoints_2dim(func, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace)
+"""
+* func: (Real,Real) -> ℝ-vector space
+"""
+function fittingcontrolpoints(func, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace)
     n1, n2 = dim(P1), dim(P2)
     A1, A2 = innerproduct_I(P1), innerproduct_I(P2)
     A = [A1[i1, j1] * A2[i2, j2] for i1 in 1:n1, i2 in 1:n2, j1 in 1:n1, j2 in 1:n2]
@@ -343,7 +346,10 @@ function fittingcontrolpoints_2dim(func, P1::AbstractBSplineSpace, P2::AbstractB
     return reshape(inv(_A) * _b, n1, n2)
 end
 
-function fittingcontrolpoints_3dim(func, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, P3::AbstractBSplineSpace)
+"""
+* func: (Real,Real,Real) -> ℝ-vector space
+"""
+function fittingcontrolpoints(func, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, P3::AbstractBSplineSpace)
     n1, n2, n3 = dim(P1), dim(P2), dim(P3)
     A1, A2, A3 = innerproduct_I(P1), innerproduct_I(P2), innerproduct_I(P3)
     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]
@@ -357,17 +363,10 @@ 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::Array{<:AbstractBSplineSpace,1})
-    # TODO: currently, this function only supports for 1-dim and 2-dim B-spline manifold.
-    d = length(Ps)
-    if d == 1
-        return fittingcontrolpoints_1dim(func, Ps[1])
-    elseif d == 2
-        return fittingcontrolpoints_2dim(func, Ps[1], Ps[2])
-    elseif d == 3
-        return fittingcontrolpoints_3dim(func, Ps[1], Ps[2], Ps[3])
-    else
-        error("fittingcontrolpoints supports only 1 to 3 dimensions.")
-    end
+function fittingcontrolpoints(func, Ps::AbstractVector{<:AbstractBSplineSpace})
+    # 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...)
 end
diff --git a/src/_Integral.jl b/src/_Integral.jl
index 21822208..444ebb5e 100644
--- a/src/_Integral.jl
+++ b/src/_Integral.jl
@@ -14,6 +14,9 @@ function integrate(func, I1::ClosedInterval{<:Real}, nip1, nodes1, weights1)
     return S1 * width(I1) / 2
 end
 
+"""
+Integrate on rectangular region.
+"""
 function integrate(func, I1::ClosedInterval{<:Real}, I2::ClosedInterval{<:Real}, nip1, nip2, nodes1, nodes2, weights1, weights2)
     dnodes1 = (width(I1) * nodes1 .+ sum(extrema(I1))) / 2
     dnodes2 = (width(I2) * nodes2 .+ sum(extrema(I2))) / 2
@@ -34,6 +37,9 @@ function integrate(func, I1::ClosedInterval{<:Real}, I2::ClosedInterval{<:Real},
     return S1 * width(I1) * width(I2) / 4
 end
 
+"""
+Integrate on rectangular solid
+"""
 function integrate(func, I1::ClosedInterval{<:Real}, I2::ClosedInterval{<:Real}, I3::ClosedInterval{<:Real}, nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3)
     dnodes1 = (width(I1) * nodes1 .+ sum(extrema(I1))) / 2
     dnodes2 = (width(I2) * nodes2 .+ sum(extrema(I2))) / 2
diff --git a/src/_Knots.jl b/src/_Knots.jl
index db2f564d..71f4d1e3 100644
--- a/src/_Knots.jl
+++ b/src/_Knots.jl
@@ -7,12 +7,12 @@ k=(k_1,\dots,k_l)
 ```
 """
 struct Knots
-    vector::Array{Float64,1}
+    vector::Vector{Float64}
     function Knots(vector::AbstractVector{<:Real})
         return new(sort(convert(Vector{Float64}, vector)))
     end
 end
-function Knots(vector::Array{Any,1})
+function Knots(vector::AbstractVector)
     if isempty(vector)
         return Knots(Float64[])
     else
@@ -32,7 +32,7 @@ function Base.:*(p₊::Integer, k::Knots)
     elseif p₊ > 0
         return sum(k for _ in 1:p₊)
     else
-        error("Polynominal degree p₊ must be non-negative.")
+        throw(DomainError(p₊, "additional degree of polynominal must be non-negative"))
     end
 end
 
@@ -47,6 +47,13 @@ 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)
 
+@doc raw"""
+Check a inclusive relation ship ``k\subset k'``.
+```math
+(1,2,3) \subseteq (1,2,3,4)
+(1,2,5) \nsubseteq (1,2,3,4)
+```
+"""
 function Base.issubset(k::Knots, k′::Knots)
     K′ = copy(k′.vector)
     for kᵢ in k.vector
diff --git a/src/_Refinement.jl b/src/_Refinement.jl
index 5739615b..ddb5aa19 100644
--- a/src/_Refinement.jl
+++ b/src/_Refinement.jl
@@ -9,7 +9,7 @@ B_{(i,p,k)} = \sum_{j}A_{i,j}B_{(j,p',k')}
 Assumption:
 * ``P ⊆ P′``
 """
-function _changebasis_R(P::AbstractBSplineSpace, P′::AbstractBSplineSpace)::Array{Float64,2}
+function _changebasis_R(P::AbstractBSplineSpace, P′::AbstractBSplineSpace)::Matrix{Float64}
     p = degree(P)
     k = knots(P)
     p′ = degree(P′)
@@ -140,7 +140,7 @@ B_{(i,p,k)} = \sum_{j}A_{i,j}B_{(j,p',k')}
 Assumption:
 * ``P ⊑ P′``
 """
-function _changebasis_I(P::AbstractBSplineSpace, P′::AbstractBSplineSpace)::Array{Float64,2}
+function _changebasis_I(P::AbstractBSplineSpace, P′::AbstractBSplineSpace)::Matrix{Float64}
     I = bsplineunity(P)
     p = degree(P)
     k = knots(P)
@@ -170,7 +170,7 @@ function changebasis(P::AbstractBSplineSpace, P′::AbstractBSplineSpace)
     elseif P ⊑ P′
         return _changebasis_I(P, P′)
     else
-        error("𝒫[p,k] ⊄ 𝒫[p′,k′] and 𝒫[p,k] ⋢ 𝒫[p′,k′]")
+        throw(DomainError((P, P′),"𝒫[p,k] ⊆ 𝒫[p′,k′] or 𝒫[p,k] ⊑ 𝒫[p′,k′] must hold."))
     end
 end
 
@@ -214,12 +214,12 @@ function refinement(M::AbstractBSplineManifold; p₊::Union{Nothing,AbstractVect
     if isnothing(p₊)
         p₊ = zeros(Int, d)
     elseif length(Ps) ≠ length(p₊)
-        error("dimension does not match")
+        throw(DimensionMismatch())
     end
     if isnothing(k₊)
         k₊ = zeros(Knots, d)
     elseif length(Ps) ≠ length(k₊)
-        error("dimension does not match")
+        throw(DimensionMismatch())
     end
 
     Ps′ = [(P = Ps[i];


Test patch:
diff --git a/test/runtests.jl b/test/runtests.jl
index 74fbb03f..e57993ff 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -360,7 +360,7 @@ using Random
         a_ref = M′.controlpoints
 
         a_tmp = fittingcontrolpoints(M, [P1′, P2′])
-        a_fit = reshape(transpose(hcat(reshape(a_tmp, prod(size(a_tmp)))...)), size(a_tmp)..., 3)
+        a_fit = BasicBSpline._arrayofvector2array(a_tmp)
 
         @test norm(a_fit - a_ref) < ε
     end
@@ -390,7 +390,7 @@ using Random
         a_ref = M′.controlpoints
 
         a_tmp = fittingcontrolpoints(M, [P1′, P2′])
-        a_fit = reshape(transpose(hcat(reshape(a_tmp, prod(size(a_tmp)))...)), size(a_tmp)..., 3)
+        a_fit = BasicBSpline._arrayofvector2array(a_tmp)
 
         @test norm(a_fit - a_ref) < ε
     end
@@ -427,7 +427,7 @@ using Random
         a_ref = M′.controlpoints
 
         a_tmp = fittingcontrolpoints(M, [P1′, P2′, P3′])
-        a_fit = reshape(transpose(hcat(reshape(a_tmp, prod(size(a_tmp)))...)), size(a_tmp)..., 4)
+        a_fit = BasicBSpline._arrayofvector2array(a_tmp)
 
         @test norm(a_fit - a_ref) < ε
     end
@@ -464,7 +464,7 @@ using Random
         a_ref = M′.controlpoints
 
         a_tmp = fittingcontrolpoints(M, [P1′, P2′, P3′])
-        a_fit = reshape(transpose(hcat(reshape(a_tmp, prod(size(a_tmp)))...)), size(a_tmp)..., 4)
+        a_fit = BasicBSpline._arrayofvector2array(a_tmp)
 
         @test norm(a_fit - a_ref) < 1.0e-6
     end