Skip to content

Commit 09d78a1

Browse files
dlfivefiftydlfivefifty
committed
Cone -> Conic
1 parent aee1d90 commit 09d78a1

File tree

3 files changed

+207
-59
lines changed

3 files changed

+207
-59
lines changed

‎Project.toml

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -26,14 +26,14 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
2626

2727
[compat]
2828
ApproxFun = "0.11"
29-
ApproxFunBase = "≥ 0.0.3"
29+
ApproxFunBase = "0.1"
3030
BandedMatrices = "0.9"
31-
BlockArrays = "0.8"
32-
DomainSets = "≥ 0.0.1"
31+
BlockArrays = "0.9"
32+
DomainSets = "0.0.2"
3333
FastGaussQuadrature = "≥ 0.3.0"
3434
FastTransforms = "≥ 0.4.1"
3535
FillArrays = "≥ 0.5"
36-
InfiniteArrays = "0.0.3"
36+
InfiniteArrays = "0.1"
3737
LazyArrays = "0.8"
3838
RecipesBase = "≥ 0.5.0"
3939
SpecialFunctions = "≥ 0.6.0"

‎src/Cone/Cone.jl

Lines changed: 157 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,153 @@
1-
export Cone, DuffyCone, LegendreCone
1+
export Cone, DuffyCone, LegendreCone, Conic, DuffyConic, LegendreConic
2+
3+
###
4+
# Conic
5+
###
6+
7+
struct Conic <: Domain{Vec{3,Float64}} end
8+
struct DuffyConic <: Space{Conic,Float64} end
9+
struct LegendreConic <: Space{Conic,Float64} end
10+
11+
@containsconstants DuffyConic
12+
@containsconstants LegendreConic
13+
14+
spacescompatible(::DuffyConic, ::DuffyConic) = true
15+
spacescompatible(::LegendreConic, ::LegendreConic) = true
16+
17+
domain(::DuffyConic) = Conic()
18+
domain(::LegendreConic) = Conic()
19+
20+
tensorizer(K::DuffyConic) = DiskTensorizer()
21+
tensorizer(K::LegendreConic) = DiskTensorizer()
22+
23+
rectspace(::DuffyConic) = NormalizedJacobi(0,1,Segment(1,0))*Fourier()
24+
25+
conemap(t,x,y) = Vec(t, t*x, t*y)
26+
conemap(txy) = ((t,x,y) = txy; conemap(t,x,y))
27+
conicpolar(t,θ) = conemap(t, cos(θ), sin(θ))
28+
conicpolar(tθ) = ((t,θ) = tθ; conicpolar(t,θ))
29+
30+
points(::Conic, n) = conicpolar.(points(rectspace(DuffyConic()), n))
31+
checkpoints(::Conic) = conicpolar.(checkpoints(rectspace(DuffyConic())))
32+
33+
34+
plan_transform(S::DuffyConic, n::AbstractVector) =
35+
TransformPlan(S, plan_transform(rectspace(S),n), Val{false})
36+
plan_itransform(S::DuffyConic, n::AbstractVector) =
37+
ITransformPlan(S, plan_itransform(rectspace(S),n), Val{false})
38+
39+
*(P::TransformPlan{<:Any,<:DuffyConic}, v::AbstractArray) = P.plan*v
40+
*(P::ITransformPlan{<:Any,<:DuffyConic}, v::AbstractArray) = P.plan*v
41+
42+
evaluate(cfs::AbstractVector, S::DuffyConic, txy) = evaluate(cfs, rectspace(S), ipolar(txy[Vec(2,3)]))
43+
evaluate(cfs::AbstractVector, S::LegendreConic, txy) = evaluate(coefficients(cfs, S, DuffyConic()), DuffyConic(), txy)
44+
45+
46+
# function transform(::DuffyConic, v)
47+
# n = (1 + isqrt(1+8*length(v))) ÷ 4
48+
# F = copy(reshape(v,n,2n-1))
49+
# Pt = plan_transform(NormalizedJacobi(0,1,Segment(1,0)), n)
50+
# Pθ = plan_transform(Fourier(), 2n-1)
51+
# for j = 1:size(F,2)
52+
# F[:,j] = Pt*F[:,j]
53+
# end
54+
# for k = 1:size(F,1)
55+
# F[k,:] = Pθ * F[k,:]
56+
# end
57+
# fromtensor(rectspace(DuffyConic()), F)
58+
# end
59+
60+
function evaluate(cfs::AbstractVector, S::DuffyConic, txy::Vec{3})
61+
t,x,y = txy
62+
@assert t ≈ sqrt(x^2+y^2)
63+
θ = atan(y,x)
64+
Fun(rectspace(S), cfs)(t,θ)
65+
end
66+
67+
68+
function _coefficients(triangleplan, v::AbstractVector{T}, ::DuffyConic, ::LegendreConic) where T
69+
F = totensor(rectspace(DuffyConic()), v)
70+
F = pad(F, :, 2size(F,1)-1)
71+
Fc = F[:,1:2:end]
72+
c_execute_tri_lo2hi(triangleplan, Fc)
73+
F[:,1:2:end] .= Fc
74+
# ignore first column
75+
Fc[:,2:end] .= F[:,2:2:end]
76+
c_execute_tri_lo2hi(triangleplan, Fc)
77+
F[:,2:2:end] .= Fc[:,2:end]
78+
fromtensor(LegendreConic(), F)
79+
end
80+
81+
function _coefficients(triangleplan, v::AbstractVector{T}, ::LegendreConic, ::DuffyConic) where T
82+
F = totensor(LegendreConic(), v)
83+
F = pad(F, :, 2size(F,1)-1)
84+
Fc = F[:,1:2:end]
85+
c_execute_tri_hi2lo(triangleplan, Fc)
86+
F[:,1:2:end] .= Fc
87+
# ignore first column
88+
Fc[:,2:end] .= F[:,2:2:end]
89+
c_execute_tri_hi2lo(triangleplan, Fc)
90+
F[:,2:2:end] .= Fc[:,2:end]
91+
fromtensor(DuffyConic(), F)
92+
end
93+
94+
function coefficients(cfs::AbstractVector{T}, CD::DuffyConic, ZD::LegendreConic) where T
95+
c = totensor(rectspace(DuffyConic()), cfs) # TODO: wasteful
96+
n = size(c,1)
97+
_coefficients(c_plan_rottriangle(n, zero(T), zero(T), zero(T)), cfs, CD, ZD)
98+
end
99+
100+
function coefficients(cfs::AbstractVector{T}, ZD::LegendreConic, CD::DuffyConic) where T
101+
c = tridevec(cfs) # TODO: wasteful
102+
n = size(c,1)
103+
_coefficients(c_plan_rottriangle(n, zero(T), zero(T), zero(T)), cfs, ZD, CD)
104+
end
105+
106+
struct LegendreConicTransformPlan{DUF,CHEB}
107+
duffyplan::DUF
108+
triangleplan::CHEB
109+
end
110+
111+
function LegendreConicTransformPlan(S::LegendreConic, v::AbstractVector{T}) where T
112+
D = plan_transform(DuffyConic(),v)
113+
F = totensor(rectspace(DuffyConic()), D*v) # wasteful, just use to figure out `size(F,1)`
114+
P = c_plan_rottriangle(size(F,1), zero(T), zero(T), zero(T))
115+
LegendreConicTransformPlan(D,P)
116+
end
117+
118+
119+
*(P::LegendreConicTransformPlan, v::AbstractVector) = _coefficients(P.triangleplan, P.duffyplan*v, DuffyConic(), LegendreConic())
120+
121+
plan_transform(K::LegendreConic, v::AbstractVector) = LegendreConicTransformPlan(K, v)
122+
123+
struct LegendreConicITransformPlan{DUF,CHEB}
124+
iduffyplan::DUF
125+
triangleplan::CHEB
126+
end
127+
128+
function LegendreConicITransformPlan(S::LegendreConic, v::AbstractVector{T}) where T
129+
F = fromtensor(LegendreConic(), v) # wasteful, just use to figure out `size(F,1)`
130+
P = c_plan_rottriangle(size(F,1), zero(T), zero(T), zero(T))
131+
D = plan_itransform(DuffyConic(), _coefficients(P, v, S, DuffyConic()))
132+
LegendreConicITransformPlan(D,P)
133+
end
134+
135+
136+
*(P::LegendreConicITransformPlan, v::AbstractVector) = P.iduffyplan*_coefficients(P.triangleplan, v, LegendreConic(), DuffyConic())
137+
138+
139+
plan_itransform(K::LegendreConic, v::AbstractVector) = LegendreConicITransformPlan(K, v)
140+
141+
142+
143+
###
144+
# Cone
145+
###
2146

