diff --git a/Project.toml b/Project.toml
index 3bd935d..80683af 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Cartan"
 uuid = "e24af43f-9fd5-4672-9264-a75f3ae40eb2"
 authors = ["Michael Reed"]
-version = "0.3.4"
+version = "0.3.5"
 
 [deps]
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
diff --git a/src/Cartan.jl b/src/Cartan.jl
index 40266a3..0bd43bd 100644
--- a/src/Cartan.jl
+++ b/src/Cartan.jl
@@ -71,7 +71,7 @@ end
 #    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)
+    TensorField(GridFrameBundle(id,PointArray(0,dom,met)),cod)
 end
 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)
@@ -242,7 +242,7 @@ find_gs(::Any, rest) = find_gs(rest)
 
 (m::TensorField{B,F,N,<:SimplexFrameBundle} where {B,F,N})(i::ImmersedTopology) = TensorField(PointCloud(m)(i),fiber(m)[vertices(i)])
 (m::TensorField{B,F,N,<:GridFrameBundle} where {B,F,N})(i::ImmersedTopology) = TensorField(base(m)(i),fiber(m))
-for fun ∈ (:Open,:Mirror,:Clamped,:Torus,:Ribbon,:Wing,:Mobius,:Klein,:Cone,:Drum,:Sphere,:Geographic,:Hopf)
+for fun ∈ (:Open,:Mirror,:Clamped,:Torus,:Ribbon,:Wing,:Mobius,:Klein,:Cone,:Polar,:Sphere,:Geographic,:Hopf)
     for top ∈ (Symbol(fun,:Topology),)
         @eval begin
             $top(m::TensorField{B,F,N,<:GridFrameBundle} where {B,F,N}) = TensorField($top(base(m)),fiber(m))
@@ -278,7 +278,13 @@ 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,t1 = points(m),(@inbounds t[1])
+(m::Grid{1})(t::Chain) = linterp(m,t)
+(m::Grid{1})(t::AbstractFloat) = linterp(m,t)
+(m::IntervalMap)(t::Chain) = linterp(m,t)
+(m::IntervalMap)(t::AbstractFloat) = linterp(m,t)
+function linterp(m,t)
+    p,f,t1 = points(m),fiber(m),(@inbounds t[1])
+    isnan(t1) && (return zero(fibertype(m))/0)
     i = searchsortedfirst(p,@inbounds t1)-1
     if !isopen(m)
         q = immersion(m)
@@ -298,11 +304,11 @@ function (m::IntervalMap)(t); p,t1 = points(m),(@inbounds t[1])
     elseif iszero(i) || i==length(p)
         return zero(fibertype(m))
     end
-    linterp(t1,p[i],p[i+1],m.cod[i],m.cod[i+1])
+    linterp(t1,p[i],p[i+1],f[i],f[i+1])
 end
-function (m::IntervalMap)(t::Vector,d=diff(m.cod)./diff(m.dom))
+#=function (m::IntervalMap)(t::Vector,d=diff(m.cod)./diff(m.dom))
     [parametric(i,m,d) for i ∈ t]
-end
+end=#
 function parametric(t,m,d=diff(codomain(m))./diff(domain(m)))
     p = points(m)
     i = searchsortedfirst(p,t)-1
@@ -333,17 +339,26 @@ function (m::TensorField{B,F,2,<:FiberProductBundle} where {B,F})(t::Real)
 end
 
 #(m::TensorField)(t::TensorField) = TensorField(base(t),m.(fiber(t)))
-#(m::GridFrameBundle)(t::TensorField) = GridFrameBundle(base(t),TensorField(base(m),fiber(m)).(fiber(t)))
+#(m::GridFrameBundle)(t::TensorField) = GridFrameBundle(PointArray(points(t),m.(fiber(t))),immersion(m))
+(X::VectorField{B,F,N} where {B,F})(Y::VectorField{B,<:Chain{V,1,T,N} where {V,T},1} where B) where N = TensorField(base(Y),X.(fiber(Y)))
+(m::GridFrameBundle{N})(t::VectorField{B,<:Chain{V,1,T,N} where {V,T},1} where B) where N = GridFrameBundle(PointArray(points(t),m.(fiber(t))),immersion(m))
 
-function (m::TensorField{B,F,N,<:SimplexFrameBundle} where {B,N})(t) where F
-    j = findfirst(t,domain(m)); iszero(j) && (return zero(F))
+(m::SimplexFrameBundle)(t::Chain) = sinterp(m,t)
+(m::TensorField{B,F,N,<:SimplexFrameBundle} where {B,F,N})(t::Chain) = sinterp(m,t)
+function sinterp(m,t::Chain{V}) where V
+    j = findfirst(t,coordinates(m))
+    iszero(j) && (return zero(fibertype(m)))
     i = immersion(m)[j]
-    Chain(codomain(m)[i])⋅(Chain(points(domain(m))[i])\t)
+    Chain{V}(fiber(m)[i])⋅(Chain{V}(points(m)[i])\t)
 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,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])
+(m::Grid{2})(t::Chain) = bilinterp(m,t)
+(m::Grid{2})(x::AbstractFloat,y::AbstractFloat) = bilinterp(m,Chain(x,y))
+(m::TensorField{B,F,N,<:RealSpace{2}} where {B,F,N})(t::Chain) = bilinterp(m,t)
+(m::TensorField{B,F,N,<:RealSpace{2}} where {B,F,N})(x,y) = bilinterp(m,Chain(x,y))
+function bilinterp(m,t::Chain{V,G,T,2} where {G,T}) where V
+    x,y,f,t1,t2 = @inbounds (points(m).v[1],points(m).v[2],fiber(m),t[1],t[2])
+    (isnan(t1) || isnan(t2)) && (return zero(fibertype(m))/0)
     i,j = searchsortedfirst(x,t1)-1,searchsortedfirst(y,t2)-1
     if !isopen(m)
         q = immersion(m)
@@ -351,7 +366,7 @@ function (m::TensorField{B,F,N,<:RealSpace{2}} where {B,N})(t::Chain{V,G,T,2} wh
             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)
+            return zero(fibertype(m))
         else
             i1,i2,j1,j2 = (
                 i01 && !iq1,i02 && !iq2,
@@ -363,18 +378,23 @@ function (m::TensorField{B,F,N,<:RealSpace{2}} where {B,N})(t::Chain{V,G,T,2} wh
             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 ???
+        # elseif condition creates 1 allocation, as opposed to 0 ???
+        return zero(fibertype(m))
     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])
+    #f1 = linterp(t[1],x[i],x[i+1],f[i,j],f[i+1,j])
+    #f2 = linterp(t[1],x[i],x[i+1],f[i,j+1],f[i+1,j+1])
     #linterp(t[2],y[j],y[j+1],f1,f2)
     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])
+        f[i,j],f[i+1,j],f[i,j+1],f[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,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])
+(m::Grid{3})(t::Chain) = trilinterp(m,t)
+(m::Grid{3})(x::AbstractFloat,y::AbstractFloat,z::AbstractFloat) = trilinterp(m,Chain(x,y,z))
+(m::TensorField{B,F,N,<:RealSpace{3}} where {B,F,N})(t::Chain) = trilinterp(m,t)
+(m::TensorField{B,F,N,<:RealSpace{3}} where {B,F,N})(x,y,z) = trilinterp(m,Chain(x,y,z))
+function trilinterp(m,t::Chain{V,G,T,3} where {G,T}) where V
+    x,y,z,f,t1,t2,t3 = @inbounds (points(m).v[1],points(m).v[2],points(m).v[3],fiber(m),t[1],t[2],t[3])
+    (isnan(t1) || isnan(t2) || isnan(t3)) && (return zero(fibertype(m))/0)
     i,j,k = (
         searchsortedfirst(x,t1)-1,
         searchsortedfirst(y,t2)-1,
@@ -386,7 +406,7 @@ function (m::TensorField{B,F,N,<:RealSpace{3}} where {B,N})(t::Chain{V,G,T,3} wh
             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)
+            return zero(fibertype(m))
         else
             i1,i2,j1,j2,k1,k2 = (
                 i01 && !iq1,i02 && !iq2,
@@ -400,23 +420,28 @@ function (m::TensorField{B,F,N,<:RealSpace{3}} where {B,N})(t::Chain{V,G,T,3} wh
             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 ???
+        # elseif condition creates 1 allocation, as opposed to 0 ???
+        return zero(fibertype(m))
     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])
+    #f1 = linterp(t[1],x[i],x[i+1],f[i,j,k],f[i+1,j,k])
+    #f2 = linterp(t[1],x[i],x[i+1],f[i,j+1,k],f[i+1,j+1,k])
     #g1 = linterp(t[2],y[j],y[j+1],f1,f2)
-    #f3 = linterp(t[1],x[i],x[i+1],m.cod[i,j,k+1],m.cod[i+1,j,k+1])
-    #f4 = linterp(t[1],x[i],x[i+1],m.cod[i,j+1,k+1],m.cod[i+1,j+1,k+1])
+    #f3 = linterp(t[1],x[i],x[i+1],f[i,j,k+1],f[i+1,j,k+1])
+    #f4 = linterp(t[1],x[i],x[i+1],f[i,j+1,k+1],f[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(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])
+        f[i,j,k],f[i+1,j,k],f[i,j+1,k],f[i+1,j+1,k],
+        f[i,j,k+1],f[i+1,j,k+1],f[i,j+1,k+1],f[i+1,j+1,k+1])
 end
 
+(m::Grid{4})(t::Chain) = quadlinterp(m,t)
+(m::Grid{4})(x::AbstractFloat,y::AbstractFloat,z::AbstractFloat,w::AbstractFloat) = quadlinterp(m,Chain(x,y,z,w))
+(m::TensorField{B,F,N,<:RealSpace{4}} where {B,F,N})(t::Chain) = quadlinterp(m,t)
 (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,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])
+function (m)(t::Chain{V,G,T,4} where {G,T}) where V
+    x,y,z,w,f,t1,t2,t3,t4 = @inbounds (points(m).v[1],points(m).v[2],points(m).v[3],points(m).v[4],fiber(m),t[1],t[2],t[3],t[4])
+    (isnan(t1) || isnan(t2) || isnan(t3) ||isnan(t4)) && (return zero(fibertype(m))/0)
     i,j,k,l = (
         searchsortedfirst(x,t1)-1,
         searchsortedfirst(y,t2)-1,
@@ -430,7 +455,7 @@ function (m::TensorField{B,F,N,<:RealSpace{4}} where {B,N})(t::Chain{V,G,T,4} wh
             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)
+            return zero(fibertype(m))
         else
             i1,i2,j1,j2,k1,k2,l1,l2 = (
                 i01 && !iq1,i02 && !iq2,
@@ -446,33 +471,38 @@ function (m::TensorField{B,F,N,<:RealSpace{4}} where {B,N})(t::Chain{V,G,T,4} wh
             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 ???
+        # elseif condition creates 1 allocation, as opposed to 0 ???
+        return zero(fibertype(m))
     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])
+    #f1 = linterp(t[1],x[i],x[i+1],f[i,j,k,l],f[i+1,j,k,l])
+    #f2 = linterp(t[1],x[i],x[i+1],f[i,j+1,k,l],f[i+1,j+1,k,l])
     #g1 = linterp(t[2],y[j],y[j+1],f1,f2)
-    #f3 = linterp(t[1],x[i],x[i+1],m.cod[i,j,k+1,l],m.cod[i+1,j,k+1,l])
-    #f4 = linterp(t[1],x[i],x[i+1],m.cod[i,j+1,k+1,l],m.cod[i+1,j+1,k+1,l])
+    #f3 = linterp(t[1],x[i],x[i+1],f[i,j,k+1,l],f[i+1,j,k+1,l])
+    #f4 = linterp(t[1],x[i],x[i+1],f[i,j+1,k+1,l],f[i+1,j+1,k+1,l])
     #g2 = linterp(t[2],y[j],y[j+1],f3,f4)
     #h1 = linterp(t[3],z[k],z[k+1],g1,g2)
-    #f5 = linterp(t[1],x[i],x[i+1],m.cod[i,j,k,l+1],m.cod[i+1,j,k,l+1])
-    #f6 = linterp(t[1],x[i],x[i+1],m.cod[i,j+1,k,l+1],m.cod[i+1,j+1,k,l+1])
+    #f5 = linterp(t[1],x[i],x[i+1],f[i,j,k,l+1],f[i+1,j,k,l+1])
+    #f6 = linterp(t[1],x[i],x[i+1],f[i,j+1,k,l+1],f[i+1,j+1,k,l+1])
     #g3 = linterp(t[2],y[j],y[j+1],f5,f6)
-    #f7 = linterp(t[1],x[i],x[i+1],m.cod[i,j,k+1,l+1],m.cod[i+1,j,k+1,l+1])
-    #f8 = linterp(t[1],x[i],x[i+1],m.cod[i,j+1,k+1,l+1],m.cod[i+1,j+1,k+1,l+1])
+    #f7 = linterp(t[1],x[i],x[i+1],f[i,j,k+1,l+1],f[i+1,j,k+1,l+1])
+    #f8 = linterp(t[1],x[i],x[i+1],f[i,j+1,k+1,l+1],f[i+1,j+1,k+1,l+1])
     #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(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],
-        m.cod[i,j,k+1,l+1],m.cod[i+1,j,k+1,l+1],m.cod[i,j+1,k+1,l+1],m.cod[i+1,j+1,k+1,l+1])
+        f[i,j,k,l],f[i+1,j,k,l],f[i,j+1,k,l],f[i+1,j+1,k,l],
+        f[i,j,k+1,l],f[i+1,j,k+1,l],f[i,j+1,k+1,l],f[i+1,j+1,k+1,l],
+        f[i,j,k,l+1],f[i+1,j,k,l+1],f[i,j+1,k,l+1],f[i+1,j+1,k,l+1],
+        f[i,j,k+1,l+1],f[i+1,j,k+1,l+1],f[i,j+1,k+1,l+1],f[i+1,j+1,k+1,l+1])
 end
 
+(m::Grid{5})(t::Chain) = quintlinterp(m,t)
+(m::Grid{5})(x::AbstractFloat,y::AbstractFloat,z::AbstractFloat,w::AbstractFloat,v::AbstractFloat) = quintlinterp(m,Chain(x,y,z,w,v))
+(m::TensorField{B,F,N,<:RealSpace{5}} where {B,F,N})(t::Chain) = quintlinterp(m,t)
 (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,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])
