Merge pull request #15 from ioannisPApapadopoulos/jp/assembly

one cell disk now working

src/diskelement.jl

@@ -68,7 +68,7 @@ function ldiv(C::ContinuousZernikeElementMode{T}, f::AbstractQuasiVector) where
     end
 
     # Seems to cache on its own, no need to memoize unlike annulus
-    c = Z0 \ f̃ # Zernike transform
+    c = pad(Z0[:,Block.(1:200)] \ f̃, axes(Z0,2)) # Zernike transform
     c̃ = ModalTrav(paddeddata(c))
     N = size(c̃.matrix, 1) # degree
 

src/finitecontinuous.jl

@@ -79,36 +79,42 @@ function _getFs(N::Int, points::AbstractVector{T}, via_Jacobi::Bool) where T
         append!(ρs, [α / β])
     end
 
-    # If all ρs are the same, then we can reduce the computational
-    # overhead by only considering one hierarchy of semiclassical
-    # Jacobi polynomials for all the cells.
-    if all(ρs .≈ ρs[1])
-        ρs = [ρs[1]]
-        κ = K
-        same_ρs = true
-    end
+    if !isempty(ρs)
+        # If all ρs are the same, then we can reduce the computational
+        # overhead by only considering one hierarchy of semiclassical
+        # Jacobi polynomials for all the cells.
+        if all(ρs .≈ ρs[1])
+            ρs = [ρs[1]]
+            κ = K
+            same_ρs = true
+        end
+
+        # Semiclassical Jacobi parameter t
+        ts = inv.(one(T) .- ρs.^2)
 
-    # Semiclassical Jacobi parameter t
-    ts = inv.(one(T) .- ρs.^2)
+        # Use broadcast notation to compute all the lowering matrices across all
+        # intervals and Fourier modes simultaneously.
 
-    # Use broadcast notation to compute all the lowering matrices across all
-    # intervals and Fourier modes simultaneously.
+        if via_Jacobi
+            Ls = _ann2element_via_Jacobi.(ts)
+            cst = [[sum.(SemiclassicalJacobiWeight.(t,1,1,0:ms[end])) for t in ts],
+                    [sum.(SemiclassicalJacobiWeight.(t,2,0,0:ms[end])) for t in ts],
+                    [sum.(SemiclassicalJacobiWeight.(t,0,2,0:ms[end])) for t in ts]]
+        else
+            Ls = _ann2element_via_lowering.(ts)
+            cst = [[sum.(SemiclassicalJacobiWeight.(t,a,a,0:ms[end])) for t in ts] for a in 1:-1:0]
+        end
 
-    if via_Jacobi
-        Ls = _ann2element_via_Jacobi.(ts)
-        cst = [[sum.(SemiclassicalJacobiWeight.(t,1,1,0:ms[end])) for t in ts],
-                [sum.(SemiclassicalJacobiWeight.(t,2,0,0:ms[end])) for t in ts],
-                [sum.(SemiclassicalJacobiWeight.(t,0,2,0:ms[end])) for t in ts]]
+        # Use broadcast notation to compute all the derivative matrices across
+        # all the intervals and Fourier modes simultaneously
+        Z = ZernikeAnnulus{T}.(ρs,1,1)
+        Ds = (Z .\ (Laplacian.(axes.(Z,1)).*Weighted.(Z)))
+        Ds = [Ds[i].ops for i in 1:K+1-κ];
     else
-        Ls = _ann2element_via_lowering.(ts)
-        cst = [[sum.(SemiclassicalJacobiWeight.(t,a,a,0:ms[end])) for t in ts] for a in 1:-1:0]
+        ts = []
+        cst = [[]]
     end
 
-    # Use broadcast notation to compute all the derivative matrices across
-    # all the intervals and Fourier modes simultaneously
-    Z = ZernikeAnnulus{T}.(ρs,1,1)
-    Ds = (Z .\ (Laplacian.(axes.(Z,1)).*Weighted.(Z)))
-    Ds = [Ds[i].ops for i in 1:K+1-κ];
 
     # Loop over the Fourier modes
     Fs = []
@@ -128,7 +134,7 @@ function _getFs(N::Int, points::AbstractVector{T}, via_Jacobi::Bool) where T
 
         D = NTuple{K+1-κ, AbstractMatrix}([Ds[i][m+1] for i in 1:K+1-κ])
 
