From 05e71bd81c13de4c655e8352d81b212fb8c09915 Mon Sep 17 00:00:00 2001 From: jghoskins <32715900+jghoskins@users.noreply.github.com> Date: Mon, 8 Sep 2025 19:46:00 -0500 Subject: [PATCH 1/2] Create hkdiffgreen_ff.m farfield --- chunkie/+chnk/+flex2d/hkdiffgreen_ff.m | 192 +++++++++++++++++++++++++ 1 file changed, 192 insertions(+) create mode 100644 chunkie/+chnk/+flex2d/hkdiffgreen_ff.m diff --git a/chunkie/+chnk/+flex2d/hkdiffgreen_ff.m b/chunkie/+chnk/+flex2d/hkdiffgreen_ff.m new file mode 100644 index 0000000..b1d489f --- /dev/null +++ b/chunkie/+chnk/+flex2d/hkdiffgreen_ff.m @@ -0,0 +1,192 @@ +function [val,grad,hess,der3,der4,der5] = hkdiffgreen_ff(k,src,targ,opts) + +prefac = 1i/(2*sqrt(2*pi))*exp(-1i*pi/4); + +[~,ns] = size(src); +[~,nt] = size(targ); + +xs = repmat(src(1,:),nt,1); +ys = repmat(src(2,:),nt,1); + +xt = repmat(targ(1,:).',1,ns); +yt = repmat(targ(2,:).',1,ns); +rt = sqrt(xt.^2+yt.^2); + +xt = xt./rt; +yt = yt./rt; + +dp = xt.*xs+yt.*ys; + +g0 = prefac*exp(-1i*dp*k); + + +% evaluate potential and derivatives + +if nargout > 0 + val = g0; +end +if nargout > 1 + grad(:,:,1) = -1i*k*g0.*xt; + grad(:,:,2) = -1i*k*g0.*yt; + grad = -grad; +end +if nargout > 2 + hess(:,:,1) = (-1i*k)^2*g0.*xt.^2; + hess(:,:,2) = (-1i*k)^2*g0.*xt.*yt; + hess(:,:,3) = (-1i*k)^2*g0.*yt.^2; +end +if nargout > 3 + der3(:,:,1) = (-1i*k)^3*g0.*(xt.^3); + der3(:,:,2) = (-1i*k)^3*g0.*(xt.^2).*yt; + der3(:,:,3) = (-1i*k)^3*g0.*(yt.^2).*xt; + der3(:,:,4) = (-1i*k)^3*g0.*(yt.^3); + der3 = -der3; +end + +if nargout > 4 + der4(:,:,1) = (-1i*k)^4*g0.*(xt.^4); + der4(:,:,2) = (-1i*k)^4*g0.*(xt.^3).*(yt); + der4(:,:,3) = (-1i*k)^4*g0.*(xt.^2).*(yt.^2); + der4(:,:,4) = (-1i*k)^4*g0.*(xt).*(yt.^3); + der4(:,:,5) = (-1i*k)^4*g0.*(yt.^4); +end + +if nargout > 5 + der5(:,:,1) = (-1i*k)^5*g0.*(xt.^5); + der5(:,:,2) = (-1i*k)^5*g0.*(xt.^4).*(yt); + der5(:,:,3) = (-1i*k)^5*g0.*(xt.^3).*(yt.^2); + der5(:,:,4) = (-1i*k)^5*g0.*(xt.^2).*(yt.^3); + der5(:,:,5) = (-1i*k)^5*g0.*(xt.^1).*(yt.^4); + der5(:,:,6) = (-1i*k)^5*g0.*(yt.^5); + der5 = -der5; +end + +end + +function [g0,g1,g21,g321,g4321,g54321] = diff_h0k0_and_rders(k,r,r2logrfac) +% g0 = g +% g1 = g' +% g21 = g'' - g'/r +% g321 = g''' - 3*g''/r + 3g'/r^2 = g''' - 3*g21/r +% g4321 = g'''' - 6*g'''/r + 15*g''/r^2 - 15*g'/r^3 +% = g''''- 6 g321/r - 3 g21/r^2 +% g54321 = g''''' - 10g''''/r + 45g'''/r^2 -105g''/r^3 + 105g'/r^4 +% = g5 - 10g4321/r - 15 g321/r^2 + +g0 = zeros(size(r)); +g1 = zeros(size(r)); +g321 = zeros(size(r)); +g4321 = zeros(size(r)); +g54321 = zeros(size(r)); +g21 = zeros(size(r)); + +io4 = 1i*0.25; +o2p = 1/(2*pi); + +isus = abs(k)*r < 1; +%isus = false(size(r)); +%isus = true(size(r)); + +% straightforward formulas for sufficiently large + +rnot = r(~isus); + +kr = k*rnot; + +h0 = besselh(0,1,kr); +h1 = besselh(1,1,kr); +h0i = besselh(0,1,1i*kr); +h1i = besselh(1,1,1i*kr); + +rm1 = 1./rnot; +rm2 = rm1.*rm1; +rm3 = rm1.*rm2; +rm4 = rm1.*rm3; + +r2fac = (1-r2logrfac)*k*k*0.5*o2p; +logr = log(rnot); +g0(~isus) = io4*(h0 - h0i) - r2fac*rnot.*rnot.*logr; +g1(~isus) = io4*(-k*h1 + 1i*k*h1i) - r2fac*(rnot+2*rnot.*logr); +g21(~isus) = io4*(-k*k*h0 + k*h1.*rm1 - k*k*h0i - 1i*k*h1i.*rm1) ... + - r2fac*(3+2*logr) - g1(~isus).*rm1; +g321(~isus) = io4*(k*k*h0.*rm1 + k*(k*k-2*rm2).*h1 ... + + k*k*h0i.*rm1 - 1i*k*(-k*k-2*rm2).*h1i) - r2fac*2*rm1 - 3*g21(~isus).*rm1; +g4321(~isus) = io4*(k*(3*rm2-k*k).*(2*h1.*rm1-k*h0) - ... + 1i*k*(3*rm2+k*k).*(2*h1i.*rm1-1i*k*h0i)) + r2fac*2*rm2 - 6*g321(~isus).*rm1 ... + - 3*g21(~isus).*rm2; +g54321(~isus) = io4*(k*( (12*k*rm3-2*k^3*rm1).*h0 + ... + (-24*rm4+7*k*k*rm2-k^4).*h1) - 1i*k*( (12*1i*k*rm3+1i*2*k^3*rm1).*h0i + ... + (-24*rm4-7*k*k*rm2-k^4).*h1i)) - r2fac*4*rm3 - 10*g4321(~isus).*rm1 - ... + 15*g321(~isus).*rm2; + +% manually cancel when small + +rsus = r(isus); +rm1 = 1./rsus; +rm2 = rm1.*rm1; +rm3 = rm1.*rm2; +rm4 = rm1.*rm3; +rm5 = rm1.*rm4; + +r2 = rsus.*rsus; +r4 = r2.*r2; + +nterms = 7; + +% relevant parts of hankel function represented as power series +[cf1,cf2] = chnk.flex2d.besselkikdiff_etc_pscoefs(nterms); +cf1(1) = cf1(1)*r2logrfac; +kpow = k*k*(k.^(4*(0:(nterms-1)))); % k^2, k^6, k^10, ... +cf1 = cf1(:).*kpow(:); cf2 = cf2(:).*kpow(:); + + +jdiff = chnk.flex2d.pseval(cf1,r4).*r2; +f = chnk.flex2d.pseval(cf2,r4).*r2; + +% differentiate power series to get derivatives +fac = 2+4*(0:(nterms-1)).'; +jdiffd1 = chnk.flex2d.pseval(cf1.*fac,r4).*rsus; +fd1 = chnk.flex2d.pseval(cf2.*fac,r4).*rsus; +d = fac.*(fac-1)-fac; +jdiffd21 = chnk.flex2d.pseval(cf1.*d,r4); +fd21 = chnk.flex2d.pseval(cf2.*d,r4); +% third deriv and higher the first term is dead (this is subtle but true +% for the 21, 321, 4321, etc derivative) +cf1 = cf1(2:end); cf2 = cf2(2:end); fac = fac(2:end); +d = fac.*(fac-1).*(fac-2) - 3*(fac.*(fac-1)-fac); +jdiffd321 = chnk.flex2d.pseval(cf1.*d,r4).*rsus.*r2; +fd321 = chnk.flex2d.pseval(cf2.*d,r4).*rsus.*r2; +% g4321 = g'''' - 6*g'''/r + 15*g''/r^2 - 15*g'/r^3 +d = fac.*(fac-1).*(fac-2).*(fac-3)-6*fac.*(fac-1).*(fac-2)+15*(fac.*(fac-1)-fac); +jdiffd4321 = chnk.flex2d.pseval(cf1.*d,r4).*r2; +fd4321 = chnk.flex2d.pseval(cf2.*d,r4).*r2; +% g54321 = g''''' - 10g''''/r + 45g'''/r^2 -105g''/r^3 + 105g'/r^4 +d = fac.*(fac-1).*(fac-2).*(fac-3).*(fac-4) - ... + 10*fac.*(fac-1).*(fac-2).*(fac-3) + 45*fac.*(fac-1).*(fac-2) - ... + 105*(fac.*(fac-1)-fac); +jdiffd54321 = chnk.flex2d.pseval(cf1.*d,r4).*rsus; +fd54321 = chnk.flex2d.pseval(cf2.*d,r4).*rsus; + +% combine to get derivative of i/4 H - K/2pi +gam = 0.57721566490153286060651209; +logr = log(rsus) + log(k) - log(2) + gam; +const1 = (1-r2logrfac)*o2p*(log(k)+gam-log(2))*k*k/2; + +j0 = besselj(0,k*rsus); +j1 = besselj(1,k*rsus); +j2 = besselj(2,k*rsus); +j3 = besselj(3,k*rsus); +j4 = besselj(4,k*rsus); +j5 = besselj(5,k*rsus); +g0(isus) = io4*j0 - o2p*(f + logr.*jdiff) + const1*r2; +g1(isus) = io4*(-k*j1) - o2p*(fd1 + logr.*jdiffd1 + jdiff.*rm1) + 2*const1*rsus; +g21(isus) = io4*(k*k*j2) - o2p*(fd21 + logr.*jdiffd21 + 2*jdiffd1.*rm1 - 2*jdiff.*rm2); +g321(isus) = io4*(-k^3*j3) - o2p*(fd321 + logr.*jdiffd321 + 3*jdiffd21.*rm1 - ... + 6*jdiffd1.*rm2 + 8*jdiff.*rm3); +g4321(isus) = io4*(k^4*j4)- o2p*(fd4321 + logr.*jdiffd4321 + 4*jdiffd321.*rm1 - ... + 12*jdiffd21.*rm2 + 32*jdiffd1.*rm3 - 48*jdiff.*rm4); +g54321(isus) = io4*(-k^5*j5) - o2p*(fd54321 + logr.*jdiffd54321 + ... + 5*jdiffd4321.*rm1 - 20*jdiffd321.*rm2 + 80*jdiffd21.*rm3 - 240*jdiffd1.*rm4 + 384*jdiff.*rm5); + +end + From 924115b6c4257b410c47b7e1033a0b5f2dd5faac Mon Sep 17 00:00:00 2001 From: peter-nekrasov Date: Mon, 8 Sep 2025 23:51:46 -0500 Subject: [PATCH 2/2] far field kernels and demos --- chunkie/+chnk/+flex2d/kern.m | 112 +++++++++ .../demo_clamped_plate_scatter_far_field.m | 183 +++++++++++++++ .../demo/demo_free_plate_scatter_far_field.m | 218 ++++++++++++++++++ .../demo_supported_plate_scatter_far_field.m | 215 +++++++++++++++++ 4 files changed, 728 insertions(+) create mode 100644 chunkie/demo/demo_clamped_plate_scatter_far_field.m create mode 100644 chunkie/demo/demo_free_plate_scatter_far_field.m create mode 100644 chunkie/demo/demo_supported_plate_scatter_far_field.m diff --git a/chunkie/+chnk/+flex2d/kern.m b/chunkie/+chnk/+flex2d/kern.m index 7eefe85..aa1f32c 100644 --- a/chunkie/+chnk/+flex2d/kern.m +++ b/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 diff --git a/chunkie/demo/demo_clamped_plate_scatter_far_field.m b/chunkie/demo/demo_clamped_plate_scatter_far_field.m new file mode 100644 index 0000000..4d3b611 --- /dev/null +++ b/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') \ No newline at end of file diff --git a/chunkie/demo/demo_free_plate_scatter_far_field.m b/chunkie/demo/demo_free_plate_scatter_far_field.m new file mode 100644 index 0000000..5d71d51 --- /dev/null +++ b/chunkie/demo/demo_free_plate_scatter_far_field.m @@ -0,0 +1,218 @@ +%DEMO_FREE_PLATE_SCATTER_FAR_FIELD +% +% Define an exterior scattering problem on a starfish-shaped domain and +% and evaluate the far field pattern +% + +% 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 + +% defining free plate kernels + +fkern1 = @(s,t) chnk.flex2d.kern(zk, s, t, 'free_plate', nu); % build the desired kernel +double = @(s,t) chnk.lap2d.kern(s,t,'d'); +hilbert = @(s,t) chnk.lap2d.kern(s,t,'hilb'); + +opts = []; +opts.sing = 'log'; + +opts2 = []; +opts2.sing = 'pv'; + +% building system matrix + +start = tic; +sysmat1 = chunkermat(chnkr,fkern1, opts); +D = chunkermat(chnkr, double, opts); +H = chunkermat(chnkr, hilbert, opts2); + +sysmat = zeros(2*chnkr.npt); +sysmat(1:2:end,1:2:end) = sysmat1(1:4:end,1:2:end) + sysmat1(3:4:end,1:2:end)*H - 2*((1+nu)/2)^2*D*D; +sysmat(2:2:end,1:2:end) = sysmat1(2:4:end,1:2:end) + sysmat1(4:4:end,1:2:end)*H; +sysmat(1:2:end,2:2:end) = sysmat1(1:4:end,2:2:end) + sysmat1(3:4:end,2:2:end); +sysmat(2:2:end,2:2:end) = sysmat1(2:4:end,2:2:end) + sysmat1(4:4:end,2:2:end); + +D = [-1/2 + (1/8)*(1+nu).^2, 0; 0, 1/2]; % jump matrix +D = kron(eye(chnkr.npt), D); + +sys = D + sysmat; +t1 = toc(start); +fprintf('%5.2e s : time to assemble matrix\n',t1) + +% building RHS + +[~, ~, hess, third] = planewave(kvec, chnkr.r); + +nx = chnkr.n(1,:).'; +ny = chnkr.n(2,:).'; + +dx = chnkr.d(1,:).'; +dy = chnkr.d(2,:).'; + +ds = sqrt(dx.*dx+dy.*dy); +taux = (dx./ds); % normalization +tauy = (dy./ds); + +kappa = signed_curvature(chnkr); +kappa = kappa(:)'; + +firstbc = -((hess(:,1).*(nx.*nx) + hess(:,2).*(2*nx.*ny) + hess(:,3).*(ny.*ny)) + ... + nu.*(hess(:,1).*(taux.*taux) + hess(:,2).*(2*taux.*tauy) + hess(:,3).*(tauy.*tauy))); + +secondbc = -((third(:,1).*(nx.*nx.*nx) + third(:,2).*(3*nx.*nx.*ny) +... + third(:,3).*(3*nx.*ny.*ny) + third(:,4).*(ny.*ny.*ny)) + ... + (2-nu).*(third(:,1).*(taux.*taux.*nx) + third(:,2).*(taux.*taux.*ny + 2*taux.*tauy.*nx) +... + third(:,3).*(2*taux.*tauy.*ny+ tauy.*tauy.*nx) +... + + third(:,4).*(tauy.*tauy.*ny)) + ... + (1-nu).*kappa'.*(hess(:,1).*taux.*taux + hess(:,2).*(2*taux.*tauy) + hess(:,3).*tauy.*tauy+... + -(hess(:,1).*nx.*nx + hess(:,2).*(2*nx.*ny) + hess(:,3).*ny.*ny))); + +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,500); 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, 'free_plate_eval', nu); + +dens_comb = zeros(3*chnkr.npt,1); +dens_comb(1:3:end) = sol(1:2:end); +dens_comb(2:3:end) = H*sol(1:2:end); +dens_comb(3:3:end) = sol(2:2:end); + +start1 = tic; +uscat = chunkerkerneval(chnkr, ikern,dens_comb,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,"Free 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, 'free_plate_eval_ff', nu); + +start1 = tic; +uscat = chunkerkerneval(chnkr, ikern_ff,dens_comb,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,'Free 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') \ No newline at end of file diff --git a/chunkie/demo/demo_supported_plate_scatter_far_field.m b/chunkie/demo/demo_supported_plate_scatter_far_field.m new file mode 100644 index 0000000..eefbb4d --- /dev/null +++ b/chunkie/demo/demo_supported_plate_scatter_far_field.m @@ -0,0 +1,215 @@ +%DEMO_SUPPORTED_PLATE_SCATTER_FAR_FIELD +% +% Define an exterior scattering problem on a starfish-shaped domain and +% evaluate the far field pattern +% +% this demonstration requires advanced usage in which additional geometric +% info is passed to the kernels using chunker data. this demo also avoids +% precision loss in the kernel evaluators by using the "native" quadrature +% on certain components of the matrix + +% 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 = 5; +amp = 0.25; +start = tic; chnkr = chunkerfunc(@(t) starfish(t,narms,amp),cparams,pref); +chnkr = makedatarows(chnkr,2); +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 + +% calculating curvature info + +kappa = signed_curvature(chnkr); +kp = arclengthder(chnkr,kappa); +kpp = arclengthder(chnkr,kp); + +% supported plate kernels expect (d/ds) kappa in the first data row +% and (d^2/ds^2) kappa in the second data row + +chnkr.data(1,:,:) = kp; +chnkr.data(2,:,:) = kpp; + +% defining supported plate kernels + +fkern1 = @(s,t) chnk.flex2d.kern(zk, s, t, 'supported_plate_log',nu); % build the desired kernel +fkern2 = @(s,t) chnk.flex2d.kern(zk, s, t, 'supported_plate_smooth',nu); % build the desired kernel + +opts = []; +opts.sing = 'log'; + +opts2 = []; +opts2.quad = 'native'; +opts2.sing = 'smooth'; + +% building system matrix + +start = tic; +M = chunkermat(chnkr,fkern1, opts); +M2 = chunkermat(chnkr,fkern2, opts2); + +c0 = (nu - 1)*(nu + 3)*(2*nu - 1)/(2*(3 - nu)); + +M(2:2:end,1:2:end) = M(2:2:end,1:2:end) + M2 + c0.*kappa(:).^2.*eye(chnkr.npt); +M = M - 0.5*eye(2*chnkr.npt); + +sys = M; +t1 = toc(start); +fprintf('%5.2e s : time to assemble matrix\n',t1) + +% building RHS + +[val, ~, hess] = planewave(kvec, chnkr.r); + +nx = chnkr.n(1,:).'; +ny = chnkr.n(2,:).'; + +dx = chnkr.d(1,:).'; +dy = chnkr.d(2,:).'; + +ds = sqrt(dx.*dx+dy.*dy); +taux = (dx./ds); % normalization +tauy = (dy./ds); + +firstbc = -val ; + +secondbc = -((hess(:, 1).*(nx.*nx) + hess(:, 2).*(2*nx.*ny) + hess(:, 3).*(ny.*ny))+... + nu*(hess(:, 1).*(taux.*taux) + hess(:, 2).*(2*taux.*tauy) + hess(:, 3).*(tauy.*tauy))); + +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,500); 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, 'supported_plate_eval', nu); + +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,"Supported 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, 'supported_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,'Supported 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') \ No newline at end of file