Skip to content

Commit ea84cf3

Browse files
dlfivefiftydlfivefifty
committed
v0.0.1
1 parent 6c6ffd3 commit ea84cf3

File tree

4 files changed

+250
-30
lines changed

4 files changed

+250
-30
lines changed

‎Project.toml

Lines changed: 8 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -29,16 +29,17 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
2929

3030
[compat]
3131
ApproxFun = "0.11"
32+
ApproxFunOrthogonalPolynomials = "0.2"
3233
ApproxFunBase = "0.1"
3334
BandedMatrices = "0.10"
3435
BlockArrays = "0.9"
35-
DomainSets = "0.0.2, 0.1"
36-
FastGaussQuadrature = "≥ 0.3.0"
37-
FastTransforms = "≥ 0.4.1"
38-
FillArrays = "≥ 0.5"
36+
DomainSets = "0.1"
37+
FastGaussQuadrature = "0.3"
38+
FastTransforms = "0.5"
39+
FillArrays = "0.6"
3940
InfiniteArrays = "0.1"
4041
LazyArrays = "0.10"
41-
RecipesBase = "≥ 0.5.0"
42-
SpecialFunctions = "≥ 0.6.0"
43-
StaticArrays = "≥ 0.3.0"
42+
RecipesBase = "0.5, 0.6, 0.7"
43+
SpecialFunctions = "0.6, 0.7"
44+
StaticArrays = "0.11"
4445
julia = "1"

‎README.md

Lines changed: 19 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -1,28 +1,26 @@
1-
# DiskFun
1+
# MultivariateOrthogonalPolynomials.jl
22

