Merge pull request #12 from ioannisPApapadopoulos/jp/assembly

Jp/assembly

Project.toml

@@ -24,7 +24,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
-AnnuliOrthogonalPolynomials = "0.0.1, 0.0.2"
+AnnuliOrthogonalPolynomials = "0.0.2"
 BandedMatrices = "0.17, 1"
 BlockArrays = "0.16"
 BlockBandedMatrices = "0.12"

src/RadialPiecewisePolynomials.jl

@@ -10,7 +10,7 @@ import ContinuumArrays: Weight, grid, ℵ₁, ℵ₀, @simplify, ProjectionFacto
 import Base: in, axes, getindex, broadcasted, tail, +, -, *, /, \, convert, OneTo, show, summary, ==, oneto, diff
 import SemiclassicalOrthogonalPolynomials: HalfWeighted
 import MultivariateOrthogonalPolynomials: BlockOneTo, ModalInterlace, BlockRange1, Plan, ModalTrav, ZernikeITransform
-import ClassicalOrthogonalPolynomials: checkpoints, ShuffledR2HC, TransformFactorization, ldiv, paddeddata, jacobimatrix, orthogonalityweight, SetindexInterlace, pad, blockedrange
+import ClassicalOrthogonalPolynomials: checkpoints, ShuffledR2HC, TransformFactorization, ldiv, paddeddata, jacobimatrix, orthogonalityweight, SetindexInterlace, pad, blockedrange, Clenshaw, recurrencecoefficients, _p0, colsupport
 import LinearAlgebra: eigvals, eigen, isapprox, SymTridiagonal, norm, factorize
 import AnnuliOrthogonalPolynomials: factorize, ZernikeAnnulusITransform
 import LazyArrays: Vcat

src/annuluselement.jl

@@ -86,6 +86,11 @@ end
 axes(Z::ContinuousZernikeAnnulusElementMode) = (Inclusion(annulus(first(Z.points), last(Z.points))), oneto(∞))
 ==(P::ContinuousZernikeAnnulusElementMode, Q::ContinuousZernikeAnnulusElementMode) = P.points == Q.points && P.m == Q.m && P.j == Q.j
 
+function show(io::IO, C::ContinuousZernikeAnnulusElementMode)
+    points = C.points
+    print(io, "ContinuousZernikeAnnulusElementMode: $points.")
+end
+
 function getindex(Z::ContinuousZernikeAnnulusElementMode{T}, xy::StaticVector{2}, j::Int)::T where {T}
     p = Z.points
     α, β = convert(T, p[1]), convert(T, p[2])
@@ -168,11 +173,11 @@ end
 ###
 
 function _sum_semiclassicaljacobiweight(t::T, a::Number, b::Number, c::Number) where T
-    # (t,a,b,c) = map(big, map(float, (t,a,b,c)))
-    # return convert(T, t^c * beta(1+a,1+b) * _₂F₁general2(1+a,-c,2+a+b,1/t))
+    (t,a,b,c) = map(big, map(float, (t,a,b,c)))
+    return convert(T, t^c * beta(1+a,1+b) * _₂F₁general2(1+a,-c,2+a+b,1/t))
     # Working much better thanks to Timon Gutleb's efforts,
     # see https://github.com/JuliaApproximation/SemiclassicalOrthogonalPolynomials.jl/pull/90
-    convert(T, first(sum.(SemiclassicalJacobiWeight.(t,a,b,c:c))))
+    # convert(T, first(sum.(SemiclassicalJacobiWeight.(t,a,b,c:c))))
 end
 
 @simplify function *(A::QuasiAdjoint{<:Any,<:ContinuousZernikeAnnulusElementMode}, B::ContinuousZernikeAnnulusElementMode)
@@ -209,22 +214,60 @@ end
         m₀ = convert(T,π) / ( t^(one(T) + m) ) * jw
         m₀ = m == 0 ? m₀ : m₀ / T(2)
 
