Merge pull request #147 from fastalgorithms/for-saw

For saw

Author: askhamwhat

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);

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,42 @@
 % 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);
+
+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;
+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 +100,58 @@
     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
+                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})];
+                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

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));
@@ -83,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 of the solution at targets close 
+%             to the corner
 %
 % Examples:
 %   sysmat = chunkermat(chnkr,kern); % standard options
@@ -135,6 +139,7 @@
 isrcip = true;
 rcip_ignore = [];
 nsub = 40;
+rcip_savedepth = 10;
 adaptive_correction = false;
 sing = 'log';
 
@@ -162,6 +167,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 +474,8 @@
     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)
         clist = chnkobj.vstruc{ivert}{1};
@@ -512,11 +522,14 @@
         optsrcip = opts;
         optsrcip.nonsmoothonly = false;
         optsrcip.corrections = false;
+        optsrcip.rcip_savedepth = rcip_savedepth;
 
-        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)
 
@@ -574,6 +587,11 @@
             vsysmat = [vsysmat;sysmat_tmp(:)];
         end    
     end
+
+    if nargout > 2
+        varargout{2} = rcipsav;
+    end
+
 
 end