From 283d0dc60fbca73d00895a6ec59454fe9f278bbb Mon Sep 17 00:00:00 2001 From: Travis Askham Date: Tue, 25 Jul 2023 10:14:59 -0400 Subject: [PATCH 1/2] documentation updates to rcip files, working verson of interpolation --- chunkie/+chnk/+rcip/IPinit.m | 6 ++ chunkie/+chnk/+rcip/Rcompchunk.m | 99 ++++++++++++++++++++++++++------ chunkie/+chnk/+rcip/SchurBana.m | 35 ++++++++++- chunkie/+chnk/+rcip/myinv.m | 1 + chunkie/+chnk/+rcip/setup.m | 17 +++++- chunkie/@chunker/merge.m | 27 +++++++-- chunkie/chunkermat.m | 11 ++-- devtools/test/rcipTest.m | 65 ++++++++++++++++----- 8 files changed, 213 insertions(+), 48 deletions(-) diff --git a/chunkie/+chnk/+rcip/IPinit.m b/chunkie/+chnk/+rcip/IPinit.m index c261171..c2828d4 100644 --- a/chunkie/+chnk/+rcip/IPinit.m +++ b/chunkie/+chnk/+rcip/IPinit.m @@ -1,10 +1,16 @@ function [IP,IPW]=IPinit(T,W) + %CHNK.RCIP.IPinit + % % construct the prolongation matrix IP that maps function values % on n_{gl} Gauss-Legendre nodes on [-1,1] to function values at % 2n_{gl} Gauss-Legendre, with shifted and scaled n_{gl} % Gauss-Legendre nodes on each subinterval [-1,0], [0,1], respectively. % % IPW is the weighted prolongation matrix acted on the left side. + % + + % author: Johan Helsing (part of RCIP tutorial) + ngl = length(T); A=ones(ngl); AA=ones(2*ngl,ngl); diff --git a/chunkie/+chnk/+rcip/Rcompchunk.m b/chunkie/+chnk/+rcip/Rcompchunk.m index c58483d..c9b3445 100644 --- a/chunkie/+chnk/+rcip/Rcompchunk.m +++ b/chunkie/+chnk/+rcip/Rcompchunk.m @@ -1,13 +1,11 @@ -function [R]=Rcompchunk(chnkr,iedgechunks,fkern,ndim, ... - Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,... - glxs,sbcmat,lvmat,u,opts) +function [R,rcipsav]=Rcompchunk(chnkr,iedgechunks,fkern,ndim, ... + Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,starL1,circL1,... + sbclmat,sbcrmat,lvmat,rvmat,u,opts) %CHNK.RCIP.Rcompchunk carry out the forward recursion for computing % the preconditioner R where geometry is described as a chunker % % This routine is not intended to be user-callable % -% Adapted from Shidong Jiang's RCIP implementation -% % Function is passed as a handle, number of equations is given by % ndim % @@ -15,32 +13,61 @@ % % Note that matrix must be scaled to have identity on the diagonal, % will not work with scaled version of identity - +% + +% author: Shidong Jiang +% modified: Jeremy Hoskins, Manas Rachh + + k = chnkr.k; dim = chnkr.dim; +nedge = size(iedgechunks,2); +glxs = chnkr.tstor; glws = chnkr.wstor; -if nargin < 14 +% return what's needed to interpolate from coarse + +rcipsav = []; +rcipsav.k = k; +rcipsav.ndim = ndim; +rcipsav.nedge = nedge; +rcipsav.Pbc = Pbc; +rcipsav.PWbc = PWbc; +rcipsav.starL = starL; +rcipsav.starL1 = starL1; +rcipsav.starS = starS; +rcipsav.circL = circL; +rcipsav.circL1 = circL1; +rcipsav.circS = circS; +rcipsav.ilist = ilist; +rcipsav.nsub = nsub; + +if (nargin < 15 || isempty(sbclmat) || isempty(sbcrmat) || ... + isempty(lvmat) || isempty(rvmat) || isempty(u)) [sbclmat,sbcrmat,lvmat,rvmat,u] = chnk.rcip.shiftedlegbasismats(k); end -if nargin < 17 +if nargin < 20 opts = []; end -nedge = size(iedgechunks,2); +savedeep = false; +if isfield(opts,'savedeep') + savedeep = opts.savedeep; +end -ileftright = zeros(nedge,1); -nextchunk = zeros(nedge,1); +rcipsav.savedeep = savedeep; + +if savedeep + rcipsav.R = cell(nsub+1,1); + rcipsav.MAT = cell(nsub,1); + rcipsav.chnkrlocals = cell(nsub,1); +end + + +% grab only those kernels relevant to this vertex -km1 = k-1; -rcs = zeros(km1,dim,nedge); -dcs = zeros(k,dim,nedge); -d2cs = zeros(k,dim,nedge); -dscal = zeros(nedge,1); -d2scal = zeros(nedge,1); -ctr = zeros(dim,nedge); if(size(fkern)==1) fkernlocal = fkern; else @@ -55,6 +82,20 @@ end +rcipsav.fkernlocal = fkernlocal; + +% get coefficients of recentered edge chunks and figure out orientation + +km1 = k-1; +rcs = zeros(km1,dim,nedge); +dcs = zeros(k,dim,nedge); +d2cs = zeros(k,dim,nedge); +dscal = zeros(nedge,1); +d2scal = zeros(nedge,1); +ctr = zeros(dim,nedge); + +ileftright = zeros(nedge,1); +nextchunk = zeros(nedge,1); for i = 1:nedge ic = iedgechunks(1,i); @@ -92,10 +133,19 @@ end end +rcipsav.ctr = ctr; +rcipsav.rcs = rcs; +rcipsav.dcs = dcs; +rcipsav.d2cs = d2cs; +rcipsav.dscal = dscal; +rcipsav.d2scal = d2scal; +rcipsav.ileftright = ileftright; +rcipsav.glxs = glxs; +rcipsav.glws = glws; pref = []; pref.k = k; -pref.nchmax = 5; +pref.nchstor = 5; R = []; @@ -108,6 +158,9 @@ ts = cell(nedge,1); chnkrlocal(1,nedge) = chunker(); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% begin recursion proper + h0=ones(nedge,1); for level=1:nsub h = h0/2^(nsub-level); @@ -175,8 +228,16 @@ if level==1 % Dumb and lazy initializer for R, for now %R=eye(nR); R = inv(MAT(starL,starL)); + if savedeep + rcipsav.R{1} = R; + end end R=chnk.rcip.SchurBana(Pbc,PWbc,MAT,R,starL,circL,starS,circS); + if savedeep + rcipsav.R{level+1} = R; + rcipsav.MAT{level} = MAT(starL,circL); + rcipsav.chnkrlocals{level} = merge(chnkrlocal); + end end end diff --git a/chunkie/+chnk/+rcip/SchurBana.m b/chunkie/+chnk/+rcip/SchurBana.m index 37643b8..8d02b71 100644 --- a/chunkie/+chnk/+rcip/SchurBana.m +++ b/chunkie/+chnk/+rcip/SchurBana.m @@ -1,10 +1,36 @@ function A=SchurBana(P,PW,K,A,starL,circL,starS,circS) - % use matrix block inversion formula to recursively compute the + %CHNK.RCIP.SCHURBANA block inverse used in RCIP + % + % uses a matrix block inversion formula to recursively compute the % preconditioner R. + % + % R_{i} = [PW^T 0] [A^(-1) U]^(-1) [P 0] + % [ 0 I] [V D] [0 I] + % + % 2N x 3N 3N x 3N 3N x 2N + % + % where A = R_{i-1}, U = K(bad,good), V = K(good,bad), D = K(good,good) + % and bad refers to the two close panels on a type b mesh and good the + % larger far panel. % + % This routine uses the Schur-Banachiewicz speedup as suggested by + % Helsing (see e.g. the RCIP Tutorial eq 31): + % + % R_{i} = [M1 M2] + % [M3 M4] + % + % M1 = PW^T A P + PW^T A U (D - VAU)^(-1) VAP + % M2 = -PW^TAU(D-VAU)^(-1) + % M3 = -(D-VAU)^(-1)VAP + % M4 = (D-VAU)^(-1) + % + % this reduces a 3Nx3N matrix inverse to a NxN matrix inverse and + % several mat-mats + % % inputs: % P, PW - nontrivial part (i.e., other than the identity matrices) - % prolongation and weighted prolongation matrices + % prolongation and weighted prolongation matrices from type c to + % type b meshes % K - the system matrix on a type-b mesh along each edge. % the type-b mesh contains three chunks, with two small chunks % close to the corner and one big chunk (twice of the small chunk @@ -34,10 +60,13 @@ % output: % A - R_{i}, the ith iteration of R. % + + % author: Johan Helsing (part of the RCIP tutorial) + VA=K(circL,starL)*A; PTA=PW'*A; PTAU=PTA*K(starL,circL); - DVAUI=chnk.rcip.myinv(K(circL,circL)-VA*K(starL,circL)); + DVAUI=inv(K(circL,circL)-VA*K(starL,circL)); DVAUIVAP=DVAUI*(VA*P); A(starS,starS)=PTA*P+PTAU*DVAUIVAP; A(circS,circS)=DVAUI; diff --git a/chunkie/+chnk/+rcip/myinv.m b/chunkie/+chnk/+rcip/myinv.m index f3fe439..7ebd1d2 100644 --- a/chunkie/+chnk/+rcip/myinv.m +++ b/chunkie/+chnk/+rcip/myinv.m @@ -1,5 +1,6 @@ function Mi=myinv(M) % *** equation (9) in Henderson & Searle *** + np=length(M)/2; A=M(1:np,1:np); U=M(1:np,np+1:2*np); diff --git a/chunkie/+chnk/+rcip/setup.m b/chunkie/+chnk/+rcip/setup.m index 6939737..35d3323 100644 --- a/chunkie/+chnk/+rcip/setup.m +++ b/chunkie/+chnk/+rcip/setup.m @@ -1,5 +1,7 @@ - - function [Pbc,PWbc,starL,circL,starS,circS,ilist] = setup(ngl,ndim,nedge,isstart) +function [Pbc,PWbc,starL,circL,starS,circS,ilist,starL1,circL1] = ... + setup(ngl,ndim,nedge,isstart) + %CHNK.RCIP.SETUP + % % setup for the RCIP forward recursion % inputs: % ngl - number of Gauss-Legendre nodes @@ -15,6 +17,7 @@ % Pbc - prolongation matrix % PWbc - weighted prolongation matrix % starL, circL - bad and good indices for the local system matrix + % starL1, circL1 - bad and good indices for arrays of nodes, normals, etc % starS, circS - bad and good indices for the preconditioner R % [T,W] = lege.exps(ngl); @@ -27,20 +30,30 @@ % circL - good indices for the system matrix M starL = []; circL = []; + starL1 = []; + circL1 = []; indg1 = 2*ngl*ndim + (1:ngl*ndim); indb1 = 1:2*ngl*ndim; + indg11 = 2*ngl+ (1:ngl); + indb11 = 1:2*ngl; indg0 = 1:ngl*ndim; indb0 = ngl*ndim + (1:2*ngl*ndim); + indg01 = 1:ngl; + indb01 = ngl + (1:2*ngl); for i=1:nedge if isstart(i) starL = [starL indb1+3*(i-1)*ngl*ndim]; circL = [circL indg1+3*(i-1)*ngl*ndim]; + starL1 = [starL1 indb11+3*(i-1)*ngl]; + circL1 = [circL1 indg11+3*(i-1)*ngl]; ilist(:,i) = [1, 2]; else starL = [starL indb0+3*(i-1)*ngl*ndim]; circL = [circL indg0+3*(i-1)*ngl*ndim]; + starL1 = [starL1 indb01+3*(i-1)*ngl]; + circL1 = [circL1 indg01+3*(i-1)*ngl]; ilist(:,i) = [2, 3]; end end diff --git a/chunkie/@chunker/merge.m b/chunkie/@chunker/merge.m index a500ae5..b7af2a0 100644 --- a/chunkie/@chunker/merge.m +++ b/chunkie/@chunker/merge.m @@ -1,19 +1,36 @@ -function chnkrout = merge(chnkrs) +function chnkrout = merge(chnkrs,pref) +%MERGE combine array of chunker objects into one chunker +% +% input: +% chnkrs - array of chunker objects, must all have same order chunks +% pref - optional, chunkerpref object +% output: +% chnkrout - chunker object containing all nodes in chunker array. the +% ordering of nodes in chnkrout has the nodes from chnkrs(1) first, +% then chnkrs(2), etc. Adjacency information is updated as +% appropriate to the new indices. chunker data rows are copied as +% well. if chnkrs have different numbers of data rows, then those +% with fewer data rows are padded with zeros on merge. +% if isempty(chnkrs) chnkrout = chunker(); return end assert(isa(chnkrs,'chunker'), 'input must be of chunker type'); +if nargin < 2 + pref = []; +end -chnkrout = chunker(); -%chnkrout = chnkrs(1); +% mandatory setting +pref.k = chnkrs(1).k; +chnkrout = chunker(pref); for i = 1:length(chnkrs) chnkrtemp = chnkrs(i); - assert(chnkrtemp.dim == chnkrout.dim,... - 'chunkers to merge must be in same dimension'); + assert(chnkrtemp.dim == chnkrout.dim && chnkrtemp.k == chnkrout.k,... + 'chunkers to merge must be in same dimension and same order'); nch = chnkrtemp.nch; nchold = chnkrout.nch; istart = nchold+1; diff --git a/chunkie/chunkermat.m b/chunkie/chunkermat.m index fb03517..8b922e2 100644 --- a/chunkie/chunkermat.m +++ b/chunkie/chunkermat.m @@ -435,6 +435,7 @@ if(icgrph && isrcip) + [sbclmat,sbcrmat,lvmat,rvmat,u] = chnk.rcip.shiftedlegbasismats(k); nch_all = horzcat(chnkobj.echnks.nch); [~,nv] = size(chnkobj.verts); ngl = chnkrs(1).k; @@ -465,7 +466,7 @@ if(norm(opdims_test(:)-ndim)>=0.5) fprintf('in chunkermat: rcip: opdims did not match up for for vertex =%d\n',ivert) - fprtinf('returning without doing any rcip correction\m'); + fprintf('returning without doing any rcip correction\n'); break end starind = zeros(1,2*ngl*ndim*nedge); @@ -480,13 +481,13 @@ end - [Pbc,PWbc,starL,circL,starS,circS,ilist] = chnk.rcip.setup(ngl,ndim, ... - nedge,isstart); + [Pbc,PWbc,starL,circL,starS,circS,ilist,starL1,circL1] = ... + chnk.rcip.setup(ngl,ndim,nedge,isstart); % this might need to be fixed in triple junction case R = chnk.rcip.Rcompchunk(chnkrs,iedgechunks,kern,ndim, ... - Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,... - glxs); + Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,starL1,circL1,... + sbclmat,sbcrmat,lvmat,rvmat,u,opts); sysmat_tmp = inv(R) - eye(2*ngl*nedge*ndim); if (~nonsmoothonly) diff --git a/devtools/test/rcipTest.m b/devtools/test/rcipTest.m index c517f37..3a80856 100644 --- a/devtools/test/rcipTest.m +++ b/devtools/test/rcipTest.m @@ -13,10 +13,7 @@ % set wave number zk = 1.1; - - - -nch = 5*ones(1,ncurve); +nch = 8*ones(1,ncurve); a = -1.0; b = 1.0; @@ -70,7 +67,7 @@ fkern = @(s,t) -2*chnk.helm2d.kern(zk,s,t,'D'); -%% +% % sources @@ -86,6 +83,10 @@ ts = 0.0+2*pi*rand(1,nt); targets = [cos(ts);sin(ts)]; targets = 0.2*targets; +targets(:,1) = [0.95;0]; +targets(:,2) = [0,0.36]; +targets(:,3) = [-0.95;0]; +targets(:,2) = [0,0.36]; scatter(sources(1,:),sources(2,:),'o'); scatter(targets(1,:),targets(2,:),'x'); @@ -133,6 +134,8 @@ ncorner = 2; corners= cell(1,ncorner); R = cell(1,ncorner); +rcipsav = cell(1,ncorner); +starinds = cell(1,ncorner); corners{1}.clist = [1,2]; @@ -147,10 +150,11 @@ ndim = 1; -nsub = 100; +nsub = 40; opts = []; + for icorner=1:ncorner clist = corners{icorner}.clist; isstart = corners{icorner}.isstart; @@ -170,30 +174,63 @@ end end + starinds{icorner} = starind; - [Pbc,PWbc,starL,circL,starS,circS,ilist] = chnk.rcip.setup(ngl,ndim, ... + [Pbc,PWbc,starL,circL,starS,circS,ilist,starL1,circL1] = chnk.rcip.setup(ngl,ndim, ... nedge,isstart); opts_use = []; iedgechunks = corners{icorner}.iedgechunks; - tic; R{icorner} = chnk.rcip.Rcompchunk(chnkr,iedgechunks,fkern,ndim, ... - Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,... - glxs); + optsrcip = []; + optsrcip.savedeep = true; + tic; [R{icorner},rcipsav{icorner}] = chnk.rcip.Rcompchunk(chnkr,iedgechunks,fkern,ndim, ... + Pbc,PWbc,nsub,starL,circL,starS,circS,ilist,starL1,circL1,..., + [],[],[],[],[],optsrcip); toc sysmat(starind,starind) = inv(R{icorner}) - eye(2*ngl*nedge*ndim); end sysmat = sysmat + eye(np); -sol = gmres(sysmat,ubdry,np,eps*20,np); +[sol,flag,relres,iter] = gmres(sysmat,ubdry,np,eps*20,np); + +% + +% interpolate to fine grid + +ndepth = 10; +cor = cell(1,ncorner); + +for icorner = 1:ncorner + solhat = sol(starinds{icorner}); + [solhatinterp,srcinfo,wts] = chnk.rcip.rhohatInterp(solhat,rcipsav{icorner},ndepth); -opts.usesmooth=true; + targtemp = targets(:,:) - rcipsav{icorner}.ctr(:,1); + targinfo = []; + targinfo.r = targtemp; + + cmat = fkern(srcinfo,targinfo); + mu = solhatinterp(:).*wts(:); + cor{icorner} = cmat*mu; +end + +% + +soltemp = sol; +for icorner = 1:ncorner + soltemp(starinds{icorner}) = 0; +end + + +opts = []; opts.verb=false; -opts.quadkgparams = {'RelTol',1e-16,'AbsTol',1.0e-16}; -start=tic; Dsol = chunkerkerneval(chnkrtotal,fkern,sol,targets,opts); +start=tic; Dsol = chunkerkerneval(chnkrtotal,fkern,soltemp,targets,opts); t1 = toc(start); fprintf('%5.2e s : time to eval at targs (slow, adaptive routine)\n',t1) +for icorner = 1:ncorner + Dsol = Dsol + cor{icorner}; +end % wchnkr = weights(chnkrtotal); From 3dc989fb7f1a97e1bc7aeef40b534495871ba664 Mon Sep 17 00:00:00 2001 From: Travis Askham Date: Tue, 31 Oct 2023 13:28:34 -0400 Subject: [PATCH 2/2] adding RCIP density interpolation --- chunkie/+chnk/+rcip/rhohatInterp.m | 140 +++++++++++++++++++++++++++++ 1 file changed, 140 insertions(+) create mode 100644 chunkie/+chnk/+rcip/rhohatInterp.m diff --git a/chunkie/+chnk/+rcip/rhohatInterp.m b/chunkie/+chnk/+rcip/rhohatInterp.m new file mode 100644 index 0000000..a573015 --- /dev/null +++ b/chunkie/+chnk/+rcip/rhohatInterp.m @@ -0,0 +1,140 @@ +function [rhohatinterp,srcinfo,wts] = rhohatInterp(rhohat,rcipsav,ndepth) +%CHNK.RCIP.RHOHATINTERP interpolate the (weight-corrected) coarse level +% density rhohat to the requested depth using the backward recursion +% +% When using an RCIP preconditioner (in the style of eq 34 of the RCIP +% tutorial), the resulting coarse level density is +% accurate for interactions separated from the vertex +% *at the coarse scale*. +% By running the recursion for the RCIP preconditioner backwards, an +% equivalent density can be reconstructed on the fine mesh which is +% appropriate for closer interactions (at a distance well-separated from +% the vertex *at the finest level in the mesh*). +% +% - Let Gamma_i be the portion of the boundary in the vicinity of +% the vertex at level i (with level nsub the coarse level portion of the +% boundary and level 1 the lowest level). +% - Let rhohat_i be the weight corrected density on a type c mesh on +% Gamma_i +% - Let the matrix I+K_i be a discretization of the integral operator on a +% type b mesh on Gamma_i +% - Let starL denote the "bad indices" of K_i and circL the "good indices"; +% the bad correspond to the pairs of small close panels on a type b mesh +% - Let starS denote the "bad indices" of R_i and circS the "good indices"; +% the bad correspond to the close panels on a type c mesh +% - Let P be the non-trivial part of an interpolator from type c to type b +% meshes. +% +% Then, to high precision +% +% rhohat_{i-1} = R_{i-1}[ [P 0] R_i^(-1) rhohat_i - +% (I+K_i)(starL,circL)rhohat_i(circS)] +% +% INPUT +% +% rhohat - entries of density rhohat corresponding to RCIP block in matrix. +% These come from the two nearest chunks to the corner and +% should be ordered so that rhohat(starS) and rhohat(circS) +% correspond to the density at the bad and good points (near and +% far chunk) in the correct order (accounting for chunk +% orientation) +% rcipsav - the quantities obtained from a previous call to Rcompchunk +% ndepth - (default: rcip.nsub) compute up to rhohat_{nsub-ndepth+1} +% +% OUTPUT +% +% rhohatinterp - array of interpolated values, +% [rhohat_nsub(circS); rhohat_{nsub-1}(circS); rhohat_{nsub-2}(circS) ... +% rhohat_1(circS); rhohat_1(starS)] +% note that these are mostly the "good" indices except at +% the lowest level (naturally) +% srcinfo - struct, containing points (srcinfo.r), normals (srcinfo.n), etc +% in the same order as rhohatinterp. +% wts - a set of weights for integrating functions sampled at these +% points + + +rhohatinterp = []; +srcinfo = []; +wts = []; + +k = rcipsav.k; +ndim = rcipsav.ndim; +nedge = rcipsav.nedge; +Pbc = rcipsav.Pbc; +PWbc = rcipsav.PWbc; +starL = rcipsav.starL; +starL1 = rcipsav.starL1; +starS = rcipsav.starS; +circL = rcipsav.circL; +circL1 = rcipsav.circL1; +circS = rcipsav.circS; +ilist = rcipsav.ilist; +nsub = rcipsav.nsub; + +if nargin < 3 + ndepth = nsub; +end + +if ndepth > nsub + msg = "depth requested deeper than RCIP recursion performed\n " + ... + "going to depth %d instead"; + warning(msg,nsub); + ndepth = nsub; +end + +savedeep = rcipsav.savedeep; + +if savedeep + + % all relevant quantities are stored, just run backward recursion + + rhohat0 = rhohat; + rhohatinterp = [rhohatinterp; rhohat0(circS)]; + r = []; + d = []; + d2 = []; + n = []; + h = []; + cl = rcipsav.chnkrlocals{nsub}; + wt = weights(cl); + r = [r, cl.rstor(:,circL1)]; + d = [d, cl.dstor(:,circL1)]; + d2 = [d2, cl.d2stor(:,circL1)]; + n = [n, cl.nstor(:,circL1)]; + wts = [wts; wt(circL1(:))]; + + R0 = rcipsav.R{nsub+1}; + for i = 1:ndepth + R1 = rcipsav.R{nsub-i+1}; + MAT = rcipsav.MAT{nsub-i+1}; + rhotemp = R0\rhohat0; + rhohat0 = R1*(Pbc*rhotemp(starS) - MAT*rhohat0(circS)); + if i == ndepth + rhohatinterp = [rhohatinterp; rhohat0]; + r = [r, cl.rstor(:,starL1)]; + d = [d, cl.dstor(:,starL1)]; + d2 = [d2, cl.d2stor(:,starL1)]; + n = [n, cl.nstor(:,starL1)]; + wt = weights(cl); + + wts = [wts; wt(starL1(:))]; + else + cl = rcipsav.chnkrlocals{nsub-i}; + rhohatinterp = [rhohatinterp; rhohat0(circS)]; + r = [r, cl.rstor(:,circL1)]; + d = [d, cl.dstor(:,circL1)]; + d2 = [d2, cl.d2stor(:,circL1)]; + n = [n, cl.nstor(:,circL1)]; + wt = weights(cl); + wts = [wts; wt(circL1(:))]; + end + R0 = R1; + end + + srcinfo.r = r; + srcinfo.d = d; + srcinfo.d2 = d2; + srcinfo.n = n; + +end