Improved complexification

chunkie/+chnk/chunkerkerneval_smooth.m

@@ -62,7 +62,7 @@
 % these are automatically ignored by the fmm and 
 % in chunkermat corrections
 if ~(strcmpi(imethod,'fmm') && isempty(flag))
-    flagslf = chnk.flagself(targinfo.r, chnkr.r, 1e-14*diam);
+    flagslf = chnk.flagself(real(targinfo.r), real(chnkr.r), 1e-14*diam);
     if isempty(flagslf)
         selfzero = sparse(opdims(1)*size(targinfo.r(:,:),2), ...
             opdims(2)*chnkr.npt);
@@ -141,7 +141,7 @@
     wts = wts(:);
 
     if strcmpi(imethod,'flam')
-        xflam1 = chnkr.r(:,:);
+        xflam1 = real(chnkr.r(:,:));
         xflam1 = repmat(xflam1,opdims(2),1);
         xflam1 = reshape(xflam1,chnkr.dim,numel(xflam1)/chnkr.dim);
 

chunkie/@chunker/flagnear.m

@@ -68,8 +68,8 @@
 
             % get all pairwise distances
             itarg = inbor(inbor > npt)-npt;
-            srcs = reshape(chnkr.rstor(:,isrc),dim,1,length(isrc));
-            targs = reshape(pts(:,itarg),dim,length(itarg),1);
+            srcs = reshape(real(chnkr.rstor(:,isrc)),dim,1,length(isrc));
+            targs = reshape(real(pts(:,itarg)),dim,length(itarg),1);
             dists = reshape(sqrt(sum((bsxfun(@minus,targs,srcs)).^2,1)), ...
                 length(itarg),length(isrc));
 

chunkie/@chunkgraph/chunkgraph.m

@@ -304,7 +304,7 @@
                         end
 
                         vs = chunkends(chnkr,[1,chnkr.nch]);
-                        vs = vs(:, [1,4]);
+                        vs = real(vs(:, [1,4]));
 
                     else
                         msg = "CHUNKGRAPH:CONSTRUCTOR: invalid edge descriptor." + ...

chunkie/@chunkgraph/findunbounded.m

@@ -36,8 +36,8 @@
         for jj =1:numel(edges)
             echnk = cgrph.echnks(abs(edges(jj)));     
             [allrends,alltends] = chunkends(echnk);
-            ts1 = alltends(:,1,:);
-            ts2 = alltends(:,2,:);
+            ts1 = real(alltends(:,1,:));
+            ts2 = real(alltends(:,2,:));
             angs1 = atan2(ts1(2,:),ts1(1,:));
             angs2 = atan2(ts2(2,:),ts2(1,:));
             angdiffs = angs2-angs1;

chunkie/@chunkgraph/vertextract.m

@@ -56,7 +56,7 @@
 end
 
 % compute the corresponding angles and sort them
-ds = [deplus,deminus];
+ds = real([deplus,deminus]);
 angs = atan2(ds(2,:),ds(1,:));
 [angs,isrtededges] = sort(angs);
 

chunkie/chunkerkerneval.m

@@ -162,6 +162,14 @@
     error("CHUNKERKERNEVAL: input 4 is not a supported type");
 end
 
+if ~isreal(chnkobj.r) && (opts_use.accel || opts_use.forcefmm)
+    warning('fmm not supported for complex chunkers')
+end
+
+if ~isreal(targinfo.r) && (opts_use.accel || opts_use.forcefmm || opts_use.flam)
+    warning('accel not supported for complex targets')
+end
+
 if icgrph && mk > 1
     ids = chunkgraphinregion(chnkobj,targinfo.r);
 else
@@ -282,7 +290,9 @@
 
 if opts_use.forcepquad
     optsflag = []; optsflag.fac = opts_use.fac;
-    flag = flagnear(chnkr0,targinfo.r,optsflag);
+    chnkrflag = chunkerpoints(struct('r',real(chnkr0.r), ... 
+        'd',real(chnkr0.d),'d2',real(chnkr0.d2)));
+    flag = flagnear(chnkrflag,targinfo.r,optsflag);
     fints(irow0) = fints(irow0) + chnk.chunkerkerneval_smooth(chnkr0,kern0,opdims0,dens0,targinfo0, ...
         flag,opts_use);
 
