deprecated StaticArrays

chakravala · Aug 28, 2020 · 31657a8 · 31657a8
1 parent 10090dd
commit 31657a8
Show file tree

Hide file tree

Showing 5 changed files with 121 additions and 98 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,12 +1,11 @@
 name = "Adapode"
 uuid = "0750cfb5-909a-49d7-927e-29b6595444bf"
 authors = ["Michael Reed"]
-version = "0.2.4"
+version = "0.2.5"
 
 [deps]
 DirectSum = "22fd7b30-a8c0-5bf2-aabe-97783860d07c"
 Grassmann = "4df31cd9-4c27-5bea-88d0-e6a7146666d8"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 AbstractTensors = "a8e43f4a-99b7-5565-8bf1-0165161caaea"
@@ -16,7 +15,6 @@ julia = "1"
 AbstractTensors = "0.5"
 DirectSum = "0.7"
 Grassmann = "0.6"
-StaticArrays = "0"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

diff --git a/src/Adapode.jl b/src/Adapode.jl
@@ -3,13 +3,13 @@ module Adapode
 #   This file is part of Aadapode.jl. It is licensed under the GPL license
 #   Adapode Copyright (C) 2019 Michael Reed
 
-using StaticArrays, SparseArrays, LinearAlgebra
+using SparseArrays, LinearAlgebra
 using AbstractTensors, DirectSum, Grassmann
 import Grassmann: value, vector, valuetype
 import Base: @pure
-import AbstractTensors: SVector, MVector, SizedVector
+import AbstractTensors: Values, Variables, FixedVector
 
-export SVector, odesolve
+export Values, odesolve
 export initmesh, pdegrad
 
 include("constants.jl")
@@ -42,7 +42,7 @@ function weights(c,fx)
     end
     return cfx
 end
-@pure shift(::Val{m},::Val{k}=Val(1)) where {m,k} = SVector{m,Int}(k:m+k-1)
+@pure shift(::Val{m},::Val{k}=Val(1)) where {m,k} = Values{m,Int}(k:m+k-1)
 @pure shift(M::Val{m},i) where {m,a} = ((shift(M,Val{0}()).+i).%(m+1)).+1
 explicit(x,h,c,fx) = (l=length(c);x+weights(h*c,l≠length(fx) ? fx[shift(Val(l))] : fx))
 explicit(x,h,c,fx,i) = (l=length(c);weights(h*c,l≠length(fx) ? fx[shift(Val(l),i)] : fx))
@@ -54,7 +54,7 @@ function butcher(x,f,h,v::Val{N}=Val(4),a::Val{A}=Val(false)) where {N,A}
     b = butcher(v,a)
     n = length(b)-A
     T = typeof(x)
-    fx = T<:Vector ? SizedVector{n,T}(undef) : MVector{n,T}(undef)
+    fx = T<:Vector ? FixedVector{n,T}(undef) : Variables{n,T}(undef)
     @inbounds fx[1] = f(x)
     for k ∈ 2:n
         @inbounds fx[k] = f(explicit(x,h,b[k-1],fx))
@@ -87,7 +87,7 @@ function initsteps(x0,t,tmin,tmax)
     return x
 end
 function initsteps(x0,t,tmin,tmax,f,m::Val{o},B=Val(4)) where o
-    initsteps(x0,t,tmin,tmax), MVector{o+1,typeof(x0)}(undef)
+    initsteps(x0,t,tmin,tmax), Variables{o+1,typeof(x0)}(undef)
 end
 
 function initsteps!(x,f,fx,t,B=Val(4))

diff --git a/src/constants.jl b/src/constants.jl
@@ -2,70 +2,98 @@
 #   This file is part of Adapode.jl. It is licensed under the GPL license
 #   Adapode Copyright (C) 2019 Michael Reed
 
-const CB = SVector( # Fixed
-    SVector( # Euler
-        SVector{0,Float64}(),
-        SVector(1)),
-    SVector( # midpoint
-        SVector(0.5),
-        SVector(0,1)),
-    SVector( # Kutta
-        SVector(0.5),
-        SVector(-1,2),
-        SVector((1,4,1)./6)),
-    SVector( # Runge-Kutta
-        SVector(0.5),
-        SVector(0,0.5),
-        SVector(0,0,1),
-        SVector((1,2,2,1)./6)))
+const CB = Values( # Fixed
+    Values( # Euler
+        Values{0,Float64}(),
+        Values(1)),
+    Values( # midpoint
+        Values(1/2),
+        Values(0,1)),
+    Values( # Kutta
+        Values(1/2),
+        Values(-1,2),
+        Values((1,4,1)./6)),
+    Values( # Runge-Kutta
+        Values(0.5),
+        Values(0,0.5),
+        Values(0,0,1),
+        Values((1,2,2,1)./6)))
 
