upgraded from TensorFields to Cartan

Author: chakravala

Project.toml

@@ -1,20 +1,20 @@
 name = "Adapode"
 uuid = "0750cfb5-909a-49d7-927e-29b6595444bf"
 authors = ["Michael Reed"]
-version = "0.2.9"
+version = "0.3"
 
 [deps]
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+Cartan = "e24af43f-9fd5-4672-9264-a75f3ae40eb2"
 Grassmann = "4df31cd9-4c27-5bea-88d0-e6a7146666d8"
-TensorFields = "86e2b4fd-d9c8-44dc-a03f-e0a387f3b4e6"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [compat]
 julia = "1"
 Requires = "1"
-TensorFields = "0.1"
-Grassmann = "0.8"
+Grassmann = "0.8,0.9"
+Cartan = "0.3"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

README.md

@@ -12,7 +12,7 @@
 [![Gitter](https://badges.gitter.im/Grassmann-jl/community.svg)](https://gitter.im/Grassmann-jl/community?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge)
 [![Build Status](https://travis-ci.org/chakravala/Adapode.jl.svg?branch=master)](https://travis-ci.org/chakravala/Adapode.jl)
 
-This Julia project originally started as a FORTRAN 95 project called [adapode](https://github.com/chakravala/adapode).
+This Julia project originally started as a FORTRAN 95 project called [adapode](https://github.com/chakravala/adapode) and evolved with [Grassmann.jl](https://github.com/chakravala/Grassmann.jl) and [Cartan.jl](https://github.com/chakravala/Cartan.jl).
 
 ```julia
 using Grassmann, Adapode, Makie
@@ -46,61 +46,62 @@ L2Projector(initmesh(0:1/5:1)[3],x->2x[2]*sin(2π*x[2])+3)
 
 Partial differential equations can also be assembled with various additional methods:
 ```julia
-PoissonSolver(p,e,t,c,f,κ,gD=1,gN=0) = mesh(t,color=solvepoisson(t,e,c,f,κ,gD,gN))
 function solvepoisson(t,e,c,f,κ,gD=0,gN=0)
     m = volumes(t)
     b = assemblefunction(t,f,m)
     A = assemblestiffness(t,c,m)
     R,r = assemblerobin(e,κ,gD,gN)
-    return (A+R)\(b+r)
+    return TensorField(t,(A+R)\(b+r))
 end
-function solveSD(t,e,c,f,δ,κ,gD=0,gN=0)
+function solvetransportdiffusion(tf,eκ,c,δ,gD=0,gN=0)
+    t,f,e,κ = base(tf),fiber(tf),base(eκ),fiber(eκ)
     m = volumes(t)
     g = gradienthat(t,m)
     A = assemblestiffness(t,c,m,g)
-    b = means(t,f)
+    b = means(immersion(t),f)
     C = assembleconvection(t,b,m,g)
     Sd = assembleSD(t,sqrt(δ)*b,m,g)
     R,r = assemblerobin(e,κ,gD,gN)
-    return (A+R-C'+Sd)\r
+    return TensorField(t,(A+R-C'+Sd)\r)
 end
 function solvetransport(t,e,c,ϵ=0.1)
     m = volumes(t)
     g = gradienthat(t,m)
     A = assemblestiffness(t,ϵ,m,g)
     b = assembleload(t,m)
     C = assembleconvection(t,c,m,g)
-    return solvedirichlet(A+C,b,e)
+    return TensorField(t,solvedirichlet(A+C,b,e))
 end
 ```
-Such modular methods can work on input meshes of any dimension.
+Such modular methods can work with a `TensorField` of any dimension.
 The following examples are based on trivially generated 1 dimensional domains:
 ```Julia
 function BackwardEulerHeat1D()
     x,m = 0:0.01:1,100; p,e,t = initmesh(x)
     T = range(0,0.5,length=m+1) # time grid
     ξ = 0.5.-abs.(0.5.-x) # initial condition
-    A = assemblestiffness(t) # assemble(t,1,2x)
-    M,b = assemblemassfunction(t,2x).+assemblerobin(e,1e6,0,0)
+    A = assemblestiffness(p(t)) # assemble(p(t),1,2x)
+    M,b = assemblemassfunction(p(t),2x).+assemblerobin(p(e),1e6,0,0)
     h = Float64(T.step); LHS = M+h*A # time step
     for l ∈ 1:m
         ξ = LHS\(M*ξ+h*b); l%10==0 && println(l*h)
     end
-    mesh(t,color=ξ)
+    lines(TensorField(p(t),ξ))
 end
 function PoissonAdaptive(g,p,e,t,c=1,a=0,f=1)
     ϵ = 1.0
+    pt,pe = p(t),p(e)
     while ϵ > 5e-5 && length(t) < 10000
-        m = volumes(t)
-        h = gradienthat(t,m)
-        A,M,b = assemble(t,c,a,f,m,h)
-        ξ = solvedirichlet(A+M,b,e)
-        η = jumps(t,c,a,f,ξ,m,h)
-        display(mesh(t,color=ξ,shading=false))
+        m = volumes(pt)
+        h = gradienthat(pt,m)
+        A,M,b = assemble(pt,c,a,f,m,h)
+        ξ = solvedirichlet(A+M,b,pe)
+        η = jumps(pt,c,a,f,ξ,m,h)
+        display(lines(TensorField(pt,ξ)))
         if typeof(g)<:AbstractRange
-            scatter!(p,ξ,markersize=0.01)
+            #scatter!(p,ξ,markersize=0.01)
         else
-            wireframe!(t,color=(:red,0.6),linewidth=3)
+            #wireframe!(t,color=(:red,0.6),linewidth=3)
         end
         ϵ = sqrt(norm(η)^2/length(η))
         println(t,", ϵ=$ϵ, α=$(ϵ/maximum(η))"); sleep(0.5)
@@ -110,5 +111,5 @@ function PoissonAdaptive(g,p,e,t,c=1,a=0,f=1)
 end
 PoissonAdaptive(refinemesh(0:0.25:1)...,1,0,x->exp(-100abs2(x[2]-0.5)))
 ```
-More general problems for finite element boundary value problems are also enabled by mesh representations imported from external sources.
+More general problems for finite element boundary value problems are also enabled by mesh representations imported from external sources and managed by `Cartan` via `Grassmann` algebra.
 These methods can automatically generalize to higher dimensional manifolds and are compatible with discrete differential geometry.

src/Adapode.jl

@@ -18,7 +18,7 @@ module Adapode
 #  |___|___||_____||___._||   __||_____||_____||_____|
 #                         |__|
 
-using SparseArrays, LinearAlgebra, Grassmann, TensorFields, Requires
+using SparseArrays, LinearAlgebra, Grassmann, Cartan, Requires
 import Grassmann: value, vector, valuetype, tangent, list
 import Grassmann: Values, Variables, FixedVector
 import Grassmann: Scalar, GradedVector, Bivector, Trivector
@@ -78,12 +78,12 @@ function explicit(x,h,c,fx,i)
 end
 function heun(x,f::Function,h)
     hfx = h*f(x)
-    fiber(x)+(hfx+h*f((base(x)+h)↦(fiber(x)+hfx)))/2
+    fiber(x)+(hfx+h*f((point(x)+h)↦(fiber(x)+hfx)))/2
 end
 function heun(x,f::TensorField,t)
     fx = f[t.i]
     hfx = t.h*fx
-    fiber(x)+(hfx+t.h*f(base(fx)↦(fiber(x)+hfx)))/2
+    fiber(x)+(hfx+t.h*f(point(fx)↦(fiber(x)+hfx)))/2
 end
 
 @pure butcher(::Val{N},::Val{A}=Val(false)) where {N,A} = A ? CBA[N] : CB[N]
@@ -94,7 +94,7 @@ function butcher(x::Section{B,F},f,h,v::Val{N}=Val(4),a::Val{A}=Val(false)) wher
     fx = F<:Vector ? FixedVector{n,F}(undef) : Variables{n,F}(undef)
     @inbounds fx[1] = f(x)
     for k ∈ 2:n
-        @inbounds fx[k] = f((base(x)+h*sum(b[k-1]))↦explicit(x,h,b[k-1],fx))
+        @inbounds fx[k] = f((point(x)+h*sum(b[k-1]))↦explicit(x,h,b[k-1],fx))
     end
     return fx
 end
@@ -108,7 +108,7 @@ end # more accurate compared with CAB[k] methods
 function predictcorrect(x,f,fx,t,k::Val{m}=Val(4)) where m
     iszero(t.s) && initsteps!(x,f,fx,t)
     @inbounds xti = x[t.i]
-    xi,tn = fiber(xti),base(xti)+t.h
+    xi,tn = fiber(xti),point(xti)+t.h
     xn = multistep!(xti,f,fx,t,k)
     t.s = (t.s%(m+1))+1; t.i += 1
     xn = multistep!(tn↦(xi+xn),f,fx,t,k,Val(true))
@@ -117,7 +117,7 @@ end
 function predictcorrect(x,f,fx,t,::Val{1})
     @inbounds xti = x[t.i]
     t.i += 1
-    fiber(xti)+t.h*f((base(xti)+t.h)↦(fiber(xti)+t.h*f(xti)))
+    fiber(xti)+t.h*f((point(xti)+t.h)↦(fiber(xti)+t.h*f(xti)))
 end
 
 initsteps(x0,t,tmax,m) = initsteps(init(x0),t,tmax,m)
@@ -142,7 +142,7 @@ end
 function predictcorrect2(x,f,fx,t,k::Val{m}=Val(4)) where m
     iszero(t.s) && initsteps!(x,f,fx,t)
     @inbounds xti = x[t.i]
-    xi,tn = fiber(xti),base(xti)+t.h
+    xi,tn = fiber(xti),point(xti)+t.h
     xn = multistep2!(xti,f,fx,t,k)
     t.s = (t.s%m)+1; t.i += 1
     xn = multistep2!(tn↦(xi+xn),f,fx,t,k,Val(true))
@@ -152,7 +152,7 @@ end
 function predictcorrect2(x,f,fx,t,::Val{1})
     @inbounds xti = x[t.i]
     t.i += 1
-    fiber(xti)+t.h*f((base(xti)+t.h)↦xti)
+    fiber(xti)+t.h*f((point(xti)+t.h)↦xti)
 end
 
 initsteps2(x0,t,tmax,f,m,B) = initsteps2(init(x0),t,tmax,f,m,B)
@@ -165,7 +165,7 @@ function initsteps!(x,f,fx,t,B=Val(4))
     @inbounds xi = x[t.i]
     for j ∈ 1:m
         @inbounds fx[j] = f(xi)
-        xi = (base(xi)+t.h) ↦ explicit(xi,f,t.h,B)
+        xi = (point(xi)+t.h) ↦ explicit(xi,f,t.h,B)
         x[t.i+j] = xi
     end
     t.s = 1+m
@@ -175,18 +175,18 @@ end
 function explicit!(x,f,t,B=Val(5))
     resize!(x,t.i,10000)
     @inbounds xti = x[t.i]
-    xi,tn = fiber(xti)
+    xi = fiber(xti)
     fx,b = butcher(xti,f,t.h,B,Val(true)),butcher(B,Val(true))
     t.i += 1
-    @inbounds x[t.i] = (base(xti)+t.h) ↦ explicit(xti,t.h,b[end-1],fx)
+    @inbounds x[t.i] = (point(xti)+t.h) ↦ explicit(xti,t.h,b[end-1],fx)
     @inbounds t.e = maximum(abs.(t.h*value(b[end]⋅fx)))
 end
 
 function predictcorrect!(x,f,fx,t,B::Val{m}=Val(4)) where m
     resize!(x,t.i+m,10000)
     iszero(t.s) && initsteps!(x,f,fx,t)
     @inbounds xti = x[t.i]
-    xi,tn = fiber(xti),base(xti)+t.h
+    xi,tn = fiber(xti),point(xti)+t.h
     p = xi + multistep!(xti,f,fx,t,B)
     t.s = (t.s%(m+1))+1
     c = xi + multistep!(tn↦p,f,fx,t,B,Val(true))
@@ -196,7 +196,7 @@ function predictcorrect!(x,f,fx,t,B::Val{m}=Val(4)) where m
 end
 function predictcorrect!(x,f,fx,t,::Val{1})
     @inbounds xti = x[t.i]
-    xi,tn = fiber(xti),base(xti)+t.h
+    xi,tn = fiber(xti),point(xti)+t.h
     p = xi + t.h*f(xti)
     c = xi + t.h*f(tn↦p)
     t.i += 1
@@ -298,7 +298,7 @@ function timeloop!(x,t,tmax,::Val{m}=Val(1)) where m
         t.s = 0
     end
     iszero(t.s) && checkstep!(t)
-    d = tmax-time(x[t.i])
+    d = tmax-point(x[t.i])
     d ≤ t.h && (t.h = d)
     done = d ≤ t.hmax
     done && truncate!(x,t.i-1)
@@ -310,6 +310,6 @@ time(x) = x[1]
 Base.resize!(x,i,h) = length(x)<i+1 && resize!(x,i+h)
 truncate!(x,i) = length(x)>i+1 && resize!(x,i)
 
-show_progress(x,t,b) = t.i%75000 == 11 && println(time(x[t.i])," out of ",b)
+show_progress(x,t,b) = t.i%75000 == 11 && println(point(x[t.i])," out of ",b)
 
 end # module

src/element.jl

@@ -16,15 +16,15 @@ export assemble, assembleglobal, assemblestiffness, assembleconvection, assemble
 export assemblemass, assemblefunction, assemblemassfunction, assembledivergence
 export assemblemassincidence, asssemblemassnodes, assemblenodes
 export assembleload, assemblemassload, assemblerobin, edges, edgesindices, neighbors
-export solvepoisson, solveSD, solvetransport, solvedirichlet, adaptpoisson
+export solvepoisson, solvetransport, solvetransportdiffusion, solvedirichlet, adaptpoisson
 export gradienthat, gradientCR, gradient, interp, nedelec, nedelecmean, jumps
 export submesh, detsimplex, iterable, callable, value, edgelengths, laplacian
 export boundary, interior, trilength, trinormals, incidence, degrees
-import Grassmann: norm, column, columns, points, pointset, edges
+import Grassmann: norm, column, columns
 using Base.Threads
 
-import TensorFields: iterpts, iterable, callable, revrot
-import TensorFields: gradienthat, laplacian, gradient, assemblelocal!
+import Cartan: points, pointset, edges, iterpts, iterable, callable, revrot
+import Cartan: gradienthat, laplacian, gradient, assemblelocal!
 
 trilength(rc) = value.(abs.(value(rc)))
 function trinormals(t)
@@ -52,8 +52,15 @@ function assembleglobal(M,t,m::T,c::C,g::F) where {T<:AbstractVector,C<:Abstract
     end
     return A
 end
+function assembleglobal(M,X::SimplexFrameBundle,m::T,c::C,g::F) where {T<:AbstractVector,C<:AbstractVector,F<:AbstractVector}
+    np = length(points(X)); A = spzeros(np,np); t = immersion(X)
+    for k ∈ 1:length(t)
+        assemblelocal!(A,M(c[k],g[k],Val(mdims(Manifold(X)))),m[k],value(t[k]))
+    end
+    return A
+end
 
-import TensorFields: weights, degrees, assembleincidence, incidence
+import Cartan: weights, degrees, assembleincidence, incidence
 
 assemblemassincidence(t,f,m=volumes(t),l=m) = assemblemassincidence(t,iterpts(t,f),iterable(t,m),iterable(t,l))
 function assemblemassincidence(t,f::F,m::V,l::T) where {F<:AbstractVector,V<:AbstractVector,T<:AbstractVector}
@@ -66,14 +73,25 @@ function assemblemassincidence(t,f::F,m::V,l::T) where {F<:AbstractVector,V<:Abs
     end
     return M,b
 end
+function assemblemassincidence(X::SimplexFrameBundle,f::F,m::V,l::T) where {F<:AbstractVector,V<:AbstractVector,T<:AbstractVector}
+    np,n,t = length(points(X)),Val(mdims(Manifold(X))),immersion(X)
+    M,b,v = spzeros(np,np), zeros(np), f
+    for k ∈ 1:length(t)
+        tk = value(t[k])
+        assemblelocal!(M,mass(nothing,nothing,n),m[k],tk)
+        b[tk] .+= v[tk]*l[k]
+    end
+    return M,b
+end
 
 assemblefunction(t,f,m=volumes(t),d=degrees(t,m)) = assembleincidence(t,iterpts(t,f)/mdims(t),m)
 assemblenodes(t,f,m=volumes(t),d=degrees(t,m)) = assembleincidence(t,iterpts(t,f)./d,m)
 assemblemassload(t,m=volumes(t),l=m,d=degrees(t)) = assemblemassincidence(t,inv.(d),m,l)
 assemblemassfunction(t,f,m=volumes(t),l=m) = assemblemassincidence(t,iterpts(t,f)/mdims(t),iterable(t,m),iterable(t,l))
+assemblemassfunction(tf::TensorField,m=volumes(base(tf)),l=m) = assemblemassfunction(base(tf),fiber(tf),m,l)
 assemblemassnodes(t,f,m=volumes(t),l=m,d=degrees(t)) = assemblemassincidence(t,iterpts(t,f)./d,iterable(t,m),iterable(t,l))
 
-import TensorFields: assembleload, interp, pretni
+import Cartan: assembleload, interp, pretni
 
 #mass(a,b,::Val{N}) where N = (ones(SMatrix{N,N,Int})+I)/Int(factorial(N+1)/factorial(N-1))
 mass(a,b,::Val{N}) where N = (x=Submanifold(N)(∇);outer(x,x)+I)/Int(factorial(N+1)/factorial(N-1))
@@ -82,6 +100,8 @@ assemblemass(t,m=volumes(t)) = assembleglobal(mass,t,iterpts(t,m))
 stiffness(c,g::Float64,::Val{2}) = (cg = c*g^2; Chain(Chain(cg,-cg),Chain(-cg,cg)))
 stiffness(c,g,::Val{N}) where N = Chain{Submanifold(N),1}(map.(*,c,value(g).⋅Ref(g)))
 assemblestiffness(t,c=1,m=volumes(t),g=gradienthat(t,m)) = assembleglobal(stiffness,t,m,iterable(c isa Real ? t : means(t),c),g)
+assemblestiffness(t::SimplexFrameBundle,c=1,m=volumes(t),g=gradienthat(t,m)) = assembleglobal(stiffness,t,m,iterable(c isa Real ? immersion(t) : means(t),c),g)
+
 # iterable(means(t),c) # mapping of c.(means(t))
 
 #=function sonicstiffness(c,g,::Val{N}) where N
@@ -133,18 +153,19 @@ function solvepoisson(t,e,c,f,κ,gD=0,gN=0)
     b = assemblenodes(t,f,m)
     A = assemblestiffness(t,c,m)
     R,r = assemblerobin(e,κ,gD,gN)
-    return (A+R)\(b+r)
+    return TensorField(t,(A+R)\(b+r))
 end
 
-function solveSD(t,e,c,f,δ,κ,gD=0,gN=0)
+function solvetransportdiffusion(tf,eκ,c,δ,gD=0,gN=0)
+    t,f,e,κ = base(tf),fiber(tf),base(eκ),fiber(eκ)
     m = volumes(t)
     g = gradienthat(t,m)
     A = assemblestiffness(t,c,m,g)
-    b = means(t,f)
+    b = means(immersion(t),f)
     C = assembleconvection(t,b,m,g)
     Sd = assembleSD(t,sqrt(δ)*b,m,g)
     R,r = assemblerobin(e,κ,gD,gN)
-    return (A+R-C'+Sd)\r
+    return TensorField(t,(A+R-C'+Sd)\r)
 end
 
 function solvetransport(t,e,c,ϵ=0.1)
@@ -153,7 +174,7 @@ function solvetransport(t,e,c,ϵ=0.1)
     A = assemblestiffness(t,ϵ,m,g)
     b = assembleload(t,m)
     C = assembleconvection(t,c,m,g)
-    return solvedirichlet(A+C,b,e)
+    TensorField(t,solvedirichlet(A+C,b,e))
 end
 
 function adaptpoisson(g,p,e,t,c=1,a=0,f=1,κ=1e6,gD=0,gN=0)
@@ -171,8 +192,12 @@ function adaptpoisson(g,p,e,t,c=1,a=0,f=1,κ=1e6,gD=0,gN=0)
     return g,p,e,t
 end
 
-solvedirichlet(M,b,e::ChainBundle) = solvedirichlet(M,b,pointset(e))
-solvedirichlet(M,b,e::ChainBundle,u) = solvedirichlet(M,b,pointset(e),u)
+#solvedirichlet(M,b,e::ChainBundle) = solvedirichlet(M,b,pointset(e))
+#solvedirichlet(M,b,e::ChainBundle,u) = solvedirichlet(M,b,pointset(e),u)
+solvedirichlet(M,b,e::SimplexFrameBundle) = solvedirichlet(M,b,vertices(e))
+solvedirichlet(M,b,e::SimplexManifold) = solvedirichlet(M,b,vertices(e))
+solvedirichlet(M,b,e::SimplexFrameBundle,u) = solvedirichlet(M,b,vertices(e),u)
+solvedirichlet(M,b,e::SimplexManifold,u) = solvedirichlet(M,b,vertices(e),u)
 function solvedirichlet(A,b,fixed,boundary)
     neq = length(b)
     free,ξ = interior(fixed,neq),zeros(eltype(b),neq)
@@ -187,15 +212,15 @@ function solvedirichlet(M,b,fixed)
     return ξ
 end
 
-import TensorFields: interior
+import Cartan: interior
 
 solvehomogenous(e,M,b) = solvedirichlet(M,b,e)
 export solvehomogenous, solveboundary
 const solveboundary = solvedirichlet # deprecate
 const edgelengths = volumes # deprecate
 const boundary = pointset # deprecate
 
-import TensorFields: facesindices, edgesindices, neighbor, neighbors, centroidvectors
+import Cartan: facesindices, edgesindices, neighbor, neighbors, centroidvectors
 
 #=
 facesindices(t,cols=columns(t)) = mdims(t) == 3 ? edgesindices(t,cols) : throw(error())