changes behavior of chunkermatapply and building corrections matrix

- removes dval from chunkermatapply interface for consistency with chunkermat
- adds a warning about fmms
- tries to speed up kerneval_smooth in the most common situation
- the previous version of corrections was slow, this avoids some
slow indexing that was happening.
- this adds a faster lege.barywts option when nodes are known

chunkie/+chnk/+quadggq/buildmattd.m

@@ -1,4 +1,4 @@
-function [spmat] = buildmattd(chnkr,kern,opdims,type,auxquads,ilist)
+function [spmat] = buildmattd(chnkr,kern,opdims,type,auxquads,ilist,corrections)
 %CHNK.QUADGGQ.BUILDMAT build matrix for given kernel and chnkr 
 % description of boundary, using special quadrature for self
 % and neighbor panels.
@@ -15,6 +15,9 @@
 if (nargin <6)
     ilist = [];
 end
+if nargin < 7
+    corrections = false;
+end
 
 k = chnkr.k;
 nch = chnkr.nch;
@@ -54,6 +57,17 @@
 
 % do smooth weight for all
 %sysmat = chnk.quadnative.buildmat(chnkr,kern,opdims,1:nch,1:nch,wts);
+
+if corrections
+    wtss = chnkr.wts;
+    wtss = repmat(wtss(:).',opdims(2),1); wtss = reshape(wtss,opdims(2)*k,nch);
+    indd = kron(eye(k),true(opdims(1),opdims(2)));
+    indd = indd(:) > 0;
+else
+    wtss = [];
+    indd = [];
+end
+
 mmat = k*nch*opdims(1); nmat = k*nch*opdims(2);
 
 nnz = k*nch*opdims(1)*k*3*opdims(2);
@@ -89,7 +103,7 @@
         % skip construction if both chunks are in the "bad" chunk list
         else
             submat = chnk.quadggq.nearbuildmat(r,d,n,d2,data,ibefore,j, ...
-                kern,opdims,xs1,wts1,ainterp1kron,ainterp1);
+                kern,opdims,xs1,wts1,ainterp1kron,ainterp1,corrections,wtss);
 
             imat = 1 + (ibefore-1)*k*opdims(1);
             induse = ict:ict+nnz1-1;
@@ -106,7 +120,7 @@
         % skip construction if both chunks are in the "bad" chunk list
       else      
         submat = chnk.quadggq.nearbuildmat(r,d,n,d2,data,iafter,j, ...
-            kern,opdims,xs1,wts1,ainterp1kron,ainterp1);
+            kern,opdims,xs1,wts1,ainterp1kron,ainterp1,corrections,wtss);
 
 
         imat = 1 + (iafter-1)*k*opdims(1);
@@ -124,7 +138,7 @@
       % skip construction if the chunk is in the "bad" chunk list  
     else
       submat = chnk.quadggq.diagbuildmat(r,d,n,d2,data,j,kern,opdims,...
-          xs0,wts0,ainterps0kron,ainterps0);
+          xs0,wts0,ainterps0kron,ainterps0,corrections,wtss,indd);
 
       imat = 1 + (j-1)*k*opdims(1);
 

chunkie/+chnk/+quadggq/diagbuildmat.m

@@ -1,5 +1,5 @@
 function submat = diagbuildmat(r,d,n,d2,data,i,fkern,opdims,...
-			 xs0,whts0,ainterps0kron,ainterps0)
+			 xs0,whts0,ainterps0kron,ainterps0,corrections,wtss,indd)
 %CHNK.QUADGGQ.DIAGBUILDMAT                  
 
 % grab specific boundary data
@@ -11,6 +11,14 @@
 else
     dd = data(:,:,i);
 end
+
+if nargin < 13
+    corrections = false;
+    wtss = [];
+    indd = [];
+end
+
+
 % interpolate boundary info
 
 % get relevant coefficients
@@ -23,13 +31,13 @@
 d2fine = cell(k,1);
 dsdt = cell(k,1);
 
-for i = 1:k
-    rfine{i} = (ainterps0{i}*(rs.')).';
-    dfine{i} = (ainterps0{i}*(ds.')).';
-    d2fine{i} = (ainterps0{i}*(d2s.')).';
-    dfinenrm = sqrt(sum(dfine{i}.^2,1));
-    nfine{i} = [dfine{i}(2,:); -dfine{i}(1,:)]./dfinenrm;
-    dsdt{i} = (dfinenrm(:)).*whts0{i};
+for j = 1:k
+    rfine{j} = (ainterps0{j}*(rs.')).';
+    dfine{j} = (ainterps0{j}*(ds.')).';
+    d2fine{j} = (ainterps0{j}*(d2s.')).';
+    dfinenrm = sqrt(sum(dfine{j}.^2,1));
+    nfine{j} = [dfine{j}(2,:); -dfine{j}(1,:)]./dfinenrm;
+    dsdt{j} = (dfinenrm(:)).*whts0{j};
 end
 
 srcinfo = [];
@@ -54,6 +62,21 @@
     smatbigi*ainterps0kron{j};
 end
 
+if corrections
+    srcinfo.r = rs;
+    srcinfo.d = ds;
+    srcinfo.d2 = d2s;
+    srcinfo.n = ns;
+
+    targinfo.r = rs; targinfo.d = ds; targinfo.d2 = d2s; targinfo.n = ns;
+    targinfo.data = dd;
+
+    sc = fkern(srcinfo,targinfo);
+    sc(indd) = 0;
+
+    wtsi = wtss(:,i);
+    submat = submat - sc.*(wtsi(:).');
+
 end
 
 

chunkie/+chnk/+quadggq/nearbuildmat.m

@@ -1,5 +1,5 @@
 function submat = nearbuildmat(r,d,n,d2,data,i,j,fkern,opdims,...
-				      xs1,whts1,ainterp1kron,ainterp1)
+				      xs1,whts1,ainterp1kron,ainterp1,corrections,wtss)
 %CHNKR.QUADGGQ.NEARBUILDMAT
 
 % grab specific boundary data
@@ -12,6 +12,12 @@
 else
     dd = data(:,:,i);
 end
+
+if nargin < 14
+    corrections = false;
+    wtss = [];
+end
+
 % interpolate boundary info
 
 % get relevant coefficients
@@ -73,5 +79,12 @@
 smatbig = fkern(srcinfo,targinfo);
 submat = smatbig*diag(dsdtndim2)*ainterp1kron;
 
+if corrections
+    srcinfo = []; srcinfo.r = rs; srcinfo.d = ds; 
+    srcinfo.d2 = d2s; srcinfo.n = ns; 
+
+    wtsj = wtss(:,j);
+    submat = submat - fkern(srcinfo,targinfo).*(wtsj(:).');
+
  end
 

chunkie/+lege/barywts.m

@@ -1,4 +1,4 @@
-function w = barywts(k)
+function w = barywts(k,x)
 %LEGE.BARYWTS
 %
 % returns the weights for barycentric Lagrange interpolation
@@ -16,7 +16,9 @@
 %
 % where x(j) is the jth Legendre node of order k
 
-x = lege.exps(k);
+if nargin < 2
+    x = lege.exps(k);
+end
 
 xx = bsxfun(@minus,x,x.');
 xx(1:k+1:end) = 1;