diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 40a7191..688e9fd 100644 --- a/.github/workflows/ci.yml +++ b/.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 \ No newline at end of file diff --git a/Project.toml b/Project.toml index 9ae6d07..d91dcac 100644 --- a/Project.toml +++ b/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] diff --git a/src/annuluselement.jl b/src/annuluselement.jl index 6a54aee..6ad635c 100644 --- a/src/annuluselement.jl +++ b/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 ### diff --git a/src/continuouszernike.jl b/src/continuouszernike.jl index a29bcce..ae88110 100644 --- a/src/continuouszernike.jl +++ b/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 diff --git a/src/zernikebasis.jl b/src/zernikebasis.jl index fc46b4d..1a4ee17 100644 --- a/src/zernikebasis.jl +++ b/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