diff --git a/chunkie/+chnk/+flam/kernbyindexr.m b/chunkie/+chnk/+flam/kernbyindexr.m
index 8f6b1eb..b280657 100644
--- a/chunkie/+chnk/+flam/kernbyindexr.m
+++ b/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);
 
diff --git a/chunkie/+chnk/+flam/proxy_square_pts.m b/chunkie/+chnk/+flam/proxy_square_pts.m
index 1d7bb7f..3f1d362 100644
--- a/chunkie/+chnk/+flam/proxy_square_pts.m
+++ b/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
\ No newline at end of file
diff --git a/chunkie/+chnk/+helm2d/fmm.m b/chunkie/+chnk/+helm2d/fmm.m
index 4d50789..d7eceb0 100644
--- a/chunkie/+chnk/+helm2d/fmm.m
+++ b/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);
diff --git a/chunkie/+chnk/+helm2d/kern.m b/chunkie/+chnk/+helm2d/kern.m
index a90e636..d697efd 100644
--- a/chunkie/+chnk/+helm2d/kern.m
+++ b/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;
 
diff --git a/chunkie/+chnk/+quadnative/buildmat.m b/chunkie/+chnk/+quadnative/buildmat.m
index c75873a..7dbbd01 100644
--- a/chunkie/+chnk/+quadnative/buildmat.m
+++ b/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;
diff --git a/chunkie/@chunkgraph/chunkgraph.m b/chunkie/@chunkgraph/chunkgraph.m
index f42edea..7a57465 100644
--- a/chunkie/@chunkgraph/chunkgraph.m
+++ b/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)
diff --git a/chunkie/@kernel/helm2d.m b/chunkie/@kernel/helm2d.m
index 64bed00..73cc8ce 100644
--- a/chunkie/@kernel/helm2d.m
+++ b/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);
 
diff --git a/chunkie/@kernel/kernel.m b/chunkie/@kernel/kernel.m
index f824811..30803a6 100644
--- a/chunkie/@kernel/kernel.m
+++ b/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',
diff --git a/chunkie/FLAM b/chunkie/FLAM
index 2c9361e..a83698e 160000
--- a/chunkie/FLAM
+++ b/chunkie/FLAM
@@ -1 +1 @@
-Subproject commit 2c9361e7bc8fb8c3255f33275c0f2f7b7bbb7572
+Subproject commit a83698efdd58293c0e4c0c52f59498ffb9da9a99
diff --git a/chunkie/chunkerinterior.m b/chunkie/chunkerinterior.m
index 9563985..cde9605 100644
--- a/chunkie/chunkerinterior.m
+++ b/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')
diff --git a/chunkie/chunkerkerneval.m b/chunkie/chunkerkerneval.m
index 24d068c..fb10f27 100644
--- a/chunkie/chunkerkerneval.m
+++ b/chunkie/chunkerkerneval.m
@@ -1,17 +1,20 @@
-function fints = chunkerkerneval(chnkr,kern,dens,targs,opts)
+function fints = chunkerkerneval(chnkobj,kern,dens,targobj,opts)
 %CHUNKERKERNEVAL compute the convolution of the integral kernel with
 % the density defined on the chunk geometry. 
 %
 % Syntax: fints = chunkerkerneval(chnkr,kern,dens,targs,opts)
 %
 % Input:
-%   chnkr - chunker object description of curve
+%   chnkobj - chunker object or chunkgraph object description of curve
 %   kern - kernel class object or kernel function taking inputs 
 %                      kern(srcinfo,targinfo) 
 %   dens - density on boundary, should have size opdims(2) x k x nch
 %          where k = chnkr.k, nch = chnkr.nch, where opdims is the 
 %           size of kern for a single src,targ pair
-%   targs - targ(1:2,i) gives the coords of the ith target
+%   targobj - object describing the target points, can be specified as
+%       * array of points
+%       * chunker object
+%       * chunkgraph object
 %
 % Optional input:
 %   opts - structure for setting various parameters
@@ -54,6 +57,16 @@
     kerneval = kern;
 end
 
+% Assign appropriate object to chnkr
+if class(chnkobj) == "chunker"
+   chnkr = chnkobj;
+elseif class(chnkobj) == "chunkgraph"
+   chnkr = merge(chnkobj.echnks);
+else
+    msg = "Unsupported object in chunkerkerneval";
+    error(msg)
+end
+
 srcinfo = []; targinfo = [];
 srcinfo.r = chnkr.r(:,1); srcinfo.d = chnkr.d(:,1); 
 srcinfo.n = chnkr.n(:,1); srcinfo.d2 = chnkr.d2(:,1);
