Laplace axisymmetric kernels with all Fourier modes

chunkie/+chnk/+axissymlap2d/g0funcall.m

@@ -62,19 +62,7 @@
             iffwd = 0;
         end
 
-    end
-
-    %disp(['inside g0mall, x = ' num2str(x)])
-    %if (iffwd == 0)
-    %    disp(['running backward recursion']);
-    %end
-
-    %if (iffwd == 1)
-    %    disp(['running backward recursion']);
-    %end
-
-    %return
-
+    end    
 
 
     gvals = zeros(maxm+1,1);
@@ -177,7 +165,6 @@
     dernext = 0;
     der = 1;
 
-
     % run the downward recurrence
     for j = 1:(nterms-maxm+1)
         i = nterms-j+1;

chunkie/+chnk/+axissymlap2d/g0funcall_vec.m

@@ -0,0 +1,57 @@
+function [gval, gdz, gdr, gdrp, gdrpr, gdzz, gdrz, gdrpz] = g0funcall_vec(r, rp, dr, z, zp, dz, m)
+%
+% chnk.axissymlap2d.g0funcall evaluates a collection of axisymmetric Laplace
+% Green's functions, defined by the expression:
+%
+%     gfunc(n) = pi*rp * \int_0^{2\pi} 1/|x - x'| e^(-i n t) dt
+%
+% it is assumed that x = (x,0,z) otherwise the integral above should pick up a
+% phase factor out front of exp(i*n*phi), where phi is the azimuthal coordinate
+% of x in cylindrical coordinates.
+%
+% The extra factor of rp (and maybe pi?) out front makes subsequent interfacing
+% with RCIP slightly easier. Modes 0 through maxm are returned, with gval(1) =
+% mode 0 and gval(maxm+1) = mode maxm. The function is even, so g_{-n} = g_n.
+%
+% The above scaling should be consistent with what is in
+% chnk.axissymlap2d.gfunc, which is for merely the zero-mode
+% 
+
+    r = reshape(r,[1,size(r,1),size(r,2)]);
+    rp = reshape(rp,[1,size(rp,1),size(rp,2)]);
+    dr = reshape(dr,[1,size(dr,1),size(dr,2)]);
+    z = reshape(z,[1,size(z,1),size(z,2)]);
+    zp = reshape(zp,[1,size(zp,1),size(zp,2)]);
+    dz = reshape(dz,[1,size(dz,1),size(dz,2)]);
+
+    t = (dz.^2+dr.^2)./(2.*r.*rp);
+    chi = t+1;
+
+    [qm, qmd, qmdd] = chnk.axissymlap2d.qleg_half_miller_vec(t,m);
+
+    gval = 2*pi*sqrt(rp./r).*qm;
+    gdz  = 2*pi*sqrt(rp./r).*qmd ...
+          ./(rp.*r).*dz;
+
+    rfac = -r/2.*qm+(-(1+t).*r+rp).*qmd;
+    gdrp  = 2*pi*sqrt(rp./r)./(rp.*r) ...
+            .*rfac;
+
+    rfac = -rp/2.*qm+(-(1+t).*rp+r).*qmd;
+    gdr  = 2*pi*sqrt(rp./r)./(rp.*r) ...
+           .*rfac;
+
+    rfac = 1./(rp.*r).*qmd + (dz./(rp.*r)).^2.*qmdd;
+    gdzz = 2*pi*sqrt(rp./r).*rfac;
+
+    rfac = -3./(2*r.^2.*rp).*qmd + (-chi./(r.^2.*rp) + 1./(r.*rp.^2)).*qmdd;
+    gdrz = 2*pi*sqrt(rp./r).*dz.*rfac;
+
+    rfac = -3./(2*rp.^2.*r).*qmd + (-chi./(rp.^2.*r) + 1./(rp.*r.^2)).*qmdd;
+    gdrpz = 2*pi*sqrt(rp./r).*dz.*rfac;
+
+    rfac = 1./(4*r.*rp).*qm ...
+        + (2*chi./(rp.*r) - 3./(2*r.^2) - 3./(2*rp.^2)).*qmd ...
+        + (-chi./r+1./rp).*(-chi./rp + 1./r).*qmdd;
+    gdrpr = 2*pi*sqrt(rp./r).*rfac;
+end

chunkie/+chnk/+axissymlap2d/gfunc.m