-constants(a,b,c) = SVector(a...,b,b-c)
+constants(a,b,c) = Values(a...,b,b-c)
 
-const CBA = SVector( # Adaptive
-    constants(SVector( # Heun-Euler
-            SVector(1)),
-        SVector(1,0),
-        SVector(0.5,0.5)),
-    constants(SVector( # Bogacki-Shampine
-            SVector(0.5),
-            SVector(0,0.75),
-            SVector(2/9,1/3,4/9)),
-        SVector(7/24, 1/4, 1/3, 1/8),
-        SVector(2/9,1/3,4/9,0)),
-    constants(SVector( # Fehlberg
-            SVector(0.25),
-            SVector(0.09375,0.28125),
-            SVector(1932/2197,-7200/2197,7296/2197),
-            SVector(439/216,-8,3680/513,-845/4104),
-            SVector(-8/27,2,3544/2565,1859/4104,-11/40)),
-        SVector(16/135,0,6656/12825,28561/56430,-0.18,2/55),
-        SVector(25/216,0,1408/2565,2197/4104,-0.2,0)),
-    constants(SVector( # Cash-Karp
-            SVector(0.2),
-            SVector(3/40,9/40),
-            SVector(3/40,-9/10,6/5),
-            SVector(-11/54,5/2,-70/27,35/27),
-            SVector(1631/55296,175/512,575/13824,44275/110592,253/4096)),
-        SVector(2825/27648,0,18575/48384,13525/55296,277/14336,0.25),
-        SVector(37/378,0,250/621,125/594,0,512/1771)),
-    constants(SVector( # Dormand-Prince
-            SVector(0.2),
-            SVector(3/40,9/40),
-            SVector(44/45,-56/15,32/9),
-            SVector(19372/6561,-25360/2187,64448/6561,-212/729),
-            SVector(9017/3168,-355/33,46732/5247,49/176,-5103/18656),
-            SVector(35/384,0,500/1113,125/192,-2187/6784,11/84)),
-        SVector(35/384,0,500/1113,125/192,-2187/6784,11/84,0),
-        SVector(5179/57600,0,7571/16695,393/640,-92097/339200,187/2100,1/40)))
+const CBA = Values( # Adaptive
+    constants(Values( # Heun-Euler
+            Values(1)),
+        Values(1,0),
+        Values(1,1)./2),
+    constants(Values( # Bogacki-Shampine
+            Values(1/2),
+            Values(0,3/4),
+            Values(2,3,4)./9),
+        Values(7/24, 1/4, 1/3, 1/8),
+        Values(2/9,1/3,4/9,0)),
+    constants(Values( # Fehlberg
+            Values(1/4),
+            Values(0.09375,0.28125),
+            Values(1932/2197,-7200/2197,7296/2197),
+            Values(439/216,-8,3680/513,-845/4104),
+            Values(-8/27,2,3544/2565,1859/4104,-11/40)),
+        Values(16/135,0,6656/12825,28561/56430,-0.18,2/55),
+        Values(25/216,0,1408/2565,2197/4104,-0.2,0)),
+    constants(Values( # Cash-Karp
+            Values(1/5),
+            Values(3/40,9/40),
+            Values(3/40,-9/10,6/5),
+            Values(-11/54,5/2,-70/27,35/27),
+            Values(1631/55296,175/512,575/13824,44275/110592,253/4096)),
+        Values(2825/27648,0,18575/48384,13525/55296,277/14336,0.25),
+        Values(37/378,0,250/621,125/594,0,512/1771)),
+    constants(Values( # Dormand-Prince
+            Values(1/5),
+            Values(3/40,9/40),
+            Values(44/45,-56/15,32/9),
+            Values(19372/6561,-25360/2187,64448/6561,-212/729),
+            Values(9017/3168,-355/33,46732/5247,49/176,-5103/18656),
+            Values(35/384,0,500/1113,125/192,-2187/6784,11/84)),
+        Values(35/384,0,500/1113,125/192,-2187/6784,11/84,0),
+        Values(5179/57600,0,7571/16695,393/640,-92097/339200,187/2100,1/40)))
 
-const CAB = SVector(SVector(1), # Adams-Bashforth
-    SVector(-1,3)./2,
-    SVector(5,-16,23)./12,
-    SVector(-9,37,-59,55)./24,
-    SVector(251,-1274,2616,-2774,1901)./720)
+const CAB = Values(Values(1), # Adams-Bashforth
+    Values(-1,3)./2,
+    Values(5,-16,23)./12,
+    Values(-9,37,-59,55)./24,
+    Values(251,-1274,2616,-2774,1901)./720)
 