@@ -86,61 +99,77 @@
 if isfield(opts,'fac'); opts_use.fac = opts.fac; end
 if isfield(opts,'eps'); opts_use.eps = opts.eps; end
 
-[dim,~] = size(targs);
+% Assign appropriate object to targinfo
+targinfo = [];
+if isa(targobj, "chunker")
+    targinfo.r = targobj.r(:,:);
+    targinfo.d = targobj.d(:,:);
+    targinfo.d2 = targobj.d2(:,:);
+    targinfo.n = targobj.n(:,:);
+elseif isa(targobj, "chunkgraph")
+    targinfo.r = targobj.r(:,:);
+    targinfo.d = targobj.d(:,:);
+    targinfo.d2 = targobj.d2(:,:);
+    targinfo.n = targobj.n(:,:);
+else
+    targinfo.r = targobj;
+end
+
+
+
+[dim,~] = size(targinfo.r);
 
 if (dim ~= 2); warning('only dimension two tested'); end
 
 if opts_use.forcesmooth
-    fints = chunkerkerneval_smooth(chnkr,kern,opdims,dens,targs, ...
+    fints = chunkerkerneval_smooth(chnkr,kern,opdims,dens,targinfo, ...
         [],opts_use);
     return
 end
 
 if opts_use.forceadap
     fints = chunkerkerneval_adap(chnkr,kern,opdims,dens, ...
-        targs,[],opts_use);
+        targinfo,[],opts_use);
     return
 end
 
 
 if opts_use.forcepquad
     optsflag = []; optsflag.fac = opts_use.fac;
-    flag = flagnear(chnkr,targs,optsflag);
-    fints = chunkerkerneval_smooth(chnkr,kern,opdims,dens,targs, ...
+    flag = flagnear(chnkr,targinfo.r,optsflag);
+    fints = chunkerkerneval_smooth(chnkr,kern,opdims,dens,targinfo, ...
         flag,opts_use);
 
     fints = fints + chunkerkerneval_ho(chnkr,kern,opdims,dens, ...
-        targs,flag,opts_use);
+        targinfo,flag,opts_use);
 
     return
 end
 
 % smooth for sufficiently far, adaptive otherwise
 
-%optsflag = []; optsflag.fac = opts_use.fac;
 rho = 1.8;
 optsflag = [];  optsflag.rho = rho;
-flag = flagnear_rectangle(chnkr,targs,optsflag);
+flag = flagnear_rectangle(chnkr,targinfo.r,optsflag);
 
 npoly = chnkr.k*2;
 nlegnew = chnk.ellipse_oversample(rho,npoly,opts_use.eps);
 nlegnew = max(nlegnew,chnkr.k);
 
 [chnkr2,dens2] = upsample(chnkr,nlegnew,dens);
-fints = chunkerkerneval_smooth(chnkr2,kern,opdims,dens2,targs, ...
+fints = chunkerkerneval_smooth(chnkr2,kern,opdims,dens2,targinfo, ...
     flag,opts_use);
 
 fints = fints + chunkerkerneval_adap(chnkr,kern,opdims,dens, ...
-        targs,flag,opts_use);
-
+        targinfo,flag,opts_use);
 
 end
 
 function fints = chunkerkerneval_ho(chnkr,kern,opdims,dens, ...
-    targs,flag,opts)
+    targinfo,flag,opts)
 
 % target
-[~,nt] = size(targs);
+[~,nt] = size(targinfo.r);
 fints = zeros(opdims(1)*nt,1);
 k = size(chnkr.r,2);
 [t,w] = lege.exps(2*k);
@@ -155,8 +184,22 @@
 % interpolation matrix
 intp = lege.matrin(k,t);          % interpolation from k to 2*k
 intp_ab = lege.matrin(k,[-1;1]);  % interpolation from k to end points
+
+targs = targinfo.r;
+
 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
 for j=1:size(chnkr.r,3)
     [ji] = find(flag(:,j));
     if(~isempty(ji))
@@ -174,7 +217,7 @@
 
 
 function fints = chunkerkerneval_smooth(chnkr,kern,opdims,dens, ...
-    targs,flag,opts)
+    targinfo,flag,opts)
 
 if isa(kern,'kernel')
     kerneval = kern.eval;
@@ -203,12 +246,10 @@
 dens = reshape(dens,opdims(2),k,nch);
 
 [~,w] = lege.exps(k);
