Skip to content

Commit efd46fe

Browse files
dlfivefiftydlfivefifty
committed
Disk transform
1 parent ef8b46f commit efd46fe

File tree

3 files changed

+151
-15
lines changed

3 files changed

+151
-15
lines changed

‎src/Disk/Disk.jl

Lines changed: 66 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,30 @@
1-
export ChebyshevDisk
1+
export ChebyshevDisk, ZernikeDisk
2+
3+
struct ZernikeDisk{V,T} <: Space{Disk{V},T}
4+
domain::Disk{V}
5+
end
6+
7+
ZernikeDisk(d::Disk{V}) where V = ZernikeDisk{V,complex(eltype(V))}(d)
8+
ZernikeDisk() = ZernikeDisk(Disk())
9+
10+
@containsconstants ZernikeDisk
11+
12+
points(K::ZernikeDisk, n::Integer) =
13+
fromcanonical.(Ref(K), points(ChebyshevDisk(), n))
14+
215

316
struct ChebyshevDisk{V,T} <: Space{Disk{V},T}
417
domain::Disk{V}
518
end
619

7-
ChebyshevDisk(d::Disk{V}) where V = ChebyshevDisk{V,complex(eltype(V))}(d)
20+
@containsconstants ChebyshevDisk
21+
22+
ChebyshevDisk(d::Disk{V}) where V = ChebyshevDisk{V,eltype(V)}(d)
823
ChebyshevDisk() = ChebyshevDisk(Disk())
924

10-
rectspace(_) = Chebyshev() * Laurent()
25+
rectspace(_) = Chebyshev() * Fourier()
26+
27+
spacescompatible(K1::ChebyshevDisk, K2::ChebyshevDisk) = domainscompatible(K1, K2)
1128

1229
function points(S::ChebyshevDisk, N)
1330
pts = points(rectspace(S), N)
@@ -66,3 +83,49 @@ fromrectcfs(S, v) = fromtensor(S, checkerboard(totensor(rectspace(S),v)))
6683

6784
evaluate(cfs::AbstractVector, S::ChebyshevDisk, x) = evaluate(torectcfs(S,cfs), rectspace(S), ipolar(x))
6885

86+
87+
88+
89+
struct FastZernikeDiskTransformPlan{DUF,CHEB}
90+
cxfplan::DUF
91+
disk2cxf::CHEB
92+
end
93+
94+
function FastZernikeDiskTransformPlan(S::ZernikeDisk, v::AbstractVector{T}) where T
95+
# n = floor(Integer,sqrt(2length(v)) + 1/2)
96+
# v = Array{T}(undef, sum(1:n))
97+
FastZernikeDiskTransformPlan(plan_transform(ChebyshevDisk(),v), CDisk2CxfPlan(n))
98+
end
99+
100+
function *(P::FastZernikeDiskTransformPlan, v)
101+
v̂ = P.cxfplan*v
102+
F = tridevec_trans(v̂)
103+
F̌ = P.disk2cheb \ F
104+
trivec(tridenormalize!(F̌))
105+
end
106+
107+
plan_transform(K::ZernikeDisk, v::AbstractVector) = FastZernikeDiskTransformPlan(K, v)
108+
109+
struct FastZernikeDiskITransformPlan{DUF,CHEB}
110+
icxfplan::DUF
111+
disk2cxf::CHEB
112+
end
113+
114+
function FastZernikeDiskITransformPlan(S::ZernikeDisk, v::AbstractVector{T}) where T
115+
# n = floor(Integer,sqrt(2length(v)) + 1/2)
116+
# v = Array{T}(undef, sum(1:n))
117+
FastZernikeDiskITransformPlan(plan_itransform(ChebyshevDisk(), v), CDisk2CxfPlan(n))
118+
end
119+
120+
function *(P::FastZernikeDiskITransformPlan, v)
121+
# n = floor(Integer,sqrt(2length(v)) + 1/2)
122+
# v = pad(v, sum(1:n))
123+
F̌ = trinormalize!(tridevec(v))
124+
F = P.disk2cxf * F̌
125+
v̂ = trivec(transpose(F))
126+
P.icxfplan*v̂
127+
end
128+
129+
130+
plan_itransform(K::ZernikeDisk, v::AbstractVector) = FastZernikeDiskITransformPlan(K, v)
131+
Base.sum(f::Fun{<:ZernikeDisk}) = Fun(f,ZernikeDisk()).coefficients[1]/π

‎src/c_transforms.jl

Lines changed: 38 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -6,21 +6,28 @@ const PlanPtr = Ptr{ft_plan_struct}
66
function c_plan_tri2cheb end
77
function c_tri2cheb end
88
function c_cheb2tri end
9+
function c_plan_disk2cxf end
10+
function c_disk2cxf end
11+
function c_cxf2disk end
912

