Merge pull request #145 from peter-nekrasov/farfield_flex

far field kernels and demos

Author: jghoskins

chunkie/+chnk/+flex2d/kern.m

@@ -228,6 +228,37 @@
 
 end
 
+% clamped plate kernels for far field evaluation
+if strcmpi(type, 'clamped_plate_eval_ff')
+
+    submat = zeros(nt,2*ns);
+
+    srcnorm = srcinfo.n;
+    srctang = srcinfo.d;
+    nx = repmat(srcnorm(1,:),nt,1);
+    ny = repmat(srcnorm(2,:),nt,1);
+    dx = repmat(srctang(1,:),nt,1);
+    dy = repmat(srctang(2,:),nt,1);
+    ds = sqrt(dx.*dx+dy.*dy);
+
+    taux = dx./ds;
+    tauy = dy./ds;
+
+    [~, ~, hess, third] = chnk.flex2d.hkdiffgreen_ff(zk, src, targ);           % Hankel part
+
+    K1 = -(1/(2*zk^2).*(third(:, :, 1).*(nx.*nx.*nx) + third(:, :, 2).*(3*nx.*nx.*ny) +...
+       third(:, :, 3).*(3*nx.*ny.*ny) + third(:, :, 4).*(ny.*ny.*ny)) ) - ...
+       (3/(2*zk^2).*(third(:, :, 1).*(nx.*taux.*taux) + third(:, :, 2).*(2*nx.*taux.*tauy + ny.*taux.*taux) +...
+       third(:, :, 3).*(nx.*tauy.*tauy + 2*ny.*taux.*tauy) + third(:, :, 4).*(ny.*tauy.*tauy)));  % G_{ny ny ny} + 3G_{ny tauy tauy}
+
+    K2 =  -(1/(2*zk^2).*(hess(:, :, 1).*(nx.*nx) + hess(:, :, 2).*(2*nx.*ny) + hess(:, :, 3).*(ny.*ny)))+...
+          (1/(2*zk^2).*(hess(:, :, 1).*(taux.*taux) + hess(:, :, 2).*(2*taux.*tauy) + hess(:, :, 3).*(tauy.*tauy))); % -G_{ny ny}  + G_{tauy tauy}
+
+    submat(:,1:2:end) = K1;
+    submat(:,2:2:end) = K2;
+
+end
+
 
 %%% FREE PLATE KERNELS
 
@@ -425,6 +456,35 @@
 
 end
 
+% free plate kernels for far field evaluation
+if strcmpi(type, 'free_plate_eval_ff')           
+   srcnorm = srcinfo.n;
+   srctang = srcinfo.d;
+   nu = varargin{1};
+
+   [val,grad] = chnk.flex2d.hkdiffgreen_ff(zk,src,targ); 
+   nx = repmat(srcnorm(1,:),nt,1);
+   ny = repmat(srcnorm(2,:),nt,1);
+
+   dx = repmat(srctang(1,:),nt,1);
+   dy = repmat(srctang(2,:),nt,1);
+
+   ds = sqrt(dx.*dx+dy.*dy);
+
+   taux = dx./ds; 
+   tauy = dy./ds;
+  
+   K1 = (-1/(2*zk^2).*(grad(:, :, 1).*(nx) + grad(:, :, 2).*ny)); 
+   K1H = ((1 + nu)/2).*(-1/(2*zk^2).*(grad(:, :, 1).*(taux) + grad(:, :, 2).*tauy));                    % G_{tauy}
+   K2 = 1/(2*zk^2).*val;
+
+   submat = zeros(nt,3*ns);
+   submat(:,1:3:end) = K1;
+   submat(:,2:3:end) = K1H;
+   submat(:,3:3:end) = K2;
+
+end
+
 %%% SUPPORTED PLATE KERNELS
 
 % boundary conditions applied to a point source
@@ -727,6 +787,58 @@
 
 end
 
