axisymmetric laplace kernels

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

chunkie/+chnk/+axissymlap2d/gfunc.m

@@ -0,0 +1,18 @@
+function [gval,gdz,gdr,gdrp] = gfunc(r,rp,dr,z,zp,dz)
+    domega3 = 2*pi;
+    n       = 3;
+
+    t = (dz.^2+dr.^2)./(2.*r.*rp);
+    [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;
+
+    rfac = -r*(n-2)/2.*q0+(-(1+t).*r+rp).*q0d;
+    gdrp  = -domega3*(rp./r).^((n-2)/2)./(rp.*r) ...
+          .*rfac;
+    rfac = -rp*(n-2)/2.*q0+(-(1+t).*rp+r).*q0d;
+    gdr  = -domega3*(rp./r).^((n-2)/2)./(rp.*r) ...
+          .*rfac;
+end

chunkie/+chnk/+axissymlap2d/gfunc_on_axis.m

@@ -0,0 +1,18 @@
+function [gval,gdz,gdr,gdrp] = gfunc_on_axis(r,rp,dr,z,zp,dz)
+    domega3 = 2*pi;
+    n       = 3;
+
+    t = (dz.^2+dr.^2)./(2.*r.*rp);
+    [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;
+
+    rfac = -r*(n-2)/2.*q0+(-(1+t).*r+rp).*q0d;
+    gdrp  = -domega3*(rp./r).^((n-2)/2)./(rp.*r) ...
+          .*rfac;
+    rfac = -rp*(n-2)/2.*q0+(-(1+t).*rp+r).*q0d;
+    gdr  = -domega3*(rp./r).^((n-2)/2)./(rp.*r) ...
+          .*rfac;
+end

chunkie/+chnk/+axissymlap2d/green.m

@@ -0,0 +1,36 @@
+function [val, grad] = green(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
+%
+% Since the kernels are not translationally invariant in r, the size
+% of the gradient is 3, for d/dr, d/dr', and d/dz
+%
+
+[~, ns] = size(src);
+[~, nt] = size(targ);
+
+gtmp = zeros(nt, ns, 3);
+
+rt = repmat(targ(1,:).',1,ns);
+rs = repmat(src(1,:),nt,1);
+dz = repmat(src(2,:),nt,1)-repmat(targ(2,:).',1,ns);
+r  = (rt + origin(1));
+rp = (rs + origin(1));
+dr = (rs-rt);
+z  = zeros(size(rt));
+zp = zeros(size(rt));
+[g,gdz,gdr,gdrp] = chnk.axissymlap2d.gfunc(r,rp,dr,z,zp,dz);
+
+val = g;
+gtmp(:,:,1) = gdr;
+gtmp(:,:,2) = gdrp;
+gtmp(:,:,3) = gdz;
+grad = gtmp;
+
+

chunkie/+chnk/+axissymlap2d/kern.m

@@ -0,0 +1,81 @@
+function submat = kern(srcinfo, targinfo, origin, type)
+%CHNK.AXISSYMLAP2D.KERN axissymmetric Laplace layer potential kernels in 2D
+% 
+% Syntax: submat = chnk.axissymlap2d.kern(srcinfo,targingo,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(src, targ, origin);
+    nx = repmat(srcnorm(1,:), nt, 1);
+    ny = repmat(srcnorm(2,:), nt, 1);
+    % Due to lack of translation invariance in r, no sign flip needed, 
+    % as gradient is computed with repsect to r'
+    submat = (grad(:,:,2).*nx + grad(:,:,3).*ny);
+end
+
+if strcmpi(type, 'sprime')
+    targnorm = targinfo.n;
+    [~, grad] = chnk.axissymlap2d.green(src, targ, origin);
+
+    nx = repmat((targnorm(1,:)).',1,ns);
+    ny = repmat((targnorm(2,:)).',1,ns);
+    submat = (grad(:,:,1).*nx - grad(:,:,3).*ny);
+
+end
+
+if strcmpi(type, 's')
+    submat = chnk.axissymlap2d.green(src, targ, origin);
+
+end
+
+

chunkie/+chnk/+axissymlap2d/qeval0_far.m

@@ -0,0 +1,24 @@
+function [q0] = qeval0_far(x)
+
+        coefs = ...
+       [2.2214414690791831235E0,
+        0,                    
+        0.41652027545234683566E0,                    
+        0,                    
+        0.22778452563800217575E0,                    
+        0,                    
+        0.15660186137612649583E0,                    
+        0,                    
+        0.11928657409509635424E0]; 
+
+
+        t = 1./x;
+        t0 = 1./sqrt(x);
+        q0 = zeros(size(t));
+
+        for i=1:9
+                q0 = q0 + t0*coefs(i);
+                t0 = t0.*t;        
+        end
+
+end

chunkie/+chnk/+axissymlap2d/qeval0_near.m

@@ -0,0 +1,34 @@
+function [q0] = qeval0_near(t)
+
+        coefs = ...
+             [1.7328679513998632735E0,
+            -0.5000000000000000000E0,
+           -0.09160849392498290919E0,
+           0.062500000000000000000E0,
+           0.019905513916401443210E0,
+          -0.017578125000000000000E0,
+          -0.006097834693194945559E0,
+           0.006103515625000000000E0,
+          0.0021674343381175963469E0,
+         -0.0023365020751953125000E0,
+         -0.0008357538695841108955E0,
+          0.0009462833404541015625E0,
+          0.0003390851133264318725E0,
+        -0.00039757043123245239258E0,
+        -0.00014242013770157621830E0,
+         0.00017140153795480728149E0,
+         0.00006133159221782055041E0,
+        -0.00007532294148404616863E0];
+
+
+        b = log(t);
+        v = 1;
+
+        q0 = zeros(size(t));
+
+        for i=1:9
+                q0 = q0 + v*coefs(2*i-1)+v.*b*coefs(2*i);
+                v  = v.*t;        
+        end
+
+end

chunkie/+chnk/+axissymlap2d/qeval0der_near.m

@@ -0,0 +1,36 @@
+function [q0d] = qeval0der_near(t)
+
+        coefs = ...
+             [1.7328679513998632735d0,
+            -0.5000000000000000000d0,
+           -0.09160849392498290919d0,
+           0.062500000000000000000d0,
+           0.019905513916401443210d0,
+          -0.017578125000000000000d0,
+          -0.006097834693194945559d0,
+           0.006103515625000000000d0,
+          0.0021674343381175963469d0,
+         -0.0023365020751953125000d0,
+         -0.0008357538695841108955d0,
+          0.0009462833404541015625d0,
+          0.0003390851133264318725d0,
+        -0.00039757043123245239258d0,
+        -0.00014242013770157621830d0,
+         0.00017140153795480728149d0,
+         0.00006133159221782055041d0,
+        -0.00007532294148404616863d0];
+
+        b = log(t);
+        v = 1;
+
+        q0d = coefs(2)./t;
+
+        for i=2:9
+                q0d = q0d + (i-1)* ...
+                   (v*coefs(2*i-1)+v.*b*coefs(2*i)) ...
+                       +v*coefs(2*i);
+                v  = v.*t;        
+        end
+
+
+end

chunkie/+chnk/+axissymlap2d/qeval1_far.m

@@ -0,0 +1,24 @@
+function [q1] = qeval1_far(x)
+
+        coefs = ...
+       [0.55536036726979578088E0,                    
+        0,                    
+        0.26032517215771677229E0,                    
+        0,                    
+        0.17083839422850163181E0,                    
+        0,                    
+        0.12723901236810277786E0,                    
+        0,                    
+        0.10139358798083190111E0];  
+
+
+        t = 1./x;
+        t0 = 1./sqrt(x).^3;
+        q1 = zeros(size(t));
+
+        for i=1:9
+                q1 = q1 + t0*coefs(i);
+                t0 = t0.*t;        
+        end
+
+end

chunkie/+chnk/+axissymlap2d/qeval1_near.m

@@ -0,0 +1,34 @@
+function [q1] = qeval1_near(t)
+
+        coefs = ...
+            [-0.2671320486001367265d0,
+            -0.5000000000000000000d0,
+             0.5248254817749487276d0,
+           -0.18750000000000000000d0,
+           -0.04098835652733573868d0,
+           0.029296875000000000000d0, 
+           0.009513531070472923783d0,
+          -0.008544921875000000000d0,
+         -0.0029774361551467310174d0, 
+          0.0030040740966796875000d0, 
+          0.0010682069932178195667d0,
+         -0.0011565685272216796875d0,
+         -0.0004138797762860294822d0, 
+         0.00046985596418380737305d0, 
+          0.00016838776925222835322d0,
+        -0.00019777100533246994019d0,
+        -0.00007084821236213522235d0, 
+         0.00008536600034858565778d0];
+
+
+        b = log(t);
+        v = 1;
+
+        q1 = zeros(size(t));
+
+        for i=1:9
+                q1 = q1 + v*coefs(2*i-1)+v.*b*coefs(2*i);
+                v  = v.*t;        
+        end
+
+end

chunkie/+chnk/+axissymlap2d/qget_zero.m

@@ -0,0 +1,16 @@
+function [qm,qmp] = qget_zero(x)
+
+        a = sqrt(2./(x+1));
+        [fF,fE] = chnk.axissymlap2d.gaus_agm(x);
+
+        q0 = a.*fF;
+        q1 = x.*q0-2*fE./(a);
+
+        qa = q0;
+        qb = q1;
+
+        qmm = 0;
+        qm  = q0;
+        qmp = q1;
+
+end