Skip to content

Commit e4044b3

Browse files
dlfivefiftydlfivefifty
committed
add coordinate, disk tests pass
1 parent 25ccdb4 commit e4044b3

File tree

9 files changed

+82
-89
lines changed

9 files changed

+82
-89
lines changed

‎README.md

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ julia> using MultivariateOrthogonalPolynomials, StaticArrays, LinearAlgebra
1212
julia> P = JacobiTriangle()
1313
JacobiTriangle(0, 0, 0)
1414

15-
julia> x,y = first.(axes(P,1)), last.(axes(P,1));
15+
julia> x,y = coordinates(P);
1616

1717
julia> u = P * (P \ (exp.(x) .* cos.(y))) # Expand in Triangle OPs
1818
JacobiTriangle(0, 0, 0) * [1.3365085377830084, 0.5687967596428205, -0.22812040274224554, 0.07733064070637755, 0.016169744493985644, -0.08714886622738759, 0.00338435674992512, 0.01220019521126353, -0.016867598915573725, 0.003930461395801074 … ]

‎examples/diskheat.jl

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -3,8 +3,7 @@ pyplot() # pyplot supports disks
33

44
Z = Zernike(1)
55
W = Weighted(Z)
6-
xy = axes(W,1)
7-
x,y = first.(xy),last.(xy)
6+
x,y = coordinates(W)
87
Δ = Z \ Laplacian(xy) * W
98
S = Z \ W
109

‎examples/diskhelmholtz.jl

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -14,9 +14,8 @@ pyplot()
1414

1515
Z = Zernike(1)
1616
W = Weighted(Z) # w*Z
17-
xy = axes(Z, 1);
18-
x, y = first.(xy), last.(xy);
19-
Δ = Z \ (Laplacian(xy) * W)
17+
x, y = coordinates(W)
18+
Δ = Z \ laplacian(W)
2019
S = Z \ W # identity
2120

2221

‎examples/multivariatelanczos.jl

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -3,8 +3,7 @@ using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, Test
33
P = Legendre()
44
P² = RectPolynomial(P,P)
55
p₀ = expand(P², 𝐱 -> 1)
6-
𝐱 = axes(P²,1)
7-
x,y = first.(𝐱),last.(𝐱)
6+
x,y = coordinates(P²)
87
w = P²/P²\ (x-y).^2
98

109
w .* p₀

‎src/MultivariateOrthogonalPolynomials.jl

Lines changed: 16 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@ import Base.Broadcast: Broadcasted, broadcasted, DefaultArrayStyle
1111
import DomainSets: boundary, EuclideanDomain
1212

1313
import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle, domain
14-
import ContinuumArrays: @simplify, Weight, weight, grid, plotgrid, TransformFactorization, ExpansionLayout, plotvalues, unweighted, plan_transform, checkpoints, transform_ldiv, AbstractBasisLayout, basis_axes, Inclusion, grammatrix, weaklaplacian, layout_broadcasted
14+
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
1515

1616
import ArrayLayouts: MemoryLayout, sublayout, sub_materialize
1717
import BlockArrays: block, blockindex, BlockSlice, viewblock, blockcolsupport, AbstractBlockStyle, BlockStyle
@@ -33,7 +33,21 @@ export MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial,
3333
Zernike, ZernikeWeight, zerniker, zernikez,
3434
AngularMomentum,
3535
RadialCoordinate, Weighted, Block, jacobimatrix, KronPolynomial, RectPolynomial,
36-
grammatrix, oneto
36+
grammatrix, oneto, coordinates, Laplacian, AbsLaplacian
37+
38+
39+
laplacian_axis(::Inclusion{<:SVector{2}}, A; dims...) = diff(A, (2,0); dims...) + diff(A, (0, 2); dims...)
40+
abslaplacian_axis(::Inclusion{<:SVector{2}}, A; dims...) = -(diff(A, (2,0); dims...) + diff(A, (0, 2); dims...))
41+
coordinates(P) = (first.(axes(P,1)), last.(axes(P,1)))
42+
43+
function diff_layout(::AbstractBasisLayout, a, (k,j)::NTuple{2,Int}; dims...)
44+
(k < 0 || j < 0) && throw(ArgumentError("order must be non-negative"))
45+
k == j == 0 && return a
46+
((k,j) == (1,0) || (k,j) == (0,1)) && return diff(a, Val((k,j)); dims...)
47+
k ≥ j && diff(diff(a, (1,0)), (k-1,j))
48+
diff(diff(a, (0,1)), (k,j-1))
49+
end
50+
3751

