From 7503c39d9a64fddb79c3fb8ade81ad70018dcbf5 Mon Sep 17 00:00:00 2001
From: Travis Askham <askhamwhat@gmail.com>
Date: Mon, 6 Jan 2025 14:20:47 -0500
Subject: [PATCH 1/8] small local fixes

---
 chunkie/chunkerfunc.m        | 13 +++++++------
 chunkie/chunkerkernevalmat.m | 15 ++++++++-------
 chunkie/chunkermat.m         | 14 +++++++++++---
 3 files changed, 26 insertions(+), 16 deletions(-)

diff --git a/chunkie/chunkerfunc.m b/chunkie/chunkerfunc.m
index dfc27949..961c5549 100644
--- a/chunkie/chunkerfunc.m
+++ b/chunkie/chunkerfunc.m
@@ -155,7 +155,7 @@
     tsplits = [ta;tb];
 end
    
-tsplits = sort(unique(tsplits),'ascend');
+tsplits = sort(uniquetol(tsplits,eps),'ascend');
 lab = length(tsplits);
 if (lab-1 > nchmax)
     error(['CHUNKERFUNC: nchmax exceeded in chunkerfunc on initial splits.\n ',...
@@ -184,14 +184,15 @@
 
 maxiter_res=nchmax-nch;
 
+xmin =  Inf;
+xmax = -Inf;
+ymin =  Inf;
+ymax = -Inf;
+
 rad_curr = 0;
 for ijk = 1:maxiter_res
 
 %       loop through all existing chunks, if resolved store, if not split
-    xmin =  Inf;
-    xmax = -Inf;
-    ymin =  Inf;
-    ymax = -Inf;
     
     ifdone=1;
     for ich=1:nchnew
@@ -224,7 +225,7 @@
             
             resol_speed_test = err1>eps;
             if nout < 2
-                resol_speed_test = err1>eps*k;
+                resol_speed_test = err1*(b-a) > eps*k;
             end
             
             xmax = max(xmax,max(r(1,:)));
diff --git a/chunkie/chunkerkernevalmat.m b/chunkie/chunkerkernevalmat.m
index 7a5f5934..ae2ed5f1 100644
--- a/chunkie/chunkerkernevalmat.m
+++ b/chunkie/chunkerkernevalmat.m
@@ -1,4 +1,4 @@
-    function mat = chunkerkernevalmat(chnkr,kern,targobj,opts)
+function [mat,varargout] = chunkerkernevalmat(chnkr,kern,targobj,opts)
 %CHUNKERKERNEVALMAT compute the matrix which maps density values on 
 % the chunk geometry to the value of the convolution of the given
 % integral kernel with the density at the specified target points
@@ -138,8 +138,6 @@
 optsadap = []; 
 optsadap.eps = eps;
 
-
-
 if forcesmooth
     mat = chunkerkernevalmat_smooth(chnkr,ftmp,opdims,targinfo, ...
         [],optssmooth);
@@ -165,11 +163,14 @@
 
 % smooth for sufficiently far, adaptive otherwise
 
-% TODO: change to chunkerkerneval system, need routine to generate
-% upsampling matrix.
+rho = 1.8;
+optsflag = [];  optsflag.rho = rho;
+flag = flagnear_rectangle(chnkr,targinfo.r,optsflag);
+
+npoly = chnkr.k*2;
+nlegnew = chnk.ellipse_oversample(rho,npoly,eps);
+nlegnew = max(nlegnew,chnkr.k);
 
-optsflag = []; optsflag.fac = fac;
-flag = flagnear(chnkr,targinfo.r,optsflag);
 spmat = chunkerkernevalmat_adap(chnkr,ftmp,opdims, ...
         targinfo,flag,optsadap);
 
diff --git a/chunkie/chunkermat.m b/chunkie/chunkermat.m
index ad49416f..ef5d6389 100644
--- a/chunkie/chunkermat.m
+++ b/chunkie/chunkermat.m
@@ -411,7 +411,9 @@
 
         % mark off the near and self interactions
         for ich = 1:chnkr.nch
-            for jch = [ich,chnkr.adj(1,ich),chnkr.adj(2,ich)]
+            jlist = [ich,chnkr.adj(1,ich),chnkr.adj(2,ich)];
+            jlist = jlist(jlist > 0);
+            for jch = jlist
                 flag((jch - 1)*chnkr.k+(1:chnkr.k), ich) = 0;
             end
         end
@@ -453,6 +455,7 @@
     npt_all = horzcat(chnkobj.echnks.npt);
     [~,nv] = size(chnkobj.verts);
     ngl = chnkrs(1).k;
+    rcipsav = cell(nv,1);
     
     for ivert=1:nv
         clist = chnkobj.vstruc{ivert}{1};
@@ -500,10 +503,11 @@
         optsrcip.nonsmoothonly = false;
         optsrcip.corrections = false;
 
-        R = chnk.rcip.Rcompchunk(chnkrs,iedgechunks,kern,ndim, ...
+        [R,rcipsav{ivert}] = chnk.rcip.Rcompchunk(chnkrs,iedgechunks,kern,ndim, ...
             Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,starL1,circL1,... 
             sbclmat,sbcrmat,lvmat,rvmat,u,optsrcip);
-       
+        rcipsav{ivert}.starind = starind;
+        
         sysmat_tmp = inv(R) - eye(2*ngl*nedge*ndim);
         if (~nonsmoothonly)
             
@@ -561,6 +565,10 @@
             vsysmat = [vsysmat;sysmat_tmp(:)];
         end    
     end
+
+    if nargout > 2
+        varargout{2} = rcipsav;
+    end
     
 end
 

From dbf992d55698db0759752c264a7917068575f35d Mon Sep 17 00:00:00 2001
From: Travis Askham <askhamwhat@gmail.com>
Date: Mon, 2 Jun 2025 16:10:01 -0400
Subject: [PATCH 2/8] update rcip with save depth option

---
 chunkie/+chnk/+rcip/Rcompchunk.m   |  26 +++----
 chunkie/+chnk/+rcip/rhohatInterp.m | 110 ++++++++++++++++-------------
 chunkie/chunkermat.m               |  11 ++-
 3 files changed, 84 insertions(+), 63 deletions(-)

diff --git a/chunkie/+chnk/+rcip/Rcompchunk.m b/chunkie/+chnk/+rcip/Rcompchunk.m
index 0d2de0f3..d3af4f76 100644
--- a/chunkie/+chnk/+rcip/Rcompchunk.m
+++ b/chunkie/+chnk/+rcip/Rcompchunk.m
@@ -53,19 +53,18 @@
     opts = [];
 end
 
-savedeep = false;
-if isfield(opts,'savedeep')
-    savedeep = opts.savedeep;
+savedepth = 10;
+if isfield(opts,'rcip_savedepth')
+    savedepth = opts.rcip_savedepth;
 end
+savedepth = max(savedepth,0);
+savedepth = min(savedepth,nsub);
 
-rcipsav.savedeep = savedeep;
-
-if savedeep
-    rcipsav.R = cell(nsub+1,1);
-    rcipsav.MAT = cell(nsub,1);
-    rcipsav.chnkrlocals = cell(nsub,1);
-end
+rcipsav.savedepth = savedepth;
 
+rcipsav.R = cell(nsub+1,1);
+rcipsav.MAT = cell(nsub,1);
+rcipsav.chnkrlocals = cell(nsub,1);
 
 % grab only those kernels relevant to this vertex
 
@@ -260,12 +259,15 @@
     if level==1    %  Dumb and lazy initializer for R, for now
   %R=eye(nR); 
         R = inv(MAT(starL,starL));
-        if savedeep
+        if level >= nsub-savedepth+1
             rcipsav.R{1} = R;
         end
     end
+    if savedepth < nsub && level == nsub-savedepth+1
+        rcipsav.R{level} = R;
+    end
     R=chnk.rcip.SchurBana(Pbc,PWbc,MAT,R,starL,circL,starS,circS);   
-    if savedeep
+    if level >= nsub-savedepth+1
         rcipsav.R{level+1} = R;
         rcipsav.MAT{level} = MAT(starL,circL);
         rcipsav.chnkrlocals{level} = merge(chnkrlocal);
diff --git a/chunkie/+chnk/+rcip/rhohatInterp.m b/chunkie/+chnk/+rcip/rhohatInterp.m
index 1ed8b26f..b10a58b2 100644
--- a/chunkie/+chnk/+rcip/rhohatInterp.m
+++ b/chunkie/+chnk/+rcip/rhohatInterp.m
@@ -3,9 +3,9 @@
 % density rhohat to the requested depth using the backward recursion
 %
 % When using an RCIP preconditioner (in the style of eq 34 of the RCIP 
-% tutorial), the resulting coarse level density is 
-% accurate for interactions separated from the vertex 
-% *at the coarse scale*. 
+% tutorial), the resulting coarse level density is accurate for
+% interactions separated from the vertex *at the coarse scale*. 
+%
 % By running the recursion for the RCIP preconditioner backwards, an 
 % equivalent density can be reconstructed on the fine mesh which is 
 % appropriate for closer interactions (at a distance well-separated from
@@ -53,26 +53,40 @@
 % wts - a set of weights for integrating functions sampled at these
 %           points
 
-				% author: Travis Askham
+% author: Travis Askham
   
-rhohatinterp = [];
-srcinfo = [];
-wts = [];
-
-k = rcipsav.k;
-ndim = rcipsav.ndim;
 nedge = rcipsav.nedge;
 Pbc = rcipsav.Pbc;
-PWbc = rcipsav.PWbc;
-starL = rcipsav.starL;
-starL1 = rcipsav.starL1;
-starS = rcipsav.starS;
-circL = rcipsav.circL;
-circL1 = rcipsav.circL1;
-circS = rcipsav.circS;
-ilist = rcipsav.ilist;
+
+starL1 = sort(rcipsav.starL1);
+starS = sort(rcipsav.starS);
+
+circL1 = sort(rcipsav.circL1);
+circS = sort(rcipsav.circS);
+
+nrho = numel([circS,starS]);
+ndens = numel(rhohat)/nrho;
+rhohat = reshape(rhohat,[nrho, ndens]);
+
 nsub = rcipsav.nsub;
 
+rhohatinterp = cell(nedge,1);
+srcinfo = cell(nedge,1);
+wts = cell(nedge,1);
+
+% figure out which edge the indices belong to 
+% THIS ASSUMES WE HAVE THE SAME ORDER ON EDGES
+circSedge = cell(nedge,1); ncS = numel(circS)/nedge;
+circL1edge = cell(nedge,1); ncL1 = numel(circL1)/nedge;
+starSedge = cell(nedge,1); nsS = numel(starS)/nedge;
+starL1edge = cell(nedge,1); nsL1 = numel(starL1)/nedge;
+for j = 1:nedge
+    circSedge{j} = circS((j-1)*ncS+1:j*ncS);
+    circL1edge{j} = circL1((j-1)*ncL1+1:j*ncL1);
+    starSedge{j} = starS((j-1)*nsS+1:j*nsS);
+    starL1edge{j} = starL1((j-1)*nsL1+1:j*nsL1);
+end
+
 if nargin < 3
     ndepth = nsub;
 end
@@ -84,58 +98,54 @@
     ndepth = nsub;
 end
 
-savedeep = rcipsav.savedeep;
+savedepth = rcipsav.savedepth;
 
-if savedeep
+if ndepth <= savedepth
     
     % all relevant quantities are stored, just run backward recursion
     
     rhohat0 = rhohat;
-    rhohatinterp = [rhohatinterp; rhohat0(circS)];
-    r = [];
-    d = [];
-    d2 = [];
-    n = [];
-    h = [];
     cl = rcipsav.chnkrlocals{nsub};
     wt = weights(cl);
-    r = [r, cl.rstor(:,circL1)];
-    d = [d, cl.dstor(:,circL1)];
-    d2 = [d2, cl.d2stor(:,circL1)];
-    n = [n, cl.nstor(:,circL1)];
-    wts = [wts; wt(circL1(:))];
+
+    for j = 1:nedge
+        rhohatinterp{j} = rhohat0(circSedge{j},:);
+        srcinfo{j}.r = cl.rstor(:,circL1edge{j});
+        srcinfo{j}.d = cl.dstor(:,circL1edge{j});
+        srcinfo{j}.d2 = cl.d2stor(:,circL1edge{j});
+        srcinfo{j}.n = cl.nstor(:,circL1edge{j});
+        wts{j} = wt(circL1edge{j});
+    end
     
     R0 = rcipsav.R{nsub+1};
     for i = 1:ndepth
         R1 = rcipsav.R{nsub-i+1};
         MAT = rcipsav.MAT{nsub-i+1};
         rhotemp = R0\rhohat0;
-        rhohat0 = R1*(Pbc*rhotemp(starS) - MAT*rhohat0(circS));
+        rhohat0 = R1*(Pbc*rhotemp(starS,:) - MAT*rhohat0(circS,:));
         if i == ndepth
-            rhohatinterp = [rhohatinterp; rhohat0];
-            r = [r, cl.rstor(:,starL1)];
-            d = [d, cl.dstor(:,starL1)];
-            d2 = [d2, cl.d2stor(:,starL1)];
-            n = [n, cl.nstor(:,starL1)];
             wt = weights(cl);
-
-            wts = [wts; wt(starL1(:))];
+            for j = 1:nedge
+                rhohatinterp{j} = [rhohatinterp{j}; rhohat0([starSedge{j} circSedge{j}],:)];
+                srcinfo{j}.r = [srcinfo{j}.r, cl.rstor(:,starL1edge{j})];
+                srcinfo{j}.d = [srcinfo{j}.d, cl.dstor(:,starL1edge{j})];
+                srcinfo{j}.d2 = [srcinfo{j}.d2, cl.d2stor(:,starL1edge{j})];
+                srcinfo{j}.n = [srcinfo{j}.n, cl.nstor(:,starL1edge{j})];                
+                wts{j} = [wts{j}, wt(starL1edge{j})];
+            end
         else
             cl = rcipsav.chnkrlocals{nsub-i};
-            rhohatinterp = [rhohatinterp; rhohat0(circS)];
-            r = [r, cl.rstor(:,circL1)];
-            d = [d, cl.dstor(:,circL1)];
-            d2 = [d2, cl.d2stor(:,circL1)];
-            n = [n, cl.nstor(:,circL1)];
             wt = weights(cl);
-            wts = [wts; wt(circL1(:))];
+            for j = 1:nedge
+                rhohatinterp{j} = [rhohatinterp{j}; rhohat0(circSedge{j},:)];
+                srcinfo{j}.r = [srcinfo{j}.r, cl.rstor(:,circL1edge{j})];
+                srcinfo{j}.d = [srcinfo{j}.d, cl.dstor(:,circL1edge{j})];
+                srcinfo{j}.d2 = [srcinfo{j}.d2, cl.d2stor(:,circL1edge{j})];
+                srcinfo{j}.n = [srcinfo{j}.n, cl.nstor(:,circL1edge{j})];
+                wts{j} = [wts{j}, wt(circL1edge{j})];
+            end
         end
         R0 = R1;
     end
     
-    srcinfo.r = r;
-    srcinfo.d = d;
-    srcinfo.d2 = d2;
-    srcinfo.n = n;
-
 end
diff --git a/chunkie/chunkermat.m b/chunkie/chunkermat.m
index 362d02f6..6cd5ee17 100644
--- a/chunkie/chunkermat.m
+++ b/chunkie/chunkermat.m
@@ -65,6 +65,7 @@
 %                      corrections for near corners if input chnkobj is
 %                      of type chunkergraph
 %           opts.rcip_ignore = [], list of vertices to ignore in rcip
+%           opts.rcip_savedepth = (10), depth to save rcip info
 %           opts.nsub_or_tol = (40) specify the level of refinements in rcip
 %                    or a tolerance where the number of levels is given by
 %                    ceiling(log_{2}(1/tol^2));
@@ -135,6 +136,7 @@
 isrcip = true;
 rcip_ignore = [];
 nsub = 40;
+rcip_savedepth = 10;
 adaptive_correction = false;
 sing = 'log';
 
@@ -162,6 +164,9 @@
 if(isfield(opts,'rcip'))
     isrcip = opts.rcip;
 end
+if(isfield(opts,'rcip_savedepth'))
+    rcip_savedepth = opts.rcip_savedepth;
+end
 if(isfield(opts,'rcip_ignore'))
     rcip_ignore = opts.rcip_ignore;
     if ~isrcip && isempty(rcip_ignore)
@@ -466,6 +471,7 @@
     npt_all = horzcat(chnkobj.echnks.npt);
     [~,nv] = size(chnkobj.verts);
     ngl = chnkrs(1).k;
+
     rcipsav = cell(nv,1);
     
     for ivert=setdiff(1:nv,rcip_ignore)
@@ -513,12 +519,14 @@
         optsrcip = opts;
         optsrcip.nonsmoothonly = false;
         optsrcip.corrections = false;
+        optsrcip.rcip_savedepth = rcip_savedepth;
 
         [R,rcipsav{ivert}] = chnk.rcip.Rcompchunk(chnkrs,iedgechunks,kern,ndim, ...
             Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,starL1,circL1,... 
             sbclmat,sbcrmat,lvmat,rvmat,u,optsrcip);
+
         rcipsav{ivert}.starind = starind;
-        
+
         sysmat_tmp = inv(R) - eye(2*ngl*nedge*ndim);
         if (~nonsmoothonly)
             
@@ -580,6 +588,7 @@
     if nargout > 2
         varargout{2} = rcipsav;
     end
+
     
 end
 

From df3065e2164d1c7811dd0d791b788769e1105274 Mon Sep 17 00:00:00 2001
From: Travis Askham <askhamwhat@gmail.com>
Date: Fri, 27 Jun 2025 21:22:08 -0400
Subject: [PATCH 3/8] fix behavior of last panel depending on orientation

---
 chunkie/+chnk/+rcip/rhohatInterp.m | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/chunkie/+chnk/+rcip/rhohatInterp.m b/chunkie/+chnk/+rcip/rhohatInterp.m
index b10a58b2..1cabc642 100644
--- a/chunkie/+chnk/+rcip/rhohatInterp.m
+++ b/chunkie/+chnk/+rcip/rhohatInterp.m
@@ -74,6 +74,8 @@
 srcinfo = cell(nedge,1);
 wts = cell(nedge,1);
 
+ileftright = rcipsav.ileftright
+
 % figure out which edge the indices belong to 
 % THIS ASSUMES WE HAVE THE SAME ORDER ON EDGES
 circSedge = cell(nedge,1); ncS = numel(circS)/nedge;
@@ -126,7 +128,11 @@
         if i == ndepth
             wt = weights(cl);
             for j = 1:nedge
-                rhohatinterp{j} = [rhohatinterp{j}; rhohat0([starSedge{j} circSedge{j}],:)];
+                if ileftright(j) == 1
+                    rhohatinterp{j} = [rhohatinterp{j}; rhohat0([circSedge{j} starSedge{j}],:)];
+                else
+                    rhohatinterp{j} = [rhohatinterp{j}; rhohat0([starSedge{j} circSedge{j}],:)];
+                end
                 srcinfo{j}.r = [srcinfo{j}.r, cl.rstor(:,starL1edge{j})];
                 srcinfo{j}.d = [srcinfo{j}.d, cl.dstor(:,starL1edge{j})];
                 srcinfo{j}.d2 = [srcinfo{j}.d2, cl.d2stor(:,starL1edge{j})];

From 83fcb25194b595a75ec265e4aaa9b8459830b579 Mon Sep 17 00:00:00 2001
From: Travis Askham <askhamwhat@gmail.com>
Date: Tue, 7 Oct 2025 09:37:40 -0400
Subject: [PATCH 4/8] bug fix

---
 chunkie/chunkerkernevalmat.m | 127 ++++++++++++++++++++++++++++++++---
 1 file changed, 119 insertions(+), 8 deletions(-)

diff --git a/chunkie/chunkerkernevalmat.m b/chunkie/chunkerkernevalmat.m
index 6df13c6b..20d0a408 100644
--- a/chunkie/chunkerkernevalmat.m
+++ b/chunkie/chunkerkernevalmat.m
@@ -283,13 +283,13 @@
 
 % smooth for sufficiently far, adaptive otherwise
 
-rho = 1.8;
-optsflag = [];  optsflag.rho = rho;
-flag = flagnear_rectangle(chnkr,targinfo.r,optsflag);
+%rho = 1.8;
+%optsflag = [];  optsflag.rho = rho;
+%flag = flagnear_rectangle(chnkr0,targinfo0.r,optsflag);
 
-npoly = chnkr.k*2;
-nlegnew = chnk.ellipse_oversample(rho,npoly,eps);
-nlegnew = max(nlegnew,chnkr.k);
+%npoly = chnkr.k*2;
+%nlegnew = chnk.ellipse_oversample(rho,npoly,eps);
+%nlegnew = max(nlegnew,chnkr.k);
 
 
 spmat = chunkerkernevalmat_adap(chnkr0,kern0,opdims0, ...
@@ -299,8 +299,8 @@
 if corrections
     % TODO: find more elegant solution that avoids building a dense flag matrix
     flaginv = ~flag;
-    mat0 = spmat - chunkerkernevalmat_smooth(chnkr0,kern0,opdims0,targinfo0, ...
-        flaginv,opts);
+    mat0 = spmat - chunkerkernevalmat_smooth2(chnkr0,kern0,opdims0,targinfo0, ...
+        flag,opts);
     mat(irow0,icol0) = sparse(mat0);
     continue
 end
@@ -367,6 +367,117 @@
 
 end
 
+function mat = chunkerkernevalmat_smooth2(chnkr,kern,opdims, ...
+    targinfo,flag,opts)
+
+if isa(kern,'kernel')
+    kerneval = kern.eval;
+else
+    kerneval = kern;
+end
+
+k = chnkr.k;
+nch = chnkr.nch;
+
+if nargin < 5
+    flag = [];
+end
+if nargin < 6
+    opts = [];
+end
+
+% Extract target info
+targs = targinfo.r;
+[~,nt] = size(targs);
+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
+datat = [];
+if isfield(targinfo, 'data')
+    datat = targinfo.data;
+end
+
+
+% using adaptive quadrature
+
+
+if isempty(flag)
+
+
+    mat = sparse(opdims(1)*nt,opdims(2)*chnkr.npt);
+    
+else
+    is = zeros(nnz(flag)*opdims(1)*opdims(2)*k,1);
+    js = is;
+    vs = is;
+    istart = 1;
+
+    r = chnkr.r;
+    d = chnkr.d;
+    n = chnkr.n;
+    d2 = chnkr.d2;
+    data = chnkr.data;
+    for i = 1:nch
+        jmat = 1 + (i-1)*k*opdims(2);
+        jmatend = i*k*opdims(2);
+                        
+        [ji] = find(flag(:,i));
+        datat2 = [];
+        if ~isempty(datat)
+            datat2 = datat(:,ji);
+        end
+        srcinfo = [];
+        srcinfo.r = r(:,:,i);
+        srcinfo.d = d(:,:,i);
+        srcinfo.n = n(:,:,i);
+        srcinfo.d2 = d2(:,:,i);
+
+        if ~isempty(data); srcinfo.data = data(:,:,i); end
+
+        targinfo = [];
+        targinfo.r = targs(:,ji);
+        targinfo.d = targd(:,ji);
+        targinfo.n = targn(:,ji);
+        targinfo.d2 = targd2(:,ji);
+        targinfo.data = datat2;
+
+        mat1 = kerneval(srcinfo,targinfo);
+        wts1 = chnkr.wts(:,i);
+        wts2 = repmat( (wts1(:)).', opdims(2), 1);
+        wts2 = ( wts2(:) ).';
+        mat1 = mat1.*wts2;
+
+        js1 = jmat:jmatend;
+        js1 = repmat( (js1(:)).',opdims(1)*numel(ji),1);
+                
+        indji = (ji-1)*opdims(1);
+        indji = repmat( (indji(:)).', opdims(1),1) + ( (1:opdims(1)).');
+        indji = indji(:);
+        
+        indji = repmat(indji,1,opdims(2)*k);
+        
+        iend = istart+numel(mat1)-1;
+        is(istart:iend) = indji(:);
+        js(istart:iend) = js1(:);
+        vs(istart:iend) = mat1(:);
+        istart = iend+1;
+    end
+    mat = sparse(is,js,vs,opdims(1)*nt,opdims(2)*chnkr.npt);
+    
+end
+
+end
+
 function mat = chunkerkernevalmat_adap(chnkr,kern,opdims, ...
     targinfo,flag,opts)
 

From 36e940806f2000ebe024e03b1b02e3a0305c91b7 Mon Sep 17 00:00:00 2001
From: Travis Askham <askhamwhat@gmail.com>
Date: Mon, 27 Oct 2025 11:50:37 -0400
Subject: [PATCH 5/8] revert

---
 chunkie/chunkerkernevalmat.m | 1 -
 1 file changed, 1 deletion(-)

diff --git a/chunkie/chunkerkernevalmat.m b/chunkie/chunkerkernevalmat.m
index 47362714..b8d5c8ce 100644
--- a/chunkie/chunkerkernevalmat.m
+++ b/chunkie/chunkerkernevalmat.m
@@ -289,7 +289,6 @@
 
 if corrections
     mat0 = spmat - chunkerkernevalmat_smooth_corrections(chnkr0,kern0,opdims0,targinfo0, ...
-
         flag,opts);
     mat(irow0,icol0) = sparse(mat0);
     continue

From 5306e7be32174f2a17a9eb7ab9f4aef1fcccf443 Mon Sep 17 00:00:00 2001
From: Travis Askham <askhamwhat@gmail.com>
Date: Mon, 27 Oct 2025 14:59:43 -0400
Subject: [PATCH 6/8] suppress a print and update test

---
 chunkie/+chnk/+rcip/rhohatInterp.m |  2 +-
 devtools/test/rcipTest.m           | 66 ++++++++++++++++++++++++------
 2 files changed, 54 insertions(+), 14 deletions(-)

diff --git a/chunkie/+chnk/+rcip/rhohatInterp.m b/chunkie/+chnk/+rcip/rhohatInterp.m
index 1cabc642..6c544fb4 100644
--- a/chunkie/+chnk/+rcip/rhohatInterp.m
+++ b/chunkie/+chnk/+rcip/rhohatInterp.m
@@ -74,7 +74,7 @@
 srcinfo = cell(nedge,1);
 wts = cell(nedge,1);
 
-ileftright = rcipsav.ileftright
+ileftright = rcipsav.ileftright;
 
 % figure out which edge the indices belong to 
 % THIS ASSUMES WE HAVE THE SAME ORDER ON EDGES
diff --git a/devtools/test/rcipTest.m b/devtools/test/rcipTest.m
index 4592c21d..d1f484b1 100644
--- a/devtools/test/rcipTest.m
+++ b/devtools/test/rcipTest.m
@@ -1,8 +1,6 @@
 rcipTest0();
 
-
 function rcipTest0()
-
 %RCIPTEST
 %
 % This file tests the rcip routines for solving the exterior dirichlet 
@@ -54,12 +52,22 @@ function rcipTest0()
 fcurve = cell(1,ncurve);
 figure(1)
 clf; hold on;
+chnkrcell = cell(ncurve,1);
 for icurve = 1:ncurve
     fcurve{icurve} = @(t) circulararc(t,cpars{icurve});
     chnkr(icurve) = chunkerfuncuni(fcurve{icurve},nch(icurve),cparams{icurve},pref);
     plot(chnkr(icurve)); quiver(chnkr(icurve));
+    chnkrcell{icurve} = chnkr(icurve);
 end
 
+cends = chunkends(chnkrcell{icurve});
+vert = [cends(:,1,1), cends(:,2,end)];
+
+cg = chunkgraph(vert,[2 1; 1 2],chnkrcell); % equivalent chunkgraph
+
+%%
+
+
 iedgechunks = [1 2; 1 chnkr(2).nch];
 isstart = [1 0];
 nedge = 2;
@@ -85,7 +93,7 @@ function rcipTest0()
 ts = 0.0+2*pi*rand(1,nt);
 targets = [cos(ts);sin(ts)];
 targets = 0.2*targets;
-targets(:,1) = [0.95;0];
+targets(:,1) = [0.999999;0];
 targets(:,2) = [0,0.36];
 targets(:,3) = [-0.95;0];
 targets(:,2) = [0,0.36];
@@ -154,8 +162,6 @@ function rcipTest0()
 
 nsub = 40;
 
-opts = [];
-
 
 for icorner=1:ncorner
     clist = corners{icorner}.clist;
@@ -180,11 +186,11 @@ function rcipTest0()
     
     [Pbc,PWbc,starL,circL,starS,circS,ilist,starL1,circL1] = chnk.rcip.setup(ngl,ndim, ...
       nedge,isstart);
-    opts_use = [];
     
     iedgechunks = corners{icorner}.iedgechunks;
     optsrcip = [];
-    optsrcip.savedeep = true;
+    optsrcip.nsub = nsub;
+    optsrcip.rcip_savedepth = Inf;
     tic; [R{icorner},rcipsav{icorner}] = chnk.rcip.Rcompchunk(chnkr,iedgechunks,fkern,ndim, ...
         Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,starL1,circL1,...,
         [],[],[],[],[],optsrcip);
@@ -194,18 +200,48 @@ function rcipTest0()
 end
 sysmat = sysmat + eye(np);
 
+opts = [];
+opts.nsub = nsub; 
+opts.rcip_savedepth = Inf;
+[sysmat2,~,rcipsav] = chunkermat(cg,fkern,opts);
+
+sysmat2 = sysmat2 + eye(np);
+
 [sol,flag,relres,iter] = gmres(sysmat,ubdry,np,eps*20,np);
+[sol2,flag,relres,iter] = gmres(sysmat2,ubdry,np,eps*20,np);
 
 %
 
 % interpolate to fine grid
 
-ndepth = 10;
+ndepth = 20;
 cor = cell(1,ncorner);
 
 for icorner = 1:ncorner
-    solhat = sol(starinds{icorner});
-    [solhatinterp,srcinfo,wts] = chnk.rcip.rhohatInterp(solhat,rcipsav{icorner},ndepth);
+    starindtmp = rcipsav{icorner}.starind;
+    solhat = sol(starindtmp);
+    [solhatinterpcell,srcinfocell,wtscell] = chnk.rcip.rhohatInterp(solhat,rcipsav{icorner},ndepth);
+
+    rtmp = [];
+    dtmp = [];
+    d2tmp = [];
+    ntmp = [];
+    datatmp = [];
+    solhatinterp = [];
+    wts = [];
+    for j = 1:length(srcinfocell)
+        rtmp = [rtmp, srcinfocell{j}.r];
+        dtmp = [dtmp, srcinfocell{j}.d];
+        d2tmp = [d2tmp, srcinfocell{j}.d2];
+        ntmp = [ntmp, srcinfocell{j}.n];
+        if isfield(srcinfocell{j},'data')
+            datatmp = [datatmp, srcinfocell{j}.data];
+        end
+        solhatinterp = [solhatinterp; solhatinterpcell{j}(:)];
+        wts = [wts; wtscell{j}(:)];
+
+    end
+    srcinfo = struct('r',rtmp,'d',dtmp,'d2',d2tmp,'n',ntmp,'data',datatmp);
 
     targtemp = targets(:,:) - rcipsav{icorner}.ctr(:,1);
     targinfo = [];
@@ -242,6 +278,13 @@ function rcipTest0()
 fprintf('relative frobenius error %5.2e\n',relerr);
 fprintf('relative l_inf/l_1 error %5.2e\n',relerr2);
 
+assert(relerr < 1e-10)
+assert(relerr2 < 1e-10)
+
+
+
+end
+
 
 %%
 
@@ -255,9 +298,6 @@ function rcipTest0()
 
 
 
-end
-
-
 function [r,d,d2] = circulararc(t,cpars)
 %%circulararc
 % return position, first and second derivatives of a circular arc

From 068126a37ca6527516bb68879c26e3b7110fb076 Mon Sep 17 00:00:00 2001
From: Manas Rachh <mrachh@flatironinstitute.org>
Date: Tue, 28 Oct 2025 12:13:31 +0530
Subject: [PATCH 7/8] updating docs for chunkermat for rcipsav

---
 chunkie/chunkermat.m | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/chunkie/chunkermat.m b/chunkie/chunkermat.m
index 6cd5ee17..15b90cd0 100644
--- a/chunkie/chunkermat.m
+++ b/chunkie/chunkermat.m
@@ -84,6 +84,9 @@
 % Optional output
 %   opts - with the updated opts structure which stores the relevant
 %          quantities in opts.auxquads.<opts.quad><opts.type>
+%   rcipsav - precomputed structure of rcip data at corners
+%             for subsequent postprocessing solution at targets closed 
+%             to the corner
 %
 % Examples:
 %   sysmat = chunkermat(chnkr,kern); % standard options

From 65b298c911193a38a5a149add4156f73b203da40 Mon Sep 17 00:00:00 2001
From: Travis Askham <askhamwhat@gmail.com>
Date: Tue, 28 Oct 2025 16:05:28 -0400
Subject: [PATCH 8/8] Update chunkermat.m

---
 chunkie/chunkermat.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/chunkie/chunkermat.m b/chunkie/chunkermat.m
index 15b90cd0..1c48c6f3 100644
--- a/chunkie/chunkermat.m
+++ b/chunkie/chunkermat.m
@@ -85,7 +85,7 @@
 %   opts - with the updated opts structure which stores the relevant
 %          quantities in opts.auxquads.<opts.quad><opts.type>
 %   rcipsav - precomputed structure of rcip data at corners
-%             for subsequent postprocessing solution at targets closed 
+%             for subsequent postprocessing of the solution at targets close 
 %             to the corner
 %
 % Examples: