add matrix version of layer potential evaluator

add some doc to starfish code and uniform chunker

chunkie/chunkerfuncuni.m

@@ -0,0 +1,115 @@
+function chnkr = chunkerfuncuni(fcurve,nch,cparams,pref)
+%CHUNKERFUNC create a chunker object corresponding to a parameterized curve
+%
+% Syntax: chnkr = chunkerfunc(fcurve,cparams,pref)
+%
+% Input: 
+%   fcurve - function handle of the form
+%               [r,d,d2] = fcurve(t)
+%            where r, d, d2 are size [dim,size(t)] arrays describing
+%            position, first derivative, and second derivative of a curve
+%            in dim dimensions parameterized by t.
+%
+% Optional input:
+%   nch - number of chunks to use (16)
+%	cparams - curve parameters structure (defaults)
+%       cparams.ta = left end of t interval (0)
+%       cparams.tb = right end of t interval (2*pi)
+%   pref - chunkerpref object or structure (defaults)
+%       pref.nchmax - maximum number of chunks (10000)
+%       pref.k - number of Legendre nodes on chunks (16)
+%
+% Examples:
+%   chnkr = chunkerfuncuni(@(t) starfish(t)); % chunk up starfish w/ standard
+%                                        % options
+%   pref = []; pref.k = 30; 
+%   cparams = []; cparams.eps = 1e-3;
+%   chnkr = chunkerfunc(@(t) starfish(t),cparams,pref); % change up options
+%   
+% see also CHUNKERPOLY, CHUNKERPREF, CHUNKER
+
+% author: Travis Askham ([email protected])
+%
+
+
+if nargin < 2
+    nch = 16;
+end
+if nargin < 3
+    cparams = [];
+end
+if nargin < 4
+    pref = chunkerpref();
+else
+    pref = chunkerpref(pref);
+end
+
+
+ta = 0.0; tb = 2*pi; ifclosed=true;
+chsmall = Inf; nover = 0;
+eps = 1.0e-6;
+lvlr = 'a'; maxchunklen = Inf; lvlrfac = 2.0;
+
+if isfield(cparams,'ta')
+    ta = cparams.ta;
+end	 
+if isfield(cparams,'tb')
+    tb = cparams.tb;
+end	 
+if isfield(cparams,'ifclosed')
+    ifclosed = cparams.ifclosed;
+end	 
+
+ts = linspace(ta,tb,nch+1);
+
+k = pref.k;
+
+dim = checkcurveparam(fcurve,ta);
+pref.dim = dim;
+nout = 3;
+out = cell(nout,1);
+
+[xs,ws] = lege.exps(k);
+
+%       . . . start chunking
+
+ab = zeros(2,nch);
+adjs = zeros(2,nch);
+ab(1,:) = ts(1:end-1);
+ab(2,:) = ts(2:end);
+
+adjs(1,:) = (1:nch) - 1;
+adjs(2,:) = (1:nch) + 1;
+
+if ifclosed
+    adjs(1,1)=nch;
+    adjs(2,nch)=1;
+else
+    adjs(1,1)=-1;
+    adjs(2,nch)=-1;
+end
+
+
+%       up to here, everything has been done in parameter space, [ta,tb]
+%       . . . finally evaluate the k nodes on each chunk, along with 
+%       derivatives and chunk lengths
+
+chnkr = chunker(pref); % empty chunker
+chnkr = chnkr.addchunk(nch);
+
+
+for i = 1:nch
+    a=ab(1,i);
+    b=ab(2,i);
+
+    ts = a + (b-a)*(xs+1)/2;
+    [out{:}] = fcurve(ts);
+    chnkr.r(:,:,i) = reshape(out{1},dim,k);
+    chnkr.d(:,:,i) = reshape(out{2},dim,k);
+    chnkr.d2(:,:,i) = reshape(out{3},dim,k);
+    chnkr.h(i) = (b-a)/2;
+end
+
+chnkr.adj = adjs(:,1:nch);
+
+end

chunkie/chunkerkernevalmat.m