+function quintlinterp(m,t::Chain{V,G,T,5} where {G,T}) where V
+    x,y,z,w,v,f,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],fiber(m),t[1],t[2],t[3],t[4],t[5])
+    (isnan(t1) || isnan(t2) || isnan(t3) || isnan(t4) || isnan(t5)) && (return zero(fibertype(m))/0)
     i,j,k,l,o = (
         searchsortedfirst(x,t1)-1,
         searchsortedfirst(y,t2)-1,
@@ -488,7 +518,7 @@ function (m::TensorField{B,F,N,<:RealSpace{5}} where {B,N})(t::Chain{V,G,T,5} wh
             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)
+            return zero(fibertype(m))
         else
             i1,i2,j1,j2,k1,k2,l1,l2,o1,o2 = (
                 i01 && !iq1,i02 && !iq2,
@@ -506,17 +536,18 @@ function (m::TensorField{B,F,N,<:RealSpace{5}} where {B,N})(t::Chain{V,G,T,5} wh
             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 ???
+        # elseif condition creates 1 allocation, as opposed to 0 ???
+        return zero(fibertype(m))
     end
     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],
-        m.cod[i,j,k+1,l+1,o],m.cod[i+1,j,k+1,l+1,o],m.cod[i,j+1,k+1,l+1,o],m.cod[i+1,j+1,k+1,l+1,o],
-        m.cod[i,j,k,l,o+1],m.cod[i+1,j,k,l,o+1],m.cod[i,j+1,k,l,o+1],m.cod[i+1,j+1,k,l,o+1],
-        m.cod[i,j,k+1,l,o+1],m.cod[i+1,j,k+1,l,o+1],m.cod[i,j+1,k+1,l,o+1],m.cod[i+1,j+1,k+1,l,o+1],
-        m.cod[i,j,k,l+1,o+1],m.cod[i+1,j,k,l+1,o+1],m.cod[i,j+1,k,l+1,o+1],m.cod[i+1,j+1,k,l+1,o+1],
-        m.cod[i,j,k+1,l+1,o+1],m.cod[i+1,j,k+1,l+1,o+1],m.cod[i,j+1,k+1,l+1,o+1],m.cod[i+1,j+1,k+1,l+1,o+1])
+        f[i,j,k,l,o],f[i+1,j,k,l,o],f[i,j+1,k,l,o],f[i+1,j+1,k,l,o],
+        f[i,j,k+1,l,o],f[i+1,j,k+1,l,o],f[i,j+1,k+1,l,o],f[i+1,j+1,k+1,l,o],
+        f[i,j,k,l+1,o],f[i+1,j,k,l+1,o],f[i,j+1,k,l+1,o],f[i+1,j+1,k,l+1,o],
+        f[i,j,k+1,l+1,o],f[i+1,j,k+1,l+1,o],f[i,j+1,k+1,l+1,o],f[i+1,j+1,k+1,l+1,o],
+        f[i,j,k,l,o+1],f[i+1,j,k,l,o+1],f[i,j+1,k,l,o+1],f[i+1,j+1,k,l,o+1],
+        f[i,j,k+1,l,o+1],f[i+1,j,k+1,l,o+1],f[i,j+1,k+1,l,o+1],f[i+1,j+1,k+1,l,o+1],
+        f[i,j,k,l+1,o+1],f[i+1,j,k,l+1,o+1],f[i,j+1,k,l+1,o+1],f[i+1,j+1,k,l+1,o+1],
+        f[i,j,k+1,l+1,o+1],f[i+1,j,k+1,l+1,o+1],f[i,j+1,k+1,l+1,o+1],f[i+1,j+1,k+1,l+1,o+1])
 end
 
 valmat(t::Values{N,<:Vector},s=size(t[1])) where N = [Values((t[q][i] for q ∈ OneTo(N))...) for i ∈ OneTo(s[1])]
@@ -544,6 +575,18 @@ end
 @inline Base.:<<(a::GlobalFiber,b::GlobalFiber) = contraction(b,~a)
 @inline Base.:>>(a::GlobalFiber,b::GlobalFiber) = contraction(~a,b)
 @inline Base.:<(a::GlobalFiber,b::GlobalFiber) = contraction(b,a)
+Base.sign(a::TensorField) = TensorField(base(a), sign.(fiber(Real(a))))
+Base.inv(a::TensorField{B,<:Real} where B) = TensorField(base(a), inv.(fiber(a)))
+Base.inv(a::TensorField{B,<:Complex} where B) = TensorField(base(a), inv.(fiber(a)))
+Base.:/(a::TensorField,b::TensorField{B,<:Real} where B) = TensorField(base(a), fiber(a)./fiber(b))
+Base.:/(a::TensorField,b::TensorField{B,<:Complex} where B) = TensorField(base(a), fiber(a)./fiber(b))
+Grassmann.eigen(t::TensorField,i::Val) = TensorField(base(t), eigen.(fiber(t),i))
+Grassmann.eigen(t::TensorField,i::Int) = TensorField(base(t), eigen.(fiber(t),Val(i)))
+Grassmann.eigvals(t::TensorField,i::Val) = TensorField(base(t), eigvals.(fiber(t),i))
+Grassmann.eigvals(t::TensorField,i::Int) = TensorField(base(t), eigvals.(fiber(t),Val(i)))
+Grassmann.eigvecs(t::TensorField,i::Val) = TensorField(base(t), eigvecs.(fiber(t),i))
+Grassmann.eigvecs(t::TensorField,i::Int) = TensorField(base(t), eigvecs.(fiber(t),Val(i)))
+Grassmann.eigpolys(t::TensorField,G::Val) = TensorField(base(t), eigpolys.(fiber(t),G))
 for type ∈ (:TensorField,)
     for (op,mop) ∈ ((:*,:wedgedot_metric),(:wedgedot,:wedgedot_metric),(:veedot,:veedot_metric),(:⋅,:contraction_metric),(:contraction,:contraction_metric),(:>,:contraction_metric),(:⊘,:⊘),(:>>>,:>>>),(:/,:/),(:^,:^))
         let bop = op ∈ (:*,:>,:>>>,:/,:^) ? :(Base.$op) : :(Grassmann.$op)
@@ -560,10 +603,10 @@ end
 for fun ∈ (:exp,:exp2,:exp10,:log,:log2,:log10,:sinh,:cosh,:abs,:sqrt,:cbrt,:cos,:sin,:tan,:cot,:sec,:csc,:asec,:acsc,:sech,:csch,:asech,:tanh,:coth,:asinh,:acosh,:atanh,:acoth,:asin,:acos,:atan,:acot,:sinc,:cosc,:cis,:abs2,:inv)
     @eval Base.$fun(t::TensorField) = TensorField(domain(t), $fun.(codomain(t),ref(metricextensor(t))))
 end
-for fun ∈ (:reverse,:involute,:clifford,:even,:odd,:scalar,:vector,:bivector,:volume,:value,:complementleft,:realvalue,:imagvalue,:outermorphism,:Outermorphism,:DiagonalOperator,:TensorOperator)
+for fun ∈ (:reverse,:involute,:clifford,:even,:odd,:scalar,:vector,:bivector,:volume,:value,:complementleft,:realvalue,:imagvalue,:outermorphism,:Outermorphism,:DiagonalOperator,:TensorOperator,:eigen,:eigvecs,:eigvals,:eigvalsreal,:eigvalscomplex,:eigvecsreal,:eigvecscomplex,:eigpolys)
     @eval Grassmann.$fun(t::TensorField) = TensorField(domain(t), $fun.(codomain(t)))
 end
-for fun ∈ (:⋆,:angle,:radius,:complementlefthodge,:pseudoabs,:pseudoabs2,:pseudoexp,:pseudolog,:pseudoinv,:pseudosqrt,:pseudocbrt,:pseudocos,:pseudosin,:pseudotan,:pseudocosh,:pseudosinh,:pseudotanh,:metric)
+for fun ∈ (:⋆,:angle,:radius,:complementlefthodge,:pseudoabs,:pseudoabs2,:pseudoexp,:pseudolog,:pseudoinv,:pseudosqrt,:pseudocbrt,:pseudocos,:pseudosin,:pseudotan,:pseudocosh,:pseudosinh,:pseudotanh,:metric,:unit)
     @eval Grassmann.$fun(t::TensorField) = TensorField(domain(t), $fun.(codomain(t),ref(metricextensor(t))))
 end
 for fun ∈ (:sum,:prod)
@@ -723,46 +766,62 @@ function __init__()
                 Makie.lines!(x,Real.(getindex.(y,i));args...)
             end
         end
-        function Makie.lines(t::IntervalMap{B,<:Endomorphism},f::Function=speed;args...) where B<:Coordinate{<:AbstractReal}
+        function Makie.lines(t::IntervalMap{B,<:TensorOperator},f::Function=speed;args...) where B<:Coordinate{<:AbstractReal}
             display(Makie.lines(getindex.(t,1),f;args...))
             for i ∈ 2:mdims(eltype(codomain(t)))
                 Makie.lines!(getindex.(t,i),f;args...)
             end
         end
-        function Makie.lines!(t::IntervalMap{B,<:Endomorphism},f::Function=speed;args...) where B<:Coordinate{<:AbstractReal}
+        function Makie.lines!(t::IntervalMap{B,<:TensorOperator},f::Function=speed;args...) where B<:Coordinate{<:AbstractReal}
             display(Makie.lines!(getindex.(t,1),f;args...))
             for i ∈ 2:mdims(eltype(codomain(t)))
                 Makie.lines!(getindex.(t,i),f;args...)
             end
         end
-        function Makie.arrows(M::VectorField,t::TensorField{B,<:Endomorphism,N,<:GridFrameBundle} where B;args...) where N
+        function Makie.arrows(M::VectorField,t::TensorField{B,<:TensorOperator,N,<:GridFrameBundle} where B;args...) where N
             Makie.arrows(TensorField(fiber(M),fiber(t));args...)
         end
-        function Makie.arrows!(M::VectorField,t::TensorField{B,<:Endomorphism,N,<:GridFrameBundle} where B;args...) where N
+        function Makie.arrows!(M::VectorField,t::TensorField{B,<:TensorOperator,N,<:GridFrameBundle} where B;args...) where N
             Makie.arrows!(TensorField(fiber(M),fiber(t)))
         end
-        function Makie.arrows(t::TensorField{<:Coordinate{<:Chain},<:Endomorphism,N,<:GridFrameBundle};args...) where N
+        function Makie.arrows(t::VectorField{<:Coordinate{<:Chain},<:TensorOperator,2,<:RealSpace{2}};args...)
             display(Makie.arrows(getindex.(t,1);args...))
             for i ∈ 2:mdims(eltype(codomain(t)))
                 Makie.arrows!(getindex.(t,i);args...)
             end
         end
-        function Makie.arrows!(t::TensorField{<:Coordinate{<:Chain},<:Endomorphism,N,<:GridFrameBundle};args...) where N
+        function Makie.arrows!(t::VectorField{<:Coordinate{<:Chain},<:TensorOperator,2,<:RealSpace{2}};args...)
             display(Makie.arrows!(getindex.(t,1);args...))
             for i ∈ 2:mdims(eltype(codomain(t)))
                 Makie.arrows!(getindex.(t,i);args...)
             end
         end
-        function Makie.arrows(t::VectorField{<:Coordinate{<:Chain},<:Endomorphism,2,<:RealSpace{2}};args...)
-            display(Makie.arrows(getindex.(t,1);args...))
+        for (fun,fun!) ∈ ((:arrows,:arrows!),(:streamplot,:streamplot!))
+            @eval begin
+                function Makie.$fun(t::TensorField{<:Coordinate{<:Chain},<:TensorOperator,N,<:GridFrameBundle};args...) where N
+                    display(Makie.$fun(getindex.(t,1);args...))
+                    for i ∈ 2:mdims(eltype(codomain(t)))
+                        Makie.$fun!(getindex.(t,i);args...)
+                    end
+                end
+                function Makie.$fun!(t::TensorField{<:Coordinate{<:Chain},<:TensorOperator,N,<:GridFrameBundle};args...) where N
+                    display(Makie.$fun!(getindex.(t,1);args...))
+                    for i ∈ 2:mdims(eltype(codomain(t)))
+                        Makie.$fun!(getindex.(t,i);args...)
+                    end
+                end
+            end
+        end
+        function Makie.streamplot(t::TensorField{<:Coordinate{<:Chain},<:TensorOperator,2,<:RealSpace{2}},m::U;args...) where U<:Union{<:VectorField,<:Function}
+            display(Makie.streamplot(getindex.(t,1),m;args...))
             for i ∈ 2:mdims(eltype(codomain(t)))
-                Makie.arrows!(getindex.(t,i);args...)
+                Makie.streamplot!(getindex.(t,i),m;args...)
             end
         end
-        function Makie.arrows!(t::VectorField{<:Coordinate{<:Chain},<:Endomorphism,2,<:RealSpace{2}};args...)
-            display(Makie.arrows!(getindex.(t,1);args...))
+        function Makie.streamplot!(t::TensorField{<:Coordinate{<:Chain},<:TensorOperator,2,<:RealSpace{2}},m::U;args...) where U<:Union{<:VectorField,<:Function}
+            display(Makie.streamplot!(getindex.(t,1),m;args...))
             for i ∈ 2:mdims(eltype(codomain(t)))
-                Makie.arrows!(getindex.(t,i);args...)
+                Makie.streamplot!(getindex.(t,i),m;args...)
             end
         end
         Makie.volume(t::VolumeGrid;args...) = Makie.volume(points(t).v...,Real.(codomain(t));args...)
@@ -807,6 +866,8 @@ function __init__()
         end
         Makie.wireframe(t::SurfaceGrid;args...) = Makie.wireframe(graph(t);args...)
         Makie.wireframe!(t::SurfaceGrid;args...) = Makie.wireframe!(graph(t);args...)
+        point2chain(x,V=Submanifold(2)) = Chain(x[1],x[2])
+        chain3vec(x) = Makie.Vec3(x[1],x[2],x[3])
         for fun ∈ (:streamplot,:streamplot!)
             @eval begin
                 Makie.$fun(f::Function,t::Rectangle;args...) = Makie.$fun(f,t.v...;args...)
@@ -815,13 +876,32 @@ function __init__()
                 Makie.$fun(m::ScalarMap,dims...;args...) = Makie.$fun(gradient_fast(m),dims...;args...)
                 Makie.$fun(m::VectorField{R,F,1,<:SimplexFrameBundle} where {R,F},dims...;args...) = Makie.$fun(p->Makie.Point(m(Chain(one(eltype(p)),p.data...))),dims...;args...)
                 Makie.$fun(m::VectorField{<:Coordinate{<:Chain},F,N,<:RealSpace} where {F,N};args...) = Makie.$fun(p->Makie.Point(m(Chain(p.data...))),points(m).v...;args...)