3852
include("ModalInterlace.jl")
3953
include("clenshawkron.jl")

‎test/test_disk.jl

Lines changed: 37 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -69,12 +69,11 @@ import ForwardDiff: hessian
6969
end
7070
end
7171

72-
@testset "Jacobi matrices" begin
72+
@testset "Jacobi matrices" begin
7373
# Setup
7474
α = 10 * rand()
7575
Z = Zernike(α)
76-
xy = axes(Z, 1)
77-
x, y = first.(xy), last.(xy)
76+
x, y = coordinates(Z)
7877
n = 150
7978

8079
# X tests
@@ -95,7 +94,7 @@ import ForwardDiff: hessian
9594
# Y tests
9695
JY = zeros(n,n)
9796
for j = 1:n
98-
JY[1:n,j] = (Z \ (y .* Z[:,j]))[1:n]
97+
JY[1:n,j] = (Z \ (y .* Z[:,j]))[1:n]
9998
end
10099
Y = Z \ (y .* Z)
101100
# The Zernike Jacobi matrices are symmetric for this normalization
@@ -112,8 +111,8 @@ import ForwardDiff: hessian
112111
@test (X*Y)[Block.(1:6),Block.(1:6)] ≈ (X[Block.(1:10),Block.(1:10)]*Y[Block.(1:10),Block.(1:10)])[Block.(1:6),Block.(1:6)]
113112

114113
# Addition of Jacobi matrices
115-
@test (X+Y)[Block.(1:6),Block.(1:6)] ≈ (X[Block.(1:6),Block.(1:6)]+Y[Block.(1:6),Block.(1:6)])
116-
@test (Y+Y)[Block.(1:6),Block.(1:6)] ≈ (Y[Block.(1:6),Block.(1:6)]+Y[Block.(1:6),Block.(1:6)])
114+
@test (X+Y)[Block.(1:6),Block.(1:6)] ≈ X[Block.(1:6),Block.(1:6)]+Y[Block.(1:6),Block.(1:6)]
115+
@test (Y+Y)[Block.(1:6),Block.(1:6)] ≈ Y[Block.(1:6),Block.(1:6)]+Y[Block.(1:6),Block.(1:6)]
117116

118117
# for now, reject non-zero first parameter options
119118
@test_throws ErrorException("Implement for non-zero first basis parameter.") jacobimatrix(Val(1),Zernike(1,1))
@@ -128,7 +127,7 @@ import ForwardDiff: hessian
128127
end
129128

130129
Z = Zernike(a,b);
131-
xy = axes(Z,1); x,y = first.(xy),last.(xy);
130+
x,y = coordinates(Z)
132131
u = Z * (Z \ exp.(x .* cos.(y)))
133132
@test u[SVector(0.1,0.2)] ≈ exp(0.1cos(0.2))
134133
end
@@ -172,7 +171,7 @@ import ForwardDiff: hessian
172171
ur = r -> r^m*f(r^2)
173172
@test ur(r) ≈ (1-r^2) * zerniker(ℓ, m, b, r)
174173

175-
D = Derivative(axes(Chebyshev(),1))
174+
D = Derivative(Chebyshev())
176175
D1 = Normalized(Jacobi(0, m+1)) \ (D * (HalfWeighted{:a}(Normalized(Jacobi(1, m)))))
177176
D2 = HalfWeighted{:b}(Normalized(Jacobi(1, m))) \ (D * (HalfWeighted{:b}(Normalized(Jacobi(0, m+1)))))
178177

@@ -219,12 +218,11 @@ import ForwardDiff: hessian
219218
# @test lap(u, xy...) ≈ Zernike(1)[xy,15] * (-4) * 5 * 1
220219

221220
WZ = Weighted(Zernike(1)) # Zernike(1) weighted by (1-r^2)
222-
Δ = Laplacian(axes(WZ,1))
221+
Δ = Laplacian(WZ)
223222
Δ_Z = Zernike(1) \ (Δ * WZ)
224223
@test exp.(Δ_Z)[1:10,1:10] == exp.(Δ_Z[1:10,1:10])
225224

