clean up chunkerinterior file, simplify pquadTest

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);
+in(~ingood) = inpolygon(pts(1,~ingood),pts(2,~ingood),xv,yv);

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