-[~,nt] = size(targs);
+[~,nt] = size(targinfo.r);
 
 fints = zeros(opdims(1)*nt,1);
 
-targinfo = []; targinfo.r = targs;
-
 % assume smooth weights are good enough
 
 % Sequence of checks, first see ifflam is set as it supercedes
@@ -272,24 +313,6 @@
     end
 else
 
-%     % hack to ignore interactions: create sparse array with very small
-%     % number instead of zero in ignored entries. kernbyindexr overwrites
-%     %   with that small number
-%     [i,j] = find(flag);
-%     
-%     % targets correspond to multiple outputs (if matrix valued kernel)
-%     inew = (opdims(1)*(i(:)-1)).' + (1:opdims(1)).'; inew = inew(:);
-%     jnew = (j(:)).' + 0*(1:opdims(1)).'; jnew = jnew(:);
-%     
-%     % source chunks have multiple points and can be multi-dimensional
-%     jnew = (opdims(2)*chnkr.k*(jnew(:)-1)).' + (1:(chnkr.k*opdims(2))).';
-%     jnew = jnew(:);
-%     inew = (inew(:)).' + 0*(1:(chnkr.k*opdims(2))).';
-%     inew = inew(:);
-%     
-%     mm = nt*opdims(1); nn = chnkr.npt*opdims(2);
-%     v = 1e-300*ones(length(inew),1); sp = sparse(inew,jnew,v,mm,nn);
-
     wts = chnkr.wts;
     wts = wts(:);
     
@@ -297,14 +320,29 @@
         xflam1 = chnkr.r(:,:);
         xflam1 = repmat(xflam1,opdims(2),1);
         xflam1 = reshape(xflam1,chnkr.dim,numel(xflam1)/chnkr.dim);
-        targsflam = repmat(targs(:,:),opdims(1),1);
-        targsflam = reshape(targsflam,chnkr.dim,numel(targsflam)/chnkr.dim);
-        matfun = @(i,j) chnk.flam.kernbyindexr(i,j,targs,chnkr,kerneval, ...
-            opdims);
+
+        targinfo_flam = [];
+        targinfo_flam.r = repelem(targinfo.r(:,:),1,opdims(1));
+        if isfield(targinfo, 'd')
+            targinfo_flam.d = repelem(targinfo.d(:,:),1,opdims(1));
+        end
+        
+        if isfield(targinfo, 'd2')
+            targinfo_flam.d2 = repelem(targinfo.d2(:,:),1,opdims(1));
+        end
+        
+        if isfield(targinfo, 'n')
+            targinfo_flam.n = repelem(targinfo.n(:,:),1,opdims(1));
+        end
+
+% TODO: Pull through data?
+
+        matfun = @(i,j) chnk.flam.kernbyindexr(i, j, targinfo_flam, ...,
+                           chnkr, kerneval, opdims);
     
 
         width = max(abs(max(chnkr)-min(chnkr)))/3;
-        tmax = max(targs(:,:),[],2); tmin = min(targs(:,:),[],2);
+        tmax = max(targinfo.r(:,:),[],2); tmin = min(targinfo.r(:,:),[],2);
         wmax = max(abs(tmax-tmin));
         width = max(width,wmax/3);  
         npxy = chnk.flam.nproxy_square(kerneval,width);
@@ -314,16 +352,16 @@
             ctr,chnkr,kerneval,opdims,pr,ptau,pw,pin);
 
         optsifmm=[]; optsifmm.Tmax=Inf;