3147
struct Cone <: Domain{Vec{3,Float64}} end
148+
149+
in(x::Vec{3}, d::Cone) = 0 ≤ x[1] ≤ 1 && x[2]^2 + x[3]^2 ≤ x[1]^2
150+
4151
struct DuffyCone <: Space{Cone,Float64} end
5152
struct LegendreCone <: Space{Cone,Float64} end
6153

@@ -13,16 +160,16 @@ spacescompatible(::LegendreCone, ::LegendreCone) = true
13160
domain(::DuffyCone) = Cone()
14161
domain(::LegendreCone) = Cone()
15162

16-
tensorizer(K::DuffyCone) = DiskTensorizer()
17-
tensorizer(K::LegendreCone) = DiskTensorizer()
163+
tensorizer(K::DuffyCone) = BallTensorizer()
164+
tensorizer(K::LegendreCone) = BallTensorizer()
18165

19-
rectspace(::DuffyCone) = NormalizedJacobi(0,1,Segment(1,0))*Fourier()
166+
rectspace(::DuffyCone) = NormalizedJacobi(0,1,Segment(1,0))*ZernikeDisk()
20167

21-
conepolar(t,θ) = Vec(t, t*cos(θ), t*sin(θ))
22-
conepolar(tθ) = ((t,θ) = tθ; conepolar(t,θ))
23168

