From aabbb0947088ebd81101da91d478c91c05ebf3b2 Mon Sep 17 00:00:00 2001 From: Travis Askham Date: Sat, 19 Oct 2024 17:47:14 -0400 Subject: [PATCH] clean up chunkerinterior file, simplify pquadTest --- chunkie/chunkerinterior.m | 37 +----------------------------- devtools/test/pquadTest.m | 48 +++++++++------------------------------ 2 files changed, 12 insertions(+), 73 deletions(-) diff --git a/chunkie/chunkerinterior.m b/chunkie/chunkerinterior.m index 91c880f..c5878ad 100644 --- a/chunkie/chunkerinterior.m +++ b/chunkie/chunkerinterior.m @@ -130,41 +130,6 @@ 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 - -% looking for closest point might fail if an inaccurate panel for a target -% is not the closest panel to that target -% here we look for other boundary pieces that are close each panel - -%flag2 = flagnear_rectangle(chnkr,chnkr.r(:,:),optsflag); - -% ignore self and adjacent -% inds = 1:chnkr.npt; -% slf = repmat(1:chnkr.nch,chnkr.k,1); slf = slf(:); -% left = chnkr.adj(1,slf); right = chnkr.adj(2,slf); -% indsl = inds(left > 0); left = left(left> 0); -% indsr = inds(right > 0); right = right(right >0); -% ii = [inds(:); indsl(:); indsr(:)]; jj = [slf(:); left(:); right(:)]; -% linind = sub2ind([chnkr.npt,chnkr.nch],ii,jj); -%flag2(linind) = 0; - -% for any target that was in the inaccurate region of a panel, also need to -% check any panel that is similarly close to that panel for closest normal -% for j = 1:chnkr.nch -% ibdrypt = find(flag2(:,j)); -% ich = unique(slf(ibdrypt)); -% ipt = find(flag(:,j)); -% -% ii = repmat(ipt(:),length(ich),1); -% jj = repmat(ich(:).',length(ipt),1); jj = jj(:); -% linind = sub2ind(size(flag),ii,jj); -% flag(linind) = 1; -% end % do Cauchy integral with oversampling @@ -229,4 +194,4 @@ yv = [yv; nan; yt(:)]; end -in(~ingood) = inpolygon(pts(1,~ingood),pts(2,~ingood),xv,yv); \ No newline at end of file +in(~ingood) = inpolygon(pts(1,~ingood),pts(2,~ingood),xv,yv); diff --git a/devtools/test/pquadTest.m b/devtools/test/pquadTest.m index aa6d9a1..e72a1d0 100644 --- a/devtools/test/pquadTest.m +++ b/devtools/test/pquadTest.m @@ -23,44 +23,15 @@ zk = norm(kvec); -% define geometry and boundary conditions -% (vertices defined by another function) -verts = chnk.demo.barbell(2.0,2.0,1.0,1.0); -nv = size(verts,2); -edgevals = randn(1,nv); - -% parameters for curve rounding/chunking routine to oversample boundary -cparams = []; -cparams.widths = 0.1*ones(size(verts,2),1);% width to cut about each corner -cparams.eps = 1e-12; % tolerance at which to resolve curve -cparams.rounded = true; - -% call smoothed-polygon chunking routine -% a smoothed version of edgevals is returned in -% chnkr.data -chnkr = chunkerpoly(verts,cparams,[],edgevals); -chnkr = refine(chnkr,struct('nover',1)); -% -% plot smoothed geometry and data - -figure(1) -clf -plot(chnkr,'-x') -hold on -quiver(chnkr) -axis equal - -figure(2) -clf -nchplot = 1:chnkr.nch; -plot3(chnkr,1) +cparams = []; cparams.maxchunklen = 2.0/real(zk); +chnkr = chunkerfunc(@(t) starfish(t),cparams); % solve and visualize the solution % build laplace dirichlet matrix % fkern = kernel('laplace','d'); -fkern = kernel('helmholtz','c',zk,[1.5,-2i]); +fkern = kernel('helmholtz','c',zk,[1,-2i]); opts = []; start = tic; C = chunkermat(chnkr,fkern,opts); t1 = toc(start); @@ -72,7 +43,7 @@ sys = -0.5*eye(chnkr.npt) + C; % rhs = chnkr.data(1,:); rhs = rhs(:); -rhs = besselh(0,zk*abs((1+1.25*1i)-(chnkr.r(1,:)+1i*chnkr.r(2,:)))); rhs = rhs(:); +rhs = besselh(0,zk*abs((2+2.25*1i)-(chnkr.r(1,:)+1i*chnkr.r(2,:)))); rhs = rhs(:); start = tic; sol = gmres(sys,rhs,[],1e-14,100); t1 = toc(start); fprintf('%5.2e s : time for dense gmres\n',t1) @@ -96,7 +67,7 @@ opts.forcepquad = true; opts.side = 'i'; % 'i' for interior, 'e' for exterior, for positively oriented curve. start = tic; -Csolpquad = chunkerkerneval(chnkr,fkern,sol,targets(:,in),opts); +tic; Csolpquad = chunkerkerneval(chnkr,fkern,sol,targets(:,in),opts); % below check slp and dlp % fkernd = kernel('helmholtz','d',zk); % fkerns = kernel('helmholtz','s',zk); @@ -106,14 +77,17 @@ t1 = toc(start); fprintf('%5.2e s : time for kerneval (Helsing-Ojala for near)\n',t1); +soltrue = besselh(0,zk*abs((2+2.25*1i)-(targets(1,in)+1i*targets(2,in)))); + start = tic; Csol = chunkerkerneval(chnkr,fkern,sol,targets(:,in)); t1 = toc(start); fprintf('%5.2e s : time for kerneval (adaptive for near)\n',t1); % Compare with reference solution Dsol -error = max(abs(Csol-Csolpquad))/max(abs(Csol)); -fprintf('%5.2e : Relative max error\n',error); +err = max(abs(Csol-Csolpquad))/max(abs(Csol)); +err2 = max(abs(soltrue(:)-Csolpquad(:)))/max(abs(soltrue(:))); +fprintf('%5.2e : Relative max error\n',err); % -assert(error < 1e-10) +assert(err < 1e-10) % profile viewer \ No newline at end of file