-        normalize_constants = [[cst[k][i][m+1] for k in 1:lastindex(cst)] for i in 1:K+1-κ]
+        normalize_constants = Vector{T}[T[cst[k][i][m+1] for k in 1:lastindex(cst)] for i in 1:K+1-κ]
 
         # Construct the structs for each Fourier mode seperately
         append!(Fs, [FiniteContinuousZernikeMode(M, points, m, j, L₁₁, L₀₁, L₁₀, D, normalize_constants, via_Jacobi, same_ρs, N)])
@@ -201,14 +207,14 @@ end
 #     zero_dirichlet_bcs!.(Fs, Δ)
 # end
 
-function zero_dirichlet_bcs!(F::FiniteContinuousZernike{T}, Δ::AbstractVector{<:AbstractMatrix{T}}) where T
+function zero_dirichlet_bcs!(F::FiniteContinuousZernike{T}, Δ::AbstractVector{<:AbstractMatrix}) where T
     @assert length(Δ) == 2*F.N-1
     # @assert Δ[1] isa LinearAlgebra.Symmetric{T, <:ArrowheadMatrix{T}}
     Fs = F.Fs #_getFs(F.N, F.points)
     zero_dirichlet_bcs!.(Fs, Δ)
 end
 
-function zero_dirichlet_bcs!(F::FiniteContinuousZernike{T}, Mf::AbstractVector{<:PseudoBlockVector{T}}) where T
+function zero_dirichlet_bcs!(F::FiniteContinuousZernike{T}, Mf::AbstractVector{<:PseudoBlockVector}) where T
     @assert length(Mf) == 2*F.N-1
     Fs = F.Fs #_getFs(F.N, F.points)
     zero_dirichlet_bcs!.(Fs, Mf)
