Repository: hyrodium/BasicBSpline.jl
Pull request #84 (2021-03-09 10:00:01)

Feature/moreperformance


Patch:
diff --git a/src/BasicBSpline.jl b/src/BasicBSpline.jl
index 0daedd7f..5c41afa1 100644
--- a/src/BasicBSpline.jl
+++ b/src/BasicBSpline.jl
@@ -29,7 +29,6 @@ include("_FastBSplineManifold.jl")
 include("_Refinement.jl")
 include("_Integral.jl")
 include("_Fitting.jl")
-include("_Geometry.jl")
 
 if VERSION < v"1.2.0"
     <=(x) = Fix2(<=, x)
diff --git a/src/_BSplineSpace.jl b/src/_BSplineSpace.jl
index ec11ac9c..0a392667 100644
--- a/src/_BSplineSpace.jl
+++ b/src/_BSplineSpace.jl
@@ -126,5 +126,6 @@ end
 function _knotindex(P,t)
     p = degree(P)
     k = knots(P)
-    return _knotindex(k[1+p:end-p], t)+p
+    l = length(k)
+    return _knotindex(view(k.vector, 1+p:l-p), t) + p
 end
diff --git a/src/_Fitting.jl b/src/_Fitting.jl
index b4bd12bc..12ecd5f9 100644
--- a/src/_Fitting.jl
+++ b/src/_Fitting.jl
@@ -9,12 +9,12 @@ Assumption:
 * ``i ≤ j``
 * ``1 ≤ n ≤ p-j+i+1``
 """
-function _b_b_int(P::AbstractBSplineSpace, i, j, n, nip, nodes, weights)
+function _b_b_int(P::AbstractBSplineSpace, i, j, n, gl)
     p = degree(P)
     I = knots(P)[j+n-1]..knots(P)[j+n]
 
     f(t) = bsplinebasis(i, P, t) * bsplinebasis(j, P, t)
-    return integrate(f, I, nip, nodes, weights)
+    return integrate(f, I, gl)
 end
 
 @doc raw"""