@@ -301,7 +311,9 @@
 
 rho = 1.8;
 optsflag = [];  optsflag.rho = rho;
-flag = flagnear_rectangle(chnkr0,targinfo0.r,optsflag);
+chnkrflag = chunkerpoints(struct('r',real(chnkr0.r), ...
+    'd',real(chnkr0.d),'d2',real(chnkr0.d2)));
+flag = flagnear_rectangle(chnkrflag,targinfo0.r,optsflag);
 
 npoly = chnkr0.k*2;
 nlegnew = chnk.ellipse_oversample(rho,npoly,opts_use.eps);

chunkie/chunkerkernevalmat.m

@@ -212,7 +212,7 @@
 optsadap = []; 
 optsadap.eps = eps;
 
-if corrections
+if corrections || nonsmoothonly
     mat = sparse(nout,icollocs(end)-1);
 else
     mat = zeros(nout,icollocs(end)-1);
@@ -267,7 +267,9 @@
 % TODO: switch to oversampling version with bounding ellipse like chunkerkerneval
 
 optsflag = []; optsflag.fac = fac;
-flag = flagnear(chnkr0,targinfo0.r,optsflag);
+chnkrflag = chunkerpoints(struct('r',real(chnkr0.r), ... 
+    'd',real(chnkr0.d),'d2',real(chnkr0.d2)));
+flag = flagnear(chnkrflag,real(targinfo0.r),optsflag);
 
 if forcepquad
     spmat = chunkerkernevalmat_pquad(chnkr0,kern0,opdims0, ...

devtools/test/complexificationTest.m

@@ -16,8 +16,6 @@ function complexificationTest0()
 zk1 = zks(1);
 zk2 = zks(2);
 
-close all
-
 a = 3;
 b = 1/a/2;
 t0 =-5;
@@ -34,6 +32,11 @@ function complexificationTest0()
 cparams.ifclosed = 0;
 chnkr = chunkerfuncuni(f, nch, cparams);
 
+rends = chunkends(chnkr,[1,chnkr.nch]);
+verts = rends(:,[1,4]);
+fchnk{1} = chnkr;
+cgrph = chunkgraph(real(verts),[1;2],fchnk);
+
 
 ddiff = kernel('helmdiff', 'd', zks);
 sdiff = (-1)*kernel('helmdiff', 's', zks);
@@ -43,13 +46,13 @@ function complexificationTest0()
 K = kernel([ddiff, sdiff; ...
            dpdiff, spdiff]);
 
-n = 2*chnkr.npt;
-sysmat = chunkermat(chnkr, K);
+n = 2*cgrph.npt;
+sysmat = chunkermat(cgrph, K);
 sysmat = sysmat - eye(n);
 
 %% Analytic solution test
 
-r = chnkr.r;
+r = cgrph.r;
 r0 = [0.1; -3];
 s = [];
 s.r = r0;
@@ -58,27 +61,23 @@ function complexificationTest0()
 dk1 = kernel('helm', 'd', zk1);
 skp1 = kernel('helm', 'sp', zk1);
 rhs = complex(zeros(n,1));
-rhs(1:2:end) = sk1.eval(s, chnkr);
-rhs(2:2:end) = skp1.eval(s, chnkr);
+rhs(1:2:end) = sk1.eval(s, cgrph);
+rhs(2:2:end) = skp1.eval(s, cgrph);
 
 dens = sysmat \ rhs;
 
-rt = [0.3; 2];
+rt = [[0.3; 2] [-0.5;0.001]];
 targ = [];
 targ.r = rt;
 Keval = kernel([dk1, (-1)*sk1]);
 
 opts = [];
-opts.forcesmooth = true;
+% opts.forcesmooth = true;
 opts.accel = false;
-pot = chunkerkerneval(chnkr, Keval, dens, targ, opts);
+pot = chunkerkerneval(cgrph, Keval, dens, targ, opts);
 
 uex = sk1.eval(s, targ);
-assert(norm(uex - pot) < 1e-10); 
-
-
-
-
+assert(norm(uex - pot) < 1e-12); 
 
 end