diff --git a/.readthedocs.yaml b/.readthedocs.yaml new file mode 100644 index 0000000..8325091 --- /dev/null +++ b/.readthedocs.yaml @@ -0,0 +1,32 @@ +# .readthedocs.yaml +# Read the Docs configuration file +# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details + +# Required +version: 2 + +# Set the OS, Python version and other tools you might need +build: + os: ubuntu-22.04 + tools: + python: "3.12" + # You can also specify other tool versions: + # nodejs: "19" + # rust: "1.64" + # golang: "1.19" + +# Build documentation in the "docs/" directory with Sphinx +sphinx: + configuration: docs/conf.py + +# Optionally build your docs in additional formats such as PDF and ePub +# formats: +# - pdf +# - epub + +# Optional but recommended, declare the Python requirements required +# to build your documentation +# See https://docs.readthedocs.io/en/stable/guides/reproducible-builds.html +python: + install: + - requirements: docs/requirements.txt diff --git a/chunkie/+chnk/+elast2d/kern.m b/chunkie/+chnk/+elast2d/kern.m index f789a89..075708c 100644 --- a/chunkie/+chnk/+elast2d/kern.m +++ b/chunkie/+chnk/+elast2d/kern.m @@ -33,12 +33,17 @@ % - type: string % * type == 'S', single layer kernel % * type == 'Strac', traction of single layer kernel +% * type == 'Sgrad', gradient of single layer kernel +% kernel is 4 x 2, where entries are organized as +% Sgrad = [grad S(1,1) grad S(1,2) +% grad S(2,1) grad S(2,2)] % * type == 'D', double layer kernel % * type == 'Dalt', alternative (smooth) double layer kernel % * type == 'Daltgrad', gradient of alternative (smooth) double layer % kernel is 4 x 2, where entries are organized as % Daltgrad = [grad Dalt(1,1) grad Dalt(1,2) % grad Dalt(2,1) grad Dalt(2,2)] +% * type == 'Dalttrac', traction of alternative (smooth) double layer % % outpts % - mat: (2 nt) x (2 ns) matrix @@ -53,8 +58,8 @@ eta = mu/(2*pi*(lam+2*mu)); zeta = (lam+mu)/(pi*(lam+2*mu)); -n = size(s.r,2); -m = size(t.r,2); +n = size(s.r(:,:),2); +m = size(t.r(:,:),2); x = (t.r(1,:)).' - s.r(1,:); y = (t.r(2,:)).' - s.r(2,:); @@ -81,6 +86,37 @@ kron(rdotv./r2,[1 0; 0 1])); mat = term1+term2; end +if (strcmpi(type,'sgrad')) + mat = zeros(4*m,2*n); + tmp = beta*x./r2; + mat(1:4:end,1:2:end) = tmp; + mat(3:4:end,2:2:end) = tmp; + tmp = beta*y./r2; + mat(2:4:end,1:2:end) = tmp; + mat(4:4:end,2:2:end) = tmp; + + % x der of x^2/ r^2 + tmp = gamma*(2*r2.*x - x.^2.*(2*x))./r4; + mat(1:4:end,1:2:end) = mat(1:4:end,1:2:end) + tmp; + % y der of x^2/ r^2 + tmp = gamma*(-x.^2.*(2*y))./r4; + mat(2:4:end,1:2:end) = mat(2:4:end,1:2:end) + tmp; + % x der of xy/ r^2 + tmp = gamma*(r2.*y-x.*y.*(2*x))./r4; + mat(3:4:end,1:2:end) = mat(3:4:end,1:2:end) + tmp; + mat(1:4:end,2:2:end) = mat(1:4:end,2:2:end) + tmp; + % y der of xy/ r^2 + tmp = gamma*(r2.*x-x.*y.*(2*y))./r4; + mat(4:4:end,1:2:end) = mat(4:4:end,1:2:end) + tmp; + mat(2:4:end,2:2:end) = mat(2:4:end,2:2:end) + tmp; + % x der of y^2/ r^2 + tmp = gamma*(- y.^2.*(2*x))./r4; + mat(3:4:end,2:2:end) = mat(3:4:end,2:2:end) + tmp; + % y der of y^2/ r^2 + tmp = gamma*(2*r2.*y-y.^2.*(2*y))./r4; + mat(4:4:end,2:2:end) = mat(4:4:end,2:2:end) + tmp; + +end if (strcmpi(type,'d')) dirx = s.n(1,:); dirx = dirx(:).'; diry = s.n(2,:); diry = diry(:).'; @@ -151,4 +187,59 @@ end +if (strcmpi(type,'dalttrac')) + dirx = s.n(1,:); dirx = dirx(:).'; + diry = s.n(2,:); diry = diry(:).'; + rdotv = x.*dirx + y.*diry; + r6 = r4.*r2; + + matg = zeros(4*m,2*n); + + % i = 1, j = 1, l = 1 + aij_xl = -zeta*(-4*x.^3.*rdotv./r6+(2*x.*rdotv+x.^2.*dirx)./r4) ... + -2*eta*(dirx./r2 - 2*rdotv.*x./r4); + matg(1:4:end,1:2:end) = aij_xl; + + % i = 1, j = 1, l = 2 + aij_xl = -zeta*(-4*x.^2.*y.*rdotv./r6+(x.^2.*diry)./r4) ... + -2*eta*(diry./r2 - 2*rdotv.*y./r4); + matg(2:4:end,1:2:end) = aij_xl; + + % i = 2, j = 1, l = 1 + aij_xl = -zeta*(-4*x.^2.*y.*rdotv./r6+(y.*rdotv+x.*y.*dirx)./r4); + matg(3:4:end,1:2:end) = aij_xl; + + % i = 2, j = 1, l = 2 + aij_xl = -zeta*(-4*x.*y.^2.*rdotv./r6+(x.*rdotv+x.*y.*diry)./r4); + matg(4:4:end,1:2:end) = aij_xl; + + % i = 1, j = 2, l = 1 + aij_xl = -zeta*(-4*x.^2.*y.*rdotv./r6+(y.*rdotv+x.*y.*dirx)./r4); + matg(1:4:end,2:2:end) = aij_xl; + + % i = 1, j = 2, l = 2 + aij_xl = -zeta*(-4*x.*y.^2.*rdotv./r6+(x.*rdotv + x.*y.*diry)./r4); + matg(2:4:end,2:2:end) = aij_xl; + + % i = 2, j = 2, l = 1 + aij_xl = -zeta*(-4*y.^2.*x.*rdotv./r6+(y.^2.*dirx)./r4) ... + -2*eta*(dirx./r2 - 2*rdotv.*x./r4); + matg(3:4:end,2:2:end) = aij_xl; + + % i = 2, j = 2, l = 2 + aij_xl = -zeta*(-4*y.^3.*rdotv./r6+(2*y.*rdotv+y.^2.*diry)./r4) ... + -2*eta*(diry./r2 - 2*rdotv.*y./r4); + matg(4:4:end,2:2:end) = aij_xl; + + mat = zeros(2*m,2*n); + n1 = t.n(1,:); n1 = n1(:); n2 = t.n(2,:); n2 = n2(:); + tmp = matg(1:4:end,:)+matg(4:4:end,:); + mat(1:2:end,:) = lam*n1.*tmp; + mat(2:2:end,:) = lam*n2.*tmp; + tmp = mu*(matg(2:4:end,:) + matg(3:4:end,:)); + mat(1:2:end,:) = mat(1:2:end,:) + tmp.*n2 + 2*mu*matg(1:4:end,:).*n1; + mat(2:2:end,:) = mat(2:2:end,:) + tmp.*n1 + 2*mu*matg(4:4:end,:).*n2; + +end + end diff --git a/chunkie/+chnk/+flam/kernbyindexr.m b/chunkie/+chnk/+flam/kernbyindexr.m index 8f6b1eb..b280657 100644 --- a/chunkie/+chnk/+flam/kernbyindexr.m +++ b/chunkie/+chnk/+flam/kernbyindexr.m @@ -1,4 +1,4 @@ -function mat = kernbyindexr(i,j,targs,chnkr,kern,opdims,spmat) +function mat = kernbyindexr(i,j,targobj,chnkr,kern,opdims,spmat) %% evaluate system matrix by entry index utility function for % general kernels, with replacement for specific entries and the % ability to add a low-rank modification, rectangular version @@ -17,6 +17,10 @@ % % i - array of row indices to compute % j - array of col indices to compute +% targobj - object describing the target points, can be specified as +% * array of points +% * chunker object +% * chunkgraph object % chnkr - chunker object describing boundary % whts - smooth integration weights on chnkr % kern - kernel function of the form kern(s,t,stau,ttau) where s and t @@ -42,15 +46,32 @@ [juni,~,ijuni] = unique(jpts); % matrix-valued entries of kernel for unique points - -ri = targs(:,iuni); rj = chnkr.r(:,juni); +rj = chnkr.r(:,juni); dj = chnkr.d(:,juni); d2j = chnkr.d2(:,juni); nj = chnkr.n(:,juni); -%di = bsxfun(@rdivide,di,sqrt(sum(di.^2,1))); -%dj = bsxfun(@rdivide,dj,sqrt(sum(dj.^2,1))); srcinfo = []; srcinfo.r = rj; srcinfo.d = dj; srcinfo.d2 = d2j; srcinfo.n = nj; -targinfo = []; targinfo.r = ri; +% Assign appropriate object to targinfo +targinfo = []; +if isa(targobj, "chunker") || isa(targobj, "chunkgraph") + targinfo.r = targobj.r(:,iuni); + targinfo.d = targobj.d(:,iuni); + targinfo.d2 = targobj.d2(:,iuni); + targinfo.n = targobj.n(:,iuni); +elseif isstruct(targobj) + if isfield(targobj,'d') + targinfo.d = targobj.d(:,iuni); + end + if isfield(targobj,'d2') + targinfo.d = targobj.d(:,iuni); + end + if isfield(targobj,'n') + targinfo.n = targobj.n(:,iuni); + end + targinfo.r = targobj.r(:,iuni); +else + targinfo.r = targobj(:,iuni); +end matuni = kern(srcinfo,targinfo); diff --git a/chunkie/+chnk/+flam/proxy_square_pts.m b/chunkie/+chnk/+flam/proxy_square_pts.m index 1d7bb7f..3f1d362 100644 --- a/chunkie/+chnk/+flam/proxy_square_pts.m +++ b/chunkie/+chnk/+flam/proxy_square_pts.m @@ -1,4 +1,4 @@ -function [pr,ptau,pw,pin] = proxy_square_pts(porder) +function [pr,ptau,pw,pin] = proxy_square_pts(porder, opts) %PROXY_SQUARE_PTS return the proxy points, unit tangents, weights, and % function handle for determining if points are within proxy surface. % This function is for the proxy surface around a unit box centered at @@ -16,12 +16,31 @@ porder = 64; end +if nargin < 2 + opts = []; +end + +iflege = 0; +if isfield(opts, 'iflege') + iflege = opts.iflege; +end + + assert(mod(porder,4) == 0, ... 'number of proxy points on square should be multiple of 4') po4 = porder/4; -pts = -1.5+3*(0:(po4-1))*1.0/po4; +if ~iflege + pts = -1.5+3*(0:(po4-1))*1.0/po4; + pw = ones(porder,1)*(4.0*3.0/porder); +else + [xleg, wleg] = lege.exps(po4); + pts = xleg.'*1.5; + wleg = wleg*1.5; + pw = [wleg; wleg; wleg; wleg]; +end + one = ones(size(pts)); pr = [pts, one*1.5, -pts, -1.5*one; @@ -30,8 +49,6 @@ ptau = [one, one*0, -one, one*0; one*0, one, one*0, -one]; -pw = ones(porder,1)*(4.0*3.0/porder); - pin = @(x) max(abs(x),[],1) < 1.5; end \ No newline at end of file diff --git a/chunkie/+chnk/+helm2d/fmm.m b/chunkie/+chnk/+helm2d/fmm.m index 4d50789..d7eceb0 100644 --- a/chunkie/+chnk/+helm2d/fmm.m +++ b/chunkie/+chnk/+helm2d/fmm.m @@ -38,6 +38,8 @@ % note that targinfo must contain targ.n % type = 'dprime' normal derivative of double layer, % note that targinfo must contain targ.n +% type = 'cprime' normal derivative of combined layer, +% note that targinfo must contain targ.n % sigma - density % varargin{1} - coef in the combined layer formula, otherwise % does nothing @@ -61,7 +63,7 @@ case {'d', 'dprime'} srcuse.dipstr = sigma(:).'; srcuse.dipvec = srcinfo.n(1:2,:); - case 'c' + case {'c', 'cprime'} coefs = varargin{1}; srcuse.charges = coefs(2)*sigma(:).'; srcuse.dipstr = coefs(1)*sigma(:).'; @@ -83,7 +85,7 @@ switch lower(type) case {'s', 'd', 'c'} varargout{1} = U.pottarg.'; - case {'sprime', 'dprime'} + case {'sprime', 'dprime', 'cprime'} if ( ~isfield(targinfo, 'n') ) error('CHUNKIE:helm2d:fmm:normals', ... 'Targets require normal info when evaluating Helmholtz kernel ''%s''.', type); diff --git a/chunkie/+chnk/+helm2d/kern.m b/chunkie/+chnk/+helm2d/kern.m index a90e636..d697efd 100644 --- a/chunkie/+chnk/+helm2d/kern.m +++ b/chunkie/+chnk/+helm2d/kern.m @@ -133,14 +133,14 @@ % normal derivative of combined field if strcmpi(type,'cprime') coef = ones(2,1); - if(nargin == 5); coef = varargin{1}; end; + if(nargin == 5); coef = varargin{1}; end targnorm = targinfo.n; srcnorm = srcinfo.n; - submat = zeros(nt,ns); + % Get gradient and hessian info - [submats,grad,hess] = chnk.helm2d.green(zk,src,targ); + [~,grad,hess] = chnk.helm2d.green(zk,src,targ); nxsrc = repmat(srcnorm(1,:),nt,1); nysrc = repmat(srcnorm(2,:),nt,1); @@ -161,7 +161,7 @@ % Dirichlet and neumann data corresponding to combined field if strcmpi(type,'c2trans') coef = ones(2,1); - if(nargin == 5); coef = varargin{1}; end; + if(nargin == 5); coef = varargin{1}; end targnorm = targinfo.n; srcnorm = srcinfo.n; diff --git a/chunkie/+chnk/+helm2d/transmission_helper.m b/chunkie/+chnk/+helm2d/transmission_helper.m index 553dd34..b199e36 100644 --- a/chunkie/+chnk/+helm2d/transmission_helper.m +++ b/chunkie/+chnk/+helm2d/transmission_helper.m @@ -151,7 +151,7 @@ c1 = coefs(d1); c2 = coefs(d2); ind1 = sum(nchs(1:i-1))*ngl*2+(1:2:2*nchs(i)*ngl); - ind2 = sum(nchs(1:i-1))*ngl*2+(1:2:2*nchs(i)*ngl); + ind2 = sum(nchs(1:i-1))*ngl*2+(2:2:2*nchs(i)*ngl); targnorm = chnkrs(i).n; nx = targnorm(1,:); nx = nx(:); ny = targnorm(2,:); ny = ny(:); 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 262180b..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,15 +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(abs(ds).^2,1)); -%taus = bsxfun(@rdivide,ds,dsnrms); -%dtnrms = sqrt(sum(abs(dt).^2,1)); -%taut = bsxfun(@rdivide,dt,dtnrms); - -ws = kron(hs(:),wts(:)); +dsnrms = sqrt(sum(ds.^2,1)); +ws = repmat(wts(:),length(j),1); dsdt = dsnrms(:).*ws; 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/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 9533aeb..438cd45 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,24 +13,78 @@ % % 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 + +if(size(fkern)==1) + fkernlocal = fkern; +else + fkernlocal(nedge,nedge) = kernel(); + for i=1:nedge + ici = iedgechunks(1,i); + for j=1:nedge + icj = iedgechunks(1,j); + fkernlocal(i,j) = fkern(ici,icj); + end + end + +end + +rcipsav.fkernlocal = fkernlocal; + +% get coefficients of recentered edge chunks and figure out orientation km1 = k-1; rcs = zeros(km1,dim,nedge); @@ -42,6 +94,8 @@ d2scal = zeros(nedge,1); ctr = zeros(dim,nedge); +ileftright = zeros(nedge,1); +nextchunk = zeros(nedge,1); for i = 1:nedge ic = iedgechunks(1,i); @@ -52,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; @@ -62,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; @@ -72,17 +125,26 @@ 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 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 = []; @@ -125,6 +187,9 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% begin recursion proper + h0=ones(nedge,1); for level=1:nsub h = h0/2^(nsub-level); @@ -161,7 +226,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; @@ -194,8 +259,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/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/+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/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 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/+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/chunk_nearparam.m b/chunkie/+chnk/chunk_nearparam.m new file mode 100644 index 0000000..efdf2fd --- /dev/null +++ b/chunkie/+chnk/chunk_nearparam.m @@ -0,0 +1,223 @@ +function [ts,rs,ds,d2s,dist2s] = chunk_nearparam(rval,pts,opts,t,u) +%CHNK.CHUNK_NEARPARAM Find nearest point on a single chunk for given +% reference point. This routine is not intended to be user-callable +% +% Syntax: [rn,dn,d2n,dist,tn,ichn] = nearest(rval,pts,opts,t,u) +% +% Input: +% rval - dim x k array of chunk positions +% pts - the reference point whose distance to curve is being computed +% +% Optional input: +% u - the matrix created by lege.exps mapping point values to coefs for +% order chnkr.k +% opts - options structure +% opts.thresh - threshold for newton (1e-9) +% opts.nitermax - maximum iterations for newton (200) +% +% Output: +% ts - ts(j) in [-1,1] is preimage on chunk that is closest to pts(:,j) +% rs - rs(:,j) chunk coordinates at ts(j) +% + +% author: Travis Askham (askhamwhat@gmail.com) + +maxnewt = 15; +thresh0 = 1.0d-14; + +if nargin < 4 || isempty(t) || nargin < 5 || isempty(u) + [t,~,u] = lege.exps(chnkr.k); +end + +if nargin < 3 + opts = []; +end + +if isfield(opts,'nitermax') + maxnewt = opts.nitermax; +end +if isfield(opts,'thresh') + 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); + +k = size(rval,2); + +rc = u*(rval.'); +drc = lege.derpol(rc); drc = [drc;zeros(1,dim)]; +d2rc = lege.derpol(drc); d2rc = [d2rc;zeros(1,dim)]; +cfs = [rc.';drc.';d2rc.'].'; + +dist2all = reshape(sum( abs(reshape(pts,dim,1,npts) ... + - reshape(rval,dim,k,1)).^2, 1),k,npts); +[dist2all,ipt] = min(dist2all,[],1); + +rs = zeros(dim,npts); +ds = zeros(dim,npts); +d2s = zeros(dim,npts); +dist2s = zeros(npts,1); +ts = zeros(npts,1); + +thresh = thresh0 * (k^2*sum(abs(drc(:))) + k*sum(abs(rc(:)))); + +for i = 1:npts + + ref = pts(:,i); + + % closest point on grid + + t0 = t(ipt(i)); + + all0 = lege.exev(t0,cfs); + r0 = all0(1:dim); + d0 = all0(dim+1:2*dim); + d20 = all0(2*dim+1:end); + + ts(i) = t(ipt(i)); + rs(:,i) = r0; + ds(:,i) = d0; + d2s(:,i) = d20; + dist2s(i) = dist2all(i); + + % try newton + + rdiff0 = r0(:) - ref; + dprime = rdiff0.'*d0(:); + dprime2 = d0(:).'*d0(:) + d0(:).'*d20(:); + + newtsuccess = false; + + for iii = 1:maxnewt + + % Newton step + + deltat = - dprime/dprime2; + + t0 = t0+deltat; + + if (t0 > 1.0) + t0 = 1; + end + if (t0 < -1.0) + t0 = -1; + end + + all0 = (lege.exev(t0,cfs)); + r0 = all0(1:dim); + d0 = all0(dim+1:2*dim); + d20 = all0(2*dim+1:end); + + rdiff0 = r0(:)-ref(:); + dprime = rdiff0.'*d0(:); + dprime2 = d0(:).'*d0(:) + d0(:).'*d20(:); + + if abs(dprime) < thresh + newtsuccess = true; + break; + end + + if t0 == 1 && dprime < 0 + newtsuccess = true; + break; + end + if t0 == -1 && dprime > 0 + newtsuccess = true; + break; + end + + end + + dist2newt = sum(abs(rdiff0.^2)); + + if dist2newt > dist2all(i) + newtsuccess = false; + else + ts(i) = t0; + rs(:,i) = r0; + ds(:,i) = d0; + d2s(:,i) = d20; + dist2s(i) = dist2newt; + end + + if ~newtsuccess + % try levenberg if Newton fails + t0 = t(ipt(i)); + + all0 = lege.exev(t0,cfs); + r0 = all0(1:dim); + d0 = all0(dim+1:2*dim); + d20 = all0(2*dim+1:end); + + rdiff0 = r0(:) - ref; + dprime = rdiff0.'*d0(:); + dprime2 = d0(:).'*d0(:); + lam = dprime2; + + dist0 = sum(abs(rdiff0).^2); + + for iii = 1:maxnewt + + % Levenberg step + + deltat = - dprime/(dprime2+lam); + + t1 = t0+deltat; + + if (t1 > 1.0) + t1 = 1; + end + if (t1 < -1.0) + t1 = -1; + end + + all1 = (lege.exev(t1,cfs)); + r1 = all1(1:dim); + d1 = all1(dim+1:2*dim); + d21 = all1(2*dim+1:end); + + rdiff1 = r1(:)-ref(:); + dist1 = sum(abs(rdiff1).^2); + + if dist1 > dist0 + lam = lam*2; + else + t0 = t1; + r0 = r1; + d0 = d1; + d20 = d21; + rdiff0 = rdiff1; + dprime = rdiff0.'*d0(:); + dprime2 = d0(:).'*d0(:) + d0(:).'*d20(:); + dist0 = dist1; + lam = lam/3; + if abs(dprime) < thresh + break; + end + + if t0 == 1 && dprime < 0 + break; + end + if t0 == -1 && dprime > 0 + break; + end + end + end + + ts(i) = t0; + rs(:,i) = r0; + ds(:,i) = d0; + d2s(:,i) = d20; + dist2s(i) = dist0; + end + +end + diff --git a/chunkie/+chnk/ellipse_oversample.m b/chunkie/+chnk/ellipse_oversample.m new file mode 100644 index 0000000..5e50137 --- /dev/null +++ b/chunkie/+chnk/ellipse_oversample.m @@ -0,0 +1,577 @@ +function nleg = ellipse_oversample(rho,npoly,tol) +%CHNK.ELLIPSE_OVERSAMPLE get oversampling factors based on bernstein +% ellipses +% +% input +% +% rho - bernstein ellipse parameter +% npoly - degree of polynomial*log + polynomial to integrate +% tol - tolerance for integral +% +% output +% +% nleg - the order of the legendre rule to use to get sufficient accuracy +% outside of the image of the bernstein ellipse under the boundary +% parameterization. + +[iprecs,npolys,rhos,nlegtab] = loadtab(); + +precs = 10.^(-iprecs); + +if (tol < min(precs)) + warning('requested tolerance out of range, oversampling may be insufficient'); +end +if (rho < min(rhos)) + warning('requested bernstein parameter out of range, oversampling may be insufficient'); +end +if (npoly > max(npolys)) + warning('requested poly degree out of range, oversampling may be insufficient'); +end + +iprec = min(nnz(precs>tol)+1,length(iprecs)); +ipol = min(nnz(npolys < npoly) + 1,length(npolys)); +irho = max(length(rhos)-nnz(rhos > rho),1); + +nleg = nlegtab(irho,ipol,iprec); + + +end + + +function [iprecs,npolys,rhos,nlegtab] = loadtab() +iprecs = zeros(5,1); +npolys = zeros(20,1); +rhos = zeros(5,1); +nlegtab = zeros(5,20,5); + iprecs( 1) = 3; + iprecs( 2) = 6; + iprecs( 3) = 9; + iprecs( 4) = 12; + iprecs( 5) = 15; + npolys( 1) = 4; + npolys( 2) = 8; + npolys( 3) = 12; + npolys( 4) = 16; + npolys( 5) = 20; + npolys( 6) = 24; + npolys( 7) = 28; + npolys( 8) = 32; + npolys( 9) = 36; + npolys( 10) = 40; + npolys( 11) = 44; + npolys( 12) = 48; + npolys( 13) = 52; + npolys( 14) = 56; + npolys( 15) = 60; + npolys( 16) = 64; + npolys( 17) = 68; + npolys( 18) = 72; + npolys( 19) = 76; + npolys( 20) = 80; + rhos( 1) =0.12000000000000000000E+01; + rhos( 2) =0.14000000000000000000E+01; + rhos( 3) =0.16000000000000000000E+01; + rhos( 4) =0.18000000000000000000E+01; + rhos( 5) =0.20000000000000000000E+01; + nlegtab( 1, 1, 1) = 14; + nlegtab( 2, 1, 1) = 9; + nlegtab( 3, 1, 1) = 7; + nlegtab( 4, 1, 1) = 6; + nlegtab( 5, 1, 1) = 5; + nlegtab( 1, 2, 1) = 14; + nlegtab( 2, 2, 1) = 9; + nlegtab( 3, 2, 1) = 7; + nlegtab( 4, 2, 1) = 6; + nlegtab( 5, 2, 1) = 6; + nlegtab( 1, 3, 1) = 14; + nlegtab( 2, 3, 1) = 9; + nlegtab( 3, 3, 1) = 8; + nlegtab( 4, 3, 1) = 7; + nlegtab( 5, 3, 1) = 7; + nlegtab( 1, 4, 1) = 14; + nlegtab( 2, 4, 1) = 9; + nlegtab( 3, 4, 1) = 8; + nlegtab( 4, 4, 1) = 8; + nlegtab( 5, 4, 1) = 8; + nlegtab( 1, 5, 1) = 14; + nlegtab( 2, 5, 1) = 10; + nlegtab( 3, 5, 1) = 9; + nlegtab( 4, 5, 1) = 9; + nlegtab( 5, 5, 1) = 8; + nlegtab( 1, 6, 1) = 14; + nlegtab( 2, 6, 1) = 10; + nlegtab( 3, 6, 1) = 10; + nlegtab( 4, 6, 1) = 9; + nlegtab( 5, 6, 1) = 9; + nlegtab( 1, 7, 1) = 14; + nlegtab( 2, 7, 1) = 11; + nlegtab( 3, 7, 1) = 10; + nlegtab( 4, 7, 1) = 10; + nlegtab( 5, 7, 1) = 10; + nlegtab( 1, 8, 1) = 14; + nlegtab( 2, 8, 1) = 11; + nlegtab( 3, 8, 1) = 11; + nlegtab( 4, 8, 1) = 10; + nlegtab( 5, 8, 1) = 10; + nlegtab( 1, 9, 1) = 14; + nlegtab( 2, 9, 1) = 12; + nlegtab( 3, 9, 1) = 11; + nlegtab( 4, 9, 1) = 11; + nlegtab( 5, 9, 1) = 11; + nlegtab( 1, 10, 1) = 14; + nlegtab( 2, 10, 1) = 12; + nlegtab( 3, 10, 1) = 12; + nlegtab( 4, 10, 1) = 11; + nlegtab( 5, 10, 1) = 11; + nlegtab( 1, 11, 1) = 14; + nlegtab( 2, 11, 1) = 13; + nlegtab( 3, 11, 1) = 12; + nlegtab( 4, 11, 1) = 12; + nlegtab( 5, 11, 1) = 12; + nlegtab( 1, 12, 1) = 14; + nlegtab( 2, 12, 1) = 13; + nlegtab( 3, 12, 1) = 12; + nlegtab( 4, 12, 1) = 12; + nlegtab( 5, 12, 1) = 12; + nlegtab( 1, 13, 1) = 15; + nlegtab( 2, 13, 1) = 13; + nlegtab( 3, 13, 1) = 13; + nlegtab( 4, 13, 1) = 13; + nlegtab( 5, 13, 1) = 12; + nlegtab( 1, 14, 1) = 15; + nlegtab( 2, 14, 1) = 14; + nlegtab( 3, 14, 1) = 13; + nlegtab( 4, 14, 1) = 13; + nlegtab( 5, 14, 1) = 13; + nlegtab( 1, 15, 1) = 15; + nlegtab( 2, 15, 1) = 14; + nlegtab( 3, 15, 1) = 14; + nlegtab( 4, 15, 1) = 13; + nlegtab( 5, 15, 1) = 13; + nlegtab( 1, 16, 1) = 16; + nlegtab( 2, 16, 1) = 14; + nlegtab( 3, 16, 1) = 14; + nlegtab( 4, 16, 1) = 14; + nlegtab( 5, 16, 1) = 13; + nlegtab( 1, 17, 1) = 16; + nlegtab( 2, 17, 1) = 15; + nlegtab( 3, 17, 1) = 14; + nlegtab( 4, 17, 1) = 14; + nlegtab( 5, 17, 1) = 14; + nlegtab( 1, 18, 1) = 16; + nlegtab( 2, 18, 1) = 15; + nlegtab( 3, 18, 1) = 15; + nlegtab( 4, 18, 1) = 14; + nlegtab( 5, 18, 1) = 14; + nlegtab( 1, 19, 1) = 17; + nlegtab( 2, 19, 1) = 15; + nlegtab( 3, 19, 1) = 15; + nlegtab( 4, 19, 1) = 15; + nlegtab( 5, 19, 1) = 14; + nlegtab( 1, 20, 1) = 17; + nlegtab( 2, 20, 1) = 16; + nlegtab( 3, 20, 1) = 15; + nlegtab( 4, 20, 1) = 15; + nlegtab( 5, 20, 1) = 15; + nlegtab( 1, 1, 2) = 31; + nlegtab( 2, 1, 2) = 18; + nlegtab( 3, 1, 2) = 13; + nlegtab( 4, 1, 2) = 11; + nlegtab( 5, 1, 2) = 10; + nlegtab( 1, 2, 2) = 31; + nlegtab( 2, 2, 2) = 18; + nlegtab( 3, 2, 2) = 13; + nlegtab( 4, 2, 2) = 11; + nlegtab( 5, 2, 2) = 10; + nlegtab( 1, 3, 2) = 31; + nlegtab( 2, 3, 2) = 18; + nlegtab( 3, 3, 2) = 14; + nlegtab( 4, 3, 2) = 12; + nlegtab( 5, 3, 2) = 11; + nlegtab( 1, 4, 2) = 31; + nlegtab( 2, 4, 2) = 18; + nlegtab( 3, 4, 2) = 14; + nlegtab( 4, 4, 2) = 13; + nlegtab( 5, 4, 2) = 12; + nlegtab( 1, 5, 2) = 31; + nlegtab( 2, 5, 2) = 18; + nlegtab( 3, 5, 2) = 15; + nlegtab( 4, 5, 2) = 13; + nlegtab( 5, 5, 2) = 12; + nlegtab( 1, 6, 2) = 31; + nlegtab( 2, 6, 2) = 19; + nlegtab( 3, 6, 2) = 15; + nlegtab( 4, 6, 2) = 14; + nlegtab( 5, 6, 2) = 13; + nlegtab( 1, 7, 2) = 31; + nlegtab( 2, 7, 2) = 19; + nlegtab( 3, 7, 2) = 16; + nlegtab( 4, 7, 2) = 15; + nlegtab( 5, 7, 2) = 14; + nlegtab( 1, 8, 2) = 31; + nlegtab( 2, 8, 2) = 19; + nlegtab( 3, 8, 2) = 16; + nlegtab( 4, 8, 2) = 15; + nlegtab( 5, 8, 2) = 15; + nlegtab( 1, 9, 2) = 31; + nlegtab( 2, 9, 2) = 20; + nlegtab( 3, 9, 2) = 17; + nlegtab( 4, 9, 2) = 16; + nlegtab( 5, 9, 2) = 15; + nlegtab( 1, 10, 2) = 31; + nlegtab( 2, 10, 2) = 20; + nlegtab( 3, 10, 2) = 18; + nlegtab( 4, 10, 2) = 17; + nlegtab( 5, 10, 2) = 16; + nlegtab( 1, 11, 2) = 31; + nlegtab( 2, 11, 2) = 21; + nlegtab( 3, 11, 2) = 18; + nlegtab( 4, 11, 2) = 17; + nlegtab( 5, 11, 2) = 17; + nlegtab( 1, 12, 2) = 31; + nlegtab( 2, 12, 2) = 21; + nlegtab( 3, 12, 2) = 19; + nlegtab( 4, 12, 2) = 18; + nlegtab( 5, 12, 2) = 18; + nlegtab( 1, 13, 2) = 31; + nlegtab( 2, 13, 2) = 21; + nlegtab( 3, 13, 2) = 19; + nlegtab( 4, 13, 2) = 19; + nlegtab( 5, 13, 2) = 18; + nlegtab( 1, 14, 2) = 31; + nlegtab( 2, 14, 2) = 22; + nlegtab( 3, 14, 2) = 20; + nlegtab( 4, 14, 2) = 19; + nlegtab( 5, 14, 2) = 19; + nlegtab( 1, 15, 2) = 31; + nlegtab( 2, 15, 2) = 22; + nlegtab( 3, 15, 2) = 20; + nlegtab( 4, 15, 2) = 20; + nlegtab( 5, 15, 2) = 20; + nlegtab( 1, 16, 2) = 31; + nlegtab( 2, 16, 2) = 23; + nlegtab( 3, 16, 2) = 21; + nlegtab( 4, 16, 2) = 20; + nlegtab( 5, 16, 2) = 20; + nlegtab( 1, 17, 2) = 31; + nlegtab( 2, 17, 2) = 23; + nlegtab( 3, 17, 2) = 21; + nlegtab( 4, 17, 2) = 21; + nlegtab( 5, 17, 2) = 21; + nlegtab( 1, 18, 2) = 31; + nlegtab( 2, 18, 2) = 23; + nlegtab( 3, 18, 2) = 22; + nlegtab( 4, 18, 2) = 22; + nlegtab( 5, 18, 2) = 21; + nlegtab( 1, 19, 2) = 31; + nlegtab( 2, 19, 2) = 24; + nlegtab( 3, 19, 2) = 22; + nlegtab( 4, 19, 2) = 22; + nlegtab( 5, 19, 2) = 22; + nlegtab( 1, 20, 2) = 32; + nlegtab( 2, 20, 2) = 24; + nlegtab( 3, 20, 2) = 23; + nlegtab( 4, 20, 2) = 23; + nlegtab( 5, 20, 2) = 22; + nlegtab( 1, 1, 3) = 48; + nlegtab( 2, 1, 3) = 28; + nlegtab( 3, 1, 3) = 20; + nlegtab( 4, 1, 3) = 16; + nlegtab( 5, 1, 3) = 14; + nlegtab( 1, 2, 3) = 48; + nlegtab( 2, 2, 3) = 28; + nlegtab( 3, 2, 3) = 20; + nlegtab( 4, 2, 3) = 17; + nlegtab( 5, 2, 3) = 15; + nlegtab( 1, 3, 3) = 48; + nlegtab( 2, 3, 3) = 28; + nlegtab( 3, 3, 3) = 21; + nlegtab( 4, 3, 3) = 17; + nlegtab( 5, 3, 3) = 16; + nlegtab( 1, 4, 3) = 48; + nlegtab( 2, 4, 3) = 28; + nlegtab( 3, 4, 3) = 21; + nlegtab( 4, 4, 3) = 18; + nlegtab( 5, 4, 3) = 16; + nlegtab( 1, 5, 3) = 48; + nlegtab( 2, 5, 3) = 28; + nlegtab( 3, 5, 3) = 22; + nlegtab( 4, 5, 3) = 19; + nlegtab( 5, 5, 3) = 17; + nlegtab( 1, 6, 3) = 48; + nlegtab( 2, 6, 3) = 28; + nlegtab( 3, 6, 3) = 22; + nlegtab( 4, 6, 3) = 19; + nlegtab( 5, 6, 3) = 18; + nlegtab( 1, 7, 3) = 48; + nlegtab( 2, 7, 3) = 28; + nlegtab( 3, 7, 3) = 23; + nlegtab( 4, 7, 3) = 20; + nlegtab( 5, 7, 3) = 18; + nlegtab( 1, 8, 3) = 48; + nlegtab( 2, 8, 3) = 29; + nlegtab( 3, 8, 3) = 23; + nlegtab( 4, 8, 3) = 20; + nlegtab( 5, 8, 3) = 19; + nlegtab( 1, 9, 3) = 48; + nlegtab( 2, 9, 3) = 29; + nlegtab( 3, 9, 3) = 24; + nlegtab( 4, 9, 3) = 21; + nlegtab( 5, 9, 3) = 20; + nlegtab( 1, 10, 3) = 48; + nlegtab( 2, 10, 3) = 30; + nlegtab( 3, 10, 3) = 24; + nlegtab( 4, 10, 3) = 22; + nlegtab( 5, 10, 3) = 20; + nlegtab( 1, 11, 3) = 48; + nlegtab( 2, 11, 3) = 30; + nlegtab( 3, 11, 3) = 25; + nlegtab( 4, 11, 3) = 22; + nlegtab( 5, 11, 3) = 21; + nlegtab( 1, 12, 3) = 48; + nlegtab( 2, 12, 3) = 30; + nlegtab( 3, 12, 3) = 25; + nlegtab( 4, 12, 3) = 23; + nlegtab( 5, 12, 3) = 22; + nlegtab( 1, 13, 3) = 48; + nlegtab( 2, 13, 3) = 31; + nlegtab( 3, 13, 3) = 26; + nlegtab( 4, 13, 3) = 24; + nlegtab( 5, 13, 3) = 23; + nlegtab( 1, 14, 3) = 48; + nlegtab( 2, 14, 3) = 31; + nlegtab( 3, 14, 3) = 26; + nlegtab( 4, 14, 3) = 24; + nlegtab( 5, 14, 3) = 23; + nlegtab( 1, 15, 3) = 48; + nlegtab( 2, 15, 3) = 31; + nlegtab( 3, 15, 3) = 27; + nlegtab( 4, 15, 3) = 25; + nlegtab( 5, 15, 3) = 24; + nlegtab( 1, 16, 3) = 48; + nlegtab( 2, 16, 3) = 32; + nlegtab( 3, 16, 3) = 27; + nlegtab( 4, 16, 3) = 25; + nlegtab( 5, 16, 3) = 25; + nlegtab( 1, 17, 3) = 48; + nlegtab( 2, 17, 3) = 32; + nlegtab( 3, 17, 3) = 28; + nlegtab( 4, 17, 3) = 26; + nlegtab( 5, 17, 3) = 26; + nlegtab( 1, 18, 3) = 48; + nlegtab( 2, 18, 3) = 32; + nlegtab( 3, 18, 3) = 28; + nlegtab( 4, 18, 3) = 27; + nlegtab( 5, 18, 3) = 26; + nlegtab( 1, 19, 3) = 49; + nlegtab( 2, 19, 3) = 33; + nlegtab( 3, 19, 3) = 29; + nlegtab( 4, 19, 3) = 27; + nlegtab( 5, 19, 3) = 27; + nlegtab( 1, 20, 3) = 49; + nlegtab( 2, 20, 3) = 33; + nlegtab( 3, 20, 3) = 29; + nlegtab( 4, 20, 3) = 28; + nlegtab( 5, 20, 3) = 28; + nlegtab( 1, 1, 4) = 66; + nlegtab( 2, 1, 4) = 37; + nlegtab( 3, 1, 4) = 27; + nlegtab( 4, 1, 4) = 22; + nlegtab( 5, 1, 4) = 19; + nlegtab( 1, 2, 4) = 66; + nlegtab( 2, 2, 4) = 37; + nlegtab( 3, 2, 4) = 27; + nlegtab( 4, 2, 4) = 22; + nlegtab( 5, 2, 4) = 20; + nlegtab( 1, 3, 4) = 66; + nlegtab( 2, 3, 4) = 37; + nlegtab( 3, 3, 4) = 28; + nlegtab( 4, 3, 4) = 23; + nlegtab( 5, 3, 4) = 20; + nlegtab( 1, 4, 4) = 66; + nlegtab( 2, 4, 4) = 37; + nlegtab( 3, 4, 4) = 28; + nlegtab( 4, 4, 4) = 24; + nlegtab( 5, 4, 4) = 21; + nlegtab( 1, 5, 4) = 66; + nlegtab( 2, 5, 4) = 38; + nlegtab( 3, 5, 4) = 29; + nlegtab( 4, 5, 4) = 24; + nlegtab( 5, 5, 4) = 22; + nlegtab( 1, 6, 4) = 66; + nlegtab( 2, 6, 4) = 38; + nlegtab( 3, 6, 4) = 29; + nlegtab( 4, 6, 4) = 25; + nlegtab( 5, 6, 4) = 22; + nlegtab( 1, 7, 4) = 66; + nlegtab( 2, 7, 4) = 38; + nlegtab( 3, 7, 4) = 30; + nlegtab( 4, 7, 4) = 25; + nlegtab( 5, 7, 4) = 23; + nlegtab( 1, 8, 4) = 66; + nlegtab( 2, 8, 4) = 39; + nlegtab( 3, 8, 4) = 30; + nlegtab( 4, 8, 4) = 26; + nlegtab( 5, 8, 4) = 24; + nlegtab( 1, 9, 4) = 66; + nlegtab( 2, 9, 4) = 39; + nlegtab( 3, 9, 4) = 31; + nlegtab( 4, 9, 4) = 27; + nlegtab( 5, 9, 4) = 24; + nlegtab( 1, 10, 4) = 66; + nlegtab( 2, 10, 4) = 39; + nlegtab( 3, 10, 4) = 31; + nlegtab( 4, 10, 4) = 27; + nlegtab( 5, 10, 4) = 25; + nlegtab( 1, 11, 4) = 66; + nlegtab( 2, 11, 4) = 40; + nlegtab( 3, 11, 4) = 31; + nlegtab( 4, 11, 4) = 28; + nlegtab( 5, 11, 4) = 26; + nlegtab( 1, 12, 4) = 66; + nlegtab( 2, 12, 4) = 40; + nlegtab( 3, 12, 4) = 32; + nlegtab( 4, 12, 4) = 28; + nlegtab( 5, 12, 4) = 26; + nlegtab( 1, 13, 4) = 66; + nlegtab( 2, 13, 4) = 40; + nlegtab( 3, 13, 4) = 32; + nlegtab( 4, 13, 4) = 29; + nlegtab( 5, 13, 4) = 27; + nlegtab( 1, 14, 4) = 66; + nlegtab( 2, 14, 4) = 41; + nlegtab( 3, 14, 4) = 33; + nlegtab( 4, 14, 4) = 30; + nlegtab( 5, 14, 4) = 28; + nlegtab( 1, 15, 4) = 66; + nlegtab( 2, 15, 4) = 41; + nlegtab( 3, 15, 4) = 33; + nlegtab( 4, 15, 4) = 30; + nlegtab( 5, 15, 4) = 28; + nlegtab( 1, 16, 4) = 66; + nlegtab( 2, 16, 4) = 41; + nlegtab( 3, 16, 4) = 34; + nlegtab( 4, 16, 4) = 31; + nlegtab( 5, 16, 4) = 29; + nlegtab( 1, 17, 4) = 66; + nlegtab( 2, 17, 4) = 42; + nlegtab( 3, 17, 4) = 34; + nlegtab( 4, 17, 4) = 31; + nlegtab( 5, 17, 4) = 30; + nlegtab( 1, 18, 4) = 66; + nlegtab( 2, 18, 4) = 42; + nlegtab( 3, 18, 4) = 35; + nlegtab( 4, 18, 4) = 32; + nlegtab( 5, 18, 4) = 31; + nlegtab( 1, 19, 4) = 67; + nlegtab( 2, 19, 4) = 42; + nlegtab( 3, 19, 4) = 35; + nlegtab( 4, 19, 4) = 33; + nlegtab( 5, 19, 4) = 31; + nlegtab( 1, 20, 4) = 67; + nlegtab( 2, 20, 4) = 43; + nlegtab( 3, 20, 4) = 36; + nlegtab( 4, 20, 4) = 33; + nlegtab( 5, 20, 4) = 32; + nlegtab( 1, 1, 5) = 85; + nlegtab( 2, 1, 5) = 47; + nlegtab( 3, 1, 5) = 34; + nlegtab( 4, 1, 5) = 28; + nlegtab( 5, 1, 5) = 24; + nlegtab( 1, 2, 5) = 85; + nlegtab( 2, 2, 5) = 47; + nlegtab( 3, 2, 5) = 34; + nlegtab( 4, 2, 5) = 28; + nlegtab( 5, 2, 5) = 24; + nlegtab( 1, 3, 5) = 85; + nlegtab( 2, 3, 5) = 47; + nlegtab( 3, 3, 5) = 35; + nlegtab( 4, 3, 5) = 29; + nlegtab( 5, 3, 5) = 25; + nlegtab( 1, 4, 5) = 85; + nlegtab( 2, 4, 5) = 47; + nlegtab( 3, 4, 5) = 35; + nlegtab( 4, 4, 5) = 29; + nlegtab( 5, 4, 5) = 26; + nlegtab( 1, 5, 5) = 85; + nlegtab( 2, 5, 5) = 47; + nlegtab( 3, 5, 5) = 36; + nlegtab( 4, 5, 5) = 30; + nlegtab( 5, 5, 5) = 26; + nlegtab( 1, 6, 5) = 85; + nlegtab( 2, 6, 5) = 48; + nlegtab( 3, 6, 5) = 36; + nlegtab( 4, 6, 5) = 30; + nlegtab( 5, 6, 5) = 27; + nlegtab( 1, 7, 5) = 85; + nlegtab( 2, 7, 5) = 48; + nlegtab( 3, 7, 5) = 37; + nlegtab( 4, 7, 5) = 31; + nlegtab( 5, 7, 5) = 28; + nlegtab( 1, 8, 5) = 85; + nlegtab( 2, 8, 5) = 48; + nlegtab( 3, 8, 5) = 37; + nlegtab( 4, 8, 5) = 32; + nlegtab( 5, 8, 5) = 28; + nlegtab( 1, 9, 5) = 85; + nlegtab( 2, 9, 5) = 49; + nlegtab( 3, 9, 5) = 38; + nlegtab( 4, 9, 5) = 32; + nlegtab( 5, 9, 5) = 29; + nlegtab( 1, 10, 5) = 85; + nlegtab( 2, 10, 5) = 49; + nlegtab( 3, 10, 5) = 38; + nlegtab( 4, 10, 5) = 33; + nlegtab( 5, 10, 5) = 30; + nlegtab( 1, 11, 5) = 85; + nlegtab( 2, 11, 5) = 49; + nlegtab( 3, 11, 5) = 39; + nlegtab( 4, 11, 5) = 33; + nlegtab( 5, 11, 5) = 30; + nlegtab( 1, 12, 5) = 85; + nlegtab( 2, 12, 5) = 50; + nlegtab( 3, 12, 5) = 39; + nlegtab( 4, 12, 5) = 34; + nlegtab( 5, 12, 5) = 31; + nlegtab( 1, 13, 5) = 85; + nlegtab( 2, 13, 5) = 50; + nlegtab( 3, 13, 5) = 39; + nlegtab( 4, 13, 5) = 35; + nlegtab( 5, 13, 5) = 32; + nlegtab( 1, 14, 5) = 85; + nlegtab( 2, 14, 5) = 50; + nlegtab( 3, 14, 5) = 40; + nlegtab( 4, 14, 5) = 35; + nlegtab( 5, 14, 5) = 32; + nlegtab( 1, 15, 5) = 85; + nlegtab( 2, 15, 5) = 51; + nlegtab( 3, 15, 5) = 40; + nlegtab( 4, 15, 5) = 36; + nlegtab( 5, 15, 5) = 33; + nlegtab( 1, 16, 5) = 85; + nlegtab( 2, 16, 5) = 51; + nlegtab( 3, 16, 5) = 41; + nlegtab( 4, 16, 5) = 36; + nlegtab( 5, 16, 5) = 34; + nlegtab( 1, 17, 5) = 85; + nlegtab( 2, 17, 5) = 52; + nlegtab( 3, 17, 5) = 41; + nlegtab( 4, 17, 5) = 37; + nlegtab( 5, 17, 5) = 34; + nlegtab( 1, 18, 5) = 85; + nlegtab( 2, 18, 5) = 52; + nlegtab( 3, 18, 5) = 42; + nlegtab( 4, 18, 5) = 37; + nlegtab( 5, 18, 5) = 35; + nlegtab( 1, 19, 5) = 85; + nlegtab( 2, 19, 5) = 52; + nlegtab( 3, 19, 5) = 42; + nlegtab( 4, 19, 5) = 38; + nlegtab( 5, 19, 5) = 36; + nlegtab( 1, 20, 5) = 85; + nlegtab( 2, 20, 5) = 53; + nlegtab( 3, 20, 5) = 43; + nlegtab( 4, 20, 5) = 39; + nlegtab( 5, 20, 5) = 37; + +end 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 b68418f..54c58dc 100644 --- a/chunkie/@chunker/chunker.m +++ b/chunkie/@chunker/chunker.m @@ -1,8 +1,70 @@ classdef chunker -%CHUNKER class which describes a curve divided into chunks (panels). +%CHUNKER class which describes a curve divided into chunks (or "panels"). +% % On each chunk the curve is represented by the values of its position, -% first and second derivatives in parameter space on a Legendre grid. +% first and second derivatives by scaled Legendre nodes. +% +% chunker properties: +% k - integer, number of Legendre nodes on each chunk +% nch - integer, number of chunks that make up the curve +% dim - integer, dimension of the ambient space in which the curve is +% embedded +% 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 +% 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 +% coordinate at the ith node on the jth chunk of the chunker +% n - dim x k x nch array of normals to the curve +% data - datadim x k x nch array of data attached to chunker points +% this data will be refined along with the chunker +% adj - 2 x nch integer array. adj(1,j) is i.d. of the chunk that +% precedes chunk j in parameter space. adj(2,j) is the i.d. of the +% chunk which follows. if adj(i,j) is 0 then that end of the chunk +% is a free end. if adj(i,j) is negative than that end of the chunk +% meets with other chunk ends in a vertex. the specific negative +% number acts as an i.d. of that vertex. % +% chunker methods: +% chunker(p,t,w) - construct an empty chunker with given preferences and +% precomputed Legendre nodes/weights (optional) +% obj = addchunk(obj,nchadd) - add nchadd chunks to the structure +% (initialized with zeros) +% obj = makedatarows(obj,nrows) - add nrows rows to the data storage. +% [obj,info] = sort(obj) - sort the chunks so that adjacent chunks are +% stored sequentially +% [rn,dn,d2n,dist,tn,ichn] = nearest(obj,ref,ich,opts,u,xover,aover) - +% find nearest point on chunker to ref +% obj = reverse(obj) - reverse chunk orientation +% rmin = min(obj) - minimum of coordinate values +% rmax = max(obj) - maximum of coordinate values +% wts = weights(obj) - scaled integration weights on curve +% obj.n = normals(obj) - recompute normal vectors +% onesmat = onesmat(obj) - matrix that corresponds to integration of a +% density on the curve +% rnormonesmat = normonesmat(obj) - matrix that corresponds to +% integration of dot product of normal with vector density on +% the curve +% plot(obj,varargin) - plot the chunker curve +% plot3(obj,idata,varargin) - 3D plot of the curve and one row of the +% data storage +% quiver(obj,varargin) - quiver plot of the chnkr points and normals +% scatter(obj,varargin) - scatter plot of the chnkr nodes +% tau = taus(obj) - unit tangents to curve +% obj = refine(obj,varargin) - refine the curve +% a = area(obj) - for a closed curve, area inside +% s = arclength(obj) - get values of arclength along curve +% kappa = signed_curvature(obj) - get signed curvature along curve +% rflag = datares(obj,opts) - check if data in data rows is resolved +% [rc,dc,d2c] = exps(obj) - get expansion coefficients for r,d,d2 +% ier = checkadjinfo(obj) - checks that the adjacency info of the curve +% is consistent +% [inds,adjs,info] = sortinfo(obj) - attempts to sort the curve and finds +% number of connected components, etc +% [re,taue] = chunkends(obj,ich) - get the endpoints of chunks +% flag = flagnear(obj,pts,opts) - flag points near the boundary + % author: Travis Askham (askhamwhat@gmail.com) @@ -14,7 +76,6 @@ d d2 adj - h n wts data @@ -24,7 +85,6 @@ dstor d2stor adjstor - hstor nstor wtsstor datastor @@ -49,6 +109,16 @@ methods function obj = chunker(p,t,w) + %CHUNKER construct an empty chunker with some default settings + % + % syntax: chnkr = chunker(p,t,w); + % + % optional input: + % p - struct of preferences + % p.k - integer, order to be used on chunks (16) + % p.dim - integer, dimension of ambient space (2) + % t, w - arrays of Legendre nodes and weights of order p.k + % if nargin < 1 || isempty(p) p = chunkerpref(); else @@ -74,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 = []; @@ -93,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 @@ -126,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 @@ -178,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); @@ -186,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; @@ -194,7 +255,6 @@ obj.wts = wtstemp; obj.d2 = d2temp; obj.adj = adjtemp; - obj.h = htemp; obj.data = datatemp; obj.nchstor = nchstornew; end @@ -273,6 +333,7 @@ obj.datastor = []; end + [obj2,f2] = upsample(obj,kup,f) [obj,info] = sort(obj) [rn,dn,d2n,dist,tn,ichn] = nearest(obj,ref,ich,opts,u,xover,aover) obj = reverse(obj) diff --git a/chunkie/@chunker/flagnear_rectangle_grid.m b/chunkie/@chunker/flagnear_rectangle_grid.m index fa60da3..79ea3cb 100644 --- a/chunkie/@chunker/flagnear_rectangle_grid.m +++ b/chunkie/@chunker/flagnear_rectangle_grid.m @@ -74,11 +74,12 @@ [xx,yy] = meshgrid(x(:),y(:)); nx = length(x(:)); ny = length(y(:)); +% bin sort structure for x and y vectors xmin = min(x(:)); xmax = max(x(:)); -nregx = max(round(length(x)/2),1); hx = (xmax-xmin)/nregx; +nregx = max(length(x)-1,1); hx = (xmax-xmin)/nregx; xids = round((x(:)-xmin)/hx); ymin = min(y(:)); ymax = max(y(:)); -nregy = max(round(length(y)/2),1); hy = (ymax-ymin)/nregy; +nregy = max(length(y)-1,1); hy = (ymax-ymin)/nregy; yids = round((y(:)-ymin)/hy); [xidsort,ix] = sort(xids,'ascend'); @@ -114,6 +115,11 @@ iymin = round((rymin-ymin)/hy)+1; iymax = round((rymax-ymin)/hy)+1; + ixmin = max(ixmin,1); + ixmax = min(ixmax,length(xiladr)-1); + iymin = max(iymin,1); + iymax = min(iymax,length(yiladr)-1); + ixrel = ix(xiladr(ixmin):(xiladr(ixmax+1)-1)); iyrel = iy(yiladr(iymin):(yiladr(iymax+1)-1)); diff --git a/chunkie/@chunker/max.m b/chunkie/@chunker/max.m index 39b4c42..b62fa85 100644 --- a/chunkie/@chunker/max.m +++ b/chunkie/@chunker/max.m @@ -17,4 +17,4 @@ % author: Travis Askham (askhamwhat@gmail.com) -rmax = max(reshape(obj.r,obj.dim,obj.k*obj.nch),[],2); \ No newline at end of file +rmax = max(real(reshape(obj.r,obj.dim,obj.k*obj.nch)),[],2); \ No newline at end of file diff --git a/chunkie/@chunker/merge.m b/chunkie/@chunker/merge.m index 2b71a92..052e3fb 100644 --- a/chunkie/@chunker/merge.m +++ b/chunkie/@chunker/merge.m @@ -1,22 +1,35 @@ -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 - -pref = []; +% mandatory setting pref.k = chnkrs(1).k; chnkrout = chunker(pref); -%chnkrout = chnkrs(1); - 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; @@ -29,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/min.m b/chunkie/@chunker/min.m index cf8648d..6ed706a 100644 --- a/chunkie/@chunker/min.m +++ b/chunkie/@chunker/min.m @@ -17,4 +17,4 @@ % author: Travis Askham (askhamwhat@gmail.com) -rmin = min(reshape(obj.r,obj.dim,obj.k*obj.nch),[],2); \ No newline at end of file +rmin = min(real(reshape(obj.r,obj.dim,obj.k*obj.nch)),[],2); \ No newline at end of file diff --git a/chunkie/@chunker/refine.m b/chunkie/@chunker/refine.m index 03bbafb..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 @@ -150,9 +144,6 @@ chunklens(i) = sum(chnkr.wts(:,i)); chunklens(nch) = sum(chnkr.wts(:,nch)); - - chunklens(i) = dot(dsdti,w); - chunklens(nch) = dot(dsdtnch,w); ifdone=0; @@ -168,7 +159,7 @@ % level restriction -if (strcmpi(lvlr,'a') || strcmpi(lvlr,'t')) +if (strcmpi(lvlr,'a')) for ijk = 1:maxiter_lvlr nchold=chnkr.nch; @@ -178,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 new file mode 100644 index 0000000..25d6f77 --- /dev/null +++ b/chunkie/@chunker/upsample.m @@ -0,0 +1,59 @@ +function [chnkrup,sigmaup] = upsample(chnkr,kup,sigma) +%UPSAMPLE return a chunker object which has been upsampled to a finer grid +% +% [chnkrup,sigmaup] = upsample(chnkr,kup,sigma); +% +% input: +% chnkr - chunker object +% kup - integer, order of panels for upsampled chunker +% sigma - optional, function on chunker to upsample +% +% output: +% +% chnkrup - upsampled verson of chunker object +% sigmaup - upsampled function (if sigma provided), empty otherwise +% + +k = chnkr.k; +dim = chnkr.dim; +nch = chnkr.nch; + +[~,~,u] = lege.exps(k); +[tu,wu,~,vu] = lege.exps(kup); + +upmat = vu(:,1:k)*u; + +pref = []; pref.k = kup; +chnkrup = chunker(pref,tu,wu); +chnkrup = chnkrup.addchunk(nch); + +nn = dim*nch; + +chnkrup.r = permute(reshape(upmat*reshape( ... + permute(chnkr.r,[2,1,3]),k,nn),kup,dim,nch),[2,1,3]); +chnkrup.d= permute(reshape(upmat*reshape( ... + permute(chnkr.d,[2,1,3]),k,nn),kup,dim,nch),[2,1,3]); +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.n = normals(chnkrup); +chnkrup.wts = weights(chnkrup); + +if chnkr.hasdata + ndata = size(chnkr.data,1); + nndata = ndata*nch; + chnkrup = chnkrup.makedatarows(ndata); + chnkrup.data = permute(reshape(upmat*reshape( ... + permute(chnkr.data,[2,1,3]),k,nndata),kup,ndata,nch),[2,1,3]); +end + +if nargin > 2 + dimsig = numel(sigma)/(nch*k); + sigma= reshape(sigma,dimsig,k,nch); + nnsig = dimsig*nch; + sigmaup = permute(reshape(upmat*reshape( ... + permute(sigma,[2,1,3]),k,nnsig),kup,dimsig,nch),[2,1,3]); +else + sigmaup = []; +end diff --git a/chunkie/@chunker/weights.m b/chunkie/@chunker/weights.m index 7de7a0e..fcce774 100644 --- a/chunkie/@chunker/weights.m +++ b/chunkie/@chunker/weights.m @@ -1,9 +1,10 @@ function wts = weights(chnkr) -%WEIGHTS integration weights suitable for smooth functions defined on the -% chunker object. Note that this routine is only to be used -% inside chunkerfunc and not accessed externally. +%WEIGHTS compute integration weights suitable for smooth functions defined +% on the chunker object. This routine is intended for use by chunker object +% constructors. It is not intended to be called by users, because +% such weights should already be stored in the chunker object. % -% This is merely the standard Legendre weights scaled to the chunks +% This is the standard Legendre weights scaled to the chunks % % Syntax: wts = weights(chnkr) % @@ -23,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 22fb4d4..111e7f2 100644 --- a/chunkie/@chunkgraph/balance.m +++ b/chunkie/@chunkgraph/balance.m @@ -25,7 +25,7 @@ echnks = obj.echnks; nverts = size(obj.verts,2); - + % Loop over all vertices in the cgraph structure for iii=1:nverts vedge = vstruc{iii}{1}; @@ -33,13 +33,12 @@ ifdone = false; while (~ifdone) - %%%%%%% find the arclengths of each panel incident to the edge + % find the arclengths of each panel on the edge incident to + % the vertex parcl = zeros([numel(vedge),1]); pinds = zeros([numel(vedge),1]); - [xleg16,wleg16,~,~] = lege.exps(16); - for ii=1:numel(vedge) if (sign(vsign(ii)) == -1) ds = obj.echnks(vedge(ii)).d(:,:,1); @@ -49,34 +48,39 @@ pinds(ii) = size(obj.echnks(vedge(ii)).r,3); end - h = obj.echnks(vedge(ii)).h(pinds(ii)); k = obj.echnks(vedge(ii)).k; - if (k ~=16) - [xleg,wleg,~,~] = lege.exps(k); - else - wleg = wleg16; - end - arc = sum(sqrt(ds(1,:).^2+ds(2,:).^2).*wleg'*h); + wleg = echnks(vedge(ii)).wstor; + + arc = sum(sqrt(ds(1,:).^2+ds(2,:).^2).*wleg'); parcl(ii) = arc; end - + % find the ratio of the max to min panels of all edges + % incident at the vertex and refine the panel with the + % max panel length amin = min(parcl); [amax,ind] = max(parcl); ind = ind(1); n2ref = floor(log(amax/amin)/log(2)); + if (n2ref >0) chnkr = obj.echnks(vedge(ind)); opts = []; opts.lvlrfac = 'a'; for ilev = 1:n2ref opts.splitchunks=([pinds(ind)]); + % If it is the last chunk, then update pinds to + % ensure that the last chunk is refined at + % subsequent levels if (pinds(ind)>1) pinds(ind) = pinds(ind)+1; end + % resort chnkr in case both end points need + % to be refined in the wrong order chnkr = refine(chnkr,opts); + chnkr = sort(chnkr); end obj.echnks(vedge(ind)) = chnkr; diff --git a/chunkie/@chunkgraph/build_v2emat.m b/chunkie/@chunkgraph/build_v2emat.m new file mode 100644 index 0000000..0357754 --- /dev/null +++ b/chunkie/@chunkgraph/build_v2emat.m @@ -0,0 +1,42 @@ +function A = build_v2emat(obj) +%V2EMAT get a "vertex to edge" matrix for the chunkgraph object +% +% input: obj - a chunkgraph object +% +% output: A - a nedges x nverts sparse matrix +% A(i,j) = 1 if edge i starts at vertex j +% A(i,j) = -1 if edge i ends at vertex j +% A(i,j) = 2 if edge i starts and ends at vertex j +% + +nverts = size(obj.verts,2); +nedges = size(obj.edgesendverts,2); + + +ii = zeros(2*nedges,1); +jj = zeros(2*nedges,1); +vv = zeros(2*nedges,1); +nf = 0; + +for i = 1:nedges + j1 = obj.edgesendverts(1,i); + j2 = obj.edgesendverts(2,i); + if j1 == j2 + nf = nf + 1; + ii(nf) = i; + jj(nf) = j1; + vv(nf) = 2; + else + nf = nf + 1; + ii(nf) = i; + jj(nf) = j1; + vv(nf) = -1; + nf = nf + 1; + ii(nf) = i; + jj(nf) = j2; + vv(nf) = 1; + end +end + +A = sparse(ii(1:nf),jj(1:nf),vv(1:nf),nedges,nverts); + diff --git a/chunkie/@chunkgraph/chunkgraph.m b/chunkie/@chunkgraph/chunkgraph.m index 7f83641..ddceffe 100644 --- a/chunkie/@chunkgraph/chunkgraph.m +++ b/chunkie/@chunkgraph/chunkgraph.m @@ -1,17 +1,25 @@ classdef chunkgraph %CHUNKGRAPH chunk graph class for storing complex domains % -% We describe a complex domain by its edges (smooth, regular -% curves) and vertices (free ends of edges or points where the ends of -% one or more edges meets) +% We describe a complex domain by a planar graph. Smooth boundary components +% (discretized by chunkers) specify the edges of the graph and the +% coordinates of free ends of the edges or points where the ends of multiple +% edges meet specify the vertices. The interiors of distinct edges are not +% allowed to intersect (i.e. intersecting curves should be split at the +% point where they intersect. +% +% The chunkgraph is specified by an array of chunkers called echnks (the +% edges), an array of points called verts (the vertices), and an array of +% indices specifying the vertices at the left and right ends of any % % properties(SetAccess=public) verts - edge2verts + edgesendverts echnks regions vstruc + v2emat end properties(SetAccess=public) @@ -26,16 +34,36 @@ end methods - function obj = chunkgraph(verts, edge2verts, fchnks, cparams, pref) + function obj = chunkgraph(verts, edgesendverts, fchnks, cparams) if (nargin == 0) return end if (numel(verts)==0) return end - obj.verts = verts; - obj.edge2verts = edge2verts; + obj.verts = verts; + + nverts = size(verts(:,:),2); + assert(nverts == size(edgesendverts,2),'edge specification not compatible with number of vertices'); + if (size(edgesendverts,1) ~= 2) + nedge = size(edgesendverts,1); + edgevertends_new = zeros(2,nedge); + nedge = size(edgesendverts,1); + for i = 1:nedge + edgevertends_new(1,i) = find(edgesendverts(i,:) == -1); + edgevertends_new(2,i) = find(edgesendverts(i,:) == 1); + end + edgesendverts = edgevertends_new; + assert(all(edgesendverts(:) ~= 0),'edge specification had an error'); + end + + obj.edgesendverts = edgesendverts; + obj.v2emat = build_v2emat(obj); obj.echnks = chunker.empty; + + if nargin < 3 + fchnks = []; + end if (nargin < 4) cploc = []; @@ -73,15 +101,11 @@ end - if (size(verts,2) ~= size(edge2verts,2)) - error('Incompatible vertex and edge sizes'); - end - echnks = chunker.empty(); - for i=1:size(edge2verts,1) + for i=1:size(edgesendverts,2) if (numel(fchnks)1) + rgn2 = rgns{2}; + [rgnout] = mergeregions(obj,rgnout,rgn2); + for ii=3:numel(rgns) + rgn2 = rgns{ii}; + [rgnout] = mergeregions(obj,rgnout,rgn2); + end + end + + regions = rgnout; + obj.regions = regions; end function obj = set.verts(obj,val) @@ -137,10 +228,10 @@ validateattributes(val,classes,{}) obj.verts = val; end - function obj = set.edge2verts(obj,val) + function obj = set.edgesendverts(obj,val) classes = {'numeric'}; validateattributes(val,classes,{}) - obj.edge2verts = val; + obj.edgesendverts = val; end function obj = set.echnks(obj,val) classes = {'chunker'}; @@ -204,7 +295,11 @@ sourceinfo.d2= d2s; sourceinfo.w = ws; end - end + + % defined in other files + spmat = build_v2emat(obj) + end + methods(Static) end diff --git a/chunkie/@chunkgraph/findregions.m b/chunkie/@chunkgraph/findregions.m index 5f39efc..4c87942 100644 --- a/chunkie/@chunkgraph/findregions.m +++ b/chunkie/@chunkgraph/findregions.m @@ -27,18 +27,17 @@ else vstruc = procverts(obj); end -nedge = size(obj.edge2verts,1); -e2v = obj.edge2verts; - -% each edge belongs to two regions (going in opposite directions) -edges = [1:nedge,-(1:nedge)]; +nedge = size(obj.edgesendverts,2); % if the user has provided iverts, get the reduced edge list if (nargin >1) - e2vtmp = e2v(:,iverts); - [iinds,jinds] = find(e2vtmp ~= 0); + v2etmp = obj.v2emat(:,iverts); + [iinds,jinds] = find(v2etmp ~= 0); iinds = unique(iinds); edges = [iinds,-iinds]; +else + % each edge belongs to two regions (going in opposite directions) + edges = [1:nedge,-(1:nedge)]; end regions = {}; @@ -56,8 +55,14 @@ estart = enum; ecycle = [enum]; - iv0 = find(e2v(abs(enum),:)==-sign(enum)); - ivc = find(e2v(abs(enum),:)== sign(enum)); + if enum > 0 + iv0 = obj.edgesendverts(1,abs(enum)); + ivc = obj.edgesendverts(2,abs(enum)); + else + iv0 = obj.edgesendverts(2,abs(enum)); + ivc = obj.edgesendverts(1,abs(enum)); + end + ifdone = false; while (~ifdone) @@ -72,7 +77,11 @@ esign = vstruc{ivc}{2}(1); end enum = -esign*enext; - ivc = find(e2v(abs(enum),:)== sign(enum)); + if enum > 0 + ivc = obj.edgesendverts(2,abs(enum)); + else + ivc = obj.edgesendverts(1,abs(enum)); + end if (enum == estart) ifdone = true; else diff --git a/chunkie/@chunkgraph/vertextract.m b/chunkie/@chunkgraph/vertextract.m index 6c00a8a..7721b91 100644 --- a/chunkie/@chunkgraph/vertextract.m +++ b/chunkie/@chunkgraph/vertextract.m @@ -25,10 +25,14 @@ % author: Jeremy Hoskins +irel = find(cgrph.v2emat(:,ivert) ~= 0); % extract the indices of the edges which terminate at ivert. -ieplus = find(cgrph.edge2verts(:,ivert) == 1); +ieplus = find(cgrph.edgesendverts(2,irel) == ivert); +ieplus = irel(ieplus); ieplus = ieplus(:).'; % extract the indices of the edges which begin at ivert. -ieminus = find(cgrph.edge2verts(:,ivert) == -1); +ieminus = find(cgrph.edgesendverts(1,irel) == ivert); +ieminus = irel(ieminus); ieminus = ieminus(:).'; + % for each incoming edge, get the tangent vector near the end (at the % last discretization node) @@ -50,8 +54,8 @@ [angs,isrtededges] = sort(angs); % compute inds and isgn using isrtedges -inds = [ieplus',ieminus']; -isgn = [ones(size(ieplus')),-ones(size(ieminus'))]; +inds = [ieplus,ieminus]; +isgn = [ones(size(ieplus)),-ones(size(ieminus))]; inds = inds(isrtededges); isgn = isgn(isrtededges); end diff --git a/chunkie/@kernel/elast2d.m b/chunkie/@kernel/elast2d.m index 25a9143..3940faf 100644 --- a/chunkie/@kernel/elast2d.m +++ b/chunkie/@kernel/elast2d.m @@ -8,6 +8,10 @@ % MU) constructs the single-layer elasticity kernel for traction with % Lamé parameters LAM and MU. % +% KERNEL.ELAST2D('sgrad', LAM, MU) or KERNEL.ELAST2D('straction', LAM, +% MU) constructs the gradient of the single-layer elasticity kernel with +% Lamé parameters LAM and MU. +% % KERNEL.ELAST2D('d', LAM, MU) or KERNEL.ELAST2D('double', LAM, MU) % constructs the double-layer elasticity kernel with Lamé parameters LAM % and MU. @@ -15,7 +19,16 @@ % KERNEL.ELAST2D('dalt', LAM, MU) constructs the alternative double-layer % elasticity kernel with Lamé parameters LAM and MU. % -% See also CHNK.ELAST2D.KERN. +% KERNEL.ELAST2D('dalttrac', LAM, MU) or +% KERNEL.ELAST2D('dalttraction', LAM, MU) +% constructs the traction of the alternative double-layer elasticity +% kernel with Lamé parameters LAM and MU. +% +% KERNEL.ELAST2D('daltgrad', LAM, MU) constructs the gradient of the +% alternative double-layer elasticity kernel with Lamé parameters LAM +% and MU. +% +% See also CHNK.ELAST2D.KERN if ( nargin < 1 ) error('Missing elasticity kernel type.'); @@ -41,6 +54,11 @@ obj.eval = @(s,t) chnk.elast2d.kern(lam, mu, s, t, 's'); obj.sing = 'log'; + case {'sgrad'} + obj.type = 'sgrad'; + obj.eval = @(s,t) chnk.elast2d.kern(lam, mu, s, t, 'sgrad'); + obj.sing = 'pv'; + case {'strac', 'straction'} obj.type = 'strac'; obj.eval = @(s,t) chnk.elast2d.kern(lam, mu, s, t, 'strac'); @@ -56,6 +74,11 @@ obj.eval = @(s,t) chnk.elast2d.kern(lam, mu, s, t, 'dalt'); obj.sing = 'smooth'; + case {'dalttrac','dalttraction'} + obj.type = 'dalttrac'; + obj.eval = @(s,t) chnk.elast2d.kern(lam, mu, s, t, 'dalttrac'); + obj.sing = 'hs'; + case {'daltgrad'} obj.type = 'daltgrad'; obj.eval = @(s,t) chnk.elast2d.kern(lam, mu, s, t, 'daltgrad'); @@ -67,4 +90,9 @@ end +icheck = exist(['fmm2d.' mexext], 'file'); +if icheck ~=3 + obj.fmm = []; +end + end diff --git a/chunkie/@kernel/helm2d.m b/chunkie/@kernel/helm2d.m index 64bed00..057b5d3 100644 --- a/chunkie/@kernel/helm2d.m +++ b/chunkie/@kernel/helm2d.m @@ -14,6 +14,10 @@ % parameter ETA, i.e., COEFS(1)*KERNEL.HELM2D('d', ZK) + % COEFS(2)*KERNEL.HELM2D('s', ZK). % +% KERNEL.HELM2D('cp', ZK, COEFS) or KERNEL.HELM2D('combined', ZK, COEFS) +% constructs the derivative of the combined-layer Helmholtz kernel with +% wavenumber ZK and parameter ETA, i.e., COEFS(1)*KERNEL.HELM2D('dp', ZK) + +% COEFS(2)*KERNEL.HELM2D('sp', ZK). % See also CHNK.HELM2D.KERN. if ( nargin < 1 ) @@ -66,9 +70,26 @@ obj.fmm = @(eps,s,t,sigma) chnk.helm2d.fmm(eps, zk, s, t, 'c', sigma, coefs); obj.sing = 'log'; + case {'cp', 'cprime'} + if ( nargin < 3 ) + warning('Missing combined layer coefficients. Defaulting to [1,1i].'); + coefs = [1,1i]; + end + obj.type = 'cprime'; + obj.params.coefs = coefs; + obj.eval = @(s,t) chnk.helm2d.kern(zk, s, t, 'cprime', coefs); + obj.fmm = @(eps,s,t,sigma) chnk.helm2d.fmm(eps, zk, s, t, 'cprime', sigma, coefs); + obj.sing = 'hs'; + otherwise error('Unknown Helmholtz kernel type ''%s''.', type); end +icheck = exist(['fmm2d.' mexext], 'file'); +if icheck ~=3 + obj.fmm = []; +end + + end diff --git a/chunkie/@kernel/helm2ddiff.m b/chunkie/@kernel/helm2ddiff.m index a1392b9..9dae103 100644 --- a/chunkie/@kernel/helm2ddiff.m +++ b/chunkie/@kernel/helm2ddiff.m @@ -92,4 +92,10 @@ end +icheck = exist(['fmm2d.' mexext], 'file'); +if icheck ~=3 + obj.fmm = []; +end + + end diff --git a/chunkie/@kernel/kernel.m b/chunkie/@kernel/kernel.m index 0162e92..a25b09f 100644 --- a/chunkie/@kernel/kernel.m +++ b/chunkie/@kernel/kernel.m @@ -9,6 +9,7 @@ % ---- ---- % 'laplace' ('lap', 'l') 's', 'd', 'sp', 'c' % 'helmholtz' ('helm', 'h') 's', 'd', 'sp', 'dp', 'c' +% 'cp' % 'helmholtz difference' ('helmdiff', 'hdiff') 's', 'd', 'sp', 'dp' % 'elasticity' ('elast', 'e') 's', 'strac', 'd', 'dalt' % 'stokes' ('stok', 's') 'svel', 'spres', 'strac', diff --git a/chunkie/@kernel/lap2d.m b/chunkie/@kernel/lap2d.m index 8db13bb..1eb548a 100644 --- a/chunkie/@kernel/lap2d.m +++ b/chunkie/@kernel/lap2d.m @@ -85,4 +85,9 @@ end +icheck = exist(['fmm2d.' mexext], 'file'); +if icheck ~=3 + obj.fmm = []; +end + end diff --git a/chunkie/@kernel/stok2d.m b/chunkie/@kernel/stok2d.m index c8f2327..e017192 100644 --- a/chunkie/@kernel/stok2d.m +++ b/chunkie/@kernel/stok2d.m @@ -84,4 +84,10 @@ end +icheck = exist(['fmm2d.' mexext], 'file'); +if icheck ~=3 + obj.fmm = []; +end + + end diff --git a/chunkie/FLAM b/chunkie/FLAM index 2c9361e..a83698e 160000 --- a/chunkie/FLAM +++ b/chunkie/FLAM @@ -1 +1 @@ -Subproject commit 2c9361e7bc8fb8c3255f33275c0f2f7b7bbb7572 +Subproject commit a83698efdd58293c0e4c0c52f59498ffb9da9a99 diff --git a/chunkie/chunkerfit.m b/chunkie/chunkerfit.m new file mode 100644 index 0000000..987ec68 --- /dev/null +++ b/chunkie/chunkerfit.m @@ -0,0 +1,185 @@ +function chnkr = chunkerfit(xy, opts) +%CHUNKERFIT Create a chunker by fitting a curve to a set of points. +% +% Syntax: chnkr = CHUNKERFIT(xy, opts) +% +% Input: +% xy - (2,:) array of points used to fit a curve +% +% Optional input: +% opts - options structure (defaults) +% opts.method = string ('spline') +% Method to use for curve fitting. +% opts.ifclosed = boolean (true) +% Flag for fitting a closed (true) or open (false) curve. +% opts.splitatpoints = boolean (false) +% If true, start the adaptive splitting process from the set of +% splits implied by the given points. This ensures that the given +% points will always coincide with panel endpoints in the resulting +% chunker object. If false, the adaptive splitting process will +% start from scratch. +% opts.cparams = struct +% Curve parameters structure to be passed to chunkerfunc after +% fitting is done. See chunkerfunc for parameters. +% opts.pref = chunkerpref object or struct +% Chunker preferences to be passed to chunkerfunc. See chunkerpref +% for parameters. +% +% Output: +% chnkr - chunker object discretizing the fitted curve +% +% Examples: +% r = chnk.curves.bymode(sort(2*pi*rand(20,1)), [2 0.5 0.2 0.7]); +% chnkr = chunkerfit(r); % Closed curve +% opts = []; opts.ifclosed = false; +% chnkr = chunkerfit(r(:,1:10), opts); % Open curve +% +% See also CHUNKERFUNC, CHUNKERPREF, CHUNKER. + +if ( nargin == 0 ) + test(); + return +end + +if ( size(xy, 1) ~= 2 ) + error('CHUNKIE:CHUNKERFIT:size', 'Points must be specified as a 2xN matrix.'); +end + +x = xy(1,:).'; +y = xy(2,:).'; + +% Set default options +defaults = []; +defaults.method = 'spline'; +defaults.ifclosed = true; +defaults.splitatpoints = false; +defaults.cparams = []; +defaults.pref = []; + +if ( nargin < 2 ) + opts = defaults; +else + for field = fieldnames(defaults).' + if ( ~isfield(opts, field{1}) ) + opts.(field{1}) = defaults.(field{1}); + end + end +end + +if ( ~strcmpi(opts.method, 'spline') ) + error('CHUNKIE:CHUNKERFIT:method', 'Unsupported method ''%s''.', opts.method); +end + +if ( opts.ifclosed ) + % Wrap the nodes around + x = [x; x(1)]; + y = [y; y(1)]; +end + +% Choose the parametric grid to be polygonal arc length +t = [0; cumsum(sqrt(diff(x).^2 + diff(y).^2))]; + +ppx = myspline(t, x, opts.ifclosed); +ppy = myspline(t, y, opts.ifclosed); +ppx_dx = ppdiff(ppx, 1); +ppy_dy = ppdiff(ppy, 1); +ppx_dxx = ppdiff(ppx, 2); +ppy_dyy = ppdiff(ppy, 2); + +% Pack polynomial coefficients into a tensor to vectorize ppval +cfs = zeros([6 size(ppx.coefs)]); +cfs(1,:,:) = ppx.coefs; +cfs(2,:,:) = ppy.coefs; +cfs(3,:,:) = ppx_dx.coefs; +cfs(4,:,:) = ppy_dy.coefs; +cfs(5,:,:) = ppx_dxx.coefs; +cfs(6,:,:) = ppy_dyy.coefs; +breaks = ppx.breaks.'; + + function [r, d, d2] = splinefunc(tt) + vals = myppval(breaks, cfs, tt); + xs = vals(1,:); + ys = vals(2,:); + dxs = vals(3,:); + dys = vals(4,:); + d2xs = vals(5,:); + d2ys = vals(6,:); + r = [xs(:).' ; ys(:).' ]; + d = [dxs(:).' ; dys(:).' ]; + d2 = [d2xs(:).'; d2ys(:).']; + end + +cparams = opts.cparams; +cparams.ifclosed = opts.ifclosed; +cparams.ta = t(1); +cparams.tb = t(end); +if ( opts.splitatpoints ) + cparams.tsplits = t; +end +chnkr = chunkerfunc(@splinefunc, cparams, opts.pref); + +end + +function test() + +rng(0) +r = chnk.curves.bymode(sort(2*pi*rand(20,1)), [2 0.5 0.2 0.7]); +chnkr = chunkerfit(r); + +plot(chnkr, 'bo-') +hold on +plot(r(1,:), r(2,:), 'r.', markersize=30) +hold off +shg + +end + +function y = myppval(breaks, cfs, x) + [~, idx] = histc(x, [-inf; breaks(2:end-1); inf]); + p = reshape((x-breaks(idx)).^(3:-1:0), [1 size(idx,1) size(cfs,3)]); + y = sum(p .* cfs(:,idx,:), 3); +end + +function pp = ppdiff(pp, n) + deg = pp.order-1; + D = diag(deg:-1:1, 1); + for k = 1:n + pp.coefs = pp.coefs * D^n; + end +end + +function pp = myspline(x, y, ifclosed) + n = length(y); + dx = diff(x); + dy = diff(y); + dydx = dy ./ dx; + b = zeros(n, 1); + A = spdiags([ [dx(2:n-1) ; 0 ; 0], ... + 2*[0 ; dx(2:n-1)+dx(1:n-2) ; 0], ... + [0 ; 0 ; dx(1:n-2)] ], ... + [-1 0 1], n, n); + + if ( ifclosed ) + % Periodic conditions: match first and second derivatives at the first + % and last points + A(1,[1 n]) = [1 -1]; + A(n,1:2) = dx(n-1)*[2 1]; + A(n,n-1:n) = A(n,n-1:n) + dx(1)*[1 2]; + b(2:n-1,:) = 3*(dx(2:n-1).*dydx(1:n-2) + dx(1:n-2).*dydx(2:n-1)); + b(n,:) = 3*(dx(n-1)*dydx(1) + dx(1)*dydx(n-1)); + else + % Not-a-knot conditions: match third derivatives at the first and last + % interior breaks + A(1,1:2) = [dx(2), x(3)-x(1)]; + A(n,n-1:n) = [x(n)-x(n-2), dx(n-2)]; + b(1,:) = ((dx(1)+2*(x(3)-x(1)))*dx(2)*dydx(1) + dx(1)^2*dydx(2)) / (x(3)-x(1)); + b(2:n-1,:) = 3*(dx(2:n-1).*dydx(1:n-2) + dx(1:n-2).*dydx(2:n-1)); + b(n,:) = (dx(n-1)^2*dydx(n-2) + ((2*(x(n)-x(n-2))+dx(n-1))*dx(n-2))*dydx(n-1)) / (x(n)-x(n-2)); + end + + % Solve linear system for the coefficients + c = A \ b; + c4 = (c(1:n-1)+c(2:n)-2*dydx(1:n-1)) ./ dx; + c3 = (dydx(1:n-1)-c(1:n-1))./dx - c4; + pp = mkpp(x, [c4./dx c3 c(1:n-1) y(1:n-1)]); +end diff --git a/chunkie/chunkerflam.m b/chunkie/chunkerflam.m index fb1a1e5..9d4590d 100644 --- a/chunkie/chunkerflam.m +++ b/chunkie/chunkerflam.m @@ -244,7 +244,7 @@ mmax = max([mmax,max(chnkr)],[],2); mmin = min([mmin,min(chnkr)],[],2); end - +xflam = real(xflam); width = max(mmax-mmin); chnkrtotal = merge(chnkrs); diff --git a/chunkie/chunkerfunc.m b/chunkie/chunkerfunc.m index c3f30d2..7bdc048 100644 --- a/chunkie/chunkerfunc.m +++ b/chunkie/chunkerfunc.m @@ -13,21 +13,22 @@ % [r,d] = fcurve(t); or [r,d,d2] = fcurve(t); % where d is the first derivative of r with respect to t and % d2 is the second derivative. in some situations, this will -% improve the convergence order. +% improve the convergence order and final precision. % % Optional input: % cparams - curve parameters structure (defaults) % cparams.ta = left end of t interval (0) % cparams.tb = right end of t interval (2*pi) +% cparams.tsplits = set of initial break points for discretization in +% parameter space (should be in [ta,tb]) % cparams.ifclosed = flag determining if the curve % is to be interpreted as a closed curve (true) % cparams.chsmall = max size of end intervals if % ifclosed == 0 (Inf) % cparams.nover = oversample resolved curve nover % times (0) -% cparams.eps = resolve coordinates, arclength, -% and first and second derivs of coordinates -% to this tolerance (1.0e-6) +% cparams.eps = tolerance to resolve coordinates and arclength +% density (1.0e-6) % cparams.lvlr = string, determines type of level % restriction to be enforced % lvlr = 'a' -> no chunk should have double the arc length @@ -42,6 +43,9 @@ % pref.nchmax - maximum number of chunks (10000) % pref.k - number of Legendre nodes on chunks (16) % +% Output: +% chnkr - a chunker object containing the discretization of the domain +% % Examples: % chnkr = chunkerfunc(@(t) starfish(t)); % chunk up starfish w/ standard % % options @@ -137,19 +141,42 @@ ab = zeros(2,nchmax); adjs = zeros(2,nchmax); -ab(1,1)=ta; -ab(2,1)=tb; -nch=1; + +if (isfield(cparams,'tsplits')) + tsplits = cparams.tsplits; + tsplits = [tsplits(:); ta; tb]; +else + tsplits = [ta;tb]; +end + +tsplits = sort(unique(tsplits),'ascend'); +lab = length(tsplits); +if (lab-1 > nchmax) + error(['nchmax exceeded in chunkerfunc on initial splits.\n ',... + 'try increasing nchmax in preference struct']); +end +if (any(tsplits > tb) || any(tsplits < ta)) + error(['tsplits outside of interval of definition.\n', ... + 'check definition of splits, ta and tb']); +end + +ab(1,1:(lab-1)) = tsplits(1:end-1); +ab(2,1:(lab-1)) = tsplits(2:end); + +nch=lab-1; +adjs(1,1:nch) = 0:(nch-1); +adjs(2,1:nch) = 2:(nch+1); + if ifclosed - adjs(1,1)=1; - adjs(2,1)=1; + adjs(1,1)=nch; + adjs(2,nch)=1; else adjs(1,1)=-1; - adjs(2,1)=-1; + adjs(2,nch)=-1; end nchnew=nch; -maxiter_res=nchmax; +maxiter_res=nchmax-nch; rad_curr = 0; for ijk = 1:maxiter_res @@ -493,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/chunkerinterior.m b/chunkie/chunkerinterior.m index 4468bec..4a23291 100644 --- a/chunkie/chunkerinterior.m +++ b/chunkie/chunkerinterior.m @@ -1,12 +1,18 @@ -function in = chunkerinterior(chnkr,pts,opts) -%CHUNKERINTERIOR returns an array indicating whether each point in -% pts is inside the domain. Assumes the domain is closed. +function [in] = chunkerinterior(chnkobj,ptsobj,opts) +%CHUNKERINTERIOR returns an array indicating whether each point specified +% % by pts is inside the domain. Assumes the domain is closed. % -% Syntax: in = chunkerinterior(chnkr,pts,opts) +% Syntax: in = chunkerinterior(chnkobj,pts,opts) +% in = chunkerinterior(chnkobj,{x,y},opts) % meshgrid version % % Input: -% chnkr - chunker object describing geometry -% pts - (chnkr.dim,:) array of points to test +% chnkobj - chunker object or chunkgraph object describing geometry +% ptsobj - object describing the target points, can be specified as +% * (chnkr.dim,:) array of points to test +% * {x,y} - length 2 cell array. the points checked then have the +% coordinates of a mesh grid [xx,yy] = meshgrid(x,y) +% * chunker object, in which case it uses chunker.r(:,:) +% * chunkgraph object, in which case it uses chunkgraph.r(:,:) % % Optional input: % opts - options structure with entries: @@ -19,7 +25,7 @@ % % Output: % in - logical array, if in(i) is true, then pts(:,i) is inside the -% domain +% domain or for a mesh grid [xx(i); yy(i)] is inside the domain. % % Examples: % chnkr = chunkerfunc(@(t) starfish(t)); @@ -29,16 +35,41 @@ % author: Travis Askham (askhamwhat@gmail.com) -if class(chnkr) == "chunkgraph" - chnkr = merge(chnkr.echnks); +grid = false; + +if nargin < 3 + opts = []; +end + +% Assign appropriate object to chnkr +if class(chnkobj) == "chunker" + chnkr = chnkobj; +elseif class(chnkobj) == "chunkgraph" + chnkr = merge(chnkobj.echnks); +else + msg = "Unsupported object in chunkerinterior"; + error(msg) end assert(chnkr.dim == 2,'interior only well-defined for 2D'); -if nargin < 3 - opts = []; +% Assign appropriate object to pts +if isa(ptsobj, "cell") + assert(length(ptsobj)==2,'second input should be either 2xnpts array or length 2 cell array'); + x = ptsobj{1}; + y = ptsobj{2}; + grid = true; + [xx,yy] = meshgrid(x,y); + pts = [xx(:).'; yy(:).']; +elseif isa(ptsobj, "chunker") + pts = ptsobj.r(:,:); +elseif isa(ptsobj, "chunkgraph") + pts = ptsobj.r(:,:); +else + pts = ptsobj; end + usefmm = true; if isfield(opts,'fmm') usefmm = opts.fmm; @@ -92,17 +123,34 @@ useflam_final = useflam; end +% use bernstein ellipses and rectangles to flag problematic points +% + +eps_local = 1e-3; + +rho = 1.2; +optsflag = []; optsflag.rho = rho; optsflag.occ = 5; +if grid + flag = flagnear_rectangle_grid(chnkr,x,y,optsflag); +else + flag = flagnear_rectangle(chnkr,pts,optsflag); +end + +npoly = chnkr.k; +nlegnew = chnk.ellipse_oversample(rho,npoly,eps_local); +nlegnew = max(nlegnew,chnkr.k); + +[chnkr2] = upsample(chnkr,nlegnew); icont = false; if usefmm_final try - eps_local = 1e-3; - wchnkr = chnkr.wts; - dens1_fmm = ones(chnkr.k*chnkr.nch,1).*wchnkr(:); + wchnkr = chnkr2.wts; + dens1_fmm = ones(chnkr2.k*chnkr2.nch,1).*wchnkr(:); pgt = 1; - vals1 = chnk.lap2d.fmm(eps_local,chnkr,pts,'d',dens1_fmm,pgt); + vals1 = chnk.lap2d.fmm(eps_local,chnkr2,pts,'d',dens1_fmm,pgt); catch fprintf('using fmm failed due to incompatible mex, try regenrating mex\n'); useflam_final = useflam; @@ -112,102 +160,82 @@ if ~usefmm_final || icont kernd = kernel('lap','d'); - dens1 = ones(chnkr.k,chnkr.nch); - wts = chnkr.wts; + dens1 = ones(chnkr2.k,chnkr2.nch); + opdims = [1 1]; if useflam_final xflam1 = chnkr.r(:,:); - matfun = @(i,j) chnk.flam.kernbyindexr(i,j,pts,chnkr,kernd,opdims); + matfun = @(i,j) chnk.flam.kernbyindexr(i,j,pts,chnkr2,kernd,opdims); [pr,ptau,pw,pin] = chnk.flam.proxy_square_pts(); pxyfun = @(rc,rx,cx,slf,nbr,l,ctr) chnk.flam.proxyfunr(rc,rx,slf,nbr,l, ... - ctr,chnkr,kernd,opdims,pr,ptau,pw,pin); - F = ifmm(matfun,pts,xflam1,200,1e-14,pxyfun); + ctr,chnkr2,kernd,opdims,pr,ptau,pw,pin); + F = ifmm(matfun,pts,xflam1,200,1e-6,pxyfun); vals1 = ifmm_mv(F,dens1(:),matfun); else optskerneval = []; optskerneval.usesmooth = 1; - vals1 = chunkerkerneval(chnkr,kernd,dens1,pts,optskerneval); + vals1 = chunkerkerneval(chnkr2,kernd,dens1,pts,optskerneval); end end in = abs(vals1+1) < abs(vals1); -% find nearest neighbors at certain level of refinement (here chosen -% uniformly, for simplicity) - -nt = size(pts,2); -xx = zeros(2,chnkr.npt + nt); -xx(:,1:chnkr.npt) = chnkr.r(:,:); -xx(:,chnkr.npt+1:chnkr.npt+nt) = pts; - -pt2chnk = repmat(1:chnkr.nch,chnkr.k,1); - -chunklens = zeros(chnkr.nch,1); -ws = chnkr.wts; -chunklens(:) = sum(ws,1); -lmax = max(chunklens)*3/chnkr.k; - -T = hypoct_uni(xx,lmax); -targ_dists = Inf(nt,1); -itarg_dists = (1:nt).'; - -rnorm = normals(chnkr); -normdist_flag = false(size(targ_dists)); - -for i = 1:length(T.nodes) - xi = T.nodes(i).xi; - if nnz(xi <= chnkr.npt) > 0 - isrc = xi(xi <= chnkr.npt); - inbor = [T.nodes(T.nodes(i).nbor).xi T.nodes(i).xi]; - if nnz(inbor > chnkr.npt) > 0 - - % get all pairwise distances - itarg = inbor(inbor > chnkr.npt)-chnkr.npt; - srcs = reshape(chnkr.r(:,isrc),chnkr.dim,1,length(isrc)); - targs = reshape(pts(:,itarg),chnkr.dim,length(itarg),1); - dists = reshape(sqrt(sum((bsxfun(@minus,targs,srcs)).^2,1)), ... - length(itarg),length(isrc)); - [mindists,inds] = min(dists,[],2); - isrc2 = isrc(inds); - - % if minimum distance for any targ is smaller than - % previously seen, we update - ifnew = mindists < targ_dists(itarg); - targ_dists(itarg(ifnew)) = mindists(ifnew); - itarg_dists(itarg(ifnew)) = isrc2(ifnew); - - % if vector from point to boundary point aligns - % with normal, probably on inside. - % if angle is near right angle or distance is small, - % flag it for more precise testing - rdiff = chnkr.r(:,isrc2(ifnew)) - pts(:,itarg(ifnew)); - rnorms = rnorm(:,isrc2(ifnew)); - lens = chunklens(pt2chnk(isrc2(ifnew))); - - rdrn = sum(rdiff.*rnorms,1); rdrn = rdrn(:); - lens = lens(:); - in(itarg(ifnew)) = rdrn > 0; - normdist_flag(itarg(ifnew)) = abs(rdrn) < 0.1*lens; +% for points where the integral might be inaccurate: +% find close boundary point and check normal direction + + +nnzpt = sum(flag~=0,2); +ipt = find(nnzpt); + +npts = numel(pts)/2; +distmins = inf(npts,1); +dss = zeros(2,npts); +rss = zeros(2,npts); + +k = chnkr.k; +[t,~,u] = lege.exps(k); + +for i = 1:chnkr.nch + + % check side based on closest boundary node + rval = chnkr.r(:,:,i); + dval = chnkr.d(:,:,i); + nval = chnkr.n(:,:,i); + [ji] = find(flag(:,i)); + ptsi = pts(:,ji); + nptsi = size(ptsi,2); + dist2all = reshape(sum( abs(reshape(ptsi,2,1,nptsi) ... + - reshape(rval,2,k,1)).^2, 1),k,nptsi); + [dist2all,ipti] = min(dist2all,[],1); + for j = 1:length(ji) + jj = ji(j); + if dist2all(j) < distmins(jj) + distmins(jj) = dist2all(j); + rss(:,jj) = rval(:,ipti(j)); + dss(:,jj) = dval(:,ipti(j)); end end -end -% more precise testing - -iflagged = find(normdist_flag); -[~,~,u] = lege.exps(chnkr.k); -tover = lege.exps(2*chnkr.k); -ainterpover = lege.matrin(chnkr.k,tover); -for i = 1:length(iflagged) - ii = iflagged(i); - ich = pt2chnk(itarg_dists(ii)); - ich = [ich; chnkr.adj(:,ich)]; - [rn,dn] = nearest(chnkr,pts(:,ii),ich,[],u,tover,ainterpover); - in(ii) = (rn(:)-pts(:,ii)).'*[dn(2);-dn(1)] > 0; + % if angle is small do a refined search for closest point + ptsn = nval(:,ipti); + diffs = rval(:,ipti)-ptsi; + dots = sum(ptsn.*diffs,1); + + jsus = abs(dots) < 2e-1*sqrt(dist2all); + jii = ji(jsus); + [~,rs,ds,~,dist2s] = chnk.chunk_nearparam(rval,pts(:,jii),[],t,u); + for j = 1:length(jii) + jj = jii(j); + if dist2s(j) < distmins(jj) + distmins(jj) = dist2s(j); + rss(:,jj) = rs(:,j); + dss(:,jj) = ds(:,j); + end + end end - - - +for i = 1:length(ipt) + jj = ipt(i); + in(jj) = (rss(:,jj)-pts(:,jj)).'*[dss(2,jj);-dss(1,jj)] > 0; end diff --git a/chunkie/chunkerkerneval.m b/chunkie/chunkerkerneval.m index 13f6197..cad3a9e 100644 --- a/chunkie/chunkerkerneval.m +++ b/chunkie/chunkerkerneval.m @@ -1,17 +1,20 @@ -function fints = chunkerkerneval(chnkr,kern,dens,targs,opts) +function fints = chunkerkerneval(chnkobj,kern,dens,targobj,opts) %CHUNKERKERNEVAL compute the convolution of the integral kernel with % the density defined on the chunk geometry. % % Syntax: fints = chunkerkerneval(chnkr,kern,dens,targs,opts) % % Input: -% chnkr - chunker object description of curve +% chnkobj - chunker object or chunkgraph object description of curve % kern - kernel class object or kernel function taking inputs % kern(srcinfo,targinfo) % dens - density on boundary, should have size opdims(2) x k x nch % where k = chnkr.k, nch = chnkr.nch, where opdims is the % size of kern for a single src,targ pair -% targs - targ(1:2,i) gives the coords of the ith target +% targobj - object describing the target points, can be specified as +% * array of points +% * chunker object +% * chunkgraph object % % Optional input: % opts - structure for setting various parameters @@ -54,6 +57,16 @@ kerneval = kern; end +% Assign appropriate object to chnkr +if class(chnkobj) == "chunker" + chnkr = chnkobj; +elseif class(chnkobj) == "chunkgraph" + chnkr = merge(chnkobj.echnks); +else + msg = "Unsupported object in chunkerkerneval"; + error(msg) +end + srcinfo = []; targinfo = []; srcinfo.r = chnkr.r(:,1); srcinfo.d = chnkr.d(:,1); srcinfo.n = chnkr.n(:,1); srcinfo.d2 = chnkr.d2(:,1); @@ -86,55 +99,77 @@ if isfield(opts,'fac'); opts_use.fac = opts.fac; end if isfield(opts,'eps'); opts_use.eps = opts.eps; end -[dim,~] = size(targs); +% Assign appropriate object to targinfo +targinfo = []; +if isa(targobj, "chunker") + targinfo.r = targobj.r(:,:); + targinfo.d = targobj.d(:,:); + targinfo.d2 = targobj.d2(:,:); + targinfo.n = targobj.n(:,:); +elseif isa(targobj, "chunkgraph") + targinfo.r = targobj.r(:,:); + targinfo.d = targobj.d(:,:); + targinfo.d2 = targobj.d2(:,:); + targinfo.n = targobj.n(:,:); +else + targinfo.r = targobj; +end + + + +[dim,~] = size(targinfo.r); if (dim ~= 2); warning('only dimension two tested'); end if opts_use.forcesmooth - fints = chunkerkerneval_smooth(chnkr,kern,opdims,dens,targs, ... + fints = chunkerkerneval_smooth(chnkr,kern,opdims,dens,targinfo, ... [],opts_use); return end if opts_use.forceadap fints = chunkerkerneval_adap(chnkr,kern,opdims,dens, ... - targs,[],opts_use); + targinfo,[],opts_use); return end if opts_use.forcepquad optsflag = []; optsflag.fac = opts_use.fac; - flag = flagnear(chnkr,targs,optsflag); - fints = chunkerkerneval_smooth(chnkr,kern,opdims,dens,targs, ... + flag = flagnear(chnkr,targinfo.r,optsflag); + fints = chunkerkerneval_smooth(chnkr,kern,opdims,dens,targinfo, ... flag,opts_use); fints = fints + chunkerkerneval_ho(chnkr,kern,opdims,dens, ... - targs,flag,opts_use); + targinfo,flag,opts_use); return end % smooth for sufficiently far, adaptive otherwise -%optsflag = []; optsflag.fac = opts_use.fac; -optsflag = []; -flag = flagnear_rectangle(chnkr,targs,optsflag); +rho = 1.8; +optsflag = []; optsflag.rho = rho; +flag = flagnear_rectangle(chnkr,targinfo.r,optsflag); + +npoly = chnkr.k*2; +nlegnew = chnk.ellipse_oversample(rho,npoly,opts_use.eps); +nlegnew = max(nlegnew,chnkr.k); -fints = chunkerkerneval_smooth(chnkr,kern,opdims,dens,targs, ... +[chnkr2,dens2] = upsample(chnkr,nlegnew,dens); +fints = chunkerkerneval_smooth(chnkr2,kern,opdims,dens2,targinfo, ... flag,opts_use); fints = fints + chunkerkerneval_adap(chnkr,kern,opdims,dens, ... - targs,flag,opts_use); - + targinfo,flag,opts_use); end function fints = chunkerkerneval_ho(chnkr,kern,opdims,dens, ... - targs,flag,opts) + targinfo,flag,opts) % target -[~,nt] = size(targs); +[~,nt] = size(targinfo.r); fints = zeros(opdims(1)*nt,1); k = size(chnkr.r,2); [t,w] = lege.exps(2*k); @@ -144,13 +179,26 @@ 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 intp_ab = lege.matrin(k,[-1;1]); % interpolation from k to end points + +targs = targinfo.r; + targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt); targn = zeros(chnkr.dim,nt); +if isfield(targinfo, 'd') + targd = targinfo.d; +end + +if isfield(targinfo, 'd2') + targd2 = targinfo.d2; +end + +if isfield(targinfo, 'n') + targn = targinfo.n; +end for j=1:size(chnkr.r,3) [ji] = find(flag(:,j)); if(~isempty(ji)) @@ -158,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); @@ -168,7 +216,7 @@ function fints = chunkerkerneval_smooth(chnkr,kern,opdims,dens, ... - targs,flag,opts) + targinfo,flag,opts) if isa(kern,'kernel') kerneval = kern.eval; @@ -197,12 +245,10 @@ dens = reshape(dens,opdims(2),k,nch); [~,w] = lege.exps(k); -[~,nt] = size(targs); +[~,nt] = size(targinfo.r); fints = zeros(opdims(1)*nt,1); -targinfo = []; targinfo.r = targs; - % assume smooth weights are good enough % Sequence of checks, first see ifflam is set as it supercedes @@ -225,8 +271,16 @@ end end - - +if strcmpi(imethod,'fmm') + icheck = exist(['fmm2d.' mexext], 'file'); + if icheck ~=3 + imethod = 'direct'; + fstr = ['CHUNKERKERNEVAL: forcefmm flag used but fmm2d not found\n' ... + 'Using direct computation instead\nMake sure fmm2d.' mexext ... + ' is in MATLAB path']; + warning(sprintf(fstr)) + end +end if strcmpi(imethod,'direct') % do dense version @@ -235,7 +289,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); @@ -249,7 +303,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); @@ -266,41 +320,36 @@ end else -% % hack to ignore interactions: create sparse array with very small -% % number instead of zero in ignored entries. kernbyindexr overwrites -% % with that small number -% [i,j] = find(flag); -% -% % targets correspond to multiple outputs (if matrix valued kernel) -% inew = (opdims(1)*(i(:)-1)).' + (1:opdims(1)).'; inew = inew(:); -% jnew = (j(:)).' + 0*(1:opdims(1)).'; jnew = jnew(:); -% -% % source chunks have multiple points and can be multi-dimensional -% jnew = (opdims(2)*chnkr.k*(jnew(:)-1)).' + (1:(chnkr.k*opdims(2))).'; -% jnew = jnew(:); -% inew = (inew(:)).' + 0*(1:(chnkr.k*opdims(2))).'; -% inew = inew(:); -% -% mm = nt*opdims(1); nn = chnkr.npt*opdims(2); -% v = 1e-300*ones(length(inew),1); sp = sparse(inew,jnew,v,mm,nn); - - wts = weights(chnkr); + wts = chnkr.wts; wts = wts(:); if strcmpi(imethod,'flam') xflam1 = chnkr.r(:,:); xflam1 = repmat(xflam1,opdims(2),1); xflam1 = reshape(xflam1,chnkr.dim,numel(xflam1)/chnkr.dim); - targsflam = repmat(targs(:,:),opdims(1),1); - targsflam = reshape(targsflam,chnkr.dim,numel(targsflam)/chnkr.dim); - % matfun = @(i,j) chnk.flam.kernbyindexr(i,j,targs,chnkr,wts,kern, ... - % opdims,sp); - matfun = @(i,j) chnk.flam.kernbyindexr(i,j,targs,chnkr,wts,kerneval, ... - opdims); + + targinfo_flam = []; + targinfo_flam.r = repelem(targinfo.r(:,:),1,opdims(1)); + if isfield(targinfo, 'd') + targinfo_flam.d = repelem(targinfo.d(:,:),1,opdims(1)); + end + + if isfield(targinfo, 'd2') + targinfo_flam.d2 = repelem(targinfo.d2(:,:),1,opdims(1)); + end + + if isfield(targinfo, 'n') + targinfo_flam.n = repelem(targinfo.n(:,:),1,opdims(1)); + end + +% TODO: Pull through data? + + matfun = @(i,j) chnk.flam.kernbyindexr(i, j, targinfo_flam, ..., + chnkr, kerneval, opdims); width = max(abs(max(chnkr)-min(chnkr)))/3; - tmax = max(targs(:,:),[],2); tmin = min(targs(:,:),[],2); + tmax = max(targinfo.r(:,:),[],2); tmin = min(targinfo.r(:,:),[],2); wmax = max(abs(tmax-tmin)); width = max(width,wmax/3); npxy = chnk.flam.nproxy_square(kerneval,width); @@ -310,21 +359,21 @@ ctr,chnkr,wts,kerneval,opdims,pr,ptau,pw,pin); optsifmm=[]; optsifmm.Tmax=Inf; - F = ifmm(matfun,targsflam,xflam1,200,1e-14,pxyfun,optsifmm); + F = ifmm(matfun,targinfo_flam.r,xflam1,200,1e-14,pxyfun,optsifmm); fints = ifmm_mv(F,dens(:),matfun); else wts2 = repmat(wts(:).', opdims(2), 1); sigma = wts2(:).*dens(:); - fints = kern.fmm(1e-14, chnkr, targs(:,:), sigma); + fints = kern.fmm(1e-14, chnkr, targinfo.r(:,:), sigma); end % delete interactions in flag array (possibly unstable approach) - targinfo = []; + if ~isempty(flag) 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); @@ -334,9 +383,24 @@ delsmooth = find(flag(:,i)); delsmoothrow = (opdims(1)*(delsmooth(:)-1)).' + (1:opdims(1)).'; delsmoothrow = delsmoothrow(:); - targinfo.r = targs(:,delsmooth); - kernmat = kerneval(srcinfo,targinfo); + targinfo_use = []; + targinfo_use.r = targinfo.r(:,delsmooth); + + if isfield(targinfo, 'd') + targinfo_use.d = targinfo.d(:,delsmooth); + end + + if isfield(targinfo, 'd2') + targinfo_use.d2 = targinfo.d2(:,delsmooth); + end + + if isfield(targinfo, 'n') + targinfo_use.n = targinfo.n(:,delsmooth); + end + + + kernmat = kerneval(srcinfo,targinfo_use); fints(delsmoothrow) = fints(delsmoothrow) - kernmat*densvals; end end @@ -345,7 +409,7 @@ end function fints = chunkerkerneval_adap(chnkr,kern,opdims,dens, ... - targs,flag,opts) + targinfo,flag,opts) if isa(kern,'kernel') kerneval = kern.eval; @@ -366,7 +430,23 @@ assert(numel(dens) == opdims(2)*k*nch,'dens not of appropriate size') dens = reshape(dens,opdims(2),k,nch); +% Extract target info +targs = targinfo.r; [~,nt] = size(targs); +targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt); +targn = zeros(chnkr.dim,nt); +if isfield(targinfo, 'd') + targd = targinfo.d; +end + +if isfield(targinfo, 'd2') + targd2 = targinfo.d2; +end + +if isfield(targinfo, 'n') + targn = targinfo.n; +end + fints = zeros(opdims(1)*nt,1); @@ -379,14 +459,12 @@ d = chnkr.d; n = chnkr.n; d2 = chnkr.d2; -h = chnkr.h; -targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt); -targn = zeros(chnkr.dim,nt); + 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); @@ -397,7 +475,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 c0a96c9..4169d53 100644 --- a/chunkie/chunkerkernevalmat.m +++ b/chunkie/chunkerkernevalmat.m @@ -1,4 +1,4 @@ -function mat = chunkerkernevalmat(chnkr,kern,targs,opts) +function mat = chunkerkernevalmat(chnkr,kern,targobj,opts) %CHUNKERKERNEVALMAT compute the matrix which maps density values on % the chunk geometry to the value of the convolution of the given % integral kernel with the density at the specified target points @@ -6,9 +6,21 @@ % Syntax: mat = chunkerkernevalmat(chnkr,kern,targs,opts) % % Input: -% chnkr - chunker object description of curve -% kern - integral kernel taking inputs kern(srcinfo,targinfo) -% targs - targ(1:2,i) gives the coords of the ith target +% chnkr - chunker object describing boundary, currently +% only supports chunkers, and not chunkgraphs +% kern - kernel function. By default, this should be a function handle +% accepting input of the form kern(srcinfo,targinfo), where srcinfo +% and targinfo are in the ptinfo struct format, i.e. +% ptinfo.r - positions (2,:) array +% ptinfo.d - first derivative in underlying +% parameterization (2,:) +% ptinfo.n - unit normals (2,:) +% ptinfo.d2 - second derivative in underlying +% parameterization (2,:) +% targobj - object describing the target points, can be specified as +% * array of points +% * chunker object +% * chunkgraph object % % Optional input: % opts - structure for setting various parameters @@ -38,6 +50,37 @@ % author: Travis Askham (askhamwhat@gmail.com) +% convert kernel to kernel object, put in singularity info +% opts.sing provides a default value for singularities if not +% defined for kernels + +if isa(kern,'function_handle') + kern2 = kernel(kern); + kern = kern2; +elseif isa(kern,'cell') + sz = size(kern); + kern2(sz(1),sz(2)) = kernel(); + for j = 1:sz(2) + for i = 1:sz(1) + if isa(kern{i,j},'function_handle') + kern2(i,j) = kernel(kern{i,j}); + elseif isa(kern{i,j},'kernel') + kern2(i,j) = kern{i,j}; + else + msg = "Second input is not a kernel object, function handle, " ... + + "or cell array"; + error(msg); + end + end + end + kern = kern2; + +elseif ~isa(kern,'kernel') + msg = "Second input is not a kernel object, function handle, " ... + + "or cell array"; + error(msg); +end + % determine operator dimensions using first two points @@ -48,7 +91,9 @@ targinfo.r = chnkr.r(:,2); targinfo.d = chnkr.d(:,2); targinfo.d2 = chnkr.d2(:,2); targinfo.n = chnkr.n(:,2); -ftemp = kern(srcinfo,targinfo); +ftmp = kern.eval; + +ftemp = ftmp(srcinfo,targinfo); opdims = size(ftemp); if nargin < 4 @@ -57,18 +102,34 @@ forcesmooth = false; forceadap = false; -forcepqud = false; +forcepquad = false; nonsmoothonly = false; fac = 1.0; eps = 1e-12; if isfield(opts,'forcesmooth'); forcesmooth = opts.forcesmooth; end if isfield(opts,'forceadap'); forceadap = opts.forceadap; end -if isfield(opts,'forcepquad'); forcepqud = opts.forcepquad; end +if isfield(opts,'forcepquad'); forcepquad = opts.forcepquad; end if isfield(opts,'nonsmoothonly'); nonsmoothonly = opts.nonsmoothonly; end if isfield(opts,'fac'); fac = opts.fac; end if isfield(opts,'eps'); eps = opts.eps; end -[dim,~] = size(targs); +% Assign appropriate object to targinfo +targinfo = []; +if isa(targobj, "chunker") + targinfo.r = targobj.r(:,:); + targinfo.d = targobj.d(:,:); + targinfo.d2 = targobj.d2(:,:); + targinfo.n = targobj.n(:,:); +elseif isa(targobj, "chunkgraph") + targinfo.r = targobj.r(:,:); + targinfo.d = targobj.d(:,:); + targinfo.d2 = targobj.d2(:,:); + targinfo.n = targobj.n(:,:); +else + targinfo.r = targobj; +end + +[dim,~] = size(targinfo.r); if (dim ~= 2); warning('only dimension two tested'); end @@ -77,24 +138,26 @@ optsadap = []; optsadap.eps = eps; + + if forcesmooth - mat = chunkerkernevalmat_smooth(chnkr,kern,opdims,targs, ... + mat = chunkerkernevalmat_smooth(chnkr,ftmp,opdims,targinfo, ... [],optssmooth); return end if forceadap - mat = chunkerkernevalmat_adap(chnkr,kern,opdims, ... - targs,[],optsadap); + mat = chunkerkernevalmat_adap(chnkr,ftmp,opdims, ... + targinfo,[],optsadap); return end -if forcepqud +if forcepquad optsflag = []; optsflag.fac = fac; - flag = flagnear(chnkr,targs,optsflag); - spmat = chunkerkernevalmat_ho(chnkr,kern,opdims, ... - targs,flag,optsadap); - mat = chunkerkernevalmat_smooth(chnkr,kern,opdims,targs, ... + flag = flagnear(chnkr,targinfo.r,optsflag); + spmat = chunkerkernevalmat_ho(chnkr,ftmp,opdims, ... + targinfo,flag,optsadap); + mat = chunkerkernevalmat_smooth(chnkr,ftmp,opdims,targinfo, ... flag,opts); mat = mat + spmat; return @@ -102,17 +165,20 @@ % smooth for sufficiently far, adaptive otherwise +% TODO: change to chunkerkerneval system, need routine to generate +% upsampling matrix. + optsflag = []; optsflag.fac = fac; -flag = flagnear(chnkr,targs,optsflag); -spmat = chunkerkernevalmat_adap(chnkr,kern,opdims, ... - targs,flag,optsadap); +flag = flagnear(chnkr,targinfo.r,optsflag); +spmat = chunkerkernevalmat_adap(chnkr,ftmp,opdims, ... + targinfo,flag,optsadap); if nonsmoothonly mat = spmat; return; end -mat = chunkerkernevalmat_smooth(chnkr,kern,opdims,targs, ... +mat = chunkerkernevalmat_smooth(chnkr,ftmp,opdims,targinfo, ... flag,opts); mat = mat + spmat; @@ -123,7 +189,7 @@ function mat = chunkerkernevalmat_smooth(chnkr,kern,opdims, ... - targs,flag,opts) + targinfo,flag,opts) if nargin < 6 flag = []; @@ -135,7 +201,6 @@ k = chnkr.k; nch = chnkr.nch; -targinfo = []; targinfo.r = targs; srcinfo = []; srcinfo.r = chnkr.r(:,:); srcinfo.n = chnkr.n(:,:); srcinfo.d = chnkr.d(:,:); srcinfo.d2 = chnkr.d2(:,:); @@ -163,7 +228,7 @@ end function mat = chunkerkernevalmat_adap(chnkr,kern,opdims, ... - targs,flag,opts) + targinfo,flag,opts) k = chnkr.k; nch = chnkr.nch; @@ -175,12 +240,30 @@ opts = []; end +% Extract target info +targs = targinfo.r; [~,nt] = size(targs); +targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt); +targn = zeros(chnkr.dim,nt); +if isfield(targinfo, 'd') + targd = targinfo.d; +end + +if isfield(targinfo, 'd2') + targd2 = targinfo.d2; +end + +if isfield(targinfo, 'n') + targn = targinfo.n; +end + % using adaptive quadrature if isempty(flag) + + mat = zeros(opdims(1)*nt,opdims(2)*chnkr.npt); [t,w] = lege.exps(2*k+1); @@ -191,7 +274,7 @@ n = chnkr.n; d2 = chnkr.d2; h = chnkr.h; - targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt); + for i = 1:nch jmat = 1 + (i-1)*k*opdims(2); jmatend = i*k*opdims(2); @@ -227,15 +310,12 @@ d = chnkr.d; n = chnkr.n; d2 = chnkr.d2; - h = chnkr.h; - targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt); - targn = zeros(chnkr.dim,nt); 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; @@ -260,7 +340,7 @@ end function mat = chunkerkernevalmat_ho(chnkr,kern,opdims, ... - targs,flag,opts) + targinfo,flag,opts) k = chnkr.k; nch = chnkr.nch; @@ -272,7 +352,24 @@ opts = []; end +% Extract target info +targs = targinfo.r; [~,nt] = size(targs); +targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt); +targn = zeros(chnkr.dim,nt); +if isfield(targinfo, 'd') + targd = targinfo.d; +end + +if isfield(targinfo, 'd2') + targd2 = targinfo.d2; +end + +if isfield(targinfo, 'n') + targn = targinfo.n; +end + + % using Helsing-Ojala quadrature if isempty(flag) % figure out what is this flag for in adaptive routine @@ -290,13 +387,11 @@ 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 intp_ab = lege.matrin(k,[-1;1]); % interpolation from k to end points - targd = zeros(chnkr.dim,nt); targd2 = zeros(chnkr.dim,nt); - targn = zeros(chnkr.dim,nt); + for i = 1:nch jmat = 1 + (i-1)*k*opdims(2); jmatend = i*k*opdims(2); @@ -304,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/chunkermat.m b/chunkie/chunkermat.m index 5e1897f..f3fd7bd 100644 --- a/chunkie/chunkermat.m +++ b/chunkie/chunkermat.m @@ -7,7 +7,7 @@ % Syntax: sysmat = chunkermat(chnkr,kern,opts) % % Input: -% chnkr - chunker object describing boundary +% chnkobj - chunker object describing boundary % kern - kernel function. By default, this should be a function handle % accepting input of the form kern(srcinfo,targinfo), where srcinfo % and targinfo are in the ptinfo struct format, i.e. @@ -438,6 +438,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; @@ -467,7 +468,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); @@ -482,13 +483,14 @@ end - [Pbc,PWbc,starL,circL,starS,circS,ilist] = chnk.rcip.setup(ngl,ndim, ... - nedge,isstart); - - % this might need to be fixed in triple junction case + [Pbc,PWbc,starL,circL,starS,circS,ilist,starL1,circL1] = ... + chnk.rcip.setup(ngl,ndim,nedge,isstart); + optsrcip = opts; + optsrcip.nonsmoothonly = false; + 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,optsrcip); sysmat_tmp = inv(R) - eye(2*ngl*nedge*ndim); if (~nonsmoothonly) 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/chunkie/chunkgraphinit.m b/chunkie/chunkgraphinit.m index 01ffde6..044b426 100644 --- a/chunkie/chunkgraphinit.m +++ b/chunkie/chunkgraphinit.m @@ -1,100 +1,7 @@ function [cgrph] = chunkgraphinit(verts,edge2verts,fchnks,cparams) warning('This method is deprecated. Use the chunkgraph constructor instead.'); - prefs = []; - cgrph = chunkgraph(prefs); - cgrph.verts = verts; - cgrph.edge2verts = edge2verts; - cgrph.echnks = chunker.empty; - - if (nargin < 4) - cploc = []; - cploc.ta = 0; - cploc.tb = 1; - cploc.ifclosed = 0; - cploc.nover = 1; - cploc.eps = 1.0d-10; - cploc.lvlr = 'a'; - else - cploc = cparams; - %mandatory settings - cploc.ta = 0; - cploc.tb = 1; - cploc.ifclosed = 0; - if (~isfield(cparams,'lvlr')) - cploc.lvlr = 'a'; - end - if (~isfield(cparams,'eps')) - cploc.eps = 1.0d-10; - end - if (~isfield(cparams,'nover')) - cploc.nover = 1; - end - end - - - pref = []; - pref.nchmax = 10000; - pref.k = 16; - - if (size(verts,2) ~= size(edge2verts,2)) - error('Incompatible vertex and edge sizes'); - end - - echnks = chunker.empty(); - for i=1:size(edge2verts,1) - if (numel(fchnks) 0 && mod(nin,2)==1) + if (irgn ~= 0) + disp("Warning: an unsupported geometry error has occurred"); + end + irgn = ii; + end + end + + if (irgn ~= 0) + rgnout = rgn1; + rgnout = [rgnout,rgn2(2:end)]; + rgnout{irgn} = [rgnout{irgn},rgn2{1}]; + end + + irgn1 = irgn; + + if (irgn == 0) + e1 = rgn1{1}{1}(1); + v1 = cgrph.edgesendverts(2,abs(e1)); + v1 = cgrph.verts(:,v1); + + irgn = 0; + for ii=2:numel(rgn2) + nin = pointinregion(cgrph,rgn2{ii},v1); + if (nin > 0 && mod(nin,2)==1) + if (irgn ~= 0) + disp("Warning: an unsupported geometry error has occurred"); + end + irgn = ii; + end + end + + if (irgn ~= 0) + rgnout = rgn2; + rgnout = [rgnout,rgn1(2:end)]; + rgnout{irgn} = [rgnout{irgn},rgn1{1}]; + end + end + + if (irgn1==0 && irgn==0) + rgnout = rgn1; + rgnout = [rgnout,rgn2(2:end)]; + rgnout{1} = [rgnout{1},rgn2{1}]; + end + +end + diff --git a/chunkie/pointinregion.m b/chunkie/pointinregion.m new file mode 100644 index 0000000..61afddb --- /dev/null +++ b/chunkie/pointinregion.m @@ -0,0 +1,19 @@ +function [nin] = pointinregion(cgrph,rgn,r0) + nin = 0; + for ii=1:numel(rgn) + iedges = rgn{ii}; + rs = []; + for jj=1:numel(iedges) + rtmp = sort(cgrph.echnks(abs(iedges(jj)))).r(:,:); + if (iedges(jj)<0) + rtmp = fliplr(rtmp(:,:)); + end + rs = [rs,rtmp]; + end + [in] = inpolygon(r0(1),r0(2),rs(1,:),rs(2,:)); + if (in == 1) + nin = nin + 1; + end + end +end + diff --git a/chunkie/regioninside.m b/chunkie/regioninside.m new file mode 100644 index 0000000..c49732a --- /dev/null +++ b/chunkie/regioninside.m @@ -0,0 +1,18 @@ +function [isinside] = regioninside(cgrph,rgn1,rgn2) + + e2 = rgn2{1}{1}(1); + v2 = cgrph.edgesendverts(2,abs(e2)); + v2 = cgrph.verts(:,v2); + irgn = 0; + for ii=2:numel(rgn1) + nin = pointinregion(cgrph,rgn1{ii},v2); + if (nin > 0 && mod(nin,2)==1) + if (irgn ~= 0) + disp("Warning: an unsupported geometry error has occurred"); + end + irgn = ii; + end + end + isinside = (irgn ~=0); +end + diff --git a/chunkie/roundedpolygon.m b/chunkie/roundedpolygon.m new file mode 100644 index 0000000..867e3ca --- /dev/null +++ b/chunkie/roundedpolygon.m @@ -0,0 +1,262 @@ +function chnkr = roundedpolygon(verts,cparams,pref,edgevals) +%CHUNKPOLY return a chunker object corresponding to the polygon with the +% given vertices. By default, a polygon with rounded corners is returned. +% Open and higher dimensional "polygons" are allowed. Optional dyadic +% refinement near corners of true polygons (no rounding). +% +% Syntax: chnkr = chunkpoly(verts,cparams,pref,edgevals) +% +% Input: +% verts - (dimv,nverts) array of "polygon" vertices +% in order +% +% Optional input: +% cparams - options structure +% cparams.rounded = true if corner rounding is to be used. +% false if no rounding is used (true) +% cparams.widths = radius around each corner point where either the +% rounding or dyadic refinement is applied +% (defaults to 1/10th of minimum +% of length of adjoining edges) +% cparams.autowidths = automatically compute widths (false) +% cparams.autowidthsfac = if using autowidths, set widths +% to autowidthsfac*minimum of adjoining +% edges (0.1) +% cparams.ifclosed = true, if it's a closed polygon +% false, if it's an open segment (true) +% cparams may also specify any of the values to send to the +% refine routine, otherwise it is called with defaults +% +% ~ Rounding parameters ~ +% cparams.eps - resolve curve to tolerance eps +% resolve coordinates, arclength, +% and first and second derivs of coordinates +% to this tolerance (1.0e-6) +% +% pref - chunkerpref object/ preference structure +% pref.k determines order of underlying Legendre nodes +% edgevals - very optional input. specifies constant values along each +% edge. The routine then stores a smoothed interpolant +% of these edge values on the rounded structure in the +% output chunker's data field +% +% Output: +% chnkr - chunker object corresponding to (rounded) polygon +% +% Examples: +% barbell_verts = chnk.demo.barbell(); +% chnkr = roundedpolygon(barbell_verts); % rounded "barbell" domain +% % with standard options +% cparams = []; cparams.rounded = false; +% pref = []; pref.k = 30; +% chnkr = chunkpoly(barbell_verts,cparams,pref); % not rounded +% +% See also CHUNKERFUNC, CHUNKER, CHUNKERPREF, REFINE + +[dimv,nv] = size(verts); + +assert(dimv > 1) + +if nargin < 2 + cparams = []; +end +if nargin < 3 + p = []; + p.dim = dimv; + pref = chunkerpref(p); +else + pref = chunkerpref(pref); + if pref.dim ~= dimv + warning('dimensions dont match, overwriting with vertex dim'); + pref.dim = dimv; + end +end +if nargin < 4 + edgevals = []; +end + +autowidths = false; +autowidthsfac = 0.1; +ifclosed = true; +eps = 1e-6; +nover = 0; + +if isfield(cparams,'ifclosed') + ifclosed = cparams.ifclosed; +end +if isfield(cparams,'eps') + eps = cparams.eps; +end +if isfield(cparams,'nover') + nover = cparams.nover; +end + +if (ifclosed) + verts2 = [verts(:,end), verts, verts(:,1)]; + edges2 = sqrt(sum(diff(verts2,1,2).^2,1)); +else + edges2 = [0, sqrt(sum(diff(verts,1,2).^2,1)), 0]; +end + +widths_not_set = true; + +if isfield(cparams,'widths') + widths = cparams.widths; + widths_not_set = false; +end +if isfield(cparams,'autowidths') + autowidths = cparams.autowidths; +end +if isfield(cparams,'autowidthsfac') + autowidthsfac = cparams.autowidthsfac; +end + +if (autowidths || widths_not_set) + widths = autowidthsfac*... + min(edges2(1:end-1),edges2(2:end)); +end + +if ifclosed + widths = [widths(:); widths(1)]; + verts = [verts, verts(:,1)]; + nv = size(verts,2); +end + +nedge = nv - 1; +nvals = 0; +if numel(edgevals) ~= 0 + assert(rem(numel(edgevals),nedge) == 0, ... + 'number of edge values should be multiple of number of edges'); + nvals = numel(edgevals)/nedge; + edgevals = reshape(edgevals,nvals,nedge); +end + +chnkr = chunker(pref); +chnkr = chnkr.makedatarows(nvals); +k = chnkr.k; dim = chnkr.dim; +[t,w] = lege.exps(k); +onek = ones(k,1); + + +if nvals > 0 + aint = lege.intmat(chnkr.k); +end + +for i = 1:nv-1 + % grab vertices + r1 = verts(:,i); r2 = verts(:,i+1); + w1 = widths(i); w2 = widths(i+1); + l = sqrt(sum((r1-r2).^2)); + assert(l > w1+w2+2*eps(1)*l,'widths too large for side'); + + % make chunk in middle + v = (r2-r1)/l; + ts = w1+(l-w2-w1)*(t+1)/2.0; + chnkr = chnkr.addchunk(); + nch = chnkr.nch; + if nvals > 0 + val1 = edgevals(:,i); + chnkr.data(:,:,nch) = bsxfun(@times,val1,onek(:).'); + end + chnkr.r(:,:,nch) = r1 + bsxfun(@times,v(:),(ts(:)).'); + chnkr.d(:,:,nch) = repmat(v(:),1,k); + chnkr.d2(:,:,nch) = zeros(dim,k); + 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) + % chunk up smoothed corner made by three verts + if (i==nv-1) + if nvals > 0 + val2 = edgevals(:,1); + end + r3 = verts(:,2); + else + if nvals > 0 + val2 = edgevals(:,i+1); + end + r3 = verts(:,i+2); + end + l2 = sqrt(sum((r2-r3).^2)); + v = -v; + v2 = (r3-r2)/l2; + cosphi = dot(v2,v); %cosine of angle between edges + sinphi = sqrt(1-cosphi*cosphi); + trange = w2*sqrt((1-cosphi)/2); %parameter space range of gaussian + hbell = abs(trange)/8.0; % width parameter of gaussian + m = sinphi/(1-cosphi); % slope of abs approx + cpt.ta = -trange; + cpt.tb = trange; + cpt.eps = eps; cpt.levrestr = 0; cpt.ifclosed = 0; + chnkrt = sort(chunkerfunc(@(t)fround(t,m,hbell,dim),cpt,pref)); + + % do optimal procrustes match of left and right ends + rl = fround(-trange,m,hbell,dim); rr = fround(trange,m,hbell,dim); + [um,~,vm] = svd([w2*v w2*v2]*([rl rr].')); + rotmat = um*vm.'; + + % copy in rotated and translated chunks + ncht = chnkrt.nch; + chnkr = chnkr.addchunk(ncht); + chnkr.r(:,:,nch+1:nch+ncht) = reshape(rotmat*chnkrt.r(:,:),dim,k,ncht)+r2; + chnkr.d(:,:,nch+1:nch+ncht) = reshape(rotmat*chnkrt.d(:,:),dim,k,ncht); + chnkr.d2(:,:,nch+1:nch+ncht) = reshape(rotmat*chnkrt.d2(:,:),dim,k,ncht); + 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(:)).'); + dsw = bsxfun(@times,w(:),ds); + dssums = sum(dsw,1); + dssums2 = cumsum([0,dssums(1:end-1)]); + dsint = aint*ds; + dsint = bsxfun(@plus,dsint,dssums2); + lencorner = sum(dsw(:)); + ss = -lencorner/2.0 + dsint; + ss = ss/lencorner*16; + erfss = erf(ss); + datass = reshape((val2(:)-val1(:))/2*((erfss(:)).'+1) ... + +val1(:),nvals,k,ncht); + chnkr.data(:,:,nch+1:nch+ncht) = datass; + end + end + +end + +if ifclosed + nch = chnkr.nch; + chnkr.adj(1,1) = nch; + chnkr.adj(2,nch) = 1; +end + +chnkr = chnkr.refine(cparams); + +% compute normals + +chnkr.n = normals(chnkr); + +end + +function [r,d,d2] = fround(t,m,h,dim) + +[y,dy,d2y] = chnk.spcl.absconvgauss(t,m,0.0,h); + +r = zeros(dim,length(t)); +d = zeros(dim,length(t)); +d2 = zeros(dim,length(t)); + +r(1,:) = t; +r(2,:) = y; +d(1,:) = 1.0; +d(2,:) = dy; +d2(1,:) = 0.0; +d2(2,:) = d2y; + +end diff --git a/chunkie/starfish.m b/chunkie/starfish.m index 4d812e3..dc53d90 100644 --- a/chunkie/starfish.m +++ b/chunkie/starfish.m @@ -1,19 +1,33 @@ function [r,d,d2] = starfish(t,varargin) -%STARFISH -% return position, first and second derivatives of parameterized starfish -% domain. +%STARFISH return position, first and second derivatives of a starfish +% domain with the parameterization % -% Inputs: -% t - points (in [0,2pi] to evaluate these quantities +% x(t) = x0 + (1+amp*cos(narms*(t+phi)))*cos(t) +% y(t) = y0 + (1+amp*cos(narms*(t+phi)))*sin(t) % -% Optional inputs: -% narms - integer, number of arms on starfish (5) -% amp - float, amplitude of starfish arms relative to radius of length 1 +% Syntax: [r,d,d2] = starfish(t,narms,amp,ctr,phi,scale) +% +% Input: +% t - array of points (in [0,2pi]) +% +% Optional input: +% narms - integer, number of arms on starfish (5) +% amp - float, amplitude of starfish arms relative to radius of length 1 % (0.3) -% ctr - float(2), x,y coordinates of center of starfish ( [0,0] ) -% phi - float, phase shift (0) -% scale - scaling factor (1.0) +% ctr - float(2), x0,y0 coordinates of center of starfish ( [0,0] ) +% phi - float, phase shift (0) +% scale - scaling factor (1.0) +% +% Output: +% r - 2 x numel(t) array of positions, r(:,i) = [x(t(i)); y(t(i))] +% d - 2 x numel(t) array of t derivative of r +% d2 - 2 x numel(t) array of second t derivative of r +% +% Examples: +% [r,d,d2] = starfish(t); % get default settings +% [r,d,d2] = starfish(t,narms,[],ctr,[],scale); % change some settings +% narms = 5; amp = 0.3; diff --git a/chunkie/trapperkerneval.m b/chunkie/trapperkerneval.m index d425b26..aaf7939 100644 --- a/chunkie/trapperkerneval.m +++ b/chunkie/trapperkerneval.m @@ -1,10 +1,21 @@ -function fints = trapperkerneval(trap,kern,dens,targs,opts) +function fints = trapperkerneval(trap,kern,dens,targobj,opts) %TRAPPERINTKERN compute the convolution of the integral kernel with wts = trap.wts; srcinfo = []; srcinfo.r = trap.r; srcinfo.d = trap.d; srcinfo.d2 = trap.d2; srcinfo.n = trap.n; -targinfo = []; targinfo.r = targs; + +% Assign appropriate object to targinfo +targinfo = []; +if isa(targobj, "trapper") + targinfo.r = targobj.r(:,:); + targinfo.d = targobj.d(:,:); + targinfo.d2 = targobj.d2(:,:); + targinfo.n = targobj.n(:,:); +else + targinfo.r = targobj; +end + mat = kern(srcinfo,targinfo); fints = mat*diag(wts)*dens; 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/chunkerfitTest.m b/devtools/test/chunkerfitTest.m new file mode 100644 index 0000000..a0cef29 --- /dev/null +++ b/devtools/test/chunkerfitTest.m @@ -0,0 +1,24 @@ + +%CHUNKERFITTEST + +clearvars; close all; +addpaths_loc(); + +% Sample a smooth curve at random points +rng(0) +n = 20; +tt = sort(2*pi*rand(n,1)); +r = chnk.curves.bymode(tt, [2 0.5 0.2 0.7]); + +opts = []; +opts.ifclosed = true; +opts.cparams = []; +opts.cparams.eps = 1e-6; +opts.pref = []; +opts.pref.k = 16; +chnkr = chunkerfit(r, opts); +assert(checkadjinfo(chnkr) == 0); + +opts.ifclosed = false; +chnkr = chunkerfit(r(:,1:10), opts); +assert(checkadjinfo(chnkr) == 0); diff --git a/devtools/test/chunkerinteriorTest.m b/devtools/test/chunkerinteriorTest.m index 549c53a..2538be6 100644 --- a/devtools/test/chunkerinteriorTest.m +++ b/devtools/test/chunkerinteriorTest.m @@ -57,18 +57,30 @@ fprintf('%5.2e s : time for chunkerinterior (with flam)\n',t2); fprintf('%5.2e s : time for chunkerinterior (with fmm)\n',t3); -ifail = in2(:) ~= (scal(:) < 1); assert(all(in(:) == (scal(:) < 1))); assert(all(in2(:) == (scal(:) < 1))); assert(all(in3(:) == (scal(:) < 1))); -% -% -% figure(1) -% clf -% hold off -% scatter(targs(1,in),targs(2,in),'go') -% hold on -% scatter(targs(1,~in),targs(2,~in),'rx') -% plot(chnkr,'b','LineWidth',2) -% axis equal +x1 = linspace(-2,2,100); +[xx,yy] = meshgrid(x1,x1); + +opts = []; +opts.fmm = true; +opts.flam = false; +start = tic; in4 = chunkerinterior(chnkr,{x1,x1},opts); t4 = toc(start); +fprintf('%5.2e s : time for chunkerinterior (meshgrid, with fmm)\n',t4); + +% Test targets specified as chunker +narms = 3; +amp = 0.1; +chnkr2 = chunkerfunc(@(t) 0.3*starfish(t,narms,amp),cparams,pref); + +opts = []; +opts.fmm = true; +opts.flam = false; +start = tic; in5 = chunkerinterior(chnkr,chnkr2,opts); t5 = toc(start); +fprintf('%5.2e s : time for chunkerinterior (chunker, with fmm)\n',t5); +assert(all(in5(:) == 1)); + + + diff --git a/devtools/test/chunkerkerneval_greenhelmTest.m b/devtools/test/chunkerkerneval_greenhelmTest.m index d9a9e31..e4e2cef 100644 --- a/devtools/test/chunkerkerneval_greenhelmTest.m +++ b/devtools/test/chunkerkerneval_greenhelmTest.m @@ -76,11 +76,9 @@ % test green's id -opts.forcesmooth=false; -opts.quadkgparams = {'RelTol',1.0e-13,'AbsTol',1.0e-13}; -start=tic; Du = chunkerkerneval(chnkr,kernd,densu,targets,opts); +start=tic; Du = chunkerkerneval(chnkr,kernd,densu,targets); toc(start) -start=tic; Sun = chunkerkerneval(chnkr,kerns,densun,targets,opts); +start=tic; Sun = chunkerkerneval(chnkr,kerns,densun,targets); toc(start) utarg2 = Sun-Du; diff --git a/devtools/test/chunkerkerneval_greenlapTest.m b/devtools/test/chunkerkerneval_greenlapTest.m index 9cccf63..ecf5c92 100644 --- a/devtools/test/chunkerkerneval_greenlapTest.m +++ b/devtools/test/chunkerkerneval_greenlapTest.m @@ -23,7 +23,7 @@ % sources -ns = 10; +ns = 100; ts = 0.0+2*pi*rand(ns,1); sources = starfish(ts,narms,amp); sources = 3.0*sources; @@ -31,7 +31,7 @@ % targets -nt = 100; +nt = 300; ts = 0.0+2*pi*rand(nt,1); targets = starfish(ts,narms,amp); targets = targets.*repmat(rand(1,nt),2,1); @@ -52,9 +52,9 @@ % kernel defs -kernd = @(s,t) chnk.lap2d.kern(s,t,'d'); -kerns = @(s,t) chnk.lap2d.kern(s,t,'s'); -kernsprime = @(s,t) chnk.lap2d.kern(s,t,'sprime'); +kernd = kernel('lap','d'); +kerns = kernel('lap','s'); +kernsprime = kernel('lap','sprime'); opdims = [1 1]; @@ -63,21 +63,22 @@ srcinfo = []; srcinfo.r = sources; targinfo = []; targinfo.r = chnkr.r(:,:); targinfo.d = chnkr.d(:,:); -kernmats = kerns(srcinfo,targinfo); -kernmatsprime = kernsprime(srcinfo,targinfo); +kernmats = kerns.eval(srcinfo,targinfo); +kernmatsprime = kernsprime.eval(srcinfo,targinfo); densu = kernmats*strengths; densun = kernmatsprime*strengths; % eval u at targets targinfo = []; targinfo.r = targets; -kernmatstarg = kerns(srcinfo,targinfo); +kernmatstarg = kerns.eval(srcinfo,targinfo); utarg = kernmatstarg*strengths; -% test green's id +% test green's id. we use FLAM, direct and FMM for smooth work -opts.quadkgparams = {'RelTol',1.0e-13,'AbsTol',1.0e-13}; +opts = []; +opts.flam = true; start=tic; Du = chunkerkerneval(chnkr,kernd,densu,targets,opts); toc(start) start=tic; Sun = chunkerkerneval(chnkr,kerns,densun,targets,opts); @@ -89,7 +90,41 @@ relerr = norm(utarg-utarg2,'fro')/norm(utarg,'fro'); -fprintf('relative frobenius error %5.2e\n',relerr); +fprintf('relative frobenius error, forcing flam %5.2e\n',relerr); + +assert(relerr < 1e-11); + +opts = []; +opts.accel = false; +start=tic; Du = chunkerkerneval(chnkr,kernd,densu,targets,opts); +toc(start) +start=tic; Sun = chunkerkerneval(chnkr,kerns,densun,targets,opts); +toc(start) + +utarg2 = Sun-Du; + +% + +relerr = norm(utarg-utarg2,'fro')/norm(utarg,'fro'); + +fprintf('relative frobenius error, forcing direct %5.2e\n',relerr); + +assert(relerr < 1e-11); + +opts = []; +opts.forcefmm = true; +start=tic; Du = chunkerkerneval(chnkr,kernd,densu,targets,opts); +toc(start) +start=tic; Sun = chunkerkerneval(chnkr,kerns,densun,targets,opts); +toc(start) + +utarg2 = Sun-Du; + +% + +relerr = norm(utarg-utarg2,'fro')/norm(utarg,'fro'); + +fprintf('relative frobenius error, forcing fmm %5.2e\n',relerr); assert(relerr < 1e-11); diff --git a/devtools/test/chunkermat_helm2dTest.m b/devtools/test/chunkermat_helm2dTest.m index c583647..1c1aa8c 100644 --- a/devtools/test/chunkermat_helm2dTest.m +++ b/devtools/test/chunkermat_helm2dTest.m @@ -102,7 +102,6 @@ opts.usesmooth=false; opts.verb=false; -opts.quadkgparams = {'RelTol',1e-16,'AbsTol',1.0e-16}; start=tic; Dsol = chunkerkerneval(chnkr,fkern,sol2,targets,opts); t1 = toc(start); fprintf('%5.2e s : time to eval at targs (slow, adaptive routine)\n',t1) diff --git a/devtools/test/chunkermat_helm2dfmmTest.m b/devtools/test/chunkermat_helm2dfmmTest.m deleted file mode 100644 index dcea14f..0000000 --- a/devtools/test/chunkermat_helm2dfmmTest.m +++ /dev/null @@ -1,120 +0,0 @@ - -%CHUNKERMAT_HELM2DFMMTEST -% -% test the matrix builder and the fmm for post processing and do a basic solve - -clearvars; close all; -iseed = 8675309; -rng(iseed); - -addpaths_loc(); - -cparams = []; -cparams.eps = 1.0e-10; -cparams.nover = 1; -pref = []; -pref.k = 32; -narms = 3; -amp = 0.25; -start = tic; chnkr = chunkerfunc(@(t) starfish(t,narms,amp),cparams,pref); -t1 = toc(start); - -fprintf('%5.2e s : time to build geo\n',t1) - -zk = 1.1; - -% sources - -ns = 10; -ts = 0.0+2*pi*rand(ns,1); -sources = starfish(ts,narms,amp); -sources = 3.0*sources; -strengths = randn(ns,1); - -% targets - -nt = 10000; -ts = 0.0+2*pi*rand(nt,1); -targets = starfish(ts,narms,amp); -targets = targets.*repmat(rand(1,nt),2,1)*0.5; - -% plot geo and sources - -xs = chnkr.r(1,:,:); xmin = min(xs(:)); xmax = max(xs(:)); -ys = chnkr.r(2,:,:); ymin = min(ys(:)); ymax = max(ys(:)); - -figure(1) -clf -hold off -plot(chnkr) -hold on -scatter(sources(1,:),sources(2,:),'o') -scatter(targets(1,:),targets(2,:),'x') -axis equal - -% - -kerns = @(s,t) chnk.helm2d.kern(zk,s,t,'s'); - -% eval u on bdry - -targs = chnkr.r; targs = reshape(targs,2,chnkr.k*chnkr.nch); -targstau = tangents(chnkr); -targstau = reshape(targstau,2,chnkr.k*chnkr.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; - -% - -% build helmholtz dirichlet matrix - -fkern = @(s,t) chnk.helm2d.kern(zk,s,t,'D'); -start = tic; D = chunkermat(chnkr,fkern); -t1 = toc(start); - -fprintf('%5.2e s : time to assemble matrix\n',t1) - -sys = -0.5*eye(chnkr.k*chnkr.nch) + D; - -rhs = ubdry; rhs = rhs(:); -start = tic; sol = gmres(sys,rhs,[],1e-14,100); t1 = toc(start); - -fprintf('%5.2e s : time for dense gmres\n',t1) - -start = tic; sol2 = sys\rhs; t1 = toc(start); - -fprintf('%5.2e s : time for dense backslash solve\n',t1) - -err = norm(sol-sol2,'fro')/norm(sol2,'fro'); - -fprintf('difference between direct and iterative %5.2e\n',err) - -% evaluate at targets and compare -wchnkr = chnkr.wts; - -% evaluate at targets using FMM and compare -iffmm = 0; - -if(iffmm) - sol_use = sol2.*wchnkr(:); - eps = 1e-6; - pgt = 1; - pot = chnk.helm2d.fmm(eps,zk,chnkr,targets,'D',sol_use); - - relerr = norm(utarg-pot,'fro')/(sqrt(chnkr.nch)*norm(utarg,'fro')); - relerr2 = norm(utarg-pot,'inf')/dot(abs(sol(:)),wchnkr(:)); - fprintf('relative frobenius error %5.2e\n',relerr); - fprintf('relative l_inf/l_1 error %5.2e\n',relerr2); - - assert(relerr < eps); -end - 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/chunkgrphconstructTest.m b/devtools/test/chunkgrphconstructTest.m index d829e77..3417276 100644 --- a/devtools/test/chunkgrphconstructTest.m +++ b/devtools/test/chunkgrphconstructTest.m @@ -45,60 +45,17 @@ verts = exp(1i*2*pi*(0:4)/5); verts = [real(verts);imag(verts)]; +% legacy connectivity info edge2verts = [-1, 1, 0, 0, 0; ... 0,-1, 1, 0, 0; ... 0, 0,-1, 1, 0; ... 0, 0, 0,-1, 1; ... 1, 0, 0, 0,-1]; edge2verts = sparse(edge2verts); -amp = 0.5; -frq = 6; -fchnks = {}; -for icurve = 1:size(edge2verts,1) - fchnks{icurve} = @(t) circulararc(t,cpars{1}); - fchnks{icurve} = @(t) sinearc(t,amp,frq); -end -%fchnks = {}; - -prefs = []; -% prefs.chsmall = 1d-4; -[cgrph] = chunkgraph(verts,edge2verts,fchnks,prefs); - -vstruc = procverts(cgrph); -rgns = findregions(cgrph); -cgrph = balance(cgrph); - - - -ncurve = 1; - -% 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; - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% connectivity info +edgesendverts = [1:5; [2:5 1]]; -verts = exp(1i*2*pi*(0:4)/5); -verts = [real(verts);imag(verts)]; - -edge2verts = [-1, 1, 0, 0, 0; ... - 0,-1, 1, 0, 0; ... - 0, 0,-1, 1, 0; ... - 0, 0, 0,-1, 1; ... - 1, 0, 0, 0,-1]; -edge2verts = sparse(edge2verts); amp = 0.5; frq = 6; fchnks = {}; @@ -109,86 +66,15 @@ %fchnks = {}; prefs = []; -% prefs.chsmall = 1d-4; -[cgrph] = chunkgraphinit(verts,edge2verts,fchnks,prefs); - -vstruc = procverts(cgrph); -rgns = findregions(cgrph); -cgrph = balance(cgrph); - - - -zk = 1.0; -fkern = @(s,t) -2*chnk.helm2d.kern(zk,s,t,'d'); - -opts = []; -opts.nonsmoothonly = false; -opts.rcip = true; -[sysmat] = chunkermat(cgrph,fkern,opts); -sysmat = sysmat + eye(size(sysmat,2)); - -% generate some targets... - -xs = -2:0.01:2; -ys = -2:0.01:2; -[X,Y] = meshgrid(xs,ys); -targs = [X(:).';Y(:).']; - -srcinfo = []; -srcinfo.sources = cgrph.r(:,:); -w = weights(cgrph); -n = normals(cgrph); - -% a quick hack to find the interior points - -srcinfo.dipstr = w(:).'; -srcinfo.dipvec = n(:,:); -eps = 1E-8; -pg = 0; -pgt = 1; -[U] = lfmm2d(eps,srcinfo,pg,targs,pgt); -U = U.pottarg; -inds = find(abs(U-2*pi) 0 && mod(nin,2)==1) +% if (irgn ~= 0) +% display("Warning: an unsupported geometry error has occurred"); +% end +% irgn = ii; +% end +% end +% +% if (irgn == 0) +% +% end +% +% rgnout = rgns{1}; +% if (numel(rgns)>1) +% rgn2 = rgns{2}; +% [rgnout] = mergeregions(cgrph,rgnout,rgn2); +% for ii=3:numel(rgns) +% rgn2 = rgns{ii}; +% [rgnout] = mergeregions(cgrph,rgnout,rgn2); +% end +% end +%%%% [in] = inpolygon(xquery,yquery,xverts,yverts) + + + +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 + +function [r,d,d2] = sinearc(t,amp,frq) +xs = t; +ys = amp*sin(frq*t); +xp = ones(size(t)); +yp = amp*frq*cos(frq*t); +xpp = zeros(size(t)); +ypp = -frq*frq*amp*sin(t); + +r = [(xs(:)).'; (ys(:)).']; +d = [(xp(:)).'; (yp(:)).']; +d2 = [(xpp(:)).'; (ypp(:)).']; +end diff --git a/devtools/test/chunkrgrphOpdimTest.m b/devtools/test/chunkrgrphOpdimTest.m index 4fb2652..38911c5 100644 --- a/devtools/test/chunkrgrphOpdimTest.m +++ b/devtools/test/chunkrgrphOpdimTest.m @@ -19,7 +19,7 @@ fchnks = {}; % this by default gives me straight lines prefs = struct('maxchunklen',0.5); -[cgrph] = chunkgraphinit(verts, edge2verts, fchnks, prefs); +[cgrph] = chunkgraph(verts, edge2verts, fchnks, prefs); % figure(1); clf; hold on; axis equal; axis off; % plot(cgrph); diff --git a/devtools/test/elastickernelsTest.m b/devtools/test/elastickernelsTest.m index 844771c..145df60 100644 --- a/devtools/test/elastickernelsTest.m +++ b/devtools/test/elastickernelsTest.m @@ -5,7 +5,7 @@ iseed = 1234; rng(iseed,'twister'); clc -addpaths_loc(); +%addpaths_loc(); % pde parameters @@ -19,9 +19,9 @@ chnkr = chunkerfunc(@(t) starfish(t),cparams); chnkr = chnkr.sort(); -f = [-1.3;2]; % charge strength -src = []; src.r = [2;1]; src.n = randn(2,1); -targ = []; targ.r = 0.1*randn(2,1); targ.n = randn(2,1); +f = [-1.3;2;0.8;-1.2]; % charge strength +src = []; src.r = [2 3.1; 1 -1.1]; src.n = randn(2,2); +targ = []; targ.r = 0.1*randn(2,4); targ.n = randn(2,4); %figure(1) %clf @@ -41,9 +41,9 @@ assert(gid_err < 1e-14); assert(min(pde_errs) < 1e-5); assert(min(pdedalt_errs) < 1e-5); -assert(min(div_errs) < 1e-9); -assert(min(trac_errs) < 1e-9); -assert(min(daltgrad_errs) < 1e-9); +assert(min(div_errs) < 1e-7); +assert(min(trac_errs) < 1e-7); +assert(min(daltgrad_errs) < 1e-7); % @@ -64,14 +64,58 @@ sigma = sysmat2\ubdry; dmatalt = chnk.elast2d.kern(lam,mu,t,targ,'dalt'); +dmataltgrad = chnk.elast2d.kern(lam,mu,t,targ,'daltgrad'); +dmatalttrac = chnk.elast2d.kern(lam,mu,t,targ,'dalttrac'); utargsol = dmatalt*(wts2.*sigma); +ugradtargsol = dmataltgrad*(wts2.*sigma); +tractargsol = dmatalttrac*(wts2.*sigma); -utarg = elasticlet(lam,mu,src,targ,f); +[utarg,tractarg,~,~,~,~,ugradtarg] = elasticlet(lam,mu,src,targ,f); err_dir = norm(utarg-utargsol)/norm(utarg); +err_dir_grad = norm(ugradtarg-ugradtargsol)/norm(ugradtarg); +err_dir_trac = norm(tractarg-tractargsol)/norm(tractarg); assert(err_dir < 1e-10); +assert(err_dir_grad < 1e-9); +assert(err_dir_trac < 1e-9); + +% + +sp = -0.5*(mu/(lam+2*mu)); stoktrac = -0.5*(lam+mu)/(lam+2*mu); +diag = sp+stoktrac; +fkern = kernel(); +fkern.eval = @(s,t) chnk.elast2d.kern(lam,mu,s,t,'strac'); +fkern.sing = 'pv'; +fkern.opdims = [2,2]; + +start = tic(); sysmat = chunkermat(chnkr,fkern); toc(start) +sysmat2 = sysmat + diag*eye(2*chnkr.npt); + +t = []; t.r = chnkr.r(:,:); +t.n = chnkr.n(:,:); t.d = chnkr.d(:,:); + +wts = weights(chnkr); wts = wts(:); +wts2 = [wts(:).'; wts(:).']; wts2 = wts2(:); + +[~,tracbdry] = elasticlet(lam,mu,src,t,f); + +sigma = gmres(sysmat2,tracbdry,[],1e-13,100); + +smat = chnk.elast2d.kern(lam,mu,t,targ,'s'); +utargsol = smat*(wts2.*sigma); + +utarg = elasticlet(lam,mu,src,targ,f); + +targperp = chnk.perp(targ.r); +nt = size(targ.r,2); +nspace = [repmat(eye(2),nt,1) targperp(:)]; +cfs = nspace\(utarg-utargsol); +utargsol = utargsol + nspace*cfs; + +err_neu = norm(utarg-utargsol)/norm(utarg); +assert(err_neu < 1e-10); function [pde_errs,pdedalt_errs,div_errs,trac_errs,gid_err,daltgrad_errs] = ... test_kernels(chnkr,s,targ,lam,mu,f,niter) @@ -94,8 +138,8 @@ smat = chnk.elast2d.kern(lam,mu,t,targ,'s'); dmat = chnk.elast2d.kern(lam,mu,t,targ,'d'); -uint = smat*(wts2.*trac) - dmat*(wts2.*u); -gid_err = norm(-ones(2,1)-uint./utarg); +uint = smat*(wts2.*trac) - dmat*(wts2.*u); uint = reshape(uint,size(utarg)); +gid_err = norm(-ones(size(utarg))-uint./utarg,'fro'); %fprintf('%5.2e Greens ID err\n',gid_err); @@ -106,7 +150,7 @@ div_errs = zeros(niter,1); trac_errs = zeros(niter,1); daltgrad_errs = zeros(niter,1); -is = randi(chnkr.npt); +is = randi(chnkr.npt,numel(f)/2,1); s = []; s.r = chnkr.r(:,is); s.d = chnkr.d(:,is); @@ -117,20 +161,21 @@ t.d = chnkr.d(:,it); t.n = chnkr.n(:,it); -[u00,trac00,dub00,dalt00,dalt00grad] = elasticlet(lam,mu,s,t,f); +[u00,trac00,dub00,dalt00,dalt00grad,dalt00trac,u00grad] = elasticlet(lam,mu,s,t,f); for i = 1:niter - h = 0.1^i; + hh = 0.1^i; + h = hh*ones(1,length(it)); t01 = t; t10 = t; t0m1 = t; tm10=t; t11 = t; t1m1 = t; tm1m1 = t; tm11=t; - t01.r(2) = t01.r(2)+h; - t10.r(1) = t10.r(1)+h; - t0m1.r(2) = t0m1.r(2)-h; - tm10.r(1) = tm10.r(1)-h; - t11.r(:) = t11.r(:)+h; - t1m1.r(:) = t1m1.r(:)+[h;-h]; - tm1m1.r(:) = tm1m1.r(:)-h; - tm11.r(:) = tm11.r(:)+[-h;h]; + t01.r(2,:) = t01.r(2,:)+h; + t10.r(1,:) = t10.r(1,:)+h; + t0m1.r(2,:) = t0m1.r(2,:)-h; + tm10.r(1,:) = tm10.r(1,:)-h; + t11.r = t11.r+h; + t1m1.r = t1m1.r+[h;-h]; + tm1m1.r = tm1m1.r-h; + tm11.r = tm11.r+[-h;h]; [u01,trac01,dub01,dalt01] = elasticlet(lam,mu,s,t01,f); [u10,trac10,dub10,dalt10] = elasticlet(lam,mu,s,t10,f); [u0m1,trac0m1,dub0m1,dalt0m1] = elasticlet(lam,mu,s,t0m1,f); @@ -140,32 +185,32 @@ [um1m1,tracm1m1,dubm1m1,daltm1m1] = elasticlet(lam,mu,s,tm1m1,f); [um11,tracm11,dubm11,daltm11] = elasticlet(lam,mu,s,tm11,f); - lapu = (u01+u10+u0m1+um10-4*u00)/h^2; - ux = (u10-um10)/(2*h); - uy = (u01-u0m1)/(2*h); - uxx = (u10+um10-2*u00)/h^2; - uyy = (u01+u0m1-2*u00)/h^2; - uxy = (u11-um11-u1m1+um1m1)/(4*h^2); + lapu = (u01+u10+u0m1+um10-4*u00)/hh^2; + ux = (u10-um10)/(2*hh); + uy = (u01-u0m1)/(2*hh); + uxx = (u10+um10-2*u00)/hh^2; + uyy = (u01+u0m1-2*u00)/hh^2; + uxy = (u11-um11-u1m1+um1m1)/(4*hh^2); uyx = uxy; - lapdub = (dub01+dub10+dub0m1+dubm10-4*dub00)/h^2; - dubx = (dub10-dubm10)/(2*h); - duby = (dub01-dub0m1)/(2*h); - dubxx = (dub10+dubm10-2*dub00)/h^2; - dubyy = (dub01+dub0m1-2*dub00)/h^2; - dubxy = (dub11-dubm11-dub1m1+dubm1m1)/(4*h^2); + lapdub = (dub01+dub10+dub0m1+dubm10-4*dub00)/hh^2; + dubx = (dub10-dubm10)/(2*hh); + duby = (dub01-dub0m1)/(2*hh); + dubxx = (dub10+dubm10-2*dub00)/hh^2; + dubyy = (dub01+dub0m1-2*dub00)/hh^2; + dubxy = (dub11-dubm11-dub1m1+dubm1m1)/(4*hh^2); dubyx = dubxy; - lapdalt = (dalt01+dalt10+dalt0m1+daltm10-4*dalt00)/h^2; - daltx = (dalt10-daltm10)/(2*h); - dalty = (dalt01-dalt0m1)/(2*h); - daltxx = (dalt10+daltm10-2*dalt00)/h^2; - daltyy = (dalt01+dalt0m1-2*dalt00)/h^2; - daltxy = (dalt11-daltm11-dalt1m1+daltm1m1)/(4*h^2); + lapdalt = (dalt01+dalt10+dalt0m1+daltm10-4*dalt00)/hh^2; + daltx = (dalt10-daltm10)/(2*hh); + dalty = (dalt01-dalt0m1)/(2*hh); + daltxx = (dalt10+daltm10-2*dalt00)/hh^2; + daltyy = (dalt01+dalt0m1-2*dalt00)/hh^2; + daltxy = (dalt11-daltm11-dalt1m1+daltm1m1)/(4*hh^2); daltyx = daltxy; - tracx = (trac10-tracm10)/(2*h); - tracy = (trac01-trac0m1)/(2*h); + tracx = (trac10-tracm10)/(2*hh); + tracy = (trac01-trac0m1)/(2*hh); pdeuh = mu*lapu + (lam+mu)*[uxx(1)+uxy(2);uyx(1)+uyy(2)]; %fprintf('%5.2e pde err\n',norm(pdeuh)/norm(u00)); @@ -179,8 +224,13 @@ %fprintf('%5.2e div err\n',norm(divtrach)/norm(trac00)) div_errs(i) = norm(divtrach)/norm(trac00); jact = [ux uy]; epsmat = 0.5*(jact + jact.'); + ugradh = jact.'; ugradh = ugradh(:); + fprintf('%5.2e sgrad err\n',norm(ugradh-u00grad)/norm(u00grad)); + daltjact = [daltx dalty]; daltepsmat = 0.5*(daltjact + daltjact.'); tracuh = (lam*(ux(1)+uy(2))*eye(2) + 2*mu*epsmat)*t.n; + tracdalt = (lam*(daltx(1)+dalty(2))*eye(2) + 2*mu*daltepsmat)*t.n; %fprintf('%5.2e trac err\n',norm(tracuh-trac00)/norm(trac00)); + fprintf('%5.2e dalttrac err\n',norm(tracdalt-dalt00trac)/norm(dalt00trac)); trac_errs(i) = norm(tracuh-trac00)/norm(trac00); daltgrad_errs(i) = ... @@ -191,18 +241,22 @@ end -function [u,trac,dub,dalt,daltgrad] = elasticlet(lam,mu,s,t,f) +function [u,trac,dub,dalt,daltgrad,dalttrac,ugrad] = elasticlet(lam,mu,s,t,f) mat = chnk.elast2d.kern(lam,mu,s,t,'s'); mattrac = chnk.elast2d.kern(lam,mu,s,t,'strac'); matdub = chnk.elast2d.kern(lam,mu,s,t,'d'); matdalt = chnk.elast2d.kern(lam,mu,s,t,'dalt'); matdaltgrad = chnk.elast2d.kern(lam,mu,s,t,'daltgrad'); +matdalttrac = chnk.elast2d.kern(lam,mu,s,t,'dalttrac'); +matgrad = chnk.elast2d.kern(lam,mu,s,t,'sgrad'); u = mat*f; trac = mattrac*f; dub = matdub*f; dalt = matdalt*f; daltgrad = matdaltgrad*f; +dalttrac = matdalttrac*f; +ugrad = matgrad*f; end diff --git a/devtools/test/flamutilitiesTest.m b/devtools/test/flamutilitiesTest.m index 86a7472..6e13270 100644 --- a/devtools/test/flamutilitiesTest.m +++ b/devtools/test/flamutilitiesTest.m @@ -31,7 +31,7 @@ end cparams = []; cparams.nover = 2; -[cgrph] = chunkgraphinit(verts,edge2verts,fchnks,cparams); +[cgrph] = chunkgraph(verts,edge2verts,fchnks,cparams); cgrph = balance(cgrph); @@ -104,16 +104,18 @@ fprintf('%5.2e s : time to build tridiag\n',t1) spmat = spmat + speye(npt); +inds = find(spmat); +spmat(inds) = sys(inds); % test matrix entry evaluator start = tic; -% opdims = [1 1]; -opdims = ones([2,1,1]); + opdims = [1 1]; +%opdims = ones([2,1,1]); +%sys2 = chnk.flam.kernbyindex(1:npt,1:npt,cgrph,kernd,opdims,spmat); sys2 = chnk.flam.kernbyindex(1:npt,1:npt,cgrph,kernd,opdims,spmat); - t1 = toc(start); fprintf('%5.2e s : time for mat entry eval on whole mat\n',t1) @@ -121,6 +123,8 @@ err2 = norm(sys2-sys,'fro')/norm(sys,'fro'); fprintf('%5.2e : fro error of build \n',err2) + assert(err2 < 1e-10); + xflam = cgrph.r(:,:); matfun = @(i,j) chnk.flam.kernbyindex(i,j,cgrph,kernd,opdims,spmat); @@ -161,7 +165,7 @@ fprintf('difference between fast-direct and iterative %5.2e\n',err) -% assert(err < 1e-10); + assert(err < 1e-10); % evaluate at targets and compare diff --git a/devtools/test/rcipTest.m b/devtools/test/rcipTest.m index 30db8e5..3de9a28 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 = chnkrtotal.wts; 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 9b8bdf1..0000000 --- a/devtools/test/refine_trefinementTest.m +++ /dev/null @@ -1,90 +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) - fprintf('oh no! 1 %d\n',i) - happy = false; - end - if (hself > 2*h2) - fprintf('oh no! 2 %d\n',i) - happy = false; - end -end - -fprintf('happy? %s\n',mat2str(happy)) - -% 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) - fprintf('oh no! 1 %d\n',i) - happy = false; - end - if (hself > 2*h2) - fprintf('oh no! 2 %d\n',i) - happy = false; - end -end - -fprintf('happy? %s\n',mat2str(happy)) diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 0000000..bfded92 --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,216 @@ +# Makefile for Sphinx documentation +# + +# You can set these variables from the command line. +SPHINXOPTS = +SPHINXBUILD = sphinx-build +PAPER = +BUILDDIR = _build + +# User-friendly check for sphinx-build +ifeq ($(shell which $(SPHINXBUILD) >/dev/null 2>&1; echo $$?), 1) +$(error The '$(SPHINXBUILD)' command was not found. Make sure you have Sphinx installed, then set the SPHINXBUILD environment variable to point to the full path of the '$(SPHINXBUILD)' executable. Alternatively you can add the directory with the executable to your PATH. If you don't have Sphinx installed, grab it from http://sphinx-doc.org/) +endif + +# Internal variables. +PAPEROPT_a4 = -D latex_paper_size=a4 +PAPEROPT_letter = -D latex_paper_size=letter +ALLSPHINXOPTS = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) . +# the i18n builder cannot share the environment and doctrees with the others +I18NSPHINXOPTS = $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) . + +.PHONY: help +help: + @echo "Please use \`make ' where is one of" + @echo " html to make standalone HTML files" + @echo " dirhtml to make HTML files named index.html in directories" + @echo " singlehtml to make a single large HTML file" + @echo " pickle to make pickle files" + @echo " json to make JSON files" + @echo " htmlhelp to make HTML files and a HTML help project" + @echo " qthelp to make HTML files and a qthelp project" + @echo " applehelp to make an Apple Help Book" + @echo " devhelp to make HTML files and a Devhelp project" + @echo " epub to make an epub" + @echo " latex to make LaTeX files, you can set PAPER=a4 or PAPER=letter" + @echo " latexpdf to make LaTeX files and run them through pdflatex" + @echo " latexpdfja to make LaTeX files and run them through platex/dvipdfmx" + @echo " text to make text files" + @echo " man to make manual pages" + @echo " texinfo to make Texinfo files" + @echo " info to make Texinfo files and run them through makeinfo" + @echo " gettext to make PO message catalogs" + @echo " changes to make an overview of all changed/added/deprecated items" + @echo " xml to make Docutils-native XML files" + @echo " pseudoxml to make pseudoxml-XML files for display purposes" + @echo " linkcheck to check all external links for integrity" + @echo " doctest to run all doctests embedded in the documentation (if enabled)" + @echo " coverage to run coverage check of the documentation (if enabled)" + +.PHONY: clean +clean: + rm -rf $(BUILDDIR)/* + +.PHONY: html +html: + $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html + @echo + @echo "Build finished. The HTML pages are in $(BUILDDIR)/html." + +.PHONY: dirhtml +dirhtml: + $(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml + @echo + @echo "Build finished. The HTML pages are in $(BUILDDIR)/dirhtml." + +.PHONY: singlehtml +singlehtml: + $(SPHINXBUILD) -b singlehtml $(ALLSPHINXOPTS) $(BUILDDIR)/singlehtml + @echo + @echo "Build finished. The HTML page is in $(BUILDDIR)/singlehtml." + +.PHONY: pickle +pickle: + $(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) $(BUILDDIR)/pickle + @echo + @echo "Build finished; now you can process the pickle files." + +.PHONY: json +json: + $(SPHINXBUILD) -b json $(ALLSPHINXOPTS) $(BUILDDIR)/json + @echo + @echo "Build finished; now you can process the JSON files." + +.PHONY: htmlhelp +htmlhelp: + $(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) $(BUILDDIR)/htmlhelp + @echo + @echo "Build finished; now you can run HTML Help Workshop with the" \ + ".hhp project file in $(BUILDDIR)/htmlhelp." + +.PHONY: qthelp +qthelp: + $(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) $(BUILDDIR)/qthelp + @echo + @echo "Build finished; now you can run "qcollectiongenerator" with the" \ + ".qhcp project file in $(BUILDDIR)/qthelp, like this:" + @echo "# qcollectiongenerator $(BUILDDIR)/qthelp/finufft.qhcp" + @echo "To view the help file:" + @echo "# assistant -collectionFile $(BUILDDIR)/qthelp/finufft.qhc" + +.PHONY: applehelp +applehelp: + $(SPHINXBUILD) -b applehelp $(ALLSPHINXOPTS) $(BUILDDIR)/applehelp + @echo + @echo "Build finished. The help book is in $(BUILDDIR)/applehelp." + @echo "N.B. You won't be able to view it unless you put it in" \ + "~/Library/Documentation/Help or install it in your application" \ + "bundle." + +.PHONY: devhelp +devhelp: + $(SPHINXBUILD) -b devhelp $(ALLSPHINXOPTS) $(BUILDDIR)/devhelp + @echo + @echo "Build finished." + @echo "To view the help file:" + @echo "# mkdir -p $$HOME/.local/share/devhelp/finufft" + @echo "# ln -s $(BUILDDIR)/devhelp $$HOME/.local/share/devhelp/finufft" + @echo "# devhelp" + +.PHONY: epub +epub: + $(SPHINXBUILD) -b epub $(ALLSPHINXOPTS) $(BUILDDIR)/epub + @echo + @echo "Build finished. The epub file is in $(BUILDDIR)/epub." + +.PHONY: latex +latex: + $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex + @echo + @echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex." + @echo "Run \`make' in that directory to run these through (pdf)latex" \ + "(use \`make latexpdf' here to do that automatically)." + +.PHONY: latexpdf +latexpdf: + $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex + @echo "Running LaTeX files through pdflatex..." + $(MAKE) -C $(BUILDDIR)/latex all-pdf + @echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex." + +.PHONY: latexpdfja +latexpdfja: + $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex + @echo "Running LaTeX files through platex and dvipdfmx..." + $(MAKE) -C $(BUILDDIR)/latex all-pdf-ja + @echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex." + +.PHONY: text +text: + $(SPHINXBUILD) -b text $(ALLSPHINXOPTS) $(BUILDDIR)/text + @echo + @echo "Build finished. The text files are in $(BUILDDIR)/text." + +.PHONY: man +man: + $(SPHINXBUILD) -b man $(ALLSPHINXOPTS) $(BUILDDIR)/man + @echo + @echo "Build finished. The manual pages are in $(BUILDDIR)/man." + +.PHONY: texinfo +texinfo: + $(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo + @echo + @echo "Build finished. The Texinfo files are in $(BUILDDIR)/texinfo." + @echo "Run \`make' in that directory to run these through makeinfo" \ + "(use \`make info' here to do that automatically)." + +.PHONY: info +info: + $(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo + @echo "Running Texinfo files through makeinfo..." + make -C $(BUILDDIR)/texinfo info + @echo "makeinfo finished; the Info files are in $(BUILDDIR)/texinfo." + +.PHONY: gettext +gettext: + $(SPHINXBUILD) -b gettext $(I18NSPHINXOPTS) $(BUILDDIR)/locale + @echo + @echo "Build finished. The message catalogs are in $(BUILDDIR)/locale." + +.PHONY: changes +changes: + $(SPHINXBUILD) -b changes $(ALLSPHINXOPTS) $(BUILDDIR)/changes + @echo + @echo "The overview file is in $(BUILDDIR)/changes." + +.PHONY: linkcheck +linkcheck: + $(SPHINXBUILD) -b linkcheck $(ALLSPHINXOPTS) $(BUILDDIR)/linkcheck + @echo + @echo "Link check complete; look for any errors in the above output " \ + "or in $(BUILDDIR)/linkcheck/output.txt." + +.PHONY: doctest +doctest: + $(SPHINXBUILD) -b doctest $(ALLSPHINXOPTS) $(BUILDDIR)/doctest + @echo "Testing of doctests in the sources finished, look at the " \ + "results in $(BUILDDIR)/doctest/output.txt." + +.PHONY: coverage +coverage: + $(SPHINXBUILD) -b coverage $(ALLSPHINXOPTS) $(BUILDDIR)/coverage + @echo "Testing of coverage in the sources finished, look at the " \ + "results in $(BUILDDIR)/coverage/python.txt." + +.PHONY: xml +xml: + $(SPHINXBUILD) -b xml $(ALLSPHINXOPTS) $(BUILDDIR)/xml + @echo + @echo "Build finished. The XML files are in $(BUILDDIR)/xml." + +.PHONY: pseudoxml +pseudoxml: + $(SPHINXBUILD) -b pseudoxml $(ALLSPHINXOPTS) $(BUILDDIR)/pseudoxml + @echo + @echo "Build finished. The pseudo-XML files are in $(BUILDDIR)/pseudoxml." diff --git a/docs/_static/.keep b/docs/_static/.keep new file mode 100644 index 0000000..e69de29 diff --git a/docs/_templates/.keep b/docs/_templates/.keep new file mode 100644 index 0000000..e69de29 diff --git a/docs/assets/images/scatter_star_trim.png b/docs/assets/images/scatter_star_trim.png new file mode 100644 index 0000000..a53f651 Binary files /dev/null and b/docs/assets/images/scatter_star_trim.png differ diff --git a/docs/conf.py b/docs/conf.py index c7eff8a..aa90e8e 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -55,16 +55,16 @@ # built documents. # # The short X.Y version. -version = '0.0' +version = '1.0' # The full version, including alpha/beta/rc tags. -release = '0.0.0' +release = '1.0.0' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. # # This is also used if you do content translation via gettext catalogs. # Usually you set "language" from the command line for these cases. -language = None +language = 'en' # List of patterns, relative to source directory, that match files and # directories to ignore when looking for source files. @@ -83,7 +83,7 @@ # The theme to use for HTML and HTML Help pages. See the documentation for # a list of builtin themes. # -html_theme = 'alabaster' +html_theme = 'sphinx_rtd_theme' # Theme options are theme-specific and customize the look and feel of a theme # further. For a list of options available for each theme, see the diff --git a/docs/gallery.rst b/docs/gallery.rst new file mode 100644 index 0000000..5e16af8 --- /dev/null +++ b/docs/gallery.rst @@ -0,0 +1,6 @@ +Gallery +========== + +.. toctree:: + :maxdepth: 2 + :caption: Contents: diff --git a/docs/getchunkie.rst b/docs/getchunkie.rst new file mode 100644 index 0000000..82d2c23 --- /dev/null +++ b/docs/getchunkie.rst @@ -0,0 +1,53 @@ +Get Chunkie +============ + +Installation Instructions +-------------------------- + +Clone the repository with the submodules + +.. code:: bash + + git clone --recurse-submodules https://github.com/fastalgorithms/chunkie.git + +and run startup.m in the install directory. +This will download the FLAM and fmm2d submodules, include FLAM in +the matlab path, and generate the fmm2d mex file if a fortran compiler +exists. + + +Troubleshooting +----------------- + +- The fmm2d mex installation is currently not supported on Windows, to + complete the mex installation, follow instructions on the `fmm2d documentation `_ +- fmm2d mex installation depends on gfortran. In case a compiler is not + found, the installation will be skipped. To install dependencies follow the procedure below based on your OS + + * MacOS + + Get xcode, command line tools by running + + .. code:: bash + + xcode-select --install + + Then install Homebrew from https://brew.sh, and finally install gfortran using + + .. code:: bash + + brew install gcc + + * Ubuntu linux + + .. code:: bash + + sudo apt-get install make build-essential gfortran + + * Fedora/centOS linux + + .. code:: bash + + sudo yum install make gcc gcc-c++ gcc-gfortran libgomp + +- If installing without submodules, chunkie depends on `FLAM `_, and optionally on the `fmm2d `_ repository. Parts of the library will not function without FLAM and its subdirectories included in the matlab path. diff --git a/docs/guide.rst b/docs/guide.rst new file mode 100644 index 0000000..b24e3e2 --- /dev/null +++ b/docs/guide.rst @@ -0,0 +1,25 @@ +Chunkie User Guide +=================== + +Detailed documentation of all high-level functions and +classes can be obtained by viewing the source file or by +typing "help [function name]" in a MATLAB command window: + +.. code:: matlab + + help chunker + + + +.. toctree:: + :maxdepth: 1 + :caption: Contents: + + guide/chunkers + guide/simplebvps + guide/largerproblems + guide/chunkgraphs + guide/chunkgraphbvps + guide/transmission + + diff --git a/docs/guide/chunkers.rst b/docs/guide/chunkers.rst new file mode 100644 index 0000000..d306bf5 --- /dev/null +++ b/docs/guide/chunkers.rst @@ -0,0 +1,153 @@ + +.. role:: matlab(code) + :language: matlab + +Curve Discretization with Chunkers +=================================== + +In chunkie, a smooth, regular curve is discretized by dividing +it into pieces, called "chunks", which are then represented +by polynomial interpolants at scaled Legendre nodes. This +information is stored in a :matlab:`chunker` object. + +Creating a Chunker +------------------- + +A chunker object can be obtained from a curve parameterization +using the function :matlab:`chunkerfunc`. The :matlab:`chunker` +class overloads some MATLAB commands, like :matlab:`plot` and +:matlab:`quiver`, to simplify common visualization tasks. The +code below creates a :matlab:`chunker` object for a circle and +plots the geometry and the normal vectors. + +.. include:: ../../chunkie/guide/guide_chunkers.m + :literal: + :code: matlab + :start-after: % START CIRCLE + :end-before: % END CIRCLE + +.. image:: ../../chunkie/guide/guide_chunkers_circle.png + :width: 350px + :alt: circle chunker + :align: center + +.. note:: + + The chunkerfunc routine expects a specific format for the curve + parameterization function. In particular, if :matlab:`fcurve` + is a MATLAB function defining the curve parameterization and + :matlab:`t` is a vector of points in parameter space, then + :matlab:`r = fcurve(t)` should be a matrix in which :matlab:`r(:,j)` + contains the :math:`(x,y)` coordinates of the point on the + curve corresponding :matlab:`t(j)`. + + If the output of :matlab:`fcurve` + is of the form :matlab:`[r,d] = fcurve(t)` or + :matlab:`[r,d,d2] = fcurve(t)`, then :matlab:`d` and :matlab:`d2` + are assumed to be first and second derivatives, respectively, + of the position with respect to the underlying parameterization. + Providing these values can improve the convergence rate and + precision achieved in integral equation calculations. + +Some utilities are provided that define a family of +parameterized curves: + +.. include:: ../../chunkie/guide/guide_chunkers.m + :literal: + :code: matlab + :start-after: % START MORE PARAMS + :end-before: % END MORE PARAMS + +.. |pic1| image:: ../../chunkie/guide/guide_chunkers_bymode.png + :width: 300px + :alt: by mode chunker + +.. |pic2| image:: ../../chunkie/guide/guide_chunkers_starfish.png + :width: 300px + :alt: starfish chunker + +|pic1| |pic2| + +Given a set of vertices, a rounded polygon can be defined: + +.. include:: ../../chunkie/guide/guide_chunkers.m + :literal: + :code: matlab + :start-after: % START ROUNDED POLY + :end-before: % END ROUNDED POLY + +.. image:: ../../chunkie/guide/guide_chunkers_barbell.png + :width: 500px + :alt: circle chunker + :align: center + + +Working with Chunkers +---------------------- + +Instances of :matlab:`chunker` objects can be manipulated in +several ways. +Users are free to update the position and derivative +fields of :matlab:`chunker` objects, though the software will +not check whether the user has done this in a consistent +manner. + +The example below takes the circle and random mode domains +created above and creates a new domain from them with +multiple components. The random mode domain is rotated and +the circle is shifted and its orientation is reversed, so +that the normal for the combined domain consistently points out +of the interior. + +.. include:: ../../chunkie/guide/guide_chunkers.m + :literal: + :code: matlab + :start-after: % START SHIFT AND REVERSE + :end-before: % END SHIFT AND REVERSE + +.. image:: ../../chunkie/guide/guide_chunkers_shiftandreverse.png + :width: 500px + :alt: combined chunker + :align: center + +It is also possible to find the points on the interior of a +:matlab:`chunker` object, which is convenient in many plotting +tasks: + +.. include:: ../../chunkie/guide/guide_chunkers.m + :literal: + :code: matlab + :start-after: % START INTERIOR + :end-before: % END INTERIOR + +.. image:: ../../chunkie/guide/guide_chunkers_interior.png + :width: 500px + :alt: interior points + :align: center + +.. note:: + + In the example above, it is crucial that the orientation + of the circle was reversed so that the combined domain + had consistent normals, i.e. so that the normals + all pointed out of the interior. Without doing this, + the :matlab:`chunkerinterior` function will fail. + +There is more available. The :matlab:`chunker` class documentation +gives a survey of the available methods: + +.. include:: ../../chunkie/@chunker/chunker.m + :literal: + :code: matlab + :start-after: classdef chunker + :end-before: % author + +.. note:: + + To obtain the documentation of a class method which has + overloaded the name of a MATLAB built-in, use the syntax + :matlab:`help class_name/method_name`. For example: + + .. code:: matlab + + help chunker/area diff --git a/docs/guide/chunkgraphbvps.rst b/docs/guide/chunkgraphbvps.rst new file mode 100644 index 0000000..a61a095 --- /dev/null +++ b/docs/guide/chunkgraphbvps.rst @@ -0,0 +1,25 @@ + +Boundary Value Problems on Chunkgraphs +========================================= + + +A Sound Hard Scattering Problem with Corners +--------------------------------------------- + + + +An Open Arc Problem +--------------------- + + + +Mixed Boundary Conditions +-------------------------- + + + + + + + + diff --git a/docs/guide/chunkgraphs.rst b/docs/guide/chunkgraphs.rst new file mode 100644 index 0000000..576f4b0 --- /dev/null +++ b/docs/guide/chunkgraphs.rst @@ -0,0 +1,3 @@ + +Describing More Complicated Geometries with Chunkgraphs +======================================================== diff --git a/docs/guide/largerproblems.rst b/docs/guide/largerproblems.rst new file mode 100644 index 0000000..f706bc1 --- /dev/null +++ b/docs/guide/largerproblems.rst @@ -0,0 +1,3 @@ + +Larger Boundary Value Problems +================================ diff --git a/docs/guide/simplebvps.rst b/docs/guide/simplebvps.rst new file mode 100644 index 0000000..c5d70fc --- /dev/null +++ b/docs/guide/simplebvps.rst @@ -0,0 +1,269 @@ + +.. role:: matlab(code) + :language: matlab + +Solving Standard Boundary Value Problems +========================================= + +Chunkie uses the integral equation method to solve PDEs. +In contrast with finite element methods, where the PDE is +typically converted to a weak form, the standard approach for +an integral equation method is to select an integral representation +of the solution that satisfies the PDE inside the domain +*a priori*. This cannot be done for just any PDE; typically, +integral equation methods are limited to constant coefficient, +homogeneous PDEs. But when the methods apply, they are +efficient and accurate. + +The following is a very high level discussion of the main ideas. +Some more specific examples are included further below. +For more, see, *inter alia*, [Rok1985]_, [CK2013]_, [GL1996]_, and +[Poz1992]_. + +Let :math:`\Omega` be a domain with boundary :math:`\Gamma`. +Suppose that we are interested in the solution of the boundary +value problem + +.. math:: + + \begin{align*} + \mathcal{L} u &= 0 & \textrm{ in } \Omega \; ,\\ + \mathcal{B} u &= f & \textrm{ on } \Gamma \; , + \end{align*} + +where :math:`\mathcal{L}` is a linear, constant coefficient PDE +operator and :math:`\mathcal{B}` is a linear trace operator that +restricts :math:`u`, or some linear operator applied to :math:`u`, to the +boundary. + +Suppose that :math:`G` is the Green function for :math:`\mathcal{L}` +and let :math:`K` be an integral kernel that is a linear +operator applied to :math:`G`. Then, :math:`u` can be represented +as a *layer potential*, i.e. it can be written in the form + +.. math:: + + u(x) = \int_\Gamma K(x,y) \sigma(y) \, ds(y) \; , + +which has the property that :math:`u` automatically solves the PDE +for any choice of density :math:`\sigma` (defined on the boundary +alone). The idea is to then solve for a specific :math:`\sigma` +for which the boundary condition holds: + +.. math:: + + \mathcal{B} \left [ \int_\Gamma K(\cdot,y)\sigma(y)\, ds(y) \right ] + = f \textrm{ on } \Gamma \; . + +Some of the art in the design of an integral equation method is +in choosing :math:`K` so that the above yields a reasonable boundary +integral equation for :math:`\sigma`. Fortunately, there is a +large literature on good quality integral representations for the most +common problems in classical physics (linear elasticity, electrostatics, +acoustics, and fluid flow). We include some examples below. + +This integral equation can then be discretized. Chunkie employs a +Nyström discretization method, i.e. the unknown :math:`\sigma` +is represented by its values at the nodes of a :matlab:`chunker` +object. Because :math:`K` is defined in terms of the Green function +for :math:`\mathcal{L}` it often has mild (or even strong) singularities +when restricted to :math:`\Gamma`. One of the key features of +chunkie is that the function :matlab:`layermat` will automatically return +a Nyström discretization of a layer potential using high-order +accurate quadrature rules [BGR2010]_. + +After solving the discrete system for :math:`\sigma`, the PDE solution, +:math:`u`, can then be recovered by evaluating the layer potential +representation at any points of interest. This integral is nearly singular +for points near the boundary. Chunkie provides the function +:matlab:`layereval` for doing this accurately, employing adaptive +integration when necessary. + + +A Laplace Problem +------------------ + +Consider the Laplace Neumann boundary value problem: + +.. math:: + + \begin{align*} + \Delta u &= 0 & \textrm{ in } \Omega \; ,\\ + \frac{\partial u}{\partial n} &= f & \textrm{ on } \Gamma \; , + \end{align*} + +where :math:`n` is the outward normal. +This is a model for the equilibrium heat distribution in a +homogeneous medium, where :math:`f` is a prescribed heat +flux along the boundary. + +For a solution to exist, :math:`f` must satisfy the compatibility +condition: + +.. math:: + + \int_\Gamma f(y) \, ds(y) = 0 \; , + +which is equivalent to requiring that there is no net flux +of heat into the body. +When solutions exist, they are not unique because adding a +constant to any solution gives another solution. + + +The solution of this problem by layer potentials has been +well understood for some time. The standard method is reviewed +here with a fair amount of detail. Users familiar with the +approach may simply want to skip down to the code example. + +The Green function for the Laplace equation is + +.. math:: + + G_0 (x,y) = -\frac{1}{2\pi} \log |x-y| + +The standard choice for the layer potential for this problem is +the *single layer potential* + +.. math:: + + u(x) = [S\sigma](x) := \int_\Gamma G_0(x,y) \sigma(y) ds(y) \; . + +Then, imposing the boundary conditions on this representation +results in the equation + +.. math:: + + \begin{align*} + f(x_0) &= \lim_{x\in \Omega, x\to x_0} n(x_0)\cdot \nabla_x + \int_\Gamma G_0(x,y) \sigma(y) \, ds(y) \\ + &= \frac{1}{2} \sigma(x_0) + P.V. \int_\Gamma n(x_0) \cdot \nabla_x G_0(x_0,y) + \sigma(y) \, ds(y) \\ + &=: \left [\left ( \frac{1}{2} \mathcal{I} + \mathcal{S}' \right ) \sigma + \right ] (x_0) \; , + \end{align*} + +where :math:`P.V.` indicates that the integral should be interpreted +in the principal value sense, :math:`\mathcal{I}` is the identity operator, and +the calligraphic letters indicate that these operators act +on functions defined on :math:`\Gamma`. The second equality above uses +a so-called "jump condition" for the single layer potential. + +It turns out that the equation for :math:`\sigma` above is a second-kind +integral equation, i.e. of the form "identity plus compact", but it +is not invertible. This has a simple remedy. The equation + +.. math:: + + f = \left ( \frac{1}{2} \mathcal{I} + \mathcal{S}' + \mathcal{W} \right) \sigma \; , + +where + +.. math:: + + [\mathcal{W}\sigma](x) = \int_\Gamma \sigma (y) \, ds(y) \; , + +has the same solutions as the original equation (so long as +:math:`f` satisfies the compatibility condition) and is +invertible. The operator :math:`\mathcal{W}` is straightforward +to discretize: + +.. math:: + + W = \begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{bmatrix} + \begin{bmatrix} w_1 & w_2 & \cdots & w_n \end{bmatrix} \; , + +where :math:`w_1,\ldots,w_n` are scaled integration weights on :math:`\Gamma`. +For historical reasons, the routine that returns this matrix +in chunkie is called :matlab:`onesmat`. + +To discretize the boundary integral equation, chunkie requires a +:matlab:`chunker` class object discretization of the boundary and a +:matlab:`kernel` class object representing the integral kernel. +Many of the most common kernel types are available in chunkie, +including the normal derivative of the single layer kernel, which is +called "sprime". This can be obtained via + +.. include:: ../../chunkie/guide/guide_simplebvps.m + :literal: + :code: matlab + :start-after: % START SPRIME + :end-before: % END SPRIME + +The most essential thing that an instance of the kernel class tells chunkie +is how to evaluate the kernel function. This is stored in the field +:matlab:`obj.eval`. The object can also store other information about +the kernel, such as fast algorithms for its evaluation or the type +of singularity it has. + +.. note:: + + The :matlab:`obj.eval` function in a kernel is expected to be + of the form :matlab:`eval(s,t)`, where :matlab:`s` and :matlab:`t` + are structs that specify information about the "sources" and "targets" + for which the kernel is being evaluated. This has the unfortunate feature + that :math:`x` and :math:`y` take the roles of "target" and "source", + respectively, in the math notation above, so the order is reversed. + + The structs :matlab:`s` and :matlab:`t` must have the fields + :matlab:`s.r` and :matlab:`t.r` which specify the + locations of the sources and targets, respectively. + For sources/targets on a :matlab:`chunker` + object, these will also contain fields :matlab:`s.d`, :matlab:`s.d2`, + :matlab:`s.n`, and :matlab:`s.data` (likewise for :matlab:`t`) containing + the chunk information for the corresponding points. For example, + the kernel :matlab:`kernsp` defined above is only defined for targets + on a :matlab:`chunker` discretization of a curve and the kernel evaluator + :matlab:`kernsp.eval` assumes that the field :matlab:`t.n` is populated + with the normal vectors at the targets. + +Because these standard kernels are built into chunkie, this PDE +can be solved with very little coding: + +.. include:: ../../chunkie/guide/guide_simplebvps.m + :literal: + :code: matlab + :start-after: % START LAPLACE NEUMANN + :end-before: % END LAPLACE NEUMANN + + +.. image:: ../../chunkie/guide/guide_simplebvps_laplaceneumann.png + :width: 500px + :alt: combined chunker + :align: center + + +A Helmholtz Scattering Problem +------------------------------- + + + + +A Stokes Flow Problem +---------------------- + + +References +------------ + +.. [Rok1985] Rokhlin, Vladimir. "Rapid solution of integral equations + of classical potential theory." Journal of Computational + Physics 60.2 (1985): 187-207. + + +.. [Poz1992] Pozrikidis, Constantine. *Boundary integral and singularity + methods for linearized viscous flow*. Cambridge University + Press, 1992. + + +.. [GL1996] Guenther, Ronald B., and John W. Lee. *Partial differential + equations of mathematical physics and integral equations*. + Courier Corporation, 1996. + +.. [CK2013] Colton, David, and Rainer Kress. *Integral equation + methods in scattering theory*. Society for Industrial + and Applied Mathematics, 2013. + +.. [BGR2010] Bremer, James, Zydrunas Gimbutas, and Vladimir Rokhlin. + "A nonlinear optimization procedure for generalized Gaussian + quadratures." SIAM Journal on Scientific Computing 32.4 (2010): + 1761-1788. diff --git a/docs/guide/transmission.rst b/docs/guide/transmission.rst new file mode 100644 index 0000000..aa5169d --- /dev/null +++ b/docs/guide/transmission.rst @@ -0,0 +1,3 @@ + +Application to Transmission Problems with Multiple Junctions +============================================================= diff --git a/docs/index.rst b/docs/index.rst index c1022a6..ce8bbec 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -3,14 +3,17 @@ You can adapt this file completely to your liking, but it should at least contain the root `toctree` directive. +.. role:: matlab(code) + :language: matlab chunkie -======== - -.. toctree:: - :maxdepth: 2 - :caption: Contents: +========== +.. image:: assets/images/scatter_star_trim.png + :width: 350 px + :alt: scattering by a smooth star domain + :align: center + chunkie is a MATLAB package for solving partial differential equations in two dimensional domains. The name "chunkie" derives @@ -19,27 +22,95 @@ discretize the domain boundary and the use of an integral equation method to solve the PDE. With chunkie, you can solve the Helmholtz equation on an exterior -domain in a few lines of MATLAB: - - import project - # Get your stuff done - project.do_stuff() +domain in a few lines of MATLAB. The code below defines and solves +a scattering problem, producing the image above: + +.. code:: matlab + + % planewave definitions + + kvec = 20*[1;-1.5]; + zk = norm(kvec); + planewave = @(kvec,r) exp(1i*sum(bsxfun(@times,kvec(:),r(:,:)))).'; + + % discretize domain + + chnkr = chunkerfunc(@(t) starfish(t,narms,amp),struct('maxchunklen',4/zk)); + + % build CFIE and solve + + fkern = kernel('helm','c',zk,[1,-zk*1i]); + sysmat = chunkermat(chnkr,fkern); + sysmat = 0.5*eye(chnkr.k*chnkr.nch) + sysmat; + + rhs = -planewave(kvec(:),chnkr.r(:,:)); + sol = gmres(sysmat,rhs,[],1e-13,100); + + % evaluate at targets + + x1 = linspace(-3,3,400); + [xxtarg,yytarg] = meshgrid(x1,x1); + targets = [xxtarg(:).';yytarg(:).']; + + in = chunkerinterior(chnkr,targets); + out = ~in; + + uscat = chunkerkerneval(chnkr,fkern,sol,targets(:,out)); + + uin = planewave(kvec,targets(:,out)); + utot = uscat(:)+planewave(kvec,targets(:,out)); + + % plot + + maxu = max(abs(utot(:))); + figure() + zztarg = nan(size(xxtarg)); + zztarg(out) = utot; + h=pcolor(xxtarg,yytarg,imag(zztarg)); + set(h,'EdgeColor','none') + hold on + plot(chnkr,'LineWidth',2) + axis equal tight + colormap(redblue) + caxis([-maxu,maxu]) Features -------- -- easy-to-use: includes pre-packaged solvers for Laplace, Helmholtz, - and Stokes equations and convenient tools for visualizing - domains and plotting solutions +- easy-to-use: most chunkie functionality only requires that + you specify two things, the geometry (as a "chunker" object) + and the integral kernel (as a "kernel" object). The distribution + includes pre-defined kernels for the Laplace, Helmholtz, + and Stokes PDEs and convenient tools for visualizing domains + and plotting solutions. +- capable: chunkie uses high order accurate quadrature rules + and fast algorithms to provide accurate solutions with linear + scaling. The chunkie utilities can be used to solve singular + problems on domains with corners and multiple material junctions + with relative ease. - flexible: allows you to define your own integral kernels and - use the quadrature routines in a modular fashion -- fast: works with FLAM (fast linear algebra in MATLAB) - for fast, direct integral equation solves + use the quadrature routines in a modular fashion. -Installation ------------- +How to Get Chunkie +------------------- -chunkie is pure MATLAB. You can download ... +- `Get chunkie `_ (includes installation instructions + and troubleshooting guide) +- `Source Code `_ + + +Using Chunkie +---------------- + +The `chunkie guide `_ +provides an overview of the chunkie approach to PDEs and a detailed +look at chunkie's capabilities. + +Gallery +-------- + +The `gallery `_ has examples of chunkie applications +submitted by users. Contributing ------------ @@ -51,6 +122,7 @@ We welcome your collaboration! First, see the - `Source Code `_ + License ------- @@ -58,9 +130,14 @@ The project is licensed under a modified BSD 3-clause license. -Indices and tables -================== -* :ref:`genindex` -* :ref:`modindex` -* :ref:`search` + +.. toctree:: + :maxdepth: 1 + :caption: Contents + + getchunkie + guide + gallery + + diff --git a/docs/requirements.txt b/docs/requirements.txt new file mode 100644 index 0000000..483a4e9 --- /dev/null +++ b/docs/requirements.txt @@ -0,0 +1 @@ +sphinx_rtd_theme