updated kernels

chunkie/+chnk/+flex2d/besseldiff_etc_pscoefs.m

@@ -0,0 +1,32 @@
+function [cf1,cf2] = besseldiff_etc_pscoefs(nterms)
+% powerseries coefficients for J_0 - 1 and the other series in
+% definition of Y_0
+%
+% these are both series with the terms z^2,z^4,z^6,..z^(2*nterms)
+%
+% cf1 are power series coeffs of even terms (excluding constant) 
+% for J_0(z) - 1
+%
+% cf2 are power series coeffs for the similar series
+% z^2/4 - (1+1/2)*(z^2/4)^2/(2!)^2 + (1+1/2+1/3)*(z^2/4)^3/(3!)^2 - ...
+%
+% which appears in the power series for H_0(z)
+
+sum1 = 0;
+fac1 = 1;
+pow1 = 1;
+
+cf1 = zeros(nterms,1);
+cf2 = zeros(nterms,1);
+sgn = 1;
+
+for j = 1:nterms
+    fac1 = fac1/(j*j);
+    sum1 = sum1 + 1.0/j;
+    pow1 = pow1*0.25;
+    cf1(j) = -sgn*pow1*fac1;
+    cf2(j) = sgn*sum1*pow1*fac1;
+    sgn = -sgn;
+end
+
+end

chunkie/+chnk/+flex2d/even_pseval.m

@@ -0,0 +1,17 @@
+function [v] = even_pseval(cf,x)
+% evaluates an even power series with no constant term 
+% with coefficients given in cf. 
+%
+% x can be an array
+
+x2 = x.*x;
+
+xp = x2;
+
+v = zeros(size(x));
+
+nterms = length(cf);
+for j = 1:nterms
+    v = v + cf(j)*xp;
+    xp = xp.*x2;
+end

chunkie/+chnk/+flex2d/fmm.m

@@ -0,0 +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

chunkie/+chnk/+flex2d/green.m

@@ -0,0 +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