-const CAM = SVector(SVector(1), # Adams-Moulton
-    SVector(1,1)./2,
-    SVector(-1,8,5)./12,
-    SVector(1,-5,19,9)./24,
-    SVector(-19,106,-264,646,251)./720)
+const CAM = Values(Values(1), # Adams-Moulton
+    Values(1,1)./2,
+    Values(-1,8,5)./12,
+    Values(1,-5,19,9)./24,
+    Values(-19,106,-264,646,251)./720)
+
+const Gauss = Values(
+    Values(
+        Values(1),
+        Values(
+            Values(1/3,1/3))),
+    Values(
+        Values(1,1,1)/3,
+        Values(
+            Values(1/6,1/6),
+            Values(2/3,1/6),
+            Values(1/6,2/3))),
+    Values(
+        Values(-27/48,25/48,25/48),
+        Values(
+            Values(1/3,1/3),
+            Values(1/5,1/5),
+            Values(3/5,1/5),
+            Values(1/5,3/5))),
+    Values(
+        Values(35494641/158896895,35494641/158896895,40960013/372527180,40960013/372527180,40960013/372527180),
+        Values(
+            Values(100320057/224958844,100320057/224958844),
+            Values(100320057/224958844,16300311/150784976),
+            Values(16300311/150784976,100320057/224958844),
+            Values(13196394/144102857,13196394/144102857),
+            Values(13196394/144102857,85438943/104595944),
+            Values(85438943/104595944,13196394/144102857))))
diff --git a/src/element.jl b/src/element.jl
@@ -26,7 +26,6 @@ function gradienthat(t,m=volumes(t))
     N = ndims(Manifold(t))
     if N == 2 #inv.(m)
         V = Manifold(points(t))
-        i = inv.(m)
         c = Chain{↓(V),1}.(inv.(m))
         Chain{V,1}.(-c,c)
     elseif N == 3
@@ -50,7 +49,7 @@ end
 gradientCR(t,m) = gradientCR(gradienthat(t,m))
 gradientCR(g) = gradientCR.(g)
 function gradientCR(g::Chain{V}) where V
