Jp/up

.github/workflows/ci.yml

@@ -23,12 +23,12 @@ jobs:
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
+      - uses: actions/cache@v4
         env:
           cache-name: cache-artifacts
         with:
@@ -41,6 +41,6 @@ jobs:
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v1
+      - uses: codecov/codecov-action@v3
         with:
           file: lcov.info

Project.toml

@@ -1,7 +1,7 @@
 name = "RadialPiecewisePolynomials"
 uuid = "7dab568b-3cf7-4f91-a977-b4631dfca174"
 authors = ["john.papad "]
-version = "0.1.4"
+version = "0.1.5"
 
 [deps]
 AnnuliOrthogonalPolynomials = "de1797fd-24c3-4035-91a2-b52aecdcfb01"
@@ -27,25 +27,25 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
-AnnuliOrthogonalPolynomials = "0.0.6"
-BandedMatrices = "0.17, 1"
-BlockArrays = "1.1"
+AnnuliOrthogonalPolynomials = "0.1.0"
+BandedMatrices = "1.9"
+BlockArrays = "1.5"
 BlockBandedMatrices = "0.13"
-ClassicalOrthogonalPolynomials = "0.13"
-ContinuumArrays = "0.18"
+ClassicalOrthogonalPolynomials = "0.15"
+ContinuumArrays = "0.19"
 DomainSets = "0.7"
-FastTransforms = "0.16"
+FastTransforms = "0.17"
 FillArrays = "1.5"
 HypergeometricFunctions = "0.3"
-LazyArrays = "2.1"
+LazyArrays = "2.6"
 MatrixFactorizations = "3"
 Memoization = "0.2"
-MultivariateOrthogonalPolynomials = "0.7"
-PiecewiseOrthogonalPolynomials = "0.5"
-QuasiArrays = "0.11"
-SemiclassicalOrthogonalPolynomials = "0.5"
+MultivariateOrthogonalPolynomials = "0.9"
+PiecewiseOrthogonalPolynomials = "0.5.4"
+QuasiArrays = "0.12"
+SemiclassicalOrthogonalPolynomials = "0.7"
 SpecialFunctions = "2"
-StaticArrays = "1.6"
+StaticArrays = "1.9"
 julia = "1.10"
 
 [extras]

src/annuluselement.jl

@@ -147,12 +147,12 @@ function ldiv(C::ContinuousZernikeAnnulusElementMode{T}, f::AbstractQuasiVector)
     # Truncate machine error tail
     Ñ = findall(x->abs(x) > 2*eps(T), c̃)
     c̃ = isempty(Ñ) ? Zeros{T}(3) : c̃[1:Ñ[end]+min(5, length(c̃)-Ñ[end])]
-    N = length(c̃) # degree
+    N = min(length(c̃), size(C.R,1)) # degree
 
     R̃ = view(C.R, 1:N, 1:N)
 
     # convert from ZernikeAnnulus(ρ,w_a,w_a) to hats + Bubble
-    dat = R̃[1:N,1:N] \ c̃
+    dat = R̃[1:N,1:N] \ c̃[1:N]
     cfs = T[]
     pad(append!(cfs, dat), axes(C,2))
 end
@@ -193,7 +193,7 @@ function mass_matrix(C::ContinuousZernikeAnnulusElementMode)
 
     m₀ = _mass_m₀(C, m, t)
     # TODO fix the excess zeros
