From 9cf3a62ebe203307a57cf4efefd50f7adf6df239 Mon Sep 17 00:00:00 2001
From: Hai Zhu <zhuh1106@gmail.com>
Date: Thu, 17 Oct 2024 15:11:33 -0700
Subject: [PATCH 1/8] helm kern split redo

---
 chunkie/+chnk/pquadwts.m     |  64 +++++++++++--------
 chunkie/@kernel/helm2d.m     |  30 +++++++++
 chunkie/@kernel/lap2d.m      |   9 +++
 chunkie/chunkerkerneval.m    |  68 +++++++++++++++------
 chunkie/chunkerkernevalmat.m |  37 +++++++++--
 devtools/test/pquadTest.m    | 115 +++++++++++++++++++++++++++++++++++
 6 files changed, 271 insertions(+), 52 deletions(-)
 create mode 100644 devtools/test/pquadTest.m

diff --git a/chunkie/+chnk/pquadwts.m b/chunkie/+chnk/pquadwts.m
index 59bc795..8edccfd 100644
--- a/chunkie/+chnk/pquadwts.m
+++ b/chunkie/+chnk/pquadwts.m
@@ -1,55 +1,63 @@
-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
 %
 % 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)
 
-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');
+varargout = cell(size(types));
 
-mat = (mat_ho_slp+real(mat_ho_dlp))*intp;  % depends on kern, different mat1?
+for j = 1:length(types)
+    type0 = types{j};
 
+    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;
+    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;
+    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)
@@ -66,6 +74,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 +155,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
diff --git a/chunkie/@kernel/helm2d.m b/chunkie/@kernel/helm2d.m
index 057b5d3..0ee5ec8 100644
--- a/chunkie/@kernel/helm2d.m
+++ b/chunkie/@kernel/helm2d.m
@@ -40,12 +40,28 @@
         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.function = cell(2,1);
