Patch chunkerinterior

chunkie/chunkerinterior.m

@@ -130,21 +130,15 @@
 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
+
+% do Cauchy integral with oversampling 
 
 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
@@ -182,61 +176,22 @@
 end
 in = abs(vals1+1) < abs(vals1);
 
-% 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
+% for points where integral rule is less definitive, use inpolygon
 
-    % 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
+ingood = or(abs(vals1+1) < 0.05,abs(vals1) < 0.05);
 
-for i = 1:length(ipt)
-    jj = ipt(i);
-    in(jj) = (rss(:,jj)-pts(:,jj)).'*[dss(2,jj);-dss(1,jj)] > 0;
+[inds,~,info] = sortinfo(chnkr);
+rr = chnkr.r(:,:,inds(:));
+ncomp = info.ncomp;
+nchs = info.nchs;
+ladr = cumsum([1;nchs(:)]);
+
+xv = []; yv = [];
+for j = 1:ncomp
+    xt = rr(1,:,ladr(j):(ladr(j+1)-1));
+    yt = rr(2,:,ladr(j):(ladr(j+1)-1));
+    xv = [xv; nan; xt(:)];
+    yv = [yv; nan; yt(:)];
 end
+
+in(~ingood) = inpolygon(pts(1,~ingood),pts(2,~ingood),xv,yv);

chunkie/chunkgraphinregion.m

@@ -77,33 +77,38 @@
 end
 
 nr = numel(cg.regions);
+
 for ir = 1:nr
     ncomp = numel(cg.regions{ir});
-    intmp = zeros(npts,1);
+    chnkrscomp(ncomp) = chunker(p,t,w);
     for ic = 1:ncomp
         edgelist = cg.regions{ir}{ic};
         nedge = numel(edgelist);
         chnkrs(nedge) = chunker(p,t,w);
         for ie = 1:nedge
             eid = edgelist(ie);
             if eid > 0
-                if ir == 1
-                    chnkrs(ie) = cg.echnks(eid);
-                else
-                    chnkrs(ie) = reverse(cg.echnks(eid));
-                end
+                chnkrs(ie) = cg.echnks(eid);
             else
-                if ir == 1
-                    chnkrs(ie) = reverse(cg.echnks(-eid));
-                else
-                    chnkrs(ie) = cg.echnks(-eid);
-                end
+                chnkrs(ie) = sort(reverse(cg.echnks(-eid)));
             end
         end
-        intmp = intmp + reshape(chunkerinterior(merge(chnkrs(1:nedge)),ptsobj,opts),npts,1);
+        nchs = [chnkrs(1:nedge).nch];
+        ladr = cumsum([1, nchs(:).']);
+        chnkrtmp = merge(chnkrs(1:nedge));
+        chnkrtmp.adj(1,1) = ladr(end)-1;
+        chnkrtmp.adj(2,end) = 1;
+        for jj = 1:(nedge-1)
+            ich1 = ladr(jj+1)-1; ich2 = ladr(jj+1); 
+            chnkrtmp.adj(1,ich2) = ich1;
+            chnkrtmp.adj(2,ich1) = ich2;
+        end
+        if ir ~= 1
+            chnkrtmp = reverse(chnkrtmp);
+        end
+        chnkrscomp(ic) = chnkrtmp;
     end
-
-    intmp = intmp > 0;
+    intmp = reshape(chunkerinterior(merge(chnkrscomp(1:ncomp)),ptsobj,opts),npts,1);
     if ir == 1
         ids(~intmp) = ir;
     else

devtools/test/chunkerinteriorTest.m

@@ -5,7 +5,6 @@
 clearvars; close all;
 seed = 8675309;
 rng(seed);
-addpaths_loc();
 
 doadap = false;
 
@@ -89,4 +88,41 @@
 targs = starfish(ttarg).*scal;
 
 in = chunkerinterior(chnkr,targs,struct('axissym',true));
-assert(all(in(:) == (scal(:) <= 1)));
+assert(all(in(:) == (scal(:) <= 1)));
+
+%% 
+
+% a stress test. previously this triggered a bug for the old version which
+% used the normals to determine interior points
+
+amp = 0.25;
+scale = .3;
+
+ctr = [-2;-1.6];
+chnkr_int= chunkerfunc(@(t) starfish(t,3,amp,ctr,pi/4,scale)); 
+chnkr_int = sort(reverse(chnkr_int));
+
+a = max(vecnorm(chnkr_int.r(:,:)))*1.01;
+chnkr_ext = chunkerfunc(@(t) [a*cos(t(:).'); a*sin(t(:).')]); 
+chnkr = merge([chnkr_ext,chnkr_int]);
+
+% return
+
+L = max(abs(chnkr.r),[],"all");
+x1 = linspace(-L,L,100);
+[xx,yy] = meshgrid(x1,x1);
+targs = [xx(:).'; yy(:).']; 
+xv=[chnkr_ext.r(1,:),nan,chnkr_int.r(1,:)];
+yv=[chnkr_ext.r(2,:),nan,chnkr_int.r(2,:)];
+tic; in0 = inpolygon(xx(:),yy(:),xv,yv); toc
+tic; in = chunkerinterior(chnkr,{x1,x1}); toc
+
+% clf
+% plot(chnkr,'g-o')
+% hold on
+% scatter(xx(in(:)),yy(in(:)),'bo')
+% scatter(xx(~in(:)),yy(~in(:)),'rx')
+
+assert(all(in0 == in));
+
+

devtools/test/mixedbcTest.m

@@ -93,14 +93,9 @@
 targets = zeros(2,length(xxtarg(:)));
 targets(1,:) = xxtarg(:); targets(2,:) = yytarg(:);
 
-start = tic; in1 = chunkerinterior(merge(cgrph.echnks(edir)),{xtarg,ytarg});
-in2 = chunkerinterior(reverse(merge(cgrph.echnks(eneu))),{xtarg,ytarg}); 
-t1 = toc(start);
-in = in1 & ~in2;
-out = ~in;
-
 ids= chunkgraphinregion(cgrph,{xtarg,ytarg});
-nnz(in-(ids==2))
+in = ids == 2;
+in2 = ids == 3;
 
 fprintf('%5.2e s : time to find points in domain\n',t1)
 

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)
@@ -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)
-% profile viewer
+assert(err < 1e-10)
+% profile viewer