From 585f17377a726e04381b80cf839469cafc59c66d Mon Sep 17 00:00:00 2001 From: ioannisPApapadopoulos Date: Fri, 3 Nov 2023 13:15:33 +0000 Subject: [PATCH] one cell disk now working --- src/diskelement.jl | 2 +- src/finitecontinuous.jl | 79 +++++++++++++++++++++---------------- src/finitecontinuousmode.jl | 59 +++++++++++++++++---------- src/finitezernike.jl | 9 +++-- src/plotting/plots.jl | 28 +++++++------ test/test_finiteannulus.jl | 31 +++++++++++++++ 6 files changed, 135 insertions(+), 73 deletions(-) diff --git a/src/diskelement.jl b/src/diskelement.jl index 14f7c2a..05df7f5 100644 --- a/src/diskelement.jl +++ b/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 diff --git a/src/finitecontinuous.jl b/src/finitecontinuous.jl index f1a179f..b813662 100644 --- a/src/finitecontinuous.jl +++ b/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 diff --git a/src/finitecontinuousmode.jl b/src/finitecontinuousmode.jl index bade559..2d5b9d0 100644 --- a/src/finitecontinuousmode.jl +++ b/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) diff --git a/src/finitezernike.jl b/src/finitezernike.jl index 010a4aa..163e608 100644 --- a/src/finitezernike.jl +++ b/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 diff --git a/src/plotting/plots.jl b/src/plotting/plots.jl index 85628ab..44b4a41 100644 --- a/src/plotting/plots.jl +++ b/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 diff --git a/test/test_finiteannulus.jl b/test/test_finiteannulus.jl index ef6e21c..bbf817d 100644 --- a/test/test_finiteannulus.jl +++ b/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 \ No newline at end of file