+        obj.splitinfo.function{1} = @(s,t) obj.eval(s,t) ...
+                                         + 1/(2*pi)*4*log(abs((s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)'))).*imag(obj.eval(s,t));
+        obj.splitinfo.function{2} = @(s,t) 4*imag(obj.eval(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.function = cell(3,1);
+        obj.splitinfo.function{1} = @(s,t) obj.eval(s,t) ...
+                                           + 1/(2*pi)*4*log(abs((s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)'))).*imag(obj.eval(s,t)) ... 
+                                           + 1/(2*pi)*real((s.n(1,:)+1i*s.n(2,:))./((s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)')));
+        obj.splitinfo.function{2} = @(s,t) 4*imag(obj.eval(s,t));
+        obj.splitinfo.function{3} = @(s,t) ones(size(t.r,2),size(s.r,2));
 
     case {'sp', 'sprime'}
         obj.type = 'sp';
@@ -69,6 +85,20 @@
         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'};
+        deval = @(s,t) chnk.helm2d.kern(zk, s, t, 'd', coefs(1));
+        seval = @(s,t) chnk.helm2d.kern(zk, s, t, 's', coefs(2));
+        obj.splitinfo.function = cell(3,1);
+        obj.splitinfo.function{1} = @(s,t) coefs(1) * (deval(s,t) ...
+                                                       + 1/(2*pi)*4*log(abs((s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)'))).*imag(deval(s,t)) ... 
+                                                       + 1/(2*pi)*real((s.n(1,:)+1i*s.n(2,:))./((s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)')))) ...
+                                         + coefs(2) * (seval(s,t) ...
+                                                       + 1/(2*pi)*4*log(abs((s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)'))).*imag(seval(s,t)));
+        obj.splitinfo.function{2} = @(s,t) coefs(1) * (4*imag(deval(s,t))) ...
+                                         + coefs(2) * (4*imag(seval(s,t)));
+        obj.splitinfo.function{3} = @(s,t) coefs(1)*ones(size(t.r,2),size(s.r,2));
 
     case {'cp', 'cprime'}
         if ( nargin < 3 )
diff --git a/chunkie/@kernel/lap2d.m b/chunkie/@kernel/lap2d.m
index e7317bd..cb32a96 100644
--- a/chunkie/@kernel/lap2d.m
+++ b/chunkie/@kernel/lap2d.m
@@ -42,12 +42,21 @@
         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.sing = 'log';
+        obj.splitinfo = [];
+        obj.splitinfo.type = {[1 0 0 0]};
+        obj.splitinfo.action = {'r'};
+        obj.splitinfo.function = {@(s,t) ones(size(t.r,2),size(s.r,2))};
 
     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.function = {@(s,t) ones(size(t.r,2),size(s.r,2))};
 
     case {'sp', 'sprime'}
         obj.type = 'sp';
diff --git a/chunkie/chunkerkerneval.m b/chunkie/chunkerkerneval.m
index 3383ee1..18ae5a4 100644
--- a/chunkie/chunkerkerneval.m
+++ b/chunkie/chunkerkerneval.m
@@ -112,6 +112,7 @@
 if isfield(opts,'forcesmooth'); opts_use.forcesmooth = opts.forcesmooth; end
 if isfield(opts,'forceadap'); opts_use.forceadap = opts.forceadap; end
 if isfield(opts,'forcepquad'); opts_use.forcepquad = opts.forcepquad; end
+if isfield(opts,'side'); opts_use.side = opts.side; end
 if isfield(opts,'flam')
     opts_use.flam = opts.flam;
 end
@@ -251,7 +252,12 @@
     fints(irow0) = fints(irow0) + chnk.chunkerkerneval_smooth(chnkr0,kern0,opdims0,dens0,targinfo0, ...
         flag,opts_use);
 
-    fints(irow0) = fints(irow0) + chunkerkerneval_ho(chnkr0,kern0,opdims0,dens0, ...
+    if ~isfield(opts_use,'side')
+        msg = "Error: for pquad, set opts.side to 'i' or 'e' ";
+        error(msg)  
+    end
+
+    fints(irow0) = fints(irow0) + chunkerkerneval_pquad(chnkr0,kern0,opdims0,dens0, ...
         targinfo0,flag,opts_use);
 
     return
@@ -279,9 +285,13 @@
 
 end
 
-function fints = chunkerkerneval_ho(chnkr,kern,opdims,dens, ...
+function fints = chunkerkerneval_pquad(chnkr,kern,opdims,dens, ...
     targinfo,flag,opts)
 
+if ~isa(kern,'kernel') || isempty(kern.splitinfo)
+    error('Helsing-Ojala quad only available for kernel class objects with splitinfo defined');
+end
+
 % target
 [~,nt] = size(targinfo.r);
 fints = zeros(opdims(1)*nt,1);
@@ -293,6 +303,7 @@
 d = chnkr.d;
 n = chnkr.n;
 d2 = chnkr.d2;
+wts = chnkr.wts;
 
 % interpolation matrix
 intp = lege.matrin(k,t);          % interpolation from k to 2*k
@@ -300,30 +311,49 @@
 
 targs = targinfo.r;
 
-targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt);
-targn = zeros(chnkr.dim,nt);
-if isfield(targinfo, 'd')
-    targd = targinfo.d;
-end
-
-if isfield(targinfo, 'd2')
-    targd2 = targinfo.d2;
-end
-
-if isfield(targinfo, 'n')
-    targn = targinfo.n;
-end
 for j=1:size(chnkr.r,3)
     [ji] = find(flag(:,j));
     if(~isempty(ji))
         idxjmat = (j-1)*k+(1:k);
 
+        targinfoji = [];
+        targinfoji.r = targinfo.r(:,ji);
+        if isfield(targinfo, 'd')
+            targinfoji.d = targinfo.d(:,ji);
+        end
 
-        % Helsing-Ojala (interior/exterior?)
-        mat1 = chnk.pquadwts(r,d,n,d2,ct,bw,j,targs(:,ji), ...
-            targd(:,ji),targn(:,ji),targd2(:,ji),kern,opdims,t,w,opts,intp_ab,intp); % depends on kern, different mat1?
+        if isfield(targinfo, 'd2')
+            targinfoji.d2 = targinfo.d2(:,ji);
+        end
+
+        if isfield(targinfo, 'n')
+            targinfoji.n = targinfo.n(:,ji);
+        end        
+
+        srcinfo = [];
+        srcinfo.r = r(:,:,j);
+        srcinfo.d = d(:,:,j);
+        srcinfo.d2 = d2(:,:,j);
+        srcinfo.n = n(:,:,j);
 
-        fints(ji) = fints(ji) + mat1*dens(idxjmat);
+        % Helsing-Ojala (interior/exterior?)
+        allmats = cell(size(kern.splitinfo.type));
+        [allmats{:}] = chnk.pquadwts(r,d,n,d2,wts,j,targs(:,ji), ...
+              t,w,opts,intp_ab,intp,kern.splitinfo.type);
+    
+%        fints(ji) = fints(ji) + matcfield*dens(idxjmat);
+        for l = 1:length(allmats)
+            switch kern.splitinfo.action{l}
+                case 'r'
+                    mat0 = real(allmats{l});
+                case 'i'
+                    mat0 = imag(allmats{l});
+                case 'c'
+                    mat0 = allmats{l};
+            end
+            fints(ji) = fints(ji) + (kron(mat0,ones(opdims')).* ...
+                kern.splitinfo.function{l}(srcinfo,targinfoji))*dens(idxjmat);
+        end
     end
 end
 end
diff --git a/chunkie/chunkerkernevalmat.m b/chunkie/chunkerkernevalmat.m
index 7a5f593..94a547f 100644
--- a/chunkie/chunkerkernevalmat.m
+++ b/chunkie/chunkerkernevalmat.m
@@ -155,8 +155,8 @@
 if forcepquad
     optsflag = []; optsflag.fac = fac;
     flag = flagnear(chnkr,targinfo.r,optsflag);
-    spmat = chunkerkernevalmat_ho(chnkr,ftmp,opdims, ...
-        targinfo,flag,optsadap);
+    spmat = chunkerkernevalmat_pquad(chnkr,kern,opdims, ...
+        targinfo,flag,opts);
     mat = chunkerkernevalmat_smooth(chnkr,ftmp,opdims,targinfo, ...
         flag,opts);
     mat = mat + spmat;
@@ -324,7 +324,7 @@
 
 end
 
-function mat = chunkerkernevalmat_ho(chnkr,kern,opdims, ...
+function mat = chunkerkernevalmat_pquad(chnkr,kern,opdims, ...
     targinfo,flag,opts)
 
 k = chnkr.k;
@@ -372,6 +372,7 @@
     d = chnkr.d;
     n = chnkr.n;
     d2 = chnkr.d2;
+    wts = chnkr.wts;
 
     % interpolation matrix 
     intp = lege.matrin(k,t);          % interpolation from k to 2*k
@@ -383,10 +384,34 @@
                         
         [ji] = find(flag(:,i));
 
+        targinfoji = [];
+        targinfoji.r = targinfo.r(:,ji);
+
+        srcinfo = [];
+        srcinfo.r = r(:,:,i);
+        srcinfo.d = d(:,:,i);
+        srcinfo.d2 = d2(:,:,i);
+        srcinfo.n = n(:,:,i);
+
         % Helsing-Ojala (interior/exterior?)
-        mat1 = chnk.pquadwts(r,d,n,d2,ct,bw,i,targs(:,ji), ...
-                    targd(:,ji),targn(:,ji),targd2(:,ji),kern,opdims,t,w,opts,intp_ab,intp); % depends on kern, different mat1?
-                
+        allmats = cell(size(kern.splitinfo.type));
+        [allmats{:}] = chnk.pquadwts(r,d,n,d2,wts,i,targs(:,ji), ...
+              t,w,opts,intp_ab,intp,kern.splitinfo.type);
+
+        mat1 = zeros(size(allmats{1}));
+        for l = 1:length(allmats)
+            switch kern.splitinfo.action{l}
+                case 'r'
+                    mat0 = real(allmats{l});
+                case 'i'
+                    mat0 = imag(allmats{l});
+                case 'c'
+                    mat0 = allmats{l};
+            end
+            mat1 = mat1 + (kron(mat0,ones(opdims)).* ...
+                kern.splitinfo.function{l}(srcinfo,targinfoji));
+        end
+
         js1 = jmat:jmatend;
         js1 = repmat( (js1(:)).',opdims(1)*numel(ji),1);
         
diff --git a/devtools/test/pquadTest.m b/devtools/test/pquadTest.m
new file mode 100644
index 0000000..cc41dfc
--- /dev/null
+++ b/devtools/test/pquadTest.m
@@ -0,0 +1,115 @@
+%DEMO_BARBELL_2
+%
+% Define a polygonal "barbell" shaped domain with 
+% prescribed constant boundary data on each edge. 
+% Solves the corresponding Helmholtz Dirichlet problem
+% using corner rounding and smoothing the boundary data 
+% across the rounded corners
+%
+
+clearvars; close all;
+iseed = 8675309;
+rng(iseed,'twister');
+addpaths_loc();
+
+% planewave vec
+
+kvec = 10*[1;-1.5];
+
+%
+
+zk = norm(kvec);
+
+% define geometry and boundary conditions
+% (vertices defined by another function)
+verts = chnk.demo.barbell(2.0,2.0,1.0,1.0);
+nv = size(verts,2);
+edgevals = randn(1,nv);
+
+% parameters for curve rounding/chunking routine to oversample boundary
+cparams = [];
+cparams.widths = 0.1*ones(size(verts,2),1);% width to cut about each corner
+cparams.eps = 1e-12; % tolerance at which to resolve curve
+cparams.rounded = true;
+
+% call smoothed-polygon chunking routine
+% a smoothed version of edgevals is returned in 
+% chnkr.data
+chnkr = chunkerpoly(verts,cparams,[],edgevals);
+chnkr = refine(chnkr,struct('nover',1));
+%
+% plot smoothed geometry and data
+
+figure(1)
+clf
+plot(chnkr,'-x')
+hold on
+quiver(chnkr)
+axis equal
+
+figure(2)
+clf
+nchplot = 1:chnkr.nch;
+plot3(chnkr,1)
+
+% solve and visualize the solution
+
+% build laplace dirichlet matrix
+
+% fkern = kernel('laplace','d');
+fkern = kernel('helmholtz','c',zk,[1.5,-2i]);
+opts = [];
+start = tic; C = chunkermat(chnkr,fkern,opts);
+t1 = toc(start);
+
+fprintf('%5.2e s : time to assemble matrix\n',t1)
+
+%
+
+sys = -0.5*eye(chnkr.npt) + C;
+
+% rhs = chnkr.data(1,:); rhs = rhs(:);
+rhs = besselh(0,zk*abs((1+1.25*1i)-(chnkr.r(1,:)+1i*chnkr.r(2,:)))); rhs = rhs(:);
+start = tic; sol = gmres(sys,rhs,[],1e-14,100); t1 = toc(start);
+
+fprintf('%5.2e s : time for dense gmres\n',t1)
+
+% evaluate at targets and plot
+
+rmin = min(chnkr); rmax = max(chnkr);
+nplot = 200;
+hx = (rmax(1)-rmin(1))/nplot;
+hy = (rmax(2)-rmin(2))/nplot;
+xtarg = linspace(rmin(1)+hx/2,rmax(1)-hx/2,nplot); 
+ytarg = linspace(rmin(2)+hy/2,rmax(2)-hy/2,nplot);
+[xxtarg,yytarg] = meshgrid(xtarg,ytarg);
+targets = zeros(2,length(xxtarg(:)));
+targets(1,:) = xxtarg(:); targets(2,:) = yytarg(:);
+
+start = tic; in = chunkerinterior(chnkr,targets); t1 = toc(start);
+fprintf('%5.2e s : time to find points in domain\n',t1)
+
+% compute layer potential at interior points
+opts.forcepquad = true; 
+opts.side = 'i'; % 'i' for interior, 'e' for exterior, for positively oriented curve.
+start = tic;
+Csolpquad = chunkerkerneval(chnkr,fkern,sol,targets(:,in),opts); 
+% below check slp and dlp
+% fkernd = kernel('helmholtz','d',zk);
+% fkerns = kernel('helmholtz','s',zk);
+% Ssolpquad = chunkerkerneval(chnkr,fkerns,sol,targets(:,in),opts); 
+% Dsolpquad = chunkerkerneval(chnkr,fkernd,sol,targets(:,in),opts); 
+% Csolpquad = Dsolpquad + (-zk*1i)*Ssolpquad;
+t1 = toc(start);
+fprintf('%5.2e s : time for kerneval (Helsing-Ojala for near)\n',t1);
+
+start = tic;
+Csol = chunkerkerneval(chnkr,fkern,sol,targets(:,in)); t1 = toc(start);
+fprintf('%5.2e s : time for kerneval (adaptive for near)\n',t1);
+
+% Compare with reference solution Dsol
+error = max(abs(Csol-Csolpquad))/max(abs(Csol));
+fprintf('%5.2e : Relative max error\n',error);
+% 
+
+assert(error < 1e-10)
\ No newline at end of file

From 40db4aedaf5d341ceec582076367ca0853b08f89 Mon Sep 17 00:00:00 2001
From: Hai Zhu <zhuh1106@gmail.com>
Date: Thu, 17 Oct 2024 16:04:20 -0700
Subject: [PATCH 2/8] minor performance tweak, add some comments

---
 chunkie/+chnk/pquadwts.m  | 93 +++++++++++++++++++++++++++++++++++++--
 chunkie/@kernel/lap2d.m   |  1 -
 chunkie/chunkerkerneval.m | 14 +++---
 devtools/test/pquadTest.m |  6 ++-
 4 files changed, 103 insertions(+), 11 deletions(-)

diff --git a/chunkie/+chnk/pquadwts.m b/chunkie/+chnk/pquadwts.m
index 8edccfd..7b892ae 100644
--- a/chunkie/+chnk/pquadwts.m
+++ b/chunkie/+chnk/pquadwts.m
@@ -24,6 +24,9 @@
 %   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
 %
 % Output
 %   varargout - integration matrices for specified singularity types
@@ -42,17 +45,26 @@
 
 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;
+
 for j = 1:length(types)
     type0 = types{j};
 
     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} = 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} = 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
@@ -213,3 +225,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
diff --git a/chunkie/@kernel/lap2d.m b/chunkie/@kernel/lap2d.m
index cb32a96..9f4962e 100644
--- a/chunkie/@kernel/lap2d.m
+++ b/chunkie/@kernel/lap2d.m
@@ -42,7 +42,6 @@
         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.sing = 'log';
         obj.splitinfo = [];
         obj.splitinfo.type = {[1 0 0 0]};
         obj.splitinfo.action = {'r'};
diff --git a/chunkie/chunkerkerneval.m b/chunkie/chunkerkerneval.m
index 18ae5a4..a53a4e8 100644
--- a/chunkie/chunkerkerneval.m
+++ b/chunkie/chunkerkerneval.m
@@ -344,15 +344,19 @@
 %        fints(ji) = fints(ji) + matcfield*dens(idxjmat);
         for l = 1:length(allmats)
             switch kern.splitinfo.action{l}
-                case 'r'
+                case 'r' % real part
                     mat0 = real(allmats{l});
-                case 'i'
+                case 'i' % imaginary part
                     mat0 = imag(allmats{l});
-                case 'c'
+                case 'c' % complex
                     mat0 = allmats{l};
             end
-            fints(ji) = fints(ji) + (kron(mat0,ones(opdims')).* ...
-                kern.splitinfo.function{l}(srcinfo,targinfoji))*dens(idxjmat);
+            mat0opdim = kron(mat0,ones(opdims'));
+            % this is reevaluating more or less the same kernel multiple times...
+            % think of a better kernel split format...
+            mat0xsplitfun = mat0opdim.* ...
+                            kern.splitinfo.function{l}(srcinfo,targinfoji);
+            fints(ji) = fints(ji) + mat0xsplitfun*dens(idxjmat);
         end
     end
 end
diff --git a/devtools/test/pquadTest.m b/devtools/test/pquadTest.m
index cc41dfc..891c778 100644
--- a/devtools/test/pquadTest.m
+++ b/devtools/test/pquadTest.m
@@ -7,6 +7,9 @@
 % across the rounded corners
 %
 
+% profile clear
+% profile on
+
 clearvars; close all;
 iseed = 8675309;
 rng(iseed,'twister');
@@ -112,4 +115,5 @@
 fprintf('%5.2e : Relative max error\n',error);
 % 
 
-assert(error < 1e-10)
\ No newline at end of file
+assert(error < 1e-10)
+% profile viewer
\ No newline at end of file

From ce774f7a6e3b62a03db8bf1a5cad2965c5c24761 Mon Sep 17 00:00:00 2001
From: Hai Zhu <zhuh1106@gmail.com>
Date: Fri, 18 Oct 2024 08:02:48 -0700
Subject: [PATCH 3/8] rewrite kernel split into function to efficiently reuse
 precomputed info, 1st for helmholtz combined field

---
 chunkie/@kernel/helm2d.m  | 26 +++++++++++++++-----------
 chunkie/chunkerkerneval.m |  7 ++-----
 devtools/test/pquadTest.m |  2 +-
 3 files changed, 18 insertions(+), 17 deletions(-)

diff --git a/chunkie/@kernel/helm2d.m b/chunkie/@kernel/helm2d.m
index 0ee5ec8..b72d72b 100644
--- a/chunkie/@kernel/helm2d.m
+++ b/chunkie/@kernel/helm2d.m
@@ -88,17 +88,7 @@
         obj.splitinfo = [];
         obj.splitinfo.type = {[0 0 0 0],[1 0 0 0],[0 0 -1 0]};
         obj.splitinfo.action = {'r','r','r'};
-        deval = @(s,t) chnk.helm2d.kern(zk, s, t, 'd', coefs(1));
-        seval = @(s,t) chnk.helm2d.kern(zk, s, t, 's', coefs(2));
-        obj.splitinfo.function = cell(3,1);
-        obj.splitinfo.function{1} = @(s,t) coefs(1) * (deval(s,t) ...
-                                                       + 1/(2*pi)*4*log(abs((s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)'))).*imag(deval(s,t)) ... 
-                                                       + 1/(2*pi)*real((s.n(1,:)+1i*s.n(2,:))./((s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)')))) ...
-                                         + coefs(2) * (seval(s,t) ...
-                                                       + 1/(2*pi)*4*log(abs((s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)'))).*imag(seval(s,t)));
-        obj.splitinfo.function{2} = @(s,t) coefs(1) * (4*imag(deval(s,t))) ...
-                                         + coefs(2) * (4*imag(seval(s,t)));
-        obj.splitinfo.function{3} = @(s,t) coefs(1)*ones(size(t.r,2),size(s.r,2));
+        obj.splitinfo.functions = @(s,t) helm2d_c(zk,s,t,coefs);
 
     case {'cp', 'cprime'}
         if ( nargin < 3 )
@@ -123,3 +113,17 @@
 
 
 end
+
+function f = helm2d_c(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
+
diff --git a/chunkie/chunkerkerneval.m b/chunkie/chunkerkerneval.m
index a53a4e8..550694d 100644
--- a/chunkie/chunkerkerneval.m
+++ b/chunkie/chunkerkerneval.m
@@ -341,7 +341,7 @@
         [allmats{:}] = chnk.pquadwts(r,d,n,d2,wts,j,targs(:,ji), ...
               t,w,opts,intp_ab,intp,kern.splitinfo.type);
     
-%        fints(ji) = fints(ji) + matcfield*dens(idxjmat);
+        funs = kern.splitinfo.functions(srcinfo,targinfoji);
         for l = 1:length(allmats)
             switch kern.splitinfo.action{l}
                 case 'r' % real part
@@ -352,10 +352,7 @@
                     mat0 = allmats{l};
             end
             mat0opdim = kron(mat0,ones(opdims'));
-            % this is reevaluating more or less the same kernel multiple times...
-            % think of a better kernel split format...
-            mat0xsplitfun = mat0opdim.* ...
-                            kern.splitinfo.function{l}(srcinfo,targinfoji);
+            mat0xsplitfun = mat0opdim.*funs{l};
             fints(ji) = fints(ji) + mat0xsplitfun*dens(idxjmat);
         end
     end
diff --git a/devtools/test/pquadTest.m b/devtools/test/pquadTest.m
index 891c778..aa6d9a1 100644
--- a/devtools/test/pquadTest.m
+++ b/devtools/test/pquadTest.m
@@ -80,7 +80,7 @@
 % evaluate at targets and plot
 
 rmin = min(chnkr); rmax = max(chnkr);
-nplot = 200;
+nplot = 300;
 hx = (rmax(1)-rmin(1))/nplot;
 hy = (rmax(2)-rmin(2))/nplot;
 xtarg = linspace(rmin(1)+hx/2,rmax(1)-hx/2,nplot); 

From 2a121c1f7f441f88a124ac5bc393a83539f01795 Mon Sep 17 00:00:00 2001
From: Hai Zhu <zhuh1106@gmail.com>
Date: Fri, 18 Oct 2024 09:13:06 -0700
Subject: [PATCH 4/8] add helm2d s_split and d_split

---
 chunkie/@kernel/helm2d.m | 36 ++++++++++++++++++++++++------------
 1 file changed, 24 insertions(+), 12 deletions(-)

diff --git a/chunkie/@kernel/helm2d.m b/chunkie/@kernel/helm2d.m
index b72d72b..80e3c2a 100644
--- a/chunkie/@kernel/helm2d.m
+++ b/chunkie/@kernel/helm2d.m
@@ -43,10 +43,7 @@
         obj.splitinfo = [];
         obj.splitinfo.type = {[0 0 0 0],[1 0 0 0]};
         obj.splitinfo.action = {'r','r'};
-        obj.splitinfo.function = cell(2,1);
-        obj.splitinfo.function{1} = @(s,t) obj.eval(s,t) ...
-                                         + 1/(2*pi)*4*log(abs((s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)'))).*imag(obj.eval(s,t));
-        obj.splitinfo.function{2} = @(s,t) 4*imag(obj.eval(s,t));
+        obj.splitinfo.functions = @(s,t) helm2d_s_split(zk,s,t);
 
     case {'d', 'double'}
         obj.type = 'd';
@@ -56,12 +53,7 @@
         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.function = cell(3,1);
-        obj.splitinfo.function{1} = @(s,t) obj.eval(s,t) ...
-                                           + 1/(2*pi)*4*log(abs((s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)'))).*imag(obj.eval(s,t)) ... 
-                                           + 1/(2*pi)*real((s.n(1,:)+1i*s.n(2,:))./((s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)')));
-        obj.splitinfo.function{2} = @(s,t) 4*imag(obj.eval(s,t));
-        obj.splitinfo.function{3} = @(s,t) ones(size(t.r,2),size(s.r,2));
+        obj.splitinfo.functions = @(s,t) helm2d_d_split(zk,s,t);
 
     case {'sp', 'sprime'}
         obj.type = 'sp';
@@ -88,7 +80,7 @@
         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(zk,s,t,coefs);
+        obj.splitinfo.functions = @(s,t) helm2d_c_split(zk,s,t,coefs);
 
     case {'cp', 'cprime'}
         if ( nargin < 3 )
@@ -114,7 +106,27 @@
 
 end
 
-function f = helm2d_c(zk,s,t,coefs)
+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,:)');

From 2cef27df5a7523658261ad427db3ab19975a5b98 Mon Sep 17 00:00:00 2001
From: Hai Zhu <zhuh1106@gmail.com>
Date: Fri, 18 Oct 2024 09:41:07 -0700
Subject: [PATCH 5/8] add lap2d kernel split functions

---
 chunkie/@kernel/lap2d.m                       | 25 +++++++++++++++++--
 chunkie/chunkerkernevalmat.m                  |  6 +++--
 .../chunkermat_quadadap_closetotouchingTest.m |  7 ++++--
 3 files changed, 32 insertions(+), 6 deletions(-)

diff --git a/chunkie/@kernel/lap2d.m b/chunkie/@kernel/lap2d.m
index 9f4962e..df7a5dd 100644
--- a/chunkie/@kernel/lap2d.m
+++ b/chunkie/@kernel/lap2d.m
@@ -45,7 +45,7 @@
         obj.splitinfo = [];
         obj.splitinfo.type = {[1 0 0 0]};
         obj.splitinfo.action = {'r'};
-        obj.splitinfo.function = {@(s,t) ones(size(t.r,2),size(s.r,2))};
+        obj.splitinfo.functions = @(s,t) lap2d_s_split(s,t);
 
     case {'d', 'double'}
         obj.type = 'd';
@@ -55,7 +55,7 @@
         obj.splitinfo = [];
         obj.splitinfo.type = {[0 0 -1 0]};
         obj.splitinfo.action = {'r'};
-        obj.splitinfo.function = {@(s,t) ones(size(t.r,2),size(s.r,2))};
+        obj.splitinfo.functions = @(s,t) lap2d_d_split(s,t);
 
     case {'sp', 'sprime'}
         obj.type = 'sp';
@@ -99,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);
@@ -111,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
+
diff --git a/chunkie/chunkerkernevalmat.m b/chunkie/chunkerkernevalmat.m
index 94a547f..5f163a8 100644
--- a/chunkie/chunkerkernevalmat.m
+++ b/chunkie/chunkerkernevalmat.m
@@ -399,6 +399,7 @@
               t,w,opts,intp_ab,intp,kern.splitinfo.type);
 
         mat1 = zeros(size(allmats{1}));
+        funs = kern.splitinfo.functions(srcinfo,targinfoji);
         for l = 1:length(allmats)
             switch kern.splitinfo.action{l}
                 case 'r'
@@ -408,8 +409,9 @@
                 case 'c'
                     mat0 = allmats{l};
             end
-            mat1 = mat1 + (kron(mat0,ones(opdims)).* ...
-                kern.splitinfo.function{l}(srcinfo,targinfoji));
+            mat0opdim = kron(mat0,ones(opdims));
+            mat0xsplitfun = mat0opdim.*funs{l};
+            mat1 = mat1 + mat0xsplitfun;
         end
 
         js1 = jmat:jmatend;
diff --git a/devtools/test/chunkermat_quadadap_closetotouchingTest.m b/devtools/test/chunkermat_quadadap_closetotouchingTest.m
index 450b309..3f1bafa 100644
--- a/devtools/test/chunkermat_quadadap_closetotouchingTest.m
+++ b/devtools/test/chunkermat_quadadap_closetotouchingTest.m
@@ -72,20 +72,23 @@
 % use adaptive routine to build matrix (self done by ggq, nbor by adaptive)
 
 eta = 1;
-fkern = @(s,t) chnk.lap2d.kern(s,t,'C',[1,eta]);
+fkern = kernel('laplace','c',[1,eta]);
 
 type = 'log';
 opts = []; opts.robust = true;
+opts.adaptive_correction = true;
 opdims = [1 1];
 start = tic;
-mata = chnk.quadadap.buildmat(chnkr,fkern,opdims,type,opts);
+mata = chunkermat(chnkr,fkern,opts);
 t1 = toc(start);
 fprintf('%5.2e s : time to assemble matrix (adaptive)\n',t1)
+opts.adaptive_correction = false;
 start = tic; mato = chunkermat(chnkr,fkern);
 t1 = toc(start);
 fprintf('%5.2e s : time to assemble matrix (original GGQ)\n',t1)
 start = tic; 
 opts.forcepquad = true; targs_bd = [chnkr.r(1,:);chnkr.r(2,:)];
+opts.side = 'e';
 matho = chunkerkernevalmat(chnkr,fkern,targs_bd,opts); 
 t1 = toc(start);
 fprintf('%5.2e s : time to assemble matrix (Helsing-Ojala)\n',t1)

From 0c9023d5528ea4f9b4ac4aa6f5f869e41c6aa31d Mon Sep 17 00:00:00 2001
From: Hai Zhu <zhuh1106@gmail.com>
Date: Fri, 18 Oct 2024 09:53:21 -0700
Subject: [PATCH 6/8] add comment

---
 chunkie/+chnk/pquadwts.m | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/chunkie/+chnk/pquadwts.m b/chunkie/+chnk/pquadwts.m
index 7b892ae..242c8fa 100644
--- a/chunkie/+chnk/pquadwts.m
+++ b/chunkie/+chnk/pquadwts.m
@@ -27,6 +27,8 @@
 %       [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
 %   varargout - integration matrices for specified singularity types

From 0d3617cf561e701f3562f6a0ec6b13da27df3172 Mon Sep 17 00:00:00 2001
From: Hai Zhu <zhuh1106@gmail.com>
Date: Fri, 18 Oct 2024 22:08:05 -0700
Subject: [PATCH 7/8] stokes slp velocity kernel split

---
 chunkie/@kernel/stok2d.m              | 24 ++++++++++++++++++++++++
 chunkie/chunkerkerneval.m             | 11 +++++++++--
 devtools/test/chunkermat_stok2dTest.m | 14 ++++++++++++++
 3 files changed, 47 insertions(+), 2 deletions(-)

diff --git a/chunkie/@kernel/stok2d.m b/chunkie/@kernel/stok2d.m
index 236cbdb..91ae858 100644
--- a/chunkie/@kernel/stok2d.m
+++ b/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,23 @@
 
 
 end
+
+function f = stok2d_s_split(mu,s,t)
+dist = (s.r(1,:)+1i*s.r(2,:))-(t.r(1,:)'+1i*t.r(2,:)');
+f = cell(3, 1);
+ntarg = numel(t.r(1,:));
+nsrc = numel(s.r(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) = (-(repmat(s.n(1,:),[ntarg 1]))).*real(dist)/(2*mu);     
+f{2}(1:2:end,2:2:end) = -((repmat(s.n(2,:),[ntarg 1]))).*real(dist)/(2*mu);
+f{2}(2:2:end,1:2:end) = (-(repmat(s.n(1,:),[ntarg 1]))).*imag(dist)/(2*mu); 
+f{2}(2:2:end,2:2:end) = -((repmat(s.n(2,:),[ntarg 1]))).*imag(dist)/(2*mu);
+f{3}(1:2:end,1:2:end) = - (repmat(s.n(2,:),[ntarg 1])).*real(dist)/(2*mu);
+f{3}(1:2:end,2:2:end) =   (repmat(s.n(1,:),[ntarg 1])).*real(dist)/(2*mu);
+f{3}(2:2:end,1:2:end) = - (repmat(s.n(2,:),[ntarg 1])).*imag(dist)/(2*mu);
+f{3}(2:2:end,2:2:end) =   (repmat(s.n(1,:),[ntarg 1])).*imag(dist)/(2*mu);
+end
diff --git a/chunkie/chunkerkerneval.m b/chunkie/chunkerkerneval.m
index 550694d..0cb8013 100644
--- a/chunkie/chunkerkerneval.m
+++ b/chunkie/chunkerkerneval.m
@@ -295,7 +295,8 @@
 % target
 [~,nt] = size(targinfo.r);
 fints = zeros(opdims(1)*nt,1);
-k = size(chnkr.r,2);
+k = chnkr.k;
+nch = chnkr.nch;
 [t,w] = lege.exps(2*k);
 ct = lege.exps(k);
 bw = lege.barywts(k);
@@ -305,6 +306,9 @@
 d2 = chnkr.d2;
 wts = chnkr.wts;
 
+assert(numel(dens) == opdims(2)*k*nch,'dens not of appropriate size')
+dens = reshape(dens,opdims(2),k,nch);
+
 % interpolation matrix
 intp = lege.matrin(k,t);          % interpolation from k to 2*k
 intp_ab = lege.matrin(k,[-1;1]);  % interpolation from k to end points
@@ -335,6 +339,9 @@
         srcinfo.d = d(:,:,j);
         srcinfo.d2 = d2(:,:,j);
         srcinfo.n = n(:,:,j);
+        densj = dens(:,:,j);
+        indji = (ji-1)*opdims(1);
+        ind = (indji(:)).' + (1:opdims(1)).';
 
         % Helsing-Ojala (interior/exterior?)
         allmats = cell(size(kern.splitinfo.type));
@@ -353,7 +360,7 @@
             end
             mat0opdim = kron(mat0,ones(opdims'));
             mat0xsplitfun = mat0opdim.*funs{l};
-            fints(ji) = fints(ji) + mat0xsplitfun*dens(idxjmat);
+            fints(ind(:)) = fints(ind(:)) + mat0xsplitfun*densj(:);
         end
     end
 end
diff --git a/devtools/test/chunkermat_stok2dTest.m b/devtools/test/chunkermat_stok2dTest.m
index cc18225..6e22f47 100644
--- a/devtools/test/chunkermat_stok2dTest.m
+++ b/devtools/test/chunkermat_stok2dTest.m
@@ -38,6 +38,7 @@
 ts = 0.0+2*pi*rand(nt,1);
 targets = starfish(ts,narms,amp);
 targets = targets.*repmat(rand(1,nt),2,1)*0.75;
+% targets = targets.*repmat(rand(1,nt),2,1)*0.99;
 
 plot(chnkr, 'r.'); hold on;
 plot(targets(1,:), targets(2,:), 'kx')
@@ -93,6 +94,19 @@
 t1 = toc(start);
 fprintf('%5.2e s : time to eval at targs (slow, adaptive routine)\n',t1)
 
+% test stokes slp velocity pquad...
+% targets = targets.*repmat(rand(1,nt),2,1)*0.99 seems to break below stokes kernel test...
+fkerns = kernel('stok', 'svel', mu);
+Ssol = chunkerkerneval(chnkr,fkerns,sol,targets,opts); 
+opts.forcepquad=true;
+opts.side = 'i';
+Ssol_pquad = chunkerkerneval(chnkr,fkerns,sol,targets,opts); 
+err = abs(Ssol - Ssol_pquad);
+assert(norm(err) < 1e-10);
+% figure(1),clf,
+% plot(chnkr), hold on
+% scatter(targets(1,:),targets(2,:),[],log10(err(1:2:end))); colorbar
+opts.forcepquad=false;
 %
 
 wchnkr = chnkr.wts;

From b7e0433169474656f3027a3c114517d9839baa63 Mon Sep 17 00:00:00 2001
From: Hai Zhu <zhuh1106@gmail.com>
Date: Fri, 18 Oct 2024 22:18:13 -0700
Subject: [PATCH 8/8] clean up stokes slp kernel split

---
 chunkie/@kernel/stok2d.m              | 20 ++++++++++++--------
 devtools/test/chunkermat_stok2dTest.m |  2 +-
 2 files changed, 13 insertions(+), 9 deletions(-)

diff --git a/chunkie/@kernel/stok2d.m b/chunkie/@kernel/stok2d.m
index 91ae858..b5bb102 100644
--- a/chunkie/@kernel/stok2d.m
+++ b/chunkie/@kernel/stok2d.m
@@ -195,20 +195,24 @@
 
 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) = (-(repmat(s.n(1,:),[ntarg 1]))).*real(dist)/(2*mu);     
-f{2}(1:2:end,2:2:end) = -((repmat(s.n(2,:),[ntarg 1]))).*real(dist)/(2*mu);
-f{2}(2:2:end,1:2:end) = (-(repmat(s.n(1,:),[ntarg 1]))).*imag(dist)/(2*mu); 
-f{2}(2:2:end,2:2:end) = -((repmat(s.n(2,:),[ntarg 1]))).*imag(dist)/(2*mu);
-f{3}(1:2:end,1:2:end) = - (repmat(s.n(2,:),[ntarg 1])).*real(dist)/(2*mu);
-f{3}(1:2:end,2:2:end) =   (repmat(s.n(1,:),[ntarg 1])).*real(dist)/(2*mu);
-f{3}(2:2:end,1:2:end) = - (repmat(s.n(2,:),[ntarg 1])).*imag(dist)/(2*mu);
-f{3}(2:2:end,2:2:end) =   (repmat(s.n(1,:),[ntarg 1])).*imag(dist)/(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
diff --git a/devtools/test/chunkermat_stok2dTest.m b/devtools/test/chunkermat_stok2dTest.m
index 6e22f47..0d1b2ed 100644
--- a/devtools/test/chunkermat_stok2dTest.m
+++ b/devtools/test/chunkermat_stok2dTest.m
@@ -94,7 +94,7 @@
 t1 = toc(start);
 fprintf('%5.2e s : time to eval at targs (slow, adaptive routine)\n',t1)
 
-% test stokes slp velocity pquad...
+% just test stokes slp velocity pquad...
 % targets = targets.*repmat(rand(1,nt),2,1)*0.99 seems to break below stokes kernel test...
 fkerns = kernel('stok', 'svel', mu);
 Ssol = chunkerkerneval(chnkr,fkerns,sol,targets,opts);