adding data to chunkgraphs and RCIP

chunkie/+chnk/+rcip/Rcompchunk.m

@@ -94,9 +94,16 @@
 d2scal = zeros(nedge,1);
 ctr = zeros(dim,nedge);
 
+if max([chnkr.datadim]) > 0
+datacs = zeros(k,max([chnkr.datadim]),nedge);
+end
+
 ileftright = zeros(nedge,1);
 nextchunk = zeros(nedge,1);
 
+% flag for whether the edges meeting at this vertex carry data
+hasdata = false(nedge,1);
+
 for i = 1:nedge
     ic = iedgechunks(1,i);
     ie = iedgechunks(2,i);
@@ -106,6 +113,10 @@
     d2 = chnkri.d2(:,:,ie);    
     il = chnkri.adj(1,ie);
     ir = chnkri.adj(2,ie);
+    hasdata(i) = chnkri.hasdata;
+    if hasdata(i)
+        datai = chnkri.data(:,:,ie);
+    end
     if (il > 0 && ir < 0)
         nextchunk(i) = il;
         ileftright(i) = 1;
@@ -116,7 +127,10 @@
         dcs(:,:,i) = u*(d.');
         d2cs(:,:,i) = u*(d2.');
         dscal(i) = 2;
-        d2scal(i) = 4; 
+        d2scal(i) = 4;
+        if hasdata(i)
+            datacs(:,1:chnkri.datadim,i) = u*(datai.');
+        end
     elseif (il < 0 && ir > 0)
         nextchunk(i) = ir;
         ileftright(i) = -1;
@@ -126,7 +140,10 @@
         dcs(:,:,i) = u*(d.');
         d2cs(:,:,i) = u*(d2.');
         dscal(i) = 2;
-        d2scal(i) = 4; 
+        d2scal(i) = 4;
+        if hasdata(i)
+            datacs(:,1:chnkri.datadim,i) = u*(datai.');
+        end
     else
         error('RCIP: edge chunk not adjacent to one vertex and one neighbor')
     end
@@ -142,6 +159,10 @@
 rcipsav.glxs = glxs;
 rcipsav.glws = glws;
 
+if max([chnkr.datadim]) > 0
+rcipsav.datacs = datacs;
+end
+
 pref = []; 
 pref.k = k;
 pref.nchstor = 5;
@@ -214,6 +235,10 @@
     for i=1:nedge
         chnkrlocal(i) = chnk.rcip.chunkerfunclocal(@(t) shiftedcurve(t,rcs(:,:,i),dcs(:,:,i), ...
             dscal(i),d2cs(:,:,i),d2scal(i),ileftright(i)),ts{i},pref,glxs,glws);
+        if hasdata(i)
+            chnkrlocal(i) = setdatalocal(chnkrlocal(i),ts{i},datacs(:,:,i), ...
+                ileftright(i),glxs);
+        end
     end
 
     % at the top level, append/prepend the next chunk
@@ -227,6 +252,13 @@
             chnkrlocal(i).d(:,:,nchi+1) = chnkr(ic).d(:,:,nc);                
             chnkrlocal(i).d2(:,:,nchi+1) = chnkr(ic).d2(:,:,nc); 
 
+            if hasdata(i)
+                if ~chnkrlocal(i).hasdata
+                    chnkrlocal(i) = chnkrlocal(i).makedatarows(size(chnkr(ic).data,1));
+                end
+                chnkrlocal(i).data(:,:,nchi+1) = chnkr(ic).data(:,:,nc);
+            end
+
             if ileftright(i) == -1
                 chnkrlocal(i).adj(1,nchi+1) = nchi;
                 chnkrlocal(i).adj(2,nchi+1) = -1;
@@ -291,3 +323,36 @@
 
 end
 
+function chnkr = setdatalocal(chnkr,ts,datac,ilr,xs)
+%SETDATALOCAL set data on a local chunker, mirroring shiftedcurve's
+% Legendre expansion for d/d2 (plain expansion, no vanishing-at-corner
+% factor, since data need not vanish at the vertex like r does)
+
+datadim = size(datac,2);
+nch = chnkr.nch;
+k = chnkr.k;
+nm1 = k-1;
+
+chnkr = chnkr.makedatarows(datadim);
+
+ab = zeros(2,nch);
+ab(1,:) = ts(1:end-1);
+ab(2,:) = ts(2:end);
+
+for i = 1:nch
+    a = ab(1,i);
+    b = ab(2,i);
+    tt = a + (b-a)*(xs+1)/2;
+
+    if (ilr == -1)
+        ttshift = 2*tt-1;
+    else
+        ttshift = 2*tt+1;
+    end
+    pols = lege.pols(ttshift,nm1); pols = pols.';
+
+    chnkr.data(:,:,i) = (pols*datac).';
+end
+
+end
+

chunkie/@chunker/arclengthder.m

@@ -19,7 +19,12 @@
 %
 
 dmat = lege.dermat(chnkr.k);
-ds = squeeze(sqrt(chnkr.d(1,:,:).^2+chnkr.d(2,:,:).^2));
+
+if chnkr.dim == 2
+    ds = squeeze(sqrt(chnkr.d(1,:,:).^2+chnkr.d(2,:,:).^2));
+elseif chnkr.dim == 3 
+    ds = squeeze(sqrt(chnkr.d(1,:,:).^2+chnkr.d(2,:,:).^2+chnkr.d(3,:,:).^2));
+end
 
 du = dmat*reshape(u,[chnkr.k chnkr.nch]);
 du = du ./ ds;

chunkie/@chunker/chunker.m

@@ -70,6 +70,8 @@
 %   [re,taue] = chunkends(obj,ich) - get the endpoints of chunks
 %   flag = flagnear(obj,pts,opts) - flag points near the boundary
 %   obj = obj.move(r0,r1,trotat,scale) - translate, rotate, etc
+%   du = arclengthder(chnkr,u) - returns arclength derivative of u along
+%              the boundary of the chunker
 
 
 % author: Travis Askham (askhamwhat@gmail.com)

chunkie/@chunkgraph/arclengthder.m

@@ -0,0 +1,32 @@
+function du = arclengthder(cgrph,u)
+%ARCLENGTHDER arc length derivative of function on chunkgraph
+%
+% Calculates the arclength derivative of u using Gauss-Legendre 
+%   spectral differentiation matrix. The derivative on the ith chunk is 
+%   given by du(:,i) where du is the output of this routine
+%
+% Syntax: du = arclengthder(cgrph,u)
+%
+% Input:
+%   cgrph - chunkgraph object
+%   u - function or density defined on cgrph
+%
+% Output:
+%   du - arclength derivative of u
+%
+% Examples:
+%   du = arclengthder(cgrph,u);
+%
+
+dmat = lege.dermat(cgrph.k);
+
+if cgrph.dim == 2
+    ds = squeeze(sqrt(cgrph.d(1,:,:).^2+cgrph.d(2,:,:).^2));
+elseif cgrph.dim == 3 
+    ds = squeeze(sqrt(cgrph.d(1,:,:).^2+cgrph.d(2,:,:).^2+cgrph.d(3,:,:).^2));
+end
+
+du = dmat*reshape(u,[cgrph.k cgrph.nch]);
+du = du ./ ds;
+
+end

chunkie/@chunkgraph/chunkgraph.m

@@ -32,7 +32,7 @@
 %   k - integer, number of Legendre nodes on each chunk
 %   dim - integer, dimension of the ambient space in which the curve is 
 %             embedded
-%   In the rest, nch = \sum_{j} chunkgraph.echnks(j).nch;
+%   nch - integer, total number of chunks in the chunkgraph
 %   npt - returns k*nch, the total number of points on the curve
 %   r - dim x k x nch array, r(:,i,j) gives the coordinates of the ith 
 %         node on the jth chunk of the chunkgraph
@@ -80,6 +80,8 @@
 %   rflag = datares(obj,opts) - check if data in data rows is resolved
 %   edge_reg = find_regions_of_edges(obj) - determine regions on either 
 %                side of each edge
+%   du = arclengthder(cgrph,u) - returns arclength derivative of u along
+%              the boundary of the chunkgraph
 %
 % Syntax:
 %
@@ -139,8 +141,13 @@
         k
         dim
         datadim
+        nch
     end
-
+
+    properties(SetAccess=private)
+        nedge
+    end
+
     methods
         function obj = chunkgraph(verts, edgesendverts, fchnks, cparams, pref)
             %CHUNKGRAPH constructor. Documented above.
@@ -173,6 +180,7 @@
             obj.edgesendverts = edgesendverts;
             obj.v2emat        = build_v2emat(obj);
             obj.echnks        = chunker.empty;
+            obj.nedge         = nedge;
 
             if nargin < 3
                 fchnks = [];
@@ -373,6 +381,12 @@
         function obj = set.wts(obj,val)
             obj.wts = val;
         end
+        function obj = set.data(obj,val)
+            nptvec = [obj.echnks.npt];
+            for ii = 1:length(obj.echnks)
+                obj.echnks(ii).data(:,:) = val(:,sum(nptvec(1:ii-1))+1:sum(nptvec(1:ii)));
+            end
+        end
         function r = get.r(obj)
             chnk = merge(obj.echnks);
             r = chnk.r;
@@ -414,6 +428,9 @@
         function dim = get.dim(obj)
             dim = size(obj.r,1);
         end
+        function nch = get.nch(obj)
+            nch = size(obj.r,3);
+        end        
         function datadim = get.datadim(obj)
             datadim = size(obj.data,1);
         end
@@ -461,7 +478,8 @@
         obj = plus(v,obj)
         obj = mtimes(A,obj)        
         obj = rotate(obj,theta,r0,r1)
-        obj = reflect(obj,theta,r0,r1)   
+        obj = reflect(obj,theta,r0,r1)  
+        du = arclengthder(obj,u)        
         err = chunk_fun_error(obj,fval)
     end
 

chunkie/@chunkgraph/makedatarows.m

@@ -14,12 +14,9 @@
 %
 %
     if (nrows > 0)
-        datatemp = obj.data;
-        datadimold = obj.datadim;
-        nch = sum(horzcat(obj.echnks.nch));
-        obj.data = zeros(datadimold+nrows,obj.k,nch);
-        obj.data(1:datadimold,:) = datatemp(:,:);
-        obj.hasdata = true;
+        for ii = 1:length(obj.echnks)
+            obj.echnks(ii) = makedatarows(obj.echnks(ii),nrows);
+        end
     else
         if (nrows < 0)
             warning('attempted to add negative rows, doing nothing');