Skip to content

Commit b7afee2

Browse files
jghoskinsjghoskins
committed
update
1 parent 0b14452 commit b7afee2

File tree

4 files changed

+486
-0
lines changed

4 files changed

+486
-0
lines changed

‎chunkie/+chnk/+helm2d_sheet/green.m‎

Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,52 @@
1+
function [val,grad,hess] = green(k,src,targ)
2+
%CHNK.HELM2D.GREEN evaluate the helmholtz green's function
3+
% for the given sources and targets
4+
5+
[~,ns] = size(src);
6+
[~,nt] = size(targ);
7+
8+
xs = repmat(src(1,:),nt,1);
9+
ys = repmat(src(2,:),nt,1);
10+
11+
xt = repmat(targ(1,:).',1,ns);
12+
yt = repmat(targ(2,:).',1,ns);
13+
14+
rx = xt-xs;
15+
ry = yt-ys;
16+
17+
rx2 = rx.*rx;
18+
ry2 = ry.*ry;
19+
20+
r2 = rx2+ry2;
21+
22+
r = sqrt(r2);
23+
24+
25+
if nargout > 0
26+
h0 = besselh(0,1,k*r)-4*besselj(0,k*r);
27+
val = 0.25*1i*h0;
28+
end
29+
30+
[m,n] = size(xs);
31+
32+
if nargout > 1
33+
h1 = besselh(1,1,k*r)-4*besselj(1,k*r);
34+
grad = zeros(m,n,2);
35+
36+
grad(:,:,1) = -1i*k*0.25*h1.*rx./r;
37+
grad(:,:,2) = -1i*k*0.25*h1.*ry./r;
38+
39+
end
40+
41+
if nargout > 2
42+
43+
hess = zeros(m,n,3);
44+
45+
r3 = r.^3;
46+
47+
h2 = 2*h1./(k*r)-h0;
48+
49+
hess(:,:,1) = 0.25*1i*k*((rx-ry).*(rx+ry).*h1./r3 - k*rx2.*h0./r2);
50+
hess(:,:,2) = 0.25*1i*k*k*rx.*ry.*h2./r2;
51+
hess(:,:,3) = 0.25*1i*k*((ry-rx).*(rx+ry).*h1./r3 - k*ry2.*h0./r2);
52+
end

‎chunkie/+chnk/+helm2d_sheet/kern.m‎

Lines changed: 329 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,329 @@
1+
function submat= kern(zk,srcinfo,targinfo,type,varargin)
2+
%CHNK.HELM2D.KERN standard Helmholtz layer potential kernels in 2D
3+
%
4+
% Syntax: submat = chnk.heml2d.kern(zk,srcinfo,targingo,type,varargin)
5+
%
6+
% Let x be targets and y be sources for these formulas, with
7+
% n_x and n_y the corresponding unit normals at those points
8+
% (if defined). Note that the normal information is obtained
9+
% by taking the perpendicular to the provided tangential deriviative
10+
% info and normalizing
11+
%
12+
% Kernels based on G(x,y) = i/4 H_0^{(1)}(zk |x-y|)
13+
%
14+
% D(x,y) = \nabla_{n_y} G(x,y)
15+
% S(x,y) = G(x,y)
16+
% S'(x,y) = \nabla_{n_x} G(x,y)
17+
% D'(x,y) = \nabla_{n_x} \nabla_{n_y} G(x,y)
18+
%
19+
% Input:
20+
% zk - complex number, Helmholtz wave number
21+
% srcinfo - description of sources in ptinfo struct format, i.e.
22+
% ptinfo.r - positions (2,:) array
23+
% ptinfo.d - first derivative in underlying
24+
% parameterization (2,:)
25+
% ptinfo.d2 - second derivative in underlying
26+
% parameterization (2,:)
27+
% targinfo - description of targets in ptinfo struct format,
28+
% if info not relevant (d/d2) it doesn't need to
29+
% be provided. sprime requires tangent info in
30+
% targinfo.d
31+
% type - string, determines kernel type
32+
% type == 'd', double layer kernel D
33+
% type == 's', single layer kernel S
34+
% type == 'sprime', normal derivative of single
35+
% layer S'
36+
% type == 'dprime', normal derivative of double layer D'
37+
% type == 'c', combined layer kernel coef(1) D + coef(2) S
38+
% type == 'stau', tangential derivative of single layer
39+
% type == 'all', returns all four layer potentials,
40+
% [coef(1,1)*D coef(1,2)*S; coef(2,1)*D' coef(2,2)*S']
41+
% type == 'c2trans' returns the combined field, and the
42+
% normal derivative of the combined field
43+
% [coef(1)*D + coef(2)*S; coef(1)*D' + coef(2)*S']
44+
% type == 'trans_rep' returns the potential corresponding
45+
% to the transmission representation
46+
% [coef(1)*D coef(2)*S]
47+
% type == 'trans_rep_prime' returns the normal derivative
48+
% corresponding to the transmission representation
49+
% [coef(1)*D' coef(2)*S']
50+
% type == 'trans_rep_grad' returns the gradient corresponding
51+
% to the transmission representation
52+
% [coef(1)*d_x D coef(2)*d_x S;
53+
% coef(1)*d_y D coef(2)*d_y S]
54+
%
55+
% append '_diff' to type to compute difference stably
56+
%
57+
% varargin{1} - coef: length 2 array in the combined layer
58+
% formula, 2x2 matrix for all kernels
59+
% otherwise does nothing
60+
%
61+
% Output:
62+
% submat - the evaluation of the selected kernel for the
63+
% provided sources and targets. the number of
64+
% rows equals the number of targets and the
65+
% number of columns equals the number of sources
66+
%
67+
% see also CHNK.HELM2D.GREEN
68+
69+
src = srcinfo.r(:,:);
70+
targ = targinfo.r(:,:);
71+
72+
[~,ns] = size(src);
73+
[~,nt] = size(targ);
74+
75+
% double layer
76+
if strcmpi(type,'d')
77+
srcnorm = srcinfo.n(:,:);
78+
[~,grad] = chnk.helm2d_sheet.green(zk,src,targ);
79+
nx = repmat(srcnorm(1,:),nt,1);
80+
ny = repmat(srcnorm(2,:),nt,1);
81+
submat = -(grad(:,:,1).*nx + grad(:,:,2).*ny);
82+
end
83+
84+
% normal derivative of single layer
85+
if strcmpi(type,'sprime')
86+
targnorm = targinfo.n(:,:);
87+
[~,grad] = chnk.helm2d_sheet.green(zk,src,targ);
88+
nx = repmat((targnorm(1,:)).',1,ns);
89+
ny = repmat((targnorm(2,:)).',1,ns);
90+
91+
submat = (grad(:,:,1).*nx + grad(:,:,2).*ny);
92+
end
93+
94+
% Tangential derivative of single layer
95+
if strcmpi(type,'stau')
96+
targtan = targinfo.d(:,:);
97+
[~,grad] = chnk.helm2d_sheet.green(zk,src,targ);
98+
dx = repmat((targtan(1,:)).',1,ns);
99+
dy = repmat((targtan(2,:)).',1,ns);
100+
ds = sqrt(dx.*dx+dy.*dy);
101+
submat = (grad(:,:,1).*dx + grad(:,:,2).*dy)./ds;
102+
end
103+
104+
% single layer
105+
if strcmpi(type,'s')
106+
submat = chnk.helm2d_sheet.green(zk,src,targ);
107+
end
108+
109+
% normal derivative of double layer
110+
if strcmpi(type,'dprime')
111+
targnorm = targinfo.n(:,:);
112+
srcnorm = srcinfo.n(:,:);
113+
[~,~,hess] = chnk.helm2d_sheet.green(zk,src,targ);
114+
nxsrc = repmat(srcnorm(1,:),nt,1);
115+
nysrc = repmat(srcnorm(2,:),nt,1);
116+
nxtarg = repmat((targnorm(1,:)).',1,ns);
117+
nytarg = repmat((targnorm(2,:)).',1,ns);
118+
submat = -(hess(:,:,1).*nxsrc.*nxtarg + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)...
119+
+ hess(:,:,3).*nysrc.*nytarg);
120+
end
121+
122+
123+
% Combined field
124+
if strcmpi(type,'c')
125+
srcnorm = srcinfo.n(:,:);
126+
coef = ones(2,1);
127+
if(nargin == 5); coef = varargin{1}; end
128+
[submats,grad] = chnk.helm2d_sheet.green(zk,src,targ);
129+
nx = repmat(srcnorm(1,:),nt,1);
130+
ny = repmat(srcnorm(2,:),nt,1);
131+
submatd = -(grad(:,:,1).*nx + grad(:,:,2).*ny);
132+
submat = coef(1)*submatd + coef(2)*submats;
133+
end
134+
135+
% normal derivative of combined field
136+
if strcmpi(type,'cprime')
137+
coef = ones(2,1);
138+
if(nargin == 5); coef = varargin{1}; end
139+
targnorm = targinfo.n(:,:);
140+
srcnorm = srcinfo.n(:,:);
141+
142+
143+
144+
% Get gradient and hessian info
145+
[~,grad,hess] = chnk.helm2d_sheet.green(zk,src,targ);
146+
147+
nxsrc = repmat(srcnorm(1,:),nt,1);
148+
nysrc = repmat(srcnorm(2,:),nt,1);
149+
nxtarg = repmat((targnorm(1,:)).',1,ns);
150+
nytarg = repmat((targnorm(2,:)).',1,ns);
151+
152+
% D'
153+
submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ...
154+
+ hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)...
155+
+ hess(:,:,3).*nysrc.*nytarg);
156+
% S'
157+
submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg);
158+
159+
submat = coef(1)*submatdp + coef(2)*submatsp;
160+
end
161+
162+
163+
% Dirichlet and neumann data corresponding to combined field
164+
if strcmpi(type,'c2trans')
165+
coef = ones(2,1);
166+
if(nargin == 5); coef = varargin{1}; end
167+
targnorm = targinfo.n(:,:);
168+
srcnorm = srcinfo.n(:,:);
169+
170+
% Get gradient and hessian info
171+
[submats,grad,hess] = chnk.helm2d_sheet.green(zk,src,targ);
172+
173+
nxsrc = repmat(srcnorm(1,:),nt,1);
174+
nysrc = repmat(srcnorm(2,:),nt,1);
175+
nxtarg = repmat((targnorm(1,:)).',1,ns);
176+
nytarg = repmat((targnorm(2,:)).',1,ns);
177+
178+
% D'
179+
submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ...
180+
+ hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)...
181+
+ hess(:,:,3).*nysrc.*nytarg);
182+
% S'
183+
submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg);
184+
% D
185+
submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc);
186+
187+
submat = zeros(2*nt,ns);
188+
189+
submat(1:2:2*nt,:) = coef(1)*submatd + coef(2)*submats;
190+
submat(2:2:2*nt,:) = coef(1)*submatdp + coef(2)*submatsp;
191+
192+
end
193+
194+
195+
% all kernels, [c11 D, c12 S; c21 D', c22 S']
196+
if strcmpi(type,'all')
197+
198+
targnorm = targinfo.n(:,:);
199+
srcnorm = srcinfo.n(:,:);
200+
cc = varargin{1};
201+
202+
submat = zeros(2*nt,2*ns);
203+
204+
% Get gradient and hessian info
205+
[submats,grad,hess] = chnk.helm2d.green(zk,src,targ);
206+
207+
nxsrc = repmat(srcnorm(1,:),nt,1);
208+
nysrc = repmat(srcnorm(2,:),nt,1);
209+
nxtarg = repmat((targnorm(1,:)).',1,ns);
210+
nytarg = repmat((targnorm(2,:)).',1,ns);
211+
212+
submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ...
213+
+ hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)...
214+
+ hess(:,:,3).*nysrc.*nytarg);
215+
% S'
216+
submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg);
217+
% D
218+
submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc);
219+
220+
221+
submat(1:2:2*nt,1:2:2*ns) = submatd*cc(1,1);
222+
submat(1:2:2*nt,2:2:2*ns) = submats*cc(1,2);
223+
submat(2:2:2*nt,1:2:2*ns) = submatdp*cc(2,1);
224+
submat(2:2:2*nt,2:2:2*ns) = submatsp*cc(2,2);
225+
end
226+
227+
% all kernels, [c11 D, c12 S; c21 D', c22 S'] (difference)
228+
if strcmpi(type,'all_diff')
229+
230+
targnorm = targinfo.n(:,:);
231+
srcnorm = srcinfo.n(:,:);
232+
cc = varargin{1};
233+
234+
submat = zeros(2*nt,2*ns);
235+
236+
% Get gradient and hessian info
237+
[submats,grad,hess] = chnk.helm2d_sheet.helmdiffgreen(zk,src,targ);
238+
239+
nxsrc = repmat(srcnorm(1,:),nt,1);
240+
nysrc = repmat(srcnorm(2,:),nt,1);
241+
nxtarg = repmat((targnorm(1,:)).',1,ns);
242+
nytarg = repmat((targnorm(2,:)).',1,ns);
243+
244+
submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ...
245+
+ hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)...
246+
+ hess(:,:,3).*nysrc.*nytarg);
247+
% S'
248+
submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg);
249+
% D
250+
submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc);
251+
252+
253+
submat(1:2:2*nt,1:2:2*ns) = submatd*cc(1,1);
254+
submat(1:2:2*nt,2:2:2*ns) = submats*cc(1,2);
255+
submat(2:2:2*nt,1:2:2*ns) = submatdp*cc(2,1);
256+
submat(2:2:2*nt,2:2:2*ns) = submatsp*cc(2,2);
257+
end
258+
259+
% Dirichlet data/potential correpsonding to transmission rep
260+
if strcmpi(type,'trans_rep')
261+
262+
coef = ones(2,1);
263+
if(nargin == 5); coef = varargin{1}; end;
264+
srcnorm = srcinfo.n(:,:);
265+
[submats,grad] = chnk.helm2d_sheet.green(zk,src,targ);
266+
nx = repmat(srcnorm(1,:),nt,1);
267+
ny = repmat(srcnorm(2,:),nt,1);
268+
submatd = -(grad(:,:,1).*nx + grad(:,:,2).*ny);
269+
270+
submat = zeros(1*nt,2*ns);
271+
submat(1:1:1*nt,1:2:2*ns) = coef(1)*submatd;
272+
submat(1:1:1*nt,2:2:2*ns) = coef(2)*submats;
273+
end
274+
275+
276+
% Neumann data corresponding to transmission rep
277+
if strcmpi(type,'trans_rep_prime')
278+
coef = ones(2,1);
279+
if(nargin == 5); coef = varargin{1}; end;
280+
targnorm = targinfo.n(:,:);
281+
srcnorm = srcinfo.n(:,:);
282+
283+
submat = zeros(nt,ns);
284+
285+
% Get gradient and hessian info
286+
[submats,grad,hess] = chnk.helm2d_sheet.green(zk,src,targ);
287+
288+
nxsrc = repmat(srcnorm(1,:),nt,1);
289+
nysrc = repmat(srcnorm(2,:),nt,1);
290+
nxtarg = repmat((targnorm(1,:)).',1,ns);
291+
nytarg = repmat((targnorm(2,:)).',1,ns);
292+
293+
% D'
294+
submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ...
295+
+ hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)...
296+
+ hess(:,:,3).*nysrc.*nytarg);
297+
% S'
298+
submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg);
299+
submat = zeros(nt,2*ns);
300+
submat(1:1:1*nt,1:2:2*ns) = coef(1)*submatdp;
301+
submat(1:1:1*nt,2:2:2*ns) = coef(2)*submatsp;
302+
end
303+
304+
305+
% Gradient correpsonding to transmission rep
306+
if strcmpi(type,'trans_rep_grad')
307+
coef = ones(2,1);
308+
if(nargin == 5); coef = varargin{1}; end;
309+
310+
srcnorm = srcinfo.n(:,:);
311+
312+
submat = zeros(nt,ns,6);
313+
% S
314+
[submats,grad,hess] = chnk.helm2d_sheet.green(zk,src,targ);
315+
316+
nxsrc = repmat(srcnorm(1,:),nt,1);
317+
nysrc = repmat(srcnorm(2,:),nt,1);
318+
% D
319+
submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc);
320+
321+
submat = zeros(2*nt,2*ns);
322+
323+
submat(1:2:2*nt,1:2:2*ns) = -coef(1)*(hess(:,:,1).*nxsrc + hess(:,:,2).*nysrc);
324+
submat(1:2:2*nt,2:2:2*ns) = coef(2)*grad(:,:,1);
325+
326+
submat(2:2:2*nt,1:2:2*ns) = -coef(1)*(hess(:,:,2).*nxsrc + hess(:,:,3).*nysrc);
327+
submat(2:2:2*nt,2:2:2*ns) = coef(2)*grad(:,:,2);
328+
end
329+

0 commit comments

Comments
 (0)