226-
xy = axes(WZ,1)
227-
x,y = first.(xy),last.(xy)
225+
x,y = coordinates(WZ)
228226
u = @. (1 - x^2 - y^2) * exp(x*cos(y))
229227
Δu = @. (-exp(x*cos(y)) * (4 - x*(-5 + x^2 + y^2)cos(y) + (-1 + x^2 + y^2)cos(y)^2 - 4x*y*sin(y) + x^2*(x^2 + y^2-1)*sin(y)^2))
230228
@test (WZ * (WZ \ u))[SVector(0.1,0.2)] ≈ u[SVector(0.1,0.2)]
@@ -233,12 +231,12 @@ import ForwardDiff: hessian
233231
@testset "Unweighted" begin
234232
c = [randn(100); zeros(∞)]
235233
Z = Zernike()
236-
Δ = Zernike(2) \ (Laplacian(axes(Z,1)) * Z)
234+
Δ = Zernike(2) \ (Laplacian(Z) * Z)
237235
@test tr(hessian(xy -> (Zernike{eltype(xy)}()*c)[xy], SVector(0.1,0.2))) ≈ (Zernike(2)*(Δ*c))[SVector(0.1,0.2)]
238236

239237
b = 0.2
240238
Z = Zernike(b)
241-
Δ = Zernike(b+2) \ (Laplacian(axes(Z,1)) * Z)
239+
Δ = Zernike(b+2) \ (Laplacian(Z) * Z)
242240
@test tr(hessian(xy -> (Zernike{eltype(xy)}(b)*c)[xy], SVector(0.1,0.2))) ≈ (Zernike(b+2)*(Δ*c))[SVector(0.1,0.2)]
243241
end
244242
end
@@ -325,9 +323,9 @@ end
325323
@testset "Fractional Laplacian on Unit Disk" begin
326324
@testset "Fractional Laplacian on Disk: (-Δ)^(β) == -Δ when β=1" begin
327325
WZ = Weighted(Zernike(1.))
328-
Δ = Laplacian(axes(WZ,1))
326+
Δ = Laplacian(WZ)
329327
Δ_Z = Zernike(1) \ (Δ * WZ)
330-
Δfrac = AbsLaplacian(axes(WZ,1),1.)
328+
Δfrac = AbsLaplacian(WZ,1.)
331329
Δ_Zfrac = Zernike(1) \ (Δfrac * WZ)
332330
@test Δ_Z[1:100,1:100] ≈ -Δ_Zfrac[1:100,1:100]
333331
end
@@ -338,10 +336,9 @@ end
338336
β = 1.34
339337
Z = Zernike(β)
340338
WZ = Weighted(Z)
341-
xy = axes(WZ,1)
342-
x,y = first.(xy),last.(xy)
339+
x,y = coordinates(WZ)
343340
# generate fractional Laplacian
344-
Δfrac = AbsLaplacian(axes(WZ,1),β)
341+
Δfrac = AbsLaplacian(WZ,β)
345342
Δ_Zfrac = Z \ (Δfrac * WZ)
346343
# define function whose fractional Laplacian is known
347344
u = @. (1 - x^2 - y^2).^β
@@ -355,10 +352,9 @@ end
355352
β = 2.11
356353
Z = Zernike(β)
357354
WZ = Weighted(Z)
358-
xy = axes(WZ,1)
359-
x,y = first.(xy),last.(xy)
355+
x,y = coordinates(WZ)
360356
# generate fractional Laplacian
361-
Δfrac = AbsLaplacian(axes(WZ,1),β)
357+
Δfrac = AbsLaplacian(WZ,β)
362358
Δ_Zfrac = Z \ (Δfrac * WZ)
363359
# define function whose fractional Laplacian is known
364360
u = @. (1 - x^2 - y^2).^β
@@ -371,10 +367,9 @@ end
371367
β = 3.14159
372368
Z = Zernike(β)
373369
WZ = Weighted(Z)
374-
xy = axes(WZ,1)
375-
x,y = first.(xy),last.(xy)
370+
x,y = coordinates(WZ)
376371
# generate fractional Laplacian
377-
Δfrac = AbsLaplacian(axes(WZ,1),β)
372+
Δfrac = AbsLaplacian(WZ,β)
378373
Δ_Zfrac = Z \ (Δfrac * WZ)
379374
# define function whose fractional Laplacian is known
380375
u = @. (1 - x^2 - y^2).^β
@@ -387,10 +382,9 @@ end
387382
β = 1.1
388383
Z = Zernike(β)
389384
WZ = Weighted(Z)
390-
xy = axes(WZ,1)
391-
x,y = first.(xy),last.(xy)
385+
x,y = coordinates(WZ)
392386
# generate fractional Laplacian
393-
Δfrac = AbsLaplacian(axes(WZ,1),β)
387+
Δfrac = AbsLaplacian(WZ,β)
394388
Δ_Zfrac = Z \ (Δfrac * WZ)
395389
# define function whose fractional Laplacian is known
396390
u = @. (1 - x^2 - y^2).^(β+1)
@@ -406,10 +400,9 @@ end
406400
β = 2.71999
407401
Z = Zernike(β)
408402
WZ = Weighted(Z)
409-
xy = axes(WZ,1)
410-
x,y = first.(xy),last.(xy)
403+
x,y = coordinates(WZ)
411404
# generate fractional Laplacian
412-
Δfrac = AbsLaplacian(axes(WZ,1),β)
405+
Δfrac = AbsLaplacian(WZ,β)
413406
Δ_Zfrac = Z \ (Δfrac * WZ)
414407
# define function whose fractional Laplacian is known
415408
u = @. (1 - x^2 - y^2).^(β+1)
@@ -425,10 +418,9 @@ end
425418
β = 2.71999
426419
Z = Zernike(β)
427420
WZ = Weighted(Z)
428-
xy = axes(WZ,1)
429-
x,y = first.(xy),last.(xy)
421+
x,y = coordinates(WZ)
430422
# generate fractional Laplacian
431-
Δfrac = AbsLaplacian(axes(WZ,1),β)
423+
Δfrac = AbsLaplacian(WZ,β)
432424
Δ_Zfrac = Z \ (Δfrac * WZ)
433425
# define function whose fractional Laplacian is known
434426
u = @. (1 - x^2 - y^2).^(β)*x
@@ -444,10 +436,9 @@ end
444436
β = 1.91239
445437
Z = Zernike(β)
446438
WZ = Weighted(Z)
447-
xy = axes(WZ,1)
448-
x,y = first.(xy),last.(xy)
439+
x,y = coordinates(WZ)
449440
# generate fractional Laplacian
450-
Δfrac = AbsLaplacian(axes(WZ,1),β)
441+
Δfrac = AbsLaplacian(WZ,β)
451442
Δ_Zfrac = Z \ (Δfrac * WZ)
452443
# define function whose fractional Laplacian is known
453444
u = @. (1 - x^2 - y^2).^(β)*y
@@ -464,10 +455,9 @@ end
464455
β = 1.21999
465456
Z = Zernike(β)
466457
WZ = Weighted(Z)
467-
xy = axes(WZ,1)
468-
x,y = first.(xy),last.(xy)
458+
x,y = coordinates(WZ)
469459
# generate fractional Laplacian
470-
Δfrac = AbsLaplacian(axes(WZ,1),β)
460+
Δfrac = AbsLaplacian(WZ,β)
471461
Δ_Zfrac = Z \ (Δfrac * WZ)
472462
# define function whose fractional Laplacian is known
473463
u = @. (1 - x^2 - y^2).^(β+1)*x
@@ -483,10 +473,9 @@ end
483473
β = 0.141
484474
Z = Zernike(β)
485475
WZ = Weighted(Z)
486-
xy = axes(WZ,1)
487-
x,y = first.(xy),last.(xy)
476+
x,y = coordinates(WZ)
488477
# generate fractional Laplacian
489-
Δfrac = AbsLaplacian(axes(WZ,1),β)
478+
Δfrac = AbsLaplacian(WZ,β)
490479
Δ_Zfrac = Z \ (Δfrac * WZ)
491480
# define function whose fractional Laplacian is known
492481
u = @. (1 - x^2 - y^2).^(β+1)*y
@@ -506,9 +495,9 @@ end
506495
Z = Zernike(β)
507496
WZ = Weighted(Z)
508497
xy = axes(WZ,1)
509-
x,y = first.(xy),last.(xy)
498+
x,y = coordinates(WZ)
510499
# generate fractional Laplacian
511-
Δfrac = AbsLaplacian(axes(WZ,1),β)
500+
Δfrac = AbsLaplacian(WZ,β)
512501
Δ_Zfrac = Z \ (Δfrac * WZ)
513502
# define function whose fractional Laplacian is known
514503
uexplicit = @. (1 - x^2 - y^2).^(β+1)
@@ -526,9 +515,9 @@ end
526515
Z = Zernike(β)
527516
WZ = Weighted(Z)
528517
xy = axes(WZ,1)
529-
x,y = first.(xy),last.(xy)
518+
x,y = coordinates(WZ)
530519
# generate fractional Laplacian
531-
Δfrac = AbsLaplacian(axes(WZ,1),β)
520+
Δfrac = AbsLaplacian(WZ,β)
532521
Δ_Zfrac = Z \ (Δfrac * WZ)
533522
# define function whose fractional Laplacian is known
534523
uexplicit = @. (1 - x^2 - y^2).^(β+1)*y
@@ -545,9 +534,9 @@ end
545534
Z = Zernike(β)
546535
WZ = Weighted(Z)
547536
xy = axes(WZ,1)
548-
x,y = first.(xy),last.(xy)
537+
x,y = coordinates(WZ)
549538
# generate fractional Laplacian
550-
Δfrac = AbsLaplacian(axes(WZ,1),β)
539+
Δfrac = AbsLaplacian(WZ,β)
551540
Δ_Zfrac = Z \ (Δfrac * WZ)
552541
# define function whose fractional Laplacian is known
553542
uexplicit = @. (1 - x^2 - y^2).^(β+1)*x

