diff --git a/chunkie/+chnk/+quadadap/buildmat.m b/chunkie/+chnk/+quadadap/buildmat.m index 6946598..ec42364 100644 --- a/chunkie/+chnk/+quadadap/buildmat.m +++ b/chunkie/+chnk/+quadadap/buildmat.m @@ -22,7 +22,6 @@ adj = chnkr.adj; d = chnkr.d; d2 = chnkr.d2; -h = chnkr.h; n = chnkr.n; [t,wts,u] = lege.exps(k); @@ -77,7 +76,7 @@ dt = d(:,:,ibefore); d2t = d2(:,:,ibefore); nt = n(:,:,ibefore); - submat = chnk.adapgausswts(r,d,n,d2,h,t,bw,i,rt,dt,nt,d2t, ... + submat = chnk.adapgausswts(r,d,n,d2,t,bw,i,rt,dt,nt,d2t, ... kern,opdims,t2,w2); imat = 1 + (ibefore-1)*k*opdims(1); @@ -90,7 +89,7 @@ dt = d(:,:,iafter); nt = n(:,:,iafter); d2t = d2(:,:,iafter); - submat = chnk.adapgausswts(r,d,n,d2,h,t,bw,i,rt,dt,nt,d2t, ... + submat = chnk.adapgausswts(r,d,n,d2,t,bw,i,rt,dt,nt,d2t, ... kern,opdims,t2,w2); imat = 1 + (iafter-1)*k*opdims(1); @@ -101,7 +100,7 @@ % self - submat = chnk.quadggq.diagbuildmat(r,d,n,d2,h,[],i,kern,opdims,... + submat = chnk.quadggq.diagbuildmat(r,d,n,d2,[],i,kern,opdims,... xs0,wts0,ainterps0kron,ainterps0); imat = 1 + (i-1)*k*opdims(1); @@ -163,7 +162,7 @@ d2t = d2(:,targfix); nt = n(:,targfix); - submat = chnk.adapgausswts(r,d,n,d2,h,t,bw,i,rt,dt,nt,d2t, ... + submat = chnk.adapgausswts(r,d,n,d2,t,bw,i,rt,dt,nt,d2t, ... kern,opdims,t2,w2); imats = bsxfun(@plus,(1:opdims(1)).',opdims(1)*(targfix(:)-1).'); diff --git a/chunkie/+chnk/+quadggq/buildmat.m b/chunkie/+chnk/+quadggq/buildmat.m index 9470819..2e24f04 100644 --- a/chunkie/+chnk/+quadggq/buildmat.m +++ b/chunkie/+chnk/+quadggq/buildmat.m @@ -50,7 +50,6 @@ adj = chnkr.adj; d = chnkr.d; d2 = chnkr.d2; -h = chnkr.h; n = chnkr.n; data = []; @@ -101,7 +100,7 @@ if ~isempty(ilist) && ismember(ibefore,ilist) && ismember(j,ilist) % skip construction if both chunks are in the "bad" chunk list else - submat = chnk.quadggq.nearbuildmat(r,d,n,d2,h,data,ibefore,j, ... + submat = chnk.quadggq.nearbuildmat(r,d,n,d2,data,ibefore,j, ... kern,opdims,xs1,wts1,ainterp1kron,ainterp1); imat = 1 + (ibefore-1)*k*opdims(1); @@ -115,7 +114,7 @@ if ~isempty(ilist) && ismember(iafter,ilist) && ismember(j,ilist) % skip construction if both chunks are in the "bad" chunk list else - submat = chnk.quadggq.nearbuildmat(r,d,n,d2,h,data,iafter,j, ... + submat = chnk.quadggq.nearbuildmat(r,d,n,d2,data,iafter,j, ... kern,opdims,xs1,wts1,ainterp1kron,ainterp1); imat = 1 + (iafter-1)*k*opdims(1); @@ -129,7 +128,7 @@ if ~isempty(ilist) && ismember(j,ilist) % skip construction if the chunk is in the "bad" chunk list else - submat = chnk.quadggq.diagbuildmat(r,d,n,d2,h,data,j,kern,opdims,... + submat = chnk.quadggq.diagbuildmat(r,d,n,d2,data,j,kern,opdims,... xs0,wts0,ainterps0kron,ainterps0); imat = 1 + (j-1)*k*opdims(1); diff --git a/chunkie/+chnk/+quadggq/buildmattd.m b/chunkie/+chnk/+quadggq/buildmattd.m index fd820a6..ee5e01c 100644 --- a/chunkie/+chnk/+quadggq/buildmattd.m +++ b/chunkie/+chnk/+quadggq/buildmattd.m @@ -22,7 +22,6 @@ adj = chnkr.adj; d = chnkr.d; d2 = chnkr.d2; -h = chnkr.h; n = chnkr.n; data = []; @@ -89,7 +88,7 @@ if ~isempty(ilist) && ismember(ibefore,ilist) && ismember(j,ilist) % skip construction if both chunks are in the "bad" chunk list else - submat = chnk.quadggq.nearbuildmat(r,d,n,d2,h,data,ibefore,j, ... + submat = chnk.quadggq.nearbuildmat(r,d,n,d2,data,ibefore,j, ... kern,opdims,xs1,wts1,ainterp1kron,ainterp1); imat = 1 + (ibefore-1)*k*opdims(1); @@ -106,7 +105,7 @@ if ~isempty(ilist) && ismember(iafter,ilist) && ismember(j,ilist) % skip construction if both chunks are in the "bad" chunk list else - submat = chnk.quadggq.nearbuildmat(r,d,n,d2,h,data,iafter,j, ... + submat = chnk.quadggq.nearbuildmat(r,d,n,d2,data,iafter,j, ... kern,opdims,xs1,wts1,ainterp1kron,ainterp1); @@ -124,7 +123,7 @@ if ~isempty(ilist) && ismember(j,ilist) % skip construction if the chunk is in the "bad" chunk list else - submat = chnk.quadggq.diagbuildmat(r,d,n,d2,h,data,j,kern,opdims,... + submat = chnk.quadggq.diagbuildmat(r,d,n,d2,data,j,kern,opdims,... xs0,wts0,ainterps0kron,ainterps0); imat = 1 + (j-1)*k*opdims(1); diff --git a/chunkie/+chnk/+quadggq/diagbuildmat.m b/chunkie/+chnk/+quadggq/diagbuildmat.m index 8c405ef..b8ceabf 100644 --- a/chunkie/+chnk/+quadggq/diagbuildmat.m +++ b/chunkie/+chnk/+quadggq/diagbuildmat.m @@ -1,10 +1,10 @@ -function submat = diagbuildmat(r,d,n,d2,h,data,i,fkern,opdims,... +function submat = diagbuildmat(r,d,n,d2,data,i,fkern,opdims,... xs0,whts0,ainterps0kron,ainterps0) %CHNK.QUADGGQ.DIAGBUILDMAT % grab specific boundary data -rs = r(:,:,i); ds = d(:,:,i); d2s = d2(:,:,i); hs = h(i); +rs = r(:,:,i); ds = d(:,:,i); d2s = d2(:,:,i); ns = n(:,:,i); if(isempty(data)) dd = []; @@ -29,7 +29,7 @@ d2fine{i} = (ainterps0{i}*(d2s.')).'; dfinenrm = sqrt(sum(dfine{i}.^2,1)); nfine{i} = [dfine{i}(2,:); -dfine{i}(1,:)]./dfinenrm; - dsdt{i} = (dfinenrm(:)).*whts0{i}*hs; + dsdt{i} = (dfinenrm(:)).*whts0{i}; end srcinfo = []; diff --git a/chunkie/+chnk/+quadggq/nearbuildmat.m b/chunkie/+chnk/+quadggq/nearbuildmat.m index 1bf815c..08502d1 100644 --- a/chunkie/+chnk/+quadggq/nearbuildmat.m +++ b/chunkie/+chnk/+quadggq/nearbuildmat.m @@ -1,4 +1,4 @@ -function submat = nearbuildmat(r,d,n,d2,h,data,i,j,fkern,opdims,... +function submat = nearbuildmat(r,d,n,d2,data,i,j,fkern,opdims,... xs1,whts1,ainterp1kron,ainterp1) %CHNKR.QUADGGQ.NEARBUILDMAT @@ -6,7 +6,7 @@ rs = r(:,:,j); ds = d(:,:,j); d2s = d2(:,:,j); ns = n(:,:,j); rt = r(:,:,i); dt = d(:,:,i); d2t = d2(:,:,i); nt = n(:,:,i); -hs = h(j); + if(isempty(data)) dd = []; else @@ -64,7 +64,7 @@ dfinenrm = sqrt(sum(dfine.^2,1)); %dfinenrm = dfine(1,:,:); % for complex contour, by SJ 09/30/21 -dsdt = dfinenrm(:).*whts1(:)*hs; +dsdt = dfinenrm(:).*whts1(:); dsdtndim2 = repmat(dsdt(:).',opdims(2),1); dsdtndim2 = dsdtndim2(:); diff --git a/chunkie/+chnk/+quadnative/buildmat.m b/chunkie/+chnk/+quadnative/buildmat.m index 7dbbd01..f4e9da2 100644 --- a/chunkie/+chnk/+quadnative/buildmat.m +++ b/chunkie/+chnk/+quadnative/buildmat.m @@ -21,7 +21,6 @@ r = chnkr.rstor; d = chnkr.dstor; d2 = chnkr.d2stor; -h = chnkr.hstor; n = chnkr.nstor; [dim,k,~] = size(r); @@ -34,10 +33,10 @@ srcinfo = []; srcinfo.r = rs; srcinfo.d = ds; srcinfo.d2 = d2s; srcinfo.n = ns; targinfo = []; targinfo.r = rt; targinfo.d = dt; targinfo.d2 = d2t; targinfo.n = nt; -hs = h(j); ht = h(i); + dsnrms = sqrt(sum(ds.^2,1)); -ws = kron(hs(:),wts(:)); +ws = repmat(wts(:),length(j),1); dsdt = dsnrms(:).*ws; diff --git a/chunkie/+chnk/+rcip/Rcomp.m b/chunkie/+chnk/+rcip/Rcomp.m deleted file mode 100644 index 4b75f6d..0000000 --- a/chunkie/+chnk/+rcip/Rcomp.m +++ /dev/null @@ -1,68 +0,0 @@ -function [R]=Rcomp(ngl,nedge,ndim,Pbc,PWbc,nsub,... - starL,circL,starS,circS,ilist,... - h0,isstart,fcurve,glxs,fkern) - % carry out the forward recursion for computing the preconditioner R - % in the RCIP method - % - % A more general version of Shidong's rcip version - % - % Function is passed as a handle, number of equations is given by - % ndim - % - % Kernel on input takes in arguments (chnkrlocal,ilistl); - % - % Note that matrix must be scaled to have identity on the diagonal, - % will not work with scaled version of identity - - - pref = []; - pref.k = ngl; - % size of the system matrix - nsys = 3*ngl*nedge*ndim; - - % size of the preconditioner R - nR = 2*ngl*nedge*ndim; - - ts = zeros(4,nedge); - chnkrlocal(1,nedge) = chunker(); - - for level=1:nsub - h = h0/2^(nsub-level); - - for i=1:nedge - if isstart(i) - ts(:,i) = [0, 0.5, 1, 2]*h(i); - else - ts(:,i) = -[2, 1, 0.5, 0]*h(i); - end - end - % construct local chunks around the corner - for i=1:nedge - chnkrlocal(i) = chnk.rcip.chunkerfunclocal(fcurve{i},ts(:,i),pref,glxs); - end -% figure -% clf -% plot(chnkrlocal,'r.') -% axis equal -% pause - % construct the system matrix for local chunks - if level == 1 - ilistl = []; - else - ilistl = ilist; - end - - opts = []; - % test for opdims ~= [1,1] - MAT = chunkermat(chnkrlocal,fkern,opts,ilistl); - - % - MAT = eye(nsys) + MAT; - if level==1 % Dumb and lazy initializer for R, for now - %R=eye(nR); - R = inv(MAT(starL,starL)); - end - R=chnk.rcip.SchurBana(Pbc,PWbc,MAT,R,starL,circL,starS,circS); - end - - end diff --git a/chunkie/+chnk/+rcip/Rcompchunk.m b/chunkie/+chnk/+rcip/Rcompchunk.m index e106bb1..477967c 100644 --- a/chunkie/+chnk/+rcip/Rcompchunk.m +++ b/chunkie/+chnk/+rcip/Rcompchunk.m @@ -106,7 +106,6 @@ d2 = chnkri.d2(:,:,ie); il = chnkri.adj(1,ie); ir = chnkri.adj(2,ie); - h = chnkri.h(ie); if (il > 0 && ir < 0) nextchunk(i) = il; ileftright(i) = 1; @@ -116,8 +115,8 @@ rcs(:,:,i) = sbcrmat*(r.'); dcs(:,:,i) = u*(d.'); d2cs(:,:,i) = u*(d2.'); - dscal(i) = h*2; - d2scal(i) = h^2*4; + dscal(i) = 2; + d2scal(i) = 4; elseif (il < 0 && ir > 0) nextchunk(i) = ir; ileftright(i) = -1; @@ -126,8 +125,8 @@ rcs(:,:,i) = sbclmat*(r.'); dcs(:,:,i) = u*(d.'); d2cs(:,:,i) = u*(d2.'); - dscal(i) = h*2; - d2scal(i) = h^2*4; + dscal(i) = 2; + d2scal(i) = 4; else error('RCIP: edge chunk not adjacent to one vertex and one neighbor') end @@ -197,7 +196,7 @@ chnkrlocal(i).r(:,:,nchi+1) = chnkr(ic).r(:,:,nc)-ctr(:,i); chnkrlocal(i).d(:,:,nchi+1) = chnkr(ic).d(:,:,nc); chnkrlocal(i).d2(:,:,nchi+1) = chnkr(ic).d2(:,:,nc); - chnkrlocal(i).h(nchi+1) = chnkr(ic).h(nc); + if ileftright(i) == -1 chnkrlocal(i).adj(1,nchi+1) = nchi; chnkrlocal(i).adj(2,nchi+1) = -1; diff --git a/chunkie/+chnk/+rcip/chunkerfunclocal.m b/chunkie/+chnk/+rcip/chunkerfunclocal.m index 8c21dee..f788463 100644 --- a/chunkie/+chnk/+rcip/chunkerfunclocal.m +++ b/chunkie/+chnk/+rcip/chunkerfunclocal.m @@ -71,10 +71,11 @@ ts = a + (b-a)*(xs+1)/2; [out{:}] = fcurve(ts); + h = (b-a)/2; chnkr.rstor(:,:,i) = reshape(out{1},dim,k); - chnkr.dstor(:,:,i) = reshape(out{2},dim,k); - chnkr.d2stor(:,:,i) = reshape(out{3},dim,k); - chnkr.hstor(i) = (b-a)/2; + chnkr.dstor(:,:,i) = reshape(out{2},dim,k)*h; + chnkr.d2stor(:,:,i) = reshape(out{3},dim,k)*h*h; + end chnkr.adjstor(:,1:nch) = adjs(:,1:nch); diff --git a/chunkie/+chnk/adapgausskerneval.m b/chunkie/+chnk/adapgausskerneval.m index a1b7235..63f53e3 100644 --- a/chunkie/+chnk/adapgausskerneval.m +++ b/chunkie/+chnk/adapgausskerneval.m @@ -1,4 +1,4 @@ -function [fints,maxrecs,numints,iers] = adapgausskerneval(r,d,n,d2,h,ct,bw,j,... +function [fints,maxrecs,numints,iers] = adapgausskerneval(r,d,n,d2,ct,bw,j,... dens,rt,nt,dt,d2t,kern,opdims,t,w,opts) %CHNK.ADAPGAUSSKERNEVAL adaptive integration for interaction of kernel on chunk % at targets @@ -13,7 +13,6 @@ % r - chnkr nodes % d - chnkr derivatives 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 @@ -66,7 +65,6 @@ ds = d(:,:,j); ns = n(:,:,j); d2s = d2(:,:,j); -hs = h(j); jstart = opdims(2)*k*(j-1)+1; jend = opdims(2)*k*j; densj = reshape(dens(jstart:jend),opdims(2),k); @@ -161,8 +159,6 @@ end -fints = fints*hs; - end function val = oneintp(a,b,rs,ds,ns,d2s,densj,ct,bw, ... diff --git a/chunkie/+chnk/adapgausswts.m b/chunkie/+chnk/adapgausswts.m index 6babbd5..3fd69be 100644 --- a/chunkie/+chnk/adapgausswts.m +++ b/chunkie/+chnk/adapgausswts.m @@ -1,4 +1,4 @@ -function [mat,maxrecs,numints,iers] = adapgausswts(r,d,n,d2,h,ct,bw,j,... +function [mat,maxrecs,numints,iers] = adapgausswts(r,d,n,d2,ct,bw,j,... rt,dt,nt,d2t,kern,opdims,t,w,opts) %CHNK.ADAPGAUSSWTS adaptive integration for interaction of kernel on chunk % at targets @@ -68,7 +68,6 @@ ds = d(:,:,j); ns = n(:,:,j); d2s = d2(:,:,j); -hs = h(j); stack = zeros(2,maxdepth); vals = zeros(opdims(1)*opdims(2)*k,maxdepth); @@ -167,8 +166,6 @@ end -mat = mat*hs; - end function val = oneintp(a,b,rs,ds,ns,d2s,ct,bw,rt,dt,nt,d2t,kern,opdims,t,w) diff --git a/chunkie/+chnk/pquadwts.m b/chunkie/+chnk/pquadwts.m index 01a95fe..59bc795 100644 --- a/chunkie/+chnk/pquadwts.m +++ b/chunkie/+chnk/pquadwts.m @@ -1,4 +1,4 @@ -function mat = pquadwts(r,d,n,d2,h,ct,bw,j,... +function mat = pquadwts(r,d,n,d2,ct,bw,j,... rt,dt,nt,d2t,kern,opdims,t,w,opts,intp_ab,intp) %CHNK.INTERPQUADWTS product integration for interaction of kernel on chunk % at targets @@ -14,7 +14,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 @@ -34,11 +33,11 @@ % % Helsing-Ojala (interior/exterior?) -h_i = h(j); + xlohi = intp_ab*(r(1,:,j)'+1i*r(2,:,j)'); % panel end points r_i = intp*(r(1,:,j)'+1i*r(2,:,j)'); % upsampled panel -d_i = h_i*(intp*(d(1,:,j)'+1i*d(2,:,j)')); % r' -d2_i = h_i^2*(intp*(d2(1,:,j)'+1i*d2(2,:,j)')); % r'' +d_i = (intp*(d(1,:,j)'+1i*d(2,:,j)')); % r' +d2_i = (intp*(d2(1,:,j)'+1i*d2(2,:,j)')); % r'' sp = abs(d_i); tang = d_i./sp; % speed, tangent n_i = -1i*tang; % normal cur = -real(conj(d2_i).*n_i)./sp.^2; % curvature diff --git a/chunkie/@chunker/chunker.m b/chunkie/@chunker/chunker.m index ded6cb4..54c58dc 100644 --- a/chunkie/@chunker/chunker.m +++ b/chunkie/@chunker/chunker.m @@ -12,9 +12,6 @@ % 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 chunker -% h - nch array of scaling factors for chunks. the chunk derivatives are -% scaled as if the coordinates in r(:,:,j) are sampled on an -% interval of length 2*h(j). This affects integration formulas. % d - dim x k x nch array, d(:,i,j) gives the time derivative of the % coordinate at the ith node on the jth chunk of the chunker % d2 - dim x k x nch array, d(:,i,j) gives the 2nd time derivative of the @@ -79,7 +76,6 @@ d d2 adj - h n wts data @@ -89,7 +85,6 @@ dstor d2stor adjstor - hstor nstor wtsstor datastor @@ -149,7 +144,6 @@ obj.wtsstor = zeros(k,nchstor); obj.d2stor = zeros(dim,k,nchstor); obj.adjstor = zeros(2,nchstor); - obj.hstor = zeros(nchstor,1); obj.vert = {}; obj.hasdata = false; obj.datastor = []; @@ -168,9 +162,6 @@ function adj = get.adj(obj) adj = obj.adjstor(:,1:obj.nch); end - function h = get.h(obj) - h = obj.hstor(1:obj.nch); - end function n = get.n(obj) n = obj.nstor(:,:,1:obj.nch); end @@ -201,9 +192,6 @@ function obj = set.adj(obj,val) obj.adjstor(:,1:obj.nch) = val; end - function obj = set.h(obj,val) - obj.hstor(1:obj.nch) = val; - end function k = get.k(obj) k = size(obj.rstor,2); end @@ -253,7 +241,6 @@ wtstemp = obj.wts; d2temp = obj.d2; adjtemp = obj.adj; - htemp = obj.h; datatemp = obj.data; obj.rstor = zeros(obj.dim,obj.k,nchstornew); obj.dstor = zeros(obj.dim,obj.k,nchstornew); @@ -261,7 +248,6 @@ obj.wtsstor = zeros(obj.k,nchstornew); obj.d2stor = zeros(obj.dim,obj.k,nchstornew); obj.adjstor = zeros(2,nchstornew); - obj.hstor = zeros(nchstornew,1); obj.datastor = zeros(obj.datadim,obj.k,nchstornew); obj.r = rtemp; obj.d = dtemp; @@ -269,7 +255,6 @@ obj.wts = wtstemp; obj.d2 = d2temp; obj.adj = adjtemp; - obj.h = htemp; obj.data = datatemp; obj.nchstor = nchstornew; end diff --git a/chunkie/@chunker/merge.m b/chunkie/@chunker/merge.m index ddc9587..052e3fb 100644 --- a/chunkie/@chunker/merge.m +++ b/chunkie/@chunker/merge.m @@ -42,7 +42,6 @@ ipos = chnkrtemp.adj > 0; chnkrtemp.adj(ipos) = chnkrtemp.adj(ipos)+nchold; chnkrout.adj(:,istart:iend) = chnkrtemp.adj; - chnkrout.h(istart:iend) = chnkrtemp.h; chnkrout.wts(:,istart:iend) = chnkrtemp.wts; chnkrout.n(:,:,istart:iend) = chnkrtemp.n; end diff --git a/chunkie/@chunker/refine.m b/chunkie/@chunker/refine.m index f1ce045..327ce6f 100644 --- a/chunkie/@chunker/refine.m +++ b/chunkie/@chunker/refine.m @@ -18,10 +18,7 @@ % (10*chnkr.nch) % opts.lvlr = level restriction flag ('a'), if 'a', enforce that no % two adjacent chunks should differ in length by more -% than a factor of approx 2 (see lvlrfac). if 't', -% enforce that the -% length in parameter space of two adjacent chunks -% differs by no more than a factor of approx 2. if 'n' +% than a factor of approx 2 (see lvlrfac). if 'n' % don't enforce (not recommended) % opts.lvlrfac = (2.0) factor for enforcing level restriction % opts.maxchunklen = maximum chunk length (Inf). enforce that no @@ -85,9 +82,6 @@ % splitting type stype = 'a'; -if strcmpi(lvlr,'t') - stype = 't'; -end % chunks told to split @@ -165,7 +159,7 @@ % level restriction -if (strcmpi(lvlr,'a') || strcmpi(lvlr,'t')) +if (strcmpi(lvlr,'a')) for ijk = 1:maxiter_lvlr nchold=chnkr.nch; @@ -175,45 +169,24 @@ i1=chnkr.adj(1,i); i2=chnkr.adj(2,i); - if (strcmpi(lvlr,'t')) - rlself = chnkr.h(i); - rl1 = rlself; - rl2 = rlself; - if (i1 > 0) - rl1 = chnkr.h(i1); - end - if (i2 > 0) - rl2 = chnkr.h(i2); - end - if (i1 < 0) - % meets at vertex (parameter space ref not recommended) - rl1 = min(chnkr.h(vert{-i1})); - end - if (i2 < 0) - rl2 = min(chnkr.h(vert{-i2})); - end - - else - - rlself = chunklens(i); + rlself = chunklens(i); - rl1=rlself; - rl2=rlself; + rl1=rlself; + rl2=rlself; - if (i1 > 0) - rl1 = chunklens(i1); - end - if (i2 > 0) - rl2 = chunklens(i2); - end - if (numel(vert) ~= 0) - if (i1 < 0) - rl1 = min(chunklens(vert{-i1})); - end - if (i2 < 0) - rl2 = min(chunklens(vert{-i2})); - end - end + if (i1 > 0) + rl1 = chunklens(i1); + end + if (i2 > 0) + rl2 = chunklens(i2); + end + if (numel(vert) ~= 0) + if (i1 < 0) + rl1 = min(chunklens(vert{-i1})); + end + if (i2 < 0) + rl2 = min(chunklens(vert{-i2})); + end end % only check if self is larger than either of adjacent blocks, diff --git a/chunkie/@chunker/sort.m b/chunkie/@chunker/sort.m index 6bad4c4..8e9887f 100644 --- a/chunkie/@chunker/sort.m +++ b/chunkie/@chunker/sort.m @@ -32,7 +32,6 @@ chnkr.r = chnkr.r(:,:,inds); chnkr.d = chnkr.d(:,:,inds); chnkr.d2 = chnkr.d2(:,:,inds); -chnkr.h = chnkr.h(inds); chnkr.n = chnkr.n(:,:,inds); chnkr.wts = chnkr.wts(:,inds); diff --git a/chunkie/@chunker/split.m b/chunkie/@chunker/split.m index 3f36283..daa172a 100644 --- a/chunkie/@chunker/split.m +++ b/chunkie/@chunker/split.m @@ -66,7 +66,7 @@ if strcmpi(stype1,'a') % first construct dsdt -dsdt = sqrt(sum(d.^2,1))*chnkr.hstor(ich); +dsdt = sqrt(sum(d.^2,1)); dsdt = dsdt(:); rltot = dot(dsdt,w); @@ -124,8 +124,6 @@ d2_2 = lege.exev(ts2,cd2); d2_2 = d2_2.'; -hold=chnkr.hstor(ich); - % update chnkr %i1=chnkr.adjstor(1,ich); @@ -133,19 +131,17 @@ chnkr = chnkr.addchunk(); -chnkr.hstor(ich) = hold*(t1+1)/2; -chnkr.hstor(nch+1) = hold*(1-t1)/2; +h1 = (t1+1)/2; +h2 = (1-t1)/2; chnkr.rstor(:,:,ich) = r_1; chnkr.rstor(:,:,nch+1) = r_2; -chnkr.dstor(:,:,ich) = d_1; -chnkr.wtsstor(:,ich) = (sqrt(sum(d_1.^2,1)).') .* chnkr.wstor * ... - chnkr.hstor(ich); -chnkr.dstor(:,:,nch+1) = d_2; -chnkr.wtsstor(:,nch+1) = (sqrt(sum(d_2.^2,1)).') .* chnkr.wstor * ... - chnkr.hstor(nch+1); -chnkr.d2stor(:,:,ich) = d2_1; -chnkr.d2stor(:,:,nch+1) = d2_2; +chnkr.dstor(:,:,ich) = d_1*h1; +chnkr.wtsstor(:,ich) = (sqrt(sum(d_1.^2,1)).') .* chnkr.wstor; +chnkr.dstor(:,:,nch+1) = d_2*h2; +chnkr.wtsstor(:,nch+1) = (sqrt(sum(d_2.^2,1)).') .* chnkr.wstor; +chnkr.d2stor(:,:,ich) = d2_1*h1*h1; +chnkr.d2stor(:,:,nch+1) = d2_2*h2*h2; chnkr.adjstor(2,ich)=nch+1; chnkr.adjstor(1,nch+1)=ich; diff --git a/chunkie/@chunker/upsample.m b/chunkie/@chunker/upsample.m index 5aeb177..25d6f77 100644 --- a/chunkie/@chunker/upsample.m +++ b/chunkie/@chunker/upsample.m @@ -36,7 +36,6 @@ chnkrup.d2= permute(reshape(upmat*reshape( ... permute(chnkr.d2,[2,1,3]),k,nn),kup,dim,nch),[2,1,3]); chnkrup.adj= chnkr.adj; -chnkrup.h= chnkr.h; chnkrup.n = normals(chnkrup); chnkrup.wts = weights(chnkrup); diff --git a/chunkie/@chunker/weights.m b/chunkie/@chunker/weights.m index af6afe9..fcce774 100644 --- a/chunkie/@chunker/weights.m +++ b/chunkie/@chunker/weights.m @@ -24,6 +24,6 @@ nch = chnkr.nch; w = chnkr.wstor; wts = reshape(sqrt(sum((chnkr.d).^2,1)),k,nch); - wts = wts.*bsxfun(@times,w(:),(chnkr.h(:)).'); + wts = wts.*w(:); end diff --git a/chunkie/@chunkgraph/balance.m b/chunkie/@chunkgraph/balance.m index 01e5be6..111e7f2 100644 --- a/chunkie/@chunkgraph/balance.m +++ b/chunkie/@chunkgraph/balance.m @@ -48,12 +48,11 @@ pinds(ii) = size(obj.echnks(vedge(ii)).r,3); end - h = obj.echnks(vedge(ii)).h(pinds(ii)); k = obj.echnks(vedge(ii)).k; wleg = echnks(vedge(ii)).wstor; - arc = sum(sqrt(ds(1,:).^2+ds(2,:).^2).*wleg'*h); + arc = sum(sqrt(ds(1,:).^2+ds(2,:).^2).*wleg'); parcl(ii) = arc; end diff --git a/chunkie/chunkerfunc.m b/chunkie/chunkerfunc.m index 4367785..7bdc048 100644 --- a/chunkie/chunkerfunc.m +++ b/chunkie/chunkerfunc.m @@ -520,10 +520,10 @@ for j = nout+1:3 out{j} = out{j-1}*dermat*(2/(b-a)); end + h = (b-a)/2; chnkr.rstor(:,:,i) = reshape(out{1},dim,k); - chnkr.dstor(:,:,i) = reshape(out{2},dim,k); - chnkr.d2stor(:,:,i) = reshape(out{3},dim,k); - chnkr.hstor(i) = (b-a)/2; + chnkr.dstor(:,:,i) = reshape(out{2},dim,k)*h; + chnkr.d2stor(:,:,i) = reshape(out{3},dim,k)*h*h; end chnkr.adjstor(:,1:nch) = adjs(:,1:nch); diff --git a/chunkie/chunkerfuncuni.m b/chunkie/chunkerfuncuni.m index b58765d..edfd459 100644 --- a/chunkie/chunkerfuncuni.m +++ b/chunkie/chunkerfuncuni.m @@ -123,10 +123,10 @@ for j = nout+1:3 out{j} = out{j-1}*dermat*(2/(b-a)); end + h = (b-a)/2; chnkr.rstor(:,:,i) = reshape(out{1},dim,k); - chnkr.dstor(:,:,i) = reshape(out{2},dim,k); - chnkr.d2stor(:,:,i) = reshape(out{3},dim,k); - chnkr.hstor(i) = (b-a)/2; + chnkr.dstor(:,:,i) = reshape(out{2},dim,k)*h; + chnkr.d2stor(:,:,i) = reshape(out{3},dim,k)*h*h; end chnkr.adjstor(:,1:nch) = adjs(:,1:nch); diff --git a/chunkie/chunkerintegral.m b/chunkie/chunkerintegral.m index 1a239fe..ed723ef 100644 --- a/chunkie/chunkerintegral.m +++ b/chunkie/chunkerintegral.m @@ -66,7 +66,7 @@ for i = 1:nch rci = rc(:,:,i); dci = dc(:,:,i); - fint = fint + chnk.intchunk.fhandle(f,rci,dci)*chnkr.h(i); + fint = fint + chnk.intchunk.fhandle(f,rci,dci); end else fint = 0.0; @@ -74,7 +74,7 @@ for i = 1:nch dci = dc(:,:,i); fci = fc(:,i); - fint = fint + chnk.intchunk.fcoefs(fci,dci)*chnkr.h(i); + fint = fint + chnk.intchunk.fcoefs(fci,dci); end end diff --git a/chunkie/chunkerkerneval.m b/chunkie/chunkerkerneval.m index fb10f27..5f967ce 100644 --- a/chunkie/chunkerkerneval.m +++ b/chunkie/chunkerkerneval.m @@ -179,7 +179,6 @@ d = chnkr.d; n = chnkr.n; d2 = chnkr.d2; -h = chnkr.h; % interpolation matrix intp = lege.matrin(k,t); % interpolation from k to 2*k @@ -207,7 +206,7 @@ % Helsing-Ojala (interior/exterior?) - mat1 = chnk.pquadwts(r,d,n,d2,h,ct,bw,j,targs(:,ji), ... + mat1 = chnk.pquadwts(r,d,n,d2,ct,bw,j,targs(:,ji), ... targd(:,ji),targn(:,ji),targd2(:,ji),kern,opdims,t,w,opts,intp_ab,intp); % depends on kern, different mat1? fints(ji) = fints(ji) + mat1*dens(idxjmat); @@ -282,7 +281,7 @@ for i = 1:nch densvals = dens(:,:,i); densvals = densvals(:); dsdtdt = sqrt(sum(abs(chnkr.d(:,:,i)).^2,1)); - dsdtdt = dsdtdt(:).*w(:)*chnkr.h(i); + dsdtdt = dsdtdt(:).*w(:); dsdtdt = repmat( (dsdtdt(:)).',opdims(2),1); densvals = densvals.*(dsdtdt(:)); srcinfo = []; srcinfo.r = chnkr.r(:,:,i); @@ -296,7 +295,7 @@ for i = 1:nch densvals = dens(:,:,i); densvals = densvals(:); dsdtdt = sqrt(sum(abs(chnkr.d(:,:,i)).^2,1)); - dsdtdt = dsdtdt(:).*w(:)*chnkr.h(i); + dsdtdt = dsdtdt(:).*w(:); dsdtdt = repmat( (dsdtdt(:)).',opdims(2),1); densvals = densvals.*(dsdtdt(:)); srcinfo = []; srcinfo.r = chnkr.r(:,:,i); @@ -366,7 +365,7 @@ for i = 1:nch densvals = dens(:,:,i); densvals = densvals(:); dsdtdt = sqrt(sum(abs(chnkr.d(:,:,i)).^2,1)); - dsdtdt = dsdtdt(:).*w(:)*chnkr.h(i); + dsdtdt = dsdtdt(:).*w(:); dsdtdt = repmat( (dsdtdt(:)).',opdims(2),1); densvals = densvals.*(dsdtdt(:)); srcinfo = []; srcinfo.r = chnkr.r(:,:,i); @@ -452,14 +451,12 @@ d = chnkr.d; n = chnkr.n; d2 = chnkr.d2; -h = chnkr.h; - if isempty(flag) % do all to all adaptive for i = 1:nch for j = 1:nt - fints1 = chnk.adapgausskerneval(r,d,n,d2,h,ct,bw,i,dens,targs(:,j), ... + fints1 = chnk.adapgausskerneval(r,d,n,d2,ct,bw,i,dens,targs(:,j), ... targd(:,j),targn(:,j),targd2(:,j),kerneval,opdims,t,w,opts); indj = (j-1)*opdims(1); @@ -470,7 +467,7 @@ else % do only those flagged for i = 1:nch [ji] = find(flag(:,i)); - [fints1,maxrec,numint,iers] = chnk.adapgausskerneval(r,d,n,d2,h,ct,bw,i,dens,targs(:,ji), ... + [fints1,maxrec,numint,iers] = chnk.adapgausskerneval(r,d,n,d2,ct,bw,i,dens,targs(:,ji), ... targd(:,ji),targn(:,ji),targd2(:,ji),kerneval,opdims,t,w,opts); indji = (ji-1)*opdims(1); diff --git a/chunkie/chunkerkernevalmat.m b/chunkie/chunkerkernevalmat.m index a697ac2..4169d53 100644 --- a/chunkie/chunkerkernevalmat.m +++ b/chunkie/chunkerkernevalmat.m @@ -310,13 +310,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); [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),kern,opdims,t,w,opts); js1 = jmat:jmatend; @@ -388,7 +387,6 @@ d = chnkr.d; n = chnkr.n; d2 = chnkr.d2; - h = chnkr.h; % interpolation matrix intp = lege.matrin(k,t); % interpolation from k to 2*k @@ -401,7 +399,7 @@ [ji] = find(flag(:,i)); % Helsing-Ojala (interior/exterior?) - mat1 = chnk.pquadwts(r,d,n,d2,h,ct,bw,i,targs(:,ji), ... + mat1 = chnk.pquadwts(r,d,n,d2,ct,bw,i,targs(:,ji), ... targd(:,ji),targn(:,ji),targd2(:,ji),kern,opdims,t,w,opts,intp_ab,intp); % depends on kern, different mat1? js1 = jmat:jmatend; diff --git a/chunkie/chunkerpoly.m b/chunkie/chunkerpoly.m index 69f8199..016d760 100644 --- a/chunkie/chunkerpoly.m +++ b/chunkie/chunkerpoly.m @@ -194,15 +194,15 @@ val1 = edgevals(:,i); chnkr.data(:,:,nch) = bsxfun(@times,val1,onek(:).'); end + h = (l-w2-w1)/2.0; chnkr.r(:,:,nch) = r1 + bsxfun(@times,v(:),(ts(:)).'); - chnkr.d(:,:,nch) = repmat(v(:),1,k); - chnkr.d2(:,:,nch) = zeros(dim,k); + chnkr.d(:,:,nch) = repmat(v(:),1,k)*h; + chnkr.d2(:,:,nch) = zeros(dim,k)*h*h; chnkr.adj(1,nch) = 0; chnkr.adj(2,nch) = 0; if nch > 1 chnkr.adj(1,nch) = nch-1; chnkr.adj(2,nch-1) = nch; end - chnkr.h(nch) = (l-w2-w1)/2.0; if or(i < nv-1,ifclosed) if rounded @@ -245,11 +245,9 @@ chnkr.adj(:,nch+1:nch+ncht) = chnkrt.adj+nch; chnkr.adj(2,nch) = nch+1; chnkr.adj(1,nch+1) = nch; - chnkr.h(nch+1:nch+ncht) = chnkrt.h; if nvals > 0 - ds = bsxfun(@times,reshape(sum((rotmat*chnkrt.d(:,:)).^2,1),k,ncht), ... - (chnkrt.h(:)).'); + ds = reshape(sum((rotmat*chnkrt.d(:,:)).^2,1),k,ncht); dsw = bsxfun(@times,w(:),ds); dssums = sum(dsw,1); dssums2 = cumsum([0,dssums(1:end-1)]); @@ -286,10 +284,11 @@ ncht = depth + 1; chnkrt = chunker(pref); chnkrt = chnkrt.addchunk(ncht); + hall = hadap*w2/2; hall = repmat(hall(:).',k,1); + hall=hall(:).'; chnkrt.r = bsxfun(@times,v,tadap)*w2 + r2; - chnkrt.d = repmat(v,1,k,ncht); - chnkrt.d2 = zeros(dim,k,ncht); - chnkrt.h = hadap*w2/2; + chnkrt.d(:,:) = repmat(v,1,k*ncht).*hall; + chnkrt.d2(:,:) = zeros(dim,k*ncht).*(hall.^2); chnkrt.adj = [0, 1:(ncht-1); 2:ncht, 0]; chnkrt = reverse(chnkrt); chnkrt = sort(chnkrt); @@ -303,7 +302,6 @@ chnkr.adj(2,nch) = nch+1; chnkr.adj(1,nch+1) = nch; chnkr.adj(2,nch+ncht) = 0; - chnkr.h(nch+1:nch+ncht) = chnkrt.h; if nvals > 0 chnkr.data(:,:,nch+1:nch+ncht) = repmat(val1,1,k,ncht); @@ -317,9 +315,9 @@ chnkrt = chunker(pref); chnkrt = chnkrt.addchunk(ncht); chnkrt.r = bsxfun(@times,v2,tadap)*w2 + r2; - chnkrt.d = repmat(v2,1,k,ncht); - chnkrt.d2 = zeros(dim,k,ncht); - chnkrt.h = hadap*w2/2; + chnkrt.d(:,:) = repmat(v2,1,k*ncht).*hall; + chnkrt.d2(:,:) = zeros(dim,k*ncht).*(hall.^2); + chnkrt.adj = [0, 1:(ncht-1); 2:ncht, 0]; chnkrt = sort(chnkrt); @@ -332,7 +330,6 @@ chnkr.adj(:,nch+1:nch+ncht) = chnkrt.adj+nch; chnkr.adj(2,nch) = nch+1; chnkr.adj(1,nch+1) = nch; - chnkr.h(nch+1:nch+ncht) = chnkrt.h; chnkr = chnkr.addvert([nch,nch+1]); diff --git a/devtools/test/adapgausswtsTest.m b/devtools/test/adapgausswtsTest.m index 54471c9..02d565b 100644 --- a/devtools/test/adapgausswtsTest.m +++ b/devtools/test/adapgausswtsTest.m @@ -78,7 +78,7 @@ opdims = zeros(2,1); opdims(1) = 1; opdims(2) = 1; -r = chnkr.r; d = chnkr.d; h = chnkr.h; d2 = chnkr.d2; n =chnkr.n; +r = chnkr.r; d = chnkr.d; d2 = chnkr.d2; n =chnkr.n; k = chnkr.k; [t,w] = lege.exps(k); bw = lege.barywts(k); @@ -98,7 +98,7 @@ ntimes = 100; start = tic; for i = 1:ntimes - [mat,maxrecs,numints,iers] = chnk.adapgausswts(r,d,n,d2,h,t,bw,j, ... + [mat,maxrecs,numints,iers] = chnk.adapgausswts(r,d,n,d2,t,bw,j, ... rt,dt,nt,d2t,fkern,opdims,t2,w2,opts); end t1 = toc(start); diff --git a/devtools/test/chunkermat_truepolygonTest.m b/devtools/test/chunkermat_truepolygonTest.m index e023ff7..e12be30 100644 --- a/devtools/test/chunkermat_truepolygonTest.m +++ b/devtools/test/chunkermat_truepolygonTest.m @@ -119,9 +119,7 @@ % evaluate at targets and compare -opts.usesmooth=false; opts.verb=false; -opts.quadkgparams = {'RelTol',1e-16,'AbsTol',1.0e-16}; start=tic; Ssol = chunkerkerneval(chnkr,kerns,sol,targets,opts); t1 = toc(start); fprintf('%5.2e s : time to eval at targs (slow, adaptive routine)\n',t1) diff --git a/devtools/test/rcip_analyticTest.m b/devtools/test/rcip_analyticTest.m deleted file mode 100644 index 4ee6b19..0000000 --- a/devtools/test/rcip_analyticTest.m +++ /dev/null @@ -1,340 +0,0 @@ - -%RCIPTEST -% -% This file tests the rcip routines for solving the exterior dirichlet -% problem on a domain defined by two arcs of circles meeting at two vertices - -clearvars; close all; -addpaths_loc(); - -ncurve = 2; -chnkr(1,ncurve) = chunker(); - -% set wave number -zk = 1.1; - - - - -nch = 5*ones(1,ncurve); - -a = -1.0; -b = 1.0; - -% set angles for top and bottom curve -theta = [pi/2.2 pi/2.4]; - -% set parametrization start and end points for chunkie objects -tab = [-theta(1)/2 -theta(2)/2; theta(1)/2 theta(2)/2]; - -% set parameters for curve parametrization -vert = [a b; 0 0]; -cpars = cell(1,ncurve); -cpars{1}.v0 = vert(:,1); cpars{1}.v1 = vert(:,2); -cpars{1}.theta = theta(1); cpars{1}.ifconvex = 0; cpars{1}.islocal = -1; - -cpars{2}.v0 = vert(:,1); cpars{2}.v1 = vert(:,2); -cpars{2}.theta = theta(2); cpars{2}.ifconvex = 2; cpars{2}.islocal = -1; - -% number of gauss legendre nodes -ngl = 16; -pref = []; -pref.k = ngl; - -cparams = cell(1,ncurve); -for i=1:ncurve - cparams{i}.ta = tab(1,i); - cparams{i}.tb = tab(2,i); - cparams{i}.ifclosed = false; -end - -inds = [0, cumsum(nch)]; - -[glxs,glwts,glu,glv] = lege.exps(ngl); - -fcurve = cell(1,ncurve); -figure(1) -clf; hold on; -for icurve = 1:ncurve - fcurve{icurve} = @(t) circulararc(t,cpars{icurve}); - chnkr(icurve) = chunkerfuncuni(fcurve{icurve},nch(icurve),cparams{icurve},pref); - plot(chnkr(icurve)); quiver(chnkr(icurve)); -end - -iedgechunks = [1 2; 1 chnkr(2).nch]; -isstart = [1 0]; -nedge = 2; -ndim=1; -[Pbc,PWbc,starL,circL,starS,circS,ilist] = chnk.rcip.setup(ngl,ndim, ... - nedge,isstart); - -fkern = @(s,t) -2*chnk.helm2d.kern(zk,s,t,'D'); - -%% - -% sources - -ns = 10; -ts = 0.0+2*pi*rand(1,ns); -sources = [cos(ts); sin(ts)]; -sources = 3.0*sources; -strengths = randn(ns,1); - -% targets - -nt = 20; -ts = 0.0+2*pi*rand(1,nt); -targets = [cos(ts);sin(ts)]; -targets = 0.2*targets; - -scatter(sources(1,:),sources(2,:),'o'); -scatter(targets(1,:),targets(2,:),'x'); -axis equal - - - - -chnkrtotal = merge(chnkr); -fkern = @(s,t) -2*chnk.helm2d.kern(zk,s,t,'D'); -np = chnkrtotal.k*chnkrtotal.nch; -start = tic; D = chunkermat(chnkrtotal,fkern); -t1 = toc(start); - - - -kerns = @(s,t) chnk.helm2d.kern(zk,s,t,'s'); - -% eval u on bdry - -targs = chnkrtotal.r; targs = reshape(targs,2,chnkrtotal.k*chnkrtotal.nch); -targstau = tangents(chnkrtotal); -targstau = reshape(targstau,2,chnkrtotal.k*chnkrtotal.nch); - -srcinfo = []; srcinfo.r = sources; -targinfo = []; targinfo.r = targs; -kernmats = kerns(srcinfo,targinfo); -ubdry = kernmats*strengths; - - -% eval u at targets - -targinfo = []; targinfo.r = targets; -kernmatstarg = kerns(srcinfo,targinfo); -utarg = kernmatstarg*strengths; - - - -fprintf('%5.2e s : time to assemble matrix\n',t1) - - -sysmat = D; - -RG = speye(np); -ncorner = 2; -corners= cell(1,ncorner); -R = cell(1,ncorner); - - -corners{1}.clist = [1,2]; -corners{1}.isstart = [1,0]; -corners{1}.nedge = 2; - - -corners{2}.clist = [1,2]; -corners{2}.isstart = [0 1]; -corners{2}.nedge = 2; - - -ndim = 1; - -nsub = 100; - -opts = []; - -for icorner=1:ncorner - clist = corners{icorner}.clist; - isstart = corners{icorner}.isstart; - nedge = corners{icorner}.nedge; - cparslocal = cell(1,nedge); - fcurvelocal = cell(1,nedge); - h0 = zeros(1,nedge); - starind = []; - for i=1:nedge - cparslocal{i} = cpars{clist(i)}; - cparslocal{i}.islocal = isstart(i); - fcurvelocal{i} = @(t) circulararc(t,cparslocal{i}); - if(isstart(i)) - h0(i) = chnkr(clist(i)).h(1); - starind = [starind inds(clist(i))*ngl*ndim+(1:2*ngl*ndim)]; - else - h0(i) = chnkr(clist(i)).h(end); - starind = [starind inds(clist(i)+1)*ngl*ndim-fliplr(0:2*ngl*ndim-1)]; - end - end - h0 = h0*2; - - [Pbc,PWbc,starL,circL,starS,circS,ilist] = chnk.rcip.setup(ngl,ndim, ... - nedge,isstart); - opts_use = []; - tic; R{icorner} = chnk.rcip.Rcomp(ngl,nedge,ndim,Pbc, ... - PWbc,nsub,starL,circL, ... - starS,circS,ilist,h0,isstart,fcurvelocal,glxs,fkern); - 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); - -opts.usesmooth=true; -opts.verb=false; -opts.quadkgparams = {'RelTol',1e-16,'AbsTol',1.0e-16}; -start=tic; Dsol = chunkerkerneval(chnkrtotal,fkern,sol,targets,opts); -t1 = toc(start); -fprintf('%5.2e s : time to eval at targs (slow, adaptive routine)\n',t1) - -% - -wchnkr = chnkrtotal.wts; - -relerr = norm(utarg-Dsol,'fro')/(sqrt(chnkrtotal.nch)*norm(utarg,'fro')); -relerr2 = norm(utarg-Dsol,'inf')/dot(abs(sol(:)),wchnkr(:)); -fprintf('relative frobenius error %5.2e\n',relerr); -fprintf('relative l_inf/l_1 error %5.2e\n',relerr2); - - -%% - - -%% -%---------------------------------------- -% -% -% Auxiliary routines for generating boundary -% - -function [r,d,d2] = circulararc(t,cpars) -%%circulararc -% return position, first and second derivatives of a circular arc -% that passes through points (x,y)=(a,0) and (x,y)=(b,0) with opening -% angle theta0. -% -% Inputs: -% t - paramter values (-theta0/2,theta0/2) to evaluate these quantities -% -% Outputs: -% r - coordinates -% d - first derivatives w.r.t. t -% d2 - second derivatives w.r.t. t - -v0 = cpars.v0; -v1 = cpars.v1; -theta0 = cpars.theta; -ifconvex = cpars.ifconvex; - -a = v0(1); -b = v1(1); - -islocal = -1; -if isfield(cpars,'islocal') - islocal = cpars.islocal; -end - -if islocal == -1 - cx = (a+b)/2; - r0 = (b-a)/2/sin(theta0/2); - - if ifconvex - cy = v0(2)+(b-a)/2/tan(theta0/2); - theta = 3*pi/2+t; - else - cy = v0(2)-(b-a)/2/tan(theta0/2); - theta = pi/2+t; - end - - - xs = r0*cos(theta); - ys = r0*sin(theta); - - xp = -ys; - yp = xs; - - xpp = -xs; - ypp = -ys; - - xs = cx + xs; - ys = cy + ys; - -else - - r0 = (b-a)/2/sin(theta0/2); - if ~ifconvex - if islocal == 0 - cx = r0*cos(pi/2+theta0/2); - sx = r0*sin(pi/2+theta0/2); - elseif islocal == 1 - cx = r0*cos(pi/2-theta0/2); - sx = r0*sin(pi/2-theta0/2); - end - elseif ifconvex - if islocal == 0 - cx = r0*cos(3*pi/2+theta0/2); - sx = r0*sin(3*pi/2+theta0/2); - elseif islocal == 1 - cx = r0*cos(3*pi/2-theta0/2); - sx = r0*sin(3*pi/2-theta0/2); - end - end - -% t2 = t.*t; -% t3 = t2.*t; -% t4 = t3.*t; -% t5 = t4.*t; -% t6 = t5.*t; -% t7 = t6.*t; -% t8 = t7.*t; -% t9 = t8.*t; -% t10 = t9.*t; -% -% -% ct = -t2/2 + t4/24 - t6/720 + t8/40320 - t10/3628800; -% st = t - t3/6 + t5/120 - t7/5040 + t9/362880; - -% better to use sin(t) directly instead of its series expansion because -% 1. the evaluation of sin(t) is accurate at t=0; -% 2. one needs to determine how many terms are needed in series expansion, -% which looks simple but it could be either inefficient or inaccurate -% depending on the value of t; 3. if we write "if" statements, it's -% difficult to vectorize the evaluation. - -% same reason that we should use -2*sin(t/2).^2 instead of its Taylor -% expansion, but of course one should NOT use cos(t)-1. - ct = -2*sin(t/2).^2; - st = sin(t); - - ctp = -st; - stp = cos(t); - - ctpp = -stp; - stpp = -st; - - xs = -sx*st + cx*ct; - ys = cx*st + sx*ct; - - xp = -sx*stp + cx*ctp; - yp = cx*stp + sx*ctp; - - xpp = -sx*stpp + cx*ctpp; - ypp = cx*stpp + sx*ctpp; -end - - -r = [(xs(:)).'; (ys(:)).']; -d = [(xp(:)).'; (yp(:)).']; -d2 = [(xpp(:)).'; (ypp(:)).']; - -end - - - diff --git a/devtools/test/refine_trefinementTest.m b/devtools/test/refine_trefinementTest.m deleted file mode 100644 index fde86aa..0000000 --- a/devtools/test/refine_trefinementTest.m +++ /dev/null @@ -1,83 +0,0 @@ -%REFINE_TREFINEMENTTEST -% -% get chunker description of a starfish domain. check if that domain -% satisfies a level restriction in the underlying parameter space. -% call a refinement routine to fix. - -clearvars; close all; -iseed = 8675309; -rng(iseed); - -addpaths_loc(); - -% define curve - -cparams = []; -cparams.eps = 1.0e-10; -cparams.nover = 0; -pref = []; -pref.k = 16; -narms = 5; -amp = 0.25; - -% form curve without strict enforcement of level restriction -% in h - -start = tic; chnkr = chunkerfunc(@(t) starfish(t,narms,amp),cparams,pref); -t1 = toc(start); -fprintf('%5.2e s : time to build geo\n',t1) -fprintf('originally, nch = %d\n',chnkr.nch) - -% check for violations of level restriction - -happy = true; -for i = 1:chnkr.nch - hself = chnkr.h(i); - i1 = chnkr.adj(1,i); - i2 = chnkr.adj(2,i); - h1 = hself; h2 = hself; - if (i1 > 0) - h1 = chnkr.h(i1); - end - if (i2 > 0) - h2 = chnkr.h(i2); - end - if (hself > 2*h1) - happy = false; - end - if (hself > 2*h2) - happy = false; - end -end - -% enforce level restriction using refine code with -% specific options - -opts.lvlr = 't'; % refine in t rather than arclength -% standard is 2.05 which is too loose for our purposes -opts.lvlrfac = 1.99; - -start = tic; chnkr = refine(chnkr,opts); t1 = toc(start); -fprintf('%5.2e s : time to refine geo\n',t1) -fprintf('now, nch = %d\n',chnkr.nch) - -happy = true; -for i = 1:chnkr.nch - hself = chnkr.h(i); - i1 = chnkr.adj(1,i); - i2 = chnkr.adj(2,i); - h1 = hself; h2 = hself; - if (i1 > 0) - h1 = chnkr.h(i1); - end - if (i2 > 0) - h2 = chnkr.h(i2); - end - if (hself > 2*h1) - happy = false; - end - if (hself > 2*h2) - happy = false; - end -end -