diff --git a/chunkie/+chnk/+flex2d/fmm.m b/chunkie/+chnk/+flex2d/fmm.m index b62111a..4d50789 100644 --- a/chunkie/+chnk/+flex2d/fmm.m +++ b/chunkie/+chnk/+flex2d/fmm.m @@ -1,121 +1,121 @@ -function varargout = fmm(eps, zk, srcinfo, targinfo, type, sigma, varargin) -%CHNK.HELM2D.FMM Fast multipole methods for evaluating Helmholtz layer -%potentials, their gradients, and Hessians. -% -% Syntax: -% pot = chnk.helm2d.fmm(eps, zk, srcinfo, targinfo, type, sigma, varargin) -% [pot, grad] = chnk.helm2d.fmm(eps, zk, srcinfo, targinfo, type, sigma, varargin) -% [pot, grad, hess] = chnk.helm2d.fmm(eps, zk, srcinfo, targinfo, type, sigma, varargin) -% -% Let x be targets and y be sources for these formulas, with -% n_x and n_y the corresponding unit normals at those points -% (if defined). -% -% Kernels based on G(x,y) = i/4 H_0^{(1)}(zk |x-y|) -% -% D(x,y) = \nabla_{n_y} G(x,y) -% S(x,y) = G(x,y) -% -% Note: the density should be scaled by the weights -% -% Input: -% eps - precision requested -% zk - complex number, Helmholtz wave number -% srcinfo - description of sources in ptinfo struct format, i.e. -% ptinfo.r - positions (2,:) array -% ptinfo.d - first derivative in underlying -% parameterization (2,:) -% ptinfo.d2 - second derivative in underlying -% parameterization (2,:) -% ptinfo.n - normal (2,:) array -% -% targinfo - description of targets in ptinfo struct format, -% type - string, determines kernel type -% type == 'd', double layer kernel D -% type == 's', single layer kernel S -% type == 'c', combined layer kernel coefs(1)*D + coefs(2)*S -% type = 'sprime' normal derivative of single layer, -% note that targinfo must contain targ.n -% type = 'dprime' normal derivative of double layer, -% note that targinfo must contain targ.n -% sigma - density -% varargin{1} - coef in the combined layer formula, otherwise -% does nothing -% -% Optional Output: -% pot - potential/neumann data corresponding to the kernel at the target locations -% grad - gradient at target locations -% hess - hessian at target locations - -if ( nargout == 0 ) - warning('CHUNKIE:helm2d:fmm:empty', ... - 'Nothing to compute in HELM2D.FMM. Returning empty array.'); - return -end - -srcuse = []; -srcuse.sources = srcinfo.r(1:2,:); -switch lower(type) - case {'s', 'sprime'} - srcuse.charges = sigma(:).'; - case {'d', 'dprime'} - srcuse.dipstr = sigma(:).'; - srcuse.dipvec = srcinfo.n(1:2,:); - case 'c' - coefs = varargin{1}; - srcuse.charges = coefs(2)*sigma(:).'; - srcuse.dipstr = coefs(1)*sigma(:).'; - srcuse.dipvec = srcinfo.n(1:2,:); -end - -if ( isstruct(targinfo) ) - targuse = targinfo.r(:,:); -else - targuse = targinfo; -end - -pg = 0; -pgt = min(nargout, 2); -U = hfmm2d(eps, zk, srcuse, pg, targuse, pgt); - -% Assign potentials -if ( nargout > 0 ) - switch lower(type) - case {'s', 'd', 'c'} - varargout{1} = U.pottarg.'; - case {'sprime', 'dprime'} - if ( ~isfield(targinfo, 'n') ) - error('CHUNKIE:helm2d:fmm:normals', ... - 'Targets require normal info when evaluating Helmholtz kernel ''%s''.', type); - end - varargout{1} = ( U.gradtarg(1,:).*targinfo.n(1,:) + ... - U.gradtarg(2,:).*targinfo.n(2,:) ).'; - otherwise - error('CHUNKIE:lap2d:fmm:pot', ... - 'Potentials not supported for Helmholtz kernel ''%s''.', type); - end -end - -% Assign gradients -if ( nargout > 1 ) - switch lower(type) - case {'s', 'd', 'c'} - varargout{2} = U.gradtarg; - otherwise - error('CHUNKIE:helm2d:fmm:grad', ... - 'Gradients not supported for Helmholtz kernel ''%s''.', type); - end -end - -% Assign Hessians -if ( nargout > 2 ) - switch lower(type) - case {'s', 'd', 'c'} - varargout{3} = U.hesstarg; - otherwise - error('CHUNKIE:helm2d:fmm:hess', ... - 'Hessians not supported for Helmholtz kernel ''%s''.', type); - end -end - -end +function varargout = fmm(eps, zk, srcinfo, targinfo, type, sigma, varargin) +%CHNK.HELM2D.FMM Fast multipole methods for evaluating Helmholtz layer +%potentials, their gradients, and Hessians. +% +% Syntax: +% pot = chnk.helm2d.fmm(eps, zk, srcinfo, targinfo, type, sigma, varargin) +% [pot, grad] = chnk.helm2d.fmm(eps, zk, srcinfo, targinfo, type, sigma, varargin) +% [pot, grad, hess] = chnk.helm2d.fmm(eps, zk, srcinfo, targinfo, type, sigma, varargin) +% +% Let x be targets and y be sources for these formulas, with +% n_x and n_y the corresponding unit normals at those points +% (if defined). +% +% Kernels based on G(x,y) = i/4 H_0^{(1)}(zk |x-y|) +% +% D(x,y) = \nabla_{n_y} G(x,y) +% S(x,y) = G(x,y) +% +% Note: the density should be scaled by the weights +% +% Input: +% eps - precision requested +% zk - complex number, Helmholtz wave number +% srcinfo - description of sources in ptinfo struct format, i.e. +% ptinfo.r - positions (2,:) array +% ptinfo.d - first derivative in underlying +% parameterization (2,:) +% ptinfo.d2 - second derivative in underlying +% parameterization (2,:) +% ptinfo.n - normal (2,:) array +% +% targinfo - description of targets in ptinfo struct format, +% type - string, determines kernel type +% type == 'd', double layer kernel D +% type == 's', single layer kernel S +% type == 'c', combined layer kernel coefs(1)*D + coefs(2)*S +% type = 'sprime' normal derivative of single layer, +% note that targinfo must contain targ.n +% type = 'dprime' normal derivative of double layer, +% note that targinfo must contain targ.n +% sigma - density +% varargin{1} - coef in the combined layer formula, otherwise +% does nothing +% +% Optional Output: +% pot - potential/neumann data corresponding to the kernel at the target locations +% grad - gradient at target locations +% hess - hessian at target locations + +if ( nargout == 0 ) + warning('CHUNKIE:helm2d:fmm:empty', ... + 'Nothing to compute in HELM2D.FMM. Returning empty array.'); + return +end + +srcuse = []; +srcuse.sources = srcinfo.r(1:2,:); +switch lower(type) + case {'s', 'sprime'} + srcuse.charges = sigma(:).'; + case {'d', 'dprime'} + srcuse.dipstr = sigma(:).'; + srcuse.dipvec = srcinfo.n(1:2,:); + case 'c' + coefs = varargin{1}; + srcuse.charges = coefs(2)*sigma(:).'; + srcuse.dipstr = coefs(1)*sigma(:).'; + srcuse.dipvec = srcinfo.n(1:2,:); +end + +if ( isstruct(targinfo) ) + targuse = targinfo.r(:,:); +else + targuse = targinfo; +end + +pg = 0; +pgt = min(nargout, 2); +U = hfmm2d(eps, zk, srcuse, pg, targuse, pgt); + +% Assign potentials +if ( nargout > 0 ) + switch lower(type) + case {'s', 'd', 'c'} + varargout{1} = U.pottarg.'; + case {'sprime', 'dprime'} + if ( ~isfield(targinfo, 'n') ) + error('CHUNKIE:helm2d:fmm:normals', ... + 'Targets require normal info when evaluating Helmholtz kernel ''%s''.', type); + end + varargout{1} = ( U.gradtarg(1,:).*targinfo.n(1,:) + ... + U.gradtarg(2,:).*targinfo.n(2,:) ).'; + otherwise + error('CHUNKIE:lap2d:fmm:pot', ... + 'Potentials not supported for Helmholtz kernel ''%s''.', type); + end +end + +% Assign gradients +if ( nargout > 1 ) + switch lower(type) + case {'s', 'd', 'c'} + varargout{2} = U.gradtarg; + otherwise + error('CHUNKIE:helm2d:fmm:grad', ... + 'Gradients not supported for Helmholtz kernel ''%s''.', type); + end +end + +% Assign Hessians +if ( nargout > 2 ) + switch lower(type) + case {'s', 'd', 'c'} + varargout{3} = U.hesstarg; + otherwise + error('CHUNKIE:helm2d:fmm:hess', ... + 'Hessians not supported for Helmholtz kernel ''%s''.', type); + end +end + +end diff --git a/chunkie/+chnk/+flex2d/green.m b/chunkie/+chnk/+flex2d/green.m index 04fdb2f..e65ce0b 100644 --- a/chunkie/+chnk/+flex2d/green.m +++ b/chunkie/+chnk/+flex2d/green.m @@ -1,52 +1,52 @@ -function [val,grad,hess] = green(k,src,targ) -%CHNK.HELM2D.GREEN evaluate the helmholtz green's function -% for the given sources and targets - -[~,ns] = size(src); -[~,nt] = size(targ); - -xs = repmat(src(1,:),nt,1); -ys = repmat(src(2,:),nt,1); - -xt = repmat(targ(1,:).',1,ns); -yt = repmat(targ(2,:).',1,ns); - -rx = xt-xs; -ry = yt-ys; - -rx2 = rx.*rx; -ry2 = ry.*ry; - -r2 = rx2+ry2; - -r = sqrt(r2); - - -if nargout > 0 - h0 = besselh(0,1,k*r); - val = 0.25*1i*h0; -end - -[m,n] = size(xs); - -if nargout > 1 - h1 = besselh(1,1,k*r); - grad = zeros(m,n,2); - - grad(:,:,1) = -1i*k*0.25*h1.*rx./r; - grad(:,:,2) = -1i*k*0.25*h1.*ry./r; - -end - -if nargout > 2 - - hess = zeros(m,n,3); - - r3 = r.^3; - - h2 = 2*h1./(k*r)-h0; - - hess(:,:,1) = 0.25*1i*k*((rx-ry).*(rx+ry).*h1./r3 - k*rx2.*h0./r2); - hess(:,:,2) = 0.25*1i*k*k*rx.*ry.*h2./r2; - hess(:,:,3) = 0.25*1i*k*((ry-rx).*(rx+ry).*h1./r3 - k*ry2.*h0./r2); -end +function [val,grad,hess] = green(k,src,targ) +%CHNK.HELM2D.GREEN evaluate the helmholtz green's function +% for the given sources and targets + +[~,ns] = size(src); +[~,nt] = size(targ); + +xs = repmat(src(1,:),nt,1); +ys = repmat(src(2,:),nt,1); + +xt = repmat(targ(1,:).',1,ns); +yt = repmat(targ(2,:).',1,ns); + +rx = xt-xs; +ry = yt-ys; + +rx2 = rx.*rx; +ry2 = ry.*ry; + +r2 = rx2+ry2; + +r = sqrt(r2); + + +if nargout > 0 + h0 = besselh(0,1,k*r); + val = 0.25*1i*h0; +end + +[m,n] = size(xs); + +if nargout > 1 + h1 = besselh(1,1,k*r); + grad = zeros(m,n,2); + + grad(:,:,1) = -1i*k*0.25*h1.*rx./r; + grad(:,:,2) = -1i*k*0.25*h1.*ry./r; + +end + +if nargout > 2 + + hess = zeros(m,n,3); + + r3 = r.^3; + + h2 = 2*h1./(k*r)-h0; + + hess(:,:,1) = 0.25*1i*k*((rx-ry).*(rx+ry).*h1./r3 - k*rx2.*h0./r2); + hess(:,:,2) = 0.25*1i*k*k*rx.*ry.*h2./r2; + hess(:,:,3) = 0.25*1i*k*((ry-rx).*(rx+ry).*h1./r3 - k*ry2.*h0./r2); +end diff --git a/chunkie/+chnk/+flex2d/helmdiffgreen.m b/chunkie/+chnk/+flex2d/helmdiffgreen.m index cf9489b..1b2b2e6 100644 --- a/chunkie/+chnk/+flex2d/helmdiffgreen.m +++ b/chunkie/+chnk/+flex2d/helmdiffgreen.m @@ -1,10 +1,17 @@ -function [val,grad,hess,der3,der4] = helmdiffgreen(k,src,targ) +function [val,grad,hess,der3,der4] = helmdiffgreen(k,src,targ,ifr2logr) %HELMDIFFGREEN evaluate the difference of the % Helmholtz Green function and the Laplace Green function % for the given sources and targets, i.e. % % G(x,y) = i/4 H_0^(1)(k|x-y|) + 1/(2 pi) log(|x-y|) % +% or the difference of the Helmholtz and Laplace Green funcions +% and k^2 r^2 log r/ 8 pi (a constant times the biharmonic Green function) +% i.e. +% +% G(x,y) = i/4 H_0^(1)(k|x-y|) + 1/(2 pi) log(|x-y|) + ... +% - k^2/(8*pi) |x-y|^2 log(|x-y|) +% % where H_0^(1) is the principal branch of the Hankel function % of the first kind. This routine avoids numerical cancellation % when |k||x-y| is small. @@ -18,6 +25,26 @@ % G_{x1x1x1x2}, G_{x1x1x2x2}, G_{x1x2x2x2}, G_{x2x2x2x2} % % derivatives are on the *target* variables +% +% input: +% +% src - (2,ns) array of source locations +% targ - (2,nt) array of target locations +% k - wave number, as above +% +% optional input: +% +% ifr2logr - boolean, default: false. If true, also subtract off the +% k^2/(8pi) r^2 log r kernel + +if nargin < 4 + ifr2logr = false; +end + +r2logrfac = 1; +if ifr2logr + r2logrfac = 0; +end [~,ns] = size(src); [~,nt] = size(targ); @@ -48,7 +75,7 @@ % get value and r derivatives -[g0,g1,g2,g3,g4] = diff_h0log_and_rders(k,r); +[g0,g1,g21,g3,g4] = diff_h0log_and_rders(k,r,r2logrfac); % evaluate potential and derivatives @@ -60,49 +87,57 @@ grad(:,:,2) = dy.*g1.*rm1; end if nargout > 2 - hess(:,:,1) = dx2.*g2.*rm2+g1.*(1.*rm1-dx2.*rm3); - hess(:,:,2) = dx.*dy.*(g2.*rm2-g1.*rm3); - hess(:,:,3) = dy2.*g2.*rm2+g1.*(1.*rm1-dy2.*rm3); + hess(:,:,1) = dx2.*g21.*rm2+g1.*rm1; + hess(:,:,2) = dx.*dy.*g21.*rm2; + hess(:,:,3) = dy2.*g21.*rm2+g1.*rm1; end if nargout > 3 - der3(:,:,1) = (dx3.*g3+3*dy2.*dx.*(g2.*rm1-g1.*rm2)).*rm3; - der3(:,:,2) = dx2.*dy.*(g3.*rm3-3*(g2.*rm4-g1.*rm5)) + ... - dy.*(g2.*rm2-g1.*rm3); - der3(:,:,3) = dx.*dy2.*(g3.*rm3-3*(g2.*rm4-g1.*rm5)) + ... - dx.*(g2.*rm2-g1.*rm3); - der3(:,:,4) = (dy3.*g3+3*dx2.*dy.*(g2.*rm1-g1.*rm2)).*rm3; + der3(:,:,1) = (dx3.*g3+3*dy2.*dx.*g21.*rm1).*rm3; + der3(:,:,2) = dx2.*dy.*(g3.*rm3-3*g21.*rm4) + ... + dy.*g21.*rm2; + der3(:,:,3) = dx.*dy2.*(g3.*rm3-3*g21.*rm4) + ... + dx.*g21.*rm2; + der3(:,:,4) = (dy3.*g3+3*dx2.*dy.*g21.*rm1).*rm3; end if nargout > 4 - der4(:,:,1) = (dx4.*(g4-6*g3.*rm1+15*(g2.*rm2-g1.*rm3))).*rm4 + ... - (6*dx2.*(g3-3*(g2.*rm1-g1.*rm2))).*rm3 + ... - (3*(g2-g1.*rm1)).*rm2; - der4(:,:,2) = (dx3.*dy.*(g4-6*g3.*rm1+15*(g2.*rm2-g1.*rm3))).*rm4 + ... - (3*dx.*dy.*(g3-3*(g2.*rm1-g1.*rm2))).*rm3; - der4(:,:,3) = dx2.*dy2.*(g4-6*g3.*rm1+15*g2.*rm2-15*g1.*rm3).*rm4 + ... - g3.*rm1 - 2*g2.*rm2 + 2*g1.*rm3; - der4(:,:,4) = dx.*dy3.*(g4-6*g3.*rm1+15*(g2.*rm2-g1.*rm3)).*rm4 + ... - 3*dx.*dy.*(g3-3*(g2.*rm1-g1.*rm2)).*rm3; - der4(:,:,5) = dy4.*(g4-6*g3.*rm1+15*(g2.*rm2-g1.*rm3)).*rm4 + ... - 6*dy2.*(g3-3*(g2.*rm1-g1.*rm2)).*rm3 + ... - 3*(g2-g1.*rm1).*rm2; + der4(:,:,1) = (dx4.*(g4-6*g3.*rm1+15*g21.*rm2)).*rm4 + ... + (6*dx2.*(g3-3*g21.*rm1)).*rm3 + ... + 3*g21.*rm2; + der4(:,:,2) = (dx3.*dy.*(g4-6*g3.*rm1+15*g21.*rm2)).*rm4 + ... + (3*dx.*dy.*(g3-3*g21.*rm1)).*rm3; + der4(:,:,3) = dx2.*dy2.*(g4-6*g3.*rm1+15*g21.*rm2).*rm4 + ... + g3.*rm1 - 2*g21.*rm2; + der4(:,:,4) = dx.*dy3.*(g4-6*g3.*rm1+15*g21.*rm2).*rm4 + ... + 3*dx.*dy.*(g3-3*g21.*rm1).*rm3; + der4(:,:,5) = dy4.*(g4-6*g3.*rm1+15*g21.*rm2).*rm4 + ... + 6*dy2.*(g3-3*g21.*rm1).*rm3 + ... + 3*g21.*rm2; end end -function [g0,g1,g2,g3,g4] = diff_h0log_and_rders(k,r) +function [g0,g1,g21,g3,g4] = diff_h0log_and_rders(k,r,r2logrfac) +% g0 = g +% g1 = g' +% g21 = g'' - g'/r +% +% maybe later: +% g321 = g''' - 3*g''/r + 3g'/r^2 +% g4321 = g'''' - 6*g'''/r + 15*g''/r^2 - 15*g'/r^3 g0 = zeros(size(r)); g1 = zeros(size(r)); -g2 = zeros(size(r)); g3 = zeros(size(r)); g4 = zeros(size(r)); +g21 = zeros(size(r)); io4 = 1i*0.25; o2p = 1/(2*pi); isus = abs(k)*r < 1; %isus = false(size(r)); +%isus = true(size(r)); % straightforward formulas for sufficiently large @@ -118,11 +153,14 @@ rm3 = rm1.*rm2; rm4 = rm1.*rm3; -g0(~isus) = io4*h0 + o2p*log(rnot); -g1(~isus) = -k*io4*h1 + o2p*rm1; -g2(~isus) = -k*k*io4*h0 + k*io4*h1.*rm1 - o2p*rm2; -g3(~isus) = k*k*io4*h0.*rm1 + io4*k*(k*k-2*rm2).*h1 + 2*o2p*rm3; -g4(~isus) = k*io4*(3*rm2-k*k).*(2*h1.*rm1-k*h0) - 6*o2p*rm4; +r2fac = (1-r2logrfac)*k*k*0.25*o2p; +logr = log(rnot); +g0(~isus) = io4*h0 + o2p*logr - r2fac*rnot.*rnot.*logr; +g1(~isus) = -k*io4*h1 + o2p*rm1 - r2fac*(rnot+2*rnot.*logr); +g21(~isus) = -k*k*io4*h0 + k*io4*h1.*rm1 - o2p*rm2 - r2fac*(3+2*logr) - ... + g1(~isus).*rm1; +g3(~isus) = k*k*io4*h0.*rm1 + io4*k*(k*k-2*rm2).*h1 + 2*o2p*rm3 - r2fac*2*rm1; +g4(~isus) = k*io4*(3*rm2-k*k).*(2*h1.*rm1-k*h0) - 6*o2p*rm4 + r2fac*2*rm2; % manually cancel when small @@ -133,11 +171,12 @@ rm4 = rm1.*rm3; gam = 0.57721566490153286060651209; -nterms = 14; +nterms = 20; const1 = (io4+(log(2)-gam-log(k))*o2p); % relevant parts of hankel function represented as power series [cf1,cf2] = chnk.flex2d.besseldiff_etc_pscoefs(nterms); +cf1(1) = cf1(1)*r2logrfac; kpow = (k*k).^(1:nterms); cf1 = cf1(:).*kpow(:); cf2 = cf2(:).*kpow(:); @@ -148,13 +187,15 @@ % differentiate power series to get derivatives fac = 2*(1:nterms); +d21 = fac(:).*(fac(:)-1)-fac(:); +fd21 = chnk.flex2d.even_pseval(cf2(:).*d21,rsus).*rm2; cf1 = cf1.*fac(:); cf2 = cf2.*fac(:); j0m1d1 = chnk.flex2d.even_pseval(cf1,rsus).*rm1; fd1 = chnk.flex2d.even_pseval(cf2,rsus).*rm1; cf1 = cf1.*(fac(:)-1); cf2 = cf2.*(fac(:)-1); j0m1d2 = chnk.flex2d.even_pseval(cf1,rsus).*rm2; - fd2 = chnk.flex2d.even_pseval(cf2,rsus).*rm2; + cf1 = cf1(:).*(fac(:)-2); cf1 = cf1(2:end); cf2 = cf2(:).*(fac(:)-2); cf2 = cf2(2:end); j0m1d3 = chnk.flex2d.even_pseval(cf1,rsus).*rm1; @@ -166,9 +207,10 @@ fd4 = chnk.flex2d.even_pseval(cf2,rsus).*rm2; % combine to get derivative of i/4 H + log/(2*pi) -g0(isus) = const1*(j0m1+1) - o2p*(f + logr.*j0m1); -g1(isus) = const1*j0m1d1 - o2p*(fd1 + logr.*j0m1d1 + j0m1.*rm1); -g2(isus) = const1*j0m1d2 - o2p*(fd2 + logr.*j0m1d2 + 2*j0m1d1.*rm1 - j0m1.*rm2); +r2fac = -(1-r2logrfac)*k*k*0.25; +g0(isus) = const1*(j0m1+1+r2fac*rsus.*rsus) - o2p*(f + logr.*j0m1); +g1(isus) = const1*(j0m1d1+2*r2fac*rsus) - o2p*(fd1 + logr.*j0m1d1 + j0m1.*rm1); +g21(isus) = const1*(j0m1d2-j0m1d1.*rm1) - o2p*(fd21 + logr.*(j0m1d2-j0m1d1.*rm1) + 2*j0m1d1.*rm1 - 2*j0m1.*rm2); g3(isus) = const1*j0m1d3 - o2p*(fd3 + logr.*j0m1d3 + 3*j0m1d2.*rm1 - ... 3*j0m1d1.*rm2 + 2*j0m1.*rm3); g4(isus) = const1*j0m1d4 - o2p*(fd4 + logr.*j0m1d4 + 4*j0m1d3.*rm1 - ... diff --git a/chunkie/+chnk/+flex2d/kern.m b/chunkie/+chnk/+flex2d/kern.m index 0893d6f..ae56994 100644 --- a/chunkie/+chnk/+flex2d/kern.m +++ b/chunkie/+chnk/+flex2d/kern.m @@ -1,7 +1,7 @@ function submat= kern(zk,srcinfo,targinfo,type,varargin) -%CHNK.HELM2D.KERN standard Helmholtz layer potential kernels in 2D +%CHNK.FLEX2D.KERN standard Modified biharmonic layer potential kernels in 2D % -% Syntax: submat = chnk.heml2d.kern(zk,srcinfo,targingo,type,varargin) +% Syntax: submat = chnk.flex2d.kern(zk,srcinfo,targingo,type,varargin) % % Let x be targets and y be sources for these formulas, with % n_x and n_y the corresponding unit normals at those points @@ -70,166 +70,15 @@ [~,ns] = size(src); [~,nt] = size(targ); -% double layer -if strcmpi(type,'d') - srcnorm = srcinfo.n; - [~,grad] = chnk.flex2d.green(zk,src,targ); - nx = repmat(srcnorm(1,:),nt,1); - ny = repmat(srcnorm(2,:),nt,1); - submat = -(grad(:,:,1).*nx + grad(:,:,2).*ny); -end - -% normal derivative of single layer -if strcmpi(type,'sprime') - targnorm = targinfo.n; - [~,grad] = chnk.flex2d.green(zk,src,targ); - nx = repmat((targnorm(1,:)).',1,ns); - ny = repmat((targnorm(2,:)).',1,ns); - submat = (grad(:,:,1).*nx + grad(:,:,2).*ny); -end - - -% Tangential derivative of single layer -if strcmpi(type,'stau') - targtan = targinfo.d; - [~,grad] = chnk.flex2d.green(zk,src,targ); - dx = repmat((targtan(1,:)).',1,ns); - dy = repmat((targtan(2,:)).',1,ns); - ds = sqrt(dx.*dx+dy.*dy); - submat = (grad(:,:,1).*dx + grad(:,:,2).*dy)./ds; -end - -% single layer -if strcmpi(type,'s') - submat = chnk.flex2d.green(zk,src,targ); -end - -% normal derivative of double layer -if strcmpi(type,'dprime') - targnorm = targinfo.n; - srcnorm = srcinfo.n; - [~,~,hess] = chnk.flex2d.green(zk,src,targ); - nxsrc = repmat(srcnorm(1,:),nt,1); - nysrc = repmat(srcnorm(2,:),nt,1); - nxtarg = repmat((targnorm(1,:)).',1,ns); - nytarg = repmat((targnorm(2,:)).',1,ns); - submat = -(hess(:,:,1).*nxsrc.*nxtarg + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)... - + hess(:,:,3).*nysrc.*nytarg); -end - -% Combined field -if strcmpi(type,'c') - srcnorm = srcinfo.n; - coef = ones(2,1); - if(nargin == 5); coef = varargin{1}; end - [submats,grad] = chnk.flex2d.green(zk,src,targ); - nx = repmat(srcnorm(1,:),nt,1); - ny = repmat(srcnorm(2,:),nt,1); - submatd = -(grad(:,:,1).*nx + grad(:,:,2).*ny); - submat = coef(1)*submatd + coef(2)*submats; -end - -% normal derivative of combined field -if strcmpi(type,'cprime') - coef = ones(2,1); - if(nargin == 5); coef = varargin{1}; end; - targnorm = targinfo.n; - srcnorm = srcinfo.n; - - submat = zeros(nt,ns); - - % Get gradient and hessian info - [submats,grad,hess] = chnk.flex2d.green(zk,src,targ); - - nxsrc = repmat(srcnorm(1,:),nt,1); - nysrc = repmat(srcnorm(2,:),nt,1); - nxtarg = repmat((targnorm(1,:)).',1,ns); - nytarg = repmat((targnorm(2,:)).',1,ns); - - % D' - submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ... - + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)... - + hess(:,:,3).*nysrc.*nytarg); - % S' - submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg); - - submat = coef(1)*submatdp + coef(2)*submatsp; -end - - -% Dirichlet and neumann data corresponding to combined field -if strcmpi(type,'c2trans') - coef = ones(2,1); - if(nargin == 5); coef = varargin{1}; end; - targnorm = targinfo.n; - srcnorm = srcinfo.n; - - % Get gradient and hessian info - [submats,grad,hess] = chnk.flex2d.green(zk,src,targ); - - nxsrc = repmat(srcnorm(1,:),nt,1); - nysrc = repmat(srcnorm(2,:),nt,1); - nxtarg = repmat((targnorm(1,:)).',1,ns); - nytarg = repmat((targnorm(2,:)).',1,ns); - % D' - submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ... - + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)... - + hess(:,:,3).*nysrc.*nytarg); - % S' - submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg); - % D - submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc); - submat = zeros(2*nt,ns); - - submat(1:2:2*nt,:) = coef(1)*submatd + coef(2)*submats; - submat(2:2:2*nt,:) = coef(1)*submatdp + coef(2)*submatsp; - -end - - -% all kernels, [c11 D, c12 S; c21 D', c22 S'] -if strcmpi(type,'all') - - targnorm = targinfo.n; - srcnorm = srcinfo.n; - cc = varargin{1}; - - submat = zeros(2*nt,2*ns); - - % Get gradient and hessian info - [submats,grad,hess] = chnk.flex2d.green(zk,src,targ); - - nxsrc = repmat(srcnorm(1,:),nt,1); - nysrc = repmat(srcnorm(2,:),nt,1); - nxtarg = repmat((targnorm(1,:)).',1,ns); - nytarg = repmat((targnorm(2,:)).',1,ns); - - submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ... - + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)... - + hess(:,:,3).*nysrc.*nytarg); - % S' - submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg); - % D - submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc); - - - submat(1:2:2*nt,1:2:2*ns) = submatd*cc(1,1); - submat(1:2:2*nt,2:2:2*ns) = submats*cc(1,2); - submat(2:2:2*nt,1:2:2*ns) = submatdp*cc(2,1); - submat(2:2:2*nt,2:2:2*ns) = submatsp*cc(2,2); -end - - -% clamped plate kernel for modified biharmonic problem, K21 Kernel is -% handled separately. +% clamped plate kernel for modified biharmonic problem if strcmpi(type, 'clamped-plate') srcnorm = srcinfo.n; srctang = srcinfo.d; targnorm = targinfo.n; - [hess, third, forth] = chnk.flex2d.thirdforth_derivatives(zk, src, targ); % Hankel part + [~, ~, hess, third, forth] = chnk.flex2d.helmdiffgreen(zk, src, targ); % Hankel part nx = repmat(srcnorm(1,:),nt,1); ny = repmat(srcnorm(2,:),nt,1); @@ -237,36 +86,51 @@ nytarg = repmat((targnorm(2,:)).',1,ns); zkimag = (1i)*zk; - [hessK, thirdK, forthK] = chnk.flex2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + [~, ~, hessK, thirdK, forthK] = chnk.flex2d.helmdiffgreen(zkimag, src, targ); % modified bessel K part dx = repmat(srctang(1,:),nt,1); dy = repmat(srctang(2,:),nt,1); - ds = sqrt(dx.*dx+dy.*dy); + ds = sqrt(dx.*dx+dy.*dy); - taux = dx./ds; - tauy = dy./ds; + taux = dx./ds; + tauy = dy./ds; - Kxx = 1/(2*zk^2).*(third(:, :, 1).*(nx.*nx.*nx) + third(:, :, 2).*(3*nx.*nx.*ny) +... - third(:, :, 4).*(3*nx.*ny.*ny) + third(:, :, 6).*(ny.*ny.*ny)) - ... + Kxx = -(1/(2*zk^2).*(third(:, :, 1).*(nx.*nx.*nx) + third(:, :, 2).*(3*nx.*nx.*ny) +... + third(:, :, 3).*(3*nx.*ny.*ny) + third(:, :, 4).*(ny.*ny.*ny)) - ... 1/(2*zk^2).*(thirdK(:, :, 1).*(nx.*nx.*nx) + thirdK(:, :, 2).*(3*nx.*nx.*ny) +... - thirdK(:, :, 4).*(3*nx.*ny.*ny) + thirdK(:, :, 6).*(ny.*ny.*ny)) + ... + thirdK(:, :, 3).*(3*nx.*ny.*ny) + thirdK(:, :, 4).*(ny.*ny.*ny))) - ... (3/(2*zk^2).*(third(:, :, 1).*(nx.*taux.*taux) + third(:, :, 2).*(2*nx.*taux.*tauy + ny.*taux.*taux) +... - third(:, :, 4).*(nx.*tauy.*tauy + 2*ny.*taux.*tauy) + third(:, :, 6).*(ny.*tauy.*tauy))-... + third(:, :, 3).*(nx.*tauy.*tauy + 2*ny.*taux.*tauy) + third(:, :, 4).*(ny.*tauy.*tauy))-... 3/(2*zk^2).*(thirdK(:, :, 1).*(nx.*taux.*taux) + thirdK(:, :, 2).*(2*nx.*taux.*tauy + ny.*taux.*taux) +... - thirdK(:, :, 4).*(nx.*tauy.*tauy + 2*ny.*taux.*tauy) + thirdK(:, :, 6).*(ny.*tauy.*tauy))); % G_{ny ny ny} + 3G_{ny tauy tauy} + thirdK(:, :, 3).*(nx.*tauy.*tauy + 2*ny.*taux.*tauy) + thirdK(:, :, 4).*(ny.*tauy.*tauy))); % G_{ny ny ny} + 3G_{ny tauy tauy} Kxy = -(1/(2*zk^2).*(hess(:, :, 1).*(nx.*nx) + hess(:, :, 2).*(2*nx.*ny) + hess(:, :, 3).*(ny.*ny))-... 1/(2*zk^2).*(hessK(:, :, 1).*(nx.*nx) + hessK(:, :, 2).*(2*nx.*ny) + hessK(:, :, 3).*(ny.*ny)))+... (1/(2*zk^2).*(hess(:, :, 1).*(taux.*taux) + hess(:, :, 2).*(2*taux.*tauy) + hess(:, :, 3).*(tauy.*tauy))-... 1/(2*zk^2).*(hessK(:, :, 1).*(taux.*taux) + hessK(:, :, 2).*(2*taux.*tauy) + hessK(:, :, 3).*(tauy.*tauy))); % -G_{ny ny} + G_{tauy tauy} - Kyx = 0; - Kyy = 1/(2*zk^2).*(third(:,:, 1).*(nx.*nx.*nxtarg) +third(:, :, 2).*(nx.*nx.*nytarg + 2*nx.*ny.*nxtarg) + third(:, :, 4).*(2*nx.*ny.*nytarg + ny.*ny.*nxtarg)+... - third(:, :,6).*(ny.*ny.*nytarg)) - 1/(2*zk^2).*(thirdK(:,:, 1).*(nx.*nx.*nxtarg) +thirdK(:, :, 2).*(nx.*nx.*nytarg + 2*nx.*ny.*nxtarg) + thirdK(:, :, 4).*(2*nx.*ny.*nytarg + ny.*ny.*nxtarg)+... - thirdK(:, :,6).*(ny.*ny.*nytarg)) - 1/(2*zk^2).*(third(:,:, 1).*(taux.*taux.*nxtarg) +third(:, :, 2).*(taux.*taux.*nytarg + 2*taux.*tauy.*nxtarg) + third(:, :, 4).*(2*taux.*tauy.*nytarg + tauy.*tauy.*nxtarg)+... - third(:, :,6).*(tauy.*tauy.*nytarg)) + 1/(2*zk^2).*(thirdK(:,:, 1).*(taux.*taux.*nxtarg) +thirdK(:, :, 2).*(taux.*taux.*nytarg + 2*taux.*tauy.*nxtarg) + thirdK(:, :, 4).*(2*taux.*tauy.*nytarg + tauy.*tauy.*nxtarg)+... - thirdK(:, :,6).*(tauy.*tauy.*nytarg)); + Kyx = -(1/(2*zk^2).*(forth(:, :, 1).*(nx.*nx.*nx.*nxtarg) + forth(:, :, 2).*(nx.*nx.*nx.*nytarg + 3*nx.*nx.*ny.*nxtarg) + ... + forth(:, :, 3).*(3*nx.*nx.*ny.*nytarg + 3*nx.*ny.*ny.*nxtarg) + forth(:, :, 4).*(3*nx.*ny.*ny.*nytarg +ny.*ny.*ny.*nxtarg)+... + forth(:, :, 5).*(ny.*ny.*ny.*nytarg)) - ... + 1/(2*zk^2).*(forthK(:, :, 1).*(nx.*nx.*nx.*nxtarg) + forthK(:, :, 2).*(nx.*nx.*nx.*nytarg + 3*nx.*nx.*ny.*nxtarg) + ... + forthK(:, :, 3).*(3*nx.*nx.*ny.*nytarg + 3*nx.*ny.*ny.*nxtarg) + forthK(:, :, 4).*(3*nx.*ny.*ny.*nytarg +ny.*ny.*ny.*nxtarg)+... + forthK(:, :, 5).*(ny.*ny.*ny.*nytarg))) - ... + (3/(2*zk^2).*(forth(:, :, 1).*(nx.*taux.*taux.*nxtarg)+ forth(:, :, 2).*(nx.*taux.*taux.*nytarg + 2*nx.*taux.*tauy.*nxtarg + ny.*taux.*taux.*nxtarg) +... + forth(:, :, 3).*(2*nx.*taux.*tauy.*nytarg + ny.*taux.*taux.*nytarg + nx.*tauy.*tauy.*nxtarg + 2*ny.*taux.*tauy.*nxtarg) + ... + forth(:, :, 4).*(nx.*tauy.*tauy.*nytarg +2*ny.*taux.*tauy.*nytarg + ny.*tauy.*tauy.*nxtarg) +... + forth(:, :, 5).*(ny.*tauy.*tauy.*nytarg))-... + 3/(2*zk^2).*(forthK(:, :, 1).*(nx.*taux.*taux.*nxtarg)+ forthK(:, :, 2).*(nx.*taux.*taux.*nytarg + 2*nx.*taux.*tauy.*nxtarg + ny.*taux.*taux.*nxtarg) +... + forthK(:, :, 3).*(2*nx.*taux.*tauy.*nytarg + ny.*taux.*taux.*nytarg + nx.*tauy.*tauy.*nxtarg + 2*ny.*taux.*tauy.*nxtarg) + ... + forthK(:, :, 4).*(nx.*tauy.*tauy.*nytarg +2*ny.*taux.*tauy.*nytarg + ny.*tauy.*tauy.*nxtarg) + ... + forthK(:, :, 5).*(ny.*tauy.*tauy.*nytarg))); + + Kyy = -(1/(2*zk^2).*(third(:,:, 1).*(nx.*nx.*nxtarg) +third(:, :, 2).*(nx.*nx.*nytarg + 2*nx.*ny.*nxtarg) + third(:, :, 3).*(2*nx.*ny.*nytarg + ny.*ny.*nxtarg)+... + third(:, :,4).*(ny.*ny.*nytarg)) -... + 1/(2*zk^2).*(thirdK(:,:, 1).*(nx.*nx.*nxtarg) +thirdK(:, :, 2).*(nx.*nx.*nytarg + 2*nx.*ny.*nxtarg) + thirdK(:, :, 3).*(2*nx.*ny.*nytarg + ny.*ny.*nxtarg)+... + thirdK(:, :,4).*(ny.*ny.*nytarg))) + (1/(2*zk^2).*(third(:,:, 1).*(taux.*taux.*nxtarg) +third(:, :, 2).*(taux.*taux.*nytarg + 2*taux.*tauy.*nxtarg) + third(:, :, 3).*(2*taux.*tauy.*nytarg + tauy.*tauy.*nxtarg)+... + third(:, :,4).*(tauy.*tauy.*nytarg)) - 1/(2*zk^2).*(thirdK(:,:, 1).*(taux.*taux.*nxtarg) +thirdK(:, :, 2).*(taux.*taux.*nytarg + 2*taux.*tauy.*nxtarg) + thirdK(:, :, 3).*(2*taux.*tauy.*nytarg + tauy.*tauy.*nxtarg)+... + thirdK(:, :,4).*(tauy.*tauy.*nytarg))); submat = zeros(2*nt,2*ns); @@ -277,122 +141,6 @@ submat(2:2:end,2:2:end) = Kyy; end -if strcmpi(type, 'clamped-plate K21') - srcnorm = srcinfo.n; - srctang = srcinfo.d; - targnorm = targinfo.n; - R = 0.001; - - [~, ~, forth] = chnk.flex2d.thirdforth_derivatives(zk, src, targ); % Hankel part - nx = repmat(srcnorm(1,:),nt,1); - ny = repmat(srcnorm(2,:),nt,1); - - nxtarg = repmat((targnorm(1,:)).',1,ns); - nytarg = repmat((targnorm(2,:)).',1,ns); - - zkimag = (1i)*zk; - [~, ~, forthK] = chnk.flex2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part - - dx = repmat(srctang(1,:),nt,1); - dy = repmat(srctang(2,:),nt,1); - - ds = sqrt(dx.*dx+dy.*dy); - - taux = dx./ds; - tauy = dy./ds; - - rx = targ(1,:).' - src(1,:); - ry = targ(2,:).' - src(2,:); - r2 = rx.^2 + ry.^2; - dists = sqrt(r2); - inds = and((dists < R),(dists > 0)); - - - submat = 1/(2*zk^2).*(forth(:, :, 1).*(nx.*nx.*nx.*nxtarg) + forth(:, :, 2).*(nx.*nx.*nx.*nytarg + 3*nx.*nx.*ny.*nxtarg) + ... - forth(:, :, 4).*(3*nx.*nx.*ny.*nytarg + 3*nx.*ny.*ny.*nxtarg) + forth(:, :, 6).*(3*nx.*ny.*ny.*nytarg +ny.*ny.*ny.*nxtarg)+... - forth(:, :, 8).*(ny.*ny.*ny.*nytarg)) - 1/(2*zk^2).*(forthK(:, :, 1).*(nx.*nx.*nx.*nxtarg) + forthK(:, :, 2).*(nx.*nx.*nx.*nytarg + 3*nx.*nx.*ny.*nxtarg) + ... - forthK(:, :, 4).*(3*nx.*nx.*ny.*nytarg + 3*nx.*ny.*ny.*nxtarg) + forthK(:, :, 6).*(3*nx.*ny.*ny.*nytarg +ny.*ny.*ny.*nxtarg)+... - forthK(:, :, 8).*(ny.*ny.*ny.*nytarg)) + ... - (3/(2*zk^2).*(forth(:, :, 1).*(nx.*taux.*taux.*nxtarg)+ forth(:, :, 2).*(nx.*taux.*taux.*nytarg + 2*nx.*taux.*tauy.*nxtarg + ny.*taux.*taux.*nxtarg) +... - forth(:, :, 4).*(2*nx.*taux.*tauy.*nytarg + ny.*taux.*taux.*nytarg + nx.*tauy.*tauy.*nxtarg + 2*ny.*taux.*tauy.*nxtarg) + ... - forth(:, :, 6).*(nx.*tauy.*tauy.*nytarg +2*ny.*taux.*tauy.*nytarg + ny.*tauy.*tauy.*nxtarg) + forth(:, :, 8).*(ny.*tauy.*tauy.*nytarg))-... - 3/(2*zk^2).*(forthK(:, :, 1).*(nx.*taux.*taux.*nxtarg)+ forthK(:, :, 2).*(nx.*taux.*taux.*nytarg + 2*nx.*taux.*tauy.*nxtarg + ny.*taux.*taux.*nxtarg) +... - forthK(:, :, 4).*(2*nx.*taux.*tauy.*nytarg + ny.*taux.*taux.*nytarg + nx.*tauy.*tauy.*nxtarg + 2*ny.*taux.*tauy.*nxtarg) + ... - forthK(:, :, 6).*(nx.*tauy.*tauy.*nytarg +2*ny.*taux.*tauy.*nytarg + ny.*tauy.*tauy.*nxtarg) + forthK(:, :, 8).*(ny.*tauy.*tauy.*nytarg))); - % - if any(inds(:)) % if statement for evalauation points that are closed to the boundary - temp = 1i*imag(submat(inds)); - - nx = nx(inds); ny = ny(inds); - - nxtarg = nxtarg(inds); - nytarg = nytarg(inds); - - dx = dx(inds); - dy = dy(inds); - - - ds = sqrt(dx.*dx+dy.*dy); - - - taux = dx./ds; % normalization - tauy = dy./ds; - - - rx = rx(inds); - ry = ry(inds); - r2 = r2(inds); - - - normaldot = nx.*nxtarg + ny.*nytarg; - - - rn = rx.*nx + ry.*ny; - - rntarg = rx.*nxtarg + ry.*nytarg; - - rtau = rx.*taux + ry.*tauy; - - - taudotnsrc = taux.*nx + tauy.*ny; - taudotntarg = taux.*nxtarg + tauy.*nytarg; - - [smoothfirst, smoothsecond] = chnk.flex2d.smooth_part_derivatives(zk, rx, ry); - - eulerconstantpart = (smoothfirst(:, :, 1)).*(nx.*nx.*nx.*nxtarg) +... - (smoothfirst(:, :, 2)).*(nx.*nx.*nx.*nytarg + 3*nx.*nx.*ny.*nxtarg) + ... - (smoothfirst(:, :, 3)).*(3*nx.*nx.*ny.*nytarg + 3*nx.*ny.*ny.*nxtarg) +... - (smoothfirst(:, :, 4)).*(3*nx.*ny.*ny.*nytarg +ny.*ny.*ny.*nxtarg)+... - (smoothfirst(:, :, 5)).*(ny.*ny.*ny.*nytarg) + ... - 3.*((smoothfirst(:, :, 1)).*(nx.*taux.*taux.*nxtarg) +... - (smoothfirst(:, :, 2)).*(nx.*taux.*taux.*nytarg + 2*nx.*taux.*tauy.*nxtarg + ny.*taux.*taux.*nxtarg) +... - (smoothfirst(:, :, 3)).*(2*nx.*taux.*tauy.*nytarg + ny.*taux.*taux.*nytarg + nx.*tauy.*tauy.*nxtarg + 2*ny.*taux.*tauy.*nxtarg) + ... - (smoothfirst(:, :, 4)).*(nx.*tauy.*tauy.*nytarg +2*ny.*taux.*tauy.*nytarg + ny.*tauy.*tauy.*nxtarg) +... - (smoothfirst(:, :, 5)).*(ny.*tauy.*tauy.*nytarg)); - - puresmoothpart = (smoothsecond(:, :, 1)).*(nx.*nx.*nx.*nxtarg) +... - (smoothsecond(:, :, 2)).*(nx.*nx.*nx.*nytarg + 3*nx.*nx.*ny.*nxtarg) + ... - (smoothsecond(:, :, 3)).*(3*nx.*nx.*ny.*nytarg + 3*nx.*ny.*ny.*nxtarg) +... - (smoothsecond(:, :, 4)).*(3*nx.*ny.*ny.*nytarg +ny.*ny.*ny.*nxtarg)+... - (smoothsecond(:, :, 5)).*(ny.*ny.*ny.*nytarg) + ... - 3.*((smoothsecond(:, :, 1)).*(nx.*taux.*taux.*nxtarg) +... - (smoothsecond(:, :, 2)).*(nx.*taux.*taux.*nytarg + 2*nx.*taux.*tauy.*nxtarg + ny.*taux.*taux.*nxtarg) +... - (smoothsecond(:, :, 3)).*(2*nx.*taux.*tauy.*nytarg + ny.*taux.*taux.*nytarg + nx.*tauy.*tauy.*nxtarg + 2*ny.*taux.*tauy.*nxtarg) + ... - (smoothsecond(:, :, 4)).*(nx.*tauy.*tauy.*nytarg +2*ny.*taux.*tauy.*nytarg + ny.*tauy.*tauy.*nxtarg) +... - (smoothsecond(:, :, 5)).*(ny.*tauy.*tauy.*nytarg)); - - submat(inds) = -(3/(4*pi)).*(normaldot)./r2 + (3/(2*pi)).*(rntarg.*rn)./(r2.^2) +... - (3/(2*pi)).*((rn.^2).*normaldot)./(r2.^2) - (2/pi).*((rn.^3).*rntarg)./(r2.^3) + ... - 3.*(-1/(2*pi).*(taudotntarg.*taudotnsrc)./r2 + (1/pi).*(rntarg.*rtau.*taudotnsrc)./(r2.^2) + ... - (1/(2*pi)).*(normaldot.*(rtau.^2))./(r2.^2) + (1/pi).*rtau.*rn.*taudotntarg./(r2.^2) -... - (2/pi).*(rntarg.*rn.*(rtau.^2))./(r2.^3) - (1/(4*pi)).*normaldot./(r2) + (1/(2*pi)).*(rntarg.*rn)./(r2.^2)) + ... - temp + eulerconstantpart + puresmoothpart; - - - end - - -end % free plate kernel for the modified biharmonic problem (K11 with no % hilbert transform subtraction, K12 kernel, K22 with no curvature part) K21 is @@ -494,9 +242,7 @@ coefs = varargin{1}; %R = 0.001; % radius of kernel replacement - - - [~, ~,~, ~, forth] = chnk.flex2d.helmdiffgreen(zk, src, targ); % Hankel part + [~, ~,~, ~, forth] = chnk.flex2d.helmdiffgreen(zk, src, targ, true); % Hankel part nx = repmat(srcnorm(1,:),nt,1); ny = repmat(srcnorm(2,:),nt,1); @@ -504,11 +250,7 @@ nytarg = repmat((targnorm(2,:)).',1,ns); zkimag = (1i)*zk; - [~, ~, ~,~, forthK] = chnk.flex2d.helmdiffgreen(zkimag, src, targ); % modified bessel K part - - - - + [~, ~, ~,~, forthK] = chnk.flex2d.helmdiffgreen(zkimag, src, targ, true); % modified bessel K part dx = repmat(srctang(1,:),nt,1); dy = repmat(srctang(2,:),nt,1); @@ -526,11 +268,22 @@ tauxtarg = dx1./ds1; tauytarg = dy1./ds1; + tauxtargnsrc = tauxtarg.*nx + tauytarg.*ny; + tauxtargntarg = tauxtarg.*nxtarg + tauytarg.*nytarg; + + rx = targ(1,:).' - src(1,:); ry = targ(2,:).' - src(2,:); r2 = rx.^2 + ry.^2; - dists = sqrt(r2); - %inds = and((dists <= R),(dists > 0)); + normaldot = nx.*nxtarg + ny.*nytarg; + + rn = rx.*nx + ry.*ny; + + rntarg = rx.*nxtarg + ry.*nytarg; + + + rtautarg = rx.*tauxtarg + ry.*tauytarg; + submat = -(1/(2*zk^2).*(forth(:, :, 1).*(nxtarg.*nxtarg.*nxtarg.*nx) + forth(:, :, 2).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... forth(:, :, 3).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) +... @@ -539,7 +292,7 @@ 1/(2*zk^2).*(forthK(:, :, 1).*(nxtarg.*nxtarg.*nxtarg.*nx) + forthK(:, :, 2).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... forthK(:, :, 3).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) +... forthK(:, :, 4).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... - forthK(:, :, 5).*(nytarg.*nytarg.*nytarg.*ny))) + (3/(4*pi)).*(nxtarg.*nx + nytarg.*ny)./r2 - ... + forthK(:, :, 5).*(nytarg.*nytarg.*nytarg.*ny))) - ... ((2-coefs(1))/(2*zk^2).*(forth(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg.*nx) + forth(:, :, 2).*(tauxtarg.*tauxtarg.*nxtarg.*ny + tauxtarg.*tauxtarg.*nytarg.*nx + 2*tauxtarg.*tauytarg.*nxtarg.*nx)+... forth(:, :, 3).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... forth(:, :, 4).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... @@ -547,92 +300,97 @@ (2-coefs(1))/(2*zk^2).*(forthK(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg.*nx) + forthK(:, :, 2).*(tauxtarg.*tauxtarg.*nxtarg.*ny + tauxtarg.*tauxtarg.*nytarg.*nx + 2*tauxtarg.*tauytarg.*nxtarg.*nx)+... forthK(:, :, 3).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... forthK(:, :, 4).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... - forthK(:, :, 5).*(tauytarg.*tauytarg.*nytarg.*ny))) +... - (2-coefs(1))/(4*pi).*(nxtarg.*nx + nytarg.*ny)./r2 - ... - ((2-coefs(1))/(2*pi)).*(nxtarg.*nx + nytarg.*ny).*(rx.*tauxtarg + ry.*tauytarg).*(rx.*tauxtarg + ry.*tauytarg)./(r2.^2) - ... + forthK(:, :, 5).*(tauytarg.*tauytarg.*nytarg.*ny))) + ... + 3/(2*pi).*(rn.*rntarg./(r2.^2)) + 3/(2*pi).*(rntarg.^2).*(normaldot)./(r2.^2) -... + (2/pi).*(rn.*(rntarg.^3))./(r2.^3) + ... + (2-coefs(1)).*(-1/(2*pi).*(tauxtargnsrc).*(tauxtargntarg)./(r2) +... + 1/(pi).*(rn.*rtautarg.*tauxtargntarg)./(r2.^2) + ... + 1/(pi).*(rtautarg.*tauxtargnsrc.*rntarg)./(r2.^2) -.... + 2/(pi).*(rn.*rntarg.*(rtautarg.*rtautarg))./(r2.^3) + ... + 1/(2*pi).*(rn.*rntarg)./(r2.^2)) -... ((1+coefs(1))/(4*pi)).*((taux.*tauxtarg + tauy.*tauytarg)./(r2) - 2*(rx.*tauxtarg + ry.*tauytarg).*(rx.*taux + ry.*tauy)./(r2.^2))-... (3/(4*pi)).*(nxtarg.*nx + nytarg.*ny)./r2 - (2-coefs(1))/(4*pi).*(nxtarg.*nx + nytarg.*ny)./r2 + ... (2-coefs(1))/(2*pi).*(nxtarg.*nx + nytarg.*ny).*(rx.*tauxtarg + ry.*tauytarg).*(rx.*tauxtarg + ry.*tauytarg)./(r2.^2); % third kernel with singularity subtraction and H[delta'] % - % if any(inds(:)) % if statement for evalauation points that are closed to the boundary - % temp = 1i*imag(submat(inds)); - % - % nx = nx(inds); ny = ny(inds); - % - % nxtarg = nxtarg(inds); - % nytarg = nytarg(inds); - % - % dx = dx(inds); - % dy = dy(inds); - % - % dx1 = dx1(inds); - % dy1 = dy1(inds); - % - % ds = sqrt(dx.*dx+dy.*dy); - % ds1 = sqrt(dx1.*dx1+dy1.*dy1); - % - % taux = dx./ds; % normalization - % tauy = dy./ds; - % - % tauxtarg = dx1./ds1; - % tauytarg = dy1./ds1; - % - % rx = rx(inds); - % ry = ry(inds); - % r2 = r2(inds); - % - % tauxtargnsrc = tauxtarg.*nx + tauytarg.*ny; - % tauxtargntarg = tauxtarg.*nxtarg + tauytarg.*nytarg; - % - % normaldot = nx.*nxtarg + ny.*nytarg; - % tangentdot = taux.*tauxtarg + tauy.*tauytarg; - % - % rn = rx.*nx + ry.*ny; - % - % rntarg = rx.*nxtarg + ry.*nytarg; - % - % rtau = rx.*taux + ry.*tauy; - % - % rtautarg = rx.*tauxtarg + ry.*tauytarg; - % - % - % - % [smoothfirst, smoothsecond] = chnk.flex2d.smooth_part_derivatives(zk, rx, ry); - % - % eulerconstantpart = (smoothfirst(:, :, 1)).*(nxtarg.*nxtarg.*nxtarg.*nx) +... - % (smoothfirst(:, :, 2)).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... - % (smoothfirst(:, :, 3)).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) +... - % (smoothfirst(:, :, 4)).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... - % (smoothfirst(:, :, 5)).*(nytarg.*nytarg.*nytarg.*ny)+ ... - % (2-coefs(1)).*((smoothfirst(:, :, 1)).*(tauxtarg.*tauxtarg.*nxtarg.*nx) +... - % (smoothfirst(:, :, 2)).*(tauxtarg.*tauxtarg.*nxtarg.*ny + tauxtarg.*tauxtarg.*nytarg.*nx + 2*tauxtarg.*tauytarg.*nxtarg.*nx)+... - % (smoothfirst(:, :, 3)).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... - % (smoothfirst(:, :, 4)).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... - % (smoothfirst(:, :, 5)).*(tauytarg.*tauytarg.*nytarg.*ny)); - % - % puresmoothpart = (smoothsecond(:, :, 1)).*(nxtarg.*nxtarg.*nxtarg.*nx) +... - % (smoothsecond(:, :, 2)).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... - % (smoothsecond(:, :, 3)).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) +... - % (smoothsecond(:, :, 4)).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... - % (smoothsecond(:, :, 5)).*(nytarg.*nytarg.*nytarg.*ny)+ ... - % (2-coefs(1)).*((smoothsecond(:, :, 1)).*(tauxtarg.*tauxtarg.*nxtarg.*nx) +... - % (smoothsecond(:, :, 2)).*(tauxtarg.*tauxtarg.*nxtarg.*ny + tauxtarg.*tauxtarg.*nytarg.*nx + 2*tauxtarg.*tauytarg.*nxtarg.*nx)+... - % (smoothsecond(:, :, 3)).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... - % (smoothsecond(:, :, 4)).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... - % (smoothsecond(:, :, 5)).*(tauytarg.*tauytarg.*nytarg.*ny)); + +end + + +% free plate kernel K21 for the modified biharmonic problem. This part +% handles the singularity subtraction and swap the evaluation to its +% asymptotic expansions if the targets and sources are close. + +if strcmpi(type, 'free plate K21 first part unsubtract') + srcnorm = srcinfo.n; + srctang = srcinfo.d; + targnorm = targinfo.n; + targtang = targinfo.d; + coefs = varargin{1}; + %R = 0.001; % radius of kernel replacement + + [~, ~,~, ~, forth] = chnk.flex2d.helmdiffgreen(zk, src, targ); % Hankel part + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + zkimag = (1i)*zk; + [~, ~, ~,~, forthK] = chnk.flex2d.helmdiffgreen(zkimag, src, targ); % modified bessel K part + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + dx1 = repmat((targtang(1,:)).',1,ns); + dy1 = repmat((targtang(2,:)).',1,ns); + + + ds = sqrt(dx.*dx+dy.*dy); + ds1 = sqrt(dx1.*dx1+dy1.*dy1); + + taux = dx./ds; % normalization + tauy = dy./ds; + + tauxtarg = dx1./ds1; + tauytarg = dy1./ds1; + + tauxtargnsrc = tauxtarg.*nx + tauytarg.*ny; + tauxtargntarg = tauxtarg.*nxtarg + tauytarg.*nytarg; + + + rx = targ(1,:).' - src(1,:); + ry = targ(2,:).' - src(2,:); + r2 = rx.^2 + ry.^2; + normaldot = nx.*nxtarg + ny.*nytarg; + + rn = rx.*nx + ry.*ny; + + rntarg = rx.*nxtarg + ry.*nytarg; + + + rtautarg = rx.*tauxtarg + ry.*tauytarg; + + + submat = -(1/(2*zk^2).*(forth(:, :, 1).*(nxtarg.*nxtarg.*nxtarg.*nx) + forth(:, :, 2).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... + forth(:, :, 3).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) +... + forth(:, :, 4).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... + forth(:, :, 5).*(nytarg.*nytarg.*nytarg.*ny)) -... + 1/(2*zk^2).*(forthK(:, :, 1).*(nxtarg.*nxtarg.*nxtarg.*nx) + forthK(:, :, 2).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... + forthK(:, :, 3).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) +... + forthK(:, :, 4).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... + forthK(:, :, 5).*(nytarg.*nytarg.*nytarg.*ny))) - ... + ((2-coefs(1))/(2*zk^2).*(forth(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg.*nx) + forth(:, :, 2).*(tauxtarg.*tauxtarg.*nxtarg.*ny + tauxtarg.*tauxtarg.*nytarg.*nx + 2*tauxtarg.*tauytarg.*nxtarg.*nx)+... + forth(:, :, 3).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... + forth(:, :, 4).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... + forth(:, :, 5).*(tauytarg.*tauytarg.*nytarg.*ny)) -... + (2-coefs(1))/(2*zk^2).*(forthK(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg.*nx) + forthK(:, :, 2).*(tauxtarg.*tauxtarg.*nxtarg.*ny + tauxtarg.*tauxtarg.*nytarg.*nx + 2*tauxtarg.*tauytarg.*nxtarg.*nx)+... + forthK(:, :, 3).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... + forthK(:, :, 4).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... + forthK(:, :, 5).*(tauytarg.*tauytarg.*nytarg.*ny))); + % third kernel with singularity subtraction and H[delta'] % - % submat(inds) = 3/(2*pi).*(rn.*rntarg./(r2.^2)) + 3/(2*pi).*(rntarg.^2).*(normaldot)./(r2.^2) -... - % (2/pi).*(rn.*(rntarg.^3))./(r2.^3) + ... - % (2-coefs(1)).*(-1/(2*pi).*(tauxtargnsrc).*(tauxtargntarg)./(r2) +... - % 1/(pi).*(rn.*rtautarg.*tauxtargntarg)./(r2.^2) + ... - % 1/(pi).*(rtautarg.*tauxtargnsrc.*rntarg)./(r2.^2) -.... - % 2/(pi).*(rn.*rntarg.*(rtautarg.*rtautarg))./(r2.^3) + ... - % 1/(2*pi).*(rn.*rntarg)./(r2.^2)) - 3/(4*pi).*normaldot./r2 +... - % ((2-coefs(1))/(2*pi)).*((rtautarg.*rtautarg).*(normaldot))./(r2.^2) -... - % ((2-coefs(1))/(4*pi)).*(normaldot)./r2 - ... - % ((1+coefs(1))/(4*pi)).*((tangentdot)./(r2) - 2*(rtautarg).*(rtau)./(r2.^2)) + temp + eulerconstantpart + puresmoothpart; - % end + end % kernels in K11 with hilbert transform subtractions. @@ -693,6 +451,61 @@ end +if strcmpi(type, 'free plate hilbert unsubtract') + srcnorm = srcinfo.n; + srctang = srcinfo.d; + targnorm = targinfo.n; + targtang = targinfo.d; + coefs = varargin{1}; + + [~,~,~, third, ~] = chnk.flex2d.helmdiffgreen(zk, src, targ); % Hankel part + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + zkimag = (1i)*zk; + [~,~, ~, thirdK, ~] = chnk.flex2d.helmdiffgreen(zkimag, src, targ); % modified bessel K part + [~,grad] = chnk.lap2d.green(src,targ,true); + + + + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + dx1 = repmat((targtang(1,:)).',1,ns); + dy1 = repmat((targtang(2,:)).',1,ns); + + + ds = sqrt(dx.*dx+dy.*dy); + ds1 = sqrt(dx1.*dx1+dy1.*dy1); + + taux = dx./ds; % normalization + tauy = dy./ds; + + tauxtarg = dx1./ds1; + tauytarg = dy1./ds1; + + + hilb = 2*(grad(:,:,1).*ny - grad(:,:,2).*nx); + + submat = -((1+ coefs(1))/2)*(1./(2*zk^2).*(third(:, :, 1).*(nxtarg.*nxtarg.*taux) + third(:, :, 2).*(nxtarg.*nxtarg.*tauy+ 2*nxtarg.*nytarg.*taux) +... + third(:, :, 3).*(2*nxtarg.*nytarg.*tauy +nytarg.*nytarg.*taux) +... + third(:, :, 4).*(nytarg.*nytarg.*tauy)) - ... + 1./(2*zk^2).*(thirdK(:, :, 1).*(nxtarg.*nxtarg.*taux) + thirdK(:, :, 2).*(nxtarg.*nxtarg.*tauy+ 2*nxtarg.*nytarg.*taux) +... + thirdK(:, :, 3).*(2*nxtarg.*nytarg.*tauy +nytarg.*nytarg.*taux) +... + thirdK(:, :, 4).*(nytarg.*nytarg.*tauy))) ... % + ((1+ coefs(1))/2).*0.25*hilb -... + -((1+ coefs(1))/2)*coefs(1).*(1./(2*zk^2).*(third(:, :, 1).*(tauxtarg.*tauxtarg.*taux) + third(:, :, 2).*(tauxtarg.*tauxtarg.*tauy + 2*tauxtarg.*tauytarg.*taux) +... + third(:, :, 3).*(2*tauxtarg.*tauytarg.*tauy + tauytarg.*tauytarg.*taux) +... + third(:, :, 4).*(tauytarg.*tauytarg.*tauy)) - ... + 1./(2*zk^2).*(thirdK(:, :, 1).*(tauxtarg.*tauxtarg.*taux) + thirdK(:, :, 2).*(tauxtarg.*tauxtarg.*tauy + 2*tauxtarg.*tauytarg.*taux) +... + thirdK(:, :, 3).*(2*tauxtarg.*tauytarg.*tauy + tauytarg.*tauytarg.*taux) +... + thirdK(:, :, 4).*(tauytarg.*tauytarg.*tauy))); %+ ((1+ coefs(1))/2)*coefs(1).*0.25*hilb; % hilbert subtraction + +end + % kernels in K21 that are coupled with Hilbert transforms. % (i.e beta* G_{nx nx nx tauy} + (2-nu)*beta*G_{taux taux tauy nx}) if strcmpi(type, 'free plate coupled hilbert') @@ -915,7 +728,7 @@ - [~, third, ~] = chnk.flex2d.thirdforth_derivatives(zk, src, targ); % Hankel part + [~, third, ~] = chnk.flex2d.helmdiffgreen(zk, src, targ, true); % Hankel part nx = repmat(srcnorm(1,:),nt,1); ny = repmat(srcnorm(2,:),nt,1); @@ -923,7 +736,7 @@ nytarg = repmat((targnorm(2,:)).',1,ns); zkimag = (1i)*zk; - [~, thirdK, ~] = chnk.flex2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + [~, thirdK, ~] = chnk.flex2d.helmdiffgreen(zkimag, src, targ, true); % modified bessel K part [~,grad] = chnk.lap2d.green(src,targ,true); % lap kernels for the adjoint of Hilbert transform dx = repmat(srctang(1,:),nt,1); @@ -969,17 +782,16 @@ end -if strcmpi(type, 'test kernel') % test G_{taux taux tauy nx} +if strcmpi(type, 'test kernel') % test K21 kernel srcnorm = srcinfo.n; srctang = srcinfo.d; targnorm = targinfo.n; targtang = targinfo.d; coefs = varargin{1}; - %R = 0.001; % radius of kernel replacement - [~, ~, forth] = chnk.flex2d.thirdforth_derivatives(zk, src, targ); % Hankel part + [~, ~,~, ~, forth] = chnk.flex2d.helmdiffgreen(zk, src, targ, true); % Hankel part nx = repmat(srcnorm(1,:),nt,1); ny = repmat(srcnorm(2,:),nt,1); @@ -987,7 +799,7 @@ nytarg = repmat((targnorm(2,:)).',1,ns); zkimag = (1i)*zk; - [~, ~, forthK] = chnk.flex2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + [~, ~, ~, ~, forthK] = chnk.flex2d.helmdiffgreen(zkimag, src, targ, true); % modified bessel K part @@ -1012,31 +824,60 @@ rx = targ(1,:).' - src(1,:); ry = targ(2,:).' - src(2,:); r2 = rx.^2 + ry.^2; - dists = sqrt(r2); + tauxtargnsrc = tauxtarg.*nx + tauytarg.*ny; + tauxtargntarg = tauxtarg.*nxtarg + tauytarg.*nytarg; + + normaldot = nx.*nxtarg + ny.*nytarg; + tangentdot = taux.*tauxtarg + tauy.*tauytarg; + + rn = rx.*nx + ry.*ny; + + rntarg = rx.*nxtarg + ry.*nytarg; + + rtau = rx.*taux + ry.*tauy; + + rtautarg = rx.*tauxtarg + ry.*tauytarg; %inds = and((dists <= R),(dists > 0)); - submat = 1/(2*zk^2).*(forth(:, :, 1).*(nxtarg.*nxtarg.*nxtarg.*nx) + forth(:, :, 2).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... - forth(:, :, 4).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) + forth(:, :, 6).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... - forth(:, :, 8).*(nytarg.*nytarg.*nytarg.*ny)) -... + submat = -(1/(2*zk^2).*(forth(:, :, 1).*(nxtarg.*nxtarg.*nxtarg.*nx) + forth(:, :, 2).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... + forth(:, :, 3).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) +... + forth(:, :, 4).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... + forth(:, :, 5).*(nytarg.*nytarg.*nytarg.*ny)) -... 1/(2*zk^2).*(forthK(:, :, 1).*(nxtarg.*nxtarg.*nxtarg.*nx) + forthK(:, :, 2).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... - forthK(:, :, 4).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) + forthK(:, :, 6).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... - forthK(:, :, 8).*(nytarg.*nytarg.*nytarg.*ny)) + (3/(4*pi)).*(nxtarg.*nx + nytarg.*ny)./r2 + ... + forthK(:, :, 3).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) +... + forthK(:, :, 4).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... + forthK(:, :, 5).*(nytarg.*nytarg.*nytarg.*ny))) - ... ((2-coefs(1))/(2*zk^2).*(forth(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg.*nx) + forth(:, :, 2).*(tauxtarg.*tauxtarg.*nxtarg.*ny + tauxtarg.*tauxtarg.*nytarg.*nx + 2*tauxtarg.*tauytarg.*nxtarg.*nx)+... - forth(:, :, 4).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... - forth(:, :, 6).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... - forth(:, :, 8).*(tauytarg.*tauytarg.*nytarg.*ny)) -... + forth(:, :, 3).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... + forth(:, :, 4).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... + forth(:, :, 5).*(tauytarg.*tauytarg.*nytarg.*ny)) -... (2-coefs(1))/(2*zk^2).*(forthK(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg.*nx) + forthK(:, :, 2).*(tauxtarg.*tauxtarg.*nxtarg.*ny + tauxtarg.*tauxtarg.*nytarg.*nx + 2*tauxtarg.*tauytarg.*nxtarg.*nx)+... - forthK(:, :, 4).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... - forthK(:, :, 6).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... - forthK(:, :, 8).*(tauytarg.*tauytarg.*nytarg.*ny)) +... - (2-coefs(1))/(4*pi).*(nxtarg.*nx + nytarg.*ny)./r2 - ... - ((2-coefs(1))/(2*pi)).*(nxtarg.*nx + nytarg.*ny).*(rx.*tauxtarg + ry.*tauytarg).*(rx.*tauxtarg + ry.*tauytarg)./(r2.^2)) - ... + forthK(:, :, 3).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... + forthK(:, :, 4).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... + forthK(:, :, 5).*(tauytarg.*tauytarg.*nytarg.*ny))) + ... + 3/(2*pi).*(rn.*rntarg./(r2.^2)) + 3/(2*pi).*(rntarg.^2).*(normaldot)./(r2.^2) -... + (2/pi).*(rn.*(rntarg.^3))./(r2.^3) + ... + (2-coefs(1)).*(-1/(2*pi).*(tauxtargnsrc).*(tauxtargntarg)./(r2) +... + 1/(pi).*(rn.*rtautarg.*tauxtargntarg)./(r2.^2) + ... + 1/(pi).*(rtautarg.*tauxtargnsrc.*rntarg)./(r2.^2) -.... + 2/(pi).*(rn.*rntarg.*(rtautarg.*rtautarg))./(r2.^3) + ... + 1/(2*pi).*(rn.*rntarg)./(r2.^2)) -... ((1+coefs(1))/(4*pi)).*((taux.*tauxtarg + tauy.*tauytarg)./(r2) - 2*(rx.*tauxtarg + ry.*tauytarg).*(rx.*taux + ry.*tauy)./(r2.^2))-... (3/(4*pi)).*(nxtarg.*nx + nytarg.*ny)./r2 - (2-coefs(1))/(4*pi).*(nxtarg.*nx + nytarg.*ny)./r2 + ... (2-coefs(1))/(2*pi).*(nxtarg.*nx + nytarg.*ny).*(rx.*tauxtarg + ry.*tauytarg).*(rx.*tauxtarg + ry.*tauytarg)./(r2.^2); + %+... + % 3/(2*pi).*(rn.*rntarg./(r2.^2)) + 3/(2*pi).*(rntarg.^2).*(normaldot)./(r2.^2) -... + % (2/pi).*(rn.*(rntarg.^3))./(r2.^3) + ... + % (2-coefs(1)).*(-1/(2*pi).*(tauxtargnsrc).*(tauxtargntarg)./(r2) +... + % 1/(pi).*(rn.*rtautarg.*tauxtargntarg)./(r2.^2) + ... + % 1/(pi).*(rtautarg.*tauxtargnsrc.*rntarg)./(r2.^2) -.... + % 2/(pi).*(rn.*rntarg.*(rtautarg.*rtautarg))./(r2.^3) + ... + % 1/(2*pi).*(rn.*rntarg)./(r2.^2)); % third kernel with singularity subtraction and H[delta'] - - % if any(inds(:)) % if statement for evalauation points that are closed to the boundary + % +... + %(2-coefs(1))/(4*pi).*(nxtarg.*nx + nytarg.*ny)./r2 - ... + % ((2-coefs(1))/(2*pi)).*(nxtarg.*nx + nytarg.*ny).*(rx.*tauxtarg + ry.*tauytarg).*(rx.*tauxtarg + ry.*tauytarg)./(r2.^2) + % if any(inds(:)) % if statement for evalauation points that are closed to the boundary % temp = 1i*imag(submat(inds)); % % nx = nx(inds); ny = ny(inds); @@ -1079,7 +920,7 @@ % % % - % [smoothfirst, smoothsecond] = chnk.flex2d.smooth_part_derivatives(zk, rx, ry); + % [smoothfirst, smoothsecond] = chnk.helm2d.smooth_part_derivatives(zk, rx, ry); % % eulerconstantpart = (smoothfirst(:, :, 1)).*(nxtarg.*nxtarg.*nxtarg.*nx) +... % (smoothfirst(:, :, 2)).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... @@ -1113,7 +954,9 @@ % ((2-coefs(1))/(2*pi)).*((rtautarg.*rtautarg).*(normaldot))./(r2.^2) -... % ((2-coefs(1))/(4*pi)).*(normaldot)./r2 - ... % ((1+coefs(1))/(4*pi)).*((tangentdot)./(r2) - 2*(rtautarg).*(rtau)./(r2.^2)) + temp + eulerconstantpart + puresmoothpart; - %end + % end + + end @@ -1135,21 +978,21 @@ - [~, third, ~] = chnk.flex2d.thirdforth_derivatives(zk, src, targ); % Hankel part + [~, ~, ~, third, ~] = chnk.flex2d.helmdiffgreen(zk, src, targ); % Hankel part zkimag = 1i*zk; - [~, thirdK, ~] = chnk.flex2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + [~,~, ~, thirdK, ~] = chnk.flex2d.helmdiffgreen(zkimag, src, targ); % modified bessel K part - submat = 1/(2*zk^2).*(third(:, :, 1).*(nx.*nx.*nx) + third(:, :, 2).*(3*nx.*nx.*ny) +... - third(:, :, 4).*(3*nx.*ny.*ny) + third(:, :, 6).*(ny.*ny.*ny)) - ... + submat = -(1/(2*zk^2).*(third(:, :, 1).*(nx.*nx.*nx) + third(:, :, 2).*(3*nx.*nx.*ny) +... + third(:, :, 3).*(3*nx.*ny.*ny) + third(:, :, 4).*(ny.*ny.*ny)) - ... 1/(2*zk^2).*(thirdK(:, :, 1).*(nx.*nx.*nx) + thirdK(:, :, 2).*(3*nx.*nx.*ny) +... - thirdK(:, :, 4).*(3*nx.*ny.*ny) + thirdK(:, :, 6).*(ny.*ny.*ny)) + ... + thirdK(:, :, 3).*(3*nx.*ny.*ny) + thirdK(:, :, 4).*(ny.*ny.*ny))) - ... (3/(2*zk^2).*(third(:, :, 1).*(nx.*taux.*taux) + third(:, :, 2).*(2*nx.*taux.*tauy + ny.*taux.*taux) +... - third(:, :, 4).*(nx.*tauy.*tauy + 2*ny.*taux.*tauy) + third(:, :, 6).*(ny.*tauy.*tauy))-... + third(:, :, 3).*(nx.*tauy.*tauy + 2*ny.*taux.*tauy) + third(:, :, 4).*(ny.*tauy.*tauy))-... 3/(2*zk^2).*(thirdK(:, :, 1).*(nx.*taux.*taux) + thirdK(:, :, 2).*(2*nx.*taux.*tauy + ny.*taux.*taux) +... - thirdK(:, :, 4).*(nx.*tauy.*tauy + 2*ny.*taux.*tauy) + thirdK(:, :, 6).*(ny.*tauy.*tauy))); + thirdK(:, :, 3).*(nx.*tauy.*tauy + 2*ny.*taux.*tauy) + thirdK(:, :, 4).*(ny.*tauy.*tauy))); % G_{ny ny ny} + 3G_{ny tauy tauy} end @@ -1167,10 +1010,10 @@ taux = dx./ds; tauy = dy./ds; - [hess, ~, ~] = chnk.flex2d.thirdforth_derivatives(zk, src, targ); % Hankel part + [~, ~, hess] = chnk.flex2d.helmdiffgreen(zk, src, targ); % Hankel part zkimag = 1i*zk; - [hessK, ~, ~] = chnk.flex2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + [~, ~, hessK, ~, ~] = chnk.flex2d.helmdiffgreen(zkimag, src, targ); % modified bessel K part submat = -(1/(2*zk^2).*(hess(:, :, 1).*(nx.*nx) + hess(:, :, 2).*(2*nx.*ny) + hess(:, :, 3).*(ny.*ny))-... 1/(2*zk^2).*(hessK(:, :, 1).*(nx.*nx) + hessK(:, :, 2).*(2*nx.*ny) + hessK(:, :, 3).*(ny.*ny)))+... @@ -1199,7 +1042,7 @@ end -if strcmpi(type, 'free plate eval first hilbert') % G_{tauy} +if strcmpi(type, 'free plate eval first hilbert') % G_{tauy} (supposed to coupled with the hilbert transform) srctang = srcinfo.d; coefs = varargin{1}; @@ -1239,84 +1082,4 @@ - - -% Dirichlet data/potential correpsonding to transmission rep -if strcmpi(type,'trans_rep') - - coef = ones(2,1); - if(nargin == 5); coef = varargin{1}; end; - srcnorm = srcinfo.n; - [submats,grad] = chnk.flex2d.green(zk,src,targ); - nx = repmat(srcnorm(1,:),nt,1); - ny = repmat(srcnorm(2,:),nt,1); - submatd = -(grad(:,:,1).*nx + grad(:,:,2).*ny); - - submat = zeros(1*nt,2*ns); - submat(1:1:1*nt,1:2:2*ns) = coef(1)*submatd; - submat(1:1:1*nt,2:2:2*ns) = coef(2)*submats; -end - -% Neumann data corresponding to transmission rep -if strcmpi(type,'trans_rep_prime') - coef = ones(2,1); - if(nargin == 5); coef = varargin{1}; end; - targnorm = targinfo.n; - srcnorm = srcinfo.n; - - submat = zeros(nt,ns); - - % Get gradient and hessian info - [submats,grad,hess] = chnk.flex2d.green(zk,src,targ); - - nxsrc = repmat(srcnorm(1,:),nt,1); - nysrc = repmat(srcnorm(2,:),nt,1); - nxtarg = repmat((targnorm(1,:)).',1,ns); - nytarg = repmat((targnorm(2,:)).',1,ns); - - % D' - submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ... - + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)... - + hess(:,:,3).*nysrc.*nytarg); - % S' - submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg); - submat = zeros(nt,2*ns) - submat(1:1:1*nt,1:2:2*ns) = coef(1)*submatdp; - submat(1:1:1*nt,2:2:2*ns) = coef(2)*submatsp; -end - - -% Gradient correpsonding to transmission rep -if strcmpi(type,'trans_rep_grad') - coef = ones(2,1); - if(nargin == 5); coef = varargin{1}; end; - - srcnorm = srcinfo.n; - - submat = zeros(nt,ns,6); - % S - [submats,grad,hess] = chnk.flex2d.green(zk,src,targ); - - nxsrc = repmat(srcnorm(1,:),nt,1); - nysrc = repmat(srcnorm(2,:),nt,1); - % D - submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc); - - submat = zeros(2*nt,2*ns); - - submat(1:2:2*nt,1:2:2*ns) = -coef(1)*(hess(:,:,1).*nxsrc + hess(:,:,2).*nysrc); - submat(1:2:2*nt,2:2:2*ns) = coef(2)*grad(:,:,1); - - submat(2:2:2*nt,1:2:2*ns) = -coef(1)*(hess(:,:,2).*nxsrc + hess(:,:,3).*nysrc); - submat(2:2:2*nt,2:2:2*ns) = coef(2)*grad(:,:,2); - - - - - - - - -end - - +end \ No newline at end of file diff --git a/chunkie/+chnk/+flex2d/kernel_replacement_test.m b/chunkie/+chnk/+flex2d/kernel_replacement_test.m deleted file mode 100644 index 993f468..0000000 --- a/chunkie/+chnk/+flex2d/kernel_replacement_test.m +++ /dev/null @@ -1,44 +0,0 @@ -clear -clc - - -zk = 6; -nu = 1/3; - -% t1 = 0; -% t2 = pi/10000; - -coefs = [nu;0]; - -srcinfo = []; - -srcinfo.r = [3; 0]; -srcinfo.n = [1; 0]; -srcinfo.d = [0; 1]; - - - -targinfo = []; - - -targinfo.r = [3; 1/1000]; -targinfo.n = [1; 0]; -targinfo.d = [0; 1]; - -fkern = @(s,t) chnk.helm2d.kern(zk, s, t, 'test kernel', coefs); -ans1 = fkern(srcinfo, targinfo); - -fkern1 = @(s,t) chnk.helm2d.kern(zk, s, t, 'free plate K21 first part', coefs); -ans2 = fkern1(srcinfo, targinfo); - - -abs(ans1-ans2) - - - - - - - - - diff --git a/chunkie/+chnk/+flex2d/smooth_part_derivatives.m b/chunkie/+chnk/+flex2d/smooth_part_derivatives.m deleted file mode 100644 index a041728..0000000 --- a/chunkie/+chnk/+flex2d/smooth_part_derivatives.m +++ /dev/null @@ -1,82 +0,0 @@ -function [firstpart, secondpart] = smooth_part_derivatives(k, rx, ry) -%CHNK.HELM2D.SMOOTH_PART_DERIVATIVES evaluate the derivatives of the smooth -% part of asymptotics of the modified biharmonic's green's function -% firstpart are derivatives of: -% 1/(2k^2)*((-1/(2*pi))*(log(k/2) + C)* \sum_{k = 0}^{5} -% ((-1)^k/(k!)^2)*(k*r/2)^{2k} + (1/(2*pi))*(log(k/2) + C)* \sum_{k = 0}^{5} -% (1/(k!)^2)*(k*r/2)^{2k}). (r = |x-y|) -% -% secondpart are derivatives of: -% 1/(2k^2)*(1/(2*pi)*sum_{k =1}^{5} a_{k} ((-1)^k/(k!)^2)*(k*r/2)^{2k} - -% 1/(2*pi)*sum_{k =1}^{5} a_{k} (1/(k!)^2)*(k*r/2)^{2k}) -% with a_k = sum_{m =1}^{k} (1/m) and r = |x-y|. -% -% Derivatives list in the order of: -% f_{x1x1x1y1}, f_{x1x1x1y2}, f_{x1x1x2y2}, f_{x1x2x2y2}, f_{x2x2x2y2} - -rx2 = rx.*rx; -ry2 = ry.*ry; - -r2 = rx2+ry2; - -r = sqrt(r2); - -[m,n] = size(rx); - -if nargout > 0 - firstpart = zeros(m, n, 5); - - firstpart(:, :, 1) = -(double(eulergamma) + log(k/2)).*(k^4).*(5*rx.^2 + ry.^2)./(64*pi); - - firstpart(:, :, 2) = -(double(eulergamma) + log(k/2)).*(k^4).*(rx.*ry)./(32*pi); - - firstpart(:, :, 3) = -(double(eulergamma) + log(k/2)).*(k^4).*(rx.^2 + ry.^2)./(64*pi); - - firstpart(:, :, 4) = -(double(eulergamma) + log(k/2)).*(k^4).*(rx.*ry)./(32*pi); - - firstpart(:, :, 5) = -(double(eulergamma) + log(k/2)).*(k^4).*(rx.^2 + 5*ry.^2)./(64*pi); -end - -if nargout > 1 - secondpart = zeros(m, n, 5); - - secondpart(:, :, 1) = 1/(4*pi*k^2).*(-(9*k^4)/16 - (11/13824).*(k^6).*(-288.*rx.^2 -72*r2) + ... - (25/1769472)*(k^8).*(-384*rx.^4 - 1152*(rx.^2).*r2 - 144*r2.^2) - ... - (137/884736000)*(k^10).*(-1920*(rx.^4).*r2 - 2880*(rx.^2).*(r2.^2) - 240*(r2.^3))) -... - 1/(4*pi*k^2).*(-(9*k^4)/16 + (11/13824).*(k^6).*(-288.*rx.^2 -72*r2) + ... - (25/1769472)*(k^8).*(-384*rx.^4 - 1152*(rx.^2).*r2 - 144*r2.^2) + ... - (137/884736000)*(k^10).*(-1920*(rx.^4).*r2 - 2880*(rx.^2).*(r2.^2) - 240*(r2.^3))); - - secondpart(:, :, 2) = 1/(4*pi*k^2).*(-(25/1769472)*(k^8).*(-384.*(rx.^3).*ry - 576.*rx.*r2.*ry) +... - (137/884736000)*(k^10).*(-1920*(rx.^3).*r2.*ry - 1440*rx.*(r2.^2).*ry) -... - (11/96)*(k^6).*rx.*ry) +... - 1/(4*pi*k^2).*(-(25/1769472)*(k^8).*(-384.*(rx.^3).*ry - 576.*rx.*r2.*ry) -... - (137/884736000)*(k^10).*(-1920*(rx.^3).*r2.*ry - 1440*rx.*(r2.^2).*ry) +... - (11/96)*(k^6).*rx.*ry); - - secondpart(:, :, 3) = 1/(4*pi*k^2).*(-(25/221184)*(k^8).*rx.*(-72*r2.*ry - 48*(ry.^3)) + ... - (137/884736000)*(k^10).*(rx.*(-144*(r2.^2).*ry - 192*r2.*(ry.^3)))-... - (11/96)*(k^6).*rx.*ry) + ... - 1/(4*pi*k^2).*((25/221184)*(k^8).*rx.*(-72*r2.*ry - 48*(ry.^3)) - ... - (137/884736000)*(k^10).*(rx.*(-144*(r2.^2).*ry - 192*r2.*(ry.^3))) +... - (11/96)*(k^6).*rx.*ry) ; - - secondpart(:, :, 4) = 1/(4*pi*k^2).*(-(25/221184)*(k^8).*rx.*(-72*r2.*ry - 48*(ry.^3)) + ... - (137/884736000)*(k^10).*(rx.*(-144*(r2.^2).*ry - 192*r2.*(ry.^3)))-... - (11/96)*(k^6).*rx.*ry) + ... - 1/(4*pi*k^2).*((25/221184)*(k^8).*rx.*(-72*r2.*ry - 48*(ry.^3)) - ... - (137/884736000)*(k^10).*(rx.*(-144*(r2.^2).*ry - 192*r2.*(ry.^3))) +... - (11/96)*(k^6).*rx.*ry) ; - - secondpart(:, :, 5) = 1/(4*pi*k^2).*(-(9*k^4)/16 - (11/13824).*(k^6).*(-288.*ry.^2 -72*r2) + ... - (25/1769472)*(k^8).*(-384*ry.^4 - 1152*(ry.^2).*r2 - 144*r2.^2) - ... - (137/884736000)*(k^10).*(-1920*(ry.^4).*r2 - 2880*(ry.^2).*(r2.^2) - 240*(r2.^3))) -... - 1/(4*pi*k^2).*(-(9*k^4)/16 + (11/13824).*(k^6).*(-288.*ry.^2 -72*r2) + ... - (25/1769472)*(k^8).*(-384*ry.^4 - 1152*(ry.^2).*r2 - 144*r2.^2) + ... - (137/884736000)*(k^10).*(-1920*(ry.^4).*r2 - 2880*(ry.^2).*(r2.^2) - 240*(r2.^3))); - -end - - - -end \ No newline at end of file diff --git a/chunkie/+chnk/+flex2d/test_ps_besseldiff.m b/chunkie/+chnk/+flex2d/test_ps_besseldiff.m deleted file mode 100644 index bec9af0..0000000 --- a/chunkie/+chnk/+flex2d/test_ps_besseldiff.m +++ /dev/null @@ -1,16 +0,0 @@ - -% test bessel0-1 and other ps coefficients - -clearvars; clc - -gam = 0.57721566490153286060651209; - -nterms = 12; -[cf1,cf2] = besseldiff_etc_pscoefs(nterms); - -z = 1.3 + 0.1*1i; -v1 = even_pseval(cf1,z); -abs(v1-(besselj(0,z)-1)) - -v2 = even_pseval(cf2,z); -abs(besselh(0,z)-besselj(0,z)-1i*(2/pi)*((log(z/2)+gam)*besselj(0,z) + v2)) diff --git a/chunkie/+chnk/+flex2d/transmission_helper.m b/chunkie/+chnk/+flex2d/transmission_helper.m deleted file mode 100644 index b199e36..0000000 --- a/chunkie/+chnk/+flex2d/transmission_helper.m +++ /dev/null @@ -1,227 +0,0 @@ -function [kerns,varargout] = transmission_helper(chnkobj,ks,cs,coefs,varargin) -%CHNK.HELM2D.TRANSMISSION_HELPER builds the matrix of kernels required to call -% chunkermat and generate boundary data which is appropriately scaled for the linear -% system if requested. -% -% Syntax: kerns = transmission_helper(chnkrs,ks,cs,coefs), returns a matrix of kernels -% [kerns,bdrydata] = transmission_helper(chnkrs,ks,cs,coefs,opts), returns -% a matrix of kernels and the boundary data due to point sources or -% waves as prescribed by opts -% -% Note: Here we always assume that the transmission problem is either in -% the 'TE' or 'TM' polarization, which implies that the boundary conditions -% are of the form -% [u] = f, [coef du/dn] = g -% -% where [u] denotes the jump in u across an interface. -% To obtain the TE polarization, coef \equiv 1, and for the TM polarization -% coef should be set to the square of the refractive index of the medium. -% -% These set of transmission problems do not cover all possibilities, and -% in particular exclude transmission problems with boundary conditions -% of the form -% [\beta u] = f. -% -% Input: -% chnkrs - cell array of chunker objects describing the geometry (1,ncurve) -% ks - wavenumbers in different regions of teh plane -% cs - cs(1,i), cs(2,i) denote the region number in the direction of the normal -% and opposite of the direction of the normal for curve chnkrs(i) -% coefs - denotes the scaling parameter for the jump in the neumann data, depends -% only on the region -% opts - optional structure for defining the boundary condition -% opts.bdry_data_type - 'pw' or 'point sources' -% opts.bdry_data_type = pw, then opts.dir should be set to incident -% angle, and opts.exposed_curves should be set to the curves for which -% incident field is to be set -% if opts.bdry_data_type = 'point sources' then opts.sources should be a cell -% array of size (1,ncurve), with opts.sources{i}.locs (2,m) locations -% of charges which generate the data for region i, and .charges are the -% corresponding charge strengths -% -% - - - if nargin == 4 - opts = []; - - elseif nargin == 5 - opts = varargin{1}; - else - fprintf('Invalid number of input arguments\n'); - fprintf('Returning without computing anything\n'); - % this behavior might need to be fixed - kerns = []; - varargout{2:nargout} = []; - return; - end - - - if (class(chnkobj) == "chunker") - chnkrs = chnkobj; - elseif(class(chnkobj) == "chunkgraph") - chnkrs = chnkobj.echnks; - else - msg = "Unsupported object in chunkermat"; - error(msg) - end - - - bdry_data_type = 'pw'; - direction = 0; - ncurve = length(chnkrs); - sources = cell(1,ncurve); - charges = cell(1,ncurve); - exposed_curves = ones(1,ncurve); - for i=1:ncurve - sources{i} = zeros(2,1); - charges{i} = 0 + 1j*0; - end - - if(isfield(opts,'bdry_data_type')) - bdry_data_type = opts.bdry_data_type; - end - if(isfield(opts,'dir')) - direction = opts.dir; - end - if(isfield(opts,'sources')) - sources = opts.sources; - end - if(isfield(opts,'charges')) - charges = opts.charges; - end - if(isfield(opts,'exposed_curves')) - exposed_curves = opts.exposed_curves; - end - - - - - kerns = cell(ncurve,ncurve); - - k1 = ks(cs(1,:)); - k2 = ks(cs(2,:)); - - alpha1 = 2./(coefs(cs(1,:)) + coefs(cs(2,:))); - alpha2 = 2./(1./coefs(cs(1,:)) + 1./coefs(cs(2,:))); - - -% First build system matrix without any corner corrections - cc1 = zeros(2,2); - cc2 = zeros(2,2); - for i=1:ncurve - - c1 = coefs(cs(1,i)); - c2 = coefs(cs(2,i)); - cc1(1,1) = -alpha1(i)*c1; - cc2(1,1) = -alpha1(i)*c2; - - cc1(1,2) = -alpha1(i); - cc2(1,2) = -alpha1(i); - - cc1(2,1) = alpha2(i); - cc2(2,1) = alpha2(i); - - cc1(2,2) = alpha2(i)/c1; - cc2(2,2) = alpha2(i)/c2; - - for j=1:ncurve - kerns{i,j} = @(s,t) -(chnk.helm2d.kern(k1(i),s,t,'all',cc1)- ... - chnk.helm2d.kern(k2(i),s,t,'all',cc2)); - end - end - - - - nchs = zeros(1,ncurve); - for i=1:ncurve - nchs(i) = chnkrs(i).nch; - end - - ngl = chnkrs(1).k; - np = sum(nchs(1:ncurve))*ngl; - - - if nargout > 1 - bdry_data = complex(zeros(2*np,1)); - if(strcmpi(bdry_data_type,'pw')) - for i=1:ncurve - d1 = cs(1,i); - d2 = cs(2,i); - c1 = coefs(d1); - c2 = coefs(d2); - ind1 = sum(nchs(1:i-1))*ngl*2+(1:2:2*nchs(i)*ngl); - ind2 = sum(nchs(1:i-1))*ngl*2+(2:2:2*nchs(i)*ngl); - targnorm = chnkrs(i).n; - nx = targnorm(1,:); nx = nx(:); - ny = targnorm(2,:); ny = ny(:); - if(exposed_curves(i)) - - if(exposed_curves(i) < 0) - k = k2(i); - a1 = alpha1(i); - a2 = -alpha2(i)/c2; - else - k = k1(i); - a1 = -alpha1(i); - a2 = alpha2(i)/c1; - - end - x = chnkrs(i).r(1,:); - x = x(:); - y = chnkrs(i).r(2,:); - y = y(:); - ct = cos(direction); - st = sin(direction); - u = exp(1i*k*(x*ct + y*st)); - dudn = 1i*k*(ct*nx+st*ny).*u; - bdry_data(ind1) = a1*u; - bdry_data(ind2) = a2*dudn; - - end - end - - elseif(strcmpi(bdry_data_type,'point sources')) - - for i=1:ncurve - d1 = cs(1,i); - d2 = cs(2,i); - c1 = coefs(d1); - c2 = coefs(d2); - - src1 = sources{d1}; - charges1 = charges{d1}; - charges1 = charges1(:); - src2 = sources{d2}; - charges2 = charges{d2}; - charges2 = charges2(:); - - ind1 = sum(nchs(1:i-1))*ngl*2+(1:2:2*nchs(i)*ngl); - ind2 = sum(nchs(1:i-1))*ngl*2+(2:2:2*nchs(i)*ngl); - - targnorm = chnkrs(i).n; - nx = targnorm(1,:); nx = nx.'; - ny = targnorm(2,:); ny = ny.'; - [u1,gradu1] = chnk.helm2d.green(k1(i),src1,chnkrs(i).r(:,:)); - u1 = u1*charges1; - gradu1(:,:,1) = gradu1(:,:,1)*charges1; - gradu1(:,:,2) = gradu1(:,:,2)*charges1; - dudn1 = squeeze(gradu1(:,1)).*nx + squeeze(gradu1(:,2)).*ny; - - [u2,gradu2] = chnk.helm2d.green(k2(i),src2,chnkrs(i).r(:,:)); - u2 = u2*charges2; - gradu2(:,:,1) = gradu2(:,:,1)*charges2; - gradu2(:,:,2) = gradu2(:,:,2)*charges2; - dudn2 = squeeze(gradu2(:,1)).*nx + squeeze(gradu2(:,2)).*ny; - - bdry_data(ind1) = alpha1(i)*(u1-u2); - bdry_data(ind2) = -alpha2(i)*(1/c1*dudn1 - 1/c2*dudn2); - end - - end - varargout{1} = bdry_data; - - end - - -end