-        F = ifmm(matfun,targsflam,xflam1,200,1e-14,pxyfun,optsifmm);
+        F = ifmm(matfun,targinfo_flam.r,xflam1,200,1e-14,pxyfun,optsifmm);
         fints = ifmm_mv(F,dens(:),matfun);
     else
         wts2 = repmat(wts(:).', opdims(2), 1);
         sigma = wts2(:).*dens(:);
-        fints = kern.fmm(1e-14, chnkr, targs(:,:), sigma);
+        fints = kern.fmm(1e-14, chnkr, targinfo.r(:,:), sigma);
     end
     % delete interactions in flag array (possibly unstable approach)
     
-    targinfo = [];
+    
     if ~isempty(flag)
         for i = 1:nch
             densvals = dens(:,:,i); densvals = densvals(:);
@@ -338,9 +376,24 @@
             delsmooth = find(flag(:,i)); 
             delsmoothrow = (opdims(1)*(delsmooth(:)-1)).' + (1:opdims(1)).';
             delsmoothrow = delsmoothrow(:);
-            targinfo.r = targs(:,delsmooth);
 
-            kernmat = kerneval(srcinfo,targinfo);
+            targinfo_use = [];
+            targinfo_use.r = targinfo.r(:,delsmooth);
+
+            if isfield(targinfo, 'd')
+                targinfo_use.d = targinfo.d(:,delsmooth);
+            end
+        
+            if isfield(targinfo, 'd2')
+                targinfo_use.d2 = targinfo.d2(:,delsmooth);
+            end
+        
+            if isfield(targinfo, 'n')
+                targinfo_use.n = targinfo.n(:,delsmooth);
+            end
+
+
+            kernmat = kerneval(srcinfo,targinfo_use);
             fints(delsmoothrow) = fints(delsmoothrow) - kernmat*densvals;
         end
     end    
@@ -349,7 +402,7 @@
 end
 
 function fints = chunkerkerneval_adap(chnkr,kern,opdims,dens, ...
-    targs,flag,opts)
+    targinfo,flag,opts)
 
 if isa(kern,'kernel')
     kerneval = kern.eval;
@@ -370,7 +423,23 @@
 assert(numel(dens) == opdims(2)*k*nch,'dens not of appropriate size')
 dens = reshape(dens,opdims(2),k,nch);
 
+% 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
+
 
 fints = zeros(opdims(1)*nt,1);
 
@@ -384,8 +453,8 @@
 n = chnkr.n;
 d2 = chnkr.d2;
 h = chnkr.h;
-targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt);
-targn = zeros(chnkr.dim,nt);
+
+
 
 if isempty(flag) % do all to all adaptive
     for i = 1:nch
diff --git a/chunkie/chunkerkernevalmat.m b/chunkie/chunkerkernevalmat.m
index c0a96c9..a697ac2 100644
--- a/chunkie/chunkerkernevalmat.m
+++ b/chunkie/chunkerkernevalmat.m
@@ -1,4 +1,4 @@
-function mat = chunkerkernevalmat(chnkr,kern,targs,opts)
+function mat = 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
@@ -6,9 +6,21 @@
 % Syntax: mat = chunkerkernevalmat(chnkr,kern,targs,opts)
 %
 % Input:
-%   chnkr - chunker object description of curve
-%   kern - integral kernel taking inputs kern(srcinfo,targinfo) 
-%   targs - targ(1:2,i) gives the coords of the ith target
+%   chnkr - chunker object describing boundary, currently
+%              only supports chunkers, and not chunkgraphs
+%   kern  - kernel function. By default, this should be a function handle
+%           accepting input of the form kern(srcinfo,targinfo), where srcinfo
+%           and targinfo are in the ptinfo struct format, i.e.
+%                ptinfo.r - positions (2,:) array
+%                ptinfo.d - first derivative in underlying
+%                     parameterization (2,:)
+%                ptinfo.n - unit normals (2,:)
+%                ptinfo.d2 - second derivative in underlying
+%                     parameterization (2,:)
+%   targobj - object describing the target points, can be specified as
+%       * array of points
+%       * chunker object
+%       * chunkgraph object
 %
 % Optional input:
 %   opts - structure for setting various parameters
@@ -38,6 +50,37 @@
 
 % author: Travis Askham (askhamwhat@gmail.com)
 
+% convert kernel to kernel object, put in singularity info 
+% opts.sing provides a default value for singularities if not 
+% defined for kernels
+
+if isa(kern,'function_handle')
+    kern2 = kernel(kern);
+    kern = kern2;
+elseif isa(kern,'cell')
+    sz = size(kern);
+    kern2(sz(1),sz(2)) = kernel();
+    for j = 1:sz(2)
+        for i = 1:sz(1)
+            if isa(kern{i,j},'function_handle')
+                kern2(i,j) = kernel(kern{i,j});
+            elseif isa(kern{i,j},'kernel')
+                kern2(i,j) = kern{i,j};
+            else
+                msg = "Second input is not a kernel object, function handle, " ...
+                    + "or cell array";
+                error(msg);
+            end
+        end
+    end
+    kern = kern2;
+    
+elseif ~isa(kern,'kernel')
+    msg = "Second input is not a kernel object, function handle, " ...
+                + "or cell array";
+    error(msg);
+end
+    
 % determine operator dimensions using first two points
 
 