+                function Makie.$fun(M::VectorField,m::VectorField{<:Coordinate{<:Chain{V}},<:Chain,2,<:RealSpace{2}};args...) where V
+                    kwargs = if haskey(args,:gridsize)
+                        wargs = Dict(args)
+                        delete!(wargs,:gridsize)
+                        (;:gridsize => (args[:gridsize]...,1),wargs...)
+                    else
+                        pairs((;:gridsize => (32,32,1),args...))
+                    end
+                    w,gs = widths(points(m)),kwargs[:gridsize]
+                    scale = 0.2sqrt(surfacearea(M)/prod(w))
+                    st = Makie.$fun(p->(z=m(Chain{V}(p[1],p[2]));Makie.Point(z[1],z[2],0)),points(m).v...,Makie.ClosedInterval(-1e-15,1e-15);arrow_size=scale*minimum(w)/minimum((gs[1],gs[2])),kwargs...)
+                    $(fun≠:streamplot ? :pl : :((fig,ax,pl))) = st
+                    pl.transformation.transform_func[] = Makie.PointTrans{3}() do p
+                        return Makie.Point(M(Chain{V}(p[1],p[2])))
+                    end
+                    jac,arr = jacobian(M),pl.plots[2]
+                    arr.rotation[] = chain3vec.(jac.(point2chain.(arr.args[1][],V)).⋅point2chain.(arr.rotation[],V))
+                    return st
+                end
             end
         end
         for fun ∈ (:arrows,:arrows!)
             @eval begin
                 Makie.$fun(t::ScalarField{<:Coordinate{<:Chain},F,2,<:RealSpace{2}} where F;args...) = Makie.$fun(Makie.Point.(fiber(graph(Real(t))))[:],Makie.Point.(fiber(normal(Real(t))))[:];args...)
-                Makie.$fun(t::VectorField{<:Coordinate{<:Chain{W,L,F,2} where {W,L,F}},<:Chain{V,G,T,2} where {V,G,T},2,<:AlignedRegion{2}};args...) = Makie.$fun(domain(t).v...,getindex.(codomain(t),1),getindex.(codomain(t),2);args...)
-                Makie.$fun(t::VectorField{<:Coordinate{<:Chain},F,2,<:RealSpace{2}} where F;args...) = Makie.$fun(Makie.Point.(points(t))[:],Makie.Point.(codomain(t))[:];args...)
+                Makie.$fun(t::VectorField{<:Coordinate{<:Chain{W,L,F,2} where {W,L,F}},<:Chain{V,G,T,2} where {V,G,T},2,<:AlignedRegion{2}};args...) = Makie.$fun(points(t).v...,getindex.(codomain(t),1),getindex.(codomain(t),2);args...)
+                Makie.$fun(t::VectorField{<:Coordinate{<:Chain},<:Chain,2,<:RealSpace{2}};args...) = Makie.$fun(Makie.Point.(points(t))[:],Makie.Point.(codomain(t))[:];args...)
                 Makie.$fun(t::VectorField{<:Coordinate{<:Chain},F,N,<:GridFrameBundle} where {F,N};args...) = Makie.$fun(Makie.Point.(points(t))[:],Makie.Point.(codomain(t))[:];args...)
                 Makie.$fun(t::VectorField,f::VectorField;args...) = Makie.$fun(Makie.Point.(fiber(t))[:],Makie.Point.(fiber(f))[:];args...)
                 Makie.$fun(t::Rectangle,f::Function;args...) = Makie.$fun(t.v...,f;args...)
@@ -891,16 +971,16 @@ function __init__()
             Makie.mesh!(GridFrameBundle(fiber(M));args...)
         end
         function Makie.mesh(M::TensorField{B,<:AbstractReal,2,<:GridFrameBundle} where B;args...)
-            Makie.mesh(GeometryBasics.Mesh(base(M));color=fiber(M)[:],args...)
+            Makie.mesh(GeometryBasics.Mesh(base(M));color=fiber(Real(M))[:],args...)
         end
         function Makie.mesh!(M::TensorField{B,<:AbstractReal,2,<:GridFrameBundle} where B;args...)
-            Makie.mesh!(GeometryBasics.Mesh(base(M));color=fiber(M)[:],args...)
+            Makie.mesh!(GeometryBasics.Mesh(base(M));color=fiber(Real(M))[:],args...)
         end
         function Makie.mesh(M::TensorField{B,<:Chain,2,<:GridFrameBundle} where B,f::TensorField;args...)
-            Makie.mesh(GridFrameBundle(fiber(M));color=fiber(f)[:],args...)
+            Makie.mesh(GridFrameBundle(fiber(M));color=fiber(Real(f))[:],args...)
         end
         function Makie.mesh!(M::TensorField{B,<:Chain,2,<:GridFrameBundle} where B,f::TensorField;args...)
-            Makie.mesh!(GridFrameBundle(fiber(M));color=fiber(f)[:],args...)
+            Makie.mesh!(GridFrameBundle(fiber(M));color=fiber(Real(f))[:],args...)
         end
         Makie.wireframe(M::TensorField{B,<:Chain,N,<:GridFrameBundle} where {B,N};args...) = Makie.wireframe(GridFrameBundle(fiber(M));args...)
         Makie.wireframe!(M::TensorField{B,<:Chain,N,<:GridFrameBundle} where {B,N};args...) = Makie.wireframe!(GridFrameBundle(fiber(M));args...)
@@ -961,7 +1041,7 @@ function __init__()
     @require GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" begin
         Base.convert(::Type{GeometryBasics.Point},t::T) where T<:LocalFiber = GeometryBasics.Point(base(t))
         #GeometryBasics.Point(t::T) where T<:TensorAlgebra = convert(GeometryBasics.Point,t)
-        function GeometryBasics.Mesh(m::GridFrameBundle{C,2} where C)
+        function GeometryBasics.Mesh(m::GridFrameBundle{2})
             nm = size(points(m))
             faces = GeometryBasics.Tesselation(GeometryBasics.Rect(0, 0, 1, 1), nm)
             uv = Chain(0.0,0.0):map(inv,Chain((nm.-1)...)):Chain(1.0,1.0)
diff --git a/src/diffgeo.jl b/src/diffgeo.jl
index 02b9f01..18ded25 100644
--- a/src/diffgeo.jl
+++ b/src/diffgeo.jl
@@ -39,37 +39,6 @@ boundlog(t::TensorField,lim=10) = TensorField(base(t), boundlog.(codomain(t),lim
 isclosed(t::IntervalMap) = norm(codomain(t)[end]-codomain(t)[1]) ≈ 0
 updatetopology(t::IntervalMap) = isclosed(t) ? TorusTopology(t) : t
 
-export Grid
-
-const Grid{N,C,PA,TA} = GridFrameBundle{C,N,PA,TA}
-Grid(v::A,t::I=OpenTopology(size(v))) where {N,T,A<:AbstractArray{T,N},I} = GridFrameBundle(0,PointArray(0,v),t)
-
-#=struct Grid{N,T,A<:AbstractArray{T,N},I<:ImmersedTopology}
-    v::A
-    t::I
-    Grid(v::A,t::I=OpenTopology(size(v))) where {N,T,A<:AbstractArray{T,N},I} = new{N,T,A,I}(v,t)
-end=#
-
-#Grid(v::A) where {N,T,A<:AbstractArray{T,N}} = Grid(v,Global{N}(InducedMetric()))
-#Grid(v::GridFrameBundle{<:Real}) = Grid(points(v))
-#Grid(v::GridFrameBundle) = Grid(points(v),fiber(metricextensor(v)))
-
-#=Base.size(m::Grid) = size(m.v)
-
-@generated function Base.getindex(g::Grid{M,T,A,<:OpenTopology} where {T,A},j::Int,::Val{N},i::Vararg{Int}) where {M,N}
-    :(Base.getindex(g.v,$([k≠N ? :(i[$k]) : :(i[$k]+j) for k ∈ 1:M]...)))
-end
-@generated function Base.getindex(g::Grid{M},j::Int,::Val{N},i::Vararg{Int}) where {M,N}
-    :(Base.getindex(g.v,Base.getindex(g.t,$([k≠N ? :(i[$k]) : :(i[$k]+j) for k ∈ 1:M]...))...))
-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 ? :(@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 ? :(@inbounds i[$k]) : :(@inbounds i[$k]+j) for k ∈ 1:M]...))...))
-end
-
 # centraldiff
 
 centraldiffdiff(f,dt,l) = centraldiff(centraldiff(f,dt,l),dt,l)
@@ -187,9 +156,9 @@ for fun ∈ (:_slow,:_fast)
     cd,grad = Symbol(:centraldiff,fun),Symbol(:gradient,fun)
     cdg,cdp,cdf = Symbol(cd,:_calc),Symbol(cd,:_points),Symbol(cd,:_fiber)
     @eval begin
-        $cdf(f,args...) = $cdg(Grid(codomain(f),immersion(f)),args...)
-        $cdp(f,args...) = $cdg(Grid(points(f),immersion(f)),args...)
-        $cdp(f::TensorField{B,F,Nf,<:RealSpace{Nf,P,<:InducedMetric} where P} where {B,F,Nf},n::Val{N},args...) where N = $cd(points(f).v[N],args...)
+        $cdf(f,args...) = $cdg(GridFrameBundle(0,PointArray(0,fiber(f)),immersion(f)),args...)
+        $cdp(f,args...) = $cdg(GridFrameBundle(0,PointArray(0,points(f)),immersion(f)),args...)
+        $cdp(f::TensorField{B,F,Nf,<:RealSpace{Nf,P,<:InducedMetric} where P} where {B,F,Nf},n::Val{N},args...) where N = $cd(points(f).v[N],subtopology(immersion(f),n),args...)
         function $grad(f::IntervalMap,d::AbstractVector=$cdp(f))
             TensorField(domain(f), $cdf(f,d))
         end
@@ -214,249 +183,258 @@ for fun ∈ (:_slow,:_fast)
             TensorField(domain(f), $cdf(f,n,d))
         end
         $grad(f::TensorField,n::Int,args...) = $grad(f,Val(n),args...)
-        $cd(f::AbstractArray,args...) = $cdg(Grid(f),args...)
-        function $cdg(f::Grid{1},dt::Real,s::Tuple=size(f))
+        $cd(f::AbstractArray,args...) = $cdg(GridFrameBundle(0,PointArray(0,f)),args...)
+        $cd(f::AbstractArray,q::QuotientTopology,args...) = $cdg(GridFrameBundle(0,PointArray(0,f),q),args...)
+        function $cdg(f::GridFrameBundle{1},dt::Real,s::Tuple=size(f))
             d = similar(points(f))
             @threads for i ∈ OneTo(s[1])
                 d[i] = $cdg(f,s,i)/$cdg(i,dt,l)
             end
             return d
         end
-        function $cdg(f::Grid{1},dt::Vector,s::Tuple=size(f))
+        function $cdg(f::GridFrameBundle{1},dt::Vector,s::Tuple=size(f))
             d = similar(points(f))
             @threads for i ∈ OneTo(s[1])
                 d[i] = $cdg(f,s,i)/dt[i]
             end
             return d
         end
-        function $cdg(f::Grid{1},s::Tuple=size(f))
+        function $cdg(f::GridFrameBundle{1},s::Tuple=size(f))
             d = similar(points(f))
             @threads for i ∈ OneTo(s[1])
                 d[i] = $cdg(f,s,i)
             end
             return d
         end
-        function $cdg(f::Grid{2},dt::AbstractMatrix,s::Tuple=size(f))
+        function $cdg(f::GridFrameBundle{2},dt::AbstractMatrix,s::Tuple=size(f))
             d = Array{Chain{Submanifold(2),1,pointtype(f),2},2}(undef,s...)
             @threads for i ∈ OneTo(s[1]); for j ∈ OneTo(s[2])
                 d[i,j] = Chain($cdg(f,s,i,j).v./dt[i,j].v)
             end end
             return d
         end
-        function $cdg(f::Grid{2},s::Tuple=size(f))
+        function $cdg(f::GridFrameBundle{2},s::Tuple=size(f))
             d = Array{Chain{Submanifold(2),1,pointtype(f),2},2}(undef,s...)
             @threads for i ∈ OneTo(s[1]); for j ∈ OneTo(s[2])
                 d[i,j] = $cdg(f,s,i,j)
             end end
             return d
         end
-        function $cdg(f::Grid{3},dt::AbstractArray{T,3} where T,s::Tuple=size(f))
+        function $cdg(f::GridFrameBundle{3},dt::AbstractArray{T,3} where T,s::Tuple=size(f))
             d = Array{Chain{Submanifold(3),1,pointtype(f),3},3}(undef,s...)
             @threads for i ∈ OneTo(s[1]); for j ∈ OneTo(s[2]); for k ∈ OneTo(s[3])
                 d[i,j,k] = Chain($cdg(f,s,i,j,k).v./dt[i,j,k].v)
             end end end
             return d
         end
-        function $cdg(f::Grid{3},s::Tuple=size(f))
+        function $cdg(f::GridFrameBundle{3},s::Tuple=size(f))
             d = Array{Chain{Submanifold(3),1,pointtype(f),3},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,i,j,k)
             end end end
             return d
         end
-        function $cdg(f::Grid{4},dt::AbstractArray{T,4} where T,s::Tuple=size(f))
+        function $cdg(f::GridFrameBundle{4},dt::AbstractArray{T,4} where T,s::Tuple=size(f))
             d = Array{Chain{Submanifold(4),1,pointtype(f),4},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] = Chain($cdg(f,s,i,j,k,l).v./dt[i,j,k,l].v)
             end end end end
             return d
         end
-        function $cdg(f::Grid{4},s::Tuple=size(f))
+        function $cdg(f::GridFrameBundle{4},s::Tuple=size(f))
             d = Array{Chain{Submanifold(4),1,pointtype(f),4},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,i,j,k,l)
             end end end end
             return d
         end
-        function $cdg(f::Grid{5},dt::AbstractArray{T,5} where T,s::Tuple=size(f))
+        function $cdg(f::GridFrameBundle{5},dt::AbstractArray{T,5} where T,s::Tuple=size(f))
             d = Array{Chain{Submanifold(5),1,pointtype(f),5},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] = Chain($cdg(f,s,i,j,k,l,o).v./dt[i,j,k,l,o].v)
             end end end end end
             return d
         end
-        function $cdg(f::Grid{5},s::Tuple=size(f))
+        function $cdg(f::GridFrameBundle{5},s::Tuple=size(f))
             d = Array{Chain{Submanifold(5),1,pointtype(f),5},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,i,j,k,l,o)
             end end end end end
             return d
         end
