chunkerflam and abs val

chunkie/+chnk/+axissymhelm2d/kern.m

@@ -40,7 +40,7 @@
 %                type == 'sprime', normal derivative of single
 %                      layer S'
 %                type == 'dprime', normal derivative of double layer D'
-%                type == 'ddiff', D_{k}' - D_{ik}', for this routine k must be real
+%                type == 'dprimediff', D_{k}' - D_{ik}', for this routine k must be real
 %                type == 'c', combined layer kernel coef(1) D + coef(2) S
 %
 %   varargin{1} - coef: length 2 array in the combined layer 

chunkie/+chnk/+flam/kernbyindex.m

@@ -37,8 +37,12 @@
 % l2scale - boolean type that determines if we should 
 %    rescale the matrix by l2scale. the default value is false. 
 
-if isa(kern,'kernel')
-    kern = kern.eval;
+if ~isa(kern,'kernel')
+    try 
+        kern = kernel(kern);
+    catch
+        error('KERNBYINDEX: fourth input kern not of supported type');
+    end
 end
 
 if nargin < 7
@@ -145,9 +149,9 @@
 
 
         if size(kern) == 1
-            matuni = kern(srcinfo,targinfo);
+            matuni = kern.eval(srcinfo,targinfo);
         else
-            matuni = kern{itrg,isrc}(srcinfo,targinfo);
+            matuni = kern(itrg,isrc).eval(srcinfo,targinfo);
         end
 
         % scale matrix by weights

chunkie/+chnk/+flam/proxyfun.m

@@ -41,13 +41,14 @@
 %    rescale the matrix by l2scale. the default value is false. 
 % 
 
-if isa(kern,'kernel')
-    kern = kern.eval;
+if ~isa(kern,'kernel')
+    try 
+        kern = kernel(kern);
+    catch
+        error('PROXYFUN: sixth input kern not of supported type');
+    end
 end
 
-if nargin < 13
-    l2scale = false;
-end
 
 % scaled proxy points and weights (no scaling necessary on tangents)
 
@@ -153,9 +154,9 @@
         opdims_trans = opdims_mat(:,j,i);
 
         if length(kern) == 1
-            matuni = kern(srcinfo,targinfo);
+            matuni = kern.eval(srcinfo,targinfo);
         else
-            matuni = kern{i,j}(srcinfo,targinfo);
+            matuni = kern(i,j).eval(srcinfo,targinfo);
         end
 
         wsrc = wts(juni);
@@ -180,9 +181,9 @@
 
         if ifaddtrans
             if length(kern) == 1
-                matuni2 = kern(targinfo,srcinfo);
+                matuni2 = kern.eval(targinfo,srcinfo);
             else
-                matuni2 = kern{j,i}(targinfo,srcinfo);
+                matuni2 = kern(j,i).eval(targinfo,srcinfo);
             end
 
             wsrc = pw(:);

chunkie/+chnk/chunkerkerneval_smooth.m

@@ -87,7 +87,7 @@
         % nothing to ignore
         for i = 1:nch
             densvals = dens(:,:,i); densvals = densvals(:);