@@ -30,7 +30,7 @@ Calculate
 \end{aligned}
 ```
 """
-function _b_b_int_R(P::AbstractBSplineSpace, i, j, nip, nodes, weights)
+function _b_b_int_R(P::AbstractBSplineSpace, i, j, gl)
     p = degree(P)
     k = knots(P)
     Δ = j - i
@@ -41,13 +41,13 @@ function _b_b_int_R(P::AbstractBSplineSpace, i, j, nip, nodes, weights)
     elseif Δ ≥ 0  # i ≤ j
         s = 0.0
         for n in 1:p-j+i+1
-            s += _b_b_int(P, i, j, n, nip, nodes, weights)
+            s += _b_b_int(P, i, j, n, gl)
         end
         return s
     else  # j < i
         s = 0.0
         for n in 1:p-i+j+1
-            s += _b_b_int(P, j, i, n, nip, nodes, weights)
+            s += _b_b_int(P, j, i, n, gl)
         end
         return s
     end
@@ -61,7 +61,7 @@ Calculate
 \end{aligned}
 ```
 """
-function _b_b_int_I(P::AbstractBSplineSpace, i, j, nip, nodes, weights)
+function _b_b_int_I(P::AbstractBSplineSpace, i, j, gl)
     p = degree(P)
     k = knots(P)
     n = dim(P)
@@ -73,13 +73,13 @@ function _b_b_int_I(P::AbstractBSplineSpace, i, j, nip, nodes, weights)
     elseif Δ ≥ 0  # i ≤ j
         s = 0.0
         for m in max(p-j+2,1):min(n-j+1,p-j+i+1)
-            s += _b_b_int(P, i, j, m, nip, nodes, weights)
+            s += _b_b_int(P, i, j, m, gl)
         end
         return s
     else  # j < i
         s = 0.0
         for m in max(p-i+2,1):min(n-i+1,p-i+j+1)
-            s += _b_b_int(P, j, i, m, nip, nodes, weights)
+            s += _b_b_int(P, j, i, m, gl)
         end
         return s
     end
@@ -89,33 +89,33 @@ function innerproduct_R(P::AbstractBSplineSpace)
     p = degree(P)
     n = dim(P)
     nip = p + 1
-    nodes, weights = gausslegendre(nip)
-    return [_b_b_int_R(P, i, j, nip, nodes, weights) for i in 1:n, j in 1:n]
+    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)
     n = dim(P)
     nip = p + 1
-    nodes, weights = gausslegendre(nip)
-    return [_b_b_int_I(P, i, j, nip, nodes, weights) for i in 1:n, j in 1:n]
+    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, nip1, nodes1, weights1)
+function _f_b_int_R(func, i1, P1::AbstractBSplineSpace, gl1)
     k1 = knots(P1)
     p1 = degree(P1)
     F(t1) = bsplinebasis(i1, P1, t1) * func(t1)
 
     f1,l1 = i1, i1+p1
 
-    S1 = integrate(F, k1[f1]..k1[f1+1], nip1, nodes1, weights1)
+    S1 = integrate(F, k1[f1]..k1[f1+1], gl1)
     for j1 in f1+1:l1
-        S1 += integrate(F, k1[j1]..k1[j1+1], nip1, nodes1, weights1)
+        S1 += integrate(F, k1[j1]..k1[j1+1], gl1)
     end
     return S1
 end
 
-function _f_b_int_I(func, i1, P1::AbstractBSplineSpace, nip1, nodes1, weights1)
+function _f_b_int_I(func, i1, P1::AbstractBSplineSpace, gl1)
     k1 = knots(P1)
     m1 = length(k1)
     p1 = degree(P1)
@@ -123,14 +123,14 @@ function _f_b_int_I(func, i1, P1::AbstractBSplineSpace, nip1, nodes1, weights1)
 
     f1,l1 = max(i1, p1+1), min(i1+p1, m1-p1-1)
 
-    S1 = integrate(F, k1[f1]..k1[f1+1], nip1, nodes1, weights1)
+    S1 = integrate(F, k1[f1]..k1[f1+1], gl1)
     for j1 in f1+1:l1
-        S1 += integrate(F, k1[j1]..k1[j1+1], nip1, nodes1, weights1)
+        S1 += integrate(F, k1[j1]..k1[j1+1], gl1)
     end
     return S1
 end
 
-function _f_b_int_R(func, i1, i2, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, nip1, nip2, nodes1, nodes2, weights1, weights2)
+function _f_b_int_R(func, i1, i2, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, gl1, gl2)
     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)
@@ -138,22 +138,22 @@ function _f_b_int_R(func, i1, i2, P1::AbstractBSplineSpace, P2::AbstractBSplineS
     f1, l1 = i1, i1+p1
     f2, l2 = i2, i2+p2
 
-    S2 = integrate(F, k1[f1]..k1[f1+1], k2[f2]..k2[f2+1], nip1, nip2, nodes1, nodes2, weights1, weights2)
+    S2 = integrate(F, k1[f1]..k1[f1+1], k2[f2]..k2[f2+1], gl1, gl2)
     for j2 in f2+1:l2
-        S2 += integrate(F, k1[f1]..k1[f1+1], k2[j2]..k2[j2+1], nip1, nip2, nodes1, nodes2, weights1, weights2)
+        S2 += integrate(F, k1[f1]..k1[f1+1], k2[j2]..k2[j2+1], gl1, gl2)
     end
     S1 = S2
     for j1 in f1+1:l1
-        S2 = integrate(F, k1[j1]..k1[j1+1], k2[f2]..k2[f2+1], nip1, nip2, nodes1, nodes2, weights1, weights2)
+        S2 = integrate(F, k1[j1]..k1[j1+1], k2[f2]..k2[f2+1], gl1, gl2)
         for j2 in f2+1:l2
-            S2 += integrate(F, k1[j1]..k1[j1+1], k2[j2]..k2[j2+1], nip1, nip2, nodes1, nodes2, weights1, weights2)
+            S2 += integrate(F, k1[j1]..k1[j1+1], k2[j2]..k2[j2+1], gl1, gl2)
         end
         S1 += S2
     end
     return S1
 end
 
-function _f_b_int_I(func, i1, i2, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, nip1, nip2, nodes1, nodes2, weights1, weights2)
+function _f_b_int_I(func, i1, i2, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, gl1, gl2)
     k1, k2 = knots(P1), knots(P2)
     m1, m2 = length(k1), length(k2)
     p1, p2 = degree(P1), degree(P2)
@@ -162,22 +162,22 @@ function _f_b_int_I(func, i1, i2, P1::AbstractBSplineSpace, P2::AbstractBSplineS
     f1, l1 = max(i1, p1+1), min(i1+p1, m1-p1-1)
     f2, l2 = max(i2, p2+1), min(i2+p2, m2-p2-1)
 
-    S2 = integrate(F, k1[f1]..k1[f1+1], k2[f2]..k2[f2+1], nip1, nip2, nodes1, nodes2, weights1, weights2)
+    S2 = integrate(F, k1[f1]..k1[f1+1], k2[f2]..k2[f2+1], gl1, gl2)
     for j2 in f2+1:l2
-        S2 += integrate(F, k1[f1]..k1[f1+1], k2[j2]..k2[j2+1], nip1, nip2, nodes1, nodes2, weights1, weights2)
+        S2 += integrate(F, k1[f1]..k1[f1+1], k2[j2]..k2[j2+1], gl1, gl2)
     end
     S1 = S2
     for j1 in f1+1:l1
-        S2 = integrate(F, k1[j1]..k1[j1+1], k2[f2]..k2[f2+1], nip1, nip2, nodes1, nodes2, weights1, weights2)
+        S2 = integrate(F, k1[j1]..k1[j1+1], k2[f2]..k2[f2+1], gl1, gl2)
         for j2 in f2+1:l2
-            S2 += integrate(F, k1[j1]..k1[j1+1], k2[j2]..k2[j2+1], nip1, nip2, nodes1, nodes2, weights1, weights2)
+            S2 += integrate(F, k1[j1]..k1[j1+1], k2[j2]..k2[j2+1], gl1, gl2)
         end
         S1 += S2
     end
     return S1
 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)
+function _f_b_int_R(func, i1, i2, i3, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, P3::AbstractBSplineSpace, gl1, gl2, gl3)
     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)
@@ -186,29 +186,29 @@ function _f_b_int_R(func, i1, i2, i3, P1::AbstractBSplineSpace, P2::AbstractBSpl
     f2, l2 = i2, i2+p2
     f3, l3 = i3, i3+p3
 
-    S3 = integrate(F, k1[f1]..k1[f1+1], k2[f2]..k2[f2+1], k3[f3]..k3[f3+1], nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3)
+    S3 = integrate(F, k1[f1]..k1[f1+1], k2[f2]..k2[f2+1], k3[f3]..k3[f3+1], gl1, gl2, gl3)
     for j3 in f3+1:l3
-        S3 += integrate(F, k1[f1]..k1[f1+1], k2[f2]..k2[f2+1], k3[j3]..k3[j3+1], nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3)
+        S3 += integrate(F, k1[f1]..k1[f1+1], k2[f2]..k2[f2+1], k3[j3]..k3[j3+1], gl1, gl2, gl3)
     end
     S2 = S3
     for j2 in f2+1:l2
-        S3 = integrate(F, k1[f1]..k1[f1+1], k2[j2]..k2[j2+1], k3[f3]..k3[f3+1], nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3)
+        S3 = integrate(F, k1[f1]..k1[f1+1], k2[j2]..k2[j2+1], k3[f3]..k3[f3+1], gl1, gl2, gl3)
         for j3 in f3+1:l3
-            S3 += integrate(F, k1[f1]..k1[f1+1], k2[j2]..k2[j2+1], k3[j3]..k3[j3+1], nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3)
+            S3 += integrate(F, k1[f1]..k1[f1+1], k2[j2]..k2[j2+1], k3[j3]..k3[j3+1], gl1, gl2, gl3)
         end
         S2 += S3
     end
     S1 = S2
     for j1 in f1+1:l1
-        S3 = integrate(F, k1[j1]..k1[j1+1], k2[f2]..k2[f2+1], k3[f3]..k3[f3+1], nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3)
+        S3 = integrate(F, k1[j1]..k1[j1+1], k2[f2]..k2[f2+1], k3[f3]..k3[f3+1], gl1, gl2, gl3)
         for j3 in f3+1:l3
-            S3 += integrate(F, k1[j1]..k1[j1+1], k2[f2]..k2[f2+1], k3[j3]..k3[j3+1], nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3)
+            S3 += integrate(F, k1[j1]..k1[j1+1], k2[f2]..k2[f2+1], k3[j3]..k3[j3+1], gl1, gl2, gl3)
         end
         S2 = S3
         for j2 in f2+1:l2
-            S3 = integrate(F, k1[j1]..k1[j1+1], k2[j2]..k2[j2+1], k3[f3]..k3[f3+1], nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3)
+            S3 = integrate(F, k1[j1]..k1[j1+1], k2[j2]..k2[j2+1], k3[f3]..k3[f3+1], gl1, gl2, gl3)
             for j3 in f3+1:l3
-                S3 += integrate(F, k1[j1]..k1[j1+1], k2[j2]..k2[j2+1], k3[j3]..k3[j3+1], nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3)
+                S3 += integrate(F, k1[j1]..k1[j1+1], k2[j2]..k2[j2+1], k3[j3]..k3[j3+1], gl1, gl2, gl3)
             end
             S2 += S3
         end
@@ -217,7 +217,7 @@ 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, nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3)
+function _f_b_int_I(func, i1, i2, i3, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace, P3::AbstractBSplineSpace, gl1, gl2, gl3)
     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)
@@ -227,29 +227,29 @@ function _f_b_int_I(func, i1, i2, i3, P1::AbstractBSplineSpace, P2::AbstractBSpl
     f2, l2 = max(i2, p2+1), min(i2+p2, m2-p2-1)
     f3, l3 = max(i3, p3+1), min(i3+p3, m3-p3-1)
 
-    S3 = integrate(F, k1[f1]..k1[f1+1], k2[f2]..k2[f2+1], k3[f3]..k3[f3+1], nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3)
+    S3 = integrate(F, k1[f1]..k1[f1+1], k2[f2]..k2[f2+1], k3[f3]..k3[f3+1], gl1, gl2, gl3)
     for j3 in f3+1:l3
-        S3 += integrate(F, k1[f1]..k1[f1+1], k2[f2]..k2[f2+1], k3[j3]..k3[j3+1], nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3)
+        S3 += integrate(F, k1[f1]..k1[f1+1], k2[f2]..k2[f2+1], k3[j3]..k3[j3+1], gl1, gl2, gl3)
     end
     S2 = S3
     for j2 in f2+1:l2
-        S3 = integrate(F, k1[f1]..k1[f1+1], k2[j2]..k2[j2+1], k3[f3]..k3[f3+1], nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3)
+        S3 = integrate(F, k1[f1]..k1[f1+1], k2[j2]..k2[j2+1], k3[f3]..k3[f3+1], gl1, gl2, gl3)
         for j3 in f3+1:l3
-            S3 += integrate(F, k1[f1]..k1[f1+1], k2[j2]..k2[j2+1], k3[j3]..k3[j3+1], nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3)
+            S3 += integrate(F, k1[f1]..k1[f1+1], k2[j2]..k2[j2+1], k3[j3]..k3[j3+1], gl1, gl2, gl3)
         end
         S2 += S3
     end
     S1 = S2
     for j1 in f1+1:l1
-        S3 = integrate(F, k1[j1]..k1[j1+1], k2[f2]..k2[f2+1], k3[f3]..k3[f3+1], nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3)
+        S3 = integrate(F, k1[j1]..k1[j1+1], k2[f2]..k2[f2+1], k3[f3]..k3[f3+1], gl1, gl2, gl3)
         for j3 in f3+1:l3
-            S3 += integrate(F, k1[j1]..k1[j1+1], k2[f2]..k2[f2+1], k3[j3]..k3[j3+1], nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3)
+            S3 += integrate(F, k1[j1]..k1[j1+1], k2[f2]..k2[f2+1], k3[j3]..k3[j3+1], gl1, gl2, gl3)
         end
         S2 = S3
         for j2 in f2+1:l2
-            S3 = integrate(F, k1[j1]..k1[j1+1], k2[j2]..k2[j2+1], k3[f3]..k3[f3+1], nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3)
+            S3 = integrate(F, k1[j1]..k1[j1+1], k2[j2]..k2[j2+1], k3[f3]..k3[f3+1], gl1, gl2, gl3)
             for j3 in f3+1:l3
-                S3 += integrate(F, k1[j1]..k1[j1+1], k2[j2]..k2[j2+1], k3[j3]..k3[j3+1], nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3)
+                S3 += integrate(F, k1[j1]..k1[j1+1], k2[j2]..k2[j2+1], k3[j3]..k3[j3+1], gl1, gl2, gl3)
             end
             S2 += S3
         end
@@ -262,8 +262,8 @@ function innerproduct_R(func, P1::AbstractBSplineSpace)
     p1 = degree(P1)
     n1 = dim(P1)
     nip1 = p1 + 1
-    nodes1, weights1 = gausslegendre(nip1)
-    b = [_f_b_int_R(func, i1, P1, nip1, nodes1, weights1) for i1 in 1:n1]
+    gl1 = GaussLegendre(nip1)
+    b = [_f_b_int_R(func, i1, P1, gl1) for i1 in 1:n1]
     return b
 end
 
@@ -271,8 +271,8 @@ function innerproduct_I(func, P1::AbstractBSplineSpace)
     p1 = degree(P1)
     n1 = dim(P1)
     nip1 = p1 + 1
-    nodes1, weights1 = gausslegendre(nip1)
-    b = [_f_b_int_I(func, i1, P1, nip1, nodes1, weights1) for i1 in 1:n1]
+    gl1 = GaussLegendre(nip1)
+    b = [_f_b_int_I(func, i1, P1, gl1) for i1 in 1:n1]
     return b
 end
 
@@ -281,9 +281,9 @@ function innerproduct_R(func, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace
     n1, n2 = dim(P1), dim(P2)
     nip1 = p1 + 1
     nip2 = p2 + 1
-    nodes1, weights1 = gausslegendre(nip1)
-    nodes2, weights2 = gausslegendre(nip2)
-    b = [_f_b_int_R(func, i1, i2, P1, P2, nip1, nip2, nodes1, nodes2, weights1, weights2) for i1 in 1:n1, i2 in 1:n2]
+    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]
     return b
 end
 
@@ -292,9 +292,9 @@ function innerproduct_I(func, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace
     n1, n2 = dim(P1), dim(P2)
     nip1 = p1 + 1
     nip2 = p2 + 1
-    nodes1, weights1 = gausslegendre(nip1)
-    nodes2, weights2 = gausslegendre(nip2)
-    b = [_f_b_int_I(func, i1, i2, P1, P2, nip1, nip2, nodes1, nodes2, weights1, weights2) for i1 in 1:n1, i2 in 1:n2]
+    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]
     return b
 end
 
@@ -304,10 +304,10 @@ function innerproduct_R(func, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace
     nip1 = p1 + 1
     nip2 = p2 + 1
     nip3 = p3 + 1
-    nodes1, weights1 = gausslegendre(nip1)
-    nodes2, weights2 = gausslegendre(nip2)
-    nodes3, weights3 = gausslegendre(nip3)
-    b = [_f_b_int_R(func, i1, i2, i3, P1, P2, P3, nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3) for i1 in 1:n1, i2 in 1:n2, i3 in 1:n3]
+    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]
     return b
 end
 
@@ -317,10 +317,10 @@ function innerproduct_I(func, P1::AbstractBSplineSpace, P2::AbstractBSplineSpace
     nip1 = p1 + 1
     nip2 = p2 + 1
     nip3 = p3 + 1
-    nodes1, weights1 = gausslegendre(nip1)
-    nodes2, weights2 = gausslegendre(nip2)
-    nodes3, weights3 = gausslegendre(nip3)
-    b = [_f_b_int_I(func, i1, i2, i3, P1, P2, P3, nip1, nip2, nip3, nodes1, nodes2, nodes3, weights1, weights2, weights3) for i1 in 1:n1, i2 in 1:n2, i3 in 1:n3]
+    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]
     return b
 end
 
diff --git a/src/_Geometry.jl b/src/_Geometry.jl
deleted file mode 100644
index 99222c7c..00000000
--- a/src/_Geometry.jl
+++ /dev/null
@@ -1,82 +0,0 @@
-# Geometry
-
-"""
-tangentvectors on B-spline manifold
-NOTICE: these methods are experimental
-"""
-tangentvectors
-
-function tangentvectors(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)
-    d̂ = size(a)[end]
-
-    S1 = Array{Float64}(undef, d̂)
-    f1 = j1-p1
-    l1 = j1
-    for xx in 1:d̂
-        S1[xx] = b1′[1]*a[f1,xx]
-        for i1 in f1+1:l1
-            S1[xx] += b1′[i1-f1+1]*a[i1,xx]
-        end
-    end
-
-    return S1
-end
-
-
-function tangentvectors(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)
-    b1′ = _bsplinebasis′(P1,t1,j1)
-    b2′ = _bsplinebasis′(P2,t2,j2)
-    d̂ = size(a)[end]
-
-    S1 = Array{Float64}(undef, d̂)
-    S2 = Array{Float64}(undef, d̂)
-    f1, f2 = j1-p1, j2-p2
-    l1, l2 = j1, j2
-
-    for xx in 1:d̂
-        S2[xx] = b2[1]*a[f1,f2,xx]
-        for i2 in f2+1:l2
-            S2[xx] += b2[i2-f2+1]*a[f1,i2,xx]
-        end
-        S1[xx] = b1′[1]*S2[xx]
-        for i1 in f1+1:l1
-            S2[xx] = b2[1]*a[i1,f2,xx]
-            for i2 in f2+1:l2
-                S2[xx] += b2[i2-f2+1]*a[i1,i2,xx]
-            end
-            S1[xx] += b1′[i1-f1+1]*S2[xx]
-        end
-    end
-    X1 = copy(S1)
-
-    for xx in 1:d̂
-        S2[xx] = b2′[1]*a[f1,f2,xx]
-        for i2 in f2+1:l2
-            S2[xx] += b2′[i2-f2+1]*a[f1,i2,xx]
-        end
-        S1[xx] = b1[1]*S2[xx]
-        for i1 in f1+1:l1
-            S2[xx] = b2′[1]*a[i1,f2,xx]
-            for i2 in f2+1:l2
-                S2[xx] += b2′[i2-f2+1]*a[i1,i2,xx]
-            end
-            S1[xx] += b1[i1-f1+1]*S2[xx]
-        end
-    end
-    X2 = copy(S1)
-
-    return (X1,X2)
-end
diff --git a/src/_Integral.jl b/src/_Integral.jl
index 444ebb5e..f3b407ef 100644
--- a/src/_Integral.jl
+++ b/src/_Integral.jl
@@ -1,9 +1,20 @@
 # Numerical Integration
 
+struct GaussLegendre
+    nip::Int  # number of integral points
+    nodes::Vector{Float64}
+    weights::Vector{Float64}
+    function GaussLegendre(nip::Int)
+        nodes, weights = gausslegendre(nip)
+        return new(nip, nodes, weights)
+    end
+end
+
 """
 Integrate on interval
 """
-function integrate(func, I1::ClosedInterval{<:Real}, nip1, nodes1, weights1)
+function integrate(func, I1::ClosedInterval{<:Real}, gl1)
+    nip1, nodes1, weights1 = gl1.nip, gl1.nodes, gl1.weights
     dnodes1 = (width(I1) * nodes1 .+ sum(extrema(I1))) / 2
 
     S1 = weights1[1] * func(dnodes1[1])
@@ -17,7 +28,9 @@ end
 """
 Integrate on rectangular region.
 """
-function integrate(func, I1::ClosedInterval{<:Real}, I2::ClosedInterval{<:Real}, nip1, nip2, nodes1, nodes2, weights1, weights2)
+function integrate(func, I1::ClosedInterval{<:Real}, I2::ClosedInterval{<:Real}, gl1, gl2)
+    nip1, nodes1, weights1 = gl1.nip, gl1.nodes, gl1.weights
+    nip2, nodes2, weights2 = gl2.nip, gl2.nodes, gl2.weights
     dnodes1 = (width(I1) * nodes1 .+ sum(extrema(I1))) / 2
     dnodes2 = (width(I2) * nodes2 .+ sum(extrema(I2))) / 2
 
@@ -40,7 +53,10 @@ 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)
+function integrate(func, I1::ClosedInterval{<:Real}, I2::ClosedInterval{<:Real}, I3::ClosedInterval{<:Real}, gl1, gl2, gl3)
+    nip1, nodes1, weights1 = gl1.nip, gl1.nodes, gl1.weights
+    nip2, nodes2, weights2 = gl2.nip, gl2.nodes, gl2.weights
+    nip3, nodes3, weights3 = gl3.nip, gl3.nodes, gl3.weights
     dnodes1 = (width(I1) * nodes1 .+ sum(extrema(I1))) / 2
     dnodes2 = (width(I2) * nodes2 .+ sum(extrema(I2))) / 2
     dnodes3 = (width(I3) * nodes3 .+ sum(extrema(I3))) / 2
diff --git a/src/_Knots.jl b/src/_Knots.jl
index 9d14aed2..7652a9d3 100644
--- a/src/_Knots.jl
+++ b/src/_Knots.jl
@@ -131,7 +131,7 @@ julia> length(k)
 Base.length(k::Knots) = length(k.vector)
 
 function ♯(k::Knots)
-    @warn "BasicBSpline.♯ is deprecated."
+    @warn "BasicBSpline.♯ is deprecated. Use `length` instead."
     return length(k::Knots)
 end
 Base.firstindex(k::Knots) = 1
@@ -160,6 +160,9 @@ Knots([1.0, 2.0, 3.0])
 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)
+Base.searchsortedfirst(k::Knots,t) = searchsortedfirst(k.vector,t)
+Base.searchsortedlast(k::Knots,t) = searchsortedlast(k.vector,t)
+Base.searchsorted(k::Knots,t) = searchsorted(k.vector,t)
 
 @doc raw"""
 Check a inclusive relation ship ``k\subset k'``.
