diff --git a/chunkie/+chnk/+quadba/buildmattd.m b/chunkie/+chnk/+quadba/buildmattd.m index dfe942e..07069d5 100644 --- a/chunkie/+chnk/+quadba/buildmattd.m +++ b/chunkie/+chnk/+quadba/buildmattd.m @@ -11,7 +11,6 @@ adj = chnkr.adj; d = chnkr.d; d2 = chnkr.d2; -h = chnkr.h; [~,~,u] = lege.exps(k); diff --git a/chunkie/+chnk/adapgausswts.m b/chunkie/+chnk/adapgausswts.m index 3fd69be..7afc58b 100644 --- a/chunkie/+chnk/adapgausswts.m +++ b/chunkie/+chnk/adapgausswts.m @@ -8,7 +8,7 @@ % WARNING: this routine currently assumes that the kernel function is % translation invariant to recenter (improves stability). % -% Syntax: [mat,maxrecs,numints,iers] = adapgausswts(r,d,d2,h,ct,bw,j, ... +% Syntax: [mat,maxrecs,numints,iers] = adapgausswts(r,d,d2,ct,bw,j, ... % rt,dt,d2t,kern,opdims,t,w,opts) % % Input: @@ -16,7 +16,6 @@ % d - chnkr derivatives at nodes % n - chnkr normals at nodes % d2 - chnkr 2nd derivatives at nodes -% h - lengths of chunks in parameter space % ct - Legendre nodes at order of chunker % bw - barycentric interpolation weights for Legendre nodes at order of % chunker diff --git a/chunkie/+chnk/chunk_nearparam.m b/chunkie/+chnk/chunk_nearparam.m index efdf2fd..9d670cd 100644 --- a/chunkie/+chnk/chunk_nearparam.m +++ b/chunkie/+chnk/chunk_nearparam.m @@ -24,6 +24,7 @@ maxnewt = 15; thresh0 = 1.0d-14; +iextra = 3; if nargin < 4 || isempty(t) || nargin < 5 || isempty(u) [t,~,u] = lege.exps(chnkr.k); @@ -40,12 +41,6 @@ thresh0 = opts.thresh; end -dist = Inf; - -if nargin < 3 - ich = 1:nch; -end - dim = size(pts,1); npts = numel(pts)/dim; pts = reshape(pts,dim,npts); @@ -92,9 +87,10 @@ rdiff0 = r0(:) - ref; dprime = rdiff0.'*d0(:); - dprime2 = d0(:).'*d0(:) + d0(:).'*d20(:); + dprime2 = (d0(:).'*d0(:) + rdiff0(:).'*d20(:)); newtsuccess = false; + ifend = 0; for iii = 1:maxnewt @@ -102,14 +98,17 @@ deltat = - dprime/dprime2; - t0 = t0+deltat; + t1 = t0 + deltat; - if (t0 > 1.0) - t0 = 1; + if (t1 > 1.0) + t1 = 1; end - if (t0 < -1.0) - t0 = -1; + if (t1 < -1.0) + t1 = -1; end + + deltat = t1-t0; + t0 = t1; all0 = (lege.exev(t0,cfs)); r0 = all0(1:dim); @@ -118,9 +117,13 @@ rdiff0 = r0(:)-ref(:); dprime = rdiff0.'*d0(:); - dprime2 = d0(:).'*d0(:) + d0(:).'*d20(:); + dprime2 = d0(:).'*d0(:) + rdiff0(:).'*d20(:); + + if min(abs(dprime),abs(deltat)) < thresh + ifend = ifend+1; + end - if abs(dprime) < thresh + if (ifend >= iextra) newtsuccess = true; break; end @@ -196,7 +199,7 @@ d20 = d21; rdiff0 = rdiff1; dprime = rdiff0.'*d0(:); - dprime2 = d0(:).'*d0(:) + d0(:).'*d20(:); + dprime2 = d0(:).'*d0(:); dist0 = dist1; lam = lam/3; if abs(dprime) < thresh diff --git a/chunkie/+chnk/funcuni.m b/chunkie/+chnk/funcuni.m index ac93bd3..b11dffb 100644 --- a/chunkie/+chnk/funcuni.m +++ b/chunkie/+chnk/funcuni.m @@ -104,16 +104,18 @@ for i = 1:nch a=ab(1,i); b=ab(2,i); + h = (b-a)/2; 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; + chnkr.d(:,:,i) = reshape(out{2},dim,k)*h; + chnkr.d2(:,:,i) = reshape(out{3},dim,k)*h*h; end chnkr.adj = adjs(:,1:nch); +chnkr.wtsstor(:,1:nch) = weights(chnkr); +chnkr.nstor(:,:,1:nch) = normals(chnkr); end diff --git a/chunkie/@chunker/arclengthdens.m b/chunkie/@chunker/arclengthdens.m index cecc8d2..470e567 100644 --- a/chunkie/@chunker/arclengthdens.m +++ b/chunkie/@chunker/arclengthdens.m @@ -1,10 +1,9 @@ function ds = arclengthdens(chnkr) %ARCLENGTHDENS arc length density on chunks % -% warning: this takes the length of the ith chunk in parameter space to -% be 2*chnkr.h(i) as opposed to 2. Thus the smooth integration rule +% The smooth integration rule % on the ith chunk is -% ds(:,i).*w*chnkr.h(i) +% ds(:,i).*w % where w are the standard Legendre weights of appropriate order and % ds is the output of this routine % diff --git a/chunkie/@chunker/arclengthfun.m b/chunkie/@chunker/arclengthfun.m index 092da6b..18dd033 100644 --- a/chunkie/@chunker/arclengthfun.m +++ b/chunkie/@chunker/arclengthfun.m @@ -13,10 +13,9 @@ istart = 1; A = lege.intmat(chnkr.k); -[~,w] = lege.exps(chnkr.k); ds = arclengthdens(chnkr); -chunklens = sum((w(:)*(chnkr.h(:).')).*ds,1); -s = (A*ds).*(chnkr.h(:).'); +chunklens = chunklen(chnkr); +s = A*ds; for i = 1:ncomp nch = nchs(i); diff --git a/chunkie/@chunker/move.m b/chunkie/@chunker/move.m index 8cbe34a..9484aef 100644 --- a/chunkie/@chunker/move.m +++ b/chunkie/@chunker/move.m @@ -19,4 +19,4 @@ obj.wts = weights(obj); -end \ No newline at end of file +end diff --git a/chunkie/@chunker/nearest.m b/chunkie/@chunker/nearest.m index d585f3e..f87c7d3 100644 --- a/chunkie/@chunker/nearest.m +++ b/chunkie/@chunker/nearest.m @@ -40,6 +40,7 @@ maxnewt = 15; thresh0 = 1.0d-9; iextra = 3; +nch = chnkr.nch; if nargin < 5 || isempty(u) [~,~,u] = lege.exps(chnkr.k); @@ -69,124 +70,15 @@ % grab chunk data ii = ich(i); r = chnkr.rstor(:,:,ii); - d = chnkr.dstor(:,:,ii); - d2 = chnkr.d2stor(:,:,ii); - h = chnkr.hstor(ii); - rc = u*(r.'); - dc = u*(d.'); - d2c = u*(d2.'); - % starting guess on over sampled grid - rover = (ainterpover*(r.')).'; - rdist0 = sqrt(sum((rover-ref(:)).^2,1)); - [disti,ind] = min(rdist0); - tt = xover(ind); - - % check the endpoints - ts = [-1;1]; - rs = (lege.exev(ts,rc)).'; - rsdist = sqrt(sum((rs-ref).^2,1)); - [dists,inds] = min(rsdist); - if dists < disti - % continue if endpoints closer - ri = rs(:,inds); - if dists < dist - dist = dists; - ds = (lege.exev(ts,dc)).'; - d2s = (lege.exev(ts,d2c)).'; - di = ds(:,inds); - d2i = d2s(:,inds); - rn = ri; dn = di; d2n = d2i; - tn = ts(inds); - ichn = ii; - end - continue; - end - - % run levenberg - - t0 = tt; - - ifend = 0; - - r0 = (lege.exev(t0,rc)).'; - - rdiff0 = ref(:) - r0(:); - - drc = lege.derpol(rc); - %d2rc = lege.derpol(drc); - thresh = thresh0; - - for iii = 1:maxnewt - %d0 = (lege.exev(t0,dc)).'*h; - d0 = (lege.exev(t0,drc)).'; % actual derivative of chunk -% d20 = (lege.exev(t0,d2rc)).'; % actual 2nd derivative of chunk - - % newton info -% -% dprime = rdiff0.'*d0(:); -% dprime2 = (d0(:).'*d0(:) + d0(:).'*d20(:)); - - % levenberg instead - - deltat = d0(:)\(rdiff0(:)); - t1 = t0+deltat; - - if (t1 > 1.0) - t1 = (1.0+t0)/2; - end - if (t1 < -1.0) - t1 = (-1.0+t0)/2; - end - - r1 = (lege.exev(t1,rc)).'; - rdiff1 = ref(:)-r1(:); - - err = min(abs(deltat),abs(sum(d0(:).*rdiff0(:)))); - if err < thresh - ifend = ifend+1; - end - - %sum(rdiff0(:).*d0(:)) - %abs(deltat) - - t0 = t1; - rdiff0 = rdiff1; - -% if (err < thresh) -% ifend = ifend + 1; -% end -% - %t0 = t1; - - if (ifend >= iextra) - break; - end - end - - % done levenberg - - ri = (lege.exev(t0,rc)).'; - disti = sqrt(sum((ri(:)-ref(:)).^2,1)); - if (ifend < iextra) - warning('newton failed in nearest'); -% deltat -% disti -% jac -% t0 -% rhs -% ifend -% sum(jac.*rhs) - else -% fprintf('success\n') - end - + [ti,ri,di,d2i,disti] = chnk.chunk_nearparam(r,ref,opts,chnkr.tstor,u); + disti = sqrt(disti); if disti < dist dist = disti; rn = ri; - dn = (lege.exev(t1,dc)).'*h; - d2n = (lege.exev(t1,d2c)).'*h*h; - tn = t1; + dn = di; + d2n = d2i; + tn = ti; ichn = ii; end end diff --git a/chunkie/@chunker/whts.m b/chunkie/@chunker/whts.m index 1da1df0..43bd30c 100644 --- a/chunkie/@chunker/whts.m +++ b/chunkie/@chunker/whts.m @@ -1,12 +1,8 @@ -function wchnk = whts(chnkr) +function wts = whts(chnkr) warning('whts is deprecated and will be removed, use weights instead'); - k = chnkr.k; - nch = chnkr.nch; - [~,w] = lege.exps(k); - wchnk = reshape(sqrt(sum((chnkr.d).^2,1)),k,nch); - wchnk = wchnk.*bsxfun(@times,w(:),(chnkr.h(:)).'); + wts = weights(chnkr); end diff --git a/chunkie/chunkerinterior.m b/chunkie/chunkerinterior.m index 4a23291..f0c29fe 100644 --- a/chunkie/chunkerinterior.m +++ b/chunkie/chunkerinterior.m @@ -99,7 +99,6 @@ chnkr.d2(:,:,istart:iend) = fliplr(chnkr.d2(:,:,1:nch)); chnkr.d2(1,:,istart:iend) = chnkr.d2(1,:,istart:iend); - chnkr.h(istart:iend) = chnkr.h(1:nch); chnkr.wts = weights(chnkr); chnkr.n = normals(chnkr); chnkr.adj(1,:) = 0:chnkr.nch-1; diff --git a/chunkie/chunkerkernevalmat.m b/chunkie/chunkerkernevalmat.m index 08ccc15..7a5f593 100644 --- a/chunkie/chunkerkernevalmat.m +++ b/chunkie/chunkerkernevalmat.m @@ -273,28 +273,13 @@ d = chnkr.d; n = chnkr.n; d2 = chnkr.d2; - h = chnkr.h; for i = 1:nch jmat = 1 + (i-1)*k*opdims(2); jmatend = i*k*opdims(2); - mat(:,jmat:jmatend) = chnk.adapgausswts(r,d,n,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; + mat(:,jmat:jmatend) = chnk.adapgausswts(r,d,n,d2,ct,bw,i,targs, ... + targd,targn,targd2,kern,opdims,t,w,opts); end else diff --git a/chunkie/chunkermat.m b/chunkie/chunkermat.m index 68da513..9aa6265 100644 --- a/chunkie/chunkermat.m +++ b/chunkie/chunkermat.m @@ -556,14 +556,13 @@ d = chnkr.d; n = chnkr.n; d2 = chnkr.d2; -h = chnkr.h; 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,n,d2,h,ct,bw,i,targs(:,ji), ... + mat1 = chnk.adapgausswts(r,d,n,d2,ct,bw,i,targs(:,ji), ... targd(:,ji),targn(:,ji),targd2(:,ji),kernev,opdims,t,w,opts); js1 = jmat:jmatend; diff --git a/devtools/test/arclengthfunTest.m b/devtools/test/arclengthfunTest.m new file mode 100644 index 0000000..1d11d70 --- /dev/null +++ b/devtools/test/arclengthfunTest.m @@ -0,0 +1,36 @@ +%ARCLENGTHFUNTEST tests the arclengthfun routine for computing distance +% along length of curve +% + +clearvars; close all; +seed = 8675309; +rng(seed); +addpaths_loc(); + +% geometry parameters and construction + + +cparams = []; +cparams.eps = 1.0e-9; +pref = []; +pref.k = 16; +narms = 0; +amp = 0.0; +start = tic; chnkr = chunkerfunc(@(t) starfish(t,narms,amp),cparams,pref); + +[s, ~, chnkr] = arclengthfun(chnkr); +% Compare to analytic arclength for circle +ts = squeeze(atan2(chnkr.r(2,:,:), chnkr.r(1,:,:))); +ts(ts<0) = ts(ts<0) + 2*pi; +assert(norm(s-ts) < 1e-12); + +% Now test two circles +chnkrs(2) = chunker(); +chnkrs(1) = chnkr; +rfac = 1.1; +chnkrs(2) = move(chnkr, [0;0], [3;0], 0, rfac); +chnkrtotal = merge(chnkrs); +[s, nchs, ~] = arclengthfun(chnkrtotal); + +assert(norm(s(:,1:nchs(1)) - ts) < 1e-12); +assert(norm(s(:,nchs(1)+1:end) - rfac*ts) < 1e-12); \ No newline at end of file diff --git a/devtools/test/chunker_nearestTest.m b/devtools/test/chunker_nearestTest.m new file mode 100644 index 0000000..54ab6c5 --- /dev/null +++ b/devtools/test/chunker_nearestTest.m @@ -0,0 +1,24 @@ +%chunker.nearest test +% + +iseed = 1234; +rng(iseed); + +% compare closest computed point on circle geometry to exact answer + +circ = chunkerfunc(@(t) chnk.curves.bymode(t,1)); +nt = 1000; +thetas = -pi+2*pi*rand(1,nt); + +% the 0.1 here keeps targets away from coordinate singularity of circle +scal = 0.1+2*rand(1,nt); +targs = [cos(thetas); sin(thetas)].*scal; + +for j = 1:nt + [rn,dn,dist,tn,ichn] = nearest(circ,targs(:,j)); + th1 = atan2(targs(2,j),targs(1,j)); + th2 = atan2(rn(2),rn(1)); + err = min(abs(th1-th2),abs(th1-th2+2*pi)); + err = min(err,abs(th1-th2-2*pi)); + assert(err < 1e-12); +end \ No newline at end of file diff --git a/devtools/test/chunkerinteriorTest.m b/devtools/test/chunkerinteriorTest.m index 2538be6..281d205 100644 --- a/devtools/test/chunkerinteriorTest.m +++ b/devtools/test/chunkerinteriorTest.m @@ -82,5 +82,11 @@ fprintf('%5.2e s : time for chunkerinterior (chunker, with fmm)\n',t5); assert(all(in5(:) == 1)); - - +% test axissym option +chnkr = chunkerfunc(@(t) starfish(t),struct('ta',-pi/2,'tb',pi/2,'ifclosed',0)); +nt = 1000; +ttarg = -pi/2+pi*rand(nt,1); scal = 2*rand(1,nt); +targs = starfish(ttarg).*scal; + +in = chunkerinterior(chnkr,targs,struct('axissym',true)); +assert(all(in(:) == (scal(:) <= 1))); \ No newline at end of file diff --git a/devtools/test/chunkerkernevalmat_greenlapTest.m b/devtools/test/chunkerkernevalmat_greenlapTest.m index ed4913f..4ffa782 100644 --- a/devtools/test/chunkerkernevalmat_greenlapTest.m +++ b/devtools/test/chunkerkernevalmat_greenlapTest.m @@ -6,7 +6,7 @@ clearvars; close all; seed = 8675309; rng(seed); -%addpaths_loc(); +addpaths_loc(); doadap = false; @@ -81,17 +81,26 @@ start=tic; Dmat = chunkerkernevalmat(chnkr,kernd,targets,opts); Du = Dmat*densu; toc(start) +opts=[]; opts.forceadap=true; +start=tic; Dmat = chunkerkernevalmat(chnkr,kernd,targets,opts); +Du2 = Dmat*densu; +toc(start) +opts=[]; start=tic; Smat = chunkerkernevalmat(chnkr,kerns,targets,opts); Sun = Smat*densun; toc(start) utarg2 = Sun-Du; +utarg3 = Sun-Du2; % relerr = norm(utarg-utarg2,'fro')/norm(utarg,'fro'); +assert(relerr < 1e-11); fprintf('relative frobenius error %5.2e\n',relerr); +relerr = norm(utarg-utarg3,'fro')/norm(utarg,'fro'); assert(relerr < 1e-11); +