Farfield flex

chunkie/+chnk/+flex2d/hkdiffgreen_ff.m

@@ -0,0 +1,192 @@
+function [val,grad,hess,der3,der4,der5] = hkdiffgreen_ff(k,src,targ,opts)
+
+prefac = 1i/(2*sqrt(2*pi))*exp(-1i*pi/4);
+
+[~,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);
+rt = sqrt(xt.^2+yt.^2);
+
+xt = xt./rt;
+yt = yt./rt;
+
+dp = xt.*xs+yt.*ys;
+
+g0 = prefac*exp(-1i*dp*k);
+
+
+%     evaluate potential and derivatives
+
+if nargout > 0
+    val = g0;  
+end
+if nargout > 1
+    grad(:,:,1) = -1i*k*g0.*xt;
+    grad(:,:,2) = -1i*k*g0.*yt;
+    grad = -grad;
+end
+if nargout > 2
+    hess(:,:,1) = (-1i*k)^2*g0.*xt.^2;
+    hess(:,:,2) = (-1i*k)^2*g0.*xt.*yt;
+    hess(:,:,3) = (-1i*k)^2*g0.*yt.^2;
+end
+if nargout > 3
+    der3(:,:,1) = (-1i*k)^3*g0.*(xt.^3);
+    der3(:,:,2) = (-1i*k)^3*g0.*(xt.^2).*yt;
+    der3(:,:,3) = (-1i*k)^3*g0.*(yt.^2).*xt;
+    der3(:,:,4) = (-1i*k)^3*g0.*(yt.^3);
+    der3 = -der3;
+end
+
+if nargout > 4
+    der4(:,:,1) = (-1i*k)^4*g0.*(xt.^4);
+    der4(:,:,2) = (-1i*k)^4*g0.*(xt.^3).*(yt);
+    der4(:,:,3) = (-1i*k)^4*g0.*(xt.^2).*(yt.^2);
+    der4(:,:,4) = (-1i*k)^4*g0.*(xt).*(yt.^3);
+    der4(:,:,5) = (-1i*k)^4*g0.*(yt.^4);
+end
+
+if nargout > 5
+    der5(:,:,1) = (-1i*k)^5*g0.*(xt.^5);
+    der5(:,:,2) = (-1i*k)^5*g0.*(xt.^4).*(yt);
+    der5(:,:,3) = (-1i*k)^5*g0.*(xt.^3).*(yt.^2);
+    der5(:,:,4) = (-1i*k)^5*g0.*(xt.^2).*(yt.^3);
+    der5(:,:,5) = (-1i*k)^5*g0.*(xt.^1).*(yt.^4);
+    der5(:,:,6) = (-1i*k)^5*g0.*(yt.^5);
+    der5 = -der5;
+end
+
+end
+
+function [g0,g1,g21,g321,g4321,g54321] = diff_h0k0_and_rders(k,r,r2logrfac)
+% g0 = g
+% g1 = g'
+% g21 = g'' - g'/r
+% g321 = g''' - 3*g''/r + 3g'/r^2 = g''' - 3*g21/r
+% g4321 = g'''' - 6*g'''/r + 15*g''/r^2 - 15*g'/r^3 
+%       = g''''- 6 g321/r - 3 g21/r^2
+% g54321 = g''''' - 10g''''/r + 45g'''/r^2 -105g''/r^3 + 105g'/r^4
+%        = g5 - 10g4321/r - 15 g321/r^2
+
+g0 = zeros(size(r));
+g1 = zeros(size(r));
+g321 = zeros(size(r));
+g4321 = zeros(size(r));
+g54321 = 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
+
+rnot = r(~isus);
+
+kr = k*rnot;
+
+h0 = besselh(0,1,kr);
+h1 = besselh(1,1,kr);
+h0i = besselh(0,1,1i*kr);
+h1i = besselh(1,1,1i*kr);
+
+rm1 = 1./rnot;
+rm2 = rm1.*rm1;
+rm3 = rm1.*rm2;
+rm4 = rm1.*rm3;
+
+r2fac = (1-r2logrfac)*k*k*0.5*o2p;
+logr = log(rnot);
+g0(~isus) = io4*(h0 - h0i) - r2fac*rnot.*rnot.*logr;
+g1(~isus) = io4*(-k*h1 + 1i*k*h1i) - r2fac*(rnot+2*rnot.*logr);
+g21(~isus) = io4*(-k*k*h0 + k*h1.*rm1 - k*k*h0i - 1i*k*h1i.*rm1) ...
+    - r2fac*(3+2*logr) - g1(~isus).*rm1;
+g321(~isus) = io4*(k*k*h0.*rm1 + k*(k*k-2*rm2).*h1 ...
+    + k*k*h0i.*rm1 - 1i*k*(-k*k-2*rm2).*h1i) - r2fac*2*rm1 - 3*g21(~isus).*rm1;
+g4321(~isus) = io4*(k*(3*rm2-k*k).*(2*h1.*rm1-k*h0) - ...
+    1i*k*(3*rm2+k*k).*(2*h1i.*rm1-1i*k*h0i)) + r2fac*2*rm2 - 6*g321(~isus).*rm1 ...
+    - 3*g21(~isus).*rm2;
+g54321(~isus) = io4*(k*( (12*k*rm3-2*k^3*rm1).*h0  + ...
+    (-24*rm4+7*k*k*rm2-k^4).*h1) - 1i*k*( (12*1i*k*rm3+1i*2*k^3*rm1).*h0i + ...
+    (-24*rm4-7*k*k*rm2-k^4).*h1i)) - r2fac*4*rm3 - 10*g4321(~isus).*rm1 - ...
+    15*g321(~isus).*rm2;
+
+% manually cancel when small
+
+rsus = r(isus);
+rm1 = 1./rsus;
+rm2 = rm1.*rm1;
+rm3 = rm1.*rm2;
+rm4 = rm1.*rm3;
+rm5 = rm1.*rm4;
+
+r2 = rsus.*rsus; 
+r4 = r2.*r2;
+
+nterms = 7;
+
+% relevant parts of hankel function represented as power series
+[cf1,cf2] = chnk.flex2d.besselkikdiff_etc_pscoefs(nterms);
+cf1(1) = cf1(1)*r2logrfac;
+kpow = k*k*(k.^(4*(0:(nterms-1)))); % k^2, k^6, k^10, ...
+cf1 = cf1(:).*kpow(:); cf2 = cf2(:).*kpow(:);
+
+
+jdiff = chnk.flex2d.pseval(cf1,r4).*r2;
+f = chnk.flex2d.pseval(cf2,r4).*r2;
+
+% differentiate power series to get derivatives
+fac = 2+4*(0:(nterms-1)).';
+jdiffd1 = chnk.flex2d.pseval(cf1.*fac,r4).*rsus;
+fd1 = chnk.flex2d.pseval(cf2.*fac,r4).*rsus;
+d = fac.*(fac-1)-fac;
+jdiffd21 = chnk.flex2d.pseval(cf1.*d,r4);
+fd21 = chnk.flex2d.pseval(cf2.*d,r4);
+% third deriv and higher the first term is dead (this is subtle but true
+% for the 21, 321, 4321, etc derivative)
+cf1 = cf1(2:end); cf2 = cf2(2:end); fac = fac(2:end);
+d = fac.*(fac-1).*(fac-2) - 3*(fac.*(fac-1)-fac);
+jdiffd321 = chnk.flex2d.pseval(cf1.*d,r4).*rsus.*r2;
+fd321 = chnk.flex2d.pseval(cf2.*d,r4).*rsus.*r2;
+% g4321 = g'''' - 6*g'''/r + 15*g''/r^2 - 15*g'/r^3 
+d = fac.*(fac-1).*(fac-2).*(fac-3)-6*fac.*(fac-1).*(fac-2)+15*(fac.*(fac-1)-fac);
+jdiffd4321 = chnk.flex2d.pseval(cf1.*d,r4).*r2;
+fd4321 = chnk.flex2d.pseval(cf2.*d,r4).*r2;
+% g54321 = g''''' - 10g''''/r + 45g'''/r^2 -105g''/r^3 + 105g'/r^4
+d = fac.*(fac-1).*(fac-2).*(fac-3).*(fac-4) - ...
+    10*fac.*(fac-1).*(fac-2).*(fac-3) + 45*fac.*(fac-1).*(fac-2) - ...
+    105*(fac.*(fac-1)-fac);
+jdiffd54321 = chnk.flex2d.pseval(cf1.*d,r4).*rsus;
+fd54321 = chnk.flex2d.pseval(cf2.*d,r4).*rsus;
+
+% combine to get derivative of i/4 H - K/2pi 
+gam = 0.57721566490153286060651209;
+logr = log(rsus) + log(k) - log(2) + gam;
+const1 = (1-r2logrfac)*o2p*(log(k)+gam-log(2))*k*k/2;
+
+j0 = besselj(0,k*rsus);
+j1 = besselj(1,k*rsus);
+j2 = besselj(2,k*rsus);
+j3 = besselj(3,k*rsus);
+j4 = besselj(4,k*rsus);
+j5 = besselj(5,k*rsus);
+g0(isus) = io4*j0 - o2p*(f  + logr.*jdiff) + const1*r2;
+g1(isus) = io4*(-k*j1) - o2p*(fd1 + logr.*jdiffd1 + jdiff.*rm1) + 2*const1*rsus;
+g21(isus) = io4*(k*k*j2) - o2p*(fd21 + logr.*jdiffd21 + 2*jdiffd1.*rm1 - 2*jdiff.*rm2);
+g321(isus) = io4*(-k^3*j3) - o2p*(fd321 + logr.*jdiffd321 + 3*jdiffd21.*rm1 - ...
+    6*jdiffd1.*rm2 + 8*jdiff.*rm3);
+g4321(isus) = io4*(k^4*j4)- o2p*(fd4321 + logr.*jdiffd4321 + 4*jdiffd321.*rm1 - ...
+    12*jdiffd21.*rm2 + 32*jdiffd1.*rm3 - 48*jdiff.*rm4);
+g54321(isus) = io4*(-k^5*j5) - o2p*(fd54321 + logr.*jdiffd54321 + ...
+    5*jdiffd4321.*rm1 - 20*jdiffd321.*rm2 + 80*jdiffd21.*rm3 - 240*jdiffd1.*rm4 + 384*jdiff.*rm5);
+
+end
+

