Skip to content

Commit e4c09f8

Browse files
dlfivefiftydlfivefifty
committed
L-shaped domain working
1 parent d082484 commit e4c09f8

File tree

2 files changed

+333
-14
lines changed

2 files changed

+333
-14
lines changed

‎examples/L-shaped Laplace.jl

Lines changed: 332 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -1,32 +1,351 @@
1-
using ApproxFun, MultivariateOrthogonalPolynomials
2-
import ApproxFun: Vec, PiecewiseSegment, ZeroOperator
1+
using ApproxFun, MultivariateOrthogonalPolynomials, Plots, SparseArrays
2+
import ApproxFun: Vec, PiecewiseSegment, ZeroOperator, Block
33
import MultivariateOrthogonalPolynomials: DirichletTriangle
44

5-
d = [Triangle(Vec(0,0), Vec(1,0), Vec(0,1)) , Triangle(Vec(1,1),Vec(1,0),Vec(0,1)) ,
6-
Triangle(Vec(1,1),Vec(0,1),Vec(0,2)) , Triangle(Vec(1,2),Vec(0,2),Vec(1,1)) ,
5+
d = [Triangle(Vec(1,1),Vec(0,1),Vec(0,2)) , Triangle(Vec(1,2),Vec(0,2),Vec(1,1)),
6+
Triangle(Vec(0,0), Vec(1,0), Vec(0,1)) , Triangle(Vec(1,1),Vec(1,0),Vec(0,1)) ,
77
Triangle(Vec(1,0), Vec(2,0), Vec(1,1)) , Triangle(Vec(2,1), Vec(1,1), Vec(2,0))]
88

9-
∂d = components(PiecewiseSegment([Vec(0,0), Vec(1,0), Vec(2,0), Vec(2,1), Vec(1,1), Vec(1,2), Vec(0,2), Vec(0,1), Vec(0,0)]))
10-
ιd = [Segment(Vec(0,1),Vec(1,0)), Segment(Vec(0,1), Vec(1,1)), Segment(Vec(1,1), Vec(2,0)),
11-
Segment(Vec(1,0), Vec(1,1)), Segment(Vec(0,2), Vec(1,1))]
129

13-
length(∂d)
10+
p = plot()
11+
for ▴ in d
12+
plot!(▴)
13+
end
14+
p
15+
16+
∂d = components(PiecewiseSegment([Vec(0,2), Vec(0,1), Vec(0,0), Vec(1,0), Vec(2,0), Vec(2,1), Vec(1,1), Vec(1,2), Vec(0,2)]))
17+
18+
for s in ∂d
19+
plot!(s)
20+
end
21+
p
22+
23+
# ιd = [Segment(Vec(0,2), Vec(1,1)), Segment(Vec(0,1), Vec(1,1)),
24+
# Segment(Vec(0,1),Vec(1,0)), Segment(Vec(1,0), Vec(1,1)),
25+
# Segment(Vec(1,1), Vec(2,0))]
26+
27+
28+
# interfaces
29+
ιd = Dict{NTuple{2,Int}, Segment{Vec{2,Int}}}()
30+
ιd[(1,2)] = Segment(Vec(0,2), Vec(1,1))
31+
ιd[(1,4)] = Segment(Vec(0,1), Vec(1,1))
32+
ιd[(3,4)] = Segment(Vec(0,1),Vec(1,0))
33+
ιd[(4,5)] = Segment(Vec(1,0), Vec(1,1))
34+
ιd[(5,6)] = Segment(Vec(1,1), Vec(2,0))
35+
36+
37+
keys(ιd)
38+
1439

1540
ds = vcat(fill.(DirichletTriangle{1,1,1}.(d),3)...) # repeat each triangle 3 times
16-
rs = [Legendre.(∂d); fill.(Legendre.(ιd),2)...; fill.(JacobiTriangle.(d),3)...]
41+
rs = [Legendre.(∂d)...,
42+
Legendre(ιd[1,4]), Legendre(ιd[4,5]), # straight interface
43+
Legendre(ιd[1,2]), Legendre(ιd[3,4]), Legendre(ιd[5,6]),
44+
Legendre(ιd[1,4]), Legendre(ιd[4,5]), # straight interface
45+
Legendre(ιd[1,2]), Legendre(ιd[3,4]), Legendre(ιd[5,6]),
46+
vcat(fill.(JacobiTriangle.(d),3)...)...] # diagonal interface
47+
48+
49+
50+
# rs = [Legendre.(∂d); fill.(Legendre.(values(ιd)),2)...; fill.(JacobiTriangle.(d),3)...]
51+
52+
ui = T -> 1 + (T-1)*3
53+
ux = T -> 2 + (T-1)*3
54+
uy = T -> 3 + (T-1)*3
55+
56+
Dx = Derivative([1,0])
57+
Dy = Derivative([0,1])
58+
59+
N = 100
60+
sprs = A -> (global N; sparse(A[Block.(1:N), Block.(1:N)]))
61+
62+
A = Matrix{SparseMatrixCSC{Float64,Int}}(undef, length(rs), length(ds))
63+
for K=1:length(rs), J=1:length(ds) # fill with zeros
64+
A[K,J] = ZeroOperator(ds[J],rs[K]) |> sprs
65+
end
66+
# add boundary conditions
67+
K = 0;
68+
K +=1; T = 1; J = ui(T); A[K,J] = (I : ds[J] → rs[K]) |> sprs
69+
K +=1; T = 3; J = ui(T); A[K,J] = (I : ds[J] → rs[K]) |> sprs
70+
K +=1; T = 3; J = ui(T); A[K,J] = (I : ds[J] → rs[K]) |> sprs
71+
K +=1; T = 5; J = ui(T); A[K,J] = (I : ds[J] → rs[K]) |> sprs
72+
K +=1; T = 6; J = ui(T); A[K,J] = (I : ds[J] → rs[K]) |> sprs
73+
K +=1; T = 6; J = ui(T); A[K,J] = (I : ds[J] → rs[K]) |> sprs
74+
K +=1; T = 2; J = ui(T); A[K,J] = (I : ds[J] → rs[K]) |> sprs
75+
K +=1; T = 2; J = ui(T); A[K,J] = (I : ds[J] → rs[K]) |> sprs
76+
# add dirichlet interface conditions
77+
for (T1,T2) in ((1,4), (4,5), (1,2), (3,4), (5,6))
78+
global K +=1;
79+
J = ui(T1); A[K,J] = (I : ds[J] → rs[K]) |> sprs
80+
J = ui(T2); A[K,J] = -(I : ds[J] → rs[K]) |> sprs
81+
end
82+
# add lr neumann
83+
K +=1;
84+
T = 1; J = uy(T); A[K,J] = (I : ds[J] → rs[K]) |> sprs
85+
T = 4; J = uy(T); A[K,J] = -(I : ds[J] → rs[K]) |> sprs
86+
# add ud neumann
87+
K +=1;
88+
T = 4; J = ux(T); A[K,J] = (I : ds[J] → rs[K]) |> sprs
89+
T = 5; J = ux(T); A[K,J] = -(I : ds[J] → rs[K]) |> sprs
90+
# add diagonal neumann
91+
for (T1,T2) in ((1,2), (3,4), (5,6))
92+
global K +=1;
93+
J = ux(T1); A[K,J] = (I : ds[J] → rs[K]) |> sprs
94+
J = uy(T1); A[K,J] = (I : ds[J] → rs[K]) |> sprs
95+
J = ux(T2); A[K,J] = -(I : ds[J] → rs[K]) |> sprs
96+
J = uy(T2); A[K,J] = -(I : ds[J] → rs[K]) |> sprs
97+
end
98+
99+
for T in 1:length(d)
100+
@show K,T
101+
global K +=1;
102+
J = ui(T); A[K,J] = (Dx : ds[J] → rs[K]) |> sprs
103+
J = ux(T); A[K,J] = -(I : ds[J] → rs[K]) |> sprs
104+
global K +=1;
105+
J = ui(T); A[K,J] = (Dy : ds[J] → rs[K]) |> sprs
106+
J = uy(T); A[K,J] = -(I : ds[J] → rs[K]) |> sprs
107+
global K +=1;
108+
J = ux(T); A[K,J] = (Dx : ds[J] → rs[K]) |> sprs
109+
J = uy(T); A[K,J] = (Dy : ds[J] → rs[K]) |> sprs
110+
end
111+
112+
M = hvcat(ntuple(_ -> size(A,2),size(A,1)), permutedims(A)...)
113+
114+
115+
rhs = vcat(coefficients.(Fun.(Ref((x,y) -> real(exp(x+im*y))), rs[1:length(∂d)], N))...)
116+
117+
rhs = vcat(coefficients.(Fun.(Ref((x,y) -> x^2), rs[1:length(∂d)], N))...)
118+
119+
120+
F = factorize(M)
121+
u_cfs = F \ pad(rhs, size(M,1))
122+
123+
u1 = Fun(ds[4], u_cfs[(4-1)*sum(1:N)+1:4*sum(1:N)])
124+
u1(0.99,0.99)
125+
126+
u1.coefficients
127+
u1(0.1,1.2)-real(exp(0.1+im*1.2))
128+
129+
plot(abs.([norm((M*u_cfs - pad(rhs, size(M,1)))[N*(K-1)+1:N*K]) for K=1:length(rs)] ); yscale=:log10)
130+
131+
length(rs)
132+
K = 30; norm((M*u_cfs - pad(rhs, size(M,1)))[N*(K-1)+1:N*K])
133+
134+
135+
U = Fun.(Ref((x,y) -> real(exp(x+im*y))), d, sum(1:N))
136+
137+
138+
u_cfs = Vector{Float64}()
139+
for T in d
140+
append!(u_cfs, pad(Fun((x,y) -> exp(x)*cos(y), DirichletTriangle{1,1,1}(T)).coefficients, sum(1:N)))
141+
append!(u_cfs, pad(Fun((x,y) -> exp(x)*cos(y), DirichletTriangle{1,1,1}(T)).coefficients, sum(1:N)))
142+
append!(u_cfs, pad(Fun((x,y) -> -exp(x)*sin(y), DirichletTriangle{1,1,1}(T)).coefficients, sum(1:N)))
143+
end
144+
145+
u_cfs
146+
147+
((M*u_cfs) - pad(rhs, size(M,1))) |> norm
148+
((M*u_cfs) - pad(rhs, size(M,1)))[1:(N*(length(rs)-6*3))] |> norm
149+
NN = (N*(length(rs)-6*3))
150+
((M*u_cfs) - pad(rhs, size(M,1)))[NN+1:NN+sum(1:N)] |> norm
151+
((M*u_cfs) - pad(rhs, size(M,1)))[NN+sum(1:N):NN+2sum(1:N)] |> norm
152+
153+
154+
T = d[1]
155+
u1 = pad(Fun((x,y) -> exp(x)*cos(y), DirichletTriangle{1,1,1}(T)).coefficients, sum(1:N))
156+
u1y = pad(Fun((x,y) -> -exp(x)*sin(y), DirichletTriangle{1,1,1}(T)).coefficients, sum(1:N))
157+
158+
159+
160+
A[19,1]*u1 + A[19,3]*u1y
161+
162+
NN
163+
164+
165+
size(M)
166+
size(u_cfs)
167+
size(rhs)
168+
169+
M*u_cfs - pad(rhs,size(M,1)) |> norm
170+
171+
(Dy*u1)(0.1,1.2)
172+
173+
Fun((x,y) -> -exp(x)*sin(y), DirichletTriangle{1,1,1}(T))(0.1,1.2)
174+
175+
176+
177+
178+
179+
A[18,1]*u_cfs[1:sum(1:N)] +
180+
A[18,2]*u_cfs[sum(1:N)+1:2sum(1:N)]
181+
å
182+
(length(rs)-6*3)
183+
184+
length(rs)
185+
186+
6*
187+
188+
189+
190+
(I : ds[1] → rs[1])[Block.(1:N), Block.(1:N)] * u_cfs[1:210]
191+
192+
17193

18194

195+
196+
197+
198+
199+
200+
201+
202+
203+
204+
205+
M[1:20,1:210]u_cfs[1:210]
206+
207+
Matrix(M[1:20,211:end]) |> norm
208+
209+
A[1,:]
210+
211+
212+
213+
214+
215+
216+
217+
218+
219+
220+
221+
222+
223+
224+
225+
226+
227+
u1.coefficients
228+
229+
230+
## Interlace operator
231+
232+
233+
234+
M
235+
19236
N,M = length(rs), length(ds)
20237
A = Matrix{Operator{Float64}}(undef, N,M)
21238
for K=1:N, J=1:M # fill with zeros
22239
A[K,J] = ZeroOperator(ds[J],rs[K])
23240
end
24241
# add boundary conditions
25-
for K = 1:length(∂d)
26-
for J = 1:length(d)
27-
if
28-
A[K,J] = I : ds[J] → rs[K]
242+
K = 0;
243+
K +=1; T = 1; J = ui(T); A[K,J] = I : ds[J] → rs[K]
244+
K +=1; T = 3; J = ui(T); A[K,J] = I : ds[J] → rs[K]
245+
K +=1; T = 3; J = ui(T); A[K,J] = I : ds[J] → rs[K]
246+
K +=1; T = 5; J = ui(T); A[K,J] = I : ds[J] → rs[K]
247+
K +=1; T = 6; J = ui(T); A[K,J] = I : ds[J] → rs[K]
248+
K +=1; T = 6; J = ui(T); A[K,J] = I : ds[J] → rs[K]
249+
K +=1; T = 2; J = ui(T); A[K,J] = I : ds[J] → rs[K]
250+
K +=1; T = 2; J = ui(T); A[K,J] = I : ds[J] → rs[K]
251+
# add dirichlet interface conditions
252+
for (T1,T2) in ((1,4), (4,5), (1,2), (3,4), (5,6))
253+
global K +=1;
254+
J = ui(T1); A[K,J] = I : ds[J] → rs[K]
255+
J = ui(T2); A[K,J] = -I : ds[J] → rs[K]
29256
end
257+
# add lr neumann
258+
K +=1;
259+
T = 1; J = uy(T); A[K,J] = I : ds[J] → rs[K]
260+
T = 4; J = uy(T); A[K,J] = -I : ds[J] → rs[K]
261+
# add up neumann
262+
K +=1;
263+
T = 4; J = ux(T); A[K,J] = I : ds[J] → rs[K]
264+
T = 5; J = ux(T); A[K,J] = -I : ds[J] → rs[K]
265+
# add diagonal neumann
266+
for (T1,T2) in ((1,2), (3,4), (5,6))
267+
global K +=1;
268+
J = ux(T1); A[K,J] = I : ds[J] → rs[K]
269+
J = uy(T1); A[K,J] = I : ds[J] → rs[K]
270+
J = ux(T2); A[K,J] = -I : ds[J] → rs[K]
271+
J = uy(T2); A[K,J] = -I : ds[J] → rs[K]
272+
end
273+
274+
for T in 1:length(d)
275+
global K +=1;
276+
J = ui(T); A[K,J] = I : ds[J] → rs[K]
277+
J = ux(T); A[K,J] = -Dx : ds[J] → rs[K]
278+
global K +=1;
279+
J = ui(T); A[K,J] = I : ds[J] → rs[K]
280+
J = uy(T); A[K,J] = -Dy : ds[J] → rs[K]
281+
global K +=1;
282+
J = ux(T); A[K,J] = Dx : ds[J] → rs[K]
283+
J = uy(T); A[K,J] = Dy : ds[J] → rs[K]
284+
end
285+
286+
L = Operator(A)
287+
288+
N = 20
289+
M = sparse(L[Block.(1:N), Block.(1:N)])
290+
291+
rhs = Fun.(Ref((x,y) -> real(exp(x+im*y))), rs[1:length(∂d)], N)
292+
293+
u_cfs = M \ pad(vcat(cfs...), size(M,1))
294+
295+
u_cfs[1:sum(1:N)]
296+
297+
Fun(rhs)
298+
299+
rangespace(L)
300+
301+
302+
J
303+
ds[7]
304+
305+
[randn(2,2) randn(2,2);
306+
randn(2,2) randn(2,2)]
307+
308+
309+
U = Fun.(Ref((x,y) -> real(exp(x+im*y))), d, sum(1:N))
310+
u_cfs = Vector{Float64}()
311+
for u in U
312+
append!(u_cfs, u.coefficients)
313+
append!(u_cfs, (Dx*u).coefficients)
314+
append!(u_cfs, (Dy*u).coefficients)
315+
end
316+
317+
M*u_cfs
318+
ds[1]
319+
rs[1]
320+
321+
322+
(I : ds[1] → rs[1])*U[1]
323+
324+
# add boundary conditions
325+
for K = 1:length(∂d)
326+
for J = 1:length(d)
327+
if domain(rs[K]) ⊆ domain(ds[J])
328+
A[K,J] = I : ds[J] → rs[K]
329+
break
330+
end
331+
end
332+
end
333+
334+
Ñ = length(∂d)
335+
for (K, J1, J2) in ((Ñ+1, 1, 4),
336+
(Ñ+2, 1, 2))
337+
A[K,J1] = I : ds[J1] → rs[K]
338+
A[K,J2] = -I : ds[J2] → rs[K]
339+
end
340+
341+
Operator(A)[1:20,1:20]
342+
343+
rs
344+
345+
346+
# add continuity
347+
for K = 1:length(ιd)
348+
30349

31350

32351
∂d[1] ⊆ d[1]

‎src/Triangle.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@ function in(p::Vec{2,Float64}, d::Triangle)
1616
0 ≤ x ≤ x + y ≤ 1
1717
end
1818

19-
issubset(a::Segment{<:Vec{2}}, b::Triangle) = error("Implement")
19+
issubset(a::Segment{<:Vec{2}}, b::Triangle) = all(in.(endpoints(a), Ref(b)))
2020

2121

2222
for op in (:-, :+)

0 commit comments

Comments
 (0)