@@ -1,4 +1,4 @@
-function [gval, gdz, gdr, gdrp] = gfunc(r, rp, dr, z, zp, dz)
+function [gval, gdz, gdr, gdrp, gdrpr, gdzz, gdrz, gdrpz] = gfunc(r, rp, dr, z, zp, dz)
 
 %
 % chnk.axissymlap2d.gfunc evaluates the zeroth order axisymmetric Laplace
@@ -14,8 +14,8 @@
     n       = 3;
 
     t = (dz.^2+dr.^2)./(2.*r.*rp);
-    [q0,q1,q0d] = chnk.axissymlap2d.qleg_half(t);
-    
+    [q0,q1,q0d] = chnk.axissymlap2d.qleg_half(t);    
+
     gval = domega3*(rp./r).^((n-2)/2).*q0;
     gdz  = domega3*(rp./r).^((n-2)/2).*q0d ...
           ./(rp.*r).*dz;
@@ -27,6 +27,27 @@
     rfac = -rp*(n-2)/2.*q0+(-(1+t).*rp+r).*q0d;
     gdr  = domega3*(rp./r).^((n-2)/2)./(rp.*r) ...
            .*rfac;
+
+
+    % Amandin hessian stuff
+    chi = 1+t;
+    q0dd = ((1/4)*q0 + 2*chi.*q0d)./(-t.^2-2*t);
+
+    rfac = 1./(rp.*r).*q0d + (dz./(rp.*r)).^2.*q0dd;
+    gdzz = 2*pi*sqrt(rp./r).*rfac;
+
+    rfac = -3./(2*r.^2.*rp).*q0d + (-chi./(r.^2.*rp) + 1./(r.*rp.^2)).*q0dd;
+    gdrz = 2*pi*sqrt(rp./r).*dz.*rfac;
+
+    rfac = -3./(2*rp.^2.*r).*q0d + (-chi./(rp.^2.*r) + 1./(rp.*r.^2)).*q0dd;
+    gdrpz = 2*pi*sqrt(rp./r).*dz.*rfac;
+
+    rfac = 1./(4*r.*rp).*q0 ...
+        + (2*chi./(rp.*r) - 3./(2*r.^2) - 3./(2*rp.^2)).*q0d ...
+        + (-chi./r+1./rp).*(-chi./rp + 1./r).*q0dd;
+    gdrpr = 2*pi*sqrt(rp./r).*rfac;
+
+
 
 
     % % check the scaling here correctly

chunkie/+chnk/+axissymlap2d/green_0th_mode.m