1013
if Libdl.find_library(libfasttransforms) ≡ libfasttransforms
1114
ft_set_threads(n::Int) = ccall((:omp_set_num_threads, libfasttransforms), Nothing, (Int, ), n)
1215
ft_set_threads(Sys.CPU_THREADS)
1316

14-
c_plan_sph2fourier(n::Int) = ccall((:ft_plan_sph2fourier, libfasttransforms), PlanPtr, (Int64, ), n)
15-
fc_sph2fourier(P::PlanPtr, A::Matrix{Float64}) = ccall((:ft_execute_sph2fourier, libfasttransforms), Nothing, (PlanPtr, Ptr{Float64}, Int64, Int64), P, A, size(A, 1), size(A, 2))
17+
c_plan_sph2fourier(n::Int) = ccall((:ft_plan_sph2fourier, libfasttransforms), PlanPtr, (Int, ), n)
18+
fc_sph2fourier(P::PlanPtr, A::Matrix{Float64}) = ccall((:ft_execute_sph2fourier, libfasttransforms), Nothing, (PlanPtr, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
1619

17-
c_plan_rottriangle(n::Int, α::Float64, β::Float64, γ::Float64) = ccall((:ft_plan_rottriangle, libfasttransforms), PlanPtr, (Int64, Float64, Float64, Float64), n, α, β, γ)
18-
c_execute_tri_hi2lo(P::PlanPtr, A::Matrix{Float64}) = ccall((:ft_execute_tri_hi2lo, libfasttransforms), Nothing, (PlanPtr, Ptr{Float64}, Int64), P, A, size(A, 2))
19-
c_execute_tri_lo2hi(P::PlanPtr, A::Matrix{Float64}) = ccall((:ft_execute_tri_lo2hi, libfasttransforms), Nothing, (PlanPtr, Ptr{Float64}, Int64), P, A, size(A, 2))
20+
c_plan_rottriangle(n::Int, α::Float64, β::Float64, γ::Float64) = ccall((:ft_plan_rottriangle, libfasttransforms), PlanPtr, (Int, Float64, Float64, Float64), n, α, β, γ)
21+
c_execute_tri_hi2lo(P::PlanPtr, A::Matrix{Float64}) = ccall((:ft_execute_tri_hi2lo, libfasttransforms), Nothing, (PlanPtr, Ptr{Float64}, Int), P, A, size(A, 2))
22+
c_execute_tri_lo2hi(P::PlanPtr, A::Matrix{Float64}) = ccall((:ft_execute_tri_lo2hi, libfasttransforms), Nothing, (PlanPtr, Ptr{Float64}, Int), P, A, size(A, 2))
2023

21-
c_plan_tri2cheb(n::Int, α::Float64, β::Float64, γ::Float64) = ccall((:ft_plan_tri2cheb, libfasttransforms), PlanPtr, (Int64, Float64, Float64, Float64), n, α, β, γ)
22-
c_tri2cheb(P::PlanPtr, A::Matrix{Float64}) = ccall((:ft_execute_tri2cheb, libfasttransforms), Nothing, (PlanPtr, Ptr{Float64}, Int64, Int64), P, A, size(A, 1), size(A, 2))
23-
c_cheb2tri(P::PlanPtr, A::Matrix{Float64}) = ccall((:ft_execute_cheb2tri, libfasttransforms), Nothing, (PlanPtr, Ptr{Float64}, Int64, Int64), P, A, size(A, 1), size(A, 2))
24+
c_plan_tri2cheb(n::Int, α::Float64, β::Float64, γ::Float64) = ccall((:ft_plan_tri2cheb, libfasttransforms), PlanPtr, (Int, Float64, Float64, Float64), n, α, β, γ)
25+
c_tri2cheb(P::PlanPtr, A::Matrix{Float64}) = ccall((:ft_execute_tri2cheb, libfasttransforms), Nothing, (PlanPtr, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
26+
c_cheb2tri(P::PlanPtr, A::Matrix{Float64}) = ccall((:ft_execute_cheb2tri, libfasttransforms), Nothing, (PlanPtr, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
27+
28+
c_plan_disk2cxf(n::Int) = ccall((:ft_plan_disk2cxf, libfasttransforms), PlanPtr, (Int,), n)
29+
c_disk2cxf(P::PlanPtr, A::Matrix{Float64}) = ccall((:ft_execute_disk2cxf, libfasttransforms), Nothing, (PlanPtr, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
30+
c_cxf2disk(P::PlanPtr, A::Matrix{Float64}) = ccall((:ft_execute_cxf2disk, libfasttransforms), Nothing, (PlanPtr, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
2431
else
2532
@warn "Cannot load FastTransforms Dylib"
2633
end
@@ -33,9 +40,8 @@ struct CTri2ChebPlan
3340
γ::Float64
3441
end
3542

36-
function CTri2ChebPlan(n::Int, α::Float64, β::Float64, γ::Float64)
43+
CTri2ChebPlan(n::Int, α::Float64, β::Float64, γ::Float64) =
3744
CTri2ChebPlan(c_plan_tri2cheb(n, α, β, γ), n, α, β, γ)
38-
end
3945

4046
function *(C::CTri2ChebPlan, A::Matrix{Float64})
4147
size(A,1) == size(A,2) == C.n || throw(ArgumentError(A))
@@ -50,3 +56,25 @@ function \(C::CTri2ChebPlan, A::Matrix{Float64})
5056
c_cheb2tri(C.plan, B)
5157
B
5258
end
59+
60+
61+
struct CDisk2CxfPlan
62+
plan::PlanPtr
63+
n::Int
64+
end
65+
66+
CDisk2CxfPlan(n::Int) = CDisk2CxfPlan(c_plan_disk2cxf(n), n)
67+
68+
function *(C::CDisk2CxfPlan, A::Matrix{Float64})
69+
(size(A,1) == C.n+1 && size(A,2) == 4C.n+1) || throw(ArgumentError(A))
70+
B = copy(A)
71+
c_disk2cxf(C.plan, B)
72+
B
73+
end
74+
75+
function \(C::CDisk2CxfPlan, A::Matrix{Float64})
76+
(size(A,1) == C.n && size(A,2) == 4C.n-3) || throw(DimensionMismatch(""))
77+
B = copy(A)
78+
c_cxf2disk(C.plan, B)
79+
B
80+
end

‎test/test_disk.jl

Lines changed: 47 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,55 @@
11
using Revise, ApproxFun, MultivariateOrthogonalPolynomials
2-
import MultivariateOrthogonalPolynomials: checkerboard, icheckerboard
2+
import MultivariateOrthogonalPolynomials: checkerboard, icheckerboard, CDisk2CxfPlan
3+
import ApproxFunBase: totensor
4+
import ApproxFunOrthogonalPolynomials: jacobip
35

46
f = (x,y) -> x*y+cos(y-0.1)+sin(x)+1;
57
ff = Fun(f, ChebyshevDisk(), 1000)
6-
@test ff(0.1,0.2) ≈ f(0.1,0.2)
78

9+
@test ff(0.1,0.2) ≈ f(0.1,0.2)
810
ff = Fun(f, ChebyshevDisk())
911
@test ff(0.1,0.2) ≈ f(0.1,0.2)
1012

13+
14+
f = Fun((x,y) -> 1, ChebyshevDisk()); n = 1;
15+
c = totensor(f.space, f.coefficients)
16+
P = CDisk2CxfPlan(n)
17+
d = P \ c
18+
@test d ≈ [1/sqrt(2)]
19+
f = Fun((x,y) -> y, ChebyshevDisk()); n = 2;
20+
c = totensor(f.space, f.coefficients)
21+
c = pad(c, n, 4n-3)
22+
P = CDisk2CxfPlan(n)
23+
d = P \ c
24+
@test d ≈ [0.0 0.5 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0]
25+
26+
f = Fun((x,y) -> x, ChebyshevDisk()); n = 2;
27+
c = totensor(f.space, f.coefficients)
28+
c = pad(c, n, 4n-3)
29+
P = CDisk2CxfPlan(n)
30+
d = P \ c
31+
@test d ≈ [0.0 0.0 0.5 0.0 0.0; 0.0 0.0 0.0 0.0 0.0]
32+
33+
34+
m,ℓ = (0,2);
35+
m,ℓ = (1,1);
36+
m,ℓ = (2,2);
37+
m,ℓ = (2,4);
38+
39+
p1 = (r,θ) -> sqrt(2ℓ+2)*r^m*jacobip((ℓ-m)÷2, m, 0, 2r^2-1)*cos(m*θ)/sqrt(2π)
40+
f = Fun((x,y) -> p1(sqrt(x^2+y^2), atan(y,x)), ChebyshevDisk());
41+
c = totensor(f.space, chop(f.coefficients,1E-12))
42+
n = size(c,1); c = sqrt(2π)*pad(c,n,4n-3)
43+
P = CDisk2CxfPlan(n)
44+
d = P \ c
45+
@test d[(ℓ-m)÷2, 2m+1] ≈ 0
46+
@test d[(ℓ-m)÷2+1, 2m+1] ≈ 1
47+
48+
49+
50+
51+
52+
53+
54+
55+

0 commit comments

Comments
 (0)