-        function $cdg(f::Grid{2},n::Val{1},dt::AbstractVector,s::Tuple=size(f))
+        function $cdg(f::GridFrameBundle{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,(@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))
+        function $cdg(f::GridFrameBundle{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,(@inbounds s[2]),n,i,j)/dt[j]
             end end
             return d
         end
-        function $cdg(f::Grid{2},n::Val{N},dt::AbstractMatrix,s::Tuple=size(f)) where N
+        function $cdg(f::GridFrameBundle{2},n::Val{N},dt::AbstractMatrix,s::Tuple=size(f)) where N
             d = Array{pointtype(f),2}(undef,s...)
             @threads for i ∈ OneTo(s[1]); for j ∈ OneTo(s[2])
                 d[i,j] = $cdg(f,s[N],n,i,j)/dt[i,j]
             end end
             return d
         end
-        function $cdg(f::Grid{2},n::Val{N},s::Tuple=size(f)) where N
+        function $cdg(f::GridFrameBundle{2},n::Val{N},s::Tuple=size(f)) where N
             d = Array{pointtype(f),2}(undef,s...)
             @threads for i ∈ OneTo(s[1]); for j ∈ OneTo(s[2])
                 d[i,j] = $cdg(f,s[N],n,i,j)
             end end
             return d
         end
-        function $cdg(f::Grid{3},n::Val{1},dt::AbstractVector,s::Tuple=size(f))
+        function $cdg(f::GridFrameBundle{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,(@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))
+        function $cdg(f::GridFrameBundle{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,(@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))
+        function $cdg(f::GridFrameBundle{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,(@inbounds s[3]),n,i,j,k)/dt[k]
             end end end
             return d
         end
-        function $cdg(f::Grid{3},n::Val{N},dt::AbstractArray,s::Tuple=size(f)) where N
+        function $cdg(f::GridFrameBundle{3},n::Val{N},dt::AbstractArray,s::Tuple=size(f)) where N
             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[N],n,i,j,k)/dt[i,j,k]
             end end end
             return d
         end
-        function $cdg(f::Grid{3},n::Val{N},s::Tuple=size(f)) where N
+        function $cdg(f::GridFrameBundle{3},n::Val{N},s::Tuple=size(f)) where N
             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[N],n,i,j,k)
             end end end
             return d
         end
-        function $cdg(f::Grid{4},n::Val{1},dt::AbstractVector,s::Tuple=size(f))
+        function $cdg(f::GridFrameBundle{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,(@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))
+        function $cdg(f::GridFrameBundle{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,(@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))
+        function $cdg(f::GridFrameBundle{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,(@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))
+        function $cdg(f::GridFrameBundle{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,(@inbounds s[4]),n,i,j,k,l)/dt[l]
             end end end end
             return d
         end
-        function $cdg(f::Grid{4},n::Val{N},dt::AbstractArray,s::Tuple=size(f)) where N
+        function $cdg(f::GridFrameBundle{4},n::Val{N},dt::AbstractArray,s::Tuple=size(f)) where N
             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[N],n,i,j,k,l)/dt[i,j,k,l]
             end end end end
             return d
         end
-        function $cdg(f::Grid{4},n::Val{N},s::Tuple=size(f)) where N
+        function $cdg(f::GridFrameBundle{4},n::Val{N},s::Tuple=size(f)) where N
             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[N],n,i,j,k,l)
             end end end end
             return d
         end
-        function $cdg(f::Grid{5},n::Val{1},dt::AbstractVector,s::Tuple=size(f))
+        function $cdg(f::GridFrameBundle{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,(@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))
+        function $cdg(f::GridFrameBundle{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,(@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))
+        function $cdg(f::GridFrameBundle{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,(@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))
+        function $cdg(f::GridFrameBundle{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,(@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))
+        function $cdg(f::GridFrameBundle{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,(@inbounds s[5]),n,i,j,k,l,o)/dt[o]
             end end end end end
             return d
         end
-        function $cdg(f::Grid{5},n::Val{N},dt::AbstractArray,s::Tuple=size(f)) where N
+        function $cdg(f::GridFrameBundle{5},n::Val{N},dt::AbstractArray,s::Tuple=size(f)) where N
             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[N],n,i,j,k,l,o)/dt[i,j,k,l,o]
             end end end end end
             return d
         end
-        function $cdg(f::Grid{5},n::Val{N},s::Tuple=size(f)) where N
+        function $cdg(f::GridFrameBundle{5},n::Val{N},s::Tuple=size(f)) where N
             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[N],n,i,j,k,l,o)
             end end end end end
             return d
         end
-        $cdg(f::Grid{1},s::Tuple,i::Int) = $cdg(f,s[1],Val(1),i)
-        @generated function $cdg(f::Grid{N},s::Tuple,i::Vararg{Int}) where N
+        $cdg(f::GridFrameBundle{1},s::Tuple,i::Int) = $cdg(f,s[1],Val(1),i)
+        @generated function $cdg(f::GridFrameBundle{N},s::Tuple,i::Vararg{Int}) where N
             :(Chain($([:($$cdg(f,s[$n],Val($n),i...)) for n ∈ list(1,N)]...)))
         end
         $cd(f::RealRegion) = ProductSpace($cd.(f.v))
-        $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))),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))
+        $cd(f::GridFrameBundle{N,Coordinate{P,G},<:PointArray{P,G,N,<:RealRegion,<:Global{N,<:InducedMetric}},<:ProductTopology}) where {N,P,G} = ProductSpace($cd.(base(base(f)).v))
+        $cd(f::GridFrameBundle{N,Coordinate{P,G},<:PointArray{P,G,N,<:RealRegion,<:Global{N,<:InducedMetric}},<:OpenTopology}) where {N,P,G} = ProductSpace($cd.(base(base(f)).v))
+        $cd(f::GridFrameBundle{N,Coordinate{P,G},<:PointArray{P,G,N,<:RealRegion,<:Global{N,<:InducedMetric}}}) where {N,P,G} = sum.(value.($cdg(f)))
+        #$cd(f::GridFrameBundle{N,Coordinate{P,G},<:PointArray{P,G,N,<:RealRegion},<:ProductTopology}) where {N,P,G} = applymetric.($cd(base(base(f))),metricextensor(f))
+        $cd(f::GridFrameBundle{N,Coordinate{P,G},<:PointArray{P,G,N,<:RealRegion},<:OpenTopology}) where {N,P,G} = applymetric.($cd(base(base(f))),metricextensor(f))
+        $cd(f::GridFrameBundle{N,Coordinate{P,G},<:PointArray{P,G,N,<:RealRegion}}) where {N,P,G} = applymetric.(sum.(value.($cdg(f))),metricextensor(f))
+        $cd(f::AbstractRange,q::OpenTopology,s::Tuple=size(f)) = $cd(f,s)
+        function $cd(f::AbstractRange,q::QuotientTopology,s::Tuple=size(f))
+            d = Vector{eltype(f)}(undef,s[1])
+            @threads for i ∈ OneTo(s[1])
+                d[i] = $cdg(i,q,step(f),s[1])
+            end
+            return d
+        end
         function $cd(f::AbstractRange,s::Tuple=size(f))
             d = Vector{eltype(f)}(undef,s[1])
             @threads for i ∈ OneTo(s[1])
@@ -482,8 +460,7 @@ applymetric(f::Chain{V,G},g::Endomorphism{W,<:Simplex} where W) where {V,G} = ap
     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]
+function centraldiff_slow_calc(f::GridFrameBundle{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])
         18f[1,n,i...]-9f[2,n,i...]+2f[3,n,i...]-11points(f)[i...]
     elseif i[N]==l
@@ -496,35 +473,35 @@ function centraldiff_slow_calc(f::Grid{M,T,PA,<:OpenTopology} where {M,T,PA},l::
         f[-2,n,i...]+8(f[1,n,i...]-f[-1,n,i...])-f[2,n,i...]
     end
 end
-function Cartan.centraldiff_slow_calc(f::Grid,l::Int,n::Val{N},i::Vararg{Int}) where N #l=size(f)[N]
+function Cartan.centraldiff_slow_calc(f::GridFrameBundle,l::Int,n::Val{N},i::Vararg{Int}) where N #l=size(f)[N]
     if isone(i[N])
-        if iszero(f.cod.r[2N-1])
+        if iszero(immersion(f).r[2N-1])
             18f[1,n,i...]-9f[2,n,i...]+2f[3,n,i...]-11points(f)[i...]
-        elseif f.cod.p[2N-1]≠2N-1
+        elseif immersion(f).p[2N-1]≠2N-1
             f[-2,n,i...]+7(f[0,n,i...]-points(f)[i...])+8(f[1,n,i...]-f[-1,n,i...])-f[2,n,i...]
         else
             (-f[-2,n,i...])+7(-f[0,n,i...]-points(f)[i...])+8(f[1,n,i...]+f[-1,n,i...])-f[2,n,i...]
         end
     elseif i[N]==l
-        if iszero(f.cod.r[2N])
+        if iszero(immersion(f).r[2N])
             11points(f)[i...]-18f[-1,n,i...]+9f[-2,n,i...]-2f[-3,n,i...]
-        elseif f.cod.p[2N]≠2N
+        elseif immersion(f).p[2N]≠2N
             f[-2,n,i...]+8(f[1,n,i...]-f[-1,n,i...])+7(points(f)[i...]-f[0,n,i...])-f[2,n,i...]
         else
             f[-2,n,i...]+8(-f[1,n,i...]-f[-1,n,i...])+7(points(f)[i...]-f[0,n,i...])+f[2,n,i...]
         end
     elseif i[N]==2
-        if iszero(f.cod.r[2N-1])
+        if iszero(immersion(f).r[2N-1])
             6f[1,n,i...]-f[2,n,i...]-3points(f)[i...]-2f[-1,n,i...]
-        elseif f.cod.p[2N-1]≠2N-1
+        elseif immersion(f).p[2N-1]≠2N-1
             f[-2,n,i...]-f[-1,n,i...]+8f[1,n,i...]-7getpoint(f,-1,n,i...)-f[2,n,i...]
         else
             (-f[-2,n,i...])+f[-1,n,i...]+8f[1,n,i...]-7getpoint(f,-1,n,i...)-f[2,n,i...]
         end
     elseif i[N]==l-1
-        if iszero(f.cod.r[2N])
+        if iszero(immersion(f).r[2N])
             3points(f)[i...]-6f[-1,n,i...]+f[-2,n,i...]+2f[1,n,i...]
-        elseif f.cod.p[2N]≠2N
+        elseif immersion(f).p[2N]≠2N
             f[-2,n,i...]+7getpoint(f,1,n,i...)-8f[-1,n,i...]+f[1,n,i...]-f[2,n,i...]
         else
             f[-2,n,i...]+7getpoint(f,1,n,i...)-8f[-1,n,i...]-f[1,n,i...]+f[2,n,i...]
@@ -580,8 +557,15 @@ function centraldiff_slow_calc(i::Int,dt::Real,l::Int)
         12dt # (8-2)*2dt
     end
 end
+function centraldiff_slow_calc(i::Int,q::QuotientTopology,dt::Real,l::Int)
+    if (i∈(1,2) && iszero(q.r[1])) || (i∈(l-1,l) && iszero(q.r[2]))
+        6dt # (8-2)*dt
+    else
+        12dt # (8-2)*2dt
+    end
+end
 
-function centraldiff_fast_calc(f::Grid{M,T,PA,<:OpenTopology} where {M,T,PA},l::Int,n::Val{N},i::Vararg{Int}) where N
+function centraldiff_fast_calc(f::GridFrameBundle{M,T,PA,<:OpenTopology} where {M,T,PA},l::Int,n::Val{N},i::Vararg{Int}) where N
     if isone(i[N]) # 4f[1,k,i...]-f[2,k,i...]-3f.v[i...]
         18f[1,n,i...]-9f[2,n,i...]+2f[3,n,i...]-11points(f)[i...]
     elseif i[N]==l # 3f.v[i...]-4f[-1,k,i...]+f[-2,k,i...]
@@ -590,19 +574,19 @@ function centraldiff_fast_calc(f::Grid{M,T,PA,<:OpenTopology} where {M,T,PA},l::
         f[1,n,i...]-f[-1,n,i...]
     end
 end
-function Cartan.centraldiff_fast_calc(f::Grid,l::Int,n::Val{N},i::Vararg{Int}) where N
+function Cartan.centraldiff_fast_calc(f::GridFrameBundle,l::Int,n::Val{N},i::Vararg{Int}) where N
     if isone(i[N])
-        if iszero(f.cod.r[2N-1])
+        if iszero(immersion(f).r[2N-1])
             18f[1,n,i...]-9f[2,n,i...]+2f[3,n,i...]-11points(f)[i...]
-        elseif f.cod.p[2N-1]≠2N-1
+        elseif immersion(f).p[2N-1]≠2N-1
             (f[1,n,i...]-points(f)[i...])+(f[0,n,i...]-f[-1,n,i...])
         else # mirror
             (f[1,n,i...]-points(f)[i...])-(f[0,n,i...]-f[-1,n,i...])
         end
     elseif i[N]==l
-        if iszero(f.cod.r[2N])
+        if iszero(immersion(f).r[2N])
             11points(f)[i...]-18f[-1,n,i...]+9f[-2,n,i...]-2f[-3,n,i...]
-        elseif f.cod.p[2N]≠2N
+        elseif immersion(f).p[2N]≠2N
             (f[1,n,i...]-f[0,n,i...])+(points(f)[i...]-f[-1,n,i...])
         else # mirror
             (f[0,n,i...]-f[1,n,i...])+(points(f)[i...]-f[-1,n,i...])
@@ -617,6 +601,14 @@ centraldiff_fast_calc(::Val{1},i::Int,d1::Real,d2::Real,l::Int) = isone(i) ? d1+
 centraldiff_fast_calc(::Val{2},i::Int,d1::Real,d2::Real,l::Int) = i==l ? d1+d2 : 2dt
 centraldiff_fast_calc(i::Int,dt::Real,l::Int) = i∈(1,l) ? 6dt : 2dt
 #centraldiff_fast_calc(i::Int,dt::Real,l::Int) = 2dt
+function centraldiff_fast_calc(i::Int,q::QuotientTopology,dt::Real,l::Int)
+    if (isone(i) && iszero(q.r[1])) || (i==l && iszero(q.r[2]))
+        6dt
+    else
+        2dt
+    end
+end
+
 
 qnumber(n,q) = (q^n-1)/(q-1)
 qfactorial(n,q) = prod(cumsum([q^k for k ∈ 0:n-1]))
@@ -736,14 +728,14 @@ cumtrapz(f::IntervalMap,j::Val{1}) = cumtrapz(f)
 cumtrapz(f::ParametricMap,j::Int) = cumtrapz(f,Val(j))
 cumtrapz(f::ParametricMap,j::Val{J}) where J = TensorField(domain(f), cumtrapz2(codomain(f),j,diff(points(f).v[J])))
 cumtrapz(f::ParametricMap{B,F,N,<:AlignedSpace{N}} where {B,F,N},j::Val{J}) where J = TensorField(domain(f), cumtrapz1(codomain(f),j,step(points(f).v[J])))
-selectzeros(n,j) = :(zeros($([i≠j ? :(s[$i]) : 1 for i ∈ 1:n]...)))
-selectzeros2(n,j) = :(zeros($([i≠j ? i<j ? :(s[$i]) : :(s[$i]-1) : 1 for i ∈ 1:n]...)))
+selectzeros(n,j) = :(zeros(T,$([i≠j ? :(s[$i]) : 1 for i ∈ 1:n]...)))
+selectzeros2(n,j) = :(zeros(T,$([i≠j ? i<j ? :(s[$i]) : :(s[$i]-1) : 1 for i ∈ 1:n]...)))
 gencat(n,j=n,cat=n≠2 ? :cat : j≠2 ? :vcat : :hcat) = :($cat($(selectzeros2(n,j)),$(j≠1 ? gencat(n,j-1) : :i);$((cat≠:cat ? () : (Expr(:kw,:dims,j),))...)))
 gencumtrapz1(n,j,h=:h,f=:f) = :(($h/2)*cumsum($(select1(n,j,:(2:$(:end)),f)).+$(select1(n,j,:(1:$(:end)-1),f)),dims=$j))
-@generated function cumtrapz1(f::Array{T,N} where T,::Val{J},h::Real,s::Tuple=size(f)) where {N,J}
+@generated function cumtrapz1(f::Array{T,N},::Val{J},h::Real,s::Tuple=size(f)) where {T,N,J}
     :(cat($(selectzeros(N,J)),$(gencumtrapz1(N,J)),dims=$J))
 end
-@generated function cumtrapz1(f::Array{T,N} where T,h::D...) where {N,D<:Real}
+@generated function cumtrapz1(f::Array{T,N},h::D...) where {T,N,D<:Real}
     Expr(:block,:(s=size(f)),
          [:(i = $(gencumtrapz1(N,j,:(h[$j]),j≠1 ? :i : :f,))) for j ∈ 1:N]...,
         gencat(N))
@@ -763,14 +755,14 @@ for N ∈ 2:5
     @eval function cumtrapz(f::ParametricMap{B,F,$N,<:AlignedSpace{$N}} where {B,F})
         TensorField(domain(f), cumtrapz1(codomain(f),$([:(step(points(f).v[$j])) for j ∈ 1:N]...)))
     end
-    @eval function cumtrapz(m::ParametricMap{B,F,$N,T} where {B,F,T},D::AbstractVector=diff.(points(m).v))
+    @eval function cumtrapz(m::ParametricMap{B,F,$N,T} where {B,F},D::AbstractVector=diff.(points(m).v)) where T
         f = codomain(m)
         s,d = size(f),D./2
         $(Expr(:block,vcat([gencumtrapz2(N,j,:(d[$j])).args for j ∈ 1:N]...)...))
         TensorField(domain(m), $(gencat(N)))
     end
     for J ∈ 1:N
-        @eval function cumtrapz2(c::Array{T,$N} where T,::Val{$J})
+        @eval function cumtrapz2(c::Array{T,$N},::Val{$J}) where T
             s,d = size(f),D/2
             $(gencumtrapz2(N,J,:d,:f))
             cat($(selectzeros(N,J)),i,dims=$J)
@@ -783,31 +775,19 @@ end
 
 # differential geometry
 
-export 𝓛, Lie, LieBracket, LieDerivative, bracket, Connection, CovariantDerivative
+import Grassmann: 𝓛, Lie, LieBracket, LieDerivative, bracket
+export 𝓛, Lie, LieBracket, LieDerivative, bracket, Connection, CovariantDerivative, action
 
-struct LieBracket end
-struct LieDerivative{X<:VectorField}
-    v::X
-end
-const 𝓛 = LieBracket()
-const Lie = LieBracket()
-
-Base.getindex(::LieBracket,X::VectorField) = LieDerivative(X)
-Base.getindex(::LieBracket,X::VectorField,Y::VectorField...) = bracket(X,Y...)
-(::LieBracket)(X::VectorField,Y::VectorField...) = bracket(X,Y...)
-(::LieBracket)(X::VectorField) = LieDerivative(X)
-(X::LieDerivative)(Y::VectorField...) = bracket(X.v,Y...)
-(X::LieDerivative)(f::ScalarField) = X.v(f)
-(X::VectorField)(f::ScalarField) = X⋅gradient(f)
-(X::GradedVector)(f::ScalarField) = X⋅gradient(f)
-LieBracket(X::VectorField,Y::VectorField...) = bracket(X,Y...)
-LieBracket(X::VectorField) = LieDerivative(X)
-bracket(X,Y) = X(Y) - Y(X)
-bracket(X,Y,Z) = X(bracket(Y,Z)) + Y(bracket(Z,X)) + Z(bracket(X,Y))
-bracket(W,X,Y,Z) = W(bracket(X,Y,Z)) + X(bracket(W,Z,Y)) + Y(bracket(W,X,Z)) + Z(bracket(W,Y,X))
-bracket(V,W,X,Y,Z) = V(bracket(W,X,Y,Z)) + W(bracket(V,X,Z,Y)) + X(bracket(V,W,Y,Z)) + Y(bracket(V,W,Z,X)) + Z(bracket(V,W,X,Y))
+(X::LieDerivative)(f::ScalarField) = action(X.v,f)
+(X::VectorField{B,<:Chain{V,1,T,N} where T,N} where B)(Y::VectorField{B,<:Chain{V,1,T,N} where T,N} where B) where {V,N} = action(X,Y)
+(X::VectorField)(f::ScalarField) = action(X,f)
+(X::GradedVector)(f::ScalarField) = action(X,f)
+(X::ScalarField)(f::ScalarField) = action(X,f)
 𝓛dot(x::Chain,y::Simplex{V}) where V = Chain{V}(Real.(x.⋅value(y)))
-function (X::VectorField)(Y::VectorField)
+action(X::VectorField,f::ScalarField) = X⋅gradient(f)
+action(X::GradedVector,f::ScalarField) = X⋅gradient(f)
+action(X::ScalarField,f::ScalarField) = X⋅gradient(f)
+function action(X::VectorField,Y::VectorField)
     TensorField(base(X),𝓛dot.(fiber(X),fiber(gradient(Y))))
 end
 
@@ -830,10 +810,15 @@ CovariantDerivative(∇::Connection,X) = ∇(X)
 
 export arclength, arctime, totalarclength, trapz, cumtrapz, linecumtrapz, psum, pcumsum
 export centraldiff, tangent, tangent_fast, unittangent, speed, normal, unitnormal
-export curvenormal, unitcurvenormal, ribbon, tangentsurface, planecurve
-export normalnorm, area, surfacearea, weingarten, gauss
-export normalnorm_slow, area_slow, surfacearea_slow, weingarten_slow, gauss_slow
-export tangent_slow, normal_slow, unittangent_slow, unitnormal_slow
+export curvenormal, unitcurvenormal, ribbon, tangentsurface, planecurve, link, linkmap
+export normalnorm, area, surfacearea, weingarten, gausssign, jacobian
+export normalnorm_slow, area_slow, surfacearea_slow, weingarten_slow, gausssign_slow
+export gaussintrinsic, gaussextrinsic, gaussintrinsicnorm, gaussextrinsicnorm
+export gaussintrinsic_slow, gausseintrinsicnorm_slow, curvatures, meancurvature
+export gaussextrinsic_slow, gaussextrinsicnorm_slow, principals, principalaxes
+export tangent_slow, normal_slow, unittangent_slow, unitnormal_slow, jacobian_slow
+export sector, sectordet, sectorintegral, sectorintegrate, linkintegral, linknumber
+export sector_slow, sectordet_slow, sectorintegral_slow, sectorintegrate_slow
 
 # use graph for IntervalMap? or RealFunction!
 tangent(f::IntervalMap) = gradient(f)
@@ -846,6 +831,8 @@ unittangent(f::VectorField,n=tangent(f)) = TensorField(domain(f), codomain(n)./n
 unittangent(f::IntervalMap) = unitgradient(f)
 unitnormal(f) = ⋆unittangent(f)
 normalnorm(f) = Real(abs(normal(f)))
+jacobian(f::VectorField) = TensorOperator(gradient(f))
+weingarten(f::VectorField) = jacobian(unitnormal(f))
 
 tangent_slow(f::IntervalMap) = gradient_slow(f)
 tangent_slow(f::ScalarField) = tangent_slow(graph(f))
@@ -857,26 +844,70 @@ unittangent_slow(f::VectorField,n=tangent_slow(f)) = TensorField(domain(f), codo
 unittangent_slow(f::IntervalMap) = unitgradient_slow(f)
 unitnormal_slow(f) = ⋆unittangent_slow(f)
 normalnorm_slow(f) = Real(abs(normal_slow(f)))
+jacobian_slow(f::VectorField) = TensorOperator(gradient_slow(f))
+weingarten_slow(f::VectorField) = jacobian_slow(unitnormal_slow(f))
 
 ribbon(f::AbstractCurve,g::Vector{<:AbstractCurve}) = TensorField(points(f)⊕LinRange(0,1,length(g)+1),hcat(fiber(f),fiber.(g)...))
 ribbon(f::AbstractCurve,g::AbstractCurve,n::Int=100) = tangentsurface(f,g-f,n)
 tangentsurface(f::AbstractCurve,g::AbstractCurve,n::Int=100) = ribbon(f,Ref(f).+(LinRange(inv(n),1,n).*Ref(g)))
 tangentsurface(f::AbstractCurve,v::Real=1,n::Int=100) = tangentsurface(f,v*tangent(f),n)
 
+sector(f::Chain{V,G,<:Real},J::Chain{V,G,<:Real}) where {V,G} = Chain{V}(f,J)
+sector(f::Real,J::Chain{W,1,<:Real} where W) = Chain(f,value(J)...)
+sector(f::Chain{V,G},J::Chain{W,1,<:Chain{V,G}}) where {W,V,G} = Chain{V}(f,value(J)...)
+sector(f::TensorField,J::TensorField) = TensorField(base(f), sector.(fiber(f),fiber(J)))
+sector(f::RealFunction) = f
+sector(f::TensorField) = TensorOperator(sector(f,gradient(f)))
+sectordet(f::RealFunction) = f
+sectordet(f::PlaneCurve) = f∧gradient(f)
+sectordet(f::VectorField) = f∧det(gradient(f))
+sectorintegral(f::RealFunction) = integral(f)
+sectorintegral(f::TensorField) = integral(sectordet(f))/mdims(fibertype(f))
+sectorintegrate(f::RealFunction) = integrate(f)
+sectorintegrate(f::TensorField) = integrate(sectordet(f))/mdims(fibertype(f))
+sector_slow(f::RealFunction) = f
+sector_slow(f::TensorField) = TensorOperator(sector(f,gradient_slow(f)))
+sectordet_slow(f::RealFunction) = f
+sectordet_slow(f::PlaneCurve) = f∧gradient_slow(f)
+sectordet_slow(f::VectorField) = f∧det(gradient_slow(f))
+sectorintegral_slow(f::RealFunction) = integral(f)
+sectorintegral_slow(f::TensorField) = integral(sectordet_slow(f))/mdims(fibertype(f))
+sectorintegrate_slow(f::RealFunction) = integrate(f)
+sectorintegrate_slow(f::TensorField) = integrate(sectordet_slow(f))/mdims(fibertype(f))
+
 area(f::VectorField) = integral(normalnorm(f))
 surfacearea(f::VectorField) = integrate(normalnorm(f))
-weingarten(f::VectorField) = gradient(unitnormal(f))
-gauss(f::VectorField) = Real(abs(det(weingarten(f))))
+principals(f::VectorField) = eigvals(shape(f))
+principals(f::VectorField,i) = eigvals(shape(f),i)
+principalaxes(f::VectorField) = eigvecs(shape(f))
+principalaxes(f::VectorField,i) = eigvecs(shape(f),i)
+curvatures(f::VectorField) = eigpolys(shape(f))
+curvatures(f::VectorField,i) = eigpolys(shape(f),i)
+meancurvature(f::VectorField) = eigpolys(shape(f),Val(1))
+gausssign(f::VectorField) = sign(sectordet(unitnormal(f)))
+gaussextrinsicnorm(f::VectorField) = norm(gaussextrinsic(f))
+gaussintrinsicnorm(f::VectorField,N=normal(f)) = norm(gaussintrinsic(f,N))
+gaussextrinsic(f::VectorField) = sectordet(unitnormal(f))
+function gaussintrinsic(f::VectorField,N=normal(f))
+    n,T = Real(abs(N)),pointtype(f)
+    Chain{Manifold(T),mdims(T)}(Real(sectordet(N/n))/n)
+end
 
 area_slow(f::VectorField) = integral(normalnorm_slow(f))
 surfacearea_slow(f::VectorField) = integrate(normalnorm_slow(f))
-weingarten_slow(f::VectorField) = gradient(unitnormal_slow(f))
-gauss_slow(f::VectorField) = Real(abs(det(weingarten_slow(f))))
+gausssign_slow(f::VectorField) = sign(sectordet_slow(unitnormal_slow(f)))
+gaussextrinsicnorm_slow(f::VectorField) = norm(gaussextrinsic_slow(f))
+gaussintrinsicnorm_slow(f::VectorField,N=normal(f)) = norm(gaussintrinsic_slow(f,N))
+gaussextrinsic_slow(f::VectorField) = sectordet_slow(unitnormal_slow(f))
+function gaussintrinsic_slow(f::VectorField,N=normal_slow(f))
+    n,T = Real(abs(N)),pointtype(f)
+    Chain{Manifold(T),mdims(T)}(Real(sectordet_slow(N/n))/n)
+end
 
-area(f::PlaneCurve) = integral(f∧gradient(f))/2
-#volume(f::VectorField) = integral(f∧det(gradient(f)))/3
-sectorarea(f::PlaneCurve) = integrate(f∧gradient(f))/2
-sectorvolume(f::VectorField) = integrate(f∧det(gradient(f)))/3
+area(f::PlaneCurve) = sectorintegral(f)
+#volume(f::VectorField) = sectorintegral(f)
+sectorarea(f::PlaneCurve) = sectorintegrate(f)
+sectorvolume(f::VectorField) = sectorintegrate(f)
 
 function speed(f::IntervalMap,d=centraldiffpoints(f),t=centraldifffiber(f,d))
     TensorField(domain(f), abs.(t))
@@ -884,7 +915,7 @@ end
 function curvenormal(f::IntervalMap,d=centraldiffpoints(f),t=centraldifffiber(f,d))
     TensorField(domain(f), centraldiff(t,d))
 end
-function unitcurvenormal(f::IntervalMap,d=centraldiffpoints(f),t=centraldifffiiber(f,d),n=centraldiff(t,d))
+function unitcurvenormal(f::IntervalMap,d=centraldiffpoints(f),t=centraldifffiber(f,d),n=centraldiff(t,d))
     TensorField(domain(f), (n./abs.(n)))
 end
 
@@ -950,6 +981,17 @@ function wronskian(f::SpaceCurve,d=centraldiffpoints(f),t=centraldifffiber(f,d),
     TensorField(domain(f), f.cod.∧t.∧n)
 end
 
+function linkmap(f::SpaceCurve,g::SpaceCurve)
+    TensorField(base(f)×base(g),[fiber(g)[j]-fiber(f)[i] for i ∈ OneTo(length(f)), j ∈ OneTo(length(g))])
+end
+link(tf::Chain,tg::Chain,f::Chain,g::Chain) = (gf = g-f; ∧(gf,tf,tg)/norm(gf)^3)
+function link(f::SpaceCurve,g::SpaceCurve)
+    tf,tg = fiber(tangent(f)),fiber(tangent(g))
+    TensorField(base(f)×base(g),[link(tf[i],tg[j],fiber(f)[i],fiber(g)[j]) for i ∈ OneTo(length(f)), j ∈ OneTo(length(g))])
+end
+linkintegral(f::SpaceCurve,g::SpaceCurve) = integral(link(f,g))/4π
+linknumber(f::SpaceCurve,g::SpaceCurve) = integrate(link(f,g))/4π
+
 function frame(f::VectorField{<:Coordinate{<:Chain},<:Chain,2,<:RealSpace{2}})#,d=centraldiff(points(f)),t=centraldiff(codomain(f),d),n=centraldiff(t,d))
     Ψ = gradient(f)
     Ψu,Ψv = getindex.(Ψ,1),getindex.(Ψ,2)
@@ -999,64 +1041,115 @@ function planecurve(κ::RealFunction,φ::Real=0.0)
     integral(Chain.(cos(int),sin(int)))
 end
 
-export surfacemetric, surfacemetricdiag, surfaceframe
+export surfacemetric, surfacemetricdiag, surfaceframe, shape
+export firstform, firstformdiag, secondform, firstsecondform, thirdform
 export applymetric, firstkind, secondkind, geodesic
 
-surfacemetric(dom::ScalarField,f::Function) = surfacemetric(TensorField(dom,f))
-function surfacemetric(t::ScalarField)
-    g = gradient(t)
-    V = Submanifold(MetricTensor([1 1; 1 1]))
+function EFG(V,dfdx,dfdy)
+    E,F,G = fiber(Real(dfdx⋅dfdx)),fiber(Real(dfdx⋅dfdy)),fiber(Real(dfdy⋅dfdy))
+    TensorOperator.(Chain{V}.(Chain{V}.(E,F),Chain{V}.(F,G)))
+end
+function EG(V,dfdx,dfdy)
+    DiagonalOperator.(Chain{V}.(fiber(Real(dfdx⋅dfdx)),fiber(Real(dfdy⋅dfdy))))
+end
+function LMN(V,n,ddfdx2,ddfdxdy,ddfdy2)
+    L,M,N = fiber(Real(ddfdx2⋅n)),fiber(Real(ddfdxdy⋅n)),fiber(Real(ddfdy2⋅n))
+    TensorOperator.(Chain{V}.(Chain{V}.(L,M),Chain{V}.(M,N)))
+end
+
+firstform(dom::ScalarField,f::Function) = firstform(TensorField(dom,f))
+function firstform(t::ScalarField,g=gradient(t),V=Submanifold(2))
     dfdx,dfdy = getindex.(g,1),getindex.(g,2)
-    E,F,G = dfdx*dfdx,dfdx*dfdy,dfdy*dfdy
-    g1,g2 = (1+E)*Λ(V).v1+F*Λ(V).v2,(1+G)*Λ(V).v2+F*Λ(V).v1
-    EFG = Outermorphism.(Chain{V}.(fiber(g1),fiber(g2)))
-    GridFrameBundle(PointArray(points(t),EFG), immersion(t))
+    E,F,G = fiber(dfdx*dfdx),fiber(dfdx*dfdy),fiber(dfdy*dfdy)
+    TensorField(base(t),TensorOperator.(Chain{V}.(Chain{V}.(E.+1,F),Chain{V}.(F,G.+1))))
 end
 
-surfacemetricdiag(dom::ScalarField,f::Function) = surfacemetricdiag(TensorField(dom,f))
-function surfacemetricdiag(t::ScalarField)
-    g = gradient(t)
-    V = Submanifold(DiagonalForm(Values(1,1)))
+firstformdiag(dom::ScalarField,f::Function) = firstformdiag(TensorField(dom,f))
+function firstformdiag(t::ScalarField,g=gradient(t),V=Submanifold(2))
     dfdx,dfdy = getindex.(g,1),getindex.(g,2)
-    g1,g2 = fiber(1+dfdx*dfdx),fiber(1+dfdy*dfdy)
-    EG = outermorphism.(DiagonalOperator.(Chain{V}.(g1,g2)))
-    GridFrameBundle(PointArray(points(t),EG),immersion(t))
+    E1,G1 = fiber(1+dfdx*dfdx),fiber(1+dfdy*dfdy)
+    TensorField(base(t),DiagonalOperator.(Chain{V}.(E1,G1)))
+end
+
+firstform(dom,f::Function) = firstform(TensorField(dom,f))
+function firstform(t,g=gradient(t),V=Submanifold(2))
+    TensorField(base(t),EFG(V,getindex.(g,1),getindex.(g,2)))
+end
+
+firstformdiag(dom,f::Function) = firstformdiag(TensorField(dom,f))
+function firstformdiag(t,g=gradien(t),V=Submanifold(2))
+    TensorField(base(t),EG(V,getindex.(g,1),getindex.(g,2)))
+end
+
+secondform(dom,f::Function) = secondform(TensorField(dom,f))
+function secondform(t,g=gradient(t),V=Submanifold(2))
+    n = ⋆unittangent(t,det(g))
+    dfdx,dfdy = getindex.(g,1),getindex.(g,2)
+    ddfdx,ddfdy2 = gradient(dfdx),gradient(dfdy,Val(2))
+    ddfdx2,ddfdxdy = getindex.(ddfdx,1),getindex.(ddfdx,2)
+    TensorField(base(t),LMN(V,n,ddfdx2,ddfdxdy,ddfdy2))
+end
+
+firstsecondform(dom,f::Function) = firstsecondform(TensorField(dom,f))
+function firstsecondform(t,g=gradient(t),V=Submanifold(2))
+    n = ⋆unittangent(t,det(g))
+    dfdx,dfdy = getindex.(g,1),getindex.(g,2)
+    ddfdx,ddfdy2 = gradient(dfdx),gradient(dfdy,Val(2))
+    ddfdx2,ddfdxdy = getindex.(ddfdx,1),getindex.(ddfdx,2)
+    return (TensorField(base(t),EFG(V,dfdx,dfdy)),
+        TensorField(base(t),LMN(V,n,ddfdx2,ddfdxdy,ddfdy2)))
+end
+
+thirdform(dom,f::Function) = thirdform(TensorField(dom,f))
+thirdform(t,V=Submanifold(2)) = firstform(t,gradient(unitnormal(t)),V)
+
+shape(dom,f::Function) = shape(TensorField(dom,f))
+function shape(t,g=gradient(t),V=Submanifold(2))
+    EFG,LMN = firstsecondform(t,g,V)
+    return inv(EFG)⋅LMN
 end
 
 surfacemetric(dom,f::Function) = surfacemetric(TensorField(dom,f))
-function surfacemetric(t)
-    g = gradient(t)
+function surfacemetric(t::DiagonalField)
+    GridFrameBundle(PointArray(points(t),fiber(outermorphism(t))),immersion(t))
+end
+function surfacemetric(t::EndomorphismField)
+    GridFrameBundle(PointArray(points(t),fiber(Outermorphism(t))),immersion(t))
+end
+function surfacemetric(t,g=gradient(t))
     V = Submanifold(MetricTensor([1 1; 1 1]))
-    dfdx,dfdy = getindex.(g,1),getindex.(g,2)
-    E,F,G = Real(dfdx⋅dfdx),Real(dfdx⋅dfdy),Real(dfdy⋅dfdy)
-    g1,g2 = E*Λ(V).v1+F*Λ(V).v2,G*Λ(V).v2+F*Λ(V).v1
-    EFG = Outermorphism.(Chain{V}.(fiber(g1),fiber(g2)))
+    EFG = Outermorphism.(fiber(firstform(t,g,V)))
     GridFrameBundle(PointArray(points(t),EFG),immersion(t))
 end
 
 surfacemetricdiag(dom,f::Function) = surfacemetricdiag(TensorField(dom,f))
-function surfacemetricdiag(t)
-    g = gradient(t)
+surfacemetricdiag(t::DiagonalField) = surfacemetric(t)
+function surfacemetricdiag(t,g=gradient(t))
     V = Submanifold(DiagonalForm(Values(1,1)))
-    dfdx,dfdy = getindex.(g,1),getindex.(g,2)
-    g1,g2 = fiber(Real(dfdx⋅dfdx)),fiber(Real(dfdy⋅dfdy))
-    EG = outermorphism.(DiagonalOperator.(Chain{V}.(g1,g2)))
+    EG = outermorphism.(fiber(firstformdiag(t,g,V)))
     GridFrameBundle(PointArray(points(t),EG),immersion(t))
 end
 
-function surfaceframe(t)
-    g = getindex.(metricextensor(t),1)
-    V = Submanifold(MetricTensor([1 1; 1 1]))
-    E,F,G = getindex.(g,1,1),getindex.(g,1,2),getindex.(g,2,2)
+surfaceframe(t::DiagonalField) = surfaceframediag(t)
+surfaceframe(t::AbstractFrameBundle) = surfaceframe(metrictensorfield(t))
+function surfaceframe(t::EndomorphismField)
+    surfaceframe(base(t),getindex.(t,1,1),getindex.(t,1,2),getindex.(t,2,2))
+end
+function surfaceframe(t,g=gradient(t))
+    dx,dy = getindex.(g,1),getindex.(g,2)
+    surfaceframe(base(t),fiber(Real(dx⋅dx)),fiber(Real(dx⋅dy)),fiber(Real(dy⋅dy)))
+end
+function surfaceframe(b,E,F,G)
     F2 = F.*F; mag,sig = sqrt.((E.*E).+F2), sign.(F2.-(E.*G))
-    TensorOperator.(Chain.(Chain.(E,F)./mag,Chain.(F,.-E)./(sig.*mag)))
+    TensorField(b,TensorOperator.(Chain.(Chain.(E,F)./mag,Chain.(F,.-E)./(sig.*mag))))
 end
-function surfaceframediag(t)
-    g = getindex.(metricextensor(t),1)
-    V = Submanifold(DiagonalForm(Values(1,1)))
-    E,G = getindex.(g,1,1),getindex.(g,2,2)
+
+surfaceframediag(t) = surfaceframediag(firstformdiag(t))
+surfaceframediag(t::AbstractFrameBundle) = surfaceframediag(metrictensorfield(t))
+function surfaceframediag(t::DiagonalField)
+    E,G = getindex.(value.(fiber(g)),1),getindex.(value.(fiber(g)),2)
     mag,sig = sqrt.(E.*E), sign.(.-(E.*G))
-    DiagonalOperator.(Chain.(E./mag,.-E./(sig.*mag)))
+    TensorField(base(t),DiagonalOperator.(Chain.(E./mag,.-E./(sig.*mag))))
 end
 
 _firstkind(dg,k,i,j) = dg[k,j][i] + dg[i,k][j] - dg[i,j][k]
@@ -1092,8 +1185,9 @@ end
     Expr(:call,:TensorOperator,Expr(:call,:Chain,[:(secondkind(ig,dg,$j)) for j ∈ list(1,mdims(fibertype(dg)))]...))
 end
 
-geodesic(Γ::TensorField) = x -> geodesic(x,Γ)
+geodesic(Γ) = x -> geodesic(x,Γ)
 geodesic(g::AbstractFrameBundle) = geodesic(secondkind(g))
+geodesic(x,Γ::Function) = (x2 = x[2]; Chain(x2,-geodesic(x2,Γ(x[1]))))
 geodesic(x,Γ::TensorField) = (x2 = x[2]; Chain(x2,-geodesic(x2,Γ(x[1]))))
 @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)]...)...)
diff --git a/src/topology.jl b/src/topology.jl
index 93e394c..afbe5a6 100644
--- a/src/topology.jl
+++ b/src/topology.jl
@@ -76,6 +76,9 @@ cross(a::ProductSpace,b::ProductSpace) = a⊕b
 @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))
 
+widths(t::AbstractVector) = t[end]-t[1]
+widths(t::ProductSpace) = widths.(t.v)
+
 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]...))]))