@@ -48,7 +91,9 @@
 targinfo.r = chnkr.r(:,2); targinfo.d = chnkr.d(:,2); 
 targinfo.d2 = chnkr.d2(:,2); targinfo.n = chnkr.n(:,2);
 
-ftemp = kern(srcinfo,targinfo);
+ftmp = kern.eval;
+
+ftemp = ftmp(srcinfo,targinfo);
 opdims = size(ftemp);
 
 if nargin < 4
@@ -57,18 +102,34 @@
 
 forcesmooth = false;
 forceadap = false;
-forcepqud = false;
+forcepquad = false;
 nonsmoothonly = false;
 fac = 1.0;
 eps = 1e-12;
 if isfield(opts,'forcesmooth'); forcesmooth = opts.forcesmooth; end
 if isfield(opts,'forceadap'); forceadap = opts.forceadap; end
-if isfield(opts,'forcepquad'); forcepqud = opts.forcepquad; end
+if isfield(opts,'forcepquad'); forcepquad = opts.forcepquad; end
 if isfield(opts,'nonsmoothonly'); nonsmoothonly = opts.nonsmoothonly; end
 if isfield(opts,'fac'); fac = opts.fac; end
 if isfield(opts,'eps'); eps = opts.eps; end
 
-[dim,~] = size(targs);
+% Assign appropriate object to targinfo
+targinfo = [];
+if isa(targobj, "chunker")
+    targinfo.r = targobj.r(:,:);
+    targinfo.d = targobj.d(:,:);
+    targinfo.d2 = targobj.d2(:,:);
+    targinfo.n = targobj.n(:,:);
+elseif isa(targobj, "chunkgraph")
+    targinfo.r = targobj.r(:,:);
+    targinfo.d = targobj.d(:,:);
+    targinfo.d2 = targobj.d2(:,:);
+    targinfo.n = targobj.n(:,:);
+else
+    targinfo.r = targobj;
+end
+
+[dim,~] = size(targinfo.r);
 
 
 if (dim ~= 2); warning('only dimension two tested'); end
@@ -77,24 +138,26 @@
 optsadap = []; 
 optsadap.eps = eps;
 
+
+
 if forcesmooth
-    mat = chunkerkernevalmat_smooth(chnkr,kern,opdims,targs, ...
+    mat = chunkerkernevalmat_smooth(chnkr,ftmp,opdims,targinfo, ...
         [],optssmooth);
     return
 end
 
 if forceadap
-    mat = chunkerkernevalmat_adap(chnkr,kern,opdims, ...
-        targs,[],optsadap);
+    mat = chunkerkernevalmat_adap(chnkr,ftmp,opdims, ...
+        targinfo,[],optsadap);
     return
 end
 
-if forcepqud
+if forcepquad
     optsflag = []; optsflag.fac = fac;
-    flag = flagnear(chnkr,targs,optsflag);
-    spmat = chunkerkernevalmat_ho(chnkr,kern,opdims, ...
-        targs,flag,optsadap);
-    mat = chunkerkernevalmat_smooth(chnkr,kern,opdims,targs, ...
+    flag = flagnear(chnkr,targinfo.r,optsflag);
+    spmat = chunkerkernevalmat_ho(chnkr,ftmp,opdims, ...
+        targinfo,flag,optsadap);
+    mat = chunkerkernevalmat_smooth(chnkr,ftmp,opdims,targinfo, ...
         flag,opts);
     mat = mat + spmat;
     return
@@ -102,17 +165,20 @@
 
 % smooth for sufficiently far, adaptive otherwise
 
+% TODO: change to chunkerkerneval system, need routine to generate
+% upsampling matrix.
+
 optsflag = []; optsflag.fac = fac;
-flag = flagnear(chnkr,targs,optsflag);
-spmat = chunkerkernevalmat_adap(chnkr,kern,opdims, ...
-        targs,flag,optsadap);
+flag = flagnear(chnkr,targinfo.r,optsflag);
+spmat = chunkerkernevalmat_adap(chnkr,ftmp,opdims, ...
+        targinfo,flag,optsadap);
 
 if nonsmoothonly
     mat = spmat;
     return;
 end
 
-mat = chunkerkernevalmat_smooth(chnkr,kern,opdims,targs, ...
+mat = chunkerkernevalmat_smooth(chnkr,ftmp,opdims,targinfo, ...
     flag,opts);
 
 mat = mat + spmat;
@@ -123,7 +189,7 @@
 
 
 function mat = chunkerkernevalmat_smooth(chnkr,kern,opdims, ...
-    targs,flag,opts)
+    targinfo,flag,opts)
 
 if nargin < 6
     flag = [];
@@ -135,7 +201,6 @@
 k = chnkr.k;
 nch = chnkr.nch;
 
-targinfo = []; targinfo.r = targs;
 srcinfo = []; srcinfo.r = chnkr.r(:,:); srcinfo.n = chnkr.n(:,:);
 srcinfo.d = chnkr.d(:,:); srcinfo.d2 = chnkr.d2(:,:);
 
@@ -163,7 +228,7 @@
 end
 
 function mat = chunkerkernevalmat_adap(chnkr,kern,opdims, ...
-    targs,flag,opts)
+    targinfo,flag,opts)
 
 k = chnkr.k;
 nch = chnkr.nch;
@@ -175,12 +240,30 @@
     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
+
 
 % using adaptive quadrature
 
 
 if isempty(flag)
+
+
     mat = zeros(opdims(1)*nt,opdims(2)*chnkr.npt);
 
     [t,w] = lege.exps(2*k+1);
@@ -191,7 +274,7 @@
     n = chnkr.n;
     d2 = chnkr.d2;
     h = chnkr.h;
-    targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt);    
+    
     for i = 1:nch
         jmat = 1 + (i-1)*k*opdims(2);
         jmatend = i*k*opdims(2);
