-
Notifications
You must be signed in to change notification settings - Fork 0
/
testLaplaceNeu.m
82 lines (69 loc) · 2.63 KB
/
testLaplaceNeu.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
% test a Laplace Neumann BVP
%
% Hai 05/18/21 updated convergence plot 12/17/21, maybe add more BVP types later
setup()
v = 1;
s.p = 16;
shape = 'star'; % or 'star', or 'star2'. test convergence on different shapes
% family of shapes...
if strcmp(shape,'ellipse')
lam = 1; ratio = 2; s.Z = @(t) lam*(sin(t) - ratio*1i*cos(t));
Np = 4:2:24;
elseif strcmp(shape,'star')
a = .3; w = 5; R = @(t) 2*(1 + a*cos(w*t-pi/2)); s.Z = @(t) R(t).*exp(1i*t); s.Z = @(t) s.Z(t-pi/2);
Np = 20:2:40;
elseif strcmp(shape,'star2')
s.Z = @(t) ((cos(t).^2+9*sin(t).^2).^(1/2)+cos(4*t).^2).*sin(t)...
-1/2*1i*((cos(t).^2+9*sin(t).^2).^(1/2)+cos(4*t).^2).*cos(t);
np = 48;
Np = 26:2:50;
end
tpan = linspace(0,pi,24+1)'; s.tpan = tpan;
s = quadr(s,[],'p','G');
% target
nx = 150; gx = ((1:nx)/nx*2)*4; ny = 150; gy = ((1:ny)/ny*2-1)*6; % set up plotting grid
[xx, yy] = meshgrid(gx,gy); zz = (xx+1i*yy);
t = [];
[IN, ON] = inpolygon(real(zz),imag(zz),real(s.x),imag(s.x));
ii = ~IN;
t.x = zz(ii(:)); % eval pts only on one side
t.nx = ones(size(t.x))/2 + 1i*ones(size(t.x));
% exact soln
alpha = linspace(-pi/3,pi/5,5);
y_source = [];
y_source.x = 3/5*(1*exp(1i*alpha).'+0.1+1i*0.3); y_source.ws = ones(size(y_source.x));
pt_force = [1;1;-1;1;1/2];
fhom = nan*xx; % exact soln
Ahom = LapAxiMat(t,y_source);
fhom(ii(:)) = Ahom*pt_force;
% h refinement, increase num of panels
err = NaN(size(Np));
for k=1:numel(Np)
% boundary discretization
np = Np(k);
tpan = linspace(0,pi,np+1)'; s.tpan = tpan;
s = quadr(s,[],'p','G');
% use grad SLP to solve a Neumann problem
A = LapGradAxiSpecialMat(s,s); % A = (-1/2*I+S')
rhs = LapGradAxiMat(s,y_source)*pt_force;
tau = A\rhs;
% use SLP to evaluate soln
u = nan*xx; % soln
u(ii(:)) = LapAxiSpecialMat(t,s)*tau;
% plot error
figure(1),clf,imagesc(gx,gy,log10(abs(u-fhom))),
colorbar, hold on, fill(real([s.x;s.x(1)]),imag([s.x;s.x(1)]),'w')
plot(s.x,'-k'), plot([s.xlo;s.xhi(end)],'.r','MarkerSize',12); caxis([-15 0]), colormap('jet')
axis equal tight
err(k) = max(abs(u(:)-fhom(:)));
end
figure(1),clf,imagesc(gx,gy,log10(abs(u-fhom))),
colorbar, hold on, fill(real([s.x;s.x(1)]),imag([s.x;s.x(1)]),'w')
plot(s.x,'-k'), plot([s.xlo;s.xhi(end)],'.r','MarkerSize',12); colormap('jet')
title('log10 err in |u|'), axis equal tight
figure(2),clf, semilogy(s.p*Np,err,'o-k'), title('BVP conv')
xlabel('$n_p$','interpreter','latex')
% fit_order = [log(Np'),ones(numel(Np),1)]\log(err');
% errfit = @(x) err(1)/Np(1)^(-2*s.p+1)*x.^(-2*s.p+1);
% p4=loglog(Np(:),err(:),'p-k','MarkerFaceColor','k');hold on, loglog(Np(:),errfit(Np(:)),'--k')
keyboard