fixing more chnkr.h issues and adding test for arclengthfun

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);
 

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
 

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
 %

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 = sum(chnkr.wts,1);
+s = A*ds;
 
 for i = 1:ncomp
     nch = nchs(i);

chunkie/@chunker/move.m

@@ -19,4 +19,4 @@
 
 obj.wts = weights(obj);
 
-end
+end

chunkie/@chunker/nearest.m

@@ -71,7 +71,6 @@
     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.');

chunkie/@chunker/whts.m

@@ -1,12 +1,12 @@
-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(:)).');
+  w = chnkr.wstor;
+  wts = reshape(sqrt(sum((chnkr.d).^2,1)),k,nch);
+  wts = wts.*w(:);
 
 end

chunkie/chunkerkernevalmat.m

@@ -273,13 +273,12 @@
     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, ...
+        mat(:,jmat:jmatend) =  chnk.adapgausswts(r,d,n,d2,ct,bw,i,targs, ...
                     targd,targd2,kern,opdims,t,w,opts);
 
         js1 = jmat:jmatend;

chunkie/chunkermat.m

@@ -573,14 +573,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;

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);