@@ -228,8 +311,6 @@
     n = chnkr.n;
     d2 = chnkr.d2;
     h = chnkr.h;
-    targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt);
-    targn = zeros(chnkr.dim,nt);
     for i = 1:nch
         jmat = 1 + (i-1)*k*opdims(2);
         jmatend = i*k*opdims(2);
@@ -260,7 +341,7 @@
 end
 
 function mat = chunkerkernevalmat_ho(chnkr,kern,opdims, ...
-    targs,flag,opts)
+    targinfo,flag,opts)
 
 k = chnkr.k;
 nch = chnkr.nch;
@@ -272,7 +353,24 @@
     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
+
+
 
 % using Helsing-Ojala quadrature
 if isempty(flag) % figure out what is this flag for in adaptive routine
@@ -295,8 +393,7 @@
     % interpolation matrix 
     intp = lege.matrin(k,t);          % interpolation from k to 2*k
     intp_ab = lege.matrin(k,[-1;1]);  % interpolation from k to end points
-    targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt);
-    targn = zeros(chnkr.dim,nt);
+    
     for i = 1:nch
         jmat = 1 + (i-1)*k*opdims(2);
         jmatend = i*k*opdims(2);
diff --git a/chunkie/chunkermat.m b/chunkie/chunkermat.m
index c60f3c8..25f8bf2 100644
--- a/chunkie/chunkermat.m
+++ b/chunkie/chunkermat.m
@@ -7,7 +7,7 @@
 % Syntax: sysmat = chunkermat(chnkr,kern,opts)
 %
 % Input:
-%   chnkr - chunker object describing boundary
+%   chnkobj - chunker object describing boundary
 %   kern  - kernel function. By default, this should be a function handle
 %           accepting input of the form kern(srcinfo,targinfo), where srcinfo
 %           and targinfo are in the ptinfo struct format, i.e.
diff --git a/chunkie/trapperkerneval.m b/chunkie/trapperkerneval.m
index d425b26..aaf7939 100644
--- a/chunkie/trapperkerneval.m
+++ b/chunkie/trapperkerneval.m
@@ -1,10 +1,21 @@
-function fints = trapperkerneval(trap,kern,dens,targs,opts)
+function fints = trapperkerneval(trap,kern,dens,targobj,opts)
 %TRAPPERINTKERN compute the convolution of the integral kernel with
 
 wts = trap.wts;
 srcinfo = []; srcinfo.r = trap.r; srcinfo.d = trap.d; srcinfo.d2 = trap.d2;
 srcinfo.n = trap.n;
-targinfo = []; targinfo.r = targs;
+
+% Assign appropriate object to targinfo
+targinfo = [];
+if isa(targobj, "trapper")
+    targinfo.r = targobj.r(:,:);
+    targinfo.d = targobj.d(:,:);
+    targinfo.d2 = targobj.d2(:,:);
+    targinfo.n = targobj.n(:,:);
+else
+    targinfo.r = targobj;
+end
+
 mat = kern(srcinfo,targinfo);
 fints = mat*diag(wts)*dens;
 