24-
points(::Cone, n) = conepolar.(points(rectspace(DuffyCone()), n))
25-
checkpoints(::Cone) = conepolar.(checkpoints(rectspace(DuffyCone())))
169+
170+
171+
points(::Cone, n) = conemap.(points(rectspace(DuffyCone()), n))
172+
checkpoints(::Cone) = conemap.(checkpoints(rectspace(DuffyCone())))
26173

27174

28175
plan_transform(S::DuffyCone, n::AbstractVector) =
@@ -53,9 +200,8 @@ evaluate(cfs::AbstractVector, S::LegendreCone, txy) = evaluate(coefficients(cfs,
53200

54201
function evaluate(cfs::AbstractVector, S::DuffyCone, txy::Vec{3})
55202
t,x,y = txy
56-
@assert t ≈ sqrt(x^2+y^2)
57-
θ = atan(y,x)
58-
Fun(rectspace(S), cfs)(t,θ)
203+
@assert x^2+y^2 ≤ t^2
204+
Fun(rectspace(S), cfs)(t,x/t,y/t)
59205
end
60206

61207

‎test/test_cone.jl

Lines changed: 46 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -1,63 +1,65 @@
11
using ApproxFun, MultivariateOrthogonalPolynomials, Test
22
import MultivariateOrthogonalPolynomials: rectspace, totensor
3+
import ApproxFunBase: plan_transform
34

4-
@testset "DuffyCone" begin
5-
f = Fun((t,x,y) -> 1, DuffyCone(), 10)
6-
@test f.coefficients ≈ [1; zeros(ncoefficients(f)-1)]
7-
@test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ 1
5+
@testset "Conic" begin
6+
@testset "DuffyConic" begin
7+
f = Fun((t,x,y) -> 1, DuffyConic(), 10)
8+
@test f.coefficients ≈ [1; zeros(ncoefficients(f)-1)]
9+
@test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ 1
810

9-
f = Fun((t,x,y) -> t, DuffyCone(), 10)
10-
@test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ sqrt(0.1^2+0.2^2)
11+
f = Fun((t,x,y) -> t, DuffyConic(), 10)
12+
@test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ sqrt(0.1^2+0.2^2)
1113

12-
f = Fun((t,x,y) -> x, DuffyCone(), 10)
14+
f = Fun((t,x,y) -> x, DuffyConic(), 10)
1315

14-
@test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ 0.1
16+
@test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ 0.1
1517

16-
f = Fun((t,x,y) -> exp(cos(t*x)*y), DuffyCone(), 1000)
17-
@test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ exp(cos(sqrt(0.1^2+0.2^2)*0.1)*0.2)
18+
f = Fun((t,x,y) -> exp(cos(t*x)*y), DuffyConic(), 1000)
19+
@test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ exp(cos(sqrt(0.1^2+0.2^2)*0.1)*0.2)
1820

19-
f = Fun((t,x,y) -> exp(cos(t*x)*y), DuffyCone())
20-
@test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ exp(cos(sqrt(0.1^2+0.2^2)*0.1)*0.2)
21+
f = Fun((t,x,y) -> exp(cos(t*x)*y), DuffyConic())
22+
@test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ exp(cos(sqrt(0.1^2+0.2^2)*0.1)*0.2)
2123

22-
m,ℓ = (1,1)
23-
f = (txy) -> ((t,x,y) = txy; θ = atan(y,x); Fun(NormalizedJacobi(0,2m+1,Segment(1,0)),[zeros(ℓ);1])(t) * 2^m * t^m * cos(m*θ))
24-
g = Fun(f, DuffyCone())
25-
t,x,y = sqrt(0.1^2+0.2^2),0.1,0.2
26-
@test g(t,x,y) ≈ f((t,x,y))
27-
end
24+
m,ℓ = (1,1)
25+
f = (txy) -> ((t,x,y) = txy; θ = atan(y,x); Fun(NormalizedJacobi(0,2m+1,Segment(1,0)),[zeros(ℓ);1])(t) * 2^m * t^m * cos(m*θ))
26+
g = Fun(f, DuffyConic())
27+
t,x,y = sqrt(0.1^2+0.2^2),0.1,0.2
28+
@test g(t,x,y) ≈ f((t,x,y))
29+
end
2830

29-
@testset "LegendreConePlan" begin
30-
m,ℓ = (1,1)
31-
f = (txy) -> ((t,x,y) = txy; θ = atan(y,x); Fun(NormalizedJacobi(0,2m+1,Segment(1,0)),[zeros(ℓ);1])(t) * 2^m * t^m * cos(m*θ))
32-
p = points(LegendreCone(), 10)
33-
P = plan_transform(LegendreCone(), f.(p))
34-
@test P.duffyplan*f.(p) ≈ Fun(f, DuffyCone()).coefficients[1:12]
35-
coefficients(g, LegendreCone())
36-
g = Fun(f, LegendreCone(), 20)
37-
t,x,y = sqrt(0.1^2+0.2^2),0.1,0.2
38-
@test g(t,x,y) ≈ f((t,x,y))
39-
end
31+
@testset "LegendreConicPlan" begin
32+
m,ℓ = (1,1)
33+
f = (txy) -> ((t,x,y) = txy; θ = atan(y,x); Fun(NormalizedJacobi(0,2m+1,Segment(1,0)),[zeros(ℓ);1])(t) * 2^m * t^m * cos(m*θ))
34+
p = points(LegendreConic(), 10)
35+
P = plan_transform(LegendreConic(), f.(p))
36+
@test P.duffyplan*f.(p) ≈ Fun(f, DuffyConic()).coefficients[1:12]
37+
g = Fun(f, LegendreConic(), 20)
38+
t,x,y = sqrt(0.1^2+0.2^2),0.1,0.2
39+
@test g(t,x,y) ≈ f((t,x,y))
40+
end
4041

41-
@testset "Legendre<>DuffyCone" begin
42-
a = randn(10)
43-
F = totensor(rectspace(DuffyCone()), a)
42+
@testset "Legendre<>DuffyConic" begin
43+
a = randn(10)
44+
F = totensor(rectspace(DuffyConic()), a)
4445

4546

46-
b = coefficients(a, LegendreCone(), DuffyCone())
47+
b = coefficients(a, LegendreConic(), DuffyConic())
4748

48-
coefficients(b, DuffyCone(), LegendreCone())
49-
end
49+
coefficients(b, DuffyConic(), LegendreConic())
50+
end
5051

51-
@testset "LegendreCone" begin
52-
f = Fun((t,x,y) -> 1, LegendreCone(), 10)
53-
@test f.coefficients ≈ [1; zeros(ncoefficients(f)-1)]
54-
@test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ 1
52+
@testset "LegendreConic" begin
53+
f = Fun((t,x,y) -> 1, LegendreConic(), 10)
54+
@test f.coefficients ≈ [1; zeros(ncoefficients(f)-1)]
55+
@test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ 1
5556

56-
f = Fun((t,x,y) -> t, LegendreCone(), 10)
57-
@test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ sqrt(0.1^2+0.2^2)
58-
f = Fun((t,x,y) -> 1+t+x+y, LegendreCone(), 10)
59-
@test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ 1+sqrt(0.1^2+0.2^2)+0.1+0.2
57+
f = Fun((t,x,y) -> t, LegendreConic(), 10)
58+
@test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ sqrt(0.1^2+0.2^2)
59+
f = Fun((t,x,y) -> 1+t+x+y, LegendreConic(), 10)
60+
@test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ 1+sqrt(0.1^2+0.2^2)+0.1+0.2
6061

6162

62-
@time Fun((t,x,y) -> 1+t+x+y, LegendreCone(), 1000)
63+
@time Fun((t,x,y) -> 1+t+x+y, LegendreConic(), 1000)
64+
end
6365
end

0 commit comments

Comments
 (0)