From 09df10911cb1ec0e671072c39c646b2d6c020d5d Mon Sep 17 00:00:00 2001
From: Michael Reed <18372368+chakravala@users.noreply.github.com>
Date: Wed, 27 Nov 2024 17:52:23 -0500
Subject: [PATCH] refined index methods

---
 Project.toml    |   2 +-
 src/Cartan.jl   | 290 ++++++++++++++++++++++++++++++------------------
 src/diffgeo.jl  |  78 +++++++++----
 src/topology.jl | 229 +++++++++++++++++++++++++-------------
 4 files changed, 394 insertions(+), 205 deletions(-)

diff --git a/Project.toml b/Project.toml
index 20ea236..7ad2b96 100644
--- a/Project.toml
+++ b/Project.toml
@@ -12,7 +12,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 [compat]
 julia = "1"
 Requires = "1"
-Grassmann = "0.8, 0.9"
+Grassmann = "0.8"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/src/Cartan.jl b/src/Cartan.jl
index dc243d3..c90fed5 100644
--- a/src/Cartan.jl
+++ b/src/Cartan.jl
@@ -30,6 +30,8 @@ import Grassmann: realvalue, imagvalue, points, metrictensor, metricextensor
 import Grassmann: Values, Variables, FixedVector, list
 import Grassmann: Scalar, GradedVector, Bivector, Trivector
 import Base: @pure, OneTo, getindex
+import LinearAlgebra: cross
+#import ElasticArrays: resize_lastdim!
 
 export Values, Derivation
 export initmesh, pdegrad, det
@@ -57,22 +59,26 @@ end
 
 struct TensorField{B,F,N,M<:AbstractFrameBundle{B,N}} <: GlobalFiber{LocalTensor{B,F},N}
     dom::M
-    cod::Array{F,N}
+    cod::Array{F,N}#ElasticArray{F,N,L,Vector{F}}
+    #function TensorField(dom::M,cod::ElasticArray{F,N,L,Vector{F}}) where {B,F,N,M<:AbstractFrameBundle{B,N},L}
     function TensorField(dom::M,cod::Array{F,N}) where {B,F,N,M<:AbstractFrameBundle{B,N}}
         new{B,F,N,M}(dom,cod)
     end
 end
 
-function TensorField(id::Int,dom::PA,cod::Array{F,N},met::GA=Global{N}(InducedMetric())) where {N,P,F,PA<:AbstractArray{P,N},GA<:AbstractArray}
+#function TensorField(dom::M,cod::Array{F,N}) where {B,F,N,M<:AbstractFrameBundle{B,N}}
+#    TensorField(dom,ElasticArray(cod))
+#end
+function TensorField(id::Int,dom::PA,cod::DenseArray{F,N},met::GA=Global{N}(InducedMetric())) where {N,P,F,PA<:AbstractArray{P,N},GA<:AbstractArray}
     TensorField(GridFrameBundle(id,dom,met),cod)
 end
-function TensorField(id::Int,dom::P,cod::Vector{F},met::G=Global{N}(InducedMetric())) where {F,P<:PointCloud,G<:AbstractVector}
+function TensorField(id::Int,dom::P,cod::DenseVector{F},met::G=Global{N}(InducedMetric())) where {F,P<:PointCloud,G<:AbstractVector}
     TensorField(SimplexFrameBundle(id,dom,met),cod)
 end
-TensorField(id::Int,dom,cod::Array,met::GlobalFiber) = TensorField(id,dom,cod,fiber(met))
+TensorField(id::Int,dom,cod::DenseArray,met::GlobalFiber) = TensorField(id,dom,cod,fiber(met))
 TensorField(dom::AbstractFrameBundle,cod::AbstractFrameBundle) = TensorField(dom,points(cod))
-TensorField(dom::AbstractArray{B,N} where B,cod::Array{F,N} where F,met::AbstractArray=Global{N}(InducedMetric())) where N = TensorField((global grid_id+=1),dom,cod,fiber(met))
-TensorField(dom::ChainBundle,cod::Vector,met::AbstractVector=Global{1}(InducedMetric())) = TensorField((global grid_id+=1),dom,cod,met)
+TensorField(dom::AbstractArray{B,N} where B,cod::DenseArray{F,N} where F,met::AbstractArray=Global{N}(InducedMetric())) where N = TensorField((global grid_id+=1),dom,cod,fiber(met))
+TensorField(dom::ChainBundle,cod::DenseVector,met::AbstractVector=Global{1}(InducedMetric())) = TensorField((global grid_id+=1),dom,cod,met)
 TensorField(a::TensorField,b::TensorField) = TensorField(fiber(a),fiber(b))
 
 #const ParametricMesh{B,F,P<:AbstractVector{<:Chain}} = TensorField{B,F,1,P}
@@ -127,7 +133,7 @@ for bundle ∈ (:TensorField,:GlobalSection)
         $bundle(dom,fun::BitArray) = $bundle(dom, Float64.(fun))
         $bundle(dom,fun::$bundle) = $bundle(dom, fiber(fun))
         $bundle(dom::$bundle,fun) = $bundle(base(dom), fun)
-        $bundle(dom::$bundle,fun::Array) = $bundle(base(dom), fun)
+        $bundle(dom::$bundle,fun::DenseArray) = $bundle(base(dom), fun)
         $bundle(dom::$bundle,fun::Function) = $bundle(base(dom), fun)
         $bundle(dom::AbstractArray,fun::AbstractRange) = $bundle(dom, collect(fun))
         $bundle(dom::AbstractArray,fun::RealRegion) = $bundle(dom, collect(fun))
@@ -164,11 +170,18 @@ iscompact(t::TensorField) = iscompact(immersion(t))
 PointCloud(t::TensorField) = PointCloud(base(t))
 Grassmann.grade(::GradedField{G}) where G = G
 
+resize(t::TensorField) = TensorField(resize(domain(t)),codomain(t))
+resize(t::GlobalSection) = GlobalSection(resize(domain(t)),codomain(t))
+
 @pure Base.eltype(::Type{<:TensorField{B,F}}) where {B,F} = LocalTensor{B,F}
 Base.getindex(m::TensorField,i::Vararg{Int}) = LocalTensor(getindex(domain(m),i...), getindex(codomain(m),i...))
 Base.getindex(m::TensorField,i::Vararg{Union{Int,Colon}}) = TensorField(domain(m)(i...), getindex(codomain(m),i...))
 #Base.setindex!(m::TensorField{B,F,1,<:Interval},s::LocalTensor,i::Vararg{Int}) where {B,F} = setindex!(codomain(m),fiber(s),i...)
 Base.setindex!(m::TensorField{B,Fm} where Fm,s::F,i::Vararg{Int}) where {B,F} = setindex!(codomain(m),s,i...)
+Base.setindex!(m::TensorField{B,Fm} where Fm,s::TensorField,::Colon,i::Int) where B = setindex!(codomain(m),codomain(s),:,i)
+Base.setindex!(m::TensorField{B,Fm} where Fm,s::TensorField,::Colon,::Colon,i::Int) where B = setindex!(codomain(m),codomain(s),:,:,i)
+Base.setindex!(m::TensorField{B,Fm} where Fm,s::TensorField,::Colon,::Colon,::Colon,i::Int) where B = setindex!(codomain(m),codomain(s),:,:,:,i)
+Base.setindex!(m::TensorField{B,Fm} where Fm,s::TensorField,::Colon,::Colon,::Colon,::Colon,i::Int) where B = setindex!(codomain(m),codomain(s),:,:,:,:,i)
 function Base.setindex!(m::TensorField{B,F,N,<:IntervalRange} where {B,F,N},s::LocalTensor,i::Vararg{Int})
     setindex!(codomain(m),fiber(s),i...)
     return s
@@ -191,6 +204,18 @@ function Base.setindex!(m::GlobalSection,s::LocalSection,i::Vararg{Int})
     return s
 end
 
+extract(x::AbstractVector,i) = (@inbounds x[i])
+extract(x::AbstractMatrix,i) = (@inbounds LocalTensor(points(x).v[end][i],x[:,i]))
+extract(x::AbstractArray{T,3} where T,i) = (@inbounds LocalTensor(points(x).v[end][i],x[:,:,i]))
+extract(x::AbstractArray{T,4} where T,i) = (@inbounds LocalTensor(points(x).v[end][i],x[:,:,:,i]))
+extract(x::AbstractArray{T,5} where T,i) = (@inbounds LocalTensor(points(x).v[end][i],x[:,:,:,:,i]))
+
+assign!(x::AbstractVector,i,s) = (@inbounds x[i] = s)
+assign!(x::AbstractMatrix,i,s) = (@inbounds x[:,i] = s)
+assign!(x::AbstractArray{T,3} where T,i,s) = (@inbounds x[:,:,i] = s)
+assign!(x::AbstractArray{T,4} where T,i,s) = (@inbounds x[:,:,:,i] = s)
+assign!(x::AbstractArray{T,5} where T,i,s) = (@inbounds x[:,:,:,:,i] = s)
+
 Base.BroadcastStyle(::Type{<:TensorField{B,F,N,P}}) where {B,F,N,P} = Broadcast.ArrayStyle{TensorField{B,F,N,P}}()
 Base.BroadcastStyle(::Type{<:GlobalSection{B,F,N,BA,FA}}) where {B,F,N,BA,FA} = Broadcast.ArrayStyle{TensorField{B,F,N,BA,FA}}()
 