@@ -233,10 +239,12 @@ function _bubble2disk_or_ann_all_modes(F::FiniteContinuousZernike, us::AbstractV
         m = ms[i]
         γs = _getγs(points, m)
         append!(γs, one(T))
-        if first(points) ≈ 0 && K > 1
+        if first(points) ≈ 0
             Us[1:Ms[i]-1,i,1] = [us[i][1]*γs[1];us[i][K+1:K:end]]
-            for k = 1:K-1 
-                Us[1:Ms[i],i,k+1] = [us[i][k];us[i][k+1]*γs[k+1];us[i][(K+k+1):K:end]] 
+            if K > 1
+                for k = 1:K-1
+                    Us[1:Ms[i],i,k+1] = [us[i][k];us[i][k+1]*γs[k+1];us[i][(K+k+1):K:end]] 
+                end
             end
         else
             for k = 1:K 
@@ -264,10 +272,11 @@ function _bubble2disk_or_ann_all_modes(F::FiniteContinuousZernike, us::AbstractV
 end
 
 
-function finite_plotvalues(F::FiniteContinuousZernike, us::AbstractVector; N=0)
+function finite_plotvalues(F::FiniteContinuousZernike, us::AbstractVector; N=0, K=0)
     T = eltype(F)
     _, Ũs = _bubble2disk_or_ann_all_modes(F, us)
-    points, K = T.(F.points), length(F.points)-1
+    points = T.(F.points)
+    K = K==0 ? length(F.points)-1 : K
     N = N == 0 ? F.N : N
     θs, rs, vals = [], [], []
     for k in 1:K
@@ -301,11 +310,11 @@ function _inf_error(K::Int, θs::AbstractVector, rs::AbstractVector, vals::Abstr
     for k = 1:K
         append!(vals_, [abs.(vals[k] - u.(RadialCoordinate.(rs[k],θs[k]')))])
     end
-    vals_, sum(maximum.(vals_))
+    vals_, norm((norm.(vals_, Inf)), Inf)
 end
 
-function inf_error(F::FiniteContinuousZernike{T}, θs::AbstractVector, rs::AbstractVector, vals::AbstractVector, u::Function) where T
-    K = lastindex(F.points)-1
+function inf_error(F::FiniteContinuousZernike{T}, θs::AbstractVector, rs::AbstractVector, vals::AbstractVector, u::Function;K=0) where T
+    K = K==0 ? lastindex(F.points)-1 : K
     _inf_error(K, θs, rs, vals, u)
 end
 

src/finitecontinuousmode.jl

@@ -86,6 +86,7 @@ end
 function _getγs(points::AbstractArray{T}, m::Int) where T
     K = length(points)-1
     first(points) > 0 && return [(one(T)-(points[k+1]/points[k+2])^2)*(points[k+1]/points[k+2])^m / (one(T)-(points[k]/points[k+1])^2) for k in 1:K-1]
+    K == 1 && return T[]
     γ = [(one(T)-(points[k+1]/points[k+2])^2)*(points[k+1]/points[k+2])^m / (one(T)-(points[k]/points[k+1])^2) for k in 2:K-1]
     return append!([(one(T)-(points[2]/points[3])^2)*(points[2]/points[3])^m / (sqrt(convert(T,2)^(m+2-iszero(m))/π) * normalizedjacobip(0, 0, m, 1.0))],γ)
 end
@@ -121,7 +122,11 @@ function ldiv(F::FiniteContinuousZernikeMode{T}, f::AbstractQuasiVector) where T
 
     bubbles = zeros(T, N-2, K)
     if first(points) ≈ 0
-        hats = vcat([fs[i][1] for i in 2:K-1], fs[end][1:2])
+        if K == 1
+            hats = [fs[1][1]]
+        else
+            hats = vcat([fs[i][1] for i in 2:K-1], fs[end][1:2])
+        end
         bubbles[:,1] = fs[1][2:N-1]
         for i in 2:K bubbles[:,i] = fs[i][3:N] end
     else
@@ -141,21 +146,25 @@ end
 function _build_top_left_block(F::FiniteContinuousZernikeMode{T}, Ms, γs::AbstractArray{T}, p::T) where T
     K = length(Ms)
     if p ≈ 0
-        a = [Ms[1][1,1]]
-        for i in 2:K
-            append!(a, [Ms[i][1,1]; Ms[i][2,2]])
-        end
+        if K > 1
+            a = [Ms[1][1,1]]
+            for i in 2:K
+                append!(a, [Ms[i][1,1]; Ms[i][2,2]])
+            end
 
-        dv = zeros(T, K)
-        dv[1] = a[1]*γs[1]^2 + a[2];
-        dv[end] = a[end];
+            dv = zeros(T, K)
+            dv[1] = a[1]*γs[1]^2 + a[2];
+            dv[end] = a[end];
 
-        for i in 1:K-2 dv[i+1] = a[2i+1]*γs[i+1]^2 + a[2i+2] end
+            for i in 1:K-2 dv[i+1] = a[2i+1]*γs[i+1]^2 + a[2i+2] end
 
-        ev = zeros(T, K-1)
-        γs = vcat(γs, one(T))
+            ev = zeros(T, K-1)
+            γs = vcat(γs, one(T))
 
-        for i in 2:K ev[i-1] = Ms[i][1,2] * γs[i] end
+            for i in 2:K ev[i-1] = Ms[i][1,2] * γs[i] end
+        else
+            dv = [Ms[1][1,1]]
+        end
     else
         a = []
         for i in 1:K
@@ -171,7 +180,9 @@ function _build_top_left_block(F::FiniteContinuousZernikeMode{T}, Ms, γs::Abstr
         for i in 1:K ev[i] = Ms[i][1,2] * γs[i] end
 
     end
-    Symmetric(BandedMatrix{T}(0=>dv, 1=>ev))
+
+    K == 1 && return Symmetric(BandedMatrix{T}(0=>dv))
+    return Symmetric(BandedMatrix{T}(0=>dv, 1=>ev))
 end
 
 # Interaction of the hats with the bubbles
@@ -185,12 +196,14 @@ function _build_second_block(F::FiniteContinuousZernikeMode{T}, Ms, γs::Abstrac
         if p ≈ 0
             append!(ev, [zeros(T, K-1)])
             dv[j][1] = Ms[1][1,j+1] * γs[1]
-            ev[j][1] = Ms[2][1,j+2]
-            for i in 2:K-1
-                dv[j][i] = Ms[i][2,j+2] * γs[i]
-                ev[j][i] = Ms[i+1][1,j+2]
+            if K > 1
+                ev[j][1] = Ms[2][1,j+2]
+                for i in 2:K-1
+                    dv[j][i] = Ms[i][2,j+2] * γs[i]
+                    ev[j][i] = Ms[i+1][1,j+2]
+                end
+                dv[j][K] = Ms[K][2,j+2]
             end
-            dv[j][K] = Ms[K][2,j+2]
         else
             append!(ev, [zeros(T, K)])
             for i in 1:K
@@ -200,6 +213,7 @@ function _build_second_block(F::FiniteContinuousZernikeMode{T}, Ms, γs::Abstrac
         end
     end
     if p ≈ 0
+        K == 1 && return [BandedMatrix{T}(0=>dv[j]) for j in 1:bs]
         return [BandedMatrix{T}(0=>dv[j], 1=>ev[j]) for j in 1:bs]
     else
         return [BandedMatrix{T}((0=>dv[j], -1=>ev[j]), (K+1, K)) for j in 1:bs]
@@ -211,7 +225,7 @@ function _build_trailing_bubbles(F::FiniteContinuousZernikeMode{T}, Ms, N::Int,
     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])
-        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])
+        Mn = vcat([Ms[1][2:N-1,2:N-1]], [reshape(view(Ms[i],3:N, 3:N)[:], N-2, N-2) for i in 2:K])
     else
         # 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]
@@ -259,6 +273,9 @@ end
     end
 end
 
+###
+# Assembly
+###
 @simplify function *(A::QuasiAdjoint{<:Any,<:FiniteContinuousZernikeMode}, B::BroadcastQuasiMatrix{<:Any, typeof(*), <:Tuple{BroadcastQuasiVector, FiniteContinuousZernikeMode}})
     λ, F = B.args
     T = promote_type(eltype(A), eltype(F))
@@ -339,7 +356,7 @@ function zero_dirichlet_bcs!(F::FiniteContinuousZernikeMode{T}, Δ::LinearAlgebr
     end
 end
 
-function zero_dirichlet_bcs!(F::FiniteContinuousZernikeMode{T}, A::Matrix{T}) where T
+function zero_dirichlet_bcs!(F::FiniteContinuousZernikeMode{T}, A::Matrix) where T
     points = F.points
     K = length(points)-1
     if first(points) > 0
@@ -352,7 +369,7 @@ function zero_dirichlet_bcs!(F::FiniteContinuousZernikeMode{T}, A::Matrix{T}) wh
     end
 end
 
-function zero_dirichlet_bcs!(F::FiniteContinuousZernikeMode{T}, Mf::PseudoBlockVector{T}) where T
+function zero_dirichlet_bcs!(F::FiniteContinuousZernikeMode{T}, Mf::PseudoBlockVector) where T
     points = F.points
     K = length(points)-1
     if !(first(points) ≈  0)

src/finitezernike.jl

@@ -313,9 +313,10 @@ function _bubble2disk_or_ann_all_modes(Z::FiniteZernikeBasis{T}, us::AbstractVec
 end
 
 
-function finite_plotvalues(Z::FiniteZernikeBasis{T}, us::AbstractVector; N=0) where T
+function finite_plotvalues(Z::FiniteZernikeBasis{T}, us::AbstractVector; N=0, K=0) where T
     _, Ũs = _bubble2disk_or_ann_all_modes(Z, us)
-    points, K = T.(Z.points), length(Z.points)-1
+    points= T.(Z.points)
+    K = K==0 ? lastindex(Z.points)-1 : K
     N = N == 0 ? Z.N : N
     Zs = Z.Zs
     θs=[]; rs=[]; vals = []   
@@ -341,8 +342,8 @@ function finite_plotvalues(Z::FiniteZernikeBasis{T}, us::AbstractVector; N=0) wh
 end
 
 ### Error collection
-function inf_error(Z::FiniteZernikeBasis{T}, θs::AbstractVector, rs::AbstractVector, vals::AbstractVector, u::Function) where T
-    K = lastindex(Z.points)-1
+function inf_error(Z::FiniteZernikeBasis{T}, θs::AbstractVector, rs::AbstractVector, vals::AbstractVector, u::Function; K=0) where T
+    K = K==0 ? lastindex(Z.points)-1 : K
     _inf_error(K, θs, rs, vals, u)
 end
 

src/plotting/plots.jl

@@ -4,13 +4,17 @@ using PyPlot, DelimitedFiles, LaTeXStrings
 ### Helper plot functions
 #######
 
-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
+function _plot(K::Int, θs::AbstractVector, rs::AbstractVector, vals::AbstractVector;ρ::T=0.0, ttl=[], vminmax=[]) where T
+    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))
+    if vminmax == []
+        vmin,vmax = minimum(minimum.(vals)), maximum(maximum.(vals))
+    else
+        vmin,vmax = vminmax[1], vminmax[2]
+    end
     norm = PyPlot.matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
 
     if ρ > 0.0
@@ -34,18 +38,18 @@ function _plot(K::Int, θs::AbstractVector, rs::AbstractVector, vals::AbstractVe
     display(gcf())
 end
 
-function plot(F::FiniteContinuousZernike{T}, θs::AbstractVector, rs::AbstractVector, vals::AbstractVector;ρ::T=0.0, ttl=[]) where T
-    K = lastindex(F.points)-1
-    _plot(K, θs, rs, vals, ρ=ρ, ttl=ttl)
+function plot(F::FiniteContinuousZernike{T}, θs::AbstractVector, rs::AbstractVector, vals::AbstractVector;ρ::T=0.0, ttl=[], vminmax=[],K=0) where T
+    K = K ==0 ? lastindex(Z.points)-1 : K
+    _plot(K, θs, rs, vals, ρ=ρ, ttl=ttl, vminmax=vminmax)
 end
 
-function plot(F::FiniteContinuousZernikeMode{T}, θs::AbstractVector, rs::AbstractVector, vals::AbstractVector;ρ::T=0.0, ttl=[]) where T
-    K = lastindex(F.points)-1
+function plot(F::FiniteContinuousZernikeMode{T}, θs::AbstractVector, rs::AbstractVector, vals::AbstractVector;ρ::T=0.0, ttl=[], K=0) where T
+    K = K ==0 ? lastindex(Z.points)-1 : K
     _plot(K, θs, rs, vals, ρ=ρ, ttl=ttl)
 end
 
-function plot(Z::FiniteZernikeBasis{T}, θs::AbstractVector, rs::AbstractVector, vals::AbstractVector;ρ::T=0.0, ttl=[]) where T
-    K = lastindex(Z.points)-1
+function plot(Z::FiniteZernikeBasis{T}, θs::AbstractVector, rs::AbstractVector, vals::AbstractVector;ρ::T=0.0, ttl=[], K=0) where T
+    K = K ==0 ? lastindex(Z.points)-1 : K
     _plot(K, θs, rs, vals, ρ=ρ, ttl=ttl)
 end
 

test/test_finiteannulus.jl

@@ -83,6 +83,27 @@ end
             # NIntegrate[Integrate[Exp[-10*(r^2*Cos[y]^2 + (r*Sin[y]-0.6)^2)]^2* (1-r^2)^2*(r^2-0.2^2)^2*r,{r,0,1}],{y,0,2 Pi}]
             @test sum(transpose.(fc) .* (M .* fc)) ≈  0.005578088780274445
             @test sum(transpose.(fc) .* (Δ .* fc)) ≈ 0.16877589535690113
+
+            # just disk element
+            N = 50; points = [0.0;1.0]
+            K = length(points)-1
+
+            F = FiniteContinuousZernike(N, points, via_Jacobi=via_Jacobi)
+            fc = F \ f0.(axes(F,1))
+            (θs, rs, vals) = finite_plotvalues(F, fc)
+            vals_, err = inf_error(F, θs, rs, vals, f0)
+            @test err < 1e-14
+            M = F' * F
+            @test length(M) == 2N-1
+            @test size(M[1]) == (K*(N/2-1), K*(N/2-1))
+            @test size(M[end]) == (2K, 2K)
+            ∇ = Derivative(axes(F,1)); Δ = (∇*F)' * (∇*F)
+            @test length(Δ) == 2N-1
+            @test size(Δ[1]) == (K*(N/2-1), K*(N/2-1))
+            @test size(Δ[end]) == (2K, 2K)
+            # NIntegrate[Integrate[Exp[-10*(r^2*Cos[y]^2 + (r*Sin[y]-0.6)^2)]^2* (1-r^2)^2*(r^2-0.2^2)^2*r,{r,0,1}],{y,0,2 Pi}]
+            @test sum(transpose.(fc) .* (M .* fc)) ≈  0.005578088780274445
+            @test sum(transpose.(fc) .* (Δ .* fc)) ≈ 0.16877589535690113
         end
     end
 
@@ -111,6 +132,16 @@ end
         @test size(M[end]) == (2K, 2K)
         # NIntegrate[Integrate[Sin[10*r*Cos[y]]^2*Sin[r^2]*r,{r,0.0,0.8}],{y,0,2 Pi}]
         @test sum(transpose.(fc) .* (M .* fc)) ≈  0.3055743848155512
+
+        # just disk element
+        F = FiniteContinuousZernike(N, [0.0;0.8])
+        fc = F \ plane_wave.(axes(F,1))
+        M = F' * (λ.(axes(F,1)) .* F)
+        @test length(M) == 2N-1
+        @test size(M[1]) == ((N/2-1), (N/2-1))
+        @test size(M[end]) == (2, 2)
+        # NIntegrate[Integrate[Sin[10*r*Cos[y]]^2*Sin[r^2]*r,{r,0.0,0.8}],{y,0,2 Pi}]
+        @test sum(transpose.(fc) .* (M .* fc)) ≈  0.3055743848155512
     end
 
 end