-        M = ApplyArray(*, L₁₁', L₁₁)
+        L = Hcat(Vcat(view(L₁₀,1:2,1:1)[1:2], Zeros{T}(∞)), Vcat(view(L₀₁,1:2,1:1)[1:2], Zeros{T}(∞)), L₁₁)
+        # TODO fix the excess zeros
+        return ApplyArray(*,Diagonal(Fill(β^2*m₀,∞)), ApplyArray(*, L', L))
 
-        a = ApplyArray(*, view(L₁₁',1:2,1:2), view(L₁₀,1:2,1))
-        b = ApplyArray(*, view(L₁₁',1:2,1:2), view(L₀₁,1:2,1))
+        # Old code - for deprecation
+        # M = ApplyArray(*, L₁₁', L₁₁)
 
-        a11 = ApplyArray(*, view(L₁₀',1,1:2)', view(L₁₀,1:2,1))
-        a12 = ApplyArray(*, view(L₁₀',1,1:2)', view(L₀₁,1:2,1))
-        a22 = ApplyArray(*, view(L₀₁',1,1:2)', view(L₀₁,1:2,1))
+        # a = ApplyArray(*, view(L₁₁',1:2,1:2), view(L₁₀,1:2,1))
+        # b = ApplyArray(*, view(L₁₁',1:2,1:2), view(L₀₁,1:2,1))
 
-        # C = Vcat(Hcat(a11[1], a12[1], [a;Zeros{T}(∞)]'), Hcat(a12[1], a22[1], [b;Zeros{T}(∞)]'))
-        C = Vcat(Hcat(a11[1], a12[1], a'), Hcat(a12[1], a22[1], b'))[1:2,1:4]
-        M = Hcat(Vcat(C[1:2,3:4]', Zeros{T}(∞,2)), M)
-        return Vcat(Hcat(β^2*m₀*C, Zeros{T}(2,∞)), β^2*m₀*M)
+        # a11 = ApplyArray(*, view(L₁₀',1,1:2)', view(L₁₀,1:2,1))
+        # a12 = ApplyArray(*, view(L₁₀',1,1:2)', view(L₀₁,1:2,1))
+        # a22 = ApplyArray(*, view(L₀₁',1,1:2)', view(L₀₁,1:2,1))
+
+        # # C = Vcat(Hcat(a11[1], a12[1], [a;Zeros{T}(∞)]'), Hcat(a12[1], a22[1], [b;Zeros{T}(∞)]'))
+        # C = Vcat(Hcat(a11[1], a12[1], a'), Hcat(a12[1], a22[1], b'))[1:2,1:4]
+        # M = Hcat(Vcat(C[1:2,3:4]', Zeros{T}(∞,2)), M)
+        # return Vcat(Hcat(β^2*m₀*C, Zeros{T}(2,∞)), β^2*m₀*M)
     end
 end
 
+@simplify function *(A::QuasiAdjoint{<:Any,<:ContinuousZernikeAnnulusElementMode}, B::BroadcastQuasiMatrix{<:Any, typeof(*), <:Tuple{BroadcastQuasiVector, <:ContinuousZernikeAnnulusElementMode}})
+    λ, C = B.args
+    T = promote_type(eltype(A), eltype(C))
+
+    @assert C.via_Jacobi == false
+    @assert A' == C
+
+    α, β = convert(T, first(C.points)), convert(T, last(C.points))
+    ρ = α / β
+    m = C.m
+
+    t = inv(one(T)-ρ^2)
+    L₁₁, L₀₁, L₁₀ = C.L₁₁, C.L₀₁, C.L₁₀
+
+    jw = C.normalize_constants[2] # _sum_semiclassicaljacobiweight(t,0,0,m)
+    m₀ = convert(T,π) / ( t^(one(T) + m) ) * jw
+    m₀ = m == 0 ? m₀ : m₀ / T(2)
+    L = Hcat(Vcat(view(L₁₀,1:2,1:1)[1:2], Zeros{T}(∞)), Vcat(view(L₀₁,1:2,1:1)[1:2], Zeros{T}(∞)), L₁₁)
+
+    # We need to compute the Jacobi matrix multiplier addition due to the
+    # variable Helmholtz coefficient λ(r²). We expand λ(r²) in chebyshevt
+    # and then use Clenshaw to compute λ(β^2*(I-X)/t) where X is the 
+    # correponding Jacobi matrix for this basis.
+    Tn = chebyshevt(C.points[1]..C.points[2])
+    u = Tn \ λ.f.(axes(Tn,1))
+    X = jacobimatrix(SemiclassicalJacobi(t, 0, 0, m))
+    W = Clenshaw(paddeddata(u), recurrencecoefficients(Tn)..., β^2*(I-X/t), _p0(Tn))
+
+    # TODO fix the excess zeros
+    return ApplyArray(*,Diagonal(Fill(β^2*m₀,∞)), ApplyArray(*, L', ApplyArray(*, W, L)))
+end
+
+
 ###
 # Gradient and L2 inner product of gradient
 ##
@@ -287,11 +330,11 @@ end
     t = inv(one(T)-ρ^2)
     jw = B.C.normalize_constants[1] # _sum_semiclassicaljacobiweight(t,1,1,m)
     m₀ = convert(T,π) / ( t^(3 + m) ) * jw
-    m₀ = m == 0 ? m₀ : m₀ / T(2)
+    m₀ = m == 0 ? m₀ : m₀ / 2
 
     Δ = - m₀ * B.C.D
     if m == 0
-        c = convert(T,2π)*(T(1) - ρ^4)
+        c = convert(T,2π)*(one(T) - ρ^4)
         C = [c -c 4*m₀*(m+1); -c c -4*m₀*(m+1)]
     else
         C = [W010(m, ρ) W_100_010(m, ρ) 4*m₀*(m+1); W_100_010(m, ρ) W100(m,ρ) -4*m₀*(m+1)]

src/finitecontinuous.jl

@@ -7,12 +7,12 @@ struct FiniteContinuousZernike{T, N<:Int, P<:AbstractVector} <: Basis{T}
     via_Jacobi::Bool
 end
 
-function FiniteContinuousZernike{T}(N::Int, points::AbstractVector, Fs::Tuple{Vararg{FiniteContinuousZernikeMode}}, via_Jacobi::Bool) where {T}
+function FiniteContinuousZernike{T}(N::Int, points::AbstractVector, Fs::Tuple{Vararg{FiniteContinuousZernikeMode}}; via_Jacobi::Bool) where {T}
     @assert length(points) > 1 && points == sort(points)
     @assert length(Fs) == 2N-1
     FiniteContinuousZernike{T, Int, typeof(points)}(N, points, Fs, via_Jacobi)
 end
-FiniteContinuousZernike(N::Int, points::AbstractVector, via_Jacobi::Bool=VIA_JACOBI) = FiniteContinuousZernike{Float64}(N, points, Tuple(_getFs(N, points, via_Jacobi)), via_Jacobi)
+FiniteContinuousZernike(N::Int, points::AbstractVector; via_Jacobi::Bool=VIA_JACOBI) = FiniteContinuousZernike{Float64}(N, points, Tuple(_getFs(N, points, via_Jacobi)), via_Jacobi=via_Jacobi)
 
 function axes(Z::FiniteContinuousZernike{T}) where T
     first(Z.points) ≈ 0 && return (Inclusion(last(Z.points)*UnitDisk{T}()), blockedrange(Vcat(length(Z.points)-1, Fill(length(Z.points) - 1, Z.N-2))))
@@ -21,6 +21,10 @@ function axes(Z::FiniteContinuousZernike{T}) where T
 end
 ==(P::FiniteContinuousZernike, Q::FiniteContinuousZernike) = P.N == Q.N && P.points == Q.points
 
+function show(io::IO, F::FiniteContinuousZernike)
+    N, points = F.N, F.points
+    print(io, "FiniteContinuousZernike at degree N=$N and endpoints $points.")
+end
 
 # Matrices for lowering to ZernikeAnnulus(1,1) via
 # the Jacobi matrix. Stable, but probably higher complexity
@@ -150,6 +154,19 @@ end
     [F' * F for F in Fs]
 end
 
+###
+# Weighted L2 inner products
+# Gives out list of weighted mass matrices of correct size
+###
+@simplify function *(A::QuasiAdjoint{<:Any,<:FiniteContinuousZernike}, λB::BroadcastQuasiMatrix{<:Any, typeof(*), <:Tuple{BroadcastQuasiVector, FiniteContinuousZernike}})
+    λ, B = λB.args
+    @assert A' == B
+    Fs = B.Fs
+    xs = [axes(F,1) for F in Fs]
+    [F' * (λ.f.(x) .* F) for (F, x) in zip(Fs, xs)]
+end
+
+
 ###
 # Gradient for constructing weak Laplacian.
 ###

src/finitecontinuousmode.jl

@@ -13,17 +13,17 @@ struct FiniteContinuousZernikeMode{T} <: Basis{T}
     b::Int # Should remove once adaptive expansion has been figured out.
 end
 
-function FiniteContinuousZernikeMode(N::Int, points::AbstractVector{T}, m::Int, j::Int, L₁₁, L₀₁, L₁₀, D, normalize_constants::AbstractVector{<:AbstractVector{<:T}}, via_Jacobi, same_ρs) where T
+function FiniteContinuousZernikeMode(N::Int, points::AbstractVector{T}, m::Int, j::Int, L₁₁, L₀₁, L₁₀, D, normalize_constants::AbstractVector{<:AbstractVector{<:T}}; via_Jacobi, same_ρs::Bool) where T
     @assert points == sort(points)
     @assert m ≥ 0
     @assert m == 0 ? j == 1 : 0 ≤ j ≤ 1
     K = first(points) ≈ 0 ? length(points)-2 : length(points) - 1
-    @assert length(L₁₁) == length(L₀₁) == length(L₁₀) == length(D) == same_ρs ? 1 : K
+    @assert length(L₁₁) == length(L₀₁) == length(L₁₀) == length(D) == (same_ρs ? 1 : K)
     @assert length(normalize_constants) ≥ 2
     FiniteContinuousZernikeMode{T}(N, points, m, j, L₁₁, L₀₁, L₁₀, D, normalize_constants, via_Jacobi, same_ρs, m+2N)
 end
 
-function FiniteContinuousZernikeMode(N::Int, points::AbstractVector{T}, m::Int, j::Int, via_Jacobi::Bool=false, same_ρs::Bool=false) where {T}
+function FiniteContinuousZernikeMode(N::Int, points::AbstractVector{T}, m::Int, j::Int; via_Jacobi::Bool=false, same_ρs::Bool=false) where {T}
     K = length(points)-1
     κ = first(points[1]) ≈ 0 ? 2 : 1
     ρs = []
@@ -48,7 +48,7 @@ function FiniteContinuousZernikeMode(N::Int, points::AbstractVector{T}, m::Int,
     D = (Z .\ (Laplacian.(axes.(Z,1)).*Weighted.(Z)))
     D = NTuple{K+1-κ, AbstractMatrix}([Ds.ops[m+1] for Ds in D])
 
-    FiniteContinuousZernikeMode(N, points, m, j, L₁₁, L₀₁, L₁₀, D, normalize_constants, via_Jacobi, same_ρs)
+    FiniteContinuousZernikeMode(N, points, m, j, L₁₁, L₀₁, L₁₀, D, normalize_constants, via_Jacobi=via_Jacobi, same_ρs=same_ρs)
 end
 
 # FiniteContinuousZernikeMode(points::AbstractVector, m::Int, j::Int, L₁₁, L₀₁, L₁₀, D, b::Int) = FiniteContinuousZernikeMode{Float64}(points, m, j, L₁₁, L₀₁, L₁₀, D, b)
@@ -61,6 +61,11 @@ function axes(Z::FiniteContinuousZernikeMode{T}) where T
 end
 ==(P::FiniteContinuousZernikeMode, Q::FiniteContinuousZernikeMode) = P.points == Q.points && P.m == Q.m && P.j == Q.j && P.b == Q.b
 
+function show(io::IO, F::FiniteContinuousZernikeMode)
+    N, points, m, j = F.N, F.points, F.m, F.j
+    print(io, "FiniteContinuousZernikeMode: N=$N, points=$points, m=$m, j=$j.")
+end
+
 function _getCs(points::AbstractVector{T}, m::Int, j::Int, b::Int, L₁₁, L₀₁, L₁₀, D, normalize_constants, via_Jacobi::Bool, same_ρs::Bool) where T
     K = length(points)-1
     if same_ρs
@@ -138,7 +143,7 @@ function _build_top_left_block(F::FiniteContinuousZernikeMode{T}, Ms, γs::Abstr
     if p ≈ 0
         a = [Ms[1][1,1]]
         for i in 2:K
-            append!(a, diag(Ms[i][1:2,1:2]))
+            append!(a, [Ms[i][1,1]; Ms[i][2,2]])
         end
 
         dv = zeros(T, K)
@@ -154,7 +159,7 @@ function _build_top_left_block(F::FiniteContinuousZernikeMode{T}, Ms, γs::Abstr
     else
         a = []
         for i in 1:K
-            append!(a, diag(Ms[i][1:2,1:2]))
+            append!(a, [Ms[i][1,1]; Ms[i][2,2]])
         end
 
         dv = zeros(T, K+1)
@@ -205,17 +210,13 @@ end
 function _build_trailing_bubbles(F::FiniteContinuousZernikeMode{T}, Ms, N::Int, bs::Int, p::T) where T
     K = length(Ms)
     if p ≈ 0
-        Mn = vcat([Ms[1][2:N-1,2:N-1]], [Ms[i][3:N, 3:N] for i in 2:K])
-    else
-        Mn = [Ms[i][3:N, 3:N] for i in 1:K]
-    end
-    if bs ==  2
-        return [Symmetric(BandedMatrix{T}(0=>view(M, band(0)), 1=>view(M, band(1)), 2=>view(M, band(2)))) for M in Mn]
-    elseif bs == 1
-        return [Symmetric(BandedMatrix{T}(0=>view(M, band(0)), 1=>view(M, band(1)))) for M in Mn]
+        # Mn = vcat([Ms[1][2:N-1,2:N-1]], [Ms[i][3:N, 3:N] for i in 2:K])
+        Mn = vcat([reshape(view(Ms[1],2:N-1,2:N-1)[:], N-2, N-2)], [reshape(view(Ms[i],3:N, 3:N)[:], N-2, N-2) for i in 2:K])
     else
-        error("Are you using _build_mass_trailing_bubbles correctly?")
+        # Mn = [Ms[i][3:N, 3:N] for i in 1:K]
+        Mn = [reshape(view(Ms[i],3:N, 3:N)[:], N-2, N-2) for i in 1:K]
     end
+    return [Symmetric(BandedMatrix{T}([i=>view(M, band(i))[:] for i in 0:bs]...)) for M in Mn]
 end
 
 function _arrow_head_matrix(F::FiniteContinuousZernikeMode, Ms, γs::AbstractArray{T}, N::Int, bs::Int, p::T) where T
@@ -258,6 +259,23 @@ end
     end
 end
 
+@simplify function *(A::QuasiAdjoint{<:Any,<:FiniteContinuousZernikeMode}, B::BroadcastQuasiMatrix{<:Any, typeof(*), <:Tuple{BroadcastQuasiVector, FiniteContinuousZernikeMode}})
+
+    λ, F = B.args
+    N, K = F.N, length(F.points)-1
+    @assert A' == F
+
+    Cs = _getCs(F)
+    xs = [axes(C,1) for C in Cs]
+    Ms = [C' * (λ.f.(x) .* C) for (C, x) in zip(Cs, xs)]
+
+    # figure out necessary bandwidth
+    # bs for number of bubbles
+    bs = min(N-2, maximum([last(colsupport(view(Ms[i], :, 1)[:]))-2 for i in 1:lastindex(Ms)]))
+    γs = _getγs(F)
+    return _arrow_head_matrix(F, Ms, γs, N, bs, first(F.points))
+end
+
 ###
 # Gradient for constructing weak Laplacian.
 ###

src/finitezernike.jl

@@ -146,7 +146,7 @@ end
     T = promote_type(eltype(FT), eltype(Z))
     F = FT.parent
 
-    @assert F.points == Z.points && F.m == Z.m && F.j == Z.j
+    @assert F.points == Z.points && F.m == Z.m && F.j == Z.j && F.via_Jacobi == false
 
     points, a, b, m = Z.points, Z.a, Z.b, Z.m
     α, β = first(points), last(points)

src/plotting/plots.jl

@@ -5,9 +5,9 @@ using PyPlot, DelimitedFiles, LaTeXStrings
 #######
 
 function _plot(K::Int, θs::AbstractVector, rs::AbstractVector, vals::AbstractVector;ρ::T=0.0, ttl=[]) where T
-    PyPlot.rc("font", family="serif", size=14)
-    rcParams = PyPlot.PyDict(PyPlot.matplotlib["rcParams"])
-    rcParams["text.usetex"] = true
+    # PyPlot.rc("font", family="serif", size=14)
+    # rcParams = PyPlot.PyDict(PyPlot.matplotlib["rcParams"])
+    # rcParams["text.usetex"] = true
     fig = plt.figure()
     ax = fig.add_subplot(111, projection="polar")
     vmin,vmax = minimum(minimum.(vals)), maximum(maximum.(vals))