diff --git a/devtools/test/chunkerinteriorTest.m b/devtools/test/chunkerinteriorTest.m
index 65a6a91..2538be6 100644
--- a/devtools/test/chunkerinteriorTest.m
+++ b/devtools/test/chunkerinteriorTest.m
@@ -67,5 +67,20 @@
 opts = [];
 opts.fmm = true;
 opts.flam = false;
-start = tic; in3 = chunkerinterior(chnkr,{x1,x1},opts); t4 = toc(start);
-fprintf('%5.2e s : time for chunkerinterior (meshgrid, with fmm)\n',t4);
\ No newline at end of file
+start = tic; in4 = chunkerinterior(chnkr,{x1,x1},opts); t4 = toc(start);
+fprintf('%5.2e s : time for chunkerinterior (meshgrid, with fmm)\n',t4);
+
+% Test targets specified as chunker
+narms = 3;
+amp = 0.1;
+chnkr2 = chunkerfunc(@(t) 0.3*starfish(t,narms,amp),cparams,pref); 
+
+opts = [];
+opts.fmm = true;
+opts.flam = false;
+start = tic; in5 = chunkerinterior(chnkr,chnkr2,opts); t5 = toc(start);
+fprintf('%5.2e s : time for chunkerinterior (chunker, with fmm)\n',t5);
+assert(all(in5(:) == 1));
+
+
+
diff --git a/devtools/test/chunkerkerneval_greenlapTest.m b/devtools/test/chunkerkerneval_greenlapTest.m
index 9cccf63..ecf5c92 100644
--- a/devtools/test/chunkerkerneval_greenlapTest.m
+++ b/devtools/test/chunkerkerneval_greenlapTest.m
@@ -23,7 +23,7 @@
 
 % sources
 
-ns = 10;
+ns = 100;
 ts = 0.0+2*pi*rand(ns,1);
 sources = starfish(ts,narms,amp);
 sources = 3.0*sources;
@@ -31,7 +31,7 @@
 
 % targets
 
-nt = 100;
+nt = 300;
 ts = 0.0+2*pi*rand(nt,1);
 targets = starfish(ts,narms,amp);
 targets = targets.*repmat(rand(1,nt),2,1);
@@ -52,9 +52,9 @@
 
 % kernel defs
 
-kernd = @(s,t) chnk.lap2d.kern(s,t,'d');
-kerns = @(s,t) chnk.lap2d.kern(s,t,'s');
-kernsprime = @(s,t) chnk.lap2d.kern(s,t,'sprime');
+kernd = kernel('lap','d');
+kerns = kernel('lap','s');
+kernsprime = kernel('lap','sprime');
 
 opdims = [1 1];
 
@@ -63,21 +63,22 @@
 srcinfo = []; srcinfo.r = sources; 
 targinfo = []; targinfo.r = chnkr.r(:,:); 
 targinfo.d = chnkr.d(:,:);
-kernmats = kerns(srcinfo,targinfo);
-kernmatsprime = kernsprime(srcinfo,targinfo);
+kernmats = kerns.eval(srcinfo,targinfo);
+kernmatsprime = kernsprime.eval(srcinfo,targinfo);
 densu = kernmats*strengths;
 densun = kernmatsprime*strengths;
 
 % eval u at targets
 
 targinfo = []; targinfo.r = targets;
-kernmatstarg = kerns(srcinfo,targinfo);
+kernmatstarg = kerns.eval(srcinfo,targinfo);
 utarg = kernmatstarg*strengths;
 
 
-% test green's id
+% test green's id. we use FLAM, direct and FMM for smooth work
 
-opts.quadkgparams = {'RelTol',1.0e-13,'AbsTol',1.0e-13};
+opts = [];
+opts.flam = true;
 start=tic; Du = chunkerkerneval(chnkr,kernd,densu,targets,opts); 
 toc(start)
 start=tic; Sun = chunkerkerneval(chnkr,kerns,densun,targets,opts); 
@@ -89,7 +90,41 @@
 
 relerr = norm(utarg-utarg2,'fro')/norm(utarg,'fro');
 