@@ -0,0 +1,53 @@
+function [val, grad, hess] = green_0th_mode(src, targ, origin)
+%
+% CHNK.AXISSYMHELM2D.GREEN evaluate the Laplace green's function
+% for the given sources and targets. 
+%
+% Note: that the first coordinate is r, and the second z.
+% The code relies on precomputed tables and hence loops are required for 
+% computing various pairwise interactions.
+% Finally, the code is not efficient in the sense that val, grad, hess 
+% are always internally computed independent of nargout
+%
+% Returns for gradient are:
+% grad = d_{r}, d_{r'}, d_{z}, d_{z'}
+%
+% Returns for hess are:
+% hess = d_{rr'}, d_{zz'}, d_{rz'}, d_{r'z}
+
+[~, ns] = size(src);
+[~, nt] = size(targ);
+
+rt = repmat(targ(1,:).',1,ns); % r
+rs = repmat(src(1,:),nt,1); % r'
+r  = (rt + origin(1));
+rp = (rs + origin(1));
+dr = rt-rs; % r - r'
+z  = repmat(targ(2,:).',1,ns);
+zp = repmat(src(2,:),nt,1);
+dz = z-zp; % z - z'
+
+[gs,gdzs,gdrs,gdrps,gdrpr,gdzz,gdrz,gdrpz] = chnk.axissymlap2d.gfunc(r,rp,dr,z,zp,dz);
+
+const = 1/(4*pi^2);
+
+val = gs*const;
+grad = zeros(nt, ns, 4);
+grad(:,:,1) = gdrs;
+grad(:,:,2) = gdrps;
+grad(:,:,3) = gdzs;
+grad(:,:,4) = -gdzs;
+grad = grad*const;
+
+hess = zeros(nt, ns, 4);
+hess(:,:,1) = gdrpr;
+hess(:,:,2) = -gdzz;
+hess(:,:,3) = -gdrz;
+hess(:,:,4) = gdrpz;
+hess = hess*const;
+
+% if anynan(hess)
+%   error('Execution stopped: NaN value detected.');
+% end
+
+end

chunkie/+chnk/+axissymlap2d/green_modal.m

@@ -0,0 +1,58 @@
+function [val, grad, hess] = green_modal(src, targ, origin, m)
+%
+% CHNK.AXISSYMHELM2D.GREEN evaluate the Laplace modal green's function
+% for the given sources and targets. 
+%
+% Note: that the first coordinate is r, and the second z.
+% The code is not efficient in the sense that val, grad, hess 
+% are always internally computed independent of nargout
+%
+% Returns for gradient are:
+% grad = d_{r}, d_{r'}, d_{z}, d_{z'}
+%
+% Returns for hess are:
+% hess = d_{rr'}, d_{zz'}, d_{rz'}, d_{r'z}
+%
+% m is the mode we return
+%
+% all_modes = true: means we return size [m*nt, ns] instead of [nt, ns]
+
+[~, ns] = size(src);
+[~, nt] = size(targ);
+
+rt = repmat(targ(1,:).',1,ns); % r
+rs = repmat(src(1,:),nt,1); % r'
+r  = (rt + origin(1));
+rp = (rs + origin(1));
+dr = rt-rs; % r - r'
+z  = repmat(targ(2,:).',1,ns);
+zp = repmat(src(2,:),nt,1);
+dz = z-zp; % z - z'
+
+[gs,gdzs,gdrs,gdrps,gdrprs,gdzzs,gdrzs,gdrpzs] = chnk.axissymlap2d.g0funcall_vec(r,rp,dr,z,zp,dz,m);
+
+
+const = 1/(4*pi^2);
+
+val = gs(m,:,:);
+val = reshape(val,[nt,ns]);
+val = const*val;
+
+grad = zeros(nt, ns, 4);
+grad(:,:,1) = gdrs(m,:,:);
+grad(:,:,2) = gdrps(m,:,:);
+grad(:,:,3) = gdzs(m,:,:);
+grad(:,:,4) = -gdzs(m,:,:);
+grad = const*grad;
+
+hess = zeros(nt, ns, 4);
+hess(:,:,1) = gdrprs(m,:,:);
+hess(:,:,2) = -gdzzs(m,:,:);
+hess(:,:,3) = -gdrzs(m,:,:);
+hess(:,:,4) = gdrpzs(m,:,:);
+hess = const*hess;
+
+
+
+
+end

chunkie/+chnk/+axissymlap2d/kern.m

@@ -74,8 +74,28 @@
 end
 
 if strcmpi(type, 's')
+<<<<<<< HEAD
+    submat = chnk.axissymlap2d.green_0th_mode(src, targ, origin);
+end
+
+if strcmpi(type, 'dprime')
+    targnorm = targinfo.n;
+    srcnorm = srcinfo.n;
+    nxt = repmat((targnorm(1,:)).',1,ns);
+    nyt = repmat((targnorm(2,:)).',1,ns);
+    nxs = repmat(srcnorm(1,:), nt, 1);
+    nys = repmat(srcnorm(2,:), nt, 1);
+
+    % hess = d_{rr'}, d_{zz'}, d_{rz'}, d_{r'z}
+    [~, ~, hess] = chnk.axissymlap2d.green_0th_mode(src, targ, origin);
+    submat = (hess(:,:,1).*nxt.*nxs + hess(:,:,3).*nys.*nxt ...
+        + hess(:,:,4).*nxs.*nyt + hess(:,:,2).*nyt.*nys);
+end
+
+=======
     submat = chnk.axissymlap2d.green(src, targ, origin);
 
 end
 
 
+>>>>>>> upstream/lap_axisym

chunkie/+chnk/+axissymlap2d/kern_0th_mode.m

@@ -0,0 +1,119 @@
+function submat = kern_0th_mode(srcinfo, targinfo, origin, type)
+%CHNK.AXISSYMLAP2D.KERN axissymmetric Laplace layer potential kernels in 2D
+% 
+% Syntax: submat = chnk.axissymlap2d.kern(srcinfo,targinfo,type)
+%
+% 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). Note that the normal information is obtained
+% by taking the perpendicular to the provided tangential deriviative
+% info and normalizing  
+%
+% Here the first and second components correspond to the r and z
+% coordinates respectively. 
+%
+% Kernels based on G(x,y) = \int_{0}^{\pi} 1/(d(t)) \, dt \, 
+% where d(t) = \sqrt(r^2 + r'^2 - 2rr' \cos(t) + (z-z')^2) with
+% x = (r,z), and y = (r',z')
+%
+% D(x,y) = \nabla_{n_y} G(x,y)
+% S(x,y) = G(x,y)
+% S'(x,y) = \nabla_{n_x} G(x,y)
+% D'(x,y) = \nabla_{n_x} \nabla_{n_y} G(x,y)
+%
+% Input:
+%   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,:)
+%   targinfo - description of targets in ptinfo struct format,
+%                if info not relevant (d/d2) it doesn't need to
+%                be provided. sprime requires tangent info in
+%                targinfo.d
+%   type - string, determines kernel type
+%                type == 'd', double layer kernel D
+%                type == 's', single layer kernel S
+%                type == 'sprime', normal derivative of single
+%                      layer S'
+%
+%
+% Output:
+%   submat - the evaluation of the selected kernel for the
+%            provided sources and targets. the number of
+%            rows equals the number of targets and the
+%            number of columns equals the number of sources  
+%
+%
+
+src = srcinfo.r; 
+targ = targinfo.r;
+
+[~, ns] = size(src);
+[~, nt] = size(targ);
+
+if strcmpi(type, 'd')
+    srcnorm = srcinfo.n;
+    [~, grad] = chnk.axissymlap2d.green_0th_mode(src, targ, origin);
+    nx = repmat(srcnorm(1,:), nt, 1);
+    ny = repmat(srcnorm(2,:), nt, 1);
+    % dr'*nr' + dz'*nz'
+    submat = -(grad(:,:,2).*nx + grad(:,:,4).*ny);
+end
+
+if strcmpi(type, 's')
+    [val, ~] = chnk.axissymlap2d.green_0th_mode(src, targ, origin);
+    submat = val;
+end
+
+if strcmpi(type, 'sprime')
+    targnorm = targinfo.n;
+    [~, grad] = chnk.axissymlap2d.green_0th_mode(src, targ, origin);
+    nx = repmat((targnorm(1,:)).',1,ns);
+    ny = repmat((targnorm(2,:)).',1,ns);
+    % dr*nr + dz*nz
+    submat = -(grad(:,:,1).*nx + grad(:,:,3).*ny);
+end
+
+if strcmpi(type, 'sgradr')
+    [~, grad] = chnk.axissymlap2d.green_0th_mode(src, targ, origin);
+    submat = grad(:,:,1);
+end
+
+if strcmpi(type, 'sgradz')
+    [~, grad] = chnk.axissymlap2d.green_0th_mode(src, targ, origin);
+    submat = grad(:,:,3);
+end
+
+if strcmpi(type, 'dgradr')
+    h = 1e-200;
+    targ_cr = targinfo;
+    targ_cr.r(1,:) = targ_cr.r(1,:) + 1i*h;
+    D_cr = kern_0th_mode(srcinfo, targ_cr, origin, 'd');
+    submat = imag(D_cr)/h;
+end
+
+if strcmpi(type, 'dgradz')
+    h = 1e-200;
+    targ_cz = targinfo;
+    targ_cz.r(2,:) = targ_cz.r(2,:) + 1i*h;
+    D_cz = kern_0th_mode(srcinfo, targ_cz, origin, 'd');
+    submat = imag(D_cz)/h;
+end
+
+if strcmpi(type, 'dprime')
+    targnorm = targinfo.n;
+    srcnorm = srcinfo.n;
+    nxt = repmat((targnorm(1,:)).',1,ns);
+    nyt = repmat((targnorm(2,:)).',1,ns);
+    nxs = repmat(srcnorm(1,:), nt, 1);
+    nys = repmat(srcnorm(2,:), nt, 1);
+
+    % hess = d_{rr'}, d_{zz'}, d_{rz'}, d_{r'z}
+    [~, ~, hess] = chnk.axissymlap2d.green_0th_mode(src, targ, origin);
+    submat = (hess(:,:,1).*nxt.*nxs + hess(:,:,3).*nys.*nxt ...
+        + hess(:,:,4).*nxs.*nyt + hess(:,:,2).*nyt.*nys);
+end
+
+end