Merge pull request #105 from haiszhu/dev-helmkernsplit-redo

laplace and helmholtz kernel split quadrature

Author: mrachh

chunkie/+chnk/pquadwts.m

@@ -1,56 +1,78 @@
-function mat = pquadwts(r,d,n,d2,ct,bw,j,...
-    rt,dt,nt,d2t,kern,opdims,t,w,opts,intp_ab,intp)
-%CHNK.INTERPQUADWTS product integration for interaction of kernel on chunk 
+function [varargout] = pquadwts(r,d,n,d2,wts,j,...
+    rt,t,w,opts,intp_ab,intp,types)
+%CHNK.pquadwts product integration for interaction of kernel on chunk 
 % at targets
 %
 % WARNING: this routine is not designed to be user-callable and assumes 
 %   a lot of precomputed values as input
 %
-% Syntax: mat = interpquadwts(r,d,d2,h,ct,bw,j, ...
-%   rt,dt,d2t,kern,opdims,t,w,opts)
+% Syntax: [varargout] = pquadwts(r,d,n,d2,wts,j,...
+%    rt,t,w,opts,intp_ab,intp,types)
 %
 % Input:
 %   r - chnkr nodes
 %   d - chnkr derivatives at nodes
 %   n - chnkr normals at nodes
 %   d2 - chnkr 2nd derivatives at nodes
-%   ct - Legendre nodes at order of chunker
-%   bw - barycentric interpolation weights for Legendre nodes at order of
-%   chunker
+%   wts - chnkr integration weights for smooth functions
 %   j - chunk of interest
-%   rt,dt,nt,d2t - position, derivative, normal,second derivative of select 
-%               target points. if any are not used by kernel (or not well
-%               defined, e.g. when not on curve), a dummy array
-%               of the appropriate size should be supplied
-%   kern - kernel function of form kern(srcinfo,targinfo)
-%   opdims - dimensions of kernel
-%   t - (Legendre) integration nodes for adaptive integration
-%   w - integration nodes for adaptive integrator (t and w not necessarily 
-%       same order as chunker order)
+%   rt - position of target points. if any are not used by kernel 
+%               (or not well defined, e.g. when not on curve), a 
+%               dummy array of the appropriate size should be supplied
+%   t - (Legendre) integration nodes
+%   w - (Legendre) integration weights
+%   opts - opts.side
+%   intp_ab - panel endpoints interpolation matrix
+%   types - specified singularity types
+%       [0 0  0 0] - smooth quadr
+%       [1 0  0 0] - log singularity
+%       [0 0 -1 0] - cauchy singularity
+%       each [a, b, c, d] 
+%       corresponds to (log(z-w))^a (log(zc-wc))^b (z-w)^c (zc-wc)^d
 %
 % Output
-%   mat - integration matrix
+%   varargout - integration matrices for specified singularity types
 %
 
 % Helsing-Ojala (interior/exterior?)
-
 xlohi = intp_ab*(r(1,:,j)'+1i*r(2,:,j)');         % panel end points
 r_i = intp*(r(1,:,j)'+1i*r(2,:,j)');              % upsampled panel
-d_i = (intp*(d(1,:,j)'+1i*d(2,:,j)'));        % r'
-d2_i = (intp*(d2(1,:,j)'+1i*d2(2,:,j)'));   % r''
+d_i = (intp*(d(1,:,j)'+1i*d(2,:,j)'));            % r'
+d2_i = (intp*(d2(1,:,j)'+1i*d2(2,:,j)'));         % r''
+wts_i = wts(:,j)';                                %
 sp = abs(d_i); tang = d_i./sp;                    % speed, tangent
 n_i = -1i*tang;                                   % normal
 cur = -real(conj(d2_i).*n_i)./sp.^2;              % curvature