-    Chain{V,1}(g.⋅SVector(
+    Chain{V,1}(g.⋅Values(
         Chain{V,1}(-1,1,1),
         Chain{V,1}(1,-1,1),
         Chain{V,1}(1,1,-1)))
@@ -59,8 +58,8 @@ end
 laplacian(t,u,m=volumes(t),g=gradienthat(t,m)) = value.(abs.(gradient(t,u,m,g)))
 gradient(t,u,m=volumes(t),g=gradienthat(t,m)) = [u[value(t[k])]⋅value(g[k]) for k ∈ 1:length(t)]
 
-for T ∈ (:SVector,:MVector)
-    @eval function assemblelocal!(M,mat::SMatrix{N,N},m,tk::$T{N}) where N
+for T ∈ (:Values,:Variables)
+    @eval function assemblelocal!(M,mat,m,tk::$T{N}) where N
         for i ∈ 1:N, j∈ 1:N
             M[tk[i],tk[j]] += mat[i,j]*m
         end
@@ -127,34 +126,31 @@ interp(t,b,d=degrees(t,b)) = assembleload(t,b,d)
 pretni(t,B::SparseMatrixCSC=incidence(t)) = assembleload(t,sparse(B'))
 pretni(t,ut,B=pretni(t)) = B*ut #interp(t,ut,B::SparseMatrixCSC) = B*ut
 
-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 = (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))
 assemblemass(t,m=volumes(t)) = assembleglobal(mass,t,iterpts(t,m))
 
-stiffness(c,g::Float64,::Val{2}) = (cg = c*g^2; SMatrix{2,2,typeof(c)}(cg,-cg,-cg,cg))
-function stiffness(c,g,::Val{N}) where N
-    A = zeros(MMatrix{N,N,typeof(c)})
-    for i ∈ 1:N, j ∈ 1:N
-        A[i,j] = c*(g[i]⋅g[j])[1]
-    end
-    return SMatrix{N,N,typeof(c)}(A)
-end
+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)
 # iterable(means(t),c) # mapping of c.(means(t))
 
-function sonicstiffness(c,g,::Val{N}) where N
+#=function sonicstiffness(c,g,::Val{N}) where N
     A = zeros(MMatrix{N,N,typeof(c)})
     for i ∈ 1:N, j ∈ 1:N
         A[i,j] = c*g[i][1]^2+g[j][2]^2
     end
     return SMatrix{N,N,typeof(c)}(A)
 end
 assemblesonic(t,c=1,m=volumes(t),g=gradienthat(t,m)) = assembleglobal(sonicstiffness,t,m,iterable(c isa Real ? t : means(t),c),g)
-# iterable(means(t),c) # mapping of c.(means(t))
+# iterable(means(t),c) # mapping of c.(means(t))=#
 
-convection(b,g,::Val{N}) where N = ones(SVector{N,Int})*column((b/N).⋅value(g))'
+#convection(b,g,::Val{N}) where N = ones(Values{N,Int})*column((b/N).⋅value(g))'
+convection(b,g,::Val{N}) where N = outer(∇,Chain(column((b/N).⋅value(g))))
 assembleconvection(t,b,m=volumes(t),g=gradienthat(t,m)) = assembleglobal(convection,t,m,b,g)
 
-SD(b,g,::Val) = (x=column(b.⋅value(g));x*x')
+#SD(b,g,::Val) = (x=column(b.⋅value(g));x*x')
+SD(b,g,::Val) = (x=Chain(column(b.⋅value(g)));outer(x,x))
 assembleSD(t,b,m=volumes(t),g=gradienthat(t,m)) = assembleglobal(SD,t,m,b,g)
 
 function assembledivergence(t,m,g)
@@ -268,9 +264,9 @@ end
 @generated function neighbors(A::SparseMatrixCSC,V,tk,k)
     N,F = ndims(Manifold(V)),(x->x>0)
     N1 = Grassmann.list(1,N)
-    x = SVector{N}([Symbol(:x,i) for i ∈ N1])
-    f = SVector{N}([:(findall($F,A[:,tk[$i]])) for i ∈ N1])
-    b = SVector{N}([Expr(:call,:neighbor,:k,x[setdiff(N1,i)]...) for i ∈ N1])
+    x = Values{N}([Symbol(:x,i) for i ∈ N1])
+    f = Values{N}([:(findall($F,A[:,tk[$i]])) for i ∈ N1])
+    b = Values{N}([Expr(:call,:neighbor,:k,x[setdiff(N1,i)]...) for i ∈ N1])
     Expr(:block,Expr(:(=),Expr(:tuple,x...),Expr(:tuple,f...)),
         Expr(:call,:(Chain{V,1}),b...))
 end
@@ -288,8 +284,8 @@ end
 function centroidvectors(t,m=means(t))
     p,nt = points(t),length(t)
     V = Manifold(p)(2,3)
-    c = Vector{SizedVector{3,Chain{V,1,Float64,2}}}(undef,nt)
-    δ = Vector{SizedVector{3,Float64}}(undef,nt)
+    c = Vector{FixedVector{3,Chain{V,1,Float64,2}}}(undef,nt)
+    δ = Vector{FixedVector{3,Float64}}(undef,nt)
     for k ∈ 1:nt
         c[k] = V.(m[k].-p[value(t[k])])
         δ[k] = value.(abs.(c[k]))
@@ -305,7 +301,8 @@ function nedelec(λ,g,v::Val{3})
     m12 = (f[3,1]-f[3,3]-2f[2,1]+f[2,3])/12
     m13 = (f[3,2]-2f[3,1]-f[2,2]+f[2,1])/12
     m23 = (f[1,2]-f[1,1]-2f[3,2]+f[3,1])/12
-    @SMatrix [m11 m12 m13; m12 m22 m23; m13 m23 m33]
+    #@SMatrix [m11 m12 m13; m12 m22 m23; m13 m23 m33]
+    Chain(Chain(m11,m12,m13),Chain(m12,m22,m23),Chain(m13,m23,m33))
 end
 
 function basisnedelec(p)
@@ -319,7 +316,7 @@ end
 function nedelecmean(t,t2e,signs,u)
     base = Grassmann.vectors(t)
     B = revrot.(base,revrot)./column(.∧(value.(base)))
-    N = basisnedelec(SVector(1,1)/3)
+    N = basisnedelec(Values(1,1)/3)
     x,y,z = columns(t2e); X,Y,Z = columns(signs,1,3)
     (u[x].*X).*(B.⋅N[1]) + (u[y].*Y).*(B.⋅N[2]) + (u[z].*Z).*(B.⋅N[3])
 end

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -2,11 +2,11 @@ using Grassmann
 using Adapode, Test
 
 basis"4"; x0 = 10.0v2+10.0v3+10.0v4
-Lorenz(x::Chain{V}) where V = Chain{V,1}(SVector{4,Float64}(
+Lorenz(x::Chain{V}) where V = Chain{V,1}(
 	1.0,
 	10.0(x[3]-x[2]),
 	x[2]*(28.0-x[4])-x[3],
-	x[2]*x[3]-(8/3)*x[4]))
+	x[2]*x[3]-(8/3)*x[4])
 for k ∈ 0:4
     @test (odesolve(Lorenz,x0,0,2π,7,Val(k),Val(4)); true)
 end