From 5775190e0d038e4fa4cb766234686ea2bc13581f Mon Sep 17 00:00:00 2001
From: Travis Askham <askhamwhat@gmail.com>
Date: Thu, 3 Jun 2021 22:06:44 -0400
Subject: [PATCH] add matrix version of layer potential evaluator add some doc
 to starfish code and uniform chunker

---
 chunkie/chunkerfuncuni.m                      | 115 +++++++++
 chunkie/chunkerkernevalmat.m                  | 243 ++++++++++++++++++
 chunkie/demo/demo_scatter.m                   |   1 +
 chunkie/starfish.m                            |  14 +
 .../test/chunkerkernevalmat_greenlapTest.m    |  97 +++++++
 5 files changed, 470 insertions(+)
 create mode 100644 chunkie/chunkerfuncuni.m
 create mode 100644 chunkie/chunkerkernevalmat.m
 create mode 100644 devtools/test/chunkerkernevalmat_greenlapTest.m

diff --git a/chunkie/chunkerfuncuni.m b/chunkie/chunkerfuncuni.m
new file mode 100644
index 0000000..b9786cb
--- /dev/null
+++ b/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 (askhamwhat@gmail.com)
+%
+
+
+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
diff --git a/chunkie/chunkerkernevalmat.m b/chunkie/chunkerkernevalmat.m
new file mode 100644
index 0000000..bf60858
--- /dev/null
+++ b/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 (askhamwhat@gmail.com)
+
+% 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
+
diff --git a/chunkie/demo/demo_scatter.m b/chunkie/demo/demo_scatter.m
index aec9ed7..68f6b58 100644
--- a/chunkie/demo/demo_scatter.m
+++ b/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;
diff --git a/chunkie/starfish.m b/chunkie/starfish.m
index 5e35283..4d812e3 100644
--- a/chunkie/starfish.m
+++ b/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;
diff --git a/devtools/test/chunkerkernevalmat_greenlapTest.m b/devtools/test/chunkerkernevalmat_greenlapTest.m
new file mode 100644
index 0000000..ed4913f
--- /dev/null
+++ b/devtools/test/chunkerkernevalmat_greenlapTest.m
@@ -0,0 +1,97 @@
+%CHUNKERKERNEVAL_GREENLAPTEST test the routines for integrating over 
+% chunks against the Green's ID for Laplace
+%
+% 
+
+clearvars; close all;
+seed = 8675309;
+rng(seed);
+%addpaths_loc();
+
+doadap = false;
+
+% geometry parameters and construction
+
+cparams = [];
+cparams.eps = 1.0e-12;
+pref = []; 
+pref.k = 16;
+narms = 5;
+amp = 0.5;
+start = tic; chnkr = chunkerfunc(@(t) starfish(t,narms,amp),cparams,pref); 
+t1 = toc(start);
+
+% sources
+
+ns = 10;
+ts = 0.0+2*pi*rand(ns,1);
+sources = starfish(ts,narms,amp);
+sources = 3.0*sources;
+strengths = randn(ns,1);
+
+% targets
+
+nt = 100;
+ts = 0.0+2*pi*rand(nt,1);
+targets = starfish(ts,narms,amp);
+targets = targets.*repmat(rand(1,nt),2,1);
+
+% plot geo and sources
+
+xs = chnkr.r(1,:,:); xmin = min(xs(:)); xmax = max(xs(:));
+ys = chnkr.r(2,:,:); ymin = min(ys(:)); ymax = max(ys(:));
+
+hold off
+plot(chnkr)
+hold on
+scatter(sources(1,:),sources(2,:),'o')
+scatter(targets(1,:),targets(2,:),'x')
+axis equal 
+
+%
+
+% 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');
+
+opdims = [1 1];
+
+% eval u and dudn on boundary
+
+srcinfo = []; srcinfo.r = sources; 
+targinfo = []; targinfo.r = chnkr.r(:,:); 
+targinfo.d = chnkr.d(:,:);
+kernmats = kerns(srcinfo,targinfo);
+kernmatsprime = kernsprime(srcinfo,targinfo);
+densu = kernmats*strengths;
+densun = kernmatsprime*strengths;
+
+% eval u at targets
+
+targinfo = []; targinfo.r = targets;
+kernmatstarg = kerns(srcinfo,targinfo);
+utarg = kernmatstarg*strengths;
+
+
+% test green's id
+
+opts=[];
+start=tic; Dmat = chunkerkernevalmat(chnkr,kernd,targets,opts); 
+Du = Dmat*densu;
+toc(start)
+start=tic; Smat = chunkerkernevalmat(chnkr,kerns,targets,opts); 
+Sun = Smat*densun;
+toc(start)
+
+utarg2 = Sun-Du;
+
+%
+
+relerr = norm(utarg-utarg2,'fro')/norm(utarg,'fro');
+
+fprintf('relative frobenius error %5.2e\n',relerr);
+
+assert(relerr < 1e-11);
+