From 01fae2fcbca3bd494e61796ed8c387ac68fe499d Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 18 Feb 2025 15:56:18 +0000 Subject: [PATCH 1/2] Vector spaces for Disks --- test/test_diskvector.jl | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) create mode 100644 test/test_diskvector.jl diff --git a/test/test_diskvector.jl b/test/test_diskvector.jl new file mode 100644 index 0000000..31f6003 --- /dev/null +++ b/test/test_diskvector.jl @@ -0,0 +1,25 @@ +using MultivariateOrthogonalPolynomials, Test, ForwardDiff +using ForwardDiff: gradient + + +k = 0; m = 0; n = 2 + +Z_x = 𝐱 -> gradient(𝐱 -> zernikez(n,m,𝐱), 𝐱)[1] +Z_y = 𝐱 -> gradient(𝐱 -> zernikez(n,m,𝐱), 𝐱)[2] + + +𝐱 = SVector(0.1,0.2) +r = norm(𝐱); ΞΈ = atan(𝐱[2], 𝐱[1]) +z = 2r^2 - 1 +zernikez(n,m,𝐱) + +W = (n,a,b) -> 2^(a+b+1)/(2n+a+b+1) * gamma(n+a+1)gamma(n+b+1)/(gamma(n+a+b+1)factorial(n)) + +@test jacobip(n,k, m,z) / sqrt(W(n,k,m)) β‰ˆ normalizedjacobip(n,k,m,z) + +r^m * cos(m*ΞΈ) * jacobip(n,k, m,z) / sqrt(W(n,k,m) / 2^(2+k+m)) + +sqrt(Ο€) * zernikez(n,m,𝐱) + +@time expand(Zernike(), Z_x) + From 962640743f9d88b79000da2d1194a8add2df2cb6 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Wed, 5 Mar 2025 15:19:45 -0600 Subject: [PATCH 2/2] Redefine zernikez --- src/MultivariateOrthogonalPolynomials.jl | 3 +- src/disk.jl | 10 ++++- test/test_diskvector.jl | 52 ++++++++++++++++++++++-- 3 files changed, 58 insertions(+), 7 deletions(-) diff --git a/src/MultivariateOrthogonalPolynomials.jl b/src/MultivariateOrthogonalPolynomials.jl index 5d1f42c..b8ffc22 100644 --- a/src/MultivariateOrthogonalPolynomials.jl +++ b/src/MultivariateOrthogonalPolynomials.jl @@ -4,7 +4,7 @@ using QuasiArrays: AbstractVector using ClassicalOrthogonalPolynomials, FastTransforms, BlockBandedMatrices, BlockArrays, DomainSets, QuasiArrays, StaticArrays, ContinuumArrays, InfiniteArrays, InfiniteLinearAlgebra, LazyArrays, SpecialFunctions, LinearAlgebra, BandedMatrices, LazyBandedMatrices, ArrayLayouts, - HarmonicOrthogonalPolynomials, RecurrenceRelationships + HarmonicOrthogonalPolynomials, RecurrenceRelationships, FillArrays import Base: axes, in, ==, +, -, /, *, ^, \, copy, copyto!, OneTo, getindex, size, oneto, all, resize!, BroadcastStyle, similar, fill!, setindex!, convert, show, summary, diff import Base.Broadcast: Broadcasted, broadcasted, DefaultArrayStyle @@ -14,6 +14,7 @@ import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle, domain import ContinuumArrays: @simplify, Weight, weight, grid, plotgrid, TransformFactorization, ExpansionLayout, plotvalues, unweighted, plan_transform, checkpoints, transform_ldiv, AbstractBasisLayout, basis_axes, Inclusion, grammatrix, weaklaplacian, layout_broadcasted, laplacian, abslaplacian, laplacian_axis, abslaplacian_axis, diff_layout, operatororder, broadcastbasis import ArrayLayouts: MemoryLayout, sublayout, sub_materialize +import FillArrays: SquareEye import BlockArrays: block, blockindex, BlockSlice, viewblock, blockcolsupport, AbstractBlockStyle, BlockStyle import BlockBandedMatrices: _BandedBlockBandedMatrix, AbstractBandedBlockBandedMatrix, _BandedMatrix, blockbandwidths, subblockbandwidths import LinearAlgebra: factorize diff --git a/src/disk.jl b/src/disk.jl index c7e4a0c..2622d3f 100644 --- a/src/disk.jl +++ b/src/disk.jl @@ -68,7 +68,7 @@ summary(io::IO, P::Zernike) = print(io, "Zernike($(P.a), $(P.b))") orthogonalityweight(Z::Zernike) = ZernikeWeight(Z.a, Z.b) -zerniker(β„“, m, a, b, r::T) where T = sqrt(convert(T,2)^(m+a+b+2-iszero(m))/Ο€) * r^m * normalizedjacobip((β„“-m) Γ· 2, b, m+a, 2r^2-1) +zerniker(β„“, m, a, b, r::T) where T = sqrt(convert(T,2)^(m+a+b+2-iszero(m))/Ο€) * r^m * normalizedjacobip(β„“, b, m+a, 2r^2-1) zerniker(β„“, m, b, r) = zerniker(β„“, m, zero(b), b, r) zerniker(β„“, m, r) = zerniker(β„“, m, zero(r), r) @@ -86,7 +86,7 @@ function getindex(Z::Zernike{T}, rΞΈ::RadialCoordinate, B::BlockIndex{1}) where β„“ = Int(block(B))-1 k = blockindex(B) m = iseven(β„“) ? k-isodd(k) : k-iseven(k) - zernikez(β„“, (isodd(k+β„“) ? 1 : -1) * m, Z.a, Z.b, rΞΈ) + zernikez((β„“-m) Γ· 2, (isodd(k+β„“) ? 1 : -1) * m, Z.a, Z.b, rΞΈ) end @@ -94,6 +94,7 @@ getindex(Z::Zernike, xy::StaticVector{2}, B::BlockIndex{1}) = Z[RadialCoordinate getindex(Z::Zernike, xy::StaticVector{2}, B::Block{1}) = [Z[xy, B[j]] for j=1:Int(B)] getindex(Z::Zernike, xy::StaticVector{2}, JR::BlockOneTo) = mortar([Z[xy,Block(J)] for J = 1:Int(JR[end])]) +basis_axes(::Inclusion{<:Any,<:UnitDisk}, v) = Zernike() ### # Jacobi matrices @@ -192,6 +193,11 @@ end # Transforms ### +function grammatrix(Z::Zernike{T}) where T + @assert Z.a == Z.b == 0 + SquareEye{T}((axes(Z,2),)) +end + function grid(S::Zernike{T}, B::Block{1}) where T N = Int(B) Γ· 2 + 1 #Β matrix rows M = 4N-3 # matrix columns diff --git a/test/test_diskvector.jl b/test/test_diskvector.jl index 31f6003..0778aef 100644 --- a/test/test_diskvector.jl +++ b/test/test_diskvector.jl @@ -1,11 +1,11 @@ -using MultivariateOrthogonalPolynomials, Test, ForwardDiff +using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, Test, ForwardDiff, StaticArrays using ForwardDiff: gradient k = 0; m = 0; n = 2 -Z_x = 𝐱 -> gradient(𝐱 -> zernikez(n,m,𝐱), 𝐱)[1] -Z_y = 𝐱 -> gradient(𝐱 -> zernikez(n,m,𝐱), 𝐱)[2] +Z_x = (n,m) -> (𝐱 -> gradient(𝐱 -> zernikez(n,m,𝐱), 𝐱)[1]) +Z_y = (n,m) -> (𝐱 -> gradient(𝐱 -> zernikez(n,m,𝐱), 𝐱)[2]) 𝐱 = SVector(0.1,0.2) @@ -21,5 +21,49 @@ r^m * cos(m*ΞΈ) * jacobip(n,k, m,z) / sqrt(W(n,k,m) / 2^(2+k+m)) sqrt(Ο€) * zernikez(n,m,𝐱) -@time expand(Zernike(), Z_x) +@time expand(Zernike(), Z_x(3,2)) +# vector OPs +o = expand(Zernike(), _ -> 1) + +v = [[o,0*o], [0*o,o]] +ip = (v,w) -> dot(v[1],w[1]) + dot(v[2],w[2]) +ip(v[1],v[2]) + +expand(Zernike(), Z_x(3,2)) + + +W_x = (n,m) -> (𝐱 -> gradient(𝐱 -> (1-norm(𝐱)^2)*zernikez(n,m,1,𝐱), 𝐱)[1]) +W_y = (n,m) -> (𝐱 -> gradient(𝐱 -> (1-norm(𝐱)^2)*zernikez(n,m,1,𝐱), 𝐱)[2]) + +βˆ‡W = (n,m) -> [expand(Zernike(), W_x(n,m)),expand(Zernike(), W_y(n,m))] + +ip(βˆ‡W(2,3), [expand(Zernike(), splat((x,y) -> 1+x+y+x^2+x*y+y^2+x^3+x^2*y+x*y^2+y^3)),expand(Zernike(), splat((x,y) -> 1+x+y+x^2+x*y+y^2+x^3+x^2*y+x*y^2+y^3))]) + +ip(,[expand(Zernike(), W_x(3,2)),expand(Zernike(), W_y(3,2))]) + +w = [expand(Zernike(), splat((x,y)->1-y^2)) expand(Zernike(), splat((x,y)->x*y)); expand(Zernike(), splat((x,y)->x*y)) expand(Zernike(), splat((x,y)->1-x^2))] + +wiW1 = (n,m) -> expand(Zernike()[:,Block.(1:20)], splat((x,y) -> [1-x^2,-x*y]' * gradient(𝐱 -> (1-norm(𝐱)^2)*zernikez(n,m,1,𝐱), SVector(x,y))/(1-x^2-y^2))) +wiW2 = (n,m) -> expand(Zernike()[:,Block.(1:20)], splat((x,y) -> [-x*y,1-y^2]' * gradient(𝐱 -> (1-norm(𝐱)^2)*zernikez(n,m,1,𝐱), SVector(x,y))/(1-x^2-y^2))) + +[wiW1(3,4),wiW2(3,4)], [expand(Zernike(), splat((x,y) -> 1+x+y+x^2+x*y+y^2+x^3+x^2*y+x*y^2+y^3)),expand(Zernike(), splat((x,y) -> 1+x+y+x^2+x*y+y^2+x^3+x^2*y+x*y^2+y^3))] + +v = [wiW1(3,4),wiW2(3,4)] +[ip(v, βˆ‡W(n,m)) for n=0:5, m=0:5] +dot(βˆ‡W(0,6)[1], βˆ‡W(4,6)[1]) +dot(βˆ‡W(0,6)[2], βˆ‡W(4,6)[2]) + + +βˆ‡W(0,6)[1][SVector(0.1,0.2)] +gradient(𝐱 -> (1-norm(𝐱)^2)*zernikez(0,6,1,𝐱), SVector(0.1,0.2)) +ip(βˆ‡W(8,4), βˆ‡W(9,4)) +[ip(βˆ‡W(8,4), βˆ‡W(n,m)) for n=0:10, m=0:6] +v = [wiW1(3,4),wiW2(3,4)] +[ip(v, βˆ‡W(n,m)) for n=0:10, m=0:6] + +zernikez(4,6,1,0.1*SVector(cos(0.2),sin(0.2))) +v[1][SVector(0.1,0.2)] + +(1-x^2) *P^(1,1) * (1-x^2) *P^(1,1) +(𝐱 -> (1-norm(𝐱)^2)*zernikez(3,4,1,𝐱))(0.1*SVector(cos(0.2),sin(0.2))) \ No newline at end of file