@@ -258,22 +283,32 @@ function quintlinterp(x,y,z,w,v,x1,x2,y1,y2,z1,z2,w1,w2,v1,v2,f11111,f21111,f121
     linterp(v,v1,v2,f1,f2)
 end
 
-reposition_odd(p,x,t) = iseven(p) ? x[end]-x[1]+t : 2x[1]-t
-reposition_even(p,x,t) = isodd(p) ? x[1]-x[end]+t : 2x[end]-t
+reposition_odd(p,x,t) = @inbounds (iseven(p) ? x[end]-x[1]+t : 2x[1]-t)
+reposition_even(p,x,t) = @inbounds (isodd(p) ? x[1]-x[end]+t : 2x[end]-t)
 @inline reposition(i1,i2,p1,p2,x,t) = i1 ? reposition_odd(p1,x,t) : i2 ? reposition_even(p2,x,t) : eltype(x)(t)
 
 (m::TensorField)(s::LocalTensor) = LocalTensor(base(s), m(fiber(s)))
-function (m::IntervalMap)(t); p = points(m)
-    i = searchsortedfirst(p,t[1])-1
+function (m::IntervalMap)(t); p,t1 = points(m),(@inbounds t[1])
+    i = searchsortedfirst(p,@inbounds t1)-1
     if !isopen(m)
         q = immersion(m)
-        if iszero(i) && !iszero(q.r[1])
-            return m(reposition_odd(q.p[1],p,t[1]))
-        elseif i==length(p) && !iszero(q.r[2])
-            return m(reposition_even(q.p[2],p,t[1]))
+        if iszero(i)
+            if iszero(@inbounds q.r[1])
+                return zero(fibertype(m))
+            else
+                return m(@inbounds reposition_odd(q.p[1],p,t1))
+            end
+        elseif i==length(p)
+            if iszero(@inbounds q.r[2])
+                return zero(fibertype(m))
+            else
+                return m(@inbounds reposition_even(q.p[2],p,t1))
+            end
         end
+    elseif iszero(i) || i==length(p)
+        return zero(fibertype(m))
     end
-    linterp(t[1],p[i],p[i+1],m.cod[i],m.cod[i+1])
+    linterp(t1,p[i],p[i+1],m.cod[i],m.cod[i+1])
 end
 function (m::IntervalMap)(t::Vector,d=diff(m.cod)./diff(m.dom))
     [parametric(i,m,d) for i ∈ t]
@@ -284,22 +319,24 @@ function parametric(t,m,d=diff(codomain(m))./diff(domain(m)))
     codomain(m)[i]+(t-p[i])*d[i]
 end
 
-(m::RectangleMap)(i::Int,j::Int=2) = isone(j) ? m[i,:] : m[:,i]
-function (m::RectangleMap)(t::AbstractFloat,j::Int=2)
+(m::RectangleMap)(t::Real) = slice(m,t)
+slice(m::RectangleMap,i::Int,j::Int=2) = isone(j) ? m[i,:] : m[:,i]
+function slice(m::RectangleMap,t::AbstractFloat,j::Int=2)
     Q,p = isone(j),points(m).v[j]
     i = searchsortedfirst(p,t)-1
     TensorField(points(m).v[Q ? 2 : 1],linterp(t,p[i],p[i+1],Q ? m.cod[i,:] : m.cod[:,i],Q ? m.cod[i+1,:] : m.cod[:,i+1]))
 end
 
-function slice(m::HyperrectangleMap,i::Int,j::Int,k::Int=3)
+(m::HyperrectangleMap)(t::Real,j::Int=3) = slice(m,t,j)
+function slice2(m::HyperrectangleMap,i::Int,j::Int,k::Int=3)
     isone(k) ? m[:,i,j] : k==2 ? m[i,:,j] : m[i,j,:]
 end
-function (m::HyperrectangleMap)(i::Int,j::Int=3)
+function slice(m::HyperrectangleMap,i::Int,j::Int=3)
     isone(j) ? m[i,:,:] : j==2 ? m[:,i,:] : m[:,:,i]
 end
 
-(m::TensorField{B,F,2,<:FiberProductBundle} where {B,F})(i::Int) = TensorField(base(base(m)),fiber(m)[:,i])
-function (m::TensorField{B,F,2,<:FiberProductBundle} where {B,F})(t::AbstractFloat)
+slice(m::TensorField{B,F,2,<:FiberProductBundle} where {B,F},i::Int) = TensorField(base(base(m)),fiber(m)[:,i])
+function (m::TensorField{B,F,2,<:FiberProductBundle} where {B,F})(t::Real)
     k = 2; p = base(m).cod.v[1]
     i = searchsortedfirst(p,t)-1
     TensorField(base(base(m)),linterp(t,p[i],p[i+1],m.cod[:,i],m.cod[:,i+1]))
@@ -315,46 +352,65 @@ function (m::TensorField{B,F,N,<:SimplexFrameBundle} where {B,N})(t) where F
 end
 
 (m::TensorField{B,F,N,<:RealSpace{2}} where {B,F,N})(x,y) = m(Chain(x,y))
-function (m::TensorField{B,F,N,<:RealSpace{2}} where {B,F,N})(t::Chain{V}) where V
-    x,y = points(m).v[1],points(m).v[2]
-    i,j = searchsortedfirst(x,t[1])-1,searchsortedfirst(y,t[2])-1
+function (m::TensorField{B,F,N,<:RealSpace{2}} where {B,N})(t::Chain{V,G,T,2} where {G,T}) where {F,V}
+    x,y,t1,t2 = @inbounds (points(m).v[1],points(m).v[2],t[1],t[2])
+    i,j = searchsortedfirst(x,t1)-1,searchsortedfirst(y,t2)-1
     if !isopen(m)
         q = immersion(m)
-        i1,i2,j1,j2 = (
-            iszero(i) && !iszero(q.r[1]),i==length(x) && !iszero(q.r[2]),
-            iszero(j) && !iszero(q.r[3]),j==length(y) && !iszero(q.r[4]))
-        if i1 || i2 || j1 || j2
-            return m(Chain{V}(
-                reposition(i1,i2,q.p[1],q.p[2],x,t[1]),
-                reposition(j1,j2,q.p[3],q.p[4],y,t[2])))
+        i01,i02,iq1,iq2,j01,j02,jq1,jq2 = @inbounds (
+            iszero(i),i==length(x),iszero(q.r[1]),iszero(q.r[2]),
+            iszero(j),j==length(y),iszero(q.r[3]),iszero(q.r[4]))
+        if (i01 && iq1) || (i02 && iq2) || (j01 && jq1) || (j02 && jq2)
+            return zero(F)
+        else
+            i1,i2,j1,j2 = (
+                i01 && !iq1,i02 && !iq2,
+                j01 && !jq1,j02 && !jq2)
+            if i1 || i2 || j1 || j2
+                return m(Chain{V}(
+                    (@inbounds reposition(i1,i2,q.p[1],q.p[2],x,t1)),
+                    (@inbounds reposition(j1,j2,q.p[3],q.p[4],y,t2))))
+            end
         end
+    elseif iszero(i) || iszero(j) || i==length(x) || j==length(y)
+        return zero(F) # elseif condition creates 1 allocation, as opposed to 0 ???
     end
     #f1 = linterp(t[1],x[i],x[i+1],m.cod[i,j],m.cod[i+1,j])
     #f2 = linterp(t[1],x[i],x[i+1],m.cod[i,j+1],m.cod[i+1,j+1])
     #linterp(t[2],y[j],y[j+1],f1,f2)
-    bilinterp(t[1],t[2],x[i],x[i+1],y[j],y[j+1],
+    bilinterp(t1,t2,x[i],x[i+1],y[j],y[j+1],
         m.cod[i,j],m.cod[i+1,j],m.cod[i,j+1],m.cod[i+1,j+1])
 end
 
 (m::TensorField{B,F,N,<:RealSpace{3}} where {B,F,N})(x,y,z) = m(Chain(x,y,z))
-function (m::TensorField{B,F,N,<:RealSpace{3}} where {B,F,N})(t::Chain{V}) where V
-    x,y,z = points(m).v[1],points(m).v[2],points(m).v[3]
+function (m::TensorField{B,F,N,<:RealSpace{3}} where {B,N})(t::Chain{V,G,T,3} where {G,T}) where {F,V}
+    x,y,z,t1,t2,t3 = @inbounds (points(m).v[1],points(m).v[2],points(m).v[3],t[1],t[2],t[3])
     i,j,k = (
-        searchsortedfirst(x,t[1])-1,
-        searchsortedfirst(y,t[2])-1,
-        searchsortedfirst(z,t[3])-1)
+        searchsortedfirst(x,t1)-1,
+        searchsortedfirst(y,t2)-1,
+        searchsortedfirst(z,t3)-1)
     if !isopen(m)
         q = immersion(m)
-        i1,i2,j1,j2,k1,k2 = (
-            iszero(i) && !iszero(q.r[1]),i==length(x) && !iszero(q.r[2]),
-            iszero(j) && !iszero(q.r[3]),j==length(y) && !iszero(q.r[4]),
-            iszero(k) && !iszero(q.r[5]),k==length(z) && !iszero(q.r[6]))
-        if i1 || i2 || j1 || j2 || k1 || k2
-            return m(Chain{V}(
-                reposition(i1,i2,q.p[1],q.p[2],x,t[1]),
-                reposition(j1,j2,q.p[3],q.p[4],y,t[2]),
-                reposition(k1,k2,q.p[5],q.p[6],z,t[3])))
+        i01,i02,iq1,iq2,j01,j02,jq1,jq2,k01,k02,kq1,kq2 = @inbounds (
+            iszero(i),i==length(x),iszero(q.r[1]),iszero(q.r[2]),
+            iszero(j),j==length(y),iszero(q.r[3]),iszero(q.r[4]),
+            iszero(k),k==length(z),iszero(q.r[5]),iszero(q.r[6]))
+        if (i01 && iq1) || (i02 && iq2) || (j01 && jq1) || (j02 && jq2) || (k01 && kq1) || (k02 && kq2)
+            return zero(F)
+        else
+            i1,i2,j1,j2,k1,k2 = (
+                i01 && !iq1,i02 && !iq2,
+                j01 && !jq1,j02 && !jq2,
+                k01 && !kq1,k02 && !kq2)
+            if i1 || i2 || j1 || j2 || k1 || k2
+                return m(Chain{V}(
+                    (@inbounds reposition(i1,i2,q.p[1],q.p[2],x,t1)),
+                    (@inbounds reposition(j1,j2,q.p[3],q.p[4],y,t2)),
+                    (@inbounds reposition(k1,k2,q.p[5],q.p[6],z,t3))))
+            end
         end
+    elseif iszero(i) || iszero(j) || iszero(k) || i==length(x) || j==length(y) || k==length(z)
+        return zero(F) # elseif condition creates 1 allocation, as opposed to 0 ???
     end
     #f1 = linterp(t[1],x[i],x[i+1],m.cod[i,j,k],m.cod[i+1,j,k])
     #f2 = linterp(t[1],x[i],x[i+1],m.cod[i,j+1,k],m.cod[i+1,j+1,k])
@@ -363,33 +419,44 @@ function (m::TensorField{B,F,N,<:RealSpace{3}} where {B,F,N})(t::Chain{V}) where
     #f4 = linterp(t[1],x[i],x[i+1],m.cod[i,j+1,k+1],m.cod[i+1,j+1,k+1])
     #g2 = linterp(t[2],y[j],y[j+1],f3,f4)
     #linterp(t[3],z[k],z[k+1],g1,g2)
-    trilinterp(t[1],t[2],t[3],x[i],x[i+1],y[j],y[j+1],z[k],z[k+1],
+    trilinterp(t1,t2,t3,x[i],x[i+1],y[j],y[j+1],z[k],z[k+1],
         m.cod[i,j,k],m.cod[i+1,j,k],m.cod[i,j+1,k],m.cod[i+1,j+1,k],
         m.cod[i,j,k+1],m.cod[i+1,j,k+1],m.cod[i,j+1,k+1],m.cod[i+1,j+1,k+1])
 end
 
 (m::TensorField{B,F,N,<:RealSpace{4}} where {B,F,N})(x,y,z,w) = m(Chain(x,y,z,w))
-function (m::TensorField{B,F,N,<:RealSpace{4}} where {B,F,N})(t::Chain{V}) where V
-    x,y,z,w = points(m).v[1],points(m).v[2],points(m).v[3],points(m).v[4]
+function (m::TensorField{B,F,N,<:RealSpace{4}} where {B,N})(t::Chain{V,G,T,4} where {G,T}) where {F,V}
+    x,y,z,w,t1,t2,t3,t4 = @inbounds (points(m).v[1],points(m).v[2],points(m).v[3],points(m).v[4],t[1],t[2],t[3],t[4])
     i,j,k,l = (
-        searchsortedfirst(x,t[1])-1,
-        searchsortedfirst(y,t[2])-1,
-        searchsortedfirst(z,t[3])-1,
-        searchsortedfirst(w,t[4])-1)
+        searchsortedfirst(x,t1)-1,
+        searchsortedfirst(y,t2)-1,
+        searchsortedfirst(z,t3)-1,
+        searchsortedfirst(w,t4)-1)
     if !isopen(m)
         q = immersion(m)
-        i1,i2,j1,j2,k1,k2,l1,l2 = (
-            iszero(i) && !iszero(q.r[1]),i==length(x) && !iszero(q.r[2]),
-            iszero(j) && !iszero(q.r[3]),j==length(y) && !iszero(q.r[4]),
-            iszero(k) && !iszero(q.r[5]),k==length(z) && !iszero(q.r[6]),
-            iszero(l) && !iszero(q.r[7]),l==length(w) && !iszero(q.r[8]))
-        if i1 || i2 || j1 || j2 || k1 || k2 || l1 || l2
-            return m(Chain{V}(
-                reposition(i1,i2,q.p[1],q.p[2],x,t[1]),
-                reposition(j1,j2,q.p[3],q.p[4],y,t[2]),
-                reposition(k1,k2,q.p[5],q.p[6],z,t[3]),
-                reposition(l1,l2,q.p[7],q.p[8],w,t[4])))
+        i01,i02,iq1,iq2,j01,j02,jq1,jq2,k01,k02,kq1,kq2,l01,l02,lq1,lq2 = @inbounds (
+            iszero(i),i==length(x),iszero(q.r[1]),iszero(q.r[2]),
+            iszero(j),j==length(y),iszero(q.r[3]),iszero(q.r[4]),
+            iszero(k),k==length(z),iszero(q.r[5]),iszero(q.r[6]),
+            iszero(l),l==length(w),iszero(q.r[7]),iszero(q.r[8]))
+        if (i01 && iq1) || (i02 && iq2) || (j01 && jq1) || (j02 && jq2) || (k01 && kq1) || (k02 && kq2) || (l01 && lq1) || (l02 && lq2)
+            return zero(F)
+        else
+            i1,i2,j1,j2,k1,k2,l1,l2 = (
+                i01 && !iq1,i02 && !iq2,
+                j01 && !jq1,j02 && !jq2,
+                k01 && !kq1,k02 && !kq2,
+                l01 && !lq1,l02 && !lq2)
+            if i1 || i2 || j1 || j2 || k1 || k2 || l1 || l2
+                return m(Chain{V}(
+                    (@inbounds reposition(i1,i2,q.p[1],q.p[2],x,t1)),
+                    (@inbounds reposition(j1,j2,q.p[3],q.p[4],y,t2)),
+                    (@inbounds reposition(k1,k2,q.p[5],q.p[6],z,t3)),
+                    (@inbounds reposition(l1,l2,q.p[7],q.p[8],w,t4))))
+            end
         end
+    elseif iszero(i) || iszero(j) || iszero(k) || iszero(l) || i==length(x) || j==length(y) || k==length(z) || l==length(w)
+        return zero(F) # elseif condition creates 1 allocation, as opposed to 0 ???
     end
     #f1 = linterp(t[1],x[i],x[i+1],m.cod[i,j,k,l],m.cod[i+1,j,k,l])
     #f2 = linterp(t[1],x[i],x[i+1],m.cod[i,j+1,k,l],m.cod[i+1,j+1,k,l])
@@ -406,7 +473,7 @@ function (m::TensorField{B,F,N,<:RealSpace{4}} where {B,F,N})(t::Chain{V}) where
     #g4 = linterp(t[2],y[j],y[j+1],f7,f8)
     #h2 = linterp(t[3],z[k],z[k+1],g3,g4)
     #linterp(t[4],w[l],w[l+1],h1,h2)
-    quadlinterp(t[1],t[2],t[3],t[4],x[i],x[i+1],y[j],y[j+1],z[k],z[k+1],w[l],w[l+1],
+    quadlinterp(t1,t2,t3,t4,x[i],x[i+1],y[j],y[j+1],z[k],z[k+1],w[l],w[l+1],
         m.cod[i,j,k,l],m.cod[i+1,j,k,l],m.cod[i,j+1,k,l],m.cod[i+1,j+1,k,l],
         m.cod[i,j,k+1,l],m.cod[i+1,j,k+1,l],m.cod[i,j+1,k+1,l],m.cod[i+1,j+1,k+1,l],
         m.cod[i,j,k,l+1],m.cod[i+1,j,k,l+1],m.cod[i,j+1,k,l+1],m.cod[i+1,j+1,k,l+1],
@@ -414,32 +481,44 @@ function (m::TensorField{B,F,N,<:RealSpace{4}} where {B,F,N})(t::Chain{V}) where
 end
 
 (m::TensorField{B,F,N,<:RealSpace{5}} where {B,F,N})(x,y,z,w,v) = m(Chain(x,y,z,w,v))
-function (m::TensorField{B,F,N,<:RealSpace{5}} where {B,F,N})(t::Chain{V}) where V
-    x,y,z,w,v = points(m).v[1],points(m).v[2],points(m).v[3],points(m).v[4],points(m).v[5]
+function (m::TensorField{B,F,N,<:RealSpace{5}} where {B,N})(t::Chain{V,G,T,5} where {G,T}) where {F,V}
+    x,y,z,w,v,t1,t2,t3,t4,t5 = @inbounds (points(m).v[1],points(m).v[2],points(m).v[3],points(m).v[4],points(m).v[5],t[1],t[2],t[3],t[4],t[5])
     i,j,k,l,o = (
-        searchsortedfirst(x,t[1])-1,
-        searchsortedfirst(y,t[2])-1,
-        searchsortedfirst(z,t[3])-1,
-        searchsortedfirst(w,t[4])-1,
-        searchsortedfirst(v,t[5])-1)
+        searchsortedfirst(x,t1)-1,
+        searchsortedfirst(y,t2)-1,
+        searchsortedfirst(z,t3)-1,
+        searchsortedfirst(w,t4)-1,
+        searchsortedfirst(v,t5)-1)
     if !isopen(m)
         q = immersion(m)
-        i1,i2,j1,j2,k1,k2,l1,l2,o1,o2 = (
-            iszero(i) && !iszero(q.r[1]),i==length(x) && !iszero(q.r[2]),
-            iszero(j) && !iszero(q.r[3]),j==length(y) && !iszero(q.r[4]),
-            iszero(k) && !iszero(q.r[5]),k==length(z) && !iszero(q.r[6]),
-            iszero(l) && !iszero(q.r[7]),l==length(w) && !iszero(q.r[8]),
-            iszero(o) && !iszero(q.r[9]),o==length(v) && !iszero(q.r[10]))
-        if i1 || i2 || j1 || j2 || k1 || k2 || l1 || l2 || o1 || o2
-            return m(Chain{V}(
-                reposition(i1,i2,q.p[1],q.p[2],x,t[1]),
-                reposition(j1,j2,q.p[3],q.p[4],y,t[2]),
-                reposition(k1,k2,q.p[5],q.p[6],z,t[3]),
-                reposition(l1,l2,q.p[7],q.p[8],w,t[4]),
-                reposition(o1,o2,q.p[9],q.p[10],v,t[5])))
+        i01,i02,iq1,iq2,j01,j02,jq1,jq2,k01,k02,kq1,kq2,l01,l02,lq1,lq2,o01,o02,oq1,oq2 = @inbounds (
+            iszero(i),i==length(x),iszero(q.r[1]),iszero(q.r[2]),
+            iszero(j),j==length(y),iszero(q.r[3]),iszero(q.r[4]),
+            iszero(k),k==length(z),iszero(q.r[5]),iszero(q.r[6]),
+            iszero(l),l==length(w),iszero(q.r[7]),iszero(q.r[8]),
+            iszero(o),o==length(v),iszero(q.r[9]),iszero(q.r[10]))
+        if (i01 && iq1) || (i02 && iq2) || (j01 && jq1) || (j02 && jq2) || (k01 && kq1) || (k02 && kq2) || (l01 && lq1) || (l02 && lq2) || (o01 && oq1) || (o02 && oq2)
+            return zero(F)
+        else
+            i1,i2,j1,j2,k1,k2,l1,l2,o1,o2 = (
+                i01 && !iq1,i02 && !iq2,
+                j01 && !jq1,j02 && !jq2,
+                k01 && !kq1,k02 && !kq2,
+                l01 && !lq1,l02 && !lq2,
+                o01 && !oq1,o02 && !oq2)
+            if i1 || i2 || j1 || j2 || k1 || k2 || l1 || l2 || o1 || o2
+                return m(Chain{V}(
+                    (@inbounds reposition(i1,i2,q.p[1],q.p[2],x,t1)),
+                    (@inbounds reposition(j1,j2,q.p[3],q.p[4],y,t2)),
+                    (@inbounds reposition(k1,k2,q.p[5],q.p[6],z,t3)),
+                    (@inbounds reposition(l1,l2,q.p[7],q.p[8],w,t4)),
+                    (@inbounds reposition(o1,o2,q.p[9],q.p[10],v,t5))))
+            end
         end
+    elseif iszero(i) || iszero(j) || iszero(k) || iszero(l) || iszero(o) || i==length(x) || j==length(y) || k==length(z) || l==length(w) || o==length(v)
+        return zero(F) # elseif condition creates 1 allocation, as opposed to 0 ???
     end
-    quintlinterp(t[1],t[2],t[3],t[4],t[5],x[i],x[i+1],y[j],y[j+1],z[k],z[k+1],w[l],w[l+1],v[o],v[o+1],
+    quintlinterp(t1,t2,t3,t4,t5,x[i],x[i+1],y[j],y[j+1],z[k],z[k+1],w[l],w[l+1],v[o],v[o+1],
         m.cod[i,j,k,l,o],m.cod[i+1,j,k,l,o],m.cod[i,j+1,k,l,o],m.cod[i+1,j+1,k,l,o],
         m.cod[i,j,k+1,l,o],m.cod[i+1,j,k+1,l,o],m.cod[i,j+1,k+1,l,o],m.cod[i+1,j+1,k+1,l,o],
         m.cod[i,j,k,l+1,o],m.cod[i+1,j,k,l+1,o],m.cod[i,j+1,k,l+1,o],m.cod[i+1,j+1,k,l+1,o],
@@ -526,25 +605,24 @@ checkdomain(a::GlobalFiber,b::GlobalFiber) = domain(a)≠domain(b) ? error("Glob
 
 Variation(dom,cod::TensorField) = TensorField(dom,cod.(dom))
 function Variation(cod::TensorField); p = points(cod).v[end]
-    TensorField(p,cod.(1:length(p)))
+    TensorField(p,slice.(Ref(cod),1:length(p)))
 end
 function Variation(cod::TensorField{B,F,2,<:FiberProductBundle} where {B,F})
     p = base(cod).cod.v[1]
-    TensorField(p,cod.(1:length(p)))
+    TensorField(p,slice.(Ref(cod),1:length(p)))
 end
 const variation = Variation
 
 alteration(dom,cod::TensorField) = TensorField(dom,cod.(dom))
 function alteration(cod::TensorField); p = points(cod).v[1]
-    TensorField(p,cod.(1:length(p),1))
+    TensorField(p,slice.(Ref(cod),1:length(p),1))
 end
 
 modification(dom,cod::TensorField) = TensorField(dom,cod.(dom))
 function modification(cod::TensorField); p = points(cod).v[2]
-    TensorField(p,cod.(1:length(p),2))
+    TensorField(p,slice.(Ref(cod),1:length(p),2))
 end
 
-
 include("diffgeo.jl")
 include("constants.jl")
 include("element.jl")
@@ -566,27 +644,27 @@ function variation(v::Variation,t,fun::Function,fun!::Function,args...)
 end
 function variation(v::TensorField,t,fun::Function,args...)
     for i ∈ 1:length(points(v).v[end])
-        display(fun(v(i),args...))
+        display(fun(slice(v,i),args...))
         sleep(t)
     end
 end
 function variation(v::TensorField,t,fun::Function,fun!::Function,args...)
-    display(fun(v(1),args...))
+    display(fun(slice(v,1),args...))
     for i ∈ 2:length(points(v).v[end])
-        fun!(v(i),args...)
+        fun!(slice(v,i),args...)
         sleep(t)
     end
 end
 function variation(v::TensorField{B,F,2,<:FiberProductBundle} where {B,F},t,fun::Function,args...)
     for i ∈ 1:length(base(v).cod.v[1])
-        display(fun(v(i),args...))
+        display(fun(slice(v,i),args...))
         sleep(t)
     end
 end
 function variation(v::TensorField{B,F,2,<:FiberProductBundle} where {B,F},t,fun::Function,fun!::Function,args...)
-    display(fun(v(1),args...))
+    display(fun(slice(v,1),args...))
     for i ∈ 2:length(base(v).cod.v[1])
-        fun!(v(i),args...)
+        fun!(slice(v,i),args...)
         sleep(t)
     end
 end
@@ -595,14 +673,14 @@ alteration(v::TensorField,fun::Function,args...) = alteration(v,0.0,fun,args...)
 alteration(v::TensorField,fun::Function,fun!::Function,args...) = alteration(v,0.0,fun,fun!,args...)
 function alteration(v::TensorField,t,fun::Function,args...)
     for i ∈ 1:length(points(v).v[1])
-        display(fun(v(i,1),args...))
+        display(fun(slice(v,i,1),args...))
         sleep(t)
     end
 end
 function alteration(v::TensorField,t,fun::Function,fun!::Function,args...)
-    display(fun(v(1,1),args...))
+    display(fun(slice(v,1,1),args...))
     for i ∈ 2:length(points(v).v[1])
-        fun!(v(i,1),args...)
+        fun!(slice(v,i,1),args...)
         sleep(t)
     end
 end
@@ -611,14 +689,14 @@ modification(v::TensorField,fun::Function,args...) = modification(v,0.0,fun,args
 modification(v::TensorField,fun::Function,fun!::Function,args...) = modification(v,0.0,fun,fun!,args...)
 function modification(v::TensorField,t,fun::Function,args...)
     for i ∈ 1:length(points(v).v[2])
-        display(fun(v(i,2),args...))
+        display(fun(slice(v,i,2),args...))
         sleep(t)
     end
 end
 function modification(v::TensorField,t,fun::Function,fun!::Function,args...)
-    display(fun(v(1,2),args...))
+    display(fun(slice(v,1,2),args...))
     for i ∈ 2:length(points(v).v[2])
-        fun!(v(i,2),args...)
+        fun!(slice(v,i,2),args...)
         sleep(t)
     end
 end
@@ -794,12 +872,12 @@ function __init__()
             alteration(M,Makie.lines!,Makie.lines!)
         end
         function linegraph(v::TensorField{B,<:Chain,3,<:GridFrameBundle} where B,args...)
-            display(Makie.lines(slice(v,1,1,1),args...))
+            display(Makie.lines(slice2(v,1,1,1),args...))
             c = (2,3),(1,3),(1,2)
             for k ∈ (1,2,3)
                 for i ∈ 1:length(points(v).v[c[k][1]])
                     for j ∈ 1:length(points(v).v[c[k][2]])
-                        Makie.lines!(slice(v,i,j,k),args...)
+                        Makie.lines!(slice2(v,i,j,k),args...)
                     end
                 end
             end
diff --git a/src/diffgeo.jl b/src/diffgeo.jl
index 3e6898d..f941bd8 100644
--- a/src/diffgeo.jl
+++ b/src/diffgeo.jl
@@ -63,10 +63,10 @@ end
 end=#
 Base.getindex(g::GridFrameBundle{C,M,<:FiberBundle,<:OpenTopology} where C,j::Int,n::Val,i::Vararg{Int}) where M = getpoint(g,j,n,i...)
 @generated function getpoint(g::GridFrameBundle{C,M,<:FiberBundle} where C,j::Int,::Val{N},i::Vararg{Int}) where {M,N}
-    :(Base.getindex(points(g),$([k≠N ? :(i[$k]) : :(i[$k]+j) for k ∈ 1:M]...)))
+    :(Base.getindex(points(g),$([k≠N ? :(@inbounds i[$k]) : :(@inbounds i[$k]+j) for k ∈ 1:M]...)))
 end
 @generated function Base.getindex(g::GridFrameBundle{C,M} where C,j::Int,n::Val{N},i::Vararg{Int}) where {M,N}
-    :(Base.getindex(points(g),Base.getindex(immersion(g),n,$([k≠N ? :(i[$k]) : :(i[$k]+j) for k ∈ 1:M]...))...))
+    :(Base.getindex(points(g),Base.getindex(immersion(g),n,$([k≠N ? :(@inbounds i[$k]) : :(@inbounds i[$k]+j) for k ∈ 1:M]...))...))
 end
 
 # centraldiff
@@ -294,14 +294,14 @@ for fun ∈ (:_slow,:_fast)
         function $cdg(f::Grid{2},n::Val{1},dt::AbstractVector,s::Tuple=size(f))
             d = Array{pointtype(f),2}(undef,s...)
             @threads for i ∈ OneTo(s[1]); for j ∈ OneTo(s[2])
-                d[i,j] = $cdg(f,s[1],n,i,j)/dt[i]
+                d[i,j] = $cdg(f,(@inbounds s[1]),n,i,j)/dt[i]
             end end
             return d
         end
         function $cdg(f::Grid{2},n::Val{2},dt::AbstractVector,s::Tuple=size(f))
             d = Array{pointtype(f),2}(undef,s...)
             @threads for i ∈ OneTo(s[1]); for j ∈ OneTo(s[2])
-                d[i,j] = $cdg(f,s[2],n,i,j)/dt[j]
+                d[i,j] = $cdg(f,(@inbounds s[2]),n,i,j)/dt[j]
             end end
             return d
         end
@@ -322,21 +322,21 @@ for fun ∈ (:_slow,:_fast)
         function $cdg(f::Grid{3},n::Val{1},dt::AbstractVector,s::Tuple=size(f))
             d = Array{pointtype(f),3}(undef,s...)
             @threads for i ∈ OneTo(s[1]); for j ∈ OneTo(s[2]); for k ∈ OneTo(s[3])
-                d[i,j,k] = $cdg(f,s[1],n,i,j,k)/dt[i]
+                d[i,j,k] = $cdg(f,(@inbounds s[1]),n,i,j,k)/dt[i]
             end end end
             return d
         end
         function $cdg(f::Grid{3},n::Val{2},dt::AbstractVector,s::Tuple=size(f))
             d = Array{pointtype(f),3}(undef,s...)
             @threads for i ∈ OneTo(s[1]); for j ∈ OneTo(s[2]); for k ∈ OneTo(s[3])
-                d[i,j,k] = $cdg(f,s[2],n,i,j,k)/dt[j]
+                d[i,j,k] = $cdg(f,(@inbounds s[2]),n,i,j,k)/dt[j]
             end end end
             return d
         end
         function $cdg(f::Grid{3},n::Val{3},dt::AbstractVector,s::Tuple=size(f))
             d = Array{pointtype(f),3}(undef,s...)
             @threads for i ∈ OneTo(s[1]); for j ∈ OneTo(s[2]); for k ∈ OneTo(s[3])
-                d[i,j,k] = $cdg(f,s[3],n,i,j,k)/dt[k]
+                d[i,j,k] = $cdg(f,(@inbounds s[3]),n,i,j,k)/dt[k]
             end end end
             return d
         end
@@ -357,28 +357,28 @@ for fun ∈ (:_slow,:_fast)
         function $cdg(f::Grid{4},n::Val{1},dt::AbstractVector,s::Tuple=size(f))
             d = Array{pointtype(f),4}(undef,s...)
             @threads for i ∈ OneTo(s[1]); for j ∈ OneTo(s[2]); for k ∈ OneTo(s[3]); for l ∈ OneTo(s[4])
-                d[i,j,k,l] = $cdg(f,s[1],n,i,j,k,l)/dt[i]
+                d[i,j,k,l] = $cdg(f,(@inbounds s[1]),n,i,j,k,l)/dt[i]
             end end end end
             return d
         end
         function $cdg(f::Grid{4},n::Val{2},dt::AbstractVector,s::Tuple=size(f))
             d = Array{pointtype(f),4}(undef,s...)
             @threads for i ∈ OneTo(s[1]); for j ∈ OneTo(s[2]); for k ∈ OneTo(s[3]); for l ∈ OneTo(s[4])
-                d[i,j,k,l] = $cdg(f,s[2],n,i,j,k,l)/dt[j]
+                d[i,j,k,l] = $cdg(f,(@inbounds s[2]),n,i,j,k,l)/dt[j]
             end end end end
             return d
         end
         function $cdg(f::Grid{4},n::Val{3},dt::AbstractVector,s::Tuple=size(f))
             d = Array{pointtype(f),4}(undef,s...)
             @threads for i ∈ OneTo(s[1]); for j ∈ OneTo(s[2]); for k ∈ OneTo(s[3]); for l ∈ OneTo(s[4])
-                d[i,j,k,l] = $cdg(f,s[3],n,i,j,k,l)/dt[k]
+                d[i,j,k,l] = $cdg(f,(@inbounds s[3]),n,i,j,k,l)/dt[k]
             end end end end
             return d
         end
         function $cdg(f::Grid{4},n::Val{4},dt::AbstractVector,s::Tuple=size(f))
             d = Array{pointtype(f),4}(undef,s...)
             @threads for i ∈ OneTo(s[1]); for j ∈ OneTo(s[2]); for k ∈ OneTo(s[3]); for l ∈ OneTo(s[4])
-                d[i,j,k,l] = $cdg(f,s[4],n,i,j,k,l)/dt[l]
+                d[i,j,k,l] = $cdg(f,(@inbounds s[4]),n,i,j,k,l)/dt[l]
             end end end end
             return d
         end
@@ -399,35 +399,35 @@ for fun ∈ (:_slow,:_fast)
         function $cdg(f::Grid{5},n::Val{1},dt::AbstractVector,s::Tuple=size(f))
             d = Array{pointtype(f),5}(undef,s...)
             @threads for i ∈ OneTo(s[1]); for j ∈ OneTo(s[2]); for k ∈ OneTo(s[3]); for l ∈ OneTo(s[4]); for o ∈ OneTo(s[5])
-                d[i,j,k,l,o] = $cdg(f,s[1],n,i,j,k,l,o)/dt[i]
+                d[i,j,k,l,o] = $cdg(f,(@inbounds s[1]),n,i,j,k,l,o)/dt[i]
             end end end end end
             return d
         end
         function $cdg(f::Grid{5},n::Val{2},dt::AbstractVector,s::Tuple=size(f))
             d = Array{pointtype(f),5}(undef,s...)
             @threads for i ∈ OneTo(s[1]); for j ∈ OneTo(s[2]); for k ∈ OneTo(s[3]); for l ∈ OneTo(s[4]); for o ∈ OneTo(s[5])
-                d[i,j,k,l,o] = $cdg(f,s[2],n,i,j,k,l,o)/dt[j]
+                d[i,j,k,l,o] = $cdg(f,(@inbounds s[2]),n,i,j,k,l,o)/dt[j]
             end end end end end
             return d
         end
         function $cdg(f::Grid{5},n::Val{3},dt::AbstractVector,s::Tuple=size(f))
             d = Array{pointtype(f),5}(undef,s...)
             @threads for i ∈ OneTo(s[1]); for j ∈ OneTo(s[2]); for k ∈ OneTo(s[3]); for l ∈ OneTo(s[4]); for o ∈ OneTo(s[5])
-                d[i,j,k,l,o] = $cdg(f,s[3],n,i,j,k,l,o)/dt[k]
+                d[i,j,k,l,o] = $cdg(f,(@inbounds s[3]),n,i,j,k,l,o)/dt[k]
             end end end end end
             return d
         end
         function $cdg(f::Grid{5},n::Val{4},dt::AbstractVector,s::Tuple=size(f))
             d = Array{pointtype(f),5}(undef,s...)
             @threads for i ∈ OneTo(s[1]); for j ∈ OneTo(s[2]); for k ∈ OneTo(s[3]); for l ∈ OneTo(s[4]); for o ∈ OneTo(s[5])
-                d[i,j,k,l,o] = $cdg(f,s[4],n,i,j,k,l,o)/dt[l]
+                d[i,j,k,l,o] = $cdg(f,(@inbounds s[4]),n,i,j,k,l,o)/dt[l]
             end end end end end
             return d
         end
         function $cdg(f::Grid{5},n::Val{5},dt::AbstractVector,s::Tuple=size(f))
             d = Array{pointtype(f),5}(undef,s...)
             @threads for i ∈ OneTo(s[1]); for j ∈ OneTo(s[2]); for k ∈ OneTo(s[3]); for l ∈ OneTo(s[4]); for o ∈ OneTo(s[5])
-                d[i,j,k,l,o] = $cdg(f,s[5],n,i,j,k,l,o)/dt[o]
+                d[i,j,k,l,o] = $cdg(f,(@inbounds s[5]),n,i,j,k,l,o)/dt[o]
             end end end end end
             return d
         end
@@ -453,9 +453,9 @@ for fun ∈ (:_slow,:_fast)
         $cd(f::GridFrameBundle{Coordinate{P,G},N,<:PointArray{P,G,N,<:RealRegion,<:Global{N,<:InducedMetric}},<:ProductTopology}) where {P,G,N} = ProductSpace($cd.(base(base(f)).v))
         $cd(f::GridFrameBundle{Coordinate{P,G},N,<:PointArray{P,G,N,<:RealRegion,<:Global{N,<:InducedMetric}},<:OpenTopology}) where {P,G,N} = ProductSpace($cd.(base(base(f)).v))
         $cd(f::GridFrameBundle{Coordinate{P,G},N,<:PointArray{P,G,N,<:RealRegion,<:Global{N,<:InducedMetric}}}) where {P,G,N} = sum.(value.($cdg(f)))
-        $cd(f::GridFrameBundle{Coordinate{P,G},N,<:PointArray{P,G,N,<:RealRegion},<:ProductTopology}) where {P,G,N} = applymetric.($cd(base(base(f))),fiber(f))
-        $cd(f::GridFrameBundle{Coordinate{P,G},N,<:PointArray{P,G,N,<:RealRegion},<:OpenTopology}) where {P,G,N} = applymetric.($cd(base(base(f))),fiber(f))
-        $cd(f::GridFrameBundle{Coordinate{P,G},N,<:PointArray{P,G,N,<:RealRegion}}) where {P,G,N} = applymetric.(sum.(value.($cdg(f))),fiber(f))
+        #$cd(f::GridFrameBundle{Coordinate{P,G},N,<:PointArray{P,G,N,<:RealRegion},<:ProductTopology}) where {P,G,N} = applymetric.($cd(base(base(f))),metricextensor(f))
+        $cd(f::GridFrameBundle{Coordinate{P,G},N,<:PointArray{P,G,N,<:RealRegion},<:OpenTopology}) where {P,G,N} = applymetric.($cd(base(base(f))),metricextensor(f))
+        $cd(f::GridFrameBundle{Coordinate{P,G},N,<:PointArray{P,G,N,<:RealRegion}}) where {P,G,N} = applymetric.(sum.(value.($cdg(f))),metricextensor(f))
         function $cd(f::AbstractRange,s::Tuple=size(f))
             d = Vector{eltype(f)}(undef,s[1])
             @threads for i ∈ OneTo(s[1])
@@ -473,7 +473,14 @@ for fun ∈ (:_slow,:_fast)
     end
 end
 
-applymetric(f::Chain{V,G},g::DiagonalOperator{V,<:Multivector}) where {V,G} = Chain{V,G}(value(f)./sqrt.(value(value(g)(Val(G)))))
+applymetric(f::Chain{V,G},g::DiagonalOperator{W,<:Multivector} where W) where {V,G} = Chain{V,G}(value(f)./sqrt.(value(value(g)(Val(G)))))
+applymetric(f::Chain{V,G},g::DiagonalOperator{W,<:Chain} where W) where {V,G} = Chain{V,G}(value(f)./sqrt.(value(value(g))))
+applymetric(f::Chain{V,G},g::Outermorphism) where {V,G} = applymetric(f,(@inbounds value(g)[1]))
+applymetric(f::Chain{V,G},g::Endomorphism{W,<:Simplex} where W) where {V,G} = applymetric(f,value(g))
+@generated function applymetric(x::Chain{V,G,T,N} where {G,T},g::Simplex) where {V,N}
+    Expr(:call,:(Chain{V}),[:(x[$k]/sqrt(g[$k,$k])) for k ∈ list(1,N)]...)
+end
+
 
 function centraldiff_slow_calc(f::Grid{M,T,PA,<:OpenTopology} where {M,T,PA},l::Int,n::Val{N},i::Vararg{Int}) where N #l=size(f)[N]
     if isone(i[N])
@@ -941,7 +948,7 @@ function planecurve(κ::RealFunction,φ::Real=0.0)
     integral(Chain.(cos(int),sin(int)))
 end
 
-export surfacemetric, surfaceframe
+export surfacemetric, surfaceframe, firstkind, secondkind, geodesicsystem, applymetric
 
 surfacemetric(dom::ScalarField,f::Function) = surfacemetric(TensorField(dom,f))
 function surfacemetric(t::ScalarField)
@@ -969,10 +976,24 @@ function surfaceframe(t)
     TensorOperator.(Chain.(Chain.(E,F)./mag,Chain.(F,.-E)./(sig.*mag)))
 end
 
-firstkind(dg,k,i,j) = dg[k,j][i] + dg[i,k][j] - dg[i,j][k]
+_firstkind(dg,k,i,j) = dg[k,j][i] + dg[i,k][j] - dg[i,j][k]
+firstkind(g::TensorField) = TensorField(base(g),TensorOperator.(firstkind.(d(g)/2)))
+@generated function firstkind(dg,i,j,k)
+    Expr(:call,:+,[:(_firstkind(dg,$l,i,j)) for l ∈ list(1,mdims(fibertype(dg)))]...)
+end
+@generated function firstkind(dg,i,j)
+    Expr(:call,:Chain,[:(firstkind(dg,i,j,$k)) for k ∈ list(1,mdims(fibertype(dg)))]...)
+end
+@generated function firstkind(dg,j)
+    Expr(:call,:Chain,[:(firstkind(dg,$i,j)) for i ∈ list(1,mdims(fibertype(dg)))]...)
+end
+@generated function firstkind(dg)
+    Expr(:call,:Chain,[:(firstkind(dg,$j)) for j ∈ list(1,mdims(fibertype(dg)))]...)
+end
+
 secondkind(g) = TensorField(base(g),TensorOperator.(secondkind.(inv.(g),d(g)/2)))
 @generated function secondkind(ig,dg,i,j,k)
-    Expr(:call,:+,[:(ig[k,$l]*firstkind(dg,$l,i,j)) for l ∈ list(1,mdims(fibertype(dg)))]...)
+    Expr(:call,:+,[:(ig[k,$l]*_firstkind(dg,$l,i,j)) for l ∈ list(1,mdims(fibertype(dg)))]...)
 end
 @generated function secondkind(ig,dg,i,j)
     Expr(:call,:Chain,[:(secondkind(ig,dg,i,j,$k)) for k ∈ list(1,mdims(fibertype(dg)))]...)
@@ -984,6 +1005,17 @@ end
     Expr(:call,:Chain,[:(secondkind(ig,dg,$j)) for j ∈ list(1,mdims(fibertype(dg)))]...)
 end
 
+geodesicsystem(x::Chain,Γ) = Chain(x,-geodesic(x,Γ))
+geodesicsystem(x::Chain,Γ,g) = Chain(applymetric(x,g),applymetric(-geodesic(x,Γ),g))
+geodesicsystem(x::Chain,Γ::TensorField) = geodesicsystem(x,Γ(x))
+geodesicsystem(x::Chain,Γ::TensorField,g::TensorField) = geodesicsystem(x,Γ(x),g(x))
+@generated function geodesic(x::Chain{V,G,T,N} where {V,G,T},Γ) where N
+    Expr(:call,:+,vcat([[:(Γ[$i,$j]*(x[$i]*x[$j])) for i ∈ list(1,N)] for j ∈ list(1,N)]...)...)
+end
+@generated function metricscale(x::Chain{V,G,T,N} where {G,T},g::Simplex) where {V,N}
+    Expr(:call,:(Chain{V}),[:(x[$k]*sqrt(g[$k,$k])) for k ∈ list(1,N)]...)
+end
+
 #export beta, betafunction
 
 function beta(a,b,n=30000)
diff --git a/src/topology.jl b/src/topology.jl
index faae902..92a0c73 100644
--- a/src/topology.jl
+++ b/src/topology.jl
@@ -13,10 +13,12 @@
 #   https://crucialflow.com
 
 export FiberProduct, FiberProductBundle, HomotopyBundle
-export LocalSection, GlobalFiber, LocalFiber
+export LocalSection, GlobalFiber, LocalFiber, localfiber, globalfiber
 export base, fiber, domain, codomain, ↦, →, ←, ↤, basetype, fibertype, graph
 export ProductSpace, RealRegion, NumberLine, Rectangle, Hyperrectangle, ⧺, ⊕
 
+resize_lastdim!(x::Vector,i) = resize!(x,i)
+
 # ProductSpace
 
 struct ProductSpace{V,T,N,M,S} <: AbstractArray{Chain{V,1,T,N},N}
@@ -42,9 +44,11 @@ Base.show(io::IO,t::RealRegion{V,T,N,<:AbstractRange} where {V,T,N}) = print(io,
 Base.iterate(t::RealRegion) = (getindex(t,1),1)
 Base.iterate(t::RealRegion,state) = (s=state+1; s≤length(t) ? (getindex(t,s),s) : nothing)
 
-@generated Base.size(m::RealRegion{V}) where V = :(($([:(size(m.v[$i])...) for i ∈ 1:mdims(V)]...),))
-@generated Base.getindex(m::RealRegion{V,T,N},i::Vararg{Int}) where {V,T,N} = :(Chain{V,1,T}(Values{N,T}($([:(m.v[$j][i[$j]]) for j ∈ 1:N]...))))
-Base.getindex(m::NumberLine{V,T},i::Int) where {V,T} = Chain{V,1,T}(Values((m.v[1][i],)))
+resize_lastdim!(m::ProductSpace,i) = (resize!(m.v[end],i); m)
+
+@generated Base.size(m::RealRegion{V}) where V = :(($([:(size(@inbounds m.v[$i])...) for i ∈ 1:mdims(V)]...),))
+@generated Base.getindex(m::RealRegion{V,T,N},i::Vararg{Int}) where {V,T,N} = :(Chain{V,1,T}(Values{N,T}($([:((@inbounds m.v[$j])[@inbounds i[$j]]) for j ∈ 1:N]...))))
+Base.getindex(m::NumberLine{V,T},i::Int) where {V,T} = Chain{V,1,T}(Values(((@inbounds m.v[1])[i],)))
 @pure Base.getindex(t::RealRegion,i::CartesianIndex) = getindex(t,i.I...)
 @pure Base.eltype(::Type{<:ProductSpace{V,T,N}}) where {V,T,N} = Chain{V,1,T,N}
 
@@ -64,17 +68,20 @@ end
 
 ⊕(a::AbstractVector{<:Real}...) = RealRegion(Values(a))
 ⊕(a::ProductSpace,b::AbstractVector{<:Real}) = RealRegion(Values(a.v...,b))
+⊕(a::ProductSpace,b::ProductSpace) = RealRegion(Values(a.v...,b.v...))
+cross(a::ProductSpace,b::AbstractVector{<:Real}) = a⊕b
+cross(a::ProductSpace,b::ProductSpace) = a⊕b
 
 @generated ⧺(a::Real...) = :(Chain($([:(a[$i]) for i ∈ 1:length(a)]...)))
 @generated ⧺(a::Complex...) = :(Chain($([:(a[$i]) for i ∈ 1:length(a)]...)))
 ⧺(a::Chain{A,G},b::Chain{B,G}) where {A,B,G} = Chain{A∪B,G}(vcat(a.v,b.v))
 
-remove(t::ProductSpace{V,T,2} where {V,T},::Val{1}) = t.v[2]
-remove(t::ProductSpace{V,T,2} where {V,T},::Val{2}) = t.v[1]
+remove(t::ProductSpace{V,T,2} where {V,T},::Val{1}) = (@inbounds t.v[2])
+remove(t::ProductSpace{V,T,2} where {V,T},::Val{2}) = (@inbounds t.v[1])
 @generated remove(t::ProductSpace{V,T,N} where {V,T},::Val{J}) where {N,J} = :(ProductSpace(domain(t).v[$(Values([i for i ∈ 1:N if i≠J]...))]))
 
 # 1
-(m::ProductSpace)(c::Colon,i::Int...) = m.v[1]
+(m::ProductSpace)(c::Colon,i::Int...) = (@inbounds m.v[1])
 (m::ProductSpace)(i::Int,c::Colon,j::Int...) = m.v[2]
 (m::ProductSpace)(i::Int,j::Int,c::Colon,k::Int...) = m.v[3]
 (m::ProductSpace)(i::Int,j::Int,k::Int,c::Colon,l::Int...) = m.v[4]
@@ -122,12 +129,15 @@ export CrossRange
 struct CrossRange <: AbstractVector{Int}
     n::Int
     m::Int
-    CrossRange(n::Int) = new(n,Int((isodd(n) ? n+1 : n)/2)-1)
+    CrossRange(n::Int) = new(n,crossrange(n))
 end
 
+crossrange(n) = Int((isodd(n) ? n+1 : n)/2)-1
+
 Base.iterate(t::CrossRange) = (getindex(t,1),1)
 Base.iterate(t::CrossRange,state) = (s=state+1; s≤length(t) ? (getindex(t,s),s) : nothing)
 
+#Base.resize!(m::CrossRange,i) = (m.n = i; m.m = crossrange(n))
 Base.size(m::CrossRange) = (length(m),)
 Base.length(m::CrossRange) = m.n
 Base.getindex(m::CrossRange,i::Int) = i≤m.m ? i+m.m : i-m.m
@@ -153,7 +163,7 @@ ProductTopology(i::Int) = ProductTopology(Values((OneTo(i),)))
 ProductTopology(i::AbstractVector) = ProductTopology(Values((i,)))
 ProductTopology(i::AbstractVector,jk::AbstractVector...) = ProductTopology(Values(i,jk...))
 
-Base.show(io::IO,t::ProductTopology{N,<:AbstractRange} where N) = print(io,'(',Values(getindex.(t.v,1)),"):(",Values(Number.(getproperty.(t.v,:step))),"):(",Values(getindex.(t.v,length.(t.v))),')')
+Base.show(io::IO,t::ProductTopology{N,<:AbstractRange} where N) = print(io,Values(getindex.(t.v,1)),':',Values(getindex.(t.v,length.(t.v))))
 
 (::Base.Colon)(min::Values{N,Int},max::Values{N,Int}) where N = ProductTopology(Colon().(value(min),value(max)))
 (::Base.Colon)(min::Values{N,Int},step::Values{N,Int},max::Values{N,Int}) where N = ProductTopology(Colon().(value(min),value(step),value(max)))
@@ -161,9 +171,17 @@ Base.show(io::IO,t::ProductTopology{N,<:AbstractRange} where N) = print(io,'(',V
 Base.iterate(t::ProductTopology) = (getindex(t,1),1)
 Base.iterate(t::ProductTopology,state) = (s=state+1; s≤length(t) ? (getindex(t,s),s) : nothing)
 
-@generated Base.size(m::ProductTopology{N}) where N = :(($([:(size(m.v[$i])...) for i ∈ 1:N]...),))
-@generated Base.getindex(m::ProductTopology{N},i::Vararg{Int}) where N = :(Values{N,Int}($([:(m.v[$j][i[$j]]) for j ∈ 1:N]...)))
-Base.getindex(m::ProductTopology{1},i::Int) = Values((m.v[1][i],))
+resize(::OneTo,i) = OneTo(i)
+resize(m::StepRange,i) = isone(m.start) ? (1:1:i) : (i:-1:1)
+resize(m::CrossRange,i) = CrossRange(i)
+resize(m::ProductTopology{1},i) = ProductTopology(@inbounds resize(m.v[1],i))
+@generated function resize(m::ProductTopology{N},i) where N
+    Expr(:call,:ProductTopology,Expr(:call,:Values,[j≠N ? :(@inbounds m.v[$j]) : :(@inbounds resize(m.v[$j],i)) for j ∈ list(1,N)]...))
+end
+
+@generated Base.size(m::ProductTopology{N}) where N = :(($([:(size(@inbounds m.v[$i])...) for i ∈ 1:N]...),))
+@generated Base.getindex(m::ProductTopology{N},i::Vararg{Int}) where N = :(Values{N,Int}($([:((@inbounds m.v[$j])[@inbounds i[$j]]) for j ∈ 1:N]...)))
+Base.getindex(m::ProductTopology{1},i::Int) = Values(((@inbounds m.v[1])[i],))
 @pure Base.getindex(t::ProductTopology,i::CartesianIndex) = getindex(t,i.I...)
 @pure Base.eltype(::Type{<:ProductTopology{N}}) where N = Values{N,Int}
 
@@ -185,6 +203,10 @@ exclude(m::ProductTopology) = m
 exclude(m::ProductTopology{N},n::Int) where N = ProductTopology(m.v[vcat([i≠n ? Values((i,)) : Values{0,Int}() for i ∈ list(1,N)]...)])
 exclude(m::ProductTopology{N},n::Int...) where N = ProductTopology(m.v[vcat([i∉n ? Values((i,)) : Values{0,Int}() for i ∈ list(1,N)]...)])
 
+cross(a::ProductTopology,b::ProductTopology) = ProductTopology(Values(a.v...,b.v...))
+cross(a::ProductTopology,b::AbstractVector{Int}) = ProductTopology(Values(a.v...,b))
+cross(a::ProductTopology,b::Int) = a × OneTo(b)
+
 # QuotientTopology
 
 struct QuotientTopology{N,L,M,O,LA<:ImmersedTopology{L,L}} <: ImmersedTopology{N,N}
@@ -192,7 +214,7 @@ struct QuotientTopology{N,L,M,O,LA<:ImmersedTopology{L,L}} <: ImmersedTopology{N
     q::Values{O,LA}
     r::Values{M,Int}
     s::Values{N,Int}
-    #QuotientTopology(p::Values{O,P},q::Values{O,LA},r::Values{M,Int},n::Values{N,Int}) where {P,L,LA<:ImmersedTopology{L,L},M,N} = QuotientTopology{N,L,M,O,P,LA}(p,q,r,n)
+    #QuotientTopology(p::Values{O,Int},q::Values{O,LA},r::Values{M,Int},n::Values{N,Int}) where {O,L,LA<:ImmersedTopology{L,L},M,N} = QuotientTopology{N,L,M,O,LA}(p,q,r,n)
 end
 
 const OpenTopology{N,L,M,LA} = QuotientTopology{N,L,M,0,LA}
@@ -201,8 +223,8 @@ QuotientTopology(n::ProductTopology) = OpenTopology(n.v)
 OpenTopology(n::ProductTopology) = OpenTopology(n.v)
 OpenTopology(n::QuotientTopology) = OpenTopology(size(n))
 OpenTopology(n::Values{N,Int}) where N = QuotientTopology(Values{0,Int}(),Values{0,Array{Values{N-1,Int},N-1}}(),zeros(Values{2N,Int}),n)
-RibbonTopology(n::Values{2,Int}) = QuotientTopology(Values(4,3),Values(ProductTopology(n[1]),ProductTopology(n[1])),Values(0,0,1,2),n)
-MobiusTopology(n::Values{2,Int}) = QuotientTopology(Values(4,3),Values(ProductTopology(n[1]:-1:1),ProductTopology(n[1]:-1:1)),Values(0,0,1,2),n)
+RibbonTopology(n::Values{2,Int}) = QuotientTopology(Values(2,1),Values(ProductTopology(n[2]),ProductTopology(n[2])),Values(1,2,0,0),n)
+MobiusTopology(n::Values{2,Int}) = QuotientTopology(Values(2,1),Values(ProductTopology(n[2]:-1:1),ProductTopology(n[2]:-1:1)),Values(1,2,0,0),n)
 WingTopology(n::Values{2,Int}) = QuotientTopology(Values(1,2),Values(ProductTopology(n[2]:-1:1),ProductTopology(n[2]:-1:1)),Values(1,2,0,0),n)
 MirrorTopology(n::Values{1,Int}) = QuotientTopology(Values((1,)),Array{Values{0,Int},0}.(Values((undef,))),Values(1,0),n)
 MirrorTopology(n::Values{2,Int}) = QuotientTopology(Values((1,)),Values((ProductTopology(n[2]),)),Values(1,0,0,0),n)
@@ -221,7 +243,7 @@ TorusTopology(n::Values{4,Int}) = QuotientTopology(Values(2,1,4,3,6,5,8,7),Value
 TorusTopology(n::Values{5,Int}) = QuotientTopology(Values(2,1,4,3,6,5,8,7,10,9),Values(ProductTopology(n[2],n[3],n[4],n[5]),ProductTopology(n[2],n[3],n[4],n[5]),ProductTopology(n[1],n[3],n[4],n[5]),ProductTopology(n[1],n[3],n[4],n[5]),ProductTopology(n[1],n[2],n[4],n[5]),ProductTopology(n[1],n[2],n[4],n[5]),ProductTopology(n[1],n[2],n[3],n[5]),ProductTopology(n[1],n[2],n[3],n[5]),ProductTopology(n[1],n[2],n[3],n[4]),ProductTopology(n[1],n[2],n[3],n[4])),Values(1,2,3,4,5,6,7,8,9,10),n)
 HopfTopology(n::Values{2,Int}) = QuotientTopology(Values(2,1,4,3),Values(ProductTopology(CrossRange(n[2])),ProductTopology(CrossRange(n[2])),ProductTopology(n[1]),ProductTopology(n[1])),Values(1,2,3,4),n)
 HopfTopology(n::Values{3,Int}) = QuotientTopology(Values(2,1,4,3,6,5),Values(ProductTopology(OneTo(n[2]),CrossRange(n[3])),ProductTopology(OneTo(n[2]),CrossRange(n[3])),ProductTopology(OneTo(n[1]),CrossRange(n[3])),ProductTopology(OneTo(n[1]),CrossRange(n[3])),ProductTopology(n[1],n[2]),ProductTopology(n[1],n[2])),Values(1,2,3,4,5,6),n)
-KleinTopology(n::Values{2,Int}) = QuotientTopology(Values(2,1,4,3),Values(ProductTopology(n[2]),ProductTopology(n[2]),ProductTopology(n[1]:-1:1),ProductTopology(n[1]:-1:1)),Values(1,2,3,4),n)
+KleinTopology(n::Values{2,Int}) = QuotientTopology(Values(2,1,4,3),Values(ProductTopology(n[2]:-1:1),ProductTopology(n[2]:-1:1),ProductTopology(1:1:n[1]),ProductTopology(1:1:n[1])),Values(1,2,3,4),n)
 ConeTopology(n::Values{2,Int}) = QuotientTopology(Values(1,4,3),Values(ProductTopology(CrossRange(n[2])),ProductTopology(n[1]),ProductTopology(n[1])),Values(1,0,2,3),n)
 DrumTopology(n::Values{2,Int}) = QuotientTopology(Values(1,2,4,3),Values(ProductTopology(CrossRange(n[2])),ProductTopology(n[2]),ProductTopology(n[1]),ProductTopology(n[1])),Values(1,2,3,4),n)
 SphereTopology(n::Values{2,Int}) = QuotientTopology(Values(1,2,4,3),Values(ProductTopology(CrossRange(n[2])),ProductTopology(CrossRange(n[2])),ProductTopology(n[1]),ProductTopology(n[1])),Values(1,2,3,4),n)
@@ -295,6 +317,15 @@ iscompact(t::QuotientTopology) = false
 iscompact(t::CompactTopology) = true
 _to_axis(f::Int) = (iseven(f) ? f : f+1)÷2
 
+function LinearAlgebra.cross(m::QuotientTopology{1},n::Int)
+    QuotientTopology(m.p,
+        Values(ProductTopology(n),ProductTopology(n)),
+        Values(m.r...,0,0),Values(m.s...,n))
+end
+function LinearAlgebra.cross(m::QuotientTopology,n::Int)
+    QuotientTopology(m.p,m.q .× n,Values(m.r...,0,0),Values(m.s...,n))
+end
+
 getlocate(i) = Values((i,))
 getlocate(a,i) = Values((i,))
 getlocate(a,i,j) = isone(a) ? Values(i,j) : Values(j,i)
@@ -306,12 +337,12 @@ function locate_fast(pr1::Int,s,i::Int)
     if isodd(pr1)
         return abs(i-1)+1
     else
-        return s[_to_axis(pr1)]-abs(i-1)
+        return @inbounds s[_to_axis(pr1)]-abs(i-1)
     end
 end
 function locate_fast(pr2::Int,n::Int,s,i::Int)
     if iseven(pr2)
-        return s[_to_axis(pr2)]+n-i
+        return @inbounds s[_to_axis(pr2)]+n-i
     else
         return (i+1)-n
     end
@@ -320,30 +351,47 @@ function locate(pr1::Int,a::Int,s,i::Int)
     if isodd(pr1)
         return abs(i-1)+1
     else
-        return s[a]-abs(i-1)
+        return @inbounds s[a]-abs(i-1)
     end
 end
 function locate(pr2::Int,a::Int,n::Int,s,i::Int)
     if iseven(pr2)
-        return s[a]+n-i
+        return @inbounds s[a]+n-i
     else
         return (i+1)-n
     end
 end
 
 function location(p,q,r::Int,s::Tuple,i::Int,jk::Int...)
-    pr = p[r]
+    pr = @inbounds p[r]
     a = _to_axis(pr)
-    getlocate(a,locate(pr,a,s,i),q[r][jk...]...)
+    getlocate(a,locate(pr,a,s,i),(@inbounds q[r])[jk...]...)
 end
 function location(p,q,r::Int,n::Int,s::Tuple,i::Int,jk::Int...)
-    pr = p[r]
+    pr = @inbounds p[r]
     a = _to_axis(pr)
-    getlocate(a,locate(pr,a,n,s,i),q[r][jk...]...)
+    getlocate(a,locate(pr,a,n,s,i),(@inbounds q[r])[jk...]...)
+end
+
+@generated function resize(m::OpenTopology{N},i) where N
+    Expr(:call,:QuotientTopology,:(m.p),:(m.q),:(m.r),
+         Expr(:call,:Values,[j≠N ? :(@inbounds m.s[$j]) : :i for j ∈ list(1,N)]...))
+end
+@generated function resize(m::QuotientTopology{N,L,M,O},i) where {N,L,M,O}
+    Expr(:call,:QuotientTopology,:(m.p),
+         Expr(:call,:Values,Expr(:tuple,[:((@inbounds m.r[$j])∉(@inbounds m.r[2N-1],@inbounds m.r[2N]) ? (@inbounds resize(m.q[$j],i)) : (@inbounds m.q[$j])) for j ∈ list(1,O)]...)),
+         :(m.r),
+         Expr(:call,:Values,[j≠N ? :(@inbounds m.s[$j]) : :i for j ∈ list(1,N)]...))
+end
+@generated function resize(m::CompactTopology{N,L,M,O},i) where {N,L,M,O}
+    Expr(:call,:QuotientTopology,:(m.p),
+         Expr(:call,:Values,[j∉(O-1,O) ? :(@inbounds resize(m.q[$j],i)) : :(@inbounds m.q[$j]) for j ∈ list(1,O)]...),
+         :(m.r),
+         Expr(:call,:Values,[j≠N ? :(@inbounds m.s[$j]) : :i for j ∈ list(1,N)]...))
 end
 
 @pure Base.eltype(::Type{<:QuotientTopology{N}}) where N = Values{N,Int}
-Base.size(m::QuotientTopology) = Tuple(m.s)
+Base.size(m::QuotientTopology) = m.s.v
 Base.iterate(t::QuotientTopology) = (getindex(t,1),1)
 Base.iterate(t::QuotientTopology,state) = (s=state+1; s≤length(t) ? (getindex(t,s),s) : nothing)
 
@@ -351,15 +399,16 @@ function Base.getindex(m::QuotientTopology{N},i::Vararg{Int,N}) where N
     N > 5 ? Values{N,Int}(i) : getindex(m,Val(0),i...)
 end
 function Base.getindex(m::QuotientTopology{1},i::Int)
-    n = size(m)[1]
+    s = size(m)
+    n = @inbounds s[1]
     ii = if (i > 1 && i < n)
         i
     elseif i < 2
-        r = m.r[1]
-        iszero(r) ? i : locate_fast(m.p[r],size(m),i)
+        r = @inbounds m.r[1]
+        iszero(r) ? i : locate_fast((@inbounds m.p[r]),s,i)
     else
-        r = m.r[2]
-        iszero(r) ? i : locate_fast(m.p[r],n,size(m),i)
+        r = @inbounds m.r[2]
+        iszero(r) ? i : locate_fast((@inbounds m.p[r]),n,s,i)
     end
     return Values{1,Int}((ii,))
 end
@@ -370,22 +419,22 @@ Base.getindex(m::QuotientTopology{N},::Val,i::Vararg{Int,N}) where N = getindex(
 Base.getindex(m::QuotientTopology{1},::Val,i::Int) = getindex(m,i)
 function Base.getindex(m::QuotientTopology{2},N::Val,i::Int,j::Int)
     s = size(m)
-    n1,n2 = s[1],s[2]
+    n1,n2 = @inbounds (s[1],s[2])
     isi,isj = (bounds(i,n1,N,Val(1)),bounds(j,n2,N,Val(2)))
     if isj && !isi
         if i < 2
-            r = m.r[1]
+            r = @inbounds m.r[1]
             !iszero(r) && (return location(m.p,m.q,r,s,i,j))
         else
-            r = m.r[2]
+            r = @inbounds m.r[2]
             !iszero(r) && (return location(m.p,m.q,r,n1,s,i,j))
         end
     elseif isi && !isj
         if j < 2
-            r = m.r[3]
+            r = @inbounds m.r[3]
             !iszero(r) && (return location(m.p,m.q,r,s,j,i))
         else
-            r = m.r[4]
+            r = @inbounds m.r[4]
             !iszero(r) && (return location(m.p,m.q,r,n2,s,j,i))
         end
     end
@@ -393,30 +442,30 @@ function Base.getindex(m::QuotientTopology{2},N::Val,i::Int,j::Int)
 end
 function Base.getindex(m::QuotientTopology{3},N::Val,i::Int,j::Int,k::Int)
     s = size(m)
-    n1,n2,n3 = s[1],s[2],s[3]
+    n1,n2,n3 = @inbounds (s[1],s[2],s[3])
     isi,isj,isk = (bounds(i,n1,N,Val(1)),bounds(j,n2,N,Val(2)),bounds(k,n3,N,Val(3)))
     if isj && isk && !isi
         if i < 2
-            r = m.r[1]
+            r = @inbounds m.r[1]
             !iszero(r) && (return location(m.p,m.q,r,s,i,j,k))
         else
-            r = m.r[2]
+            r = @inbounds m.r[2]
             !iszero(r) && (return location(m.p,m.q,r,n1,s,i,j,k))
         end
     elseif isi && isk && !isj
         if j < 2
-            r = m.r[3]
+            r = @inbounds m.r[3]
             !iszero(r) && (return location(m.p,m.q,r,s,j,i,k))
         else
-            r = m.r[4]
+            r = @inbounds m.r[4]
             !iszero(r) && (return location(m.p,m.q,r,n2,s,j,i,k))
         end
     elseif isi && isj && !isk
         if k < 2
-            r = m.r[5]
+            r = @inbounds m.r[5]
             !iszero(r) && (return location(m.p,m.q,r,s,k,i,j))
         else
-            r = m.r[6]
+            r = @inbounds m.r[6]
             !iszero(r) && (return location(m.p,m.q,r,n3,s,k,i,j))
         end
     end
@@ -424,38 +473,38 @@ function Base.getindex(m::QuotientTopology{3},N::Val,i::Int,j::Int,k::Int)
 end
 function Base.getindex(m::QuotientTopology{4},N::Val,i::Int,j::Int,k::Int,l::Int)
     s = size(m)
-    n1,n2,n3,n4 = s[1],s[2],s[3],s[4]
+    n1,n2,n3,n4 = @inbounds (s[1],s[2],s[3],s[4])
     isi,isj,isk,isl = (bounds(i,n1,N,Val(1)),bounds(j,n2,N,Val(2)),bounds(k,n3,N,Val(3)),bounds(l,n4,N,Val(4)))
     if isj && isk && isl && !isi
         if i < 2
-            r = m.r[1]
+            r = @inbounds m.r[1]
             !iszero(r) && (return location(m.p,m.q,r,s,i,j,k,l))
         else
-            r = m.r[2]
+            r = @inbounds m.r[2]
             !iszero(r) && (return location(m.p,m.q,r,n1,s,i,j,k,l))
         end
     elseif isi && isk && isl && !isj
         if j < 2
-            r = m.r[3]
+            r = @inbounds m.r[3]
             !iszero(r) && (return location(m.p,m.q,r,s,j,i,k,l))
         else
-            r = m.r[4]
+            r = @inbounds m.r[4]
             !iszero(r) && (return location(m.p,m.q,r,n2,s,j,i,k,l))
         end
     elseif isi && isj && isl && !isk
         if k < 2
-            r = m.r[5]
+            r = @inbounds m.r[5]
             !iszero(r) && (return location(m.p,m.q,r,s,k,i,j,l))
         else
-            r = m.r[6]
+            r = @inbounds m.r[6]
             !iszero(r) && (return location(m.p,m.q,r,n3,s,k,i,j,l))
         end
     elseif isi && isj && isk && !isl
         if l < 2
-            r = m.r[7]
+            r = @inbounds m.r[7]
             !iszero(r) && (return location(m.p,m.q,r,s,l,i,j,k))
         else
-            r = m.r[8]
+            r = @inbounds m.r[8]
             !iszero(r) && (return location(m.p,m.q,r,n4,s,l,i,j,k))
         end
     end
@@ -463,46 +512,46 @@ function Base.getindex(m::QuotientTopology{4},N::Val,i::Int,j::Int,k::Int,l::Int
 end
 function Base.getindex(m::QuotientTopology{5},N::Val,i::Int,j::Int,k::Int,l::Int,o::Int)
     s = size(m)
-    n1,n2,n3,n4,n5 = s[1],s[2],s[3],s[4],s[5]
+    n1,n2,n3,n4,n5 = @inbounds (s[1],s[2],s[3],s[4],s[5])
     isi,isj,isk,isl,iso = (bounds(i,n1,N,Val(1)),bounds(j,n2,N,Val(2)),bounds(k,n3,N,Val(3)),bounds(l,n4,N,Val(4)),bounds(o,n5,N,Val(5)))
     if isj && isk && isl && iso && !isi
         if i < 2
-            r = m.r[1]
+            r = @inbounds m.r[1]
             !iszero(r) && (return location(m.p,m.q,r,s,i,j,k,l,o))
         else
-            r = m.r[2]
+            r = @inbounds m.r[2]
             !iszero(r) && (return location(m.p,m.q,r,n1,s,i,j,k,l,o))
         end
     elseif isi && isk && isl && iso && !isj
         if j < 2
-            r = m.r[3]
+            r = @inbounds m.r[3]
             !iszero(r) && (return location(m.p,m.q,r,s,j,i,k,l,o))
         else
-            r = m.r[4]
+            r = @inbounds m.r[4]
             !iszero(r) && (return location(m.p,m.q,r,n2,s,j,i,k,l,o))
         end
     elseif isi && isj && isl && iso && !isk
         if k < 2
-            r = m.r[5]
+            r = @inbounds m.r[5]
             !iszero(r) && (return location(m.p,m.q,r,s,k,i,j,l,o))
         else
-            r = m.r[6]
+            r = @inbounds m.r[6]
             !iszero(r) && (return location(m.p,m.q,r,n3,s,k,i,j,l,o))
         end
     elseif isi && isj && isk && iso && !isl
         if l < 2
-            r = m.r[7]
+            r = @inbounds m.r[7]
             !iszero(r) && (return location(m.p,m.q,r,s,l,i,j,k,o))
         else
-            r = m.r[8]
+            r = @inbounds m.r[8]
             !iszero(r) && (return location(m.p,m.q,r,n4,s,l,i,j,k,o))
         end
     elseif isi && isj && isk && isl && !iso
         if o < 2
-            r = m.r[9]
+            r = @inbounds m.r[9]
             !iszero(r) && (return location(m.p,m.q,r,s,o,i,j,k,l))
         else
-            r = m.r[10]
+            r = @inbounds m.r[10]
             !iszero(r) && (return location(m.p,m.q,r,n4,s,o,i,j,k,l))
         end
     end
@@ -510,17 +559,20 @@ function Base.getindex(m::QuotientTopology{5},N::Val,i::Int,j::Int,k::Int,l::Int
 end
 
 function findface(m,r,i,vals)
-    if iszero(r[i]) || m.p[r[i]] ∉ r || m.q[r[i]][vals...] ≠ vals
+    ri = r[i]
+    if iszero(ri) || m.p[ri] ∉ r || m.q[ri][vals...] ≠ vals
         0
     else
-        m.p[r[i]]≠r[1] ? 2 : 1
+        m.p[ri]≠(@inbounds r[1]) ? 2 : 1
     end
 end
 function findface(m,r,i,vals,ex...)
-    if iszero(r[i]) || m.p[r[i]] ∉ r || exclude(m.q[r[i]],ex...)[vals...] ≠ vals
+    ri = r[i]
+    if iszero(ri) || m.p[ri] ∉ r || exclude(m.q[ri],ex...)[vals...] ≠ vals
         0
     else
-        findfirst(x->x==m.p[r[i]],r)
+        pri = m.p[ri]
+        findfirst(x->x==pri,r)
     end
 end
 
@@ -557,16 +609,27 @@ function subtopology(m::QuotientTopology,i::NTuple{N,Int},::Colon,j::NTuple{M,In
         findface(m,r,3,vals,N1),findface(m,r,4,vals,N1))
     p1z,p2z,p3z,p4z = iszero.((p1,p2,p3,p4))
     pz = !(p1z)+!(p2z)+!(p3z)+!(p4z)
-    a = Values{pz,Int}((p1z ? () : (p1,))...,(p2z ? () : (p2,))...,(p3z ? () : (p3,))...,(p4z ? () : (p4,))...)
+    vpz = if iszero(pz)
+        Values{0}
+    elseif isone(pz)
+        Values{1}
+    elseif pz==2
+        Values{2}
+    elseif pz==3
+        Values{3}
+    else
+        Values{4}
+    end
+    a = iszero(pz) ? Values{0,Int}() : vpz((p1z ? () : (p1,))...,(p2z ? () : (p2,))...,(p3z ? () : (p3,))...,(p4z ? () : (p4,))...)
     b = if iszero(pz)
-        Values{pz,Array{Values{1,Int},1}}()
-    else; Values{pz}(
+        Values{0,Array{Values{1,Int},1}}()
+    else; vpz(
         (p1z ? () : (ProductTopology(m.q[r[1]].v[M1-1]),))...,
         (p2z ? () : (ProductTopology(m.q[r[2]].v[M1-1]),))...,
         (p3z ? () : (ProductTopology(m.q[r[3]].v[N1]),))...,
         (p4z ? () : (ProductTopology(m.q[r[4]].v[N1]),))...)
     end
-    c = Values(p1z ? 0 : 1,p2z ? 0 : !(p1z)+!(p2z),p3z ? 0 : !(p1z)+!(p2z)+!(p3z),p4z ? 0 : pz)
+    c = Values{4,Int}(p1z ? 0 : 1,p2z ? 0 : !(p1z)+!(p2z),p3z ? 0 : !(p1z)+!(p2z)+!(p3z),p4z ? 0 : pz)
     QuotientTopology(a,b,c,Values(s[N1],s[M1]))
 end
 function subtopology(m::QuotientTopology,i::NTuple{N,Int},::Colon,j::NTuple{M,Int},::Colon,k::NTuple{L,Int},::Colon,l::NTuple{O,Int} where O) where {N,M,L}
@@ -794,7 +857,6 @@ metricextensor(c::Coordinate) = fiber(c)
 metrictensor(c) = InducedMetric()
 metrictensor(c::Coordinate) = TensorOperator(fiber(c)[1])
 
-
 Base.getindex(s::Coordinate,i::Int...) = getindex(s.v.first,i...)
 Base.getindex(s::Coordinate,i::Integer...) = getindex(s.v.first,i...)
 
@@ -838,6 +900,10 @@ const Section = LocalTensor
 const ↦, domain, codomain = LocalTensor, base, fiber
 ↤(F,B) = B ↦ F
 
+localfiber(x) = x
+localfiber(x::LocalTensor) = fiber(x)
+localfiber(x::LocalSection) = fiber(x)
+
 @inline Base.:<<(a::LocalFiber,b::LocalFiber) = contraction(b,~a)
 @inline Base.:>>(a::LocalFiber,b::LocalFiber) = contraction(~a,b)
 @inline Base.:<(a::LocalFiber,b::LocalFiber) = contraction(b,a)
@@ -975,13 +1041,18 @@ function deletepointcloud!(P::Int)
     nothing
 end
 
+⊕(a::PointArray,b::AbstractVector{<:Real}) = PointArray(points(a)⊕b)
+cross(a::PointArray,b::AbstractVector{<:Real}) = a⊕b
+
 Base.size(m::PointArray) = size(m.dom)
 Base.firstindex(m::PointCloud) = 1
 Base.lastindex(m::PointCloud) = length(points(m))
 Base.length(m::PointCloud) = length(points(m))
-Base.resize!(m::PointCloud,i::Int) = (resize!(points(m),i),resize!(metricextensor(m),i))
+Base.resize!(m::PointCloud,i::Int) = ((resize!(points(m),i),resize!(metricextensor(m),i)); m)
 Base.broadcast(f,t::PointArray) = PointArray(f.(points(t)),f.(metricextensor(t)))
 Base.broadcast(f,t::PointCloud) = PointCloud(f.(points(t)),f.(metricextensor(t)))
+resize_lastdim!(m::Global,i) = m
+resize_lastdim!(m::PointArray,i) = ((resize_lastdim!(m.dom,i),resize_lastdim!(m.cod,i)); m)
 
 function (m::PointArray)(i::Vararg{Union{Int,Colon}})
     pa = points(m)(i...)
@@ -1084,6 +1155,9 @@ abstract type GlobalFiber{E,N} <: FiberBundle{E,N} end
 Base.@pure isfiberbundle(::GlobalFiber) = true
 Base.@pure isfiberbundle(::Any) = false
 
+globalfiber(x) = x
+globalfiber(x::GlobalFiber) = fiber(x)
+
 topology(m::GlobalFiber) = topology(immersion(m))
 vertices(m::GlobalFiber) = vertices(immersion(m))
 iscover(m::GlobalFiber) = iscover(immersion(m))
@@ -1094,8 +1168,8 @@ arcdomain(t::GlobalFiber) = unitdomain(t)*arclength(codomain(t))
 graph(t::GlobalFiber) = graph.(t)
 
 Base.size(m::GlobalFiber) = size(m.cod)
-Base.resize!(m::GlobalFiber,i) = (resize!(domain(m),i),resize!(codomain(m),i))
-
+Base.resize!(m::GlobalFiber,i) = ((resize!(domain(m),i),resize!(codomain(m),i)); m)
+resize_lastdim!(m::GlobalFiber,i) = ((resize_lastdim!(domain(m),i),resize_lastdim!(codomain(m),i)); m)
 # AbstractFrameBundle
 
 export AbstractFrameBundle, GridFrameBundle, SimplexFrameBundle, FacetFrameBundle
@@ -1152,7 +1226,12 @@ basetype(::Type{<:GridFrameBundle{C}}) where C = basetype(C)
 fibertype(m::GridFrameBundle) = fibertype(base(m))
 fibertype(::Type{<:GridFrameBundle{C}}) where C = fibertype(C)
 
-Base.resize!(m::GridFrameBundle,i) = resize!(base(m),i)
+⊕(a::GridFrameBundle,b::AbstractVector{<:Real}) = GridFrameBundle(base(a)⊕b,immersion(a) × length(b))
+cross(a::GridFrameBundle,b::AbstractVector{<:Real}) = a⊕b
+
+resize_lastdim!(m::GridFrameBundle,i) = (resize_lastdim!(base(m),i); m)
+resize(m::GridFrameBundle) = GridFrameBundle(m.id,m.dom,resize(m.cod,size(m.dom)[end]))
+Base.resize!(m::GridFrameBundle,i) = (resize!(base(m),i); m)
 Base.broadcast(f,t::GridFrameBundle) = GridFrameBundle(f.(base(t)))
 
 (m::GridFrameBundle)(i::ImmersedTopology) = GridFrameBundle(bundle(m),base(m),i)