@@ -245,7 +248,7 @@ HopfTopology(n::Values{2,Int}) = QuotientTopology(Values(2,1,4,3),Values(Product
 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]:-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)
+PolarTopology(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{1,Int}) = TorusTopology(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)
 SphereTopology(n::Values{3,Int}) = QuotientTopology(Values(1,2,3,4,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)
@@ -281,7 +284,7 @@ HopfParameter(n::Values{2,Int}) = TensorField(GridFrameBundle(PointArray(LinRang
 HopfParameter(n::Values{3,Int}) = TensorField(GridFrameBundle(PointArray((LinRange(7π/16/n[1],7π/16,n[1])⊕LinRange(0,2π,n[2]))⊕LinRange(0,4π,n[3])),HopfTopology(n)))
 KleinParameter(n::Values{2,Int}) = TensorField(GridFrameBundle(PointArray(LinRange(0,2π,n[1])⊕LinRange(0,2π,n[2])),KleinTopology(n)))
 ConeParameter(n::Values{2,Int}) = TensorField(GridFrameBundle(PointArray(LinRange(0,1,n[1])⊕LinRange(0,2π,n[2])),ConeTopology(n)))
-DrumParameter(n::Values{2,Int}) = TensorField(GridFrameBundle(PointArray(LinRange(0,1,n[1])⊕LinRange(0,2π,n[2])),DrumTopology(n)))
+PolarParameter(n::Values{2,Int}) = TensorField(GridFrameBundle(PointArray(LinRange(0,1,n[1])⊕LinRange(0,2π,n[2])),PolarTopology(n)))
 SphereParameter(n::Values{1,Int}) = TorusParameter(n)
 SphereParameter(n::Values{2,Int}) = TensorField(GridFrameBundle(PointArray(LinRange(0,π,n[1])⊕LinRange(0,2π,n[2])),SphereTopology(n)))
 SphereParameter(n::Values{3,Int}) = TensorField(GridFrameBundle(PointArray((LinRange(0,π,n[1])⊕LinRange(0,π,n[2]))⊕LinRange(0,2π,n[3])),SphereTopology(n)))
@@ -310,7 +313,7 @@ for mod ∈ (:Topology,:Parameter)
             end
         end
     end
-    for (fun,(n,m)) ∈ ((:Ribbon,(60,20)),(:Wing,(60,20)),(:Mobius,(60,20)),(:Klein,(60,60)),(:Cone,(30,:(2n+1))),(:Drum,(30,:(2n))),(:Sphere,(30,:(2n+1))),(:Geographic,(:(2n+1),30)))
+    for (fun,(n,m)) ∈ ((:Ribbon,(60,20)),(:Wing,(60,20)),(:Mobius,(60,20)),(:Klein,(60,60)),(:Cone,(30,:(2n+1))),(:Polar,(30,:(2n))),(:Sphere,(30,:(2n+1))),(:Geographic,(:(2n+1),30)))
         for typ ∈ (Symbol(fun,mod),)
             @eval begin
                 export $typ
@@ -327,9 +330,20 @@ iscompact(t::QuotientTopology) = false
 iscompact(t::CompactTopology) = true
 _to_axis(f::Int) = (iseven(f) ? f : f+1)÷2
 
+zeroprodtop(r,n) = iszero(r) ? () : (ProductTopology(n),)
+LinearAlgebra.cross(m::OpenTopology,n::OpenTopology) = OpenTopology(m.s...,n.s...)
+LinearAlgebra.cross(m::OpenTopology{1},n::OpenTopology{1}) = OpenTopology(m.s...,n.s...)
+function LinearAlgebra.cross(m::QuotientTopology{1},n::QuotientTopology{1})
+    M,N = m.s[1],n.s[1]
+    QuotientTopology(Values(m.p...,(n.p.+2)...),
+        Values((zeroprodtop(m.r[1],N)...,zeroprodtop(m.r[2],N)...,zeroprodtop(n.r[1],M)...,zeroprodtop(n.r[2],M)...)),
+        Values(m.r...,iszero(n.r[1]) ? 0 : n.r[1]+length(m.p),iszero(n.r[2]) ? 0 : n.r[2]+length(m.p)),Values(M,N))
+end
+LinearAlgebra.cross(m::OpenTopology,n::Int) = OpenTopology(m.s...,n)
+LinearAlgebra.cross(m::OpenTopology{1},n::Int) = OpenTopology(m.s...,n)
 function LinearAlgebra.cross(m::QuotientTopology{1},n::Int)
     QuotientTopology(m.p,
-        Values(ProductTopology(n),ProductTopology(n)),
+        Values((zeroprodtop(m.r[1],n)...,zeroprodtop(m.r[2],n)...)),
         Values(m.r...,0,0),Values(m.s...,n))
 end
 function LinearAlgebra.cross(m::QuotientTopology,n::Int)
@@ -594,6 +608,21 @@ subtopology(m::QuotientTopology{2},i::NTuple{N,Int},::Colon,j::NTuple{M,Int},::C
 subtopology(m::QuotientTopology{3},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} = m
 subtopology(m::QuotientTopology{4},i::NTuple{N,Int},::Colon,j::NTuple{M,Int},::Colon,k::NTuple{L,Int},::Colon,l::NTuple{O,Int},::Colon,o::NTuple{Y,Int} where Y) where {N,M,L,O} = m
 subtopology(m::QuotientTopology{5},i::NTuple{N,Int},::Colon,j::NTuple{M,Int},::Colon,k::NTuple{L,Int},::Colon,l::NTuple{O,Int},::Colon,o::NTuple{Y,Int},::Colon,z::NTuple{Z,Int} where Z) where {N,M,L,O,Y} = m
+function subtopology(m::QuotientTopology{M},::Val{N}) where {M,N}
+    r1,r2 = m.r[2N-1],m.r[2N]
+    r1z,r2z,n = iszero(r1),iszero(r2),size(m)[N]
+    if r1z
+        if r2z
+            OpenTopology(n)
+        else
+            QuotientTopology(Values((isodd(m.p[r2]) ? 1 : 2,)),Array{Values{0,Int},0}.(Values((undef,))),Values(0,1),Values((n,)))
+        end
+    elseif r2z && !r1z
+        QuotientTopology(Values((isodd(m.p[r1]) ? 1 : 2,)),Array{Values{0,Int},0}.(Values((undef,))),Values(1,0),Values((n,)))
+    else
+        QuotientTopology(Values(isodd(m.p[r1]) ? 1 : 2,isodd(m.p[r2]) ? 1 : 2),Array{Values{0,Int},0}.(Values((undef,undef))),Values(1,2),Values((n,)))
+    end
+end
 function subtopology(m::QuotientTopology,i::NTuple{N,Int},::Colon,j::NTuple{M,Int} where M) where N
     N1 = N+1
     r,vals = m.r[Values(2N1-1,2N1)],Values(i...,j...)
@@ -917,6 +946,18 @@ 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)
+Base.sign(s::LocalTensor) = LocalTensor(base(s),sign(Real(fiber(s))))
+Base.inv(a::LocalTensor{B,<:Real} where B) = LocalTensor(base(a), inv(fiber(a)))
+Base.inv(a::LocalTensor{B,<:Complex} where B) = LocalTensor(base(a), inv(fiber(a)))
+Base.:/(a::LocalTensor,b::LocalTensor{B,<:Real} where B) = LocalTensor(base(a), fiber(a)/fiber(b))
+Base.:/(a::LocalTensor,b::LocalTensor{B,<:Complex} where B) = LocalTensor(base(a), fiber(a)/fiber(b))
+Grassmann.eigen(t::LocalTensor,i::Val) = LocalTensor(base(t), eigen(fiber(t),i))
+Grassmann.eigen(t::LocalTensor,i::Int) = LocalTensor(base(t), eigen(fiber(t),i))
+Grassmann.eigvals(t::LocalTensor,i::Val) = LocalTensor(base(t), eigvals(fiber(t),i))
+Grassmann.eigvals(t::LocalTensor,i::Int) = LocalTensor(base(t), eigvals(fiber(t),i))
+Grassmann.eigvecs(t::LocalTensor,i::Val) = LocalTensor(base(t), eigvecs(fiber(t),i))
+Grassmann.eigvecs(t::LocalTensor,i::Int) = LocalTensor(base(t), eigvecs(fiber(t),i))
+Grassmann.eigpolys(t::LocalTensor,G::Val) = LocalTensor(base(t), eigpolys(fiber(t),G))
 for type ∈ (:Coordinate,:LocalSection,:LocalTensor)
     for tensor ∈ (:Single,:Couple,:PseudoCouple,:Chain,:Spinor,:AntiSpinor,:Multivector,:DiagonalOperator,:TensorOperator,:Outermorphism)
         @eval (T::Type{<:$tensor})(s::$type) = $type(base(s), T(fiber(s)))
@@ -927,10 +968,10 @@ for type ∈ (:Coordinate,:LocalSection,:LocalTensor)
     for fun ∈ (:inv,:exp,:exp2,:exp10,:log,:log2,:log10,:sinh,:cosh,:abs,:sqrt,:cbrt,:cos,:sin,:tan,:cot,:sec,:csc,:asec,:acsc,:sech,:csch,:asech,:tanh,:coth,:asinh,:acosh,:atanh,:acoth,:asin,:acos,:atan,:acot,:sinc,:cosc,:cis,:abs2)
         @eval Base.$fun(s::$type) = $type(base(s), $fun(fiber(s),metricextensor(base(s))))
     end
-    for fun ∈ (:reverse,:involute,:clifford,:even,:odd,:scalar,:vector,:bivector,:volume,:value,:curl,:∂,:d,:complementleft,:realvalue,:imagvalue,:outermorphism,:Outermorphism,:DiagonalOperator,:TensorOperator)
+    for fun ∈ (:reverse,:involute,:clifford,:even,:odd,:scalar,:vector,:bivector,:volume,:value,:curl,:∂,:d,:complementleft,:realvalue,:imagvalue,:outermorphism,:Outermorphism,:DiagonalOperator,:TensorOperator,:eigen,:eigvecs,:eigvals,:eigvalsreal,:eigvalscomplex,:eigvecsreal,:eigvecscomplex,:eigpolys)
         @eval Grassmann.$fun(s::$type) = $type(base(s), $fun(fiber(s)))
     end
-    for fun ∈ (:⋆,:angle,:radius,:complementlefthodge,:pseudoabs,:pseudoabs2,:pseudoexp,:pseudolog,:pseudoinv,:pseudosqrt,:pseudocbrt,:pseudocos,:pseudosin,:pseudotan,:pseudocosh,:pseudosinh,:pseudotanh,:metric)
+    for fun ∈ (:⋆,:angle,:radius,:complementlefthodge,:pseudoabs,:pseudoabs2,:pseudoexp,:pseudolog,:pseudoinv,:pseudosqrt,:pseudocbrt,:pseudocos,:pseudosin,:pseudotan,:pseudocosh,:pseudosinh,:pseudotanh,:metric,:unit)
         @eval Grassmann.$fun(s::$type) = $type(base(s), $fun(fiber(s),metricextensor(base(s))))
     end
     for op ∈ (:+,:-,:&,:∧,:∨)
@@ -1054,7 +1095,9 @@ function deletepointcloud!(P::Int)
     nothing
 end
 
+⊕(a::PointArray,b::PointArray) = PointArray(points(a)⊕points(b))
 ⊕(a::PointArray,b::AbstractVector{<:Real}) = PointArray(points(a)⊕b)
+cross(a::PointArray,b::PointArray) = a⊕b
 cross(a::PointArray,b::AbstractVector{<:Real}) = a⊕b
 
 Base.size(m::PointArray) = size(m.dom)
@@ -1183,6 +1226,7 @@ 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)); m)
 resize_lastdim!(m::GlobalFiber,i) = ((resize_lastdim!(domain(m),i),resize_lastdim!(codomain(m),i)); m)
+
 # AbstractFrameBundle
 
 export AbstractFrameBundle, GridFrameBundle, SimplexFrameBundle, FacetFrameBundle
@@ -1190,8 +1234,22 @@ export IntervalRange, AlignedRegion, AlignedSpace
 
 abstract type AbstractFrameBundle{M,N} <: GlobalFiber{M,N} end
 
+base(m::AbstractFrameBundle) = m.p
+immersion(m::AbstractFrameBundle) = m.t
+coordinates(m::AbstractFrameBundle) = m
+points(m::AbstractFrameBundle) = points(base(m))
+fiber(m::AbstractFrameBundle) = fiber(base(m))
+metricextensor(m::AbstractFrameBundle) = metricextensor(base(m))
+metrictensor(m::AbstractFrameBundle) = metrictensor(base(m))
+pointtype(m::AbstractFrameBundle) = basetype(m)
+pointtype(m::Type{<:AbstractFrameBundle}) = basetype(m)
+metrictype(m::AbstractFrameBundle) = fibertype(m)
+metrictype(m::Type{<:AbstractFrameBundle}) = fibertype(m)
 Base.size(m::AbstractFrameBundle) = size(points(m))
 
+@pure isbundle(::AbstractFrameBundle) = true
+@pure isbundle(t) = false
+
 @pure Grassmann.Manifold(m::AbstractFrameBundle) = Manifold(points(m))
 @pure LinearAlgebra.rank(m::AbstractFrameBundle) = rank(points(m))
 @pure Grassmann.grade(m::AbstractFrameBundle) = grade(points(m))
@@ -1202,20 +1260,20 @@ Base.size(m::AbstractFrameBundle) = size(points(m))
 
 grid_id = 0
 
-struct GridFrameBundle{C<:Coordinate,N,PA<:FiberBundle{C,N},TA<:ImmersedTopology} <: AbstractFrameBundle{C,N}
+struct GridFrameBundle{N,C<:Coordinate,PA<:FiberBundle{C,N},TA<:ImmersedTopology} <: AbstractFrameBundle{C,N}
     id::Int
-    dom::PA
-    cod::TA
-    GridFrameBundle(id::Int,p::PA,t::TA=OpenTopology(size(p))) where {C<:Coordinate,N,PA<:FiberBundle{C,N},TA<:ImmersedTopology} = new{C,N,PA,TA}(id,p,t)
-    GridFrameBundle(p::PA,t::TA=OpenTopology(size(p))) where {C<:Coordinate,N,PA<:FiberBundle{C,N},TA<:ImmersedTopology} = new{C,N,PA,TA}((global grid_id+=1),p,t)
+    p::PA
+    t::TA
+    GridFrameBundle(id::Int,p::PA,t::TA=OpenTopology(size(p))) where {N,C<:Coordinate,PA<:FiberBundle{C,N},TA<:ImmersedTopology} = new{N,C,PA,TA}(id,p,t)
+    GridFrameBundle(p::PA,t::TA=OpenTopology(size(p))) where {N,C<:Coordinate,PA<:FiberBundle{C,N},TA<:ImmersedTopology} = new{N,C,PA,TA}((global grid_id+=1),p,t)
 end
 
-GridFrameBundle(id::Int,p::PA,g::GA) where {P,G,N,PA<:AbstractArray{P,N},GA<:AbstractArray{G,N}} = GridFrameBundle(id,PointArray(0,p,g))
-GridFrameBundle(p::PA,g::GA) where {P,G,N,PA<:AbstractArray{P,N},GA<:AbstractArray{G,N}} = GridFrameBundle((global grid_id+=1),p,g)
+GridFrameBundle(id::Int,p::PA,g::GA) where {N,P,G,PA<:AbstractArray{P,N},GA<:AbstractArray{G,N}} = GridFrameBundle(id,PointArray(0,p,g))
+GridFrameBundle(p::PA,g::GA) where {N,P,G,PA<:AbstractArray{P,N},GA<:AbstractArray{G,N}} = GridFrameBundle((global grid_id+=1),p,g)
 
-const IntervalRange{P<:Real,G,PA<:AbstractRange,GA} = GridFrameBundle{Coordinate{P,G},1,<:PointVector{P,G,PA,GA}}
-const AlignedRegion{N,P<:Chain,G<:InducedMetric,PA<:RealRegion{V,<:Real,N,<:AbstractRange} where V,GA<:Global} = GridFrameBundle{Coordinate{P,G},N,PointArray{P,G,N,PA,GA}}
-const AlignedSpace{N,P<:Chain,G<:InducedMetric,PA<:RealRegion{V,<:Real,N,<:AbstractRange} where V,GA} = GridFrameBundle{Coordinate{P,G},N,PointArray{P,G,N,PA,GA}}
+const IntervalRange{P<:Real,G,PA<:AbstractRange,GA} = GridFrameBundle{1,Coordinate{P,G},<:PointVector{P,G,PA,GA}}
+const AlignedRegion{N,P<:Chain,G<:InducedMetric,PA<:RealRegion{V,<:Real,N,<:AbstractRange} where V,GA<:Global} = GridFrameBundle{N,Coordinate{P,G},PointArray{P,G,N,PA,GA}}
+const AlignedSpace{N,P<:Chain,G<:InducedMetric,PA<:RealRegion{V,<:Real,N,<:AbstractRange} where V,GA} = GridFrameBundle{N,Coordinate{P,G},PointArray{P,G,N,PA,GA}}
 
 GridFrameBundle(id::Int,p::PA) where {N,P,PA<:AbstractArray{P,N}} = GridFrameBundle(id,PointArray(0,p,Global{N}(InducedMetric())))
 GridFrameBundle(p::PA) where {N,P,PA<:AbstractArray{P,N}} = GridFrameBundle(p,Global{N}(InducedMetric()))
@@ -1224,44 +1282,52 @@ GridFrameBundle(dom::GridFrameBundle,fun::Array) = GridFrameBundle(base(dom), fu
 GridFrameBundle(dom::GridFrameBundle,fun::Function) = GridFrameBundle(base(dom), fun)
 GridFrameBundle(dom::AbstractArray,fun::Function) = GridFrameBundle(dom, fun.(dom))
 
-base(m::GridFrameBundle) = m.dom
-points(m::GridFrameBundle) = points(base(m))
-ImmersedTopology(m::GridFrameBundle) = m.cod
-metricextensor(m::GridFrameBundle) = metricextensor(base(m))
-metrictensor(m::GridFrameBundle) = metrictensor(base(m))
-coordinates(t::GridFrameBundle) = t
-pointtype(m::GridFrameBundle) = basetype(m)
-pointtype(m::Type{<:GridFrameBundle}) = basetype(m)
-metrictype(m::GridFrameBundle) = fibertype(m)
-metrictype(m::Type{<:GridFrameBundle}) = fibertype(m)
+isopen(t::GridFrameBundle) = isopen(immersion(t))
+iscompact(t::GridFrameBundle) = iscompact(immersion(t))
+bundle(m::GridFrameBundle) = m.id
 basetype(m::GridFrameBundle) = basetype(base(m))
-basetype(::Type{<:GridFrameBundle{C}}) where C = basetype(C)
+basetype(::Type{<:GridFrameBundle{N,C} where N}) where C = basetype(C)
 fibertype(m::GridFrameBundle) = fibertype(base(m))
-fibertype(::Type{<:GridFrameBundle{C}}) where C = fibertype(C)
+fibertype(::Type{<:GridFrameBundle{N,C} where N}) where C = fibertype(C)
 
+⊕(a::GridFrameBundle{1},b::GridFrameBundle{1}) = GridFrameBundle(base(a)⊕base(b),immersion(a) × immersion(b))
 ⊕(a::GridFrameBundle,b::AbstractVector{<:Real}) = GridFrameBundle(base(a)⊕b,immersion(a) × length(b))
+cross(a::GridFrameBundle,b::GridFrameBundle) = a⊕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]))
+resize(m::GridFrameBundle) = GridFrameBundle(m.id,base(m),resize(immersion(m),size(base(m))[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)
 (m::GridFrameBundle)(i::Vararg{Union{Int,Colon}}) = GridFrameBundle(0,base(m)(i...),immersion(m)(i...))
-@pure Base.eltype(::Type{<:GridFrameBundle{C}}) where C = C
+@pure Base.eltype(::Type{<:GridFrameBundle{N,C} where N}) where C = C
 Base.getindex(m::GridFrameBundle,i::Vararg{Int}) = getindex(base(m),i...)
 Base.getindex(m::GridFrameBundle,i::Vararg{Union{Int,Colon}}) = GridFrameBundle(0,getindex(base(m),i...),immersion(m)(i...))
 Base.setindex!(m::GridFrameBundle,s,i::Vararg{Int}) = setindex!(base(m),s,i...)
 
-Base.BroadcastStyle(::Type{<:GridFrameBundle{C,N,PA}}) where {C,N,PA} = Broadcast.ArrayStyle{GridFrameBundle{C,N,PA}}()
+export Grid
+const Grid = GridFrameBundle
 
-function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{GridFrameBundle{C,N,PA}}}, ::Type{ElType}) where {C,N,PA,ElType<:Coordinate}
+#Grid(v::A,t::I=OpenTopology(size(v))) where {N,T,A<:AbstractArray{T,N},I} = GridFrameBundle(0,PointArray(0,v),t)
+
+Base.getindex(g::GridFrameBundle{M,C,<:FiberBundle,<:OpenTopology} where C,j::Int,n::Val,i::Vararg{Int}) where M = getpoint(g,j,n,i...)
+@generated function getpoint(g::GridFrameBundle{M,C,<:FiberBundle} where C,j::Int,::Val{N},i::Vararg{Int}) where {M,N}
+    :(Base.getindex(points(g),$([k≠N ? :(@inbounds i[$k]) : :(@inbounds i[$k]+j) for k ∈ list(1,M)]...)))
+end
+@generated function Base.getindex(g::GridFrameBundle{M},j::Int,n::Val{N},i::Vararg{Int}) where {M,N}
+    :(Base.getindex(points(g),Base.getindex(immersion(g),n,$([k≠N ? :(@inbounds i[$k]) : :(@inbounds i[$k]+j) for k ∈ list(1,M)]...))...))
+end
+
+Base.BroadcastStyle(::Type{<:GridFrameBundle{N,C,PA}}) where {N,C,PA} = Broadcast.ArrayStyle{GridFrameBundle{N,C,PA}}()
+
+function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{GridFrameBundle{N,C,PA}}}, ::Type{ElType}) where {N,C,PA,ElType<:Coordinate}
     ax = axes(bc)
     # Use the data type to create the output
     GridFrameBundle(similar(Array{pointtype(ElType),N}, ax), similar(Array{metrictype(ElType),N}, ax))
 end