chunkie/+chnk/+flex2d/kern.m

@@ -221,6 +221,34 @@
     submat(:,1:2:end) = K1;
     submat(:,2:2:end) = K2;
 
+% clamped plate kernels for far field evaluation
+case {'clamped_plate_eval_ff'}
+
+    submat = zeros(nt,2*ns);
+
+    srcnorm = srcinfo.n;
+    srctang = srcinfo.d;
+    nx = repmat(srcnorm(1,:),nt,1);
+    ny = repmat(srcnorm(2,:),nt,1);
+    dx = repmat(srctang(1,:),nt,1);
+    dy = repmat(srctang(2,:),nt,1);
+    ds = sqrt(dx.*dx+dy.*dy);
+
+    taux = dx./ds;
+    tauy = dy./ds;
+
+    [~, ~, hess, third] = chnk.flex2d.hkdiffgreen_ff(zk, src, targ);           % Hankel part
+
+    K1 = -(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)) ) - ...
+       (3/(2*zk^2).*(third(:, :, 1).*(nx.*taux.*taux) + third(:, :, 2).*(2*nx.*taux.*tauy + ny.*taux.*taux) +...
+       third(:, :, 3).*(nx.*tauy.*tauy + 2*ny.*taux.*tauy) + third(:, :, 4).*(ny.*tauy.*tauy)));  % G_{ny ny ny} + 3G_{ny tauy tauy}
+
+    K2 =  -(1/(2*zk^2).*(hess(:, :, 1).*(nx.*nx) + hess(:, :, 2).*(2*nx.*ny) + hess(:, :, 3).*(ny.*ny)))+...
+          (1/(2*zk^2).*(hess(:, :, 1).*(taux.*taux) + hess(:, :, 2).*(2*taux.*tauy) + hess(:, :, 3).*(tauy.*tauy))); % -G_{ny ny}  + G_{tauy tauy}
+
+    submat(:,1:2:end) = K1;
+    submat(:,2:2:end) = K2;
 
 %%% FREE PLATE KERNELS
 