+% supported plate kernels for far field evaluation
+if strcmpi(type, 'supported_plate_eval_ff')
+    srcnorm = srcinfo.n;
+    srctang = srcinfo.d;
+    srcd2 = srcinfo.d2;
+    coefs = varargin{1};
+    nu = coefs(1);
+    
+    nx = repmat(srcnorm(1,:),nt,1);
+    ny = repmat(srcnorm(2,:),nt,1);
+    
+    dx = repmat(srctang(1,:),nt,1);
+    dy = repmat(srctang(2,:),nt,1);
+    
+    ds = sqrt(dx.*dx+dy.*dy);
+    
+    taux = dx./ds;
+    tauy = dy./ds;
+    
+    dx1 = repmat(srctang(1,:),nt,1);
+    dy1 = repmat(srctang(2,:),nt,1);
+    
+    d2x1 = repmat(srcd2(1,:),nt,1);
+    d2y1 = repmat(srcd2(2,:),nt,1);
+    
+    denom = sqrt(dx1.^2+dy1.^2).^3;
+    numer = dx1.*d2y1-d2x1.*dy1;
+    
+    kappa = numer./denom; 
+    
+    kp = repmat(srcinfo.data(1,:),nt,1);
+    
+    a1 = 2-nu;
+    a2 = (-1+nu)*(7+nu)/(3 - nu);
+    a3 = (1-nu)*(3+nu)/(1+nu);
+    
+    [~, grad, hess, third] = chnk.flex2d.hkdiffgreen_ff(zk, src, targ, false); 
+    
+    K1 = -1/(2*zk^2)*(third(:,:,1).*nx.^3 + 3*third(:,:,2).*nx.^2.*ny + 3*third(:,:,3).*nx.*ny.^2 + third(:,:,4).*ny.^3) + ...
+         -a1/(2*zk^2)*(third(:,:,1).*nx.*taux.^2 + third(:,:,2).*(ny.*taux.^2 + 2*nx.*taux.*tauy) + third(:,:,3).*(nx.*tauy.^2 + 2*ny.*taux.*tauy) + third(:,:,4).*ny.*tauy.^2) + ...
+         a2*kappa./(2*zk^2).*(hess(:,:,1).*nx.^2 + 2*hess(:,:,2).*nx.*ny + hess(:,:,3).*ny.^2) + ...
+         -a3*kp./(2*zk^2).*(grad(:,:,1).*taux + grad(:,:,2).*tauy);
+
+    K2 = -1/(2*zk^2).*(grad(:,:,1).*nx + grad(:,:,2).*ny);
+
+    submat = zeros(nt,2*ns);
+
+    submat(:,1:2:end) = K1;
+    submat(:,2:2:end) = K2;
+
+end
+
 end
 
 

chunkie/demo/demo_clamped_plate_scatter_far_field.m