-function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{GridFrameBundle{C,N,PA}}}, ::Type{ElType}) where {C,N,PA,ElType}
+function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{GridFrameBundle{N,C,PA}}}, ::Type{ElType}) where {N,C,PA,ElType}
     t = find_gf(bc)
     # Use the data type to create the output
     GridFrameBundle(similar(Array{ElType,N}, axes(bc)), metricextensor(t))
@@ -1286,19 +1352,13 @@ end
 
 SimplexFrameBundle(id::Int,p,t,g) = SimplexFrameBundle(PointCloud(id,p,g),t)
 SimplexFrameBundle(id::Int,p,t) = SimplexFrameBundle(PointCloud(id,p),t)
+SimplexFrameBundle(p::P,t,g::G) where {P<:AbstractVector,G<:AbstractVector} = SimplexFrameBundle(PointCloud(p,g),t)
 #SimplexFrameBundle(p::AbstractVector,t) = SimplexFrameBundle(PointCloud(p),t)
 
 (p::PointCloud)(t::ImmersedTopology) = SimplexFrameBundle(p,t)
 PointCloud(m::SimplexFrameBundle) = m.p
-ImmersedTopology(m::SimplexFrameBundle) = m.t
-coordinates(t::SimplexFrameBundle) = t
-points(m::SimplexFrameBundle) = points(PointCloud(m))
-metricextensor(m::SimplexFrameBundle) = metricextensor(PointCloud(m))
-metrictensor(m::SimplexFrameBundle) = metrictensor(PointCloud(m))
-pointtype(m::SimplexFrameBundle) = basetype(m)
-pointtype(m::Type{<:SimplexFrameBundle}) = basetype(m)
-metrictype(m::SimplexFrameBundle) = fibertype(m)
-metrictype(m::Type{<:SimplexFrameBundle}) = fibertype(m)
+bundle(m::SimplexFrameBundle) = m.id
+deletebundle!(m::SimplexFrameBundle) = deletepointcloud!(bundle(m))
 basetype(::SimplexFrameBundle{P}) where P = pointtype(P)
 basetype(::Type{<:SimplexFrameBundle{P}}) where P = pointtype(P)
 fibertype(::SimplexFrameBundle{P}) where P = metrictype(P)