-wxp_i = w.*d_i;                                     % complex speed weights (Helsing's wzp)
+wxp_i = w.*d_i;                                   % complex speed weights (Helsing's wzp)
+
+varargout = cell(size(types));
+
+% [As, Ad, A1, A2, A3, A4] = SDspecialquad(struct('x',rt(1,:)' + 1i*rt(2,:)'),...
+%                     struct('x',r_i,'nx',n_i,'wxp',wxp_i),xlohi(1),xlohi(2),opts.side);
+[As, Ad] = SDspecialquad(struct('x',rt(1,:)' + 1i*rt(2,:)'),...
+                    struct('x',r_i,'nx',n_i,'wxp',wxp_i),xlohi(1),xlohi(2),opts.side);
+As = As*intp;
+Ad  = Ad*intp;
 
-mat_ho_slp = Sspecialquad(struct('x',rt(1,:)' + 1i*rt(2,:)'),...
-                          struct('x',r_i,'nx',n_i,'wxp',wxp_i),xlohi(1),xlohi(2),'e');
-mat_ho_dlp = Dspecialquad(struct('x',rt(1,:)' + 1i*rt(2,:)'),...
-                          struct('x',r_i,'nx',n_i,'wxp',wxp_i),xlohi(1),xlohi(2),'e');
+for j = 1:length(types)
+    type0 = types{j};
 
-mat = (mat_ho_slp+real(mat_ho_dlp))*intp;  % depends on kern, different mat1?
+    if (all(type0 == [0 0 0 0]))
+      varargout{j} = ones(size(rt,2),numel(wts_i)).*wts_i;
+    elseif (all(type0 == [1 0 0 0]))
+        % varargout{j} = Sspecialquad(struct('x',rt(1,:)' + 1i*rt(2,:)'),...
+        %       struct('x',r_i,'nx',n_i,'wxp',wxp_i),xlohi(1),xlohi(2),opts.side)*intp;
+        varargout{j} = As;
+    elseif (all(type0 == [0 0 -1 0]))
+            % varargout{j} = Dspecialquad(struct('x',rt(1,:)' + 1i*rt(2,:)'),...
+            %               struct('x',r_i,'nx',n_i,'wxp',wxp_i),xlohi(1),xlohi(2),opts.side)*intp;
+            varargout{j} = Ad;
+    else
+        error("split panel quad type " + join(string([1 2 3]),",") + " not available");
+    end
 
 end
+end
 
 function [A, A1, A2] = Sspecialquad(t,s,a,b,side)
 % SSPECIALQUAD - SLP val+grad close-eval Helsing "special quadrature" matrix
@@ -66,6 +88,7 @@
 % Efficient only if multiple targs, since O(p^3).
 % See Helsing-Ojala 2008 (special quadr Sec 5.1-2), Helsing 2009 mixed (p=16),
 % and Helsing's tutorial demo11b.m LGIcompRecFAS()
+% See https://github.com/ahbarnett/BIE2D/blob/master/panels/LapSLP_closepanel.m
 if nargin<5, side = 'i'; end     % interior or exterior
 zsc = (b-a)/2; zmid = (b+a)/2; % rescaling factor and midpoint of src segment
 y = (s.x-zmid)/zsc; x = (t.x-zmid)/zsc;  % transformed src nodes, targ pts
@@ -146,6 +169,7 @@
 % Efficient only if multiple targs, since O(p^3).
 % See Helsing-Ojala 2008 (special quadr Sec 5.1-2), Helsing 2009 mixed (p=16),
 % and Helsing's tutorial demo11b.m M1IcompRecFS()
+% See https://github.com/ahbarnett/BIE2D/blob/master/panels/LapDLP_closepanel.m
 if nargin<5, side = 'i'; end     % interior or exterior
 zsc = (b-a)/2; zmid = (b+a)/2; % rescaling factor and midpoint of src segment
 y = (s.x-zmid)/zsc; x = (t.x-zmid)/zsc;  % transformed src nodes, targ pts
@@ -203,3 +227,76 @@
     end
 end
 end
+
+function [As, A, A1, A2, A3, A4] = SDspecialquad(t,s,a,b,side)
+% S+D together... (with Sn) should be the same as requesting S or D
+% more expensive if request Dn...
+% https://github.com/ahbarnett/BIE2D/blob/master/panels/LapDLP_closepanel.m
+%%% d = 100;
+zsc = (b-a)/2; zmid = (b+a)/2; % rescaling factor and midpoint of src segment
+y = (s.x-zmid)/zsc; x = (t.x-zmid)/zsc;  % transformed src nodes, targ pts
+%figure; plot(x,'.'); hold on; plot(y,'+-'); plot([-1 1],[0 0],'ro'); % debug
+N = numel(x);                            % # of targets
+p = numel(s.x);                          % assume panel order is # nodes
+if N*p==0
+    As = 0; A = 0; A1=0; A2=0;
+    return
+end
+c = (1-(-1).^(1:p))./(1:p);              % Helsing c_k, k = 1..p.
+V = ones(p,p); for k=2:p, V(:,k) = V(:,k-1).*y; end  % Vandermonde mat @ nodes
+P = zeros(p+1,N);      % Build P, Helsing's p_k vectorized on all targs...
+d = 1.1; inr = abs(x)<=d; ifr = abs(x)>d;      % near & far treat separately
+%gam = 1i;
+gam = exp(1i*pi/4);  % smaller makes cut closer to panel. barnett 4/17/18
+if side == 'e', gam = conj(gam); end   % note gam is a phase, rots branch cut
+P(1,:) = log(gam) + log((1-x)./(gam*(-1-x)));  % init p_1 for all targs int
+
+% upwards recurrence for near targets, faster + more acc than quadr...
+% (note rotation of cut in log to -Im; so cut in x space is lower unit circle)
+if N>1 || (N==1 && inr==1) % Criterion added by Bowei Wu 03/05/15 to ensure inr not empty
+  for k=1:p 
+    P(k+1,inr) = x(inr).'.*P(k,inr) + c(k); 
+  end  % recursion for p_k
+end
+% fine quadr (no recurrence) for far targets (too inaccurate cf downwards)...
+Nf = numel(find(ifr)); wxp = s.wxp/zsc; % rescaled complex speed weights
+if Nf>0 % Criterion added by Bowei Wu 03/05/15 to ensure ifr not empty
+  P(end,ifr) = sum(((wxp.*(V(:,end).*y(:)))*ones(1,Nf))./bsxfun(@minus,y,x(ifr).'));  % int y^p/(y-x)
+  for ii = p:-1:2
+    P( ii,ifr) = (P(ii+1,ifr)-c(ii))./x(ifr).';
+  end
+end
+
+Q = zeros(p,N); % compute q_k from p_k via Helsing 2009 eqn (18)... (p even!)
+% Note a rot ang appears here too...  4/17/18
+%gam = exp(1i*pi/4); % 1i;  % moves a branch arc as in p_1
+%if side == 'e', gam = conj(gam); end   % note gam is a phase, rots branch cut
+Q(1:2:end,:) = P(2:2:end,:) - repmat(log((1-x.').*(-1-x.')),[ceil(p/2) 1]); % guessed!
+% (-1)^k, k odd, note each log has branch cut in semicircle from -1 to 1
+% Q(2:2:end,:) = P(3:2:end,:) - repmat(log((1-x.')./((-1-x.'))),[floor(p/2) 1]);  % same cut as for p_1
+Q(2:2:end,:) = P(3:2:end,:) - repmat(log(gam) + log((1-x.')./(gam*(-1-x.'))),[floor(p/2) 1]);  % same cut as for p_1
+Q = Q.*repmat(1./(1:p)',[1 N]); % k even
+As = real((V.'\Q).'.*repmat((1i*s.nx)',[N 1])*zsc)/(2*pi*abs(zsc));
+As = As*abs(zsc) - log(abs(zsc))/(2*pi)*repmat(abs(s.wxp)',[N 1]); % unscale, yuk
+warning('off','MATLAB:nearlySingularMatrix');
+% A = real((V.'\P).'*(1i/(2*pi)));         % solve for special quadr weights
+A = ((V.'\P(1:p,:)).'*(1i/(2*pi)));         % do not take real for the eval of Stokes DLP non-laplace term. Bowei 10/19/14
+%A = (P.'*inv(V))*(1i/(2*pi));   % equiv in exact arith, but not bkw stable.
+if nargout>2
+  R =  -(kron(ones(p,1),1./(1-x.')) + kron((-1).^(0:p-1).',1./(1+x.'))) +...
+      repmat((0:p-1)',[1 N]).*[zeros(1,N); P(1:p-1,:)];  % hypersingular kernel weights of Helsing 2009 eqn (14)
+  Az = (V.'\R).'*(1i/(2*pi*zsc));  % solve for targ complex-deriv mat & rescale
+  A1 = Az;
+  if nargout > 3
+    S = -(kron(ones(p,1),1./(1-x.').^2) - kron((-1).^(0:p-1).',1./(1+x.').^2))/2 +...
+        repmat((0:p-1)',[1 N]).*[zeros(1,N); R(1:p-1,:)]/2; % supersingular kernel weights
+    Azz = (V.'\S).'*(1i/(2*pi*zsc^2));
+    if nargout > 4
+      A1 = real(Az); A2 = -imag(Az);  % note sign for y-deriv from C-deriv
+      A3 = real(Azz); A4 = -imag(Azz);    
+    else
+      A1 = Az; A2 = Azz; 
+    end
+  end
+end
+end

chunkie/@kernel/helm2d.m

@@ -40,12 +40,20 @@
         obj.eval = @(s,t) chnk.helm2d.kern(zk, s, t, 's');
         obj.fmm  = @(eps,s,t,sigma) chnk.helm2d.fmm(eps, zk, s, t, 's', sigma);
         obj.sing = 'log';
+        obj.splitinfo = [];
+        obj.splitinfo.type = {[0 0 0 0],[1 0 0 0]};
+        obj.splitinfo.action = {'r','r'};
+        obj.splitinfo.functions = @(s,t) helm2d_s_split(zk,s,t);
 
     case {'d', 'double'}
         obj.type = 'd';
         obj.eval = @(s,t) chnk.helm2d.kern(zk, s, t, 'd');
         obj.fmm  = @(eps,s,t,sigma) chnk.helm2d.fmm(eps, zk, s, t, 'd', sigma);
         obj.sing = 'log';
+        obj.splitinfo = [];
+        obj.splitinfo.type = {[0 0 0 0],[1 0 0 0],[0 0 -1 0]};
+        obj.splitinfo.action = {'r','r','r'};
+        obj.splitinfo.functions = @(s,t) helm2d_d_split(zk,s,t);
 
     case {'sp', 'sprime'}
         obj.type = 'sp';
@@ -69,6 +77,10 @@
         obj.eval = @(s,t) chnk.helm2d.kern(zk, s, t, 'c', coefs);
         obj.fmm  = @(eps,s,t,sigma) chnk.helm2d.fmm(eps, zk, s, t, 'c', sigma, coefs);
         obj.sing = 'log';
+        obj.splitinfo = [];
+        obj.splitinfo.type = {[0 0 0 0],[1 0 0 0],[0 0 -1 0]};
+        obj.splitinfo.action = {'r','r','r'};
+        obj.splitinfo.functions = @(s,t) helm2d_c_split(zk,s,t,coefs);
 
     case {'cp', 'cprime'}
         if ( nargin < 3 )
@@ -93,3 +105,37 @@
 
 
 end
+
+function f = helm2d_s_split(zk,s,t)
+seval = chnk.helm2d.kern(zk, s, t, 's');
+dist = (s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)');
+logeval = log(abs(dist));
+f = cell(2,1);
+f{1} = seval + 1/(2*pi)*4*logeval.*imag(seval);
+f{2} = 4*imag(seval);
+end
+
+function f = helm2d_d_split(zk,s,t)
+deval = chnk.helm2d.kern(zk, s, t, 'd');
+dist = (s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)');
+logeval = log(abs(dist));
+cauchyeval = (s.n(1,:)+1i*s.n(2,:))./(dist);
+f = cell(3, 1);
+f{1} = deval + 1/(2*pi)*4*logeval.*imag(deval) + 1/(2*pi)*real(cauchyeval);
+f{2} = 4*imag(deval);
+f{3} = ones(size(t.r,2),size(s.r,2));
+end
+
+function f = helm2d_c_split(zk,s,t,coefs)
+seval = chnk.helm2d.kern(zk, s, t, 's');
+deval = chnk.helm2d.kern(zk, s, t, 'd');
+dist = (s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)');
+logeval = log(abs(dist));
+cauchyeval = (s.n(1,:)+1i*s.n(2,:))./(dist);
+f = cell(3, 1);
+f{1} = coefs(1) * (deval + 2/pi*logeval.*imag(deval) + 1/(2*pi)*real(cauchyeval)) + ...
+       coefs(2) * (seval + 2/pi*logeval.*imag(seval));
+f{2} = coefs(1) * (4*imag(deval)) + coefs(2) * (4*imag(seval));
+f{3} = coefs(1)*ones(size(t.r,2),size(s.r,2));
+end
+

chunkie/@kernel/lap2d.m

@@ -42,12 +42,20 @@
         obj.eval = @(s,t) chnk.lap2d.kern(s, t, 's');
         obj.fmm  = @(eps,s,t,sigma) chnk.lap2d.fmm(eps, s, t, 's', sigma);
         obj.sing = 'log';
+        obj.splitinfo = [];
+        obj.splitinfo.type = {[1 0 0 0]};
+        obj.splitinfo.action = {'r'};
+        obj.splitinfo.functions = @(s,t) lap2d_s_split(s,t);
 
     case {'d', 'double'}
         obj.type = 'd';
         obj.eval = @(s,t) chnk.lap2d.kern(s, t, 'd');
         obj.fmm  = @(eps,s,t,sigma) chnk.lap2d.fmm(eps, s, t, 'd', sigma);
         obj.sing = 'smooth';
+        obj.splitinfo = [];
+        obj.splitinfo.type = {[0 0 -1 0]};
+        obj.splitinfo.action = {'r'};
+        obj.splitinfo.functions = @(s,t) lap2d_d_split(s,t);
 
     case {'sp', 'sprime'}
         obj.type = 'sp';
@@ -91,6 +99,10 @@
         obj.eval = @(s,t) chnk.lap2d.kern(s, t, 'c', coefs);
         obj.fmm  = @(eps,s,t,sigma) chnk.lap2d.fmm(eps, s, t, 'c', sigma, coefs);
         obj.sing = 'log';
+        obj.splitinfo = [];
+        obj.splitinfo.type = {[1 0 0 0],[0 0 -1 0]};
+        obj.splitinfo.action = {'r','r'};
+        obj.splitinfo.functions = @(s,t) lap2d_c_split(s,t,coefs);
 
     otherwise
         error('Unknown Laplace kernel type ''%s''.', type);
@@ -103,3 +115,20 @@
 end
 
 end
+
+function f = lap2d_s_split(s,t)
+f = cell(1,1);
+f{1} = ones(size(t.r,2),size(s.r,2));
+end
+
+function f = lap2d_d_split(s,t)
+f = cell(1,1);
+f{1} = ones(size(t.r,2),size(s.r,2));
+end
+
+function f = lap2d_c_split(s,t,coefs)
+f = cell(2,1);
+f{1} = coefs(2)*ones(size(t.r,2),size(s.r,2));
+f{2} = coefs(1)*ones(size(t.r,2),size(s.r,2));
+end
+

chunkie/@kernel/stok2d.m

@@ -78,6 +78,10 @@
         obj.fmm = @(eps, s, t, sigma) chnk.stok2d.fmm(eps, mu, s, t, 'svel', sigma);
         obj.opdims = [2, 2];
         obj.sing = 'log';
+        obj.splitinfo = [];
+        obj.splitinfo.type = {[1 0 0 0],[0 0 -1 0],[0 0 -1 0]};
+        obj.splitinfo.action = {'r','r','i'};
+        obj.splitinfo.functions = @(s,t) stok2d_s_split(mu, s, t);
 
     case {'spres', 'spressure'}
         obj.type = 'spres';
@@ -188,3 +192,27 @@
 
 
 end
+
+function f = stok2d_s_split(mu,s,t)
+dist = (s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)');
+distr = real(dist);
+disti = imag(dist);
+f = cell(3, 1);
+ntarg = numel(t.r(1,:));
+nsrc = numel(s.r(1,:));
+sn1mat = repmat(s.n(1,:),[ntarg 1]);
+sn2mat = repmat(s.n(2,:),[ntarg 1]);
+f{1} = zeros(2*ntarg,2*nsrc);
+f{2} = zeros(2*ntarg,2*nsrc);
+f{3} = zeros(2*ntarg,2*nsrc);
+f{1}(1:2:end,1:2:end) = 1/(2*mu);
+f{1}(2:2:end,2:2:end) = 1/(2*mu);
+f{2}(1:2:end,1:2:end) = -sn1mat.*distr/(2*mu);     
+f{2}(1:2:end,2:2:end) = -sn2mat.*distr/(2*mu);
+f{2}(2:2:end,1:2:end) = -sn1mat.*disti/(2*mu); 
+f{2}(2:2:end,2:2:end) = -sn2mat.*disti/(2*mu);
+f{3}(1:2:end,1:2:end) = -sn2mat.*distr/(2*mu);
+f{3}(1:2:end,2:2:end) =  sn1mat.*distr/(2*mu);
+f{3}(2:2:end,1:2:end) = -sn2mat.*disti/(2*mu);
+f{3}(2:2:end,2:2:end) =  sn1mat.*disti/(2*mu);
+end