update in axissymlap2d greens function to include the hessians

chunkie/+chnk/+axissymhelm2d/axissymhelm2d.m

@@ -0,0 +1,108 @@
+function obj = axissymhelm2d(type, zk, coefs)
+%KERNEL.AXISSYMHELM2D   Construct the axissymmetric Helmholtz kernel.
+%   KERNEL.AXISSYMHELM2D('s', ZK) or KERNEL.AXISSYMHELM2D('single', ZK) 
+%   constructs the axissymmetric single-layer Helmholtz kernel with 
+%   wavenumber ZK.
+%
+%   KERNEL.AXISSYMHELM2D('d', ZK) or KERNEL.AXISSYMHELM2D('double', ZK) 
+%   constructs the axissymmetric double-layer Helmholtz kernel with 
+%   wavenumber ZK.
+%
+%   KERNEL.AXISSYMHELM2D('sp', ZK) or KERNEL.AXISSYMHELM2D('sprime', ZK) 
+%   constructs the normal derivative of the axissymmetric single-layer 
+%   Helmholtz kernel with wavenumber ZK.
+%
+%   KERNEL.AXISSYMHELM2D('c', ZK, COEFS) or 
+%   KERNEL.AXISSYMHELM2D('combined', ZK, COEFS)
+%   constructs the combined-layer axissymmetric Helmholtz kernel with 
+%   wavenumber ZK and parameter COEFS, 
+%   i.e., COEFS(1)*KERNEL.AXISSYMHELM2D('d', ZK) + 
+%         COEFS(2)*KERNEL.AXISSYMHELM2D('s', ZK).
+%
+% NOTES: The axissymetric kernels are currently supported only for purely
+% real or purely imaginary wave numbers
+%
+% See also CHNK.AXISSYMHELM2D.KERN.
+
+if ( nargin < 1 )
+    error('Missing Helmholtz kernel type.');
+end
+
+if ( nargin < 2 )
+    error('Missing Helmholtz wavenumber.');
+end
+
+
+zr = real(zk); zi = imag(zk);
+if abs(zr)*abs(zi) > eps
+    error('Only purely real or purely imaginary wavenumbers supported');
+end
+
+obj = kernel();
+obj.name = 'axissymhelmholtz';
+obj.params.zk = zk;
+obj.opdims = [1 1];
+
+switch lower(type)
+
+    case {'s', 'single'}
+        obj.type = 's';
+        obj.eval = @(s,t) chnk.axissymhelm2d.kern(zk, s, t, [0,0], 's');
+        obj.shifted_eval = @(s,t,o) chnk.axissymhelm2d.kern(zk, s, t, o, 's');
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    case {'d', 'double'}
+        obj.type = 'd';
+        obj.eval = @(s,t) chnk.axissymhelm2d.kern(zk, s, t, [0,0], 'd');
+        obj.shifted_eval = @(s,t,o) chnk.axissymhelm2d.kern(zk, s, t, o, 'd');
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    case {'sp', 'sprime'}
+        obj.type = 'sp';
+        obj.eval = @(s,t) chnk.axissymhelm2d.kern(zk, s, t, [0,0], 'sprime');
+        obj.shifted_eval = @(s,t,o) chnk.axissymhelm2d.kern(zk, s, t, o, 'sprime');
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    case {'c', 'combined'}
+        if ( nargin < 3 )
+            warning('Missing combined layer coefficients. Defaulting to [1,1i].');
+            coefs = [1,1i];
+        end
+        obj.type = 'c';
+        obj.params.coefs = coefs;
+        obj.eval = @(s,t) chnk.axissymhelm2d.kern(zk, s, t, [0,0], 'c', coefs);
+        obj.shifted_eval = @(s,t,o) chnk.axissymhelm2d.kern(zk, s, t, o, 'c', coefs);
+        obj.fmm = [];
+        obj.sing = 'log';
+     case {'neu_rpcomb'}
+        obj.type = 'neu_rpcomb';
+        if ( nargin < 3 )
+            warning('Missing coefficient of 1i*D. Defaulting to 1.');
+            coefs = 1;
+        end
+        obj.eval = @(s,t) chnk.axissymhelm2d.kern(zk, s, t, [0,0], 'neu_rpcomb', coefs);
+        obj.shifted_eval = @(s,t,o) chnk.axissymhelm2d.kern(zk, s, t, o, 'neu_rpcomb', coefs);
+        obj.fmm = [];
+        obj.sing = 'log';
+        obj.opdims = [3 3];
+        obj.params.c1 = -1.0/(0.5 + 0.25*1i*coefs);
+        obj.params.c2 = -1i*coefs/(0.5 + 0.25*1i*coefs);
+
+    case {'dp_diff_d0p'}
+        obj.type = 'dp_diff_d0p';
+        obj.eval = @(s,t) chnk.axissymhelm2d.kern(zk, s, t, [0,0], 'dprime') - ...
+                            chnk.axissymlap2d.kern(s, t, [0,0], 'dprime', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymhelm2d.kern(zk, s, t, o, 'sprime') - ...
+                            chnk.axissymlap2d.kern(s, t, o, 'dprime', n);
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    otherwise
+        error('Unknown axissym Helmholtz kernel type ''%s''.', type);
+
+end
+
+end

chunkie/+chnk/+axissymhelm2d/kern.m

@@ -299,6 +299,56 @@
   submat(1:3:end, 3:3:end) = c2*spikmat;
   submat(2:3:end, 1:3:end) = -sikmat;
   submat(3:3:end, 1:3:end) = -spikmat;
+
+case{'spp','sprimeprime'}
+  targnorm = targinfo.n(:,:);
+  [~,~,hess] = chnk.axissymhelm2d.green(zk, src, targ, origin);
+  nxtarg = repmat((targnorm(1,:)).',1,ns);
+  nytarg = repmat((targnorm(2,:)).',1,ns);
+  submat = hess(:,:,1).*nxtarg.*nxtarg + 2*hess(:,:,5).*nytarg.*nxtarg ...
+                                        + hess(:,:,3).*nytarg.*nytarg;
+
+  fker = @(x, s, t, rnt) fsprimeprime(x, zk, s, t, rnt, origin);
+    for j=1:ns
+        for i=1:nt
+            rt = targ(1,i) + origin(1);
+            dr = (src(1,j) - targ(1,i));
+            dz = (src(2,j) - targ(2,i));
+            r0   = sqrt(rt^2+(rt+dr)^2+dz^2);
+            alph = (dr^2+dz^2)/r0^2;
+            if alph > 2e-4 && alph < 0.2
+                [x0, w0] = get_grid(zk, rt, dr, dz);
+                fvals = fker(x0, src(:, j), targ(:,i), targnorm(:,i));
+                submat(i,j) = 2*w0.'*fvals;
+            end
+        end
+    end
+
+case {'q','quad','quadruple','quadrupole'} % q := spp'
+  srcnorm = srcinfo.n(:,:);
+  [~,~,hess] = chnk.axissymhelm2d.green(zk, src, targ, origin);
+  nxsrc = repmat(srcnorm(1,:),nt,1);
+  nysrc = repmat(srcnorm(2,:),nt,1);
+  submat = hess(:,:,2).*nxsrc.*nxsrc - 2*hess(:,:,6).*nysrc.*nxsrc ... 
+                                      + hess(:,:,3).*nysrc.*nysrc;
+
+  fker = @(x, s, t, rns) fqlp(x, zk, s, t, rns, origin);
+    for j=1:ns
+        for i=1:nt
+            rt = targ(1,i) + origin(1);
+            dr = (src(1,j) - targ(1,i));
+            dz = (src(2,j) - targ(2,i));
+            r0   = sqrt(rt^2+(rt+dr)^2+dz^2);
+            alph = (dr^2+dz^2)/r0^2;
+            if alph > 2e-4 && alph < 0.2
+                [x0, w0] = get_grid(zk, rt, dr, dz);
+                fvals = fker(x0, src(:, j), targ(:,i), srcnorm(:,i));
+                submat(i,j) = 2*w0.'*fvals;
+            end
+        end
+    end
+
+
 otherwise
     error('Unknown axissymmetric Helmholtz kernel type ''%s''.', type); 
 end
@@ -357,6 +407,36 @@
            rndt.*rnds.*(-zk^2.*r.^2 - 3*1j*zk.*r + 3)./r.^5).*exp(1j*zk*r)/4/pi.*(rs + o(1));
 end
 
+function f = fqlp (x, zk, s, t, rns, o)
+    rs = s(1); zs = s(2);
+    rt = t(1); zt = t(2);
+
+    sxhalf = sin(x/2);
+    sxhalf2 = sxhalf.*sxhalf;
+    cx = 1-2*sxhalf2;
+
+    rnds = ((rt + o(1)).*cx  - (rs + o(1))).*rns(1) + (zt - zs).*rns(2);
+
+    r = sqrt((rs-rt).^2 + (zs-zt).^2 + 4*(rs+o(1)).*(rt+o(1)).*sxhalf2);
+    f = ((1j*zk.*r-1)./r.^3 + ...
+           rnds.^2.*(-zk.*r.^2 - 3*1j*zk.*r + 3)./r.^5).*exp(1j*zk*r)/4/pi.*(rs+o(1));
+end
+
+function f = fsprimeprime (x, zk, s, t, rnt, o)
+    rs = s(1); zs = s(2);
+    rt = t(1); zt = t(2);
+
+    sxhalf = sin(x/2);
+    sxhalf2 = sxhalf.*sxhalf;
+    cx = 1-2*sxhalf2;
+
+    rndt = ((rt + o(1)) - (rs + o(1)).*cx).*rnt(1) + (zt - zs).*rnt(2);
+
+    r = sqrt((rs-rt).^2 + (zs-zt).^2 + 4*(rs+o(1)).*(rt+o(1)).*sxhalf2);
+    f = ((1j*zk.*r-1)./r.^3 + ...
+           rndt.^2.*(-zk.*r.^2 - 3*1j*zk.*r + 3)./r.^5).*exp(1j*zk*r)/4/pi.*(rs+o(1));
+end
+
 function [fkp, fik, fikp, fkdiff] = get_neu_kers(zk, cx, sxhalf2, s, t, rns, rnt, o)
 
     rs = s(1); zs = s(2);

chunkie/+chnk/+axissymlap2d/axissymlap2d.m

@@ -0,0 +1,99 @@
+function obj = axissymlap2d(type, n)
+%KERNEL.AXISSYMHELM2D   Construct the axissymmetric Helmholtz kernel.
+%   KERNEL.AXISSYMHELM2D('s', ZK) or KERNEL.AXISSYMHELM2D('single', ZK) 
+%   constructs the axissymmetric single-layer Helmholtz kernel with 
+%   wavenumber ZK.
+%
+%   KERNEL.AXISSYMHELM2D('d', ZK) or KERNEL.AXISSYMHELM2D('double', ZK) 
+%   constructs the axissymmetric double-layer Helmholtz kernel with 
+%   wavenumber ZK.
+%
+%   KERNEL.AXISSYMHELM2D('sp', ZK) or KERNEL.AXISSYMHELM2D('sprime', ZK) 
+%   constructs the normal derivative of the axissymmetric single-layer 
+%   Helmholtz kernel with wavenumber ZK.
+%
+%   KERNEL.AXISSYMHELM2D('c', ZK, COEFS) or 
+%   KERNEL.AXISSYMHELM2D('combined', ZK, COEFS)
+%   constructs the combined-layer axissymmetric Helmholtz kernel with 
+%   wavenumber ZK and parameter COEFS, 
+%   i.e., COEFS(1)*KERNEL.AXISSYMHELM2D('d', ZK) + 
+%         COEFS(2)*KERNEL.AXISSYMHELM2D('s', ZK).
+%
+% NOTES: The axissymetric kernels are currently supported only for purely
+% real or purely imaginary wave numbers
+%
+% See also CHNK.AXISSYMHELM2D.KERN.
+
+if ( nargin < 1 )
+    error('Missing Laplace kernel type.');
+end
+
+obj = kernel();
+obj.name = 'axissymlaplace';
+obj.opdims = [1 1];
+
+switch lower(type)
+
+    case {'s', 'single'}
+        obj.type = 's';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 's', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 's', n);
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    case {'d', 'double'}
+        obj.type = 'd';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'd', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'd', n);
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    case {'sp', 'sprime'}
+        obj.type = 'sp';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'sprime', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'sprime', n);
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    case {'dp', 'dprime'}
+        obj.type = 'dp';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'dprime', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'dprime', n);
+        obj.fmm = [];
+        obj.sing = 'hs';
+
+    case {'q', 'quad', 'quadruple', 'quadrupole'}
+        obj.type = 'q';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'q', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'q', n);
+        obj.fmm = [];
+        obj.sing = 'hs';
+
+    case {'spp', 'sprimeprime'}
+        obj.type = 'spp';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'spp', n);
+        obj.shifted_eval = @(s,t) chnk.axissymlap2d.kern(s, t, o, 'spp', n);
+        obj.fmm = [];
+        obj.sing = 'hs';
+
+    case {'q_sum_dp'}
+        obj.type = 'q_sum_dp';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'q_sum_dp', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'q_sum_dp', n);
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    case {'spp_sum_dp'}
+        obj.type = 'spp_sum_dp';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'spp_sum_dp', n);
+        obj.shifted_eval = @(s,t) chnk.axissymlap2d.kern(s, t, o, 'spp_sum_dp', n);
+        obj.fmm = [];
+        obj.sing = 'log';
+
+
+    otherwise
+        error('Unknown axissym Laplace kernel type ''%s''.', type);
+
+end
+
+end