@@ -412,6 +440,33 @@
    submat(:,2:3:end) = K1H;
    submat(:,3:3:end) = K2;
 
+% free plate kernels for far field evaluation
+case {'free_plate_eval_ff'}           
+   srcnorm = srcinfo.n;
+   srctang = srcinfo.d;
+   nu = varargin{1};
+
+   [val,grad] = chnk.flex2d.hkdiffgreen_ff(zk,src,targ); 
+   nx = repmat(srcnorm(1,:),nt,1);
+   ny = repmat(srcnorm(2,:),nt,1);
+
+   dx = repmat(srctang(1,:),nt,1);
+   dy = repmat(srctang(2,:),nt,1);
+
+   ds = sqrt(dx.*dx+dy.*dy);
+
+   taux = dx./ds; 
+   tauy = dy./ds;
+
+   K1 = (-1/(2*zk^2).*(grad(:, :, 1).*(nx) + grad(:, :, 2).*ny)); 
+   K1H = ((1 + nu)/2).*(-1/(2*zk^2).*(grad(:, :, 1).*(taux) + grad(:, :, 2).*tauy));                    % G_{tauy}
+   K2 = 1/(2*zk^2).*val;
+
+   submat = zeros(nt,3*ns);
+   submat(:,1:3:end) = K1;
+   submat(:,2:3:end) = K1H;
+   submat(:,3:3:end) = K2;
+
 %%% SUPPORTED PLATE KERNELS
 
 % boundary conditions applied to a point source
