Merge branch 'master' of https://github.com/fastalgorithms/chunkie

chunkie/+chnk/+flam/kernbyindexr.m

@@ -1,4 +1,4 @@
-function mat = kernbyindexr(i,j,targs,chnkr,kern,opdims,spmat)
+function mat = kernbyindexr(i,j,targobj,chnkr,kern,opdims,spmat)
 %% evaluate system matrix by entry index utility function for 
 % general kernels, with replacement for specific entries and the
 % ability to add a low-rank modification, rectangular version
@@ -17,6 +17,10 @@
 %
 % i - array of row indices to compute
 % j - array of col indices to compute
+% targobj - object describing the target points, can be specified as
+%       * array of points
+%       * chunker object
+%       * chunkgraph object
 % chnkr - chunker object describing boundary
 % whts - smooth integration weights on chnkr
 % kern - kernel function of the form kern(s,t,stau,ttau) where s and t 
@@ -42,15 +46,32 @@
 [juni,~,ijuni] = unique(jpts);
 
 % matrix-valued entries of kernel for unique points
-
-ri = targs(:,iuni); rj = chnkr.r(:,juni);
+rj = chnkr.r(:,juni);
 dj = chnkr.d(:,juni); d2j = chnkr.d2(:,juni);
 nj = chnkr.n(:,juni);
-%di = bsxfun(@rdivide,di,sqrt(sum(di.^2,1)));
-%dj = bsxfun(@rdivide,dj,sqrt(sum(dj.^2,1)));
 srcinfo = []; srcinfo.r = rj; srcinfo.d = dj; srcinfo.d2 = d2j;
 srcinfo.n = nj;
-targinfo = []; targinfo.r = ri;
+% Assign appropriate object to targinfo
+targinfo = [];
+if isa(targobj, "chunker") || isa(targobj, "chunkgraph")
+    targinfo.r = targobj.r(:,iuni);
+    targinfo.d = targobj.d(:,iuni);
+    targinfo.d2 = targobj.d2(:,iuni);
+    targinfo.n = targobj.n(:,iuni);
+elseif isstruct(targobj)
+    if isfield(targobj,'d')
+        targinfo.d = targobj.d(:,iuni);
+    end
+    if isfield(targobj,'d2')
+        targinfo.d = targobj.d(:,iuni);
+    end
+    if isfield(targobj,'n')
+        targinfo.n = targobj.n(:,iuni);
+    end
+    targinfo.r = targobj.r(:,iuni);
+else
+    targinfo.r = targobj(:,iuni);
+end
 
 matuni = kern(srcinfo,targinfo);
 

chunkie/+chnk/+flam/proxy_square_pts.m

@@ -1,4 +1,4 @@
-function [pr,ptau,pw,pin] = proxy_square_pts(porder)
+function [pr,ptau,pw,pin] = proxy_square_pts(porder, opts)
 %PROXY_SQUARE_PTS return the proxy points, unit tangents, weights, and 
 % function handle for determining if points are within proxy surface.
 % This function is for the proxy surface around a unit box centered at 
@@ -16,12 +16,31 @@
     porder = 64;
 end
 
+if nargin < 2
+    opts = [];
+end
+
+iflege = 0;
+if isfield(opts, 'iflege')
+    iflege = opts.iflege;
+end
+
+
 assert(mod(porder,4) == 0, ...
     'number of proxy points on square should be multiple of 4')
 
 po4 = porder/4;
 