-fprintf('relative frobenius error %5.2e\n',relerr);
+fprintf('relative frobenius error, forcing flam %5.2e\n',relerr);
+
+assert(relerr < 1e-11);
+
+opts = [];
+opts.accel = false;
+start=tic; Du = chunkerkerneval(chnkr,kernd,densu,targets,opts); 
+toc(start)
+start=tic; Sun = chunkerkerneval(chnkr,kerns,densun,targets,opts); 
+toc(start)
+
+utarg2 = Sun-Du;
+
+%
+
+relerr = norm(utarg-utarg2,'fro')/norm(utarg,'fro');
+
+fprintf('relative frobenius error, forcing direct %5.2e\n',relerr);
+
+assert(relerr < 1e-11);
+
+opts = [];
+opts.forcefmm = true;
+start=tic; Du = chunkerkerneval(chnkr,kernd,densu,targets,opts); 
+toc(start)
+start=tic; Sun = chunkerkerneval(chnkr,kerns,densun,targets,opts); 
+toc(start)
+
+utarg2 = Sun-Du;
+
+%
+
+relerr = norm(utarg-utarg2,'fro')/norm(utarg,'fro');
+
+fprintf('relative frobenius error, forcing fmm %5.2e\n',relerr);
 
 assert(relerr < 1e-11);
 
diff --git a/devtools/test/chunkgrphconstructTest.m b/devtools/test/chunkgrphconstructTest.m
index d829e77..6a6e50f 100644
--- a/devtools/test/chunkgrphconstructTest.m
+++ b/devtools/test/chunkgrphconstructTest.m
@@ -110,7 +110,7 @@
 
 prefs      = [];
 % prefs.chsmall = 1d-4;
-[cgrph] = chunkgraphinit(verts,edge2verts,fchnks,prefs);
+[cgrph] = chunkgraph(verts,edge2verts,fchnks,prefs);
 
 vstruc = procverts(cgrph);
 rgns = findregions(cgrph);
diff --git a/devtools/test/chunkgrphrcipTest.m b/devtools/test/chunkgrphrcipTest.m
index 154b123..391316a 100644
--- a/devtools/test/chunkgrphrcipTest.m
+++ b/devtools/test/chunkgrphrcipTest.m
@@ -62,7 +62,7 @@
 
 prefs      = [];
 % prefs.chsmall = 1d-4;
-[cgrph] = chunkgraphinit(verts,edge2verts,fchnks,prefs);
+[cgrph] = chunkgraph(verts,edge2verts,fchnks,prefs);
 
 vstruc = procverts(cgrph);
 rgns = findregions(cgrph);
@@ -110,7 +110,7 @@
 
 prefs      = [];
 % prefs.chsmall = 1d-4;
-[cgrph] = chunkgraphinit(verts,edge2verts,fchnks,prefs);
+[cgrph] = chunkgraph(verts,edge2verts,fchnks,prefs);
 
 vstruc = procverts(cgrph);
 rgns = findregions(cgrph);
diff --git a/devtools/test/chunkgrphrcipTransmissionTest.m b/devtools/test/chunkgrphrcipTransmissionTest.m
index 5547c62..7f1c4cf 100644
--- a/devtools/test/chunkgrphrcipTransmissionTest.m
+++ b/devtools/test/chunkgrphrcipTransmissionTest.m
@@ -55,7 +55,7 @@
 
 prefs      = [];
 % prefs.chsmall = 1d-4;
-[cgrph] = chunkgraphinit(verts,edge2verts,fchnks,prefs);
+[cgrph] = chunkgraph(verts,edge2verts,fchnks,prefs);
 
 
 vstruc = procverts(cgrph);
diff --git a/devtools/test/chunkrgrphOpdimTest.m b/devtools/test/chunkrgrphOpdimTest.m
index 4fb2652..38911c5 100644
--- a/devtools/test/chunkrgrphOpdimTest.m
+++ b/devtools/test/chunkrgrphOpdimTest.m
@@ -19,7 +19,7 @@
 fchnks = {}; % this by default gives me straight lines
 
 prefs = struct('maxchunklen',0.5);
-[cgrph] = chunkgraphinit(verts, edge2verts, fchnks, prefs);
+[cgrph] = chunkgraph(verts, edge2verts, fchnks, prefs);
 
 % figure(1); clf; hold on; axis equal; axis off;
 % plot(cgrph);
diff --git a/devtools/test/flamutilitiesTest.m b/devtools/test/flamutilitiesTest.m
index 0698b35..cd8d7bb 100644
--- a/devtools/test/flamutilitiesTest.m
+++ b/devtools/test/flamutilitiesTest.m
@@ -31,7 +31,7 @@
 end
 cparams = [];
 cparams.nover = 2;
-[cgrph] = chunkgraphinit(verts,edge2verts,fchnks,cparams);
+[cgrph] = chunkgraph(verts,edge2verts,fchnks,cparams);
 
 vstruc = procverts(cgrph);
 rgns = findregions(cgrph);