-            dsdtdt = sqrt(sum(abs(chnkr.d(:,:,i)).^2,1));
+            dsdtdt = sqrt(sum(chnkr.d(:,:,i).^2,1));
             dsdtdt = dsdtdt(:).*w(:);
             dsdtdt = repmat( (dsdtdt(:)).',opdims(2),1);
             densvals = densvals.*(dsdtdt(:));
@@ -108,7 +108,7 @@
         % ignore interactions in flag array
         for i = 1:nch
             densvals = dens(:,:,i); densvals = densvals(:);
-            dsdtdt = sqrt(sum(abs(chnkr.d(:,:,i)).^2,1));
+            dsdtdt = sqrt(sum(chnkr.d(:,:,i).^2,1));
             dsdtdt = dsdtdt(:).*w(:);
             dsdtdt = repmat( (dsdtdt(:)).',opdims(2),1);
             densvals = densvals.*(dsdtdt(:));
@@ -183,7 +183,7 @@
     if ~isempty(flag)
         for i = 1:nch
             densvals = dens(:,:,i); densvals = densvals(:);
-            dsdtdt = sqrt(sum(abs(chnkr.d(:,:,i)).^2,1));
+            dsdtdt = sqrt(sum(chnkr.d(:,:,i).^2,1));
             dsdtdt = dsdtdt(:).*w(:);
             dsdtdt = repmat( (dsdtdt(:)).',opdims(2),1);
             densvals = densvals.*(dsdtdt(:));

chunkie/chunkerflam.m

@@ -143,8 +143,8 @@
     error(msg)
 end
 
-for chnkr = chnkrs
-    if or(chnkr.nch < 1,chnkr.k < 1)
+for i=1:length(chnkrs)
+    if chnkrs(i).nch < 1 || chnkrs(i).k < 1
         warning('empty chunker, doing nothing')
         return
     end

devtools/test/chunkgrphrcipTransmissionTest.m

@@ -48,12 +48,15 @@
 ncurve = size(edge2verts,1);
 
 fchnks    = cell(1,ncurve);
+cparams = cell(ncurve,1);
 for icurve = 1:ncurve
     fchnks{icurve} = @(t) circulararc(t,cpars{1});
     fchnks{icurve} = @(t) sinearc(t,amp,frq);
+    cparams{icurve}.ta = 0;
+    cparams{icurve}.tb = 1;
 end
 
-[cgrph] = chunkgraph(verts,edge2verts,fchnks);
+[cgrph] = chunkgraph(verts,edge2verts,fchnks,cparams);
 
 
 vstruc = procverts(cgrph);
@@ -63,7 +66,7 @@
 quiver(cgrph);
 
 nregions = 2;
-ks = [1.1;2.1]*30;
+ks = [1.1;2.1]*10;
 coefs = [1.0;1.0];
 cs(1,1:ncurve) = 1;
 cs(2,1:ncurve) = 2;

devtools/test/complexificationTest.m

@@ -0,0 +1,95 @@
+% test solve_sommerfeld_dens and eval_sommerfeld_dens
+% and compare to complexification solution
+
+%%%%
+%
+%       .   .   .   testing a point charge below
+%
+%%%%
+
+eps = 1E-16;
+zks = pi*[1,1.3];
+zk1 = zks(1);
+zk2 = zks(2);
+
+close all
+
+a = 3;
+b = 1/a/2;
+t0 =-5;
+t1 = 5;
+
+f = @(t) flat_interface(t, a, b, t0, t1);
+nch = 50;
+xrad = -log(eps)/abs(real(zk1)) + max([abs(t0),abs(t1)]);
+xmin = -xrad;
+xmax =  xrad;
+cparams = [];
+cparams.ta = xmin;
+cparams.tb = xmax;
+cparams.ifclosed = 0;
+chnkr = chunkerfuncuni(f, nch, cparams);
+
+
+ddiff = kernel('helmdiff', 'd', zks);
+sdiff = (-1)*kernel('helmdiff', 's', zks);
+dpdiff = kernel('helmdiff', 'dp', zks);
+spdiff = (-1)*kernel('helmdiff', 'sp', zks);
+
+K = kernel([ddiff, sdiff; ...
+           dpdiff, spdiff]);
+
+n = 2*chnkr.npt;
+sysmat = chunkermat(chnkr, K);
+sysmat = sysmat - eye(n);
+
+%% Analytic solution test
+
+r = chnkr.r;
+r0 = [0.1; -3];
+s = [];
+s.r = r0;
+
+sk1 = kernel('helm', 's', zk1);
+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);
+
+dens = sysmat \ rhs;
+
+rt = [0.3; 2];
+targ = [];
+targ.r = rt;
+Keval = kernel([dk1, (-1)*sk1]);
+
+opts = [];
+opts.forcesmooth = true;
+opts.accel = false;
+pot = chunkerkerneval(chnkr, Keval, dens, targ, opts);
+
+uex = sk1.eval(s, targ);
+assert(norm(uex - pot) < 1e-10); 
+
+
+
+function [f,fd,fdd] = flat_interface(t, a, b, t0, t1)
+
+    phi   = @(t,u,v,z) u*(t-z).*erfc(u*(t-z))*v - exp(-u^2*(t-z).^2)/sqrt(pi)*v;
+    phid  = @(t,u,v,z) u*erfc(u*(t-z))*v;
+    phidd = @(t,u,v,z) -u*u*exp(-u^2*(t-z).^2)*2*v/sqrt(pi);
+    f = zeros([2,size(t)]);
+    fd = zeros([2,size(t)]);
+    fdd = zeros([2,size(t)]);
+
+    f(1,:) = t + 1i*(phi(t,a,b,t0) - phi(t,-a,b,t1)); 
+    fd(1,:)= 1 + 1i*(phid(t,a,b,t0) - phid(t,-a,b,t1));
+    fdd(1,:) = 1i*(phidd(t,a,b,t0) - phidd(t,-a,b,t1));
+
+    f(2,:) = 0;
+    fd(2,:) = 0;
+    fdd(2,:) = 0;
+
+end
+