-pts = -1.5+3*(0:(po4-1))*1.0/po4;
+if ~iflege
+    pts = -1.5+3*(0:(po4-1))*1.0/po4;
+    pw = ones(porder,1)*(4.0*3.0/porder);
+else
+    [xleg, wleg] = lege.exps(po4);
+    pts = xleg.'*1.5;
+    wleg = wleg*1.5;   
+    pw = [wleg; wleg; wleg; wleg];
+end
+
 one = ones(size(pts));
 
 pr = [pts, one*1.5, -pts, -1.5*one;
@@ -30,8 +49,6 @@
 ptau = [one, one*0, -one, one*0;
     one*0, one, one*0, -one];
 
-pw = ones(porder,1)*(4.0*3.0/porder);
-
 pin = @(x) max(abs(x),[],1) < 1.5;
 
 end

chunkie/+chnk/+helm2d/fmm.m

@@ -38,6 +38,8 @@
 %                   note that targinfo must contain targ.n
 %                type = 'dprime' normal derivative of double layer,
 %                   note that targinfo must contain targ.n
+%                type = 'cprime' normal derivative of combined layer,
+%                   note that targinfo must contain targ.n
 %   sigma - density
 %   varargin{1} - coef in the combined layer formula, otherwise
 %                does nothing
@@ -61,7 +63,7 @@
     case {'d', 'dprime'}
         srcuse.dipstr = sigma(:).';
         srcuse.dipvec = srcinfo.n(1:2,:);
-    case 'c'
+    case {'c', 'cprime'}
         coefs = varargin{1};
         srcuse.charges = coefs(2)*sigma(:).';
         srcuse.dipstr  = coefs(1)*sigma(:).';
@@ -83,7 +85,7 @@
     switch lower(type)
         case {'s', 'd', 'c'}
             varargout{1} = U.pottarg.';
-        case {'sprime', 'dprime'}
+        case {'sprime', 'dprime', 'cprime'}
             if ( ~isfield(targinfo, 'n') )
                 error('CHUNKIE:helm2d:fmm:normals', ...
                     'Targets require normal info when evaluating Helmholtz kernel ''%s''.', type);

chunkie/+chnk/+helm2d/kern.m

@@ -133,14 +133,14 @@
 % normal derivative of combined field
 if strcmpi(type,'cprime')
   coef = ones(2,1);
-  if(nargin == 5); coef = varargin{1}; end;
+  if(nargin == 5); coef = varargin{1}; end
   targnorm = targinfo.n;
   srcnorm = srcinfo.n;
 
-  submat = zeros(nt,ns);
+
 
   % Get gradient and hessian info
-  [submats,grad,hess] = chnk.helm2d.green(zk,src,targ);
+  [~,grad,hess] = chnk.helm2d.green(zk,src,targ);
 
   nxsrc = repmat(srcnorm(1,:),nt,1);
   nysrc = repmat(srcnorm(2,:),nt,1);
@@ -161,7 +161,7 @@
 % Dirichlet and neumann data corresponding to combined field
 if strcmpi(type,'c2trans') 
   coef = ones(2,1);
-  if(nargin == 5); coef = varargin{1}; end;
+  if(nargin == 5); coef = varargin{1}; end
   targnorm = targinfo.n;
   srcnorm = srcinfo.n;
 

chunkie/+chnk/+quadnative/buildmat.m

@@ -36,12 +36,7 @@
 targinfo = []; targinfo.r = rt; targinfo.d = dt; targinfo.d2 = d2t; targinfo.n = nt;
 hs = h(j); ht = h(i);
 
-dsnrms = sqrt(sum((ds).^2,1));
-%taus = bsxfun(@rdivide,ds,dsnrms);
-
-%dtnrms = sqrt(sum(abs(dt).^2,1));
-%taut = bsxfun(@rdivide,dt,dtnrms);
-
+dsnrms = sqrt(sum(ds.^2,1));
 ws = kron(hs(:),wts(:));
 
 dsdt = dsnrms(:).*ws;

chunkie/+chnk/+rcip/IPinit.m

@@ -1,10 +1,16 @@
   function [IP,IPW]=IPinit(T,W)
+  %CHNK.RCIP.IPinit 
+  %
   % construct the prolongation matrix IP that maps function values
   % on n_{gl} Gauss-Legendre nodes on [-1,1] to function values at 
   % 2n_{gl} Gauss-Legendre, with shifted and scaled n_{gl}
   % Gauss-Legendre nodes on each subinterval [-1,0], [0,1], respectively.
   %
   % IPW is the weighted prolongation matrix acted on the left side. 
+  %
+
+  % author: Johan Helsing (part of RCIP tutorial)
+
   ngl = length(T);
   A=ones(ngl);
   AA=ones(2*ngl,ngl);

chunkie/+chnk/+rcip/Rcompchunk.m

@@ -1,46 +1,73 @@
-function [R]=Rcompchunk(chnkr,iedgechunks,fkern,ndim, ...
-    Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,...
-    glxs,sbcmat,lvmat,u,opts)
+function [R,rcipsav]=Rcompchunk(chnkr,iedgechunks,fkern,ndim, ...
+    Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,starL1,circL1,...
+    sbclmat,sbcrmat,lvmat,rvmat,u,opts)
 %CHNK.RCIP.Rcompchunk carry out the forward recursion for computing
 % the preconditioner R where geometry is described as a chunker
 %
 % This routine is not intended to be user-callable 
 %
-% Adapted from Shidong Jiang's RCIP implementation
-%
 % Function is passed as a handle, number of equations is given by
 % ndim
 %
 % Kernel on input takes in arguments (chnkrlocal,ilistl);
 % 
 % Note that matrix must be scaled to have identity on the diagonal,
 % will not work with scaled version of identity
-
+%
+
+% author: Shidong Jiang
+% modified: Jeremy Hoskins, Manas Rachh
+
+
 k = chnkr.k;  
 dim = chnkr.dim;
+nedge = size(iedgechunks,2);
 
+glxs = chnkr.tstor;
 glws = chnkr.wstor;
 
-if nargin < 14
+% return what's needed to interpolate from coarse
+
+rcipsav = [];
+rcipsav.k = k;
+rcipsav.ndim = ndim;
+rcipsav.nedge = nedge;
+rcipsav.Pbc = Pbc;
+rcipsav.PWbc = PWbc;
+rcipsav.starL = starL;
+rcipsav.starL1 = starL1;
+rcipsav.starS = starS;
+rcipsav.circL = circL;
+rcipsav.circL1 = circL1;
+rcipsav.circS = circS;
+rcipsav.ilist = ilist;
+rcipsav.nsub = nsub;
+
+if (nargin < 15 || isempty(sbclmat) || isempty(sbcrmat) || ...
+        isempty(lvmat) || isempty(rvmat) || isempty(u))
     [sbclmat,sbcrmat,lvmat,rvmat,u] = chnk.rcip.shiftedlegbasismats(k); 
 end
 
-if nargin < 17
+if nargin < 20
     opts = [];
 end
 
-nedge = size(iedgechunks,2);
+savedeep = false;
+if isfield(opts,'savedeep')
+    savedeep = opts.savedeep;
+end
 
-ileftright = zeros(nedge,1);
-nextchunk = zeros(nedge,1);
+rcipsav.savedeep = savedeep;
+
+if savedeep
+    rcipsav.R = cell(nsub+1,1);
+    rcipsav.MAT = cell(nsub,1);
+    rcipsav.chnkrlocals = cell(nsub,1);
+end
+
+
+% grab only those kernels relevant to this vertex
 
-km1 = k-1;
-rcs = zeros(km1,dim,nedge);
-dcs = zeros(k,dim,nedge);
-d2cs = zeros(k,dim,nedge);
-dscal = zeros(nedge,1);
-d2scal = zeros(nedge,1);
-ctr = zeros(dim,nedge);
 if(size(fkern)==1)
     fkernlocal = fkern;
 else
@@ -55,6 +82,20 @@
 
 end
 
+rcipsav.fkernlocal = fkernlocal;
+
+% get coefficients of recentered edge chunks and figure out orientation
+
+km1 = k-1;
+rcs = zeros(km1,dim,nedge);
+dcs = zeros(k,dim,nedge);
+d2cs = zeros(k,dim,nedge);
+dscal = zeros(nedge,1);
+d2scal = zeros(nedge,1);
+ctr = zeros(dim,nedge);
+
+ileftright = zeros(nedge,1);
+nextchunk = zeros(nedge,1);
 
 for i = 1:nedge
     ic = iedgechunks(1,i);
@@ -92,10 +133,19 @@
     end
 end
 
+rcipsav.ctr = ctr;
+rcipsav.rcs = rcs;
+rcipsav.dcs = dcs;
+rcipsav.d2cs = d2cs;
+rcipsav.dscal = dscal;
+rcipsav.d2scal = d2scal;
+rcipsav.ileftright = ileftright;
+rcipsav.glxs = glxs;
+rcipsav.glws = glws;
 
 pref = []; 
 pref.k = k;
-pref.nchmax = 5;
+pref.nchstor = 5;
 
 R = [];
 
@@ -108,6 +158,9 @@
 ts = cell(nedge,1);
 chnkrlocal(1,nedge) = chunker();
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% begin recursion proper
+
 h0=ones(nedge,1);
 for level=1:nsub
     h = h0/2^(nsub-level);
@@ -176,8 +229,16 @@
     if level==1    %  Dumb and lazy initializer for R, for now
   %R=eye(nR); 
         R = inv(MAT(starL,starL));
+        if savedeep
+            rcipsav.R{1} = R;
+        end
     end
     R=chnk.rcip.SchurBana(Pbc,PWbc,MAT,R,starL,circL,starS,circS);   
+    if savedeep
+        rcipsav.R{level+1} = R;
+        rcipsav.MAT{level} = MAT(starL,circL);
+        rcipsav.chnkrlocals{level} = merge(chnkrlocal);
+    end
 end
 
 end