@@ -0,0 +1,183 @@
+%DEMO_CLAMPED_PLATE_SCATTER_FAR_FIELD
+%
+% Define an exterior scattering problem on a starfish-shaped domain and 
+% solve
+%
+
+% planewave vec
+
+kvec = [8;0];
+zk = norm(kvec);
+nu = 1/3;
+
+
+% discretize domain
+
+cparams = [];
+cparams.eps = 1.0e-6;
+cparams.nover = 0;
+cparams.maxchunklen = 4.0/zk; % setting a chunk length helps when the
+                              % frequency is known
+                              
+pref = []; 
+pref.k = 16;
+narms = 3;
+amp = 0.25;
+start = tic; chnkr = chunkerfunc(@(t) starfish(t,narms,amp),cparams,pref); 
+t1 = toc(start);
+
+fprintf('%5.2e s : time to build geo\n',t1)
+
+% plot geometry and data
+
+figure(1)
+clf
+plot(chnkr,'-x')
+hold on
+quiver(chnkr)
+axis equal
+
+% assembling system matrix
+
+fkern =  @(s,t) chnk.flex2d.kern(zk, s, t, 'clamped_plate');
+
+kappa = signed_curvature(chnkr);
+kappa = kappa(:);
+
+opts = [];
+opts.sing = 'log';
+
+start = tic;
+sys = chunkermat(chnkr,fkern, opts);
+sys = sys - 0.5*eye(2*chnkr.npt);
+sys(2:2:end,1:2:end) = sys(2:2:end,1:2:end) + kappa.*eye(chnkr.npt);
+
+t1 = toc(start);
+fprintf('%5.2e s : time to assemble matrix\n',t1)
+
+% building RHS
+
+[r1, grad] = planewave(kvec, chnkr.r);
+
+nx = chnkr.n(1,:); 
+ny = chnkr.n(2,:);
+
+normalderiv = grad(:, 1).*(nx.')+ grad(:, 2).*(ny.');                                % Dirichlet and Neumann BC(Clamped BC)                         
+
+firstbc = -r1;
+secondbc = -normalderiv;
+
+rhs = zeros(2*chnkr.npt, 1); rhs(1:2:end) = firstbc ; rhs(2:2:end) = secondbc;
+
+% Solving linear system
+
+start = tic; sol = gmres(sys,rhs,[],1e-12,100); t1 = toc(start);
+fprintf('%5.2e s : time for dense gmres\n',t1)    
+
+% evaluate at targets and plot
+
+rmin = min(chnkr); rmax = max(chnkr);
+xl = rmax(1)-rmin(1);
+yl = rmax(2)-rmin(2);
+nplot = 200;
+xtarg = linspace(rmin(1)-xl/2,rmax(1)+xl/2,nplot); 
+ytarg = linspace(rmin(2)-yl/2,rmax(2)+yl/2,nplot);
+[xxtarg,yytarg] = meshgrid(xtarg,ytarg);
+targets = zeros(2,length(xxtarg(:)));
+targets(1,:) = xxtarg(:); targets(2,:) = yytarg(:);
+
+start = tic; in = chunkerinterior(chnkr,{xtarg,ytarg}); t1 = toc(start);
+out = ~in;
+
+fprintf('%5.2e s : time to find points in domain\n',t1)
+
+ikern = @(s,t) chnk.flex2d.kern(zk, s, t, 'clamped_plate_eval');                              % build the kernel of evaluation          
+
+start1 = tic;
+uscat = chunkerkerneval(chnkr, ikern,sol, targets(:,out));
+t2 = toc(start1);
+fprintf('%5.2e s : time for kernel eval (for plotting)\n',t2)
+
+uin = planewave(kvec,targets(:,out));
+utot = uscat(:)+uin(:);
+
+maxu = max(abs(uin(:)));
+
+figure(2)
+clf
+
+t = tiledlayout(1,3,'TileSpacing','compact');
+
+nexttile
+zztarg = nan(size(xxtarg));
+zztarg(out) = uin;
+h=pcolor(xxtarg,yytarg,imag(zztarg),"FaceColor","interp");
+set(h,'EdgeColor','none')
+clim([-maxu,maxu])
+colormap(redblue);
+hold on
+plot(chnkr,'k','LineWidth',2)
+axis equal tight
+set(gca, "box","off","Xtick",[],"Ytick",[]);
+title('$u^{\textrm{inc}}$','Interpreter','latex','FontSize',12)
+
+nexttile
+zztarg = nan(size(xxtarg));
+zztarg(out) = uscat;
+h=pcolor(xxtarg,yytarg,imag(zztarg),"FaceColor","interp");
+set(h,'EdgeColor','none')
+clim([-maxu,maxu])
+colormap(redblue);
+hold on
+plot(chnkr,'k','LineWidth',2)
+axis equal tight
+set(gca, "box","off","Xtick",[],"Ytick",[]);
+title('$u^{\textrm{scat}}$','Interpreter','latex','FontSize',12)
+
+nexttile
+zztarg = nan(size(xxtarg));
+zztarg(out) = utot;
+h=pcolor(xxtarg,yytarg,imag(zztarg),"FaceColor","interp");
+set(h,'EdgeColor','none')
+clim([-maxu,maxu])
+colormap(redblue);
+hold on
+plot(chnkr,'k','LineWidth',2)
+axis equal tight
+set(gca, "box","off","Xtick",[],"Ytick",[]);
+colorbar
+title('$u^{\textrm{tot}}$','Interpreter','latex','FontSize',12)
+title(t,"Clamped Plate BCs")
+
+
+% evaluate the far field representation
+
+ts = linspace(-pi,pi,300);
+ts = ts(1:end-1);
+
+targets_ff = [cos(ts); sin(ts)];
+
+ikern_ff = @(s,t) chnk.flex2d.kern(zk, s, t, 'clamped_plate_eval_ff', nu); 
+
+start1 = tic;
+uscat = chunkerkerneval(chnkr, ikern_ff,sol,targets_ff);
+t2 = toc(start1);
+fprintf('%5.2e s : time for kernel eval (for plotting)\n',t2)
+
+figure(3)
+clf 
+
+t = tiledlayout(1,3,'TileSpacing','compact');
+title(t,'Clamped Plate Far Field Pattern')
+
+nexttile
+plot(ts,real(uscat));
+title('Real')
+
+nexttile
+plot(ts,imag(uscat));
+title('Imag')
+
+nexttile
+plot(ts,abs(uscat));
+title('Abs')