@@ -707,6 +762,57 @@
     submat(:,1:2:end) = K1;
     submat(:,2:2:end) = K2;
 
+
+% supported plate kernels for far field evaluation
+case{'supported_plate_eval_ff'}
+    srcnorm = srcinfo.n;
+    srctang = srcinfo.d;
+    srcd2 = srcinfo.d2;
+    coefs = varargin{1};
+    nu = coefs(1);
+
+    nx = repmat(srcnorm(1,:),nt,1);
+    ny = repmat(srcnorm(2,:),nt,1);
+
+    dx = repmat(srctang(1,:),nt,1);
+    dy = repmat(srctang(2,:),nt,1);
+
+    ds = sqrt(dx.*dx+dy.*dy);
+
+    taux = dx./ds;
+    tauy = dy./ds;
+
+    dx1 = repmat(srctang(1,:),nt,1);
+    dy1 = repmat(srctang(2,:),nt,1);
+
+    d2x1 = repmat(srcd2(1,:),nt,1);
+    d2y1 = repmat(srcd2(2,:),nt,1);
+
+    denom = sqrt(dx1.^2+dy1.^2).^3;
+    numer = dx1.*d2y1-d2x1.*dy1;
+
+    kappa = numer./denom; 
+
+    kp = repmat(srcinfo.data(1,:),nt,1);
+
+    a1 = 2-nu;
+    a2 = (-1+nu)*(7+nu)/(3 - nu);
+    a3 = (1-nu)*(3+nu)/(1+nu);
+
+    [~, grad, hess, third] = chnk.flex2d.hkdiffgreen_ff(zk, src, targ, false); 
+
+    K1 = -1/(2*zk^2)*(third(:,:,1).*nx.^3 + 3*third(:,:,2).*nx.^2.*ny + 3*third(:,:,3).*nx.*ny.^2 + third(:,:,4).*ny.^3) + ...
+         -a1/(2*zk^2)*(third(:,:,1).*nx.*taux.^2 + third(:,:,2).*(ny.*taux.^2 + 2*nx.*taux.*tauy) + third(:,:,3).*(nx.*tauy.^2 + 2*ny.*taux.*tauy) + third(:,:,4).*ny.*tauy.^2) + ...
+         a2*kappa./(2*zk^2).*(hess(:,:,1).*nx.^2 + 2*hess(:,:,2).*nx.*ny + hess(:,:,3).*ny.^2) + ...
+         -a3*kp./(2*zk^2).*(grad(:,:,1).*taux + grad(:,:,2).*tauy);
+
+    K2 = -1/(2*zk^2).*(grad(:,:,1).*nx + grad(:,:,2).*ny);
+
+    submat = zeros(nt,2*ns);
+
+    submat(:,1:2:end) = K1;
+    submat(:,2:2:end) = K2;
+
 end
 
 end