-    return ApplyArray(*,Diagonal(Fill(β^2*m₀,∞)), ApplyArray(*, C.R', C.R))
+    return ApplyArray(*,Diagonal(Fill(β^2*m₀,size(C.R,1))), ApplyArray(*, C.R', C.R))
 end
 
 
@@ -210,7 +210,7 @@ end
     # 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.
+    # corresponding 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))
@@ -229,7 +229,8 @@ function assembly_matrix(C::ContinuousZernikeAnnulusElementMode, Λ::AbstractMat
     m₀ = _mass_m₀(C, m, t)
 
     # TODO fix the excess zeros
-    ApplyArray(*,Diagonal(Fill(β^2*m₀,∞)), ApplyArray(*, C.R', ApplyArray(*, Λ, C.R)))
+    sz = size(C.R)
+    ApplyArray(*,Diagonal(Fill(β^2*m₀,sz)), ApplyArray(*, C.R', ApplyArray(*, view(Λ,1:sz[1],1:sz[2]), C.R)))
 end
 
 
@@ -309,8 +310,8 @@ function stiffness_matrix(C::ContinuousZernikeAnnulusElementMode)
         C = [W010(m, ρ) W_100_010(m, ρ) 4*m₀*(m+1); W_100_010(m, ρ) W100(m,ρ) -4*m₀*(m+1)]
     end
 
-    Δ = [[C[1:2,3]'; Zeros{T}(∞,2)] Δ]
-    Vcat(Hcat(C, Zeros{T}(2,∞)), Δ)
+    Δ = [[C[1:2,3]'; Zeros{T}(size(Δ,2)-1,2)] Δ]
+    Vcat(Hcat(C, Zeros{T}(2,size(Δ,2)-3)), Δ)
 end
 
 ###

src/continuouszernike.jl

@@ -25,18 +25,18 @@ end
 
 # Matrices for lowering to ZernikeAnnulus(0,0) via
 # direct lowering. Less stable, but probably lower complexity.
-function _ann2element_via_raising(t::T) where T
+function _ann2element_via_raising_N(t::T, N::Int) where T
     # {T} did not change the speed.
-    Q₀₀ = SemiclassicalJacobi{T}.(t, 0, 0, 0:∞)
-    Q₀₁ = SemiclassicalJacobi{T}.(t, 0, 1, 0:∞)
-    Q₁₀ = SemiclassicalJacobi{T}.(t, 1, 0, 0:∞)
-    Q₁₁ = SemiclassicalJacobi{T}.(t, 1, 1, 0:∞)
+    Q₀₀ = SemiclassicalJacobi{T}.(t, 0, 0, 0:N)
+    Q₀₁ = SemiclassicalJacobi{T}.(t, 0, 1, 0:N)
+    Q₁₀ = SemiclassicalJacobi{T}.(t, 1, 0, 0:N)
+    Q₁₁ = SemiclassicalJacobi{T}.(t, 1, 1, 0:N)
 
     R₁₁ = (Weighted.(Q₀₀) .\ Weighted.(Q₁₁)) / t^2
     R₀₁ = BroadcastVector{AbstractVector}((Q, P) -> (Weighted(Q) \ Weighted(P))[1:2,1] / t, Q₀₀, Q₀₁)
     R₁₀ = BroadcastVector{AbstractVector}((Q, P) -> (Weighted(Q) \ Weighted(P))[1:2,1] / t, Q₀₀, Q₁₀)
 
-    BroadcastVector{AbstractMatrix}((R11, R01, R10)->Hcat(Vcat(R10, Zeros{T}(∞)), Vcat(R01, Zeros{T}(∞)), R11), R₁₁, R₀₁, R₁₀)
+    BroadcastVector{AbstractMatrix}((R11, R01, R10)->Hcat(Vcat(R10, Zeros{T}(size(R11,1)+1)), Vcat(R01, Zeros{T}(size(R11,2))), R11), R₁₁, R₀₁, R₁₀)
 end
 
 function _getMs_ms_js(N::Int)
@@ -78,7 +78,7 @@ function _getFs(N::Int, points::AbstractVector{T}) where T
         # intervals and Fourier modes simultaneously.
 
 
-        Rs = _ann2element_via_raising.(ts)
+        Rs = _ann2element_via_raising_N.(ts, N)
         cst = [[sum.(SemiclassicalJacobiWeight.(t,a,a,0:ms[end])) for t in ts] for a in 1:-1:0]
 
         # Use broadcast notation to compute all the derivative matrices across
@@ -104,13 +104,13 @@ function _getFs(N::Int, points::AbstractVector{T}) where T
     for (M, m, j) in zip(Ms, ms, js)
         # Extract the lowering and differentiation matrices associated
         # with each Fourier mode and store in the Tuples
-        R = NTuple{K+1-κ, AbstractMatrix}([Rs[i][m+1] for i in 1:K+1-κ])
-        D = NTuple{K+1-κ, AbstractMatrix}([Ds[i][m+1] for i in 1:K+1-κ])
+        R = NTuple{K+1-κ, BandedMatrix}([BandedMatrix(view(Rs[i][m+1],1:M, 1:M), (1,2)) for i in 1:K+1-κ])
+        D = NTuple{K+1-κ, Tridiagonal}([Tridiagonal(Ds[i][m+1][1:M, 1:M]) 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-κ]
         Cs = Tuple(_getCs(points, m, j, N, R, D, normalize_constants, same_ρs))
 
-        # Construct the structs for each Fourier mode seperately
+        # Construct the structs for each Fourier mode separately
         append!(Fs, [ContinuousZernikeMode(M, points, m, j, Cs, normalize_constants, same_ρs, N)])
     end
     return Fs

src/zernikebasis.jl

@@ -156,7 +156,7 @@ function gram_matrix(C::ContinuousZernikeAnnulusElementMode, Ψ::ZernikeBasisMod
 
     if a == 0 && b == 0
         m₀ = _mass_m₀(C,m,t)
-        ApplyArray(*,Diagonal(Fill(β^2*m₀,∞)), C.R')
+        ApplyArray(*,Diagonal(Fill(β^2*m₀,size(C.R))), C.R')
     else
         error("L²-inner product between ContinuousZernikeAnnulusElementMode and ZernikeBasisModeElement not implemented for parameters Ψ.a = $a and Ψ.b = $b")
     end
@@ -176,7 +176,7 @@ function gram_matrix(C::ContinuousZernikeElementMode, Ψ::ZernikeBasisModeElemen
 
     if a == 0 && b == 0
         R = Zernike(0) \ Weighted(Zernike(1))
-        Vcat(Hcat(β^2, Zeros{T}(1,∞)), β^2*R.ops[m+1]')
+        Vcat(Hcat(β^2, Zeros{T}(1,size(R.ops[m+1],2))), β^2*R.ops[m+1]')
     else
         error("L²-inner product between ContinuousZernikeElementMode and ZernikeBasisModeElement not implemented for parameters Ψ.a = $a and Ψ.b = $b")
     end