Merge pull request #143 from fastalgorithms/quasi_periodic_shapes

Update shape of vector valued quasiperiodic kernels

Author: askhamwhat

chunkie/+chnk/+helm2d/kern.m

@@ -40,7 +40,7 @@
 %                       [coef(1,1)*D coef(1,2)*S; coef(2,1)*D' coef(2,2)*S']
 %                type == 'c2trans' returns the combined field, and the 
 %                          normal derivative of the combined field
-%                        [coef(1)*D + coef(2)*S; coef(1)*D' + coef(2)*S']
+%                        [coef(1,1)*D + coef(1,2)*S; coef(2,1)*D' + coef(2,2)*S']
 %                type == 'trans_rep' returns the potential corresponding
 %                           to the transmission representation
 %                        [coef(1)*D coef(2)*S]
@@ -274,6 +274,7 @@
 case {'c2trans', 'c2t', 'c2tr'}
   coef = ones(2,1);
   if(nargin == 5); coef = varargin{1}; end
+  if numel(coef) == 2, coef = repmat(coef(:).',2,1); end
   targnorm = targinfo.n(:,:);
   srcnorm = srcinfo.n(:,:);
 
@@ -296,8 +297,8 @@
 
   submat = zeros(2*nt,ns);
 
-  submat(1:2:2*nt,:) = coef(1)*submatd + coef(2)*submats;
-  submat(2:2:2*nt,:) = coef(1)*submatdp + coef(2)*submatsp;
+  submat(1:2:2*nt,:) = coef(1,1)*submatd + coef(1,2)*submats;
+  submat(2:2:2*nt,:) = coef(2,1)*submatdp + coef(2,2)*submatsp;
 
 
 % Dirichlet and Neumann data corresponding to combined field (difference)

chunkie/+chnk/+helm2dquas/green.m

@@ -1,5 +1,5 @@
-function [val,grad,hess] = green(src,targ,zk,kappa,d,sn,l)
-%CHNK.HELM2DQUAS.GREEN evaluate the quasiperiodic helmholtz green's function
+function [val,grad,hess] = green(src,targ,zk,kappa,d,sn,l,ising)
+%CHNK.HELM2DQUAS.GREEN evaluate the quasiperiodic Helmholtz Green's function
 % for the given sources and targets
 %
 % Input:
@@ -8,12 +8,14 @@
 %   d - period
 %   sn - precomputed lattice sum integrals 
 %       (see chnk.helm2dquas.latticecoefs)
+%   l - number of periodic copies computed explicitly
+%   ising - if set to 0, only include the periodic copies. If set to 1,
+%       include the free-space part
 %
 % see also CHNK.HELM2DQUAS.KERN
 [~,nsrc] = size(src);
 [~,ntarg] = size(targ);
 
-
 xs = repmat(src(1,:),ntarg,1);
 ys = repmat(src(2,:),ntarg,1);
 
@@ -23,6 +25,11 @@
 rx = xt-xs;
 ry = yt-ys;
 
+
+rx = rx(:);
+ry = ry(:);
+
+
 nx = fix(rx/d);
 rx = rx - nx*d;
 
@@ -34,14 +41,9 @@
 
 r = sqrt(r2);
 
-rx = rx(:);
-ry = ry(:);
-r = r(:);
-
 npt = size(r,1);
 
-ythresh = d/2;
-% ythresh = d/10;
+ythresh = 2*d/2;
 iclose = abs(ry) < ythresh;
 ifar = ~iclose;
 
@@ -52,6 +54,9 @@
 ryclose = ry(iclose);
 rclose = r(iclose);
 
+nxclose = nx(iclose);
+nptclose = size(rxclose, 1);
+
 nkappa = length(kappa);
 
 val = zeros(nkappa,npt,1);
@@ -76,55 +81,57 @@
 % beta = sqrt((xi_m.^2-zk^2));
 beta = sqrt(1i*(xi_m-zk)).*sqrt(-1i*(xi_m+zk));
 
-% fhat = exp(-beta.*sqrt(ryfar.^2))./beta.*exp(1i*xi_m.*rxfar)/2;
 fhat = exp(-beta.*sqrt(ryfar.^2) + 1i*xi_m.*rxfar)./(2*beta);
 val(:,ifar,:) = sum(fhat,3)/(d);
 if nargout > 1
 grad(:,ifar,1) = sum(1i*xi_m.*fhat,3)/d;
 grad(:,ifar,2) = sum(-beta.*(sqrt(ryfar.^2)./ryfar).*fhat,3)/d;
-% grad(:,ifar,:) =[gx,gy];
 end
 
 if nargout >2
 hess(:,ifar,1) = sum(-xi_m.^2.*fhat,3)/d;
 hess(:,ifar,2) = sum(-1i*xi_m.*beta.*(sqrt(ryfar.^2)./ryfar).*fhat,3)/d;
 hess(:,ifar,3) = sum((beta.*(sqrt(ryfar.^2)./ryfar)).^2.*fhat,3)/d;
-% hess(:,ifar,:) = [hxx,hxy,hyy];
 end
 
-
 end
 
 alpha = (exp(1i*kappa(:)*d));
 
-val_near=0;
-grad_near = zeros(1,1,2);
-hess_near = zeros(1,1,3);
+val_near= zeros(nkappa,nptclose);
+grad_near = zeros(nkappa,nptclose,2);
+hess_near = zeros(nkappa,nptclose,3);
+ls = -l:l;
 if ~isempty(rxclose)
-    for i = -l:l
+    for i = ls
+        if ising == 1
+            iuse = true(nptclose,1);
+        else
+            iuse = nxclose ~= -i;
+        end
+
         rxi = rxclose - i*d;
         if nargout>2
         [vali,gradi,hessi] = chnk.helm2d.green(zk,[0;0],[rxi.';ryclose.']);
         vali = reshape(vali,1,[],1);
         gradi = reshape(gradi,1,[],2);
         hessi = reshape(hessi,1,[],3);
-        val_near = val_near + vali.*alpha.^i;
-        grad_near = grad_near + gradi.*alpha.^i;
-        hess_near = hess_near + hessi.*alpha.^i;
+        val_near(:,iuse) = val_near(:,iuse) + vali(:,iuse).*alpha.^i;
+        grad_near(:,iuse,:) = grad_near(:,iuse,:) + gradi(:,iuse,:).*alpha.^i;
+        hess_near(:,iuse,:) = hess_near(:,iuse,:) + hessi(:,iuse,:).*alpha.^i;
         elseif nargout > 1
         [vali,gradi] = chnk.helm2d.green(zk,[0;0],[rxi.';ryclose.']);
         vali = reshape(vali,1,[],1);
         gradi = reshape(gradi,1,[],2);
-        val_near = val_near + vali.*alpha.^i;
-        grad_near = grad_near + gradi.*alpha.^i;
+        val_near(:,iuse) = val_near(:,iuse) + vali(:,iuse).*alpha.^i;
+        grad_near(:,iuse,:) = grad_near(:,iuse,:) + gradi(:,iuse,:).*alpha.^i;
         else
         vali = chnk.helm2d.green(zk,[0;0],[rxi.';ryclose.']);
         vali = reshape(vali,1,[],1);
-        val_near = val_near + vali.*alpha.^i;
+        val_near(:,iuse) = val_near(:,iuse) + vali(:,iuse).*alpha.^i;
         end
     end
-    
-    % sn = sn.';
+
     N = size(sn,2)-1;
     ns = (0:N);
     ns_use = (0:N+2);
@@ -139,25 +146,26 @@
         Js(:,i) = besselj(ns_use(i),zk*rclose);
         end
     end
+    % t1 = tic;
     eip = (rxclose+1i*ryclose)./rclose;
     eipn = reshape(eip.^ns,1,[], N+1);
     cs = (eipn+1./eipn)/2;
 
     Js = reshape(Js,1,[],N+3);
     sn = reshape(sn,nkappa, 1, N+1);
 
-    val_far = 0.25*1i*Js(:,:,1).*sn(:,:,1) + 0.5*1i*sum(sn(:,:,2:end).*Js(:,:,2:end-2).*cs(:,:,2:end),3);
+    tmp = reshape(Js(:,:,2:end-2).*cs(:,:,2:end),[],N);
+    val_far = 0.25*1i*Js(:,:,1).*sn(:,:,1) + 0.5*1i*sn(:,2:end)*tmp.';
     val(:,iclose) = val_near+val_far;
 
-    
-    
-    
     if nargout >1
         DJs = cat(3,-Js(:,:,2),.5*(Js(:,:,1:end-3)-Js(:,:,3:end-1)))*zk;
         ss = (eipn-1./eipn)/2i;
 
-        grad_far_p = 0.25*1i*DJs(:,:,1).*sn(:,:,1) + 0.5*1i*sum(sn(:,:,2:end).*DJs(:,:,2:end).*cs(:,:,2:end),3);
-        grad_far_t = (0.5*1i*sum(-reshape((1:N),1,1,[]).*sn(:,:,2:end).*Js(:,:,2:end-2).*ss(:,:,2:end),3))./rclose.';
+        tmp = reshape(DJs(:,:,2:end).*cs(:,:,2:end),[],N);
+        grad_far_p = 0.25*1i*DJs(:,:,1).*sn(:,:,1) + 0.5*1i*sn(:,2:end)*tmp.';
+        tmp = reshape(Js(:,:,2:end-2).*ss(:,:,2:end),[],N)./rclose;
+        grad_far_t = (0.5*1i*((-reshape((1:N),1,[]).*sn(:,2:end))*tmp.'));
 
         grad_far = cat(3,cs(:,:,2).*grad_far_p - ss(:,:,2).*grad_far_t, ss(:,:,2).*grad_far_p + cs(:,:,2).*grad_far_t);
         grad(:,iclose,:) = grad_near + grad_far; 
@@ -172,37 +180,85 @@
         tmp_n = rclose.^(-4).*(-ns.*ryclose.*Js(:,:,1:end-2).*(ns.*ryclose.*cs+2*rxclose.*ss)+ ...
             rclose.*ryclose.*(ryclose.*cs + 2*ns.*rxclose.*ss).*DJs+ ...
             rclose.^2.*rxclose.^2.*cs.*DDJs);
-        hess_far_xx = 0.25*1i*tmp_n(:,:,1).*sn(:,:,1)+.5*1i*sum(sn(:,:,2:end).*tmp_n(:,:,2:end),3);
+        tmp_n = reshape(tmp_n,[],N+1);
+        hess_far_xx = 0.25*1i*tmp_n(:,1).'.*sn(:,1)+.5*1i*sn(:,2:end)*tmp_n(:,2:end).';
 
         tmp_n = ns.*Js(:,:,1:end-2).*(ns.*rxclose.*ryclose.*cs+(rxclose.^2-ryclose.^2).*ss).*rclose.^(-4) ...
             +rclose.^(-3).*(-(rxclose.*ryclose.*cs+ns.*(rxclose.^2-ryclose.^2).*ss).*DJs+...
             rclose.*rxclose.*ryclose.*cs.*DDJs);
 
-
-        hess_far_xy = 0.25*1i*tmp_n(:,:,1).*sn(:,:,1)+.5*1i*sum(sn(:,:,2:end).*tmp_n(:,:,2:end),3);
+        tmp_n = reshape(tmp_n,[],N+1);
+        hess_far_xy = 0.25*1i*tmp_n(:,1).'.*sn(:,1)+.5*1i*sn(:,2:end)*tmp_n(:,2:end).';
 
         tmp_n = rclose.^(-4).*(-ns.*rxclose.*Js(:,:,1:end-2).*(ns.*rxclose.*cs-2*ryclose.*ss)+ ...
             rclose.*rxclose.*(rxclose.*cs - 2*ns.*ryclose.*ss).*DJs+ ...
             rclose.^2.*ryclose.^2.*cs.*DDJs);
-        hess_far_yy = 0.25*1i*tmp_n(:,:,1).*sn(:,:,1)+.5*1i*sum(sn(:,:,2:end).*tmp_n(:,:,2:end),3);
+        tmp_n = reshape(tmp_n,[],N+1);
+        hess_far_yy = 0.25*1i*tmp_n(:,1).'.*sn(:,1)+.5*1i*sn(:,2:end)*tmp_n(:,2:end).';
 
         hess_far = cat(3,hess_far_xx, hess_far_xy, hess_far_yy);
 
-
         hess(:,iclose,:) = hess_near + hess_far;
     end
 end
 
 quasi_phase = exp(1i*kappa(:)*nx(:).'*d);
 
-val = reshape(quasi_phase.*val,nkappa*ntarg,nsrc);
-if nargout>1
-grad = reshape(quasi_phase.*grad,nkappa*ntarg,nsrc,2);
-end
-if nargout>2
-hess = reshape(quasi_phase.*hess,nkappa*ntarg,nsrc,3);
-end
-% end
 
-end
+if nargout == 1
+    val = quasi_phase.*val;
+
+    if ising == 0
+        isub = (abs(nx(:)) > max(ls)) | ifar;
+
+        if any(isub)
+        vali = chnk.helm2d.green(zk,[0;0],[rx(isub).'+ nx(isub).'*d;ry(isub).']);
+        vali = reshape(vali,1,[],1);
+        val(:,isub,:) = val(:,isub,:) - vali;
+        end
+    end
 
+    val = reshape(val,nkappa*ntarg,nsrc);
+elseif nargout == 2
+    val = quasi_phase.*val;
+    grad = quasi_phase.*grad;
+    
+    if ising == 0
+        isub = (abs(nx(:)) > max(ls)) | ifar;
+    
+        if any(isub)
+        [vali, gradi] = chnk.helm2d.green(zk,[0;0],[rx(isub).' + nx(isub).'*d;ry(isub).']);
+        vali = reshape(vali,1,[],1);
+        gradi = reshape(gradi,1,[],2);
+        val(:,isub,:) = val(:,isub,:) - vali;
+        grad(:,isub,:) = grad(:,isub,:) - gradi;
+        end
+    end
+
+    val = reshape(val,nkappa*ntarg,nsrc);
+    grad = reshape(grad,nkappa*ntarg,nsrc,2);
+elseif nargout == 3
+    val = quasi_phase.*val;
+    grad = quasi_phase.*grad;
+    hess = quasi_phase.*hess;
+    
+    if ising == 0
+        isub = (abs(nx(:)) > max(ls)) | ifar;
+
+        if any(isub)
+        [vali, gradi, hessi] = chnk.helm2d.green(zk,[0;0],[rx(isub).' + nx(isub).'*d;ry(isub).']);
+        vali = reshape(vali,1,[],1);
+        gradi = reshape(gradi,1,[],2);
+        hessi = reshape(hessi,1,[],3);
+
+        val(:,isub,:) = val(:,isub,:) - vali;
+        grad(:,isub,:) = grad(:,isub,:) - gradi;
+        hess(:,isub,:) = hess(:,isub,:) - hessi;
+        end
+    end
+
+    val = reshape(val,nkappa*ntarg,nsrc);
+    grad = reshape(grad,nkappa*ntarg,nsrc,2);
+    hess = reshape(hess,nkappa*ntarg,nsrc,3);
+end
+end