@@ -0,0 +1,243 @@
+function mat = chunkerkernevalmat(chnkr,kern,targs,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
+%
+% Syntax: mat = chunkerkerneval(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
+%
+% Optional input:
+%   opts - structure for setting various parameters
+%       opts.forcesmooth - if = true, only use the smooth integration rule
+%                           (false)
+%       opts.forceadap - if = true, only use adaptive quadrature (false)
+%    NOTE: only one of forcesmooth or forceadap is allowed. If both 
+%           false, a hybrid algorithm is used, where the smooth rule is 
+%           applied for targets separated by opts.fac*length of chunk from
+%           a given chunk and adaptive integration is used otherwise
+%       opts.fac = the factor times the chunk length used to decide 
+%               between adaptive/smooth rule
+%       opts.eps = tolerance for adaptive integration
+%       opts.nonsmoothonly = boolean (false), if true, only compute the
+%                         entries for which a special quadrature is used
+%                         (e.g. self and neighbor interactoins) and return
+%                         in a sparse array.
+
+%
+% output:
+%   mat - (opdims(1)*nt) x opdims(2)*chnkr.npt matrix mapping values 
+%         of a density on the boundary to values of convolution at target 
+%         points
+%
+% see also CHUNKERKERNEVAL
+
+
+% author: Travis Askham ([email protected])
+
+% determine operator dimensions using first two points
+
+
+srcinfo = []; targinfo = [];
+srcinfo.r = chnkr.r(:,1); srcinfo.d = chnkr.d(:,1); 
+srcinfo.d2 = chnkr.d2(:,1);
+targinfo.r = chnkr.r(:,2); targinfo.d = chnkr.d(:,2); 
+targinfo.d2 = chnkr.d2(:,2);
+
+ftemp = kern(srcinfo,targinfo);
+opdims = size(ftemp);
+
+if nargin < 5
+    opts = [];
+end
+
+forcesmooth = false;
+forceadap = 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,'nonsmoothonly'); nonsmoothonly = opts.nonsmoothonly; end
+if isfield(opts,'fac'); fac = opts.fac; end
+if isfield(opts,'eps'); eps = opts.eps; end
+
+[dim,~] = size(targs);
+
+if (dim ~= 2); warning('only dimension two tested'); end
+
+optssmooth = []; 
+optsadap = []; 
+optsadap.eps = eps;
+
+if forcesmooth
+    mat = chunkerkernevalmat_smooth(chnkr,kern,opdims,targs, ...
+        [],optssmooth);
+    return
+end
+
+if forceadap
+    mat = chunkerkernevalmat_adap(chnkr,kern,opdims, ...
+        targs,[],optsadap);
+    return
+end
+
+% smooth for sufficiently far, adaptive otherwise
+
+optsflag = []; optsflag.fac = fac;
+flag = flagnear(chnkr,targs,optsflag);
+spmat = chunkerkernevalmat_adap(chnkr,kern,opdims, ...
+        targs,flag,optsadap);
+
+if nonsmoothonly
+    mat = spmat;
+    return;
+end
+
+mat = chunkerkernevalmat_smooth(chnkr,kern,opdims,targs, ...
+    flag,opts);
+
+mat = mat + spmat;
+
+
+end
+
+
+
+function mat = chunkerkernevalmat_smooth(chnkr,kern,opdims, ...
+    targs,flag,opts)
+
+if nargin < 6
+    flag = [];
+end
+if nargin < 7
+    opts = [];
+end
+
+k = chnkr.k;
+nch = chnkr.nch;
+
+targinfo = []; targinfo.r = targs;
+srcinfo = []; srcinfo.r = chnkr.r(:,:); 
+srcinfo.d = chnkr.d(:,:); srcinfo.d2 = chnkr.d2(:,:);
+
+mat = kern(srcinfo,targinfo);
+wts = weights(chnkr);
+wts2 = repmat( (wts(:)).', opdims(2), 1);
+wts2 = ( wts2(:) ).';
+mat = mat.*wts2;
+
+if isempty(flag)
+    % nothing to erase
+    return
+else
+    % delete interactions in flag array
+    for i = 1:nch
+        jmat = 1 + (i-1)*k*opdims(2);
+        jmatend = i*k*opdims(2);
+
+        rowkill = find(flag(:,i)); 
+        rowkill = (opdims(1)*(rowkill(:)-1)).' + (1:opdims(1)).';
+        mat(rowkill,jmat:jmatend) = 0;
+    end
+end
+
+end
+
+function mat = chunkerkernevalmat_adap(chnkr,kern,opdims, ...
+    targs,flag,opts)
+
+k = chnkr.k;
+nch = chnkr.nch;
+
+if nargin < 5
+    flag = [];
+end
+if nargin < 6
+    opts = [];
+end
+
+[~,nt] = size(targs);
+
+% using adaptive quadrature
+
+
+if isempty(flag)
+    mat = zeros(opdims(1)*nt,opdims(2)*chnkr.npt);
+
+    [t,w] = lege.exps(2*k+1);
+    ct = lege.exps(k);
+    bw = lege.barywts(k);
+    r = chnkr.r;
+    d = chnkr.d;
+    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);
+
+        mat(:,jmat:jmatend) =  chnk.adapgausswts(r,d,d2,h,ct,bw,i,targs, ...
+                    targd,targd2,kern,opdims,t,w,opts);
+
+        js1 = jmat:jmatend;
+        js1 = repmat( (js1(:)).',1,opdims(1)*numel(ji));
+
+        indji = (ji-1)*opdims(1);
+        indji = repmat( (indji(:)).', opdims(1),1) + ( (1:opdims(1)).');
+        indji = indji(:);
+        indji = repmat(indji,1,opdims(2)*k);
+
+        iend = istart+numel(mat1)-1;
+        is(istart:iend) = indji(:);
+        js(istart:iend) = js1(:);
+        vs(istart:iend) = mat1(:);
+        istart = iend+1;
+    end
+
+else
+    is = zeros(nnz(flag)*opdims(1)*opdims(2)*k,1);
+    js = is;
+    vs = is;
+    istart = 1;
+
+    [t,w] = lege.exps(2*k+1);
+    ct = lege.exps(k);
+    bw = lege.barywts(k);
+    r = chnkr.r;
+    d = chnkr.d;
+    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);
+
+        [ji] = find(flag(:,i));
+        mat1 =  chnk.adapgausswts(r,d,d2,h,ct,bw,i,targs(:,ji), ...
+                    targd(:,ji),targd2(:,ji),kern,opdims,t,w,opts);
+
+        js1 = jmat:jmatend;
+        js1 = repmat( (js1(:)).',opdims(1)*numel(ji),1);
+
+        indji = (ji-1)*opdims(1);
+        indji = repmat( (indji(:)).', opdims(1),1) + ( (1:opdims(1)).');
+        indji = indji(:);
+
+        indji = repmat(indji,1,opdims(2)*k);
+
+        iend = istart+numel(mat1)-1;
+        is(istart:iend) = indji(:);
+        js(istart:iend) = js1(:);
+        vs(istart:iend) = mat1(:);
+        istart = iend+1;
+    end
+    mat = sparse(is,js,vs,opdims(1)*nt,opdims(2)*chnkr.npt);
+
+end
+
+end
+

chunkie/demo/demo_scatter.m

@@ -24,6 +24,7 @@
 cparams.nover = 0;
 cparams.maxchunklen = 4.0/zk; % setting a chunk length helps when the
                               % frequency is known
+
 pref = []; 
 pref.k = 16;
 narms =5;

chunkie/starfish.m

@@ -1,5 +1,19 @@
 
 function [r,d,d2] = starfish(t,varargin)
+%STARFISH
+% return position, first and second derivatives of parameterized starfish
+% domain. 
+%
+% Inputs:
+% t - points (in [0,2pi] to evaluate these quantities
+%
+% Optional inputs:
+% narms - integer, number of arms on starfish (5)
+% amp - float, amplitude of starfish arms relative to radius of length 1
+%               (0.3)
+% ctr - float(2), x,y coordinates of center of starfish ( [0,0] )
+% phi - float, phase shift (0)
+% scale - scaling factor (1.0)
 
 narms = 5;
 amp = 0.3;