chunkie/+chnk/+axissymlap2d/gaus_agm.m

@@ -0,0 +1,38 @@
+function [rk0,re0] = gaus_agm(x)
+
+eps  = 1E-15;
+a    = sqrt(2./(x+1));
+delt = 1./sqrt(1-a.*a);
+aa0  = delt + sqrt(delt.*delt-1);
+bb0   = 1./(delt+sqrt(delt.*delt-1));
+a0  = ones(size(delt));
+b0  = 1./delt;
+
+
+fact = ((a0+b0)/2).^2;
+
+for i=1:1000
+    a1 = (a0+b0)/2;
+    b1 = sqrt(a0.*b0);
+
+    aa1 = (aa0+bb0)/2;
+    bb1 = sqrt(aa0.*bb0);
+    a0 = a1;
+    b0 = b1;
+    aa0 = aa1;
+    bb0 = bb1;
+
+    c0 = (a1-b1)/2;
+    fact = fact-(c0.*c0)*2^(i);
+    drel = abs(a0-b0)./abs(a0);
+    drel2= abs(aa0-bb0)./abs(aa0);
+    if (max(drel+drel2) <2*eps)
+        break
+    end
+
+end
+
+rk0 = pi./(2*aa0.*sqrt(1-a.*a));
+re0 = pi*fact./(2*a0);
+
+end