‎test/test_rect.jl

Lines changed: 3 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -26,8 +26,7 @@ using Base: oneto
2626
T,U = ChebyshevT(),ChebyshevU()
2727
T² = RectPolynomial(Fill(T, 2))
2828
T²ₙ = T²[:,Block.(Base.OneTo(5))]
29-
𝐱 = axes(T²ₙ,1)
30-
x,y = first.(𝐱),last.(𝐱)
29+
x,y = coordinates(T²ₙ)
3130
@test T²ₙ \ one.(x) == [1; zeros(14)]
3231
@test (T² \ x)[1:5] ≈[0;1;zeros(3)]
3332

@@ -50,8 +49,7 @@ using Base: oneto
5049
TU = RectPolynomial(T, U)
5150
X = jacobimatrix(Val{1}(), TU)
5251
Y = jacobimatrix(Val{2}(), TU)
53-
𝐱 = axes(TU, 1)
54-
x, y = first.(𝐱), last.(𝐱)
52+
x, y = coordinates(TU)
5553
N = 10
5654
KR = Block.(1:N)
5755
@test (TU \ (x .* TU))[KR,KR] == X[KR,KR]
@@ -159,8 +157,7 @@ using Base: oneto
159157
@testset "variable coefficients" begin
160158
T,U = ChebyshevT(), ChebyshevU()
161159
P = RectPolynomial(T, U)
162-
𝐱 = axes(P,1)
163-
x,y = first.(𝐱), last.(𝐱)
160+
x,y = coordinates(P)
164161
X = P\(x .* P)
165162
Y = P\(y .* P)
166163

