Merge pull request #58 from fastalgorithms/dev-chunkerkerneval

Dev chunkerkerneval

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/@chunkgraph/chunkgraph.m

@@ -1,9 +1,16 @@
 classdef chunkgraph
 %CHUNKGRAPH chunk graph class for storing complex domains
 %
-% We describe a complex domain by its edges (smooth, regular
-% curves) and vertices (free ends of edges or points where the ends of 
-% one or more edges meets)
+% We describe a complex domain by a planar graph. Smooth boundary components
+% (discretized by chunkers) specify the edges of the graph and the
+% coordinates of free ends of the edges or points where the ends of multiple
+% edges meet specify the vertices. The interiors of distinct edges are not 
+% allowed to intersect (i.e. intersecting curves should be split at the
+% point where they intersect. 
+%
+% The chunkgraph is specified by an array of chunkers called echnks (the
+% edges), an array of points called verts (the vertices), and an array of
+% indices specifying the vertices at the left and right ends of any 
 %
 % 
     properties(SetAccess=public)

chunkie/@kernel/helm2d.m

@@ -14,6 +14,10 @@
 %   parameter ETA, i.e., COEFS(1)*KERNEL.HELM2D('d', ZK) + 
 %   COEFS(2)*KERNEL.HELM2D('s', ZK).
 %
+%   KERNEL.HELM2D('cp', ZK, COEFS) or KERNEL.HELM2D('combined', ZK, COEFS)
+%   constructs the derivative of the combined-layer Helmholtz kernel with 
+%   wavenumber ZK and  parameter ETA, i.e., COEFS(1)*KERNEL.HELM2D('dp', ZK) + 
+%   COEFS(2)*KERNEL.HELM2D('sp', ZK).
 % See also CHNK.HELM2D.KERN.
 
 if ( nargin < 1 )
@@ -66,6 +70,17 @@
         obj.fmm  = @(eps,s,t,sigma) chnk.helm2d.fmm(eps, zk, s, t, 'c', sigma, coefs);
         obj.sing = 'log';
 
+    case {'cp', 'cprime'}
+        if ( nargin < 3 )
+            warning('Missing combined layer coefficients. Defaulting to [1,1i].');
+            coefs = [1,1i];
+        end
+        obj.type = 'cprime';
+        obj.params.coefs = coefs;
+        obj.eval = @(s,t) chnk.helm2d.kern(zk, s, t, 'cprime', coefs);
+        obj.fmm  = @(eps,s,t,sigma) chnk.helm2d.fmm(eps, zk, s, t, 'cprime', sigma, coefs);
+        obj.sing = 'hs';
+
     otherwise
         error('Unknown Helmholtz kernel type ''%s''.', type);
 

chunkie/@kernel/kernel.m

@@ -9,6 +9,7 @@
 %      ----                                         ----
 %      'laplace'    ('lap', 'l')                    's', 'd', 'sp', 'c'
 %      'helmholtz'  ('helm', 'h')                   's', 'd', 'sp', 'dp', 'c'
+%                                                   'cp'
 %      'helmholtz difference' ('helmdiff', 'hdiff') 's', 'd', 'sp', 'dp'
 %      'elasticity' ('elast', 'e')                  's', 'strac', 'd', 'dalt'
 %      'stokes'     ('stok', 's')                   'svel', 'spres', 'strac',

chunkie/chunkerinterior.m

@@ -1,16 +1,18 @@
-function [in] = chunkerinterior(chnkr,pts,opts)
+function [in] = chunkerinterior(chnkobj,ptsobj,opts)
 %CHUNKERINTERIOR returns an array indicating whether each point specified 
 % % by pts is inside the domain. Assumes the domain is closed.
 %
-% Syntax: in = chunkerinterior(chnkr,pts,opts)
-%         in = chunkerinterior(chnkr,{x,y},opts) % meshgrid version
+% Syntax: in = chunkerinterior(chnkobj,pts,opts)
+%         in = chunkerinterior(chnkobj,{x,y},opts) % meshgrid version
 %
 % Input:
-%   chnkr - chunker object describing geometry
-%   pts - (chnkr.dim,:) array of points to test
-% 
-%   {x,y} - length 2 cell array. the points checked then have the
-%       coordinates of a mesh grid [xx,yy] = meshgrid(x,y)
+%   chnkobj - chunker object or chunkgraph object describing geometry
+%   ptsobj - object describing the target points, can be specified as
+%       * (chnkr.dim,:) array of points to test
+%       * {x,y} - length 2 cell array. the points checked then have the
+%           coordinates of a mesh grid [xx,yy] = meshgrid(x,y)
+%       * chunker object, in which case it uses chunker.r(:,:) 
+%       * chunkgraph object, in which case it uses chunkgraph.r(:,:) 
 %
 % Optional input:
 %   opts - options structure with entries:
@@ -34,23 +36,39 @@
 
 grid = false;
 
-assert(chnkr.dim == 2,'interior only well-defined for 2D');
 
 if nargin < 3
     opts = [];
 end
 
-if isa(pts,"cell")
-    assert(length(pts)==2,'second input should be either 2xnpts array or length 2 cell array');
-    x = pts{1};
-    y = pts{2};
-    grid = true;
+% Assign appropriate object to chnkr
+if class(chnkobj) == "chunker"
+   chnkr = chnkobj;
+elseif class(chnkobj) == "chunkgraph"
+   chnkr = merge(chnkobj.echnks);
+else
+    msg = "Unsupported object in chunkerinterior";
+    error(msg)
 end
 
-if grid
+assert(chnkr.dim == 2,'interior only well-defined for 2D');
+
+% Assign appropriate object to pts
+if isa(ptsobj, "cell")
+    assert(length(ptsobj)==2,'second input should be either 2xnpts array or length 2 cell array');
+    x = ptsobj{1};
+    y = ptsobj{2};
+    grid = true;
     [xx,yy] = meshgrid(x,y);
     pts = [xx(:).'; yy(:).'];
-end    
+elseif isa(ptsobj, "chunker")
+    pts = ptsobj.r(:,:);
+elseif isa(ptsobj, "chunkgraph")
+    pts = ptsobj.r(:,:);
+else
+    pts = ptsobj;
+end
+
 
 usefmm = true;
 if isfield(opts,'fmm')