3-
[![Build Status](https://travis-ci.org/dlfivefifty/DiskFun.jl.svg?branch=master)](https://travis-ci.org/dlfivefifty/DiskFun.jl)
4-
5-
6-
7-
The following solves Poisson `Δu =f` with zero Dirichlet conditions
8-
on a disk
3+
[![Build Status](https://travis-ci.org/JuliaApproximation/MultivariateOrthogonalPolynomials.jl.svg?branch=master)](https://travis-ci.org/JuliaApproximation/MultivariateOrthogonalPolynomials.jl)
94

5+
This is an experimental package to add support for multivariate orthogonal polynomials on disks, spheres, triangles, and other simple
6+
geometries to [ApproxFun.jl](https://github.com/JuliaApproximation/ApproxFun.jl). At the moment it primarily supports triangles. For example,
7+
we can solve variable coefficient Helmholtz on the triangle with zero Dirichlet conditions as follows:
108
```julia
11-
d = Disk()
12-
f = Fun((x,y)->exp(-10(x+.2)^2-20(y-.1)^2),d)
13-
u = [dirichlet(d);lap(d)]\Any[0.,f]
14-
ApproxFun.plot(u) # Requires Gadfly or PyPlot
9+
using ApproxFun, MultivariateOrthogonalPolynomials
10+
x,y = Fun(Triangle())
11+
Δ = Laplacian() : TriangleWeight(1.0,1.0,1.0,JacobiTriangle(1.0,1.0,1.0))
12+
V = x*y^2
13+
L = Δ + 200^2*V
14+
u = \(L, ones(Triangle()); tolerance=1E-5)
1515
```
16+
See the examples folder for more examples, including non-zero Dirichlet conditions, Neumann conditions, and piecing together multiple triangles. In particular, the [examples](examples/triangleexamples.jl) from Olver, Townsend & Vasil 2019.
1617

1718

18-
The following solves beam equation `u_tt + Δ^2u = 0`
19-
on a disk
19+
This code relies on Slevinsky's [FastTransforms](https://github.com/MikaelSlevinsky/FastTransforms) C library for calculating transforms between values and coefficients. At the moment the path to the compiled FastTransforms library is hard coded in [c_transforms.jl](src/c_transforms.jl).
2020

21-
```julia
22-
d = Disk()
23-
u0 = Fun((x,y)->exp(-50x^2-40(y-.1)^2)+.5exp(-30(x+.5)^2-40(y+.2)^2),d)
24-
B= [dirichlet(d),neumann(d)]
25-
L = -lap(d)^2
26-
h = 0.001
27-
timeevolution(2,B,L,u0,h) # Requires GLPlot
28-
```
21+
## References
22+
23+
24+
S. Olver, A. Townsend & G.M. Vasil (2019), [A sparse spectral method on triangles](https://arxiv.org/pdf/1902.04863.pdf), arXiv:1902.04
25+
S. Olver & Y. Xuan (2019), [Orthogonal polynomials in and on a quadratic surface of revolution](https://arxiv.org/abs/1906.12305.pdf), arXiv:1906.12305
26+
G.M. Vasil, K.J. Burns, D. Lecoanet, S. Olver, B.P. Brown & J.S. Oishi (2016), [Tensor calculus in polar coordinates using Jacobi polynomials](http://arxiv.org/pdf/1509.07624.pdf), J. Comp. Phys., 325: 53–73

‎examples/helmholtz.jl

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,11 @@
1-
using ApproxFun, MultivariateOrthogonalPolynomials, BlockArrays, FillArrays
2-
using Plots
1+
using ApproxFun, MultivariateOrthogonalPolynomials, BlockArrays, FillArrays, Plots
32
import PyPlot
43
x,y = Fun(Triangle())
54
Δ = Laplacian() : TriangleWeight(1.0,1.0,1.0,JacobiTriangle(1.0,1.0,1.0))
65
V = x*y^2
76
L = Δ + 200^2*V
7+
u = L \ ones(Triangle())
8+
89
M = L[Block.(1:50),Block.(1:50)];
910

1011
PyPlot.spy(M)

‎examples/triangleexamples.jl

Lines changed: 220 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,220 @@
1+
using ApproxFun, MultivariateOrthogonalPolynomials, BlockArrays, SpecialFunctions, FillArrays, Plots
2+
import ApproxFun: blockbandwidths, Vec, PiecewiseSegment
3+
import MultivariateOrthogonalPolynomials: DirichletTriangle
4+
5+
6+
7+
#######
8+
# Example 1 Poisson
9+
#######
10+
11+
S = TriangleWeight(1,1,1,JacobiTriangle(1,1,1))
12+
Δ = Laplacian(S)
13+
14+
N = 500
15+
M = sparse(Δ[Block.(1:N), Block.(1:N)])
16+
F = lu(M)
17+
u_1 = F \ pad(coefficients(Fun(1, rangespace(Δ))), size(M,1))
18+
u_2 = F \ pad(coefficients(Fun((x,y) -> x^2 + y^2, rangespace(Δ))), size(M,1))
19+
u_3 = F \ pad(coefficients(Fun((x,y) -> x*y*(1-x-y), rangespace(Δ))), size(M,1))
20+
u_4 = F \ pad(coefficients(Fun((x,y) -> x*y*(1-x-y) *cos(300x*y), rangespace(Δ))), size(M,1))
21+
u_5 = F \ pad(coefficients(Fun((x,y) -> exp(-25((x-0.2)^2+(y-0.2)^2)), rangespace(Δ))), size(M,1))
22+
23+
gr()
24+
25+
plot(abs.(u_1)[1:10_000]; yscale=:log10, xscale=:log10, label="1")
26+
plot!(abs.(u_2)[1:10_000]; yscale=:log10, xscale=:log10, label="x^2 + y^2")
27+
plot!(abs.(u_3)[1:10_000]; yscale=:log10, xscale=:log10, label="xy(1-x-y)")
28+
plot!(abs.(u_4)[1:10_000]; yscale=:log10, xscale=:log10, label="xy(1-x-y)cos(300xy)")
29+
30+
#######
31+
# Example 2 Variable coefficient Helmholts
32+
######
33+
34+
S = TriangleWeight(1,1,1,JacobiTriangle(1,1,1))
35+
x,y = Fun(Triangle())
36+
37+
V = 1-(3*(x-1)^4 + 5*y^4)
38+
39+
40+
f = Fun((x,y) -> x*y*exp(x) , JacobiTriangle(1.,1.,1.))
41+
42+
k = 100
43+
L = Laplacian(S) + (k^2 * V)
44+
N = ceil(Int,2k)
45+
@time M̃ = L[Block.(1:N), Block.(1:N)]
46+
@time M = sparse(M̃)
47+
@time F = lu(M)
48+
@time u_c = F \ pad(coefficients(f), size(M,1))
49+
u = Fun(domainspace(L), u_c)
50+
51+
52+
#######
53+
# Example 3 Bi-harmonic
54+
#######
55+
56+
S = TriangleWeight(2,2,2, JacobiTriangle(2,2,2))
57+
Δ = Laplacian(S)
58+
L = Δ^2
59+
60+
N = 1001
61+
M = sparse(L[Block.(1:N), Block.(1:N)])
62+
F = lu(M)
63+
64+
65+
u_1 = F \ pad(coefficients(Fun(1, rangespace(L))), size(M,1))
66+
u_2 = F \ pad(coefficients(Fun((x,y) -> x^2 + y^2, rangespace(L))), size(M,1))
67+
u_3 = F \ pad(coefficients(Fun((x,y) -> x*y*(1-x-y), rangespace(L))), size(M,1))
68+
# u_4 = F \ pad(coefficients(Fun((x,y) -> x*y*(1-x-y) *cos(300x*y), rangespace(L))), size(M,1))
69+
u_4 = F \ pad(coefficients(Fun((x,y) -> exp(-25((x-0.2)^2+(y-0.2)^2)), rangespace(L))), size(M,1))
70+
71+
72+
bnrms = u -> [norm(PseudoBlockArray(u, 1:N)[Block(K)],Inf) for K=1:N]
73+
74+
plot(bnrms(u_1); xscale=:log10, yscale=:log10, linewidth=3, label="1")
75+
plot!(bnrms(u_2); xscale=:log10, yscale=:log10, linewidth=3, linestyle=:dot, label="x^2 + y^2")
76+
plot!(bnrms(u_3); xscale=:log10, yscale=:log10, linewidth=3, linestyle=:dot, label="x^2 + y^2")
77+
plot!(bnrms(u_4); xscale=:log10, yscale=:log10, linewidth=3, linestyle=:dot, label="x^2 + y^2")
78+
79+
80+
plot(abs.(coefficients(Fun((x,y) -> exp(-25((x-0.2)^2+(y-0.2)^2)), rangespace(L)))); scale=:log10, yscale=:log10)
81+
82+
#######
83+
# Example 4 Laplace's equation
84+
#######
85+
86+
S = DirichletTriangle{1,1,1}()
87+
B = Dirichlet() : S
88+
Δ = Laplacian() : S
89+
# L = [B; Δ]
90+
91+
fs = ((x,y) -> exp(x)*cos(y),
92+
(x,y) -> x^2,
93+
(x,y) -> x^3*(1-x)^3*(1-y)^3,
94+
(x,y) -> x^3*(1-x)^3*(1-y)^3 * cos(100x*y))
95+
96+
N = 1000;
97+
@time BM = sparse(B[Block.(1:N), Block.(1:N)]);
98+
@time ΔM = sparse(Δ[Block.(1:N), Block.(1:N)]);
99+
@time M = [BM; ΔM];
100+
M = [sparse(I, size(M,1), 2) M] # add tau term
101+
@time F = factorize(M);
102+
103+
# RHS
104+
k = 1
105+
∂u = Fun(fs[k], rangespace(B))
106+
r = pad(coefficients(∂u, rangespace(B)), size(Fs[N],1))
107+
@time u_cfs = F \ r
108+
109+
110+
#######
111+
# Example 5 Transport equation
112+
#######
113+
114+
S = DirichletTriangle{0,1,0}()
115+
∂S = Legendre(Segment(Vec(0.0,0), Vec(1.0,0)))
116+
B = I : S → ∂S
117+
118+
Dx = Derivative([1,0]) : S
119+
Dt = Derivative([0,1]) : S
120+
121+
122+
x,_ = Fun(∂S)
123+
u₀ = x*(1-x)*exp(x)
124+
125+
# u_t = u_x
126+
L = [B; Dx - Dt]
127+
N = 20; M = sparse(L[Block.(1:N), Block.(1:N)]);
128+
u_c = M \ pad(coefficients([u₀; 0], rangespace(L)), size(M,1))
129+
u = Fun(domainspace(L), u_c)
130+
131+
# 0.5u_t = u_x
132+
S = DirichletTriangle{0,1,1}()
133+
∂S = Legendre(Segment(Vec(0.0,0), Vec(1.0,0))) , Legendre(Segment(Vec(0.0,1.0),Vec(1.0,0)))
134+
B_x = I : S → ∂S[1]
135+
B_z = I : S → ∂S[2]
136+
Dx = Derivative([1,0]) : S
137+
Dt = Derivative([0,1]) : S
138+
139+
L = [B_x; B_z; Dx - 0.5Dt]
140+
N = 100; M = sparse(L[Block.(1:N), Block.(1:N)])
141+
u₀_x = Fun((x,_) -> x*exp(x-1), rangespace(B_x))
142+
u₀_z = Fun((x,y) -> (1-y), rangespace(B_z))
143+
u_c = M \ pad(coefficients([u₀_x; u₀_z; 0], rangespace(L)), size(M,1))
144+
u = Fun(domainspace(L), u_c)
145+
146+
# u_t + u_x = 0
147+
S = DirichletTriangle{1,1,0}()
148+
∂S = Legendre(Segment(Vec(0.0,0), Vec(1.0,0))) , Legendre(Segment(Vec(0.0,0.0),Vec(0.0,1.0)))
149+
B_x = I : S → ∂S[1]
150+
B_y = I : S → ∂S[2]
151+
Dx = Derivative([1,0]) : S
152+
Dt = Derivative([0,1]) : S
153+
154+
L = [B_x; B_y; Dt + Dx]
155+
N = 100; M = sparse(L[Block.(1:N), Block.(1:N)])
156+
u₀_x = Fun((x,_) -> (1-x)*exp(x), rangespace(B_x))
157+
u₀_y = Fun((x,y) -> (1-y), rangespace(B_y))
158+
159+
u_c = M \ pad(coefficients([u₀_x; u₀_y; 0], rangespace(L)), size(M,1))
160+
u = Fun(domainspace(L), u_c)
161+
162+
163+
#######
164+
# Example 7 Helmholtz in a polygon
165+
#######
166+
d = Triangle() ∪ Triangle(Vec(1,1),Vec(1,0),Vec(0,1)) ∪
167+
Triangle(Vec(1,1),Vec(0,1),Vec(0,2)) ∪
168+
Triangle(Vec(0,0),Vec(0,1),Vec(-1,1.5))
169+
170+
171+
f = Fun((x,y) -> cos(x*sin(y)), d)
172+
173+
∂d = PiecewiseSegment([Vec(0.,0), Vec(1.,0), Vec(1,1), Vec(0,2), Vec(0,1), Vec(-1,1.5), Vec(0,0)])
174+
175+
S = DirichletTriangle{1,1,1}.(components(d))
176+
∂S = Legendre.(components(∂d))
177+
178+
∂S12 = Legendre(Segment(Vec(1.,0),Vec(0,1)))
179+
∂S23 = Legendre(Segment(Vec(1.,1),Vec(0,1)))
180+
∂S14 = Legendre(Segment(Vec(0.,0),Vec(0,1)))
181+
182+
Dx = Derivative([1,0])
183+
Dy = Derivative([0,1])
184+
185+
k = 20
186+
187+
L = [(I : S[1] → ∂S[1]) 0 0 0 0 0 0 0 0 0 0 0;
188+
0 0 0 (I : S[2] → ∂S[2]) 0 0 0 0 0 0 0 0;
189+
0 0 0 0 0 0 (I : S[3] → ∂S[3]) 0 0 0 0 0;
190+
0 0 0 0 0 0 (I : S[3] → ∂S[4]) 0 0 0 0 0;
191+
0 0 0 0 0 0 0 0 0 (I : S[4] → ∂S[5]) 0 0;
192+
0 0 0 0 0 0 0 0 0 (I : S[4] → ∂S[6]) 0 0;
193+
(I : S[1] → ∂S12) 0 0 -(I : S[2] → ∂S12) 0 0 0 0 0 0 0 0;
194+
0 (I : S[1] → ∂S12) (I : S[1] → ∂S12) 0 -(I : S[2] → ∂S12) -(I : S[2] → ∂S12) 0 0 0 0 0 0;
195+
0 0 0 (I : S[2] → ∂S23) 0 0 -(I : S[3] → ∂S23) 0 0 0 0 0;
196+
0 0 0 0 0 (I : S[2] → ∂S23) 0 0 -(I : S[3] → ∂S23) 0 0 0;
197+
(I : S[1] → ∂S14) 0 0 0 0 0 0 0 0 -(I : S[4] → ∂S14) 0 0;
198+
0 (I : S[1] → ∂S14) 0 0 0 0 0 0 0 0 -(I : S[4] → ∂S14) 0;
199+
Dx -I 0 0 0 0 0 0 0 0 0 0;
200+
Dy 0 -I 0 0 0 0 0 0 0 0 0;
201+
0 0 0 Dx -I 0 0 0 0 0 0 0;
202+
0 0 0 Dy 0 -I 0 0 0 0 0 0;
203+
0 0 0 0 0 0 (Dx : S[3]) -(I : S[3]) 0 0 0 0;
204+
0 0 0 0 0 0 Dy 0 -I 0 0 0;
205+
0 0 0 0 0 0 0 0 0 Dx -I 0;
206+
0 0 0 0 0 0 0 0 0 (Dy : S[4]) 0 -(I : S[4])
207+
k^2*I Dx Dy 0 0 0 0 0 0 0 0 0;
208+
0 0 0 k^2*I Dx Dy 0 0 0 0 0 0;
209+
0 0 0 0 0 0 k^2*I Dx Dy 0 0 0;
210+
0 0 0 0 0 0 0 0 0 k^2*I Dx Dy
211+
]
212+
213+
214+
f = ones(∂d)
215+
216+
N = 101
217+
@time M = sparse(L[Block.(1:N), Block.(1:N)]);
218+
@time F = qr(M);
219+
@time u_c = F \ pad(coefficients(vcat(components(f)..., Zeros(18)...), rangespace(L)), size(M,1));
220+
u = Fun(domainspace(L), u_c)

0 commit comments

Comments
 (0)