‎test/test_rectdisk.jl

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,7 @@ import MultivariateOrthogonalPolynomials: dunklxu_raising, dunklxu_lowering, Ang
99
@test P ≠ DunklXuDisk(0.123)
1010

1111
xy = axes(P,1)
12-
x,y = first.(xy),last.(xy)
12+
x,y = coordinates(P)
1313
@test xy[SVector(0.1,0.2)] == SVector(0.1,0.2)
1414
@test x[SVector(0.1,0.2)] == 0.1
1515
@test y[SVector(0.1,0.2)] == 0.2
@@ -27,7 +27,7 @@ import MultivariateOrthogonalPolynomials: dunklxu_raising, dunklxu_lowering, Ang
2727
@test WP ≠ WQ
2828
@test WP == WP
2929

30-
x, y = first.(axes(P, 1)), last.(axes(P, 1))
30+
x, y = coordinates(P)
3131

3232
L = WP \ WQ
3333
R = Q \ P
@@ -39,8 +39,8 @@ import MultivariateOrthogonalPolynomials: dunklxu_raising, dunklxu_lowering, Ang
3939

4040
@test (DunklXuDisk() \ WeightedDunklXuDisk(1.0))[Block.(1:N), Block.(1:N)] ≈ (WeightedDunklXuDisk(0.0) \ WeightedDunklXuDisk(1.0))[Block.(1:N), Block.(1:N)]
4141

42-
∂x = Derivative(axes(P, 1), (1,0))
43-
∂y = Derivative(axes(P, 1), (0,1))
42+
∂x = Derivative(P, (1,0))
43+
∂y = Derivative(P, (0,1))
4444

4545
Dx = Q \ (∂x * P)
4646
Dy = Q \ (∂y * P)

0 commit comments

Comments
 (0)