@@ -1307,6 +1367,11 @@ Base.size(m::SimplexFrameBundle) = size(vertices(m))
 
 Base.broadcast(f,t::SimplexFrameBundle) = SimplexFrameBundle(f.(PointCloud(t)),ImmersedTopology(t))
 
+Base.firstindex(m::SimplexFrameBundle) = 1
+Base.lastindex(m::SimplexFrameBundle) = length(vertices(m))
+Base.length(m::SimplexFrameBundle) = length(vertices(m))
+#Base.resize!(m::SimplexFrameBundle,n::Int) = resize!(value(m),n)
+
 (m::SimplexFrameBundle)(i::ImmersedTopology) = SimplexFrameBundle(bundle(m),PointCloud(m),i)
 Base.getindex(m::SimplexFrameBundle,i::Chain{V,1}) where V = Chain{Manifold(V),1}(points(m)[value(i)])
 Base.getindex(m::SimplexFrameBundle,i::Values{N,Int}) where N = points(m)[value(i)]
@@ -1344,26 +1409,6 @@ end
 #FacetFrameBundle(id::Int,p,t) = FacetFrameBundle(PointCloud(id,p),t)
 #FacetFrameBundle(p::AbstractVector,t) = FacetFrameBundle(PointCloud(p),t)
 