@@ -169,13 +172,14 @@ Check a inclusive relation ship ``k\subset k'``.
 ```
 """
 function Base.issubset(k::Knots, k′::Knots)
-    K′ = copy(k′.vector)
-    for kᵢ in k.vector
-        i = findfirst(==(kᵢ), K′)
+    v = k′.vector
+    l = length(v)
+    i = 0
+    for kᵢ in k
+        i = findnext(==(kᵢ), v, i+1)
         if isnothing(i)
             return false
         end
-        deleteat!(K′, i)
     end
     return true
 end
@@ -201,25 +205,36 @@ julia> 𝔫(k,2.0)
 2
 ```
 """
-𝔫(k::Knots, t::Real) = count(==(t), k.vector)
+function 𝔫(k::Knots, t::Real)
+    # for small case, this is faster
+    # return count(==(t), k.vector)
+
+    # for large case, this is faster
+    return length(searchsorted(k,t))
+end
 
 """
 Find an index ``i`` such that ``k_{i} ≤ t < k_{i+1}``.
 """
 function _knotindex₊₀(k::Union{Knots, AbstractVector{<:Real}}, t::Real)
-    return findfirst(i -> k[i]≤t<k[i+1], 1:length(k)-1)
+    return searchsortedlast(k, t)
 end
 
 """
 Find an index ``i`` such that ``k_{i} < t ≤ k_{i+1}``.
 """
 function _knotindex₋₀(k::Union{Knots, AbstractVector{<:Real}}, t::Real)
-    return findfirst(i -> k[i]<t≤k[i+1], 1:length(k)-1)
+    return searchsortedfirst(k, t) - 1
 end
 
 """
 Find an index ``i`` such that ``k_{i} ≤ t < k_{i+1}`` or ``k_{i} < t = k_{i+1} = k_{\text{end}}``.
 """
 function _knotindex(k::Union{Knots, AbstractVector{<:Real}}, t::Real)
-    return findfirst(i -> (k[i]≤t<k[i+1])|(k[i]<t==k[i+1]==k[end]), 1:length(k)-1)
+    j = searchsortedlast(k, t)
+    if j < length(k)
+        return j
+    else
+        return searchsortedfirst(k, t) - 1
+    end
 end


Test patch:
diff --git a/test/runtests.jl b/test/runtests.jl
index c55cc76c..70b8b7f8 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -10,5 +10,6 @@ using Random
     include("test_FastBSplineSpace.jl")
     include("test_BSplineBasis.jl")
     include("test_BSplineManifold.jl")
+    include("test_Integral.jl")
     include("test_Fitting.jl")
 end
diff --git a/test/test_Integral.jl b/test/test_Integral.jl
new file mode 100644
index 00000000..aa1ebd20
--- /dev/null
+++ b/test/test_Integral.jl
@@ -0,0 +1,21 @@
+@testset begin
+    gl1 = BasicBSpline.GaussLegendre(1)
+    gl2 = BasicBSpline.GaussLegendre(2)
+    gl3 = BasicBSpline.GaussLegendre(3)
+    gl4 = BasicBSpline.GaussLegendre(4)
+
+    I = 2..8
+    coef = rand(10)
+
+    f(n,x) = sum(coef[i]*x^(i-1) for i in 1:n)
+    F(n,x) = sum(coef[i]*x^(i)/i for i in 1:n)
+
+    @test F(1,maximum(I)) - F(1,minimum(I)) ≈ BasicBSpline.integrate(x->f(1,x), I, gl1)
+    @test F(2,maximum(I)) - F(2,minimum(I)) ≈ BasicBSpline.integrate(x->f(2,x), I, gl1)
+    @test F(3,maximum(I)) - F(3,minimum(I)) ≈ BasicBSpline.integrate(x->f(3,x), I, gl2)
+    @test F(4,maximum(I)) - F(4,minimum(I)) ≈ BasicBSpline.integrate(x->f(4,x), I, gl2)
+    @test F(5,maximum(I)) - F(5,minimum(I)) ≈ BasicBSpline.integrate(x->f(5,x), I, gl3)
+    @test F(6,maximum(I)) - F(6,minimum(I)) ≈ BasicBSpline.integrate(x->f(6,x), I, gl3)
+    @test F(7,maximum(I)) - F(7,minimum(I)) ≈ BasicBSpline.integrate(x->f(7,x), I, gl4)
+    @test F(8,maximum(I)) - F(8,minimum(I)) ≈ BasicBSpline.integrate(x->f(8,x), I, gl4)
+end
diff --git a/test/test_Knots.jl b/test/test_Knots.jl
index 0829dce9..436781f5 100644
--- a/test/test_Knots.jl
+++ b/test/test_Knots.jl
@@ -1,46 +1,62 @@
+k = Knots([1,2,2,3])
+
 @testset "Knots" begin
     @test Knots([1,2,3]) == Knots(1:3) == Knots(1,2,3)
     @test Knots(2,4,5) == Knots([2,4,5])
     @test zero(Knots) == Knots() == Knots([])
     @test Knots(1:3) == Knots([3,2,1])
 
-    k = Knots([1,2,2,3])
     @test ♯(k) == length(k) == 4
 
     @test Knots([-1,2,3]) + 2 * Knots([2,5]) == Knots([-1,2,2,2,3,5,5])
     @test Knots([1,2,3]) + Knots([2,4,5]) == Knots([1,2,2,3,4,5])
     @test 2 * Knots([2,3]) == Knots([2,2,3,3])
 
-    @test unique(Knots([1,2,2,3])) == Knots([1,2,3])
+    @test unique(k) == Knots([1,2,3])
 
-    @test Knots([1,2,3]) ⊆ Knots([1,2,3])
-    @test Knots([1,2,3]) ⊇ Knots([1,2,3])
-    @test Knots([1,2,2,3]) ⊆ Knots([1,2,2,3,5])
-    @test Knots([1,2,2,3]) ⊈ Knots([1,2,3,5])
-    @test Knots([1,2,2,3,5]) ⊇ Knots([1,2,2,3])
-    @test Knots([1,2,3,5]) ⊉ Knots([1,2,2,3])
+    @testset "inclusive relation" begin
+        @test Knots([1,2,3])     ⊆ Knots([1,2,3])
+        @test Knots([1,2,3])     ⊇ Knots([1,2,3])
+        @test Knots([1,2,2,3])   ⊆ Knots([1,2,2,3,5])
+        @test Knots([1,2,2,3])   ⊈ Knots([1,2,3,5])
+        @test Knots([1,2,2,3,5]) ⊇ Knots([1,2,2,3])
+        @test Knots([1,2,3,5])   ⊉ Knots([1,2,2,3])
+    end
 
     @test 𝔫(k, 0.3) == 0
     @test 𝔫(k, 1.0) == 1
     @test 𝔫(k, 2.0) == 2
     @test 1 ∈ k
     @test 1.5 ∉ k
-    @test BasicBSpline._knotindex(k,1.5) == 1
-    @test BasicBSpline._knotindex(k,2.5) == 3
-    @test BasicBSpline._knotindex₋₀(k,1.5) == 1
-    @test BasicBSpline._knotindex₋₀(k,2.5) == 3
-    @test BasicBSpline._knotindex₊₀(k,1.5) == 1
-    @test BasicBSpline._knotindex₊₀(k,2.5) == 3
-
-    @test BasicBSpline._knotindex(k,1) == 1
-    @test BasicBSpline._knotindex(k,2) == 3
-    @test BasicBSpline._knotindex(k,3) == 3
-    @test BasicBSpline._knotindex₋₀(k,1) == nothing
-    @test BasicBSpline._knotindex₋₀(k,2) == 1
-    @test BasicBSpline._knotindex₋₀(k,3) == 3
-    @test BasicBSpline._knotindex₊₀(k,1) == 1
-    @test BasicBSpline._knotindex₊₀(k,2) == 3
-    @test BasicBSpline._knotindex₊₀(k,3) == nothing
+
+    @testset "knotindex for usual case" begin
+        @test BasicBSpline._knotindex(k,1.5) == 1
+        @test BasicBSpline._knotindex(k,2.5) == 3
+        @test BasicBSpline._knotindex₋₀(k,1.5) == 1
+        @test BasicBSpline._knotindex₋₀(k,2.5) == 3
+        @test BasicBSpline._knotindex₊₀(k,1.5) == 1
+        @test BasicBSpline._knotindex₊₀(k,2.5) == 3
+    end
+
+    @testset "knotindex for corner case" begin
+        @test BasicBSpline._knotindex(k,0) == 0
+        @test BasicBSpline._knotindex(k,1) == 1
+        @test BasicBSpline._knotindex(k,2) == 3
+        @test BasicBSpline._knotindex(k,3) == 3
+        @test BasicBSpline._knotindex(k,4) == 4
+
+        @test BasicBSpline._knotindex₋₀(k,0) == 0
+        @test BasicBSpline._knotindex₋₀(k,1) == 0
+        @test BasicBSpline._knotindex₋₀(k,2) == 1
+        @test BasicBSpline._knotindex₋₀(k,3) == 3
+        @test BasicBSpline._knotindex₋₀(k,4) == 4
+
+        @test BasicBSpline._knotindex₊₀(k,0) == 0
+        @test BasicBSpline._knotindex₊₀(k,1) == 1
+        @test BasicBSpline._knotindex₊₀(k,2) == 3
+        @test BasicBSpline._knotindex₊₀(k,3) == 4
+        @test BasicBSpline._knotindex₊₀(k,4) == 4
+    end
 
     @test string(k) == "Knots([1.0, 2.0, 2.0, 3.0])"
     @test string(Knots()) == "Knots([])"