-PointCloud(m::FacetFrameBundle) = m.p
-ImmersedTopology(m::FacetFrameBundle) = m.t
-coordinates(t::FacetFrameBundle) = t
-points(m::FacetFrameBundle) = points(PointCloud(m))
-metricextensor(m::FacetFrameBundle) = metricextensor(PointCloud(m))
-metrictensor(m::FacetFrameBundle) = metrictensor(PointCloud(m))
-pointtype(m::FacetFrameBundle) = basetype(m)
-pointtype(m::Type{<:FacetFrameBundle}) = basetype(m)
-metrictype(m::FacetFrameBundle) = fibertype(m)
-metrictype(m::Type{<:FacetFrameBundle}) = fibertype(m)
-basetype(::FacetFrameBundle{P}) where P = pointtype(P)
-basetype(::Type{<:FacetFrameBundle{P}}) where P = pointtype(P)
-fibertype(::FacetFrameBundle{P}) where P = metrictype(P)
-fibertype(::Type{<:FacetFrameBundle{P}}) where P = metrictype(P)
-
-Base.broadcast(f,t::FacetFrameBundle) = FacetFrameBundle(0,f.(PointCloud(t)),ImmersedTopology(t))
-
-function SimplexFrameBundle(p::P,t,g::G) where {P<:AbstractVector,G<:AbstractVector}
-    SimplexFrameBundle(PointCloud(p,g),t)
-end
 function SimplexFrameBundle(m::FacetFrameBundle)
     SimplexFrameBundle(PointCloud(m.id,point_cache[m.id],point_metric_cache[m.id]),ImmersedTopology(m))
 end
@@ -1372,23 +1417,21 @@ function FacetFrameBundle(m::SimplexFrameBundle)
     FacetFrameBundle(m.id,PointCloud(0,barycenter.(m[et]),barycenter.(getindex.(Ref(metricextensor(m)),et))),ImmersedTopology(m))
 end
 
-bundle(m::GridFrameBundle) = m.id
-bundle(m::SimplexFrameBundle) = m.id
+PointCloud(m::FacetFrameBundle) = m.p
+basetype(::FacetFrameBundle{P}) where P = pointtype(P)
+basetype(::Type{<:FacetFrameBundle{P}}) where P = pointtype(P)
+fibertype(::FacetFrameBundle{P}) where P = metrictype(P)
+fibertype(::Type{<:FacetFrameBundle{P}}) where P = metrictype(P)
+
+Base.broadcast(f,t::FacetFrameBundle) = FacetFrameBundle(0,f.(PointCloud(t)),ImmersedTopology(t))
+
 bundle(m::FacetFrameBundle) = m.id
-deletebundle!(m::SimplexFrameBundle) = deletepointcloud!(bundle(m))
 deletebundle!(m::FacetFrameBundle) = deletepointcloud!(bundle(m))
-@pure isbundle(::AbstractFrameBundle) = true
-@pure isbundle(t) = false
 #@pure ispoints(t::Submanifold{V}) where V = isbundle(V) && rank(V) == 1 && !isbundle(Manifold(V))
 #@pure ispoints(t) = isbundle(t) && rank(t) == 1 && !isbundle(Manifold(t))
 #@pure islocal(t) = isbundle(t) && rank(t)==1 && valuetype(t)==Int && ispoints(Manifold(t))
 #@pure iscell(t) = isbundle(t) && islocal(Manifold(t))
 
-Base.firstindex(m::SimplexFrameBundle) = 1
-Base.lastindex(m::SimplexFrameBundle) = length(vertices(m))
-Base.length(m::SimplexFrameBundle) = length(vertices(m))
-#Base.resize!(m::SimplexFrameBundle,n::Int) = resize!(value(m),n)
-
 @pure Base.eltype(::Type{<:FacetFrameBundle{P,G}}) where {P,G} = Coordinate{P,G}
 function Base.getindex(m::FacetFrameBundle,i::Int)
     ind = getimage(m,i)
@@ -1406,9 +1449,9 @@ end
 # FiberProductBundle
 
 struct FiberProductBundle{P,N,SA<:AbstractArray,PA<:AbstractArray} <: AbstractFrameBundle{Coordinate{P,InducedMetric},N}
-    dom::SA
-    cod::PA
-    FiberProductBundle{P}(s::SA,p::PA) where {P,M,N,SA<:AbstractArray{S,M} where S,PA<:AbstractArray{F,N} where F} = new{P,M+N,SA,PA}(s,p)
+    s::SA
+    g::PA
+    FiberProductBundle{P}(s::SA,g::PA) where {P,M,N,SA<:AbstractArray{S,M} where S,PA<:AbstractArray{F,N} where F} = new{P,M+N,SA,PA}(s,g)
 end
 
 varmanifold(N::Int) = Submanifold(N+1)(list(1,N)...)
@@ -1422,10 +1465,14 @@ function ⊕(a::SimplexFrameBundle,b::AbstractVector{B}) where B<:Real
     FiberProductBundle{P}(a,ProductSpace{W(N)}(Values((b,))))
 end
 
-Base.size(m::FiberProductBundle) = (length(base(m)),size(fiber(m))...)
+basetype(m::FiberProductBundle{P}) where P = basetype(P)
+basetype(::Type{<:FiberProductBundle{P}}) where P = basetype(P)
+fibertype(m::FiberProductBundle{P}) where P = fibertype(P)
+fibertype(::Type{<:FiberProductBundle{P}}) where P = fibertype(P)
+Base.size(m::FiberProductBundle) = (length(m.s),size(m.g)...)
 
 @pure Base.eltype(::Type{<:FiberProductBundle{P}}) where P = Coordinate{P,InducedMetric}
-Base.getindex(m::FiberProductBundle,i::Int,j::Vararg{Int}) = Coordinate(getindex(points(base(m)),i) ⧺ getindex(fiber(m),j...), InducedMetric())
+Base.getindex(m::FiberProductBundle,i::Int,j::Vararg{Int}) = Coordinate(getindex(points(m.s),i) ⧺ getindex(m.g,j...), InducedMetric())
 #=Base.setindex!(m::FiberProductBundle{P},s::P,i::Int,j::Vararg{Int}) where P = setindex!(points(m),s,i,j...)
 Base.setindex!(m::FiberProductBundle{P,G} where P,s::G,i::Int,j::Vararg{Int}) where G = setindex!(metricextensor(m),s,i,j...)
 function Base.setindex!(m::FiberProductBundle,s::Coordinate,i::Int,j::Vararg{Int})
@@ -1444,14 +1491,11 @@ end
 
 (p::FiberProduct)(t::ImmersedTopology) = HomotopyBundle(p,t)
 FiberProduct(m::HomotopyBundle) = m.p
-ImmersedTopology(m::HomotopyBundle) = m.t
-coordinates(t::HomotopyBundle) = t
-points(m::HomotopyBundle) = points(FiberProduct(m))
-fiberspace(m::HomotopyBundle) = fiberspace(FiberProduct(m))
-pointtype(m::HomotopyBundle) = basetype(m)
-pointtype(m::Type{<:HomotopyBundle}) = basetype(m)
+fiberspace(m::HomotopyBundle) = fiberspace(base(m))
 basetype(::HomotopyBundle{P}) where P = pointtype(P)
 basetype(::Type{<:HomotopyBundle{P}}) where P = pointtype(P)
+fibertype(::HomotopyBundle{P}) where P = fibertype(P)
+fibertype(::Type{<:HomotopyBundle{P}}) where P = fibertype(P)
 Base.size(m::HomotopyBundle) = size(FiberProduct(m))
 
 Base.broadcast(f,t::HomotopyBundle) = HomotopyBundle(f.(FiberProduct(t)),ImmersedTopology(t))