From 789698669d7e86b41ea05cefae3076f402b32318 Mon Sep 17 00:00:00 2001 From: jghoskins <32715900+jghoskins@users.noreply.github.com> Date: Thu, 23 May 2024 08:33:36 -0500 Subject: [PATCH] updates --- chunkie/+chnk/+helm2d/kern.m | 1368 +++++++++++++---- .../+chnk/+helm2d/smooth_part_derivatives.m | 82 + .../+chnk/+helm2d/thirdforth_derivatives.m | 178 +++ chunkie/+chnk/+lap2d/kern.m | 9 + chunkie/circle.m | 23 + .../demo/flexural/test_free_plate_circle.m | 201 +++ chunkie/ellipse.m | 16 + 7 files changed, 1582 insertions(+), 295 deletions(-) create mode 100644 chunkie/+chnk/+helm2d/smooth_part_derivatives.m create mode 100644 chunkie/+chnk/+helm2d/thirdforth_derivatives.m create mode 100644 chunkie/circle.m create mode 100644 chunkie/demo/flexural/test_free_plate_circle.m create mode 100644 chunkie/ellipse.m diff --git a/chunkie/+chnk/+helm2d/kern.m b/chunkie/+chnk/+helm2d/kern.m index d697efd..d9661ab 100644 --- a/chunkie/+chnk/+helm2d/kern.m +++ b/chunkie/+chnk/+helm2d/kern.m @@ -1,295 +1,1073 @@ -function submat= kern(zk,srcinfo,targinfo,type,varargin) -%CHNK.HELM2D.KERN standard Helmholtz layer potential kernels in 2D -% -% Syntax: submat = chnk.heml2d.kern(zk,srcinfo,targingo,type,varargin) -% -% Let x be targets and y be sources for these formulas, with -% n_x and n_y the corresponding unit normals at those points -% (if defined). Note that the normal information is obtained -% by taking the perpendicular to the provided tangential deriviative -% info and normalizing -% -% Kernels based on G(x,y) = i/4 H_0^{(1)}(zk |x-y|) -% -% D(x,y) = \nabla_{n_y} G(x,y) -% S(x,y) = G(x,y) -% S'(x,y) = \nabla_{n_x} G(x,y) -% D'(x,y) = \nabla_{n_x} \nabla_{n_y} G(x,y) -% -% Input: -% zk - complex number, Helmholtz wave number -% srcinfo - description of sources in ptinfo struct format, i.e. -% ptinfo.r - positions (2,:) array -% ptinfo.d - first derivative in underlying -% parameterization (2,:) -% ptinfo.d2 - second derivative in underlying -% parameterization (2,:) -% targinfo - description of targets in ptinfo struct format, -% if info not relevant (d/d2) it doesn't need to -% be provided. sprime requires tangent info in -% targinfo.d -% type - string, determines kernel type -% type == 'd', double layer kernel D -% type == 's', single layer kernel S -% type == 'sprime', normal derivative of single -% layer S' -% type == 'dprime', normal derivative of double layer D' -% type == 'c', combined layer kernel coef(1) D + coef(2) S -% type == 'stau', tangential derivative of single layer -% type == 'all', returns all four layer potentials, -% [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'] -% type == 'trans_rep' returns the potential corresponding -% to the transmission representation -% [coef(1)*D coef(2)*S] -% type == 'trans_rep_prime' returns the normal derivative -% corresponding to the transmission representation -% [coef(1)*D' coef(2)*S'] -% type == 'trans_rep_grad' returns the gradient corresponding -% to the transmission representation -% [coef(1)*d_x D coef(2)*d_x S; -% coef(1)*d_y D coef(2)*d_y S] -% -% varargin{1} - coef: length 2 array in the combined layer -% formula, 2x2 matrix for all kernels -% otherwise does nothing -% -% Output: -% submat - the evaluation of the selected kernel for the -% provided sources and targets. the number of -% rows equals the number of targets and the -% number of columns equals the number of sources -% -% see also CHNK.HELM2D.GREEN - -src = srcinfo.r; -targ = targinfo.r; - -[~,ns] = size(src); -[~,nt] = size(targ); - -% double layer -if strcmpi(type,'d') - srcnorm = srcinfo.n; - [~,grad] = chnk.helm2d.green(zk,src,targ); - nx = repmat(srcnorm(1,:),nt,1); - ny = repmat(srcnorm(2,:),nt,1); - submat = -(grad(:,:,1).*nx + grad(:,:,2).*ny); -end - -% normal derivative of single layer -if strcmpi(type,'sprime') - targnorm = targinfo.n; - [~,grad] = chnk.helm2d.green(zk,src,targ); - nx = repmat((targnorm(1,:)).',1,ns); - ny = repmat((targnorm(2,:)).',1,ns); - - submat = (grad(:,:,1).*nx + grad(:,:,2).*ny); -end - - -% Tangential derivative of single layer -if strcmpi(type,'stau') - targtan = targinfo.d; - [~,grad] = chnk.helm2d.green(zk,src,targ); - dx = repmat((targtan(1,:)).',1,ns); - dy = repmat((targtan(2,:)).',1,ns); - ds = sqrt(dx.*dx+dy.*dy); - submat = (grad(:,:,1).*dx + grad(:,:,2).*dy)./ds; -end - -% single layer -if strcmpi(type,'s') - submat = chnk.helm2d.green(zk,src,targ); -end - -% normal derivative of double layer -if strcmpi(type,'dprime') - targnorm = targinfo.n; - srcnorm = srcinfo.n; - [~,~,hess] = chnk.helm2d.green(zk,src,targ); - nxsrc = repmat(srcnorm(1,:),nt,1); - nysrc = repmat(srcnorm(2,:),nt,1); - nxtarg = repmat((targnorm(1,:)).',1,ns); - nytarg = repmat((targnorm(2,:)).',1,ns); - submat = -(hess(:,:,1).*nxsrc.*nxtarg + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)... - + hess(:,:,3).*nysrc.*nytarg); -end - -% Combined field -if strcmpi(type,'c') - srcnorm = srcinfo.n; - coef = ones(2,1); - if(nargin == 5); coef = varargin{1}; end - [submats,grad] = chnk.helm2d.green(zk,src,targ); - nx = repmat(srcnorm(1,:),nt,1); - ny = repmat(srcnorm(2,:),nt,1); - submatd = -(grad(:,:,1).*nx + grad(:,:,2).*ny); - submat = coef(1)*submatd + coef(2)*submats; -end - -% normal derivative of combined field -if strcmpi(type,'cprime') - coef = ones(2,1); - if(nargin == 5); coef = varargin{1}; end - targnorm = targinfo.n; - srcnorm = srcinfo.n; - - - - % Get gradient and hessian info - [~,grad,hess] = chnk.helm2d.green(zk,src,targ); - - nxsrc = repmat(srcnorm(1,:),nt,1); - nysrc = repmat(srcnorm(2,:),nt,1); - nxtarg = repmat((targnorm(1,:)).',1,ns); - nytarg = repmat((targnorm(2,:)).',1,ns); - - % D' - submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ... - + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)... - + hess(:,:,3).*nysrc.*nytarg); - % S' - submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg); - - submat = coef(1)*submatdp + coef(2)*submatsp; -end - - -% Dirichlet and neumann data corresponding to combined field -if strcmpi(type,'c2trans') - coef = ones(2,1); - if(nargin == 5); coef = varargin{1}; end - targnorm = targinfo.n; - srcnorm = srcinfo.n; - - % Get gradient and hessian info - [submats,grad,hess] = chnk.helm2d.green(zk,src,targ); - - nxsrc = repmat(srcnorm(1,:),nt,1); - nysrc = repmat(srcnorm(2,:),nt,1); - nxtarg = repmat((targnorm(1,:)).',1,ns); - nytarg = repmat((targnorm(2,:)).',1,ns); - - % D' - submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ... - + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)... - + hess(:,:,3).*nysrc.*nytarg); - % S' - submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg); - % D - submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc); - - 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; - -end - - -% all kernels, [c11 D, c12 S; c21 D', c22 S'] -if strcmpi(type,'all') - - targnorm = targinfo.n; - srcnorm = srcinfo.n; - cc = varargin{1}; - - submat = zeros(2*nt,2*ns); - - % Get gradient and hessian info - [submats,grad,hess] = chnk.helm2d.green(zk,src,targ); - - nxsrc = repmat(srcnorm(1,:),nt,1); - nysrc = repmat(srcnorm(2,:),nt,1); - nxtarg = repmat((targnorm(1,:)).',1,ns); - nytarg = repmat((targnorm(2,:)).',1,ns); - - submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ... - + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)... - + hess(:,:,3).*nysrc.*nytarg); - % S' - submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg); - % D - submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc); - - - submat(1:2:2*nt,1:2:2*ns) = submatd*cc(1,1); - submat(1:2:2*nt,2:2:2*ns) = submats*cc(1,2); - submat(2:2:2*nt,1:2:2*ns) = submatdp*cc(2,1); - submat(2:2:2*nt,2:2:2*ns) = submatsp*cc(2,2); -end - -% Dirichlet data/potential correpsonding to transmission rep -if strcmpi(type,'trans_rep') - - coef = ones(2,1); - if(nargin == 5); coef = varargin{1}; end; - srcnorm = srcinfo.n; - [submats,grad] = chnk.helm2d.green(zk,src,targ); - nx = repmat(srcnorm(1,:),nt,1); - ny = repmat(srcnorm(2,:),nt,1); - submatd = -(grad(:,:,1).*nx + grad(:,:,2).*ny); - - submat = zeros(1*nt,2*ns); - submat(1:1:1*nt,1:2:2*ns) = coef(1)*submatd; - submat(1:1:1*nt,2:2:2*ns) = coef(2)*submats; -end - -% Neumann data corresponding to transmission rep -if strcmpi(type,'trans_rep_prime') - coef = ones(2,1); - if(nargin == 5); coef = varargin{1}; end; - targnorm = targinfo.n; - srcnorm = srcinfo.n; - - submat = zeros(nt,ns); - - % Get gradient and hessian info - [submats,grad,hess] = chnk.helm2d.green(zk,src,targ); - - nxsrc = repmat(srcnorm(1,:),nt,1); - nysrc = repmat(srcnorm(2,:),nt,1); - nxtarg = repmat((targnorm(1,:)).',1,ns); - nytarg = repmat((targnorm(2,:)).',1,ns); - - % D' - submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ... - + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)... - + hess(:,:,3).*nysrc.*nytarg); - % S' - submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg); - submat = zeros(nt,2*ns) - submat(1:1:1*nt,1:2:2*ns) = coef(1)*submatdp; - submat(1:1:1*nt,2:2:2*ns) = coef(2)*submatsp; -end - - -% Gradient correpsonding to transmission rep -if strcmpi(type,'trans_rep_grad') - coef = ones(2,1); - if(nargin == 5); coef = varargin{1}; end; - - srcnorm = srcinfo.n; - - submat = zeros(nt,ns,6); - % S - [submats,grad,hess] = chnk.helm2d.green(zk,src,targ); - - nxsrc = repmat(srcnorm(1,:),nt,1); - nysrc = repmat(srcnorm(2,:),nt,1); - % D - submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc); - - submat = zeros(2*nt,2*ns); - - submat(1:2:2*nt,1:2:2*ns) = -coef(1)*(hess(:,:,1).*nxsrc + hess(:,:,2).*nysrc); - submat(1:2:2*nt,2:2:2*ns) = coef(2)*grad(:,:,1); - - submat(2:2:2*nt,1:2:2*ns) = -coef(1)*(hess(:,:,2).*nxsrc + hess(:,:,3).*nysrc); - submat(2:2:2*nt,2:2:2*ns) = coef(2)*grad(:,:,2); -end - - +function submat= kern(zk,srcinfo,targinfo,type,varargin) +%CHNK.HELM2D.KERN standard Helmholtz layer potential kernels in 2D +% +% Syntax: submat = chnk.heml2d.kern(zk,srcinfo,targingo,type,varargin) +% +% Let x be targets and y be sources for these formulas, with +% n_x and n_y the corresponding unit normals at those points +% (if defined). Note that the normal information is obtained +% by taking the perpendicular to the provided tangential deriviative +% info and normalizing +% +% Kernels based on G(x,y) = i/4 H_0^{(1)}(zk |x-y|) +% +% D(x,y) = \nabla_{n_y} G(x,y) +% S(x,y) = G(x,y) +% S'(x,y) = \nabla_{n_x} G(x,y) +% D'(x,y) = \nabla_{n_x} \nabla_{n_y} G(x,y) +% +% Input: +% zk - complex number, Helmholtz wave number +% srcinfo - description of sources in ptinfo struct format, i.e. +% ptinfo.r - positions (2,:) array +% ptinfo.d - first derivative in underlying +% parameterization (2,:) +% ptinfo.d2 - second derivative in underlying +% parameterization (2,:) +% targinfo - description of targets in ptinfo struct format, +% if info not relevant (d/d2) it doesn't need to +% be provided. sprime requires tangent info in +% targinfo.d +% type - string, determines kernel type +% type == 'd', double layer kernel D +% type == 's', single layer kernel S +% type == 'sprime', normal derivative of single +% layer S' +% type == 'dprime', normal derivative of double layer D' +% type == 'c', combined layer kernel coef(1) D + coef(2) S +% type == 'stau', tangential derivative of single layer +% type == 'all', returns all four layer potentials, +% [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'] +% type == 'trans_rep' returns the potential corresponding +% to the transmission representation +% [coef(1)*D coef(2)*S] +% type == 'trans_rep_prime' returns the normal derivative +% corresponding to the transmission representation +% [coef(1)*D' coef(2)*S'] +% type == 'trans_rep_grad' returns the gradient corresponding +% to the transmission representation +% [coef(1)*d_x D coef(2)*d_x S; +% coef(1)*d_y D coef(2)*d_y S] +% +% varargin{1} - coef: length 2 array in the combined layer +% formula, 2x2 matrix for all kernels +% otherwise does nothing +% +% Output: +% submat - the evaluation of the selected kernel for the +% provided sources and targets. the number of +% rows equals the number of targets and the +% number of columns equals the number of sources +% +% see also CHNK.HELM2D.GREEN + +src = srcinfo.r; +targ = targinfo.r; + +[~,ns] = size(src); +[~,nt] = size(targ); + +% double layer +if strcmpi(type,'d') + srcnorm = srcinfo.n; + [~,grad] = chnk.helm2d.green(zk,src,targ); + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + submat = -(grad(:,:,1).*nx + grad(:,:,2).*ny); +end + +% normal derivative of single layer +if strcmpi(type,'sprime') + targnorm = targinfo.n; + [~,grad] = chnk.helm2d.green(zk,src,targ); + nx = repmat((targnorm(1,:)).',1,ns); + ny = repmat((targnorm(2,:)).',1,ns); + submat = (grad(:,:,1).*nx + grad(:,:,2).*ny); +end + + +% Tangential derivative of single layer +if strcmpi(type,'stau') + targtan = targinfo.d; + [~,grad] = chnk.helm2d.green(zk,src,targ); + dx = repmat((targtan(1,:)).',1,ns); + dy = repmat((targtan(2,:)).',1,ns); + ds = sqrt(dx.*dx+dy.*dy); + submat = (grad(:,:,1).*dx + grad(:,:,2).*dy)./ds; +end + +% single layer +if strcmpi(type,'s') + submat = chnk.helm2d.green(zk,src,targ); +end + +% normal derivative of double layer +if strcmpi(type,'dprime') + targnorm = targinfo.n; + srcnorm = srcinfo.n; + [~,~,hess] = chnk.helm2d.green(zk,src,targ); + nxsrc = repmat(srcnorm(1,:),nt,1); + nysrc = repmat(srcnorm(2,:),nt,1); + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + submat = -(hess(:,:,1).*nxsrc.*nxtarg + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)... + + hess(:,:,3).*nysrc.*nytarg); +end + +% Combined field +if strcmpi(type,'c') + srcnorm = srcinfo.n; + coef = ones(2,1); + if(nargin == 5); coef = varargin{1}; end + [submats,grad] = chnk.helm2d.green(zk,src,targ); + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + submatd = -(grad(:,:,1).*nx + grad(:,:,2).*ny); + submat = coef(1)*submatd + coef(2)*submats; +end + +% normal derivative of combined field +if strcmpi(type,'cprime') + coef = ones(2,1); + if(nargin == 5); coef = varargin{1}; end; + targnorm = targinfo.n; + srcnorm = srcinfo.n; + + submat = zeros(nt,ns); + + % Get gradient and hessian info + [submats,grad,hess] = chnk.helm2d.green(zk,src,targ); + + nxsrc = repmat(srcnorm(1,:),nt,1); + nysrc = repmat(srcnorm(2,:),nt,1); + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + % D' + submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ... + + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)... + + hess(:,:,3).*nysrc.*nytarg); + % S' + submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg); + + submat = coef(1)*submatdp + coef(2)*submatsp; +end + + +% Dirichlet and neumann data corresponding to combined field +if strcmpi(type,'c2trans') + coef = ones(2,1); + if(nargin == 5); coef = varargin{1}; end; + targnorm = targinfo.n; + srcnorm = srcinfo.n; + + % Get gradient and hessian info + [submats,grad,hess] = chnk.helm2d.green(zk,src,targ); + + nxsrc = repmat(srcnorm(1,:),nt,1); + nysrc = repmat(srcnorm(2,:),nt,1); + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + % D' + submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ... + + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)... + + hess(:,:,3).*nysrc.*nytarg); + % S' + submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg); + % D + submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc); + + 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; + +end + + +% all kernels, [c11 D, c12 S; c21 D', c22 S'] +if strcmpi(type,'all') + + targnorm = targinfo.n; + srcnorm = srcinfo.n; + cc = varargin{1}; + + submat = zeros(2*nt,2*ns); + + % Get gradient and hessian info + [submats,grad,hess] = chnk.helm2d.green(zk,src,targ); + + nxsrc = repmat(srcnorm(1,:),nt,1); + nysrc = repmat(srcnorm(2,:),nt,1); + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ... + + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)... + + hess(:,:,3).*nysrc.*nytarg); + % S' + submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg); + % D + submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc); + + + submat(1:2:2*nt,1:2:2*ns) = submatd*cc(1,1); + submat(1:2:2*nt,2:2:2*ns) = submats*cc(1,2); + submat(2:2:2*nt,1:2:2*ns) = submatdp*cc(2,1); + submat(2:2:2*nt,2:2:2*ns) = submatsp*cc(2,2); +end + + +% clamped plate kernel for modified biharmonic problem +if strcmpi(type, 'clamped-plate') + srcnorm = srcinfo.n; + srctang = srcinfo.d; + targnorm = targinfo.n; + + [hess, third, forth] = chnk.helm2d.thirdforth_derivatives(zk, src, targ); % Hankel part + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + zkimag = (1i)*zk; + [hessK, thirdK, forthK] = chnk.helm2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + ds = sqrt(dx.*dx+dy.*dy); + + taux = dx./ds; + tauy = dy./ds; + + Kxx = 1/(2*zk^2).*(third(:, :, 1).*(nx.*nx.*nx) + third(:, :, 2).*(3*nx.*nx.*ny) +... + third(:, :, 4).*(3*nx.*ny.*ny) + third(:, :, 6).*(ny.*ny.*ny)) - ... + 1/(2*zk^2).*(thirdK(:, :, 1).*(nx.*nx.*nx) + thirdK(:, :, 2).*(3*nx.*nx.*ny) +... + thirdK(:, :, 4).*(3*nx.*ny.*ny) + thirdK(:, :, 6).*(ny.*ny.*ny)) + ... + (3/(2*zk^2).*(third(:, :, 1).*(nx.*taux.*taux) + third(:, :, 2).*(2*nx.*taux.*tauy + ny.*taux.*taux) +... + third(:, :, 4).*(nx.*tauy.*tauy + 2*ny.*taux.*tauy) + third(:, :, 6).*(ny.*tauy.*tauy))-... + 3/(2*zk^2).*(thirdK(:, :, 1).*(nx.*taux.*taux) + thirdK(:, :, 2).*(2*nx.*taux.*tauy + ny.*taux.*taux) +... + thirdK(:, :, 4).*(nx.*tauy.*tauy + 2*ny.*taux.*tauy) + thirdK(:, :, 6).*(ny.*tauy.*tauy))); % G_{ny ny ny} + 3G_{ny tauy tauy} + + Kxy = -(1/(2*zk^2).*(hess(:, :, 1).*(nx.*nx) + hess(:, :, 2).*(2*nx.*ny) + hess(:, :, 3).*(ny.*ny))-... + 1/(2*zk^2).*(hessK(:, :, 1).*(nx.*nx) + hessK(:, :, 2).*(2*nx.*ny) + hessK(:, :, 3).*(ny.*ny)))+... + (1/(2*zk^2).*(hess(:, :, 1).*(taux.*taux) + hess(:, :, 2).*(2*taux.*tauy) + hess(:, :, 3).*(tauy.*tauy))-... + 1/(2*zk^2).*(hessK(:, :, 1).*(taux.*taux) + hessK(:, :, 2).*(2*taux.*tauy) + hessK(:, :, 3).*(tauy.*tauy))); % -G_{ny ny} + G_{tauy tauy} + + Kyx = 1/(2*zk^2).*(forth(:, :, 1).*(nx.*nx.*nx.*nxtarg) + forth(:, :, 2).*(nx.*nx.*nx.*nytarg + 3*nx.*nx.*ny.*nxtarg) + ... + forth(:, :, 4).*(3*nx.*nx.*ny.*nytarg + 3*nx.*ny.*ny.*nxtarg) + forth(:, :, 6).*(3*nx.*ny.*ny.*nytarg +ny.*ny.*ny.*nxtarg)+... + forth(:, :, 8).*(ny.*ny.*ny.*nytarg)) - 1/(2*zk^2).*(forthK(:, :, 1).*(nx.*nx.*nx.*nxtarg) + forthK(:, :, 2).*(nx.*nx.*nx.*nytarg + 3*nx.*nx.*ny.*nxtarg) + ... + forthK(:, :, 4).*(3*nx.*nx.*ny.*nytarg + 3*nx.*ny.*ny.*nxtarg) + forthK(:, :, 6).*(3*nx.*ny.*ny.*nytarg +ny.*ny.*ny.*nxtarg)+... + forthK(:, :, 8).*(ny.*ny.*ny.*nytarg)) + ... + (3/(2*zk^2).*(forth(:, :, 1).*(nx.*taux.*taux.*nxtarg)+ forth(:, :, 2).*(nx.*taux.*taux.*nytarg + 2*nx.*taux.*tauy.*nxtarg + ny.*taux.*taux.*nxtarg) +... + forth(:, :, 4).*(2*nx.*taux.*tauy.*nytarg + ny.*taux.*taux.*nytarg + nx.*tauy.*tauy.*nxtarg + 2*ny.*taux.*tauy.*nxtarg) + ... + forth(:, :, 6).*(nx.*tauy.*tauy.*nytarg +2*ny.*taux.*tauy.*nytarg + ny.*tauy.*tauy.*nxtarg) + forth(:, :, 8).*(ny.*tauy.*tauy.*nytarg))-... + 3/(2*zk^2).*(forthK(:, :, 1).*(nx.*taux.*taux.*nxtarg)+ forthK(:, :, 2).*(nx.*taux.*taux.*nytarg + 2*nx.*taux.*tauy.*nxtarg + ny.*taux.*taux.*nxtarg) +... + forthK(:, :, 4).*(2*nx.*taux.*tauy.*nytarg + ny.*taux.*taux.*nytarg + nx.*tauy.*tauy.*nxtarg + 2*ny.*taux.*tauy.*nxtarg) + ... + forthK(:, :, 6).*(nx.*tauy.*tauy.*nytarg +2*ny.*taux.*tauy.*nytarg + ny.*tauy.*tauy.*nxtarg) + forthK(:, :, 8).*(ny.*tauy.*tauy.*nytarg))); + + Kyy = 1/(2*zk^2).*(third(:,:, 1).*(nx.*nx.*nxtarg) +third(:, :, 2).*(nx.*nx.*nytarg + 2*nx.*ny.*nxtarg) + third(:, :, 4).*(2*nx.*ny.*nytarg + ny.*ny.*nxtarg)+... + third(:, :,6).*(ny.*ny.*nytarg)) - 1/(2*zk^2).*(thirdK(:,:, 1).*(nx.*nx.*nxtarg) +thirdK(:, :, 2).*(nx.*nx.*nytarg + 2*nx.*ny.*nxtarg) + thirdK(:, :, 4).*(2*nx.*ny.*nytarg + ny.*ny.*nxtarg)+... + thirdK(:, :,6).*(ny.*ny.*nytarg)) - 1/(2*zk^2).*(third(:,:, 1).*(taux.*taux.*nxtarg) +third(:, :, 2).*(taux.*taux.*nytarg + 2*taux.*tauy.*nxtarg) + third(:, :, 4).*(2*taux.*tauy.*nytarg + tauy.*tauy.*nxtarg)+... + third(:, :,6).*(tauy.*tauy.*nytarg)) + 1/(2*zk^2).*(thirdK(:,:, 1).*(taux.*taux.*nxtarg) +thirdK(:, :, 2).*(taux.*taux.*nytarg + 2*taux.*tauy.*nxtarg) + thirdK(:, :, 4).*(2*taux.*tauy.*nytarg + tauy.*tauy.*nxtarg)+... + thirdK(:, :,6).*(tauy.*tauy.*nytarg)); + + submat = zeros(2*nt,2*ns); + + submat(1:2:end,1:2:end) = Kxx; + submat(1:2:end,2:2:end) = Kxy; + + submat(2:2:end,1:2:end) = Kyx; + submat(2:2:end,2:2:end) = Kyy; +end + +% free plate kernel for the modified biharmonic problem (K11 with no +% hilbert transform subtraction, K12 kernel, K22 with no curvature part) K21 is +% handled in a separate type. +if strcmpi(type, 'free plate first part') + srcnorm = srcinfo.n; + srctang = srcinfo.d; + targnorm = targinfo.n; + targtang = targinfo.d; + coefs = varargin{1}; + + + + + [hess, third, ~] = chnk.helm2d.thirdforth_derivatives(zk, src, targ); % Hankel part + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + zkimag = (1i)*zk; + [hessK, thirdK, ~] = chnk.helm2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + + + + + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + dx1 = repmat((targtang(1,:)).',1,ns); + dy1 = repmat((targtang(2,:)).',1,ns); + + + ds = sqrt(dx.*dx+dy.*dy); + ds1 = sqrt(dx1.*dx1+dy1.*dy1); + + taux = dx./ds; % normalization + tauy = dy./ds; + + tauxtarg = dx1./ds1; + tauytarg = dy1./ds1; + + + + + K11 = 1/(2*zk^2).*(third(:, :, 1).*(nxtarg.*nxtarg.*nx) + third(:, :, 2).*(nxtarg.*nxtarg.*ny + 2*nxtarg.*nytarg.*nx) +... + third(:, :, 4).*(2*nxtarg.*nytarg.*ny + nytarg.*nytarg.*nx) +... + third(:, :, 6).*(nytarg.*nytarg.*ny)) - ... + 1/(2*zk^2).*(thirdK(:, :, 1).*(nxtarg.*nxtarg.*nx) + thirdK(:, :, 2).*(nxtarg.*nxtarg.*ny + 2*nxtarg.*nytarg.*nx) +... + thirdK(:, :, 4).*(2*nxtarg.*nytarg.*ny + nytarg.*nytarg.*nx) +... + thirdK(:, :, 6).*(nytarg.*nytarg.*ny)) + ... + coefs(1)/(2*zk^2).*(third(:, :, 1).*(tauxtarg.*tauxtarg.*nx) + third(:, :, 2).*(tauxtarg.*tauxtarg.*ny + 2*tauxtarg.*tauytarg.*nx) +... + third(:, :, 4).*(2*tauxtarg.*tauytarg.*ny + tauytarg.*tauytarg.*nx) +... + third(:, :, 6).*(tauytarg.*tauytarg.*ny)) - ... + coefs(1)/(2*zk^2).*(thirdK(:, :, 1).*(tauxtarg.*tauxtarg.*nx) + thirdK(:, :, 2).*(tauxtarg.*tauxtarg.*ny + 2*tauxtarg.*tauytarg.*nx) +... + thirdK(:, :, 4).*(2*tauxtarg.*tauytarg.*ny + tauytarg.*tauytarg.*nx) +... + thirdK(:, :, 6).*(tauytarg.*tauytarg.*ny)); % first kernel with no hilbert transforms (G_{nx nx ny + nu G_{taux taux ny}). + + + + + K12 = 1/(2*zk^2).*(hess(:, :, 1).*nxtarg.*nxtarg + hess(:, :, 2).*(2*nxtarg.*nytarg) + hess(:, :, 3).*nytarg.*nytarg)-... + 1/(2*zk^2).*(hessK(:, :, 1).*nxtarg.*nxtarg + hessK(:, :, 2).*(2*nxtarg.*nytarg) + hessK(:, :, 3).*nytarg.*nytarg)+... + coefs(1)/(2*zk^2).*(hess(:, :, 1).*tauxtarg.*tauxtarg + hess(:, :, 2).*(2*tauxtarg.*tauytarg) + hess(:, :, 3).*tauytarg.*tauytarg)-... + coefs(1)/(2*zk^2).*(hessK(:, :, 1).*tauxtarg.*tauxtarg + hessK(:, :, 2).*(2*tauxtarg.*tauytarg) + ... + hessK(:, :, 3).*tauytarg.*tauytarg) ; % G_{nx nx} + nu G_{taux taux} + + K21 = 0; + + + + + K22 = -1./(2*zk^2).*(third(:, :, 1).*(nxtarg.*nxtarg.*nxtarg) + third(:, :, 2).*(3*nxtarg.*nxtarg.*nytarg) +... + third(:, :, 4).*(3*nxtarg.*nytarg.*nytarg) + third(:, :, 6).*(nytarg.*nytarg.*nytarg)) + ... + 1./(2*zk^2).*(thirdK(:, :, 1).*(nxtarg.*nxtarg.*nxtarg) + thirdK(:, :, 2).*(3*nxtarg.*nxtarg.*nytarg)+... + thirdK(:, :, 4).*(3*nxtarg.*nytarg.*nytarg) + thirdK(:, :, 6).*(nytarg.*nytarg.*nytarg)) -... + (2-coefs(1))/(2*zk^2).*(third(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg) + third(:, :, 2).*(tauxtarg.*tauxtarg.*nytarg + 2*tauxtarg.*tauytarg.*nxtarg) +... + third(:, :, 4).*(2*tauxtarg.*tauytarg.*nytarg + tauytarg.*tauytarg.*nxtarg) +... + + third(:, :, 6).*(tauytarg.*tauytarg.*nytarg)) + ... + (2-coefs(1))/(2*zk^2).*(thirdK(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg) + thirdK(:, :, 2).*(tauxtarg.*tauxtarg.*nytarg + 2*tauxtarg.*tauytarg.*nxtarg) +... + thirdK(:, :, 4).*(2*tauxtarg.*tauytarg.*nytarg + tauytarg.*tauytarg.*nxtarg) +... + + thirdK(:, :, 6).*(tauytarg.*tauytarg.*nytarg)); % G_{nx nx nx} + (2-nu) G_{taux taux nx} + + submat = zeros(2*nt,2*ns); + + submat(1:2:end,1:2:end) = K11; + submat(1:2:end,2:2:end) = K12; + + submat(2:2:end,1:2:end) = K21; + submat(2:2:end,2:2:end) = K22; + +end + + +% free plate kernel K21 for the modified biharmonic problem. This part +% handles the singularity subtraction and swap the evaluation to its +% asymptotic expansions if the targets and sources are close. + +if strcmpi(type, 'free plate K21 first part') + srcnorm = srcinfo.n; + srctang = srcinfo.d; + targnorm = targinfo.n; + targtang = targinfo.d; + coefs = varargin{1}; + R = 0.001; % radius of kernel replacement + + + + [~, ~, forth] = chnk.helm2d.thirdforth_derivatives(zk, src, targ); % Hankel part + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + zkimag = (1i)*zk; + [~, ~, forthK] = chnk.helm2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + + + + + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + dx1 = repmat((targtang(1,:)).',1,ns); + dy1 = repmat((targtang(2,:)).',1,ns); + + + ds = sqrt(dx.*dx+dy.*dy); + ds1 = sqrt(dx1.*dx1+dy1.*dy1); + + taux = dx./ds; % normalization + tauy = dy./ds; + + tauxtarg = dx1./ds1; + tauytarg = dy1./ds1; + + rx = targ(1,:).' - src(1,:); + ry = targ(2,:).' - src(2,:); + r2 = rx.^2 + ry.^2; + dists = sqrt(r2); + inds = and((dists < R),(dists > 0)); + + submat = 1/(2*zk^2).*(forth(:, :, 1).*(nxtarg.*nxtarg.*nxtarg.*nx) + forth(:, :, 2).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... + forth(:, :, 4).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) + forth(:, :, 6).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... + forth(:, :, 8).*(nytarg.*nytarg.*nytarg.*ny)) -... + 1/(2*zk^2).*(forthK(:, :, 1).*(nxtarg.*nxtarg.*nxtarg.*nx) + forthK(:, :, 2).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... + forthK(:, :, 4).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) + forthK(:, :, 6).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... + forthK(:, :, 8).*(nytarg.*nytarg.*nytarg.*ny)) + (3/(4*pi)).*(nxtarg.*nx + nytarg.*ny)./r2 + ... + ((2-coefs(1))/(2*zk^2).*(forth(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg.*nx) + forth(:, :, 2).*(tauxtarg.*tauxtarg.*nxtarg.*ny + tauxtarg.*tauxtarg.*nytarg.*nx + 2*tauxtarg.*tauytarg.*nxtarg.*nx)+... + forth(:, :, 4).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... + forth(:, :, 6).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... + forth(:, :, 8).*(tauytarg.*tauytarg.*nytarg.*ny)) -... + (2-coefs(1))/(2*zk^2).*(forthK(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg.*nx) + forthK(:, :, 2).*(tauxtarg.*tauxtarg.*nxtarg.*ny + tauxtarg.*tauxtarg.*nytarg.*nx + 2*tauxtarg.*tauytarg.*nxtarg.*nx)+... + forthK(:, :, 4).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... + forthK(:, :, 6).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... + forthK(:, :, 8).*(tauytarg.*tauytarg.*nytarg.*ny)) +... + (2-coefs(1))/(4*pi).*(nxtarg.*nx + nytarg.*ny)./r2 - ... + ((2-coefs(1))/(2*pi)).*(nxtarg.*nx + nytarg.*ny).*(rx.*tauxtarg + ry.*tauytarg).*(rx.*tauxtarg + ry.*tauytarg)./(r2.^2)) - ... + ((1+coefs(1))/(4*pi)).*((taux.*tauxtarg + tauy.*tauytarg)./(r2) - 2*(rx.*tauxtarg + ry.*tauytarg).*(rx.*taux + ry.*tauy)./(r2.^2))-... + (3/(4*pi)).*(nxtarg.*nx + nytarg.*ny)./r2 - (2-coefs(1))/(4*pi).*(nxtarg.*nx + nytarg.*ny)./r2 + ... + (2-coefs(1))/(2*pi).*(nxtarg.*nx + nytarg.*ny).*(rx.*tauxtarg + ry.*tauytarg).*(rx.*tauxtarg + ry.*tauytarg)./(r2.^2); + % third kernel with singularity subtraction and H[delta'] + + if any(inds(:)) % if statement for evalauation points that are closed to the boundary + temp = 1i*imag(submat(inds)); + + nx = nx(inds); ny = ny(inds); + + nxtarg = nxtarg(inds); + nytarg = nytarg(inds); + + dx = dx(inds); + dy = dy(inds); + + dx1 = dx1(inds); + dy1 = dy1(inds); + + ds = sqrt(dx.*dx+dy.*dy); + ds1 = sqrt(dx1.*dx1+dy1.*dy1); + + taux = dx./ds; % normalization + tauy = dy./ds; + + tauxtarg = dx1./ds1; + tauytarg = dy1./ds1; + + rx = rx(inds); + ry = ry(inds); + r2 = r2(inds); + + tauxtargnsrc = tauxtarg.*nx + tauytarg.*ny; + tauxtargntarg = tauxtarg.*nxtarg + tauytarg.*nytarg; + + normaldot = nx.*nxtarg + ny.*nytarg; + tangentdot = taux.*tauxtarg + tauy.*tauytarg; + + rn = rx.*nx + ry.*ny; + + rntarg = rx.*nxtarg + ry.*nytarg; + + rtau = rx.*taux + ry.*tauy; + + rtautarg = rx.*tauxtarg + ry.*tauytarg; + + + + [smoothfirst, smoothsecond] = chnk.helm2d.smooth_part_derivatives(zk, rx, ry); + + eulerconstantpart = (smoothfirst(:, :, 1)).*(nxtarg.*nxtarg.*nxtarg.*nx) +... + (smoothfirst(:, :, 2)).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... + (smoothfirst(:, :, 3)).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) +... + (smoothfirst(:, :, 4)).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... + (smoothfirst(:, :, 5)).*(nytarg.*nytarg.*nytarg.*ny)+ ... + (2-coefs(1)).*((smoothfirst(:, :, 1)).*(tauxtarg.*tauxtarg.*nxtarg.*nx) +... + (smoothfirst(:, :, 2)).*(tauxtarg.*tauxtarg.*nxtarg.*ny + tauxtarg.*tauxtarg.*nytarg.*nx + 2*tauxtarg.*tauytarg.*nxtarg.*nx)+... + (smoothfirst(:, :, 3)).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... + (smoothfirst(:, :, 4)).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... + (smoothfirst(:, :, 5)).*(tauytarg.*tauytarg.*nytarg.*ny)); + + puresmoothpart = (smoothsecond(:, :, 1)).*(nxtarg.*nxtarg.*nxtarg.*nx) +... + (smoothsecond(:, :, 2)).*(nxtarg.*nxtarg.*nxtarg.*ny + 3*nxtarg.*nxtarg.*nytarg.*nx) + ... + (smoothsecond(:, :, 3)).*(3*nxtarg.*nxtarg.*nytarg.*ny + 3*nxtarg.*nytarg.*nytarg.*nx) +... + (smoothsecond(:, :, 4)).*(3*nxtarg.*nytarg.*nytarg.*ny +nytarg.*nytarg.*nytarg.*nx)+... + (smoothsecond(:, :, 5)).*(nytarg.*nytarg.*nytarg.*ny)+ ... + (2-coefs(1)).*((smoothsecond(:, :, 1)).*(tauxtarg.*tauxtarg.*nxtarg.*nx) +... + (smoothsecond(:, :, 2)).*(tauxtarg.*tauxtarg.*nxtarg.*ny + tauxtarg.*tauxtarg.*nytarg.*nx + 2*tauxtarg.*tauytarg.*nxtarg.*nx)+... + (smoothsecond(:, :, 3)).*(tauxtarg.*tauxtarg.*nytarg.*ny + 2*tauxtarg.*tauytarg.*nxtarg.*ny + tauytarg.*tauytarg.*nxtarg.*nx + 2*tauxtarg.*tauytarg.*nytarg.*nx)+... + (smoothsecond(:, :, 4)).*(tauytarg.*tauytarg.*nxtarg.*ny + 2*tauxtarg.*tauytarg.*nytarg.*ny + tauytarg.*tauytarg.*nytarg.*nx) +... + (smoothsecond(:, :, 5)).*(tauytarg.*tauytarg.*nytarg.*ny)); + + submat(inds) = 3/(2*pi).*(rn.*rntarg./(r2.^2)) + 3/(2*pi).*(rntarg.^2)./(r2.^2) -... + (2/pi).*(rn.*(rntarg.^3))./(r2.^3) + ... + (2-coefs(1)).*(-1/(2*pi).*(tauxtargnsrc).*(tauxtargntarg)./(r2) +... + 1/(pi).*(rn.*rtautarg.*tauxtargntarg)./(r2.^2) + ... + 1/(pi).*(rtautarg.*tauxtargnsrc.*rntarg)./(r2.^2) -.... + 2/(pi).*(rn.*rntarg.*(rtautarg.*rtautarg))./(r2.^3) + ... + 1/(2*pi).*(rn.*rntarg)./(r2.^2)) - 3/(4*pi).*normaldot./r2 +... + ((2-coefs(1))/(2*pi)).*((rtautarg.*rtautarg).*(normaldot))./(r2.^2) -... + ((2-coefs(1))/(4*pi)).*(normaldot)./r2 - ... + ((1+coefs(1))/(4*pi)).*((tangentdot)./(r2) - 2*(rtautarg).*(rtau)./(r2.^2)) + temp + eulerconstantpart + puresmoothpart; + end +end + +% kernels in K11 with hilbert transform subtractions. +% (i.e. beta*(G_{nx nx tauy} + 1/4 H + nu*(G_{taux taux tauy} + 1/4 H)) + +if strcmpi(type, 'free plate hilbert subtract') + srcnorm = srcinfo.n; + srctang = srcinfo.d; + targnorm = targinfo.n; + targtang = targinfo.d; + coefs = varargin{1}; + + [~, third, ~] = chnk.helm2d.thirdforth_derivatives(zk, src, targ); % Hankel part + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + zkimag = (1i)*zk; + [~, thirdK, ~] = chnk.helm2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + [~,grad] = chnk.lap2d.green(src,targ,true); + + + + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + dx1 = repmat((targtang(1,:)).',1,ns); + dy1 = repmat((targtang(2,:)).',1,ns); + + + ds = sqrt(dx.*dx+dy.*dy); + ds1 = sqrt(dx1.*dx1+dy1.*dy1); + + taux = dx./ds; % normalization + tauy = dy./ds; + + tauxtarg = dx1./ds1; + tauytarg = dy1./ds1; + + + hilb = 2*(grad(:,:,1).*ny - grad(:,:,2).*nx); + + submat = ((1+ coefs(1))/2)*(1./(2*zk^2).*(third(:, :, 1).*(nxtarg.*nxtarg.*taux) + third(:, :, 2).*(nxtarg.*nxtarg.*tauy+ 2*nxtarg.*nytarg.*taux) +... + third(:, :, 4).*(2*nxtarg.*nytarg.*tauy +nytarg.*nytarg.*taux) +... + third(:, :, 6).*(nytarg.*nytarg.*tauy)) - ... + 1./(2*zk^2).*(thirdK(:, :, 1).*(nxtarg.*nxtarg.*taux) + thirdK(:, :, 2).*(nxtarg.*nxtarg.*tauy+ 2*nxtarg.*nytarg.*taux) +... + thirdK(:, :, 4).*(2*nxtarg.*nytarg.*tauy +nytarg.*nytarg.*taux) +... + thirdK(:, :, 6).*(nytarg.*nytarg.*tauy)) + 0.25*hilb) +... + ((1+ coefs(1))/2)*coefs(1).*(1./(2*zk^2).*(third(:, :, 1).*(tauxtarg.*tauxtarg.*taux) + third(:, :, 2).*(tauxtarg.*tauxtarg.*tauy + 2*tauxtarg.*tauytarg.*taux) +... + third(:, :, 4).*(2*tauxtarg.*tauytarg.*tauy + tauytarg.*tauytarg.*taux) +... + third(:, :, 6).*(tauytarg.*tauytarg.*tauy)) - ... + 1./(2*zk^2).*(thirdK(:, :, 1).*(tauxtarg.*tauxtarg.*taux) + thirdK(:, :, 2).*(tauxtarg.*tauxtarg.*tauy + 2*tauxtarg.*tauytarg.*taux) +... + thirdK(:, :, 4).*(2*tauxtarg.*tauytarg.*tauy + tauytarg.*tauytarg.*taux) +... + thirdK(:, :, 6).*(tauytarg.*tauytarg.*tauy)) + 0.25*hilb); % hilbert subtraction + +end + +% kernels in K21 that are coupled with Hilbert transforms. +% (i.e beta* G_{nx nx nx tauy} + (2-nu)*beta*G_{taux taux tauy nx}) +if strcmpi(type, 'free plate coupled hilbert') + srcnorm = srcinfo.n; + srctang = srcinfo.d; + targnorm = targinfo.n; + targtang = targinfo.d; + coefs = varargin{1}; + + [~, ~, forth] = chnk.helm2d.thirdforth_derivatives(zk, src, targ); % Hankel part + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + zkimag = (1i)*zk; + [~, ~, forthK] = chnk.helm2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + + + + + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + dx1 = repmat((targtang(1,:)).',1,ns); + dy1 = repmat((targtang(2,:)).',1,ns); + + + ds = sqrt(dx.*dx+dy.*dy); + ds1 = sqrt(dx1.*dx1+dy1.*dy1); + + taux = dx./ds; % normalization + tauy = dy./ds; + + tauxtarg = dx1./ds1; + tauytarg = dy1./ds1; + + + + submat = ((1+ coefs(1))/2).*(1/(2*zk^2).*(forth(:, :, 1).*(nxtarg.*nxtarg.*nxtarg.*taux) + ... + forth(:, :, 2).*(nxtarg.*nxtarg.*nxtarg.*tauy + 3*nxtarg.*nxtarg.*nytarg.*taux) + ... + forth(:, :, 4).*(3*nxtarg.*nxtarg.*nytarg.*tauy + 3*nxtarg.*nytarg.*nytarg.*taux) +... + forth(:, :, 6).*(3*nxtarg.*nytarg.*nytarg.*tauy + nytarg.*nytarg.*nytarg.*taux) +... + forth(:, :, 8).*(nytarg.*nytarg.*nytarg.*tauy)) - ... + 1/(2*zk^2).*(forthK(:, :, 1).*(nxtarg.*nxtarg.*nxtarg.*taux) + ... + forthK(:, :, 2).*(nxtarg.*nxtarg.*nxtarg.*tauy + 3*nxtarg.*nxtarg.*nytarg.*taux) + ... + forthK(:, :, 4).*(3*nxtarg.*nxtarg.*nytarg.*tauy + 3*nxtarg.*nytarg.*nytarg.*taux) +... + forthK(:, :, 6).*(3*nxtarg.*nytarg.*nytarg.*tauy + nytarg.*nytarg.*nytarg.*taux) +... + forthK(:, :, 8).*(nytarg.*nytarg.*nytarg.*tauy))) +... + ((2-coefs(1))/2)*(1+coefs(1)).*(1/(2*zk^2).*(forth(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg.*taux) + ... + forth(:, :, 2).*(tauxtarg.*tauxtarg.*nxtarg.*tauy + tauxtarg.*tauxtarg.*nytarg.*taux + 2*tauxtarg.*tauytarg.*nxtarg.*taux) + ... + forth(:, :, 4).*(tauxtarg.*tauxtarg.*nytarg.*tauy + 2*tauxtarg.*tauytarg.*nxtarg.*tauy + tauytarg.*tauytarg.*nxtarg.*taux + 2*tauxtarg.*tauytarg.*nytarg.*taux) +... + forth(:, :, 6).*(tauytarg.*tauytarg.*nxtarg.*tauy + 2*tauxtarg.*tauytarg.*nytarg.*tauy + tauytarg.*tauytarg.*nytarg.*taux) +... + forth(:, :, 8).*(tauytarg.*tauytarg.*nytarg.*tauy)) - ... + 1/(2*zk^2).*(forthK(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg.*taux) + ... + forthK(:, :, 2).*(tauxtarg.*tauxtarg.*nxtarg.*tauy + tauxtarg.*tauxtarg.*nytarg.*taux + 2*tauxtarg.*tauytarg.*nxtarg.*taux) + ... + forthK(:, :, 4).*(tauxtarg.*tauxtarg.*nytarg.*tauy + 2*tauxtarg.*tauytarg.*nxtarg.*tauy + tauytarg.*tauytarg.*nxtarg.*taux + 2*tauxtarg.*tauytarg.*nytarg.*taux) +... + forthK(:, :, 6).*(tauytarg.*tauytarg.*nxtarg.*tauy + 2*tauxtarg.*tauytarg.*nytarg.*tauy + tauytarg.*tauytarg.*nytarg.*taux) +... + forthK(:, :, 8).*(tauytarg.*tauytarg.*nytarg.*tauy))); + +end + +% Updated part in K21 that is not coupled with Hilbert transform. +%(i.e. (1-nu)*(-G_{nx nx ny} + G_{taux taux ny}). + +if strcmpi(type, 'free plate K21 second part') + srcnorm = srcinfo.n; + srctang = srcinfo.d; + targnorm = targinfo.n; + targtang = targinfo.d; + coefs = varargin{1}; + + + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + dx1 = repmat((targtang(1,:)).',1,ns); + dy1 = repmat((targtang(2,:)).',1,ns); + + + ds = sqrt(dx.*dx+dy.*dy); + ds1 = sqrt(dx1.*dx1+dy1.*dy1); + + taux = dx./ds; % normalization + tauy = dy./ds; + + tauxtarg = dx1./ds1; + tauytarg = dy1./ds1; + + zkimag = (1i)*zk; + + [~, third, ~] = chnk.helm2d.thirdforth_derivatives(zk, src, targ); % Hankel part + [~, thirdK, ~] = chnk.helm2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + + submat = -((1-coefs(1))/(2*zk^2).*(third(:, :, 1).*(nxtarg.*nxtarg.*nx) + third(:, :, 2).*(nxtarg.*nxtarg.*ny + 2*nxtarg.*nytarg.*nx) +... + third(:, :, 4).*(2*nxtarg.*nytarg.*ny + nytarg.*nytarg.*nx) +... + third(:, :, 6).*(nytarg.*nytarg.*ny)) - ... + (1-coefs(1))/(2*zk^2).*(thirdK(:, :, 1).*(nxtarg.*nxtarg.*nx) + thirdK(:, :, 2).*(nxtarg.*nxtarg.*ny + 2*nxtarg.*nytarg.*nx) +... + thirdK(:, :, 4).*(2*nxtarg.*nytarg.*ny + nytarg.*nytarg.*nx) +... + thirdK(:, :, 6).*(nytarg.*nytarg.*ny))) + ... + (1-coefs(1))/(2*zk^2).*(third(:, :, 1).*(tauxtarg.*tauxtarg.*nx) + third(:, :, 2).*(tauxtarg.*tauxtarg.*ny + 2*tauxtarg.*tauytarg.*nx) +... + third(:, :, 4).*(2*tauxtarg.*tauytarg.*ny + tauytarg.*tauytarg.*nx) +... + third(:, :, 6).*(tauytarg.*tauytarg.*ny)) - ... + (1-coefs(1))/(2*zk^2).*(thirdK(:, :, 1).*(tauxtarg.*tauxtarg.*nx) + thirdK(:, :, 2).*(tauxtarg.*tauxtarg.*ny + 2*tauxtarg.*tauytarg.*nx) +... + thirdK(:, :, 4).*(2*tauxtarg.*tauytarg.*ny + tauytarg.*tauytarg.*nx) +... + thirdK(:, :, 6).*(tauytarg.*tauytarg.*ny)); % (1-nu)*(- G_{nx nx ny} + G_{ny taux taux}) +end + +% Updated part in K21 that is coupled with Hilbert transform. +%(i.e. (1-nu)*(beta*(G_{taux taux tauy} + 1/4 H) - beta*(G_{nx nx tauy} + 1/4 H)). + +if strcmpi(type, 'free plate K21 hilbert part') + srcnorm = srcinfo.n; + srctang = srcinfo.d; + targnorm = targinfo.n; + targtang = targinfo.d; + coefs = varargin{1}; + + [~, third, ~] = chnk.helm2d.thirdforth_derivatives(zk, src, targ); % Hankel part + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + zkimag = (1i)*zk; + [~, thirdK, ~] = chnk.helm2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + [~,grad] = chnk.lap2d.green(src,targ,true); + + + + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + dx1 = repmat((targtang(1,:)).',1,ns); + dy1 = repmat((targtang(2,:)).',1,ns); + + + ds = sqrt(dx.*dx+dy.*dy); + ds1 = sqrt(dx1.*dx1+dy1.*dy1); + + taux = dx./ds; % normalization + tauy = dy./ds; + + tauxtarg = dx1./ds1; + tauytarg = dy1./ds1; + + + hilb = 2*(grad(:,:,1).*ny - grad(:,:,2).*nx); + + submat = (1-coefs(1)).*(((1+ coefs(1))/2).*(1./(2*zk^2).*(third(:, :, 1).*(tauxtarg.*tauxtarg.*taux) + third(:, :, 2).*(tauxtarg.*tauxtarg.*tauy + 2*tauxtarg.*tauytarg.*taux) +... + third(:, :, 4).*(2*tauxtarg.*tauytarg.*tauy + tauytarg.*tauytarg.*taux) +... + third(:, :, 6).*(tauytarg.*tauytarg.*tauy)) - ... + 1./(2*zk^2).*(thirdK(:, :, 1).*(tauxtarg.*tauxtarg.*taux) + thirdK(:, :, 2).*(tauxtarg.*tauxtarg.*tauy + 2*tauxtarg.*tauytarg.*taux) +... + thirdK(:, :, 4).*(2*tauxtarg.*tauytarg.*tauy + tauytarg.*tauytarg.*taux) +... + thirdK(:, :, 6).*(tauytarg.*tauytarg.*tauy)) + 0.25*hilb) - ... + ((1+ coefs(1))/2)*(1./(2*zk^2).*(third(:, :, 1).*(nxtarg.*nxtarg.*taux) + third(:, :, 2).*(nxtarg.*nxtarg.*tauy+ 2*nxtarg.*nytarg.*taux) +... + third(:, :, 4).*(2*nxtarg.*nytarg.*tauy +nytarg.*nytarg.*taux) +... + third(:, :, 6).*(nytarg.*nytarg.*tauy)) - ... + 1./(2*zk^2).*(thirdK(:, :, 1).*(nxtarg.*nxtarg.*taux) + thirdK(:, :, 2).*(nxtarg.*nxtarg.*tauy+ 2*nxtarg.*nytarg.*taux) +... + thirdK(:, :, 4).*(2*nxtarg.*nytarg.*tauy +nytarg.*nytarg.*taux) +... + thirdK(:, :, 6).*(nytarg.*nytarg.*tauy)) + 0.25*hilb)) ; + + + +end + +% Updated part in K22. (i.e. (1-nu)*(G_{taux taux} - G_{nx nx}) + +if strcmpi(type, 'free plate K22 second part') + srctang = srcinfo.d; + targnorm = targinfo.n; + targtang = targinfo.d; + coefs = varargin{1}; + + + + [hess, ~, ~] = chnk.helm2d.thirdforth_derivatives(zk, src, targ); % Hankel part + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + zkimag = (1i)*zk; + [hessK, ~, ~] = chnk.helm2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + + + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + dx1 = repmat((targtang(1,:)).',1,ns); + dy1 = repmat((targtang(2,:)).',1,ns); + + + ds = sqrt(dx.*dx+dy.*dy); + ds1 = sqrt(dx1.*dx1+dy1.*dy1); + + + tauxtarg = dx1./ds1; + tauytarg = dy1./ds1; + + + submat = (1-coefs(1)).*(1/(2*zk^2).*(hess(:, :, 1).*tauxtarg.*tauxtarg + hess(:, :, 2).*(2*tauxtarg.*tauytarg) + hess(:, :, 3).*tauytarg.*tauytarg)-... + 1/(2*zk^2).*(hessK(:, :, 1).*tauxtarg.*tauxtarg + hessK(:, :, 2).*(2*tauxtarg.*tauytarg) + ... + hessK(:, :, 3).*tauytarg.*tauytarg)-... + (1/(2*zk^2).*(hess(:, :, 1).*nxtarg.*nxtarg + hess(:, :, 2).*(2*nxtarg.*nytarg) + hess(:, :, 3).*nytarg.*nytarg)-... + 1/(2*zk^2).*(hessK(:, :, 1).*nxtarg.*nxtarg + hessK(:, :, 2).*(2*nxtarg.*nytarg) + hessK(:, :, 3).*nytarg.*nytarg))); + + +end + + + + + + + +if strcmpi(type, 'test kernel') % test G_{taux taux tauy nx} + srcnorm = srcinfo.n; + srctang = srcinfo.d; + targnorm = targinfo.n; + targtang = targinfo.d; + coefs = varargin{1}; + + + [hess, third, forth] = chnk.helm2d.thirdforth_derivatives(zk, src, targ); % Hankel part + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + zkimag = (1i)*zk; + [hessK, thirdK, forthK] = chnk.helm2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + [~,grad] = chnk.lap2d.green(src,targ,true); + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + dx1 = repmat((targtang(1,:)).',1,ns); + dy1 = repmat((targtang(2,:)).',1,ns); + + + ds = sqrt(dx.*dx+dy.*dy); + ds1 = sqrt(dx1.*dx1+dy1.*dy1); + + taux = dx./ds; % normalization + tauy = dy./ds; + + tauxtarg = dx1./ds1; + tauytarg = dy1./ds1; + + + hilb = 2*(grad(:,:,1).*ny - grad(:,:,2).*nx); + + submat = 1/(2*zk^2).*(forth(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg.*taux) + ... + forth(:, :, 2).*(tauxtarg.*tauxtarg.*nxtarg.*tauy + tauxtarg.*tauxtarg.*nytarg.*taux + 2*tauxtarg.*tauytarg.*nxtarg.*taux) + ... + forth(:, :, 4).*(tauxtarg.*tauxtarg.*nytarg.*tauy + 2*tauxtarg.*tauytarg.*nxtarg.*tauy + tauytarg.*tauytarg.*nxtarg.*taux + 2*tauxtarg.*tauytarg.*nytarg.*taux) +... + forth(:, :, 6).*(tauytarg.*tauytarg.*nxtarg.*tauy + 2*tauxtarg.*tauytarg.*nytarg.*tauy + tauytarg.*tauytarg.*nytarg.*taux) +... + forth(:, :, 8).*(tauytarg.*tauytarg.*nytarg.*tauy)) - ... + 1/(2*zk^2).*(forthK(:, :, 1).*(tauxtarg.*tauxtarg.*nxtarg.*taux) + ... + forthK(:, :, 2).*(tauxtarg.*tauxtarg.*nxtarg.*tauy + tauxtarg.*tauxtarg.*nytarg.*taux + 2*tauxtarg.*tauytarg.*nxtarg.*taux) + ... + forthK(:, :, 4).*(tauxtarg.*tauxtarg.*nytarg.*tauy + 2*tauxtarg.*tauytarg.*nxtarg.*tauy + tauytarg.*tauytarg.*nxtarg.*taux + 2*tauxtarg.*tauytarg.*nytarg.*taux) +... + forthK(:, :, 6).*(tauytarg.*tauytarg.*nxtarg.*tauy + 2*tauxtarg.*tauytarg.*nytarg.*tauy + tauytarg.*tauytarg.*nytarg.*taux) +... + forthK(:, :, 8).*(tauytarg.*tauytarg.*nytarg.*tauy)); +end + + + + + +% first-kernel of evaluation for the clamped plate +if strcmpi(type, 'first kernel') + 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; + + + + [~, third, ~] = chnk.helm2d.thirdforth_derivatives(zk, src, targ); % Hankel part + + zkimag = 1i*zk; + [~, thirdK, ~] = chnk.helm2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + + + + submat = 1/(2*zk^2).*(third(:, :, 1).*(nx.*nx.*nx) + third(:, :, 2).*(3*nx.*nx.*ny) +... + third(:, :, 4).*(3*nx.*ny.*ny) + third(:, :, 6).*(ny.*ny.*ny)) - ... + 1/(2*zk^2).*(thirdK(:, :, 1).*(nx.*nx.*nx) + thirdK(:, :, 2).*(3*nx.*nx.*ny) +... + thirdK(:, :, 4).*(3*nx.*ny.*ny) + thirdK(:, :, 6).*(ny.*ny.*ny)) + ... + (3/(2*zk^2).*(third(:, :, 1).*(nx.*taux.*taux) + third(:, :, 2).*(2*nx.*taux.*tauy + ny.*taux.*taux) +... + third(:, :, 4).*(nx.*tauy.*tauy + 2*ny.*taux.*tauy) + third(:, :, 6).*(ny.*tauy.*tauy))-... + 3/(2*zk^2).*(thirdK(:, :, 1).*(nx.*taux.*taux) + thirdK(:, :, 2).*(2*nx.*taux.*tauy + ny.*taux.*taux) +... + thirdK(:, :, 4).*(nx.*tauy.*tauy + 2*ny.*taux.*tauy) + thirdK(:, :, 6).*(ny.*tauy.*tauy))); + +end + +% second-kernel of evaluation for the clamped plate +if strcmpi(type, 'second kernel') + 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, ~, ~] = chnk.helm2d.thirdforth_derivatives(zk, src, targ); % Hankel part + + zkimag = 1i*zk; + [hessK, ~, ~] = chnk.helm2d.thirdforth_derivatives(zkimag, src, targ); % modified bessel K part + + submat = -(1/(2*zk^2).*(hess(:, :, 1).*(nx.*nx) + hess(:, :, 2).*(2*nx.*ny) + hess(:, :, 3).*(ny.*ny))-... + 1/(2*zk^2).*(hessK(:, :, 1).*(nx.*nx) + hessK(:, :, 2).*(2*nx.*ny) + hessK(:, :, 3).*(ny.*ny)))+... + (1/(2*zk^2).*(hess(:, :, 1).*(taux.*taux) + hess(:, :, 2).*(2*taux.*tauy) + hess(:, :, 3).*(tauy.*tauy))-... + 1/(2*zk^2).*(hessK(:, :, 1).*(taux.*taux) + hessK(:, :, 2).*(2*taux.*tauy) + hessK(:, :, 3).*(tauy.*tauy))); % -G_{ny ny} + G_{tauy tauy} +end + + +if strcmpi(type, 'free plate eval first') % G_{ny} + srcnorm = srcinfo.n; + [~,grad] = chnk.helm2d.green(zk,src,targ); % Hankel part + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + + + + zkimag = (1i)*zk; + [~,gradK] = chnk.helm2d.green(zkimag,src,targ); % modified bessel K part + + + + + submat = (-1/(2*zk^2).*(grad(:, :, 1).*(nx) + grad(:, :, 2).*ny) + ... + 1/(2*zk^2).*(gradK(:, :, 1).*(nx) + gradK(:, :, 2).*ny)); + + +end + +if strcmpi(type, 'free plate eval first hilbert') % G_{tauy} + srctang = srcinfo.d; + coefs = varargin{1}; + + [~,grad] = chnk.helm2d.green(zk,src,targ); % Hankel part + + + zkimag = (1i)*zk; + [~,gradK] = chnk.helm2d.green(zkimag,src,targ); % modified bessel K part + + dx = repmat(srctang(1,:),nt,1); + dy = repmat(srctang(2,:),nt,1); + + ds = sqrt(dx.*dx+dy.*dy); + + + taux = dx./ds; % normalization + tauy = dy./ds; + + + + submat = ((1 + coefs(1))/2).*(-1/(2*zk^2).*(grad(:, :, 1).*(taux) + grad(:, :, 2).*tauy) + ... + 1/(2*zk^2).*(gradK(:, :, 1).*(taux) + gradK(:, :, 2).*tauy)); % G_{tauy} +end + + +if strcmpi(type, 'free plate eval second') % G = 1/(2k^2) (i/4 H_0^{1} - 1/2pi K_0) + + [val,~] = chnk.helm2d.green(zk,src,targ); % Hankel part + + zkimag = (1i)*zk; + [valK,~] = chnk.helm2d.green(zkimag,src,targ); % modified bessel K part + + submat = 1/(2*zk^2).*val - 1/(2*zk^2).*valK; + +end + + + + + + +% Dirichlet data/potential correpsonding to transmission rep +if strcmpi(type,'trans_rep') + + coef = ones(2,1); + if(nargin == 5); coef = varargin{1}; end; + srcnorm = srcinfo.n; + [submats,grad] = chnk.helm2d.green(zk,src,targ); + nx = repmat(srcnorm(1,:),nt,1); + ny = repmat(srcnorm(2,:),nt,1); + submatd = -(grad(:,:,1).*nx + grad(:,:,2).*ny); + + submat = zeros(1*nt,2*ns); + submat(1:1:1*nt,1:2:2*ns) = coef(1)*submatd; + submat(1:1:1*nt,2:2:2*ns) = coef(2)*submats; +end + +% Neumann data corresponding to transmission rep +if strcmpi(type,'trans_rep_prime') + coef = ones(2,1); + if(nargin == 5); coef = varargin{1}; end; + targnorm = targinfo.n; + srcnorm = srcinfo.n; + + submat = zeros(nt,ns); + + % Get gradient and hessian info + [submats,grad,hess] = chnk.helm2d.green(zk,src,targ); + + nxsrc = repmat(srcnorm(1,:),nt,1); + nysrc = repmat(srcnorm(2,:),nt,1); + nxtarg = repmat((targnorm(1,:)).',1,ns); + nytarg = repmat((targnorm(2,:)).',1,ns); + + % D' + submatdp = -(hess(:,:,1).*nxsrc.*nxtarg ... + + hess(:,:,2).*(nysrc.*nxtarg+nxsrc.*nytarg)... + + hess(:,:,3).*nysrc.*nytarg); + % S' + submatsp = (grad(:,:,1).*nxtarg + grad(:,:,2).*nytarg); + submat = zeros(nt,2*ns) + submat(1:1:1*nt,1:2:2*ns) = coef(1)*submatdp; + submat(1:1:1*nt,2:2:2*ns) = coef(2)*submatsp; +end + + +% Gradient correpsonding to transmission rep +if strcmpi(type,'trans_rep_grad') + coef = ones(2,1); + if(nargin == 5); coef = varargin{1}; end; + + srcnorm = srcinfo.n; + + submat = zeros(nt,ns,6); + % S + [submats,grad,hess] = chnk.helm2d.green(zk,src,targ); + + nxsrc = repmat(srcnorm(1,:),nt,1); + nysrc = repmat(srcnorm(2,:),nt,1); + % D + submatd = -(grad(:,:,1).*nxsrc + grad(:,:,2).*nysrc); + + submat = zeros(2*nt,2*ns); + + submat(1:2:2*nt,1:2:2*ns) = -coef(1)*(hess(:,:,1).*nxsrc + hess(:,:,2).*nysrc); + submat(1:2:2*nt,2:2:2*ns) = coef(2)*grad(:,:,1); + + submat(2:2:2*nt,1:2:2*ns) = -coef(1)*(hess(:,:,2).*nxsrc + hess(:,:,3).*nysrc); + submat(2:2:2*nt,2:2:2*ns) = coef(2)*grad(:,:,2); + + + + + + + + +end + + diff --git a/chunkie/+chnk/+helm2d/smooth_part_derivatives.m b/chunkie/+chnk/+helm2d/smooth_part_derivatives.m new file mode 100644 index 0000000..a041728 --- /dev/null +++ b/chunkie/+chnk/+helm2d/smooth_part_derivatives.m @@ -0,0 +1,82 @@ +function [firstpart, secondpart] = smooth_part_derivatives(k, rx, ry) +%CHNK.HELM2D.SMOOTH_PART_DERIVATIVES evaluate the derivatives of the smooth +% part of asymptotics of the modified biharmonic's green's function +% firstpart are derivatives of: +% 1/(2k^2)*((-1/(2*pi))*(log(k/2) + C)* \sum_{k = 0}^{5} +% ((-1)^k/(k!)^2)*(k*r/2)^{2k} + (1/(2*pi))*(log(k/2) + C)* \sum_{k = 0}^{5} +% (1/(k!)^2)*(k*r/2)^{2k}). (r = |x-y|) +% +% secondpart are derivatives of: +% 1/(2k^2)*(1/(2*pi)*sum_{k =1}^{5} a_{k} ((-1)^k/(k!)^2)*(k*r/2)^{2k} - +% 1/(2*pi)*sum_{k =1}^{5} a_{k} (1/(k!)^2)*(k*r/2)^{2k}) +% with a_k = sum_{m =1}^{k} (1/m) and r = |x-y|. +% +% Derivatives list in the order of: +% f_{x1x1x1y1}, f_{x1x1x1y2}, f_{x1x1x2y2}, f_{x1x2x2y2}, f_{x2x2x2y2} + +rx2 = rx.*rx; +ry2 = ry.*ry; + +r2 = rx2+ry2; + +r = sqrt(r2); + +[m,n] = size(rx); + +if nargout > 0 + firstpart = zeros(m, n, 5); + + firstpart(:, :, 1) = -(double(eulergamma) + log(k/2)).*(k^4).*(5*rx.^2 + ry.^2)./(64*pi); + + firstpart(:, :, 2) = -(double(eulergamma) + log(k/2)).*(k^4).*(rx.*ry)./(32*pi); + + firstpart(:, :, 3) = -(double(eulergamma) + log(k/2)).*(k^4).*(rx.^2 + ry.^2)./(64*pi); + + firstpart(:, :, 4) = -(double(eulergamma) + log(k/2)).*(k^4).*(rx.*ry)./(32*pi); + + firstpart(:, :, 5) = -(double(eulergamma) + log(k/2)).*(k^4).*(rx.^2 + 5*ry.^2)./(64*pi); +end + +if nargout > 1 + secondpart = zeros(m, n, 5); + + secondpart(:, :, 1) = 1/(4*pi*k^2).*(-(9*k^4)/16 - (11/13824).*(k^6).*(-288.*rx.^2 -72*r2) + ... + (25/1769472)*(k^8).*(-384*rx.^4 - 1152*(rx.^2).*r2 - 144*r2.^2) - ... + (137/884736000)*(k^10).*(-1920*(rx.^4).*r2 - 2880*(rx.^2).*(r2.^2) - 240*(r2.^3))) -... + 1/(4*pi*k^2).*(-(9*k^4)/16 + (11/13824).*(k^6).*(-288.*rx.^2 -72*r2) + ... + (25/1769472)*(k^8).*(-384*rx.^4 - 1152*(rx.^2).*r2 - 144*r2.^2) + ... + (137/884736000)*(k^10).*(-1920*(rx.^4).*r2 - 2880*(rx.^2).*(r2.^2) - 240*(r2.^3))); + + secondpart(:, :, 2) = 1/(4*pi*k^2).*(-(25/1769472)*(k^8).*(-384.*(rx.^3).*ry - 576.*rx.*r2.*ry) +... + (137/884736000)*(k^10).*(-1920*(rx.^3).*r2.*ry - 1440*rx.*(r2.^2).*ry) -... + (11/96)*(k^6).*rx.*ry) +... + 1/(4*pi*k^2).*(-(25/1769472)*(k^8).*(-384.*(rx.^3).*ry - 576.*rx.*r2.*ry) -... + (137/884736000)*(k^10).*(-1920*(rx.^3).*r2.*ry - 1440*rx.*(r2.^2).*ry) +... + (11/96)*(k^6).*rx.*ry); + + secondpart(:, :, 3) = 1/(4*pi*k^2).*(-(25/221184)*(k^8).*rx.*(-72*r2.*ry - 48*(ry.^3)) + ... + (137/884736000)*(k^10).*(rx.*(-144*(r2.^2).*ry - 192*r2.*(ry.^3)))-... + (11/96)*(k^6).*rx.*ry) + ... + 1/(4*pi*k^2).*((25/221184)*(k^8).*rx.*(-72*r2.*ry - 48*(ry.^3)) - ... + (137/884736000)*(k^10).*(rx.*(-144*(r2.^2).*ry - 192*r2.*(ry.^3))) +... + (11/96)*(k^6).*rx.*ry) ; + + secondpart(:, :, 4) = 1/(4*pi*k^2).*(-(25/221184)*(k^8).*rx.*(-72*r2.*ry - 48*(ry.^3)) + ... + (137/884736000)*(k^10).*(rx.*(-144*(r2.^2).*ry - 192*r2.*(ry.^3)))-... + (11/96)*(k^6).*rx.*ry) + ... + 1/(4*pi*k^2).*((25/221184)*(k^8).*rx.*(-72*r2.*ry - 48*(ry.^3)) - ... + (137/884736000)*(k^10).*(rx.*(-144*(r2.^2).*ry - 192*r2.*(ry.^3))) +... + (11/96)*(k^6).*rx.*ry) ; + + secondpart(:, :, 5) = 1/(4*pi*k^2).*(-(9*k^4)/16 - (11/13824).*(k^6).*(-288.*ry.^2 -72*r2) + ... + (25/1769472)*(k^8).*(-384*ry.^4 - 1152*(ry.^2).*r2 - 144*r2.^2) - ... + (137/884736000)*(k^10).*(-1920*(ry.^4).*r2 - 2880*(ry.^2).*(r2.^2) - 240*(r2.^3))) -... + 1/(4*pi*k^2).*(-(9*k^4)/16 + (11/13824).*(k^6).*(-288.*ry.^2 -72*r2) + ... + (25/1769472)*(k^8).*(-384*ry.^4 - 1152*(ry.^2).*r2 - 144*r2.^2) + ... + (137/884736000)*(k^10).*(-1920*(ry.^4).*r2 - 2880*(ry.^2).*(r2.^2) - 240*(r2.^3))); + +end + + + +end \ No newline at end of file diff --git a/chunkie/+chnk/+helm2d/thirdforth_derivatives.m b/chunkie/+chnk/+helm2d/thirdforth_derivatives.m new file mode 100644 index 0000000..3342611 --- /dev/null +++ b/chunkie/+chnk/+helm2d/thirdforth_derivatives.m @@ -0,0 +1,178 @@ +function [hess, third, forth] = thirdforth_derivatives(k,src,targ) +%CHNK.HELM2D.GREEN evaluate the second, third, and forth derivatives of +% helmholtz green's function: i/4 H_0^1(k|x-y|). +% for the given sources and targets + +% the third derivatives are: (listed in order) +% G_{x1x1y1}, G_{x1x1y2}, G_{x1x2y1}, G_{x1x2y2}, G_{x2x2y1}, G_{x2x2y2} +% +% forth derivatives are: +%G_{x1x1x1y1}, G_{x1x1x1y2}, G_{x1x1x2y1}, G_{x1x1x2y2}, G_{x1x2x2y1}, +%G_{x1x2x2y2}, G_{x2x2x2y1}, G_{x2x2x2y2} + +[~,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); + +rx = xt-xs; +ry = yt-ys; + +rx2 = rx.*rx; +ry2 = ry.*ry; + +r2 = rx2+ry2; + +r = sqrt(r2); + +[m,n] = size(xs); + +if nargout > 0 + h0 = besselh(0,1,k*r); + hess = zeros(m,n,3); + h1 = besselh(1,1,k*r); + r3 = r.^3; + + h2 = 2*h1./(k*r)-h0; + + hess(:,:,1) = 0.25*1i*k*((rx-ry).*(rx+ry).*h1./r3 - k*rx2.*h0./r2); + hess(:,:,2) = 0.25*1i*k*k*rx.*ry.*h2./r2; + hess(:,:,3) = 0.25*1i*k*((ry-rx).*(rx+ry).*h1./r3 - k*ry2.*h0./r2); +end + + +if nargout > 1 + h3 = besselh(0,1,k*r) - besselh(2, 1, k*r); + h4 = h1 - besselh(3, 1, k*r); + + third = zeros(m,n,6); + + third(:,:,1) = 3i*k*0.25*h1.*(rx.^3)./(r.^5) - 3i*k*0.25*rx.*h1./(r.^3) - ... + 3i*(k^2)*0.25*(rx.^3).*h3./(2*r.^4) + 3i*(k^2)*0.25.*rx.*h3./(2*r.^2) - 1i*(k^3)*0.25*(rx.^3).*h1./(2*r.^3) -... + 1i*(k^3)*0.25.*(rx.^3).*h4./(4*r.^3); + + third(:,:,2) = 3i*k*0.25*(rx.^2).*ry.*h1./(r.^5) - 1i*0.25*k.*ry.*h1./(r.^3) - ... + (1i)*(0.25)*3*(k^2)*(rx.^2).*(ry).*h3./(2*r.^4) + 1i*(k^2)*0.25*ry.*h3./(2*r.^2) - ... + 1i*(k^3)*0.25*(rx.^2).*ry.*h1./(2*r.^3) - 1i*(k^3)*0.25*(rx.^2).*ry.*h4./(4*r.^3); + + third(:,:, 3) = 3i*k*0.25*(rx.^2).*ry.*h1./(r.^5) - 1i*0.25*k.*ry.*h1./(r.^3) - ... + (1i)*(0.25)*3*(k^2)*(rx.^2).*(ry).*h3./(2*r.^4) + 1i*(k^2)*0.25*ry.*h3./(2*r.^2) - ... + 1i*(k^3)*0.25*(rx.^2).*ry.*h1./(2*r.^3) - 1i*(k^3)*0.25*(rx.^2).*ry.*h4./(4*r.^3); + + third(:, :, 4) = -1i*k*0.25*rx.*h1./(r.^3) + 3i*k*0.25.*rx.*(ry.^2).*h1./(r.^5) + ... + 1i*(k^2).*rx.*h3./(8*r.^2) - 3i*(k^2).*rx.*(ry.^2).*h3./(8*r.^4) - ... + 1i*(k^3).*(rx).*(ry.^2).*h1./(8*r.^3) - 1i*(k^3)*rx.*(ry.^2).*h4./(16*r.^3); + + third(:, :, 5) = -1i*k*0.25.*rx.*h1./(r.^3) + 3i*k*0.25.*rx.*(ry.^2).*h1./(r.^5) + ... + 1i*(k^2).*rx.*h3./(8*r.^2) - 3i*(k^2).*rx.*(ry.^2).*h3./(8*r.^4) - ... + 1i*(k^3).*(rx).*(ry.^2).*h1./(8*r.^3) - 1i*(k^3)*rx.*(ry.^2).*h4./(16*r.^3); + + third(:, :, 6) = -3i*k*0.25.*ry.*h1./(r.^3) + 3i*0.25*k*(ry.^3).*h1./(r.^5) + ... + 3i*0.25*(k^2)*ry.*h3./(2*r.^2) - 3i*(k^2)*(0.25)*(ry.^3).*h3./(2*r.^4) - ... + 1i*(0.25)*(k^3).*(ry.^3).*h1./(2*r.^3) - 1i*(k^3)*0.25.*(ry.^3).*h4./(4*r.^3); +end + +if nargout > 2 + h5 = besselh(2, 1, k*r) - besselh(4, 1, k*r); + + forth = zeros(m,n,8); + + forth(:, :, 1) = 0.25*(1i)*(-15)*k.*(rx.^4).*h1./(r.^7) + 0.25*(1i)*18*k*(rx.^2).*h1./(r.^5)-... + 3*(1i)*(0.25)*k*h1./(r.^3) + 0.25*(1i)*(15*k^2).*(rx.^4).*h3./(2*r.^6) - ... + (1i)*(0.25)*9*(k^2)*(rx.^2).*h3./(r.^4) + (1i)*(0.25)*3*(k^2).*h3./(2*r.^2)+... + (1i)*(0.25)*(k^3)*(rx.^4).*h1./(r.^5) + (1i)*(0.25)*(k^3)*(rx.^4).*h4./(2*r.^5) + ... + -(1i)*(0.25)*(k^3)*(rx.^2).*h1./(r.^3) - (1i)*(0.25)*(k^3)*(rx.^2).*h4./(2*r.^3) +... + (1i)*(0.25)*(3*k^3)*(rx.^4).*h1./(2*r.^5) + (1i)*(0.25)*(3*k^3)*(rx.^4).*h4./(4*r.^5) -... + (1i)*(0.25)*(3*k^3)*(rx.^2).*h1./(2*r.^3) - (1i)*(0.25)*(3*k^3)*(rx.^2).*h4./(4*r.^3) +... + (1i)*(0.25)*(0.5)*(k^3)*(rx.^4).*h1./(r.^5) - (1i)*(0.25)*(0.5)*(k^3)*(rx.^2).*h1./(r.^3)-... + (1i)*(0.25)*(0.5)*(k^4)*(rx.^4).*h3./(2*r.^4) + (1i)*(0.25)*(0.5)*(k^3)*(rx.^4).*h4./(2*r.^5) -... + (1i)*(0.25)*(0.5)*(k^3)*(rx.^2).*h4./(2*r.^3) + (1i)*0.25*(0.5)*(k^4)*(-rx.^4).*h3./(4*r.^4) +... + (1i)*0.25*(0.5)*(k^4)*(rx.^4).*h5./(4*r.^4); + + forth(:, :, 2) = (0.25)*(1i)*(-15)*k.*(rx.^3).*(ry).*h1./(r.^7) + 0.25*(1i)*9*k*(rx.*ry).*h1./(r.^5) + ... + (0.25)*(1i)*(15)*(k^2).*(rx.^3).*(ry).*(h3)./(2*r.^6) - (0.25)*(1i)*9*(k^2)*(rx.*ry).*h3./(2*r.^4) -... + ((0.25)*(1i)*(k^3)*(-rx.^3).*(ry).*h1./(r.^5) + 0.25*(1i)*(k^3)*(-rx.^3).*ry.*h4./(2*r.^5)) + ... + (0.25)*(1i)*(3*k^3).*(rx.^3).*ry.*h1./(2*r.^5) + (0.25)*(1i)*(3*k^3).*(rx.^3).*ry.*h4./(4*r.^5) -... + (0.25)*(1i)*(3*k^3)*(rx.*ry).*h1./(2*r.^3) -(0.25)*(1i)*(3*k^3).*(rx.*ry).*h4./(4*r.^3) +... + (0.25)*(1i)*(k^3)*(rx.^3).*ry.*h1./(2*r.^5) - (0.25)*(1i)*(k^4).*(rx.^3).*ry.*h3./(4*r.^4) + ... + (0.25)*(1i)*(k^3)*(rx.^3).*ry.*h4./(4*r.^5) + (0.25)*(1i)*(k^4).*(rx.^3).*(-ry).*h3./(8*r.^4) + ... + (0.25)*(1i)*(k^4)*(rx.^3).*ry.*h5./(8*r.^4); + + forth(:, :, 3) = (0.25)*(1i)*(-15)*k.*(rx.^3).*(ry).*h1./(r.^7) + 0.25*(1i)*9*k*(rx.*ry).*h1./(r.^5) + ... + (0.25)*(1i)*(15)*(k^2).*(rx.^3).*(ry).*(h3)./(2*r.^6) - (0.25)*(1i)*9*(k^2)*(rx.*ry).*h3./(2*r.^4) -... + ((0.25)*(1i)*(k^3)*(-rx.^3).*(ry).*h1./(r.^5) + 0.25*(1i)*(k^3)*(-rx.^3).*ry.*h4./(2*r.^5)) + ... + (0.25)*(1i)*(3*k^3).*(rx.^3).*ry.*h1./(2*r.^5) + (0.25)*(1i)*(3*k^3).*(rx.^3).*ry.*h4./(4*r.^5) -... + (0.25)*(1i)*(3*k^3)*(rx.*ry).*h1./(2*r.^3) -(0.25)*(1i)*(3*k^3).*(rx.*ry).*h4./(4*r.^3) +... + (0.25)*(1i)*(k^3)*(rx.^3).*ry.*h1./(2*r.^5) - (0.25)*(1i)*(k^4).*(rx.^3).*ry.*h3./(4*r.^4) + ... + (0.25)*(1i)*(k^3)*(rx.^3).*ry.*h4./(4*r.^5) + (0.25)*(1i)*(k^4).*(rx.^3).*(-ry).*h3./(8*r.^4) + ... + (0.25)*(1i)*(k^4)*(rx.^3).*ry.*h5./(8*r.^4); + + forth(:, :, 4) = (0.25)*(1i)*(3*k)*(rx.^2).*h1./(r.^5) -(0.25)*(1i)*k.*h1./(r.^3) - ... + (0.25)*(1i)*15*k.*(rx.^2).*(ry.^2).*(h1)./(r.^7) + (0.25)*(1i)*(3*k).*(ry.^2).*h1./(r.^5)-... + (0.25)*(1i)*(3*k^2).*(rx.^2).*h3./(2*r.^4) + (0.25)*(1i)*(k^2).*h3./(2*r.^2) +... + (0.25)*(1i)*15*(k^2).*(rx.^2).*(ry.^2).*h3./(2*r.^6) - (0.25)*(1i)*(3*k^2)*(ry.^2).*h3./(2*r.^4) -... + ((0.25)*(1i)*(k^3).*(-rx.^2).*(ry.^2).*h1./(r.^5) + (0.25)*(1i)*(k^3).*(-rx.^2).*(ry.^2).*h4./(2*r.^5)) +... + (0.25)*(1i)*(3*k^3).*(rx.^2).*(ry.^2).*h1./(2*r.^5) + (0.25)*(1i)*(3*k^3).*(rx.^2).*(ry.^2).*h4./(4*r.^5) - ... + (0.25)*(1i)*(k^3).*(ry.^2).*h1./(2*r.^3) - (0.25)*(1i)*(k^3).*(ry.^2).*h4./(4*r.^3) +... + (0.25)*(1i)*(k^3).*(-rx.^2).*h1./(2*r.^3) + (0.25)*(1i)*(k^3).*(rx.^2).*(ry.^2).*h1./(2*r.^5)-... + (0.25)*(1i)*(k^4).*(rx.^2).*(ry.^2).*h3./(4*r.^4) - (0.25)*(1i)*(k^3).*(rx.^2).*h4./(4*r.^3) + ... + (0.25)*(1i)*(k^3).*(rx.^2).*(ry.^2).*h4./(4*r.^5) -(0.25)*(1i)*(k^4)*(rx.^2).*(ry.^2).*h3./(8*r.^4)+... + (0.25)*(1i)*(k^4).*(rx.^2).*(ry.^2).*h5./(8*r.^4); + + forth(:, :, 5) = (0.25)*(1i)*(3*k)*(rx.^2).*h1./(r.^5) -(0.25)*(1i)*k.*h1./(r.^3) - ... + (0.25)*(1i)*15*k.*(rx.^2).*(ry.^2).*(h1)./(r.^7) + (0.25)*(1i)*(3*k).*(ry.^2).*h1./(r.^5)-... + (0.25)*(1i)*(3*k^2).*(rx.^2).*h3./(2*r.^4) + (0.25)*(1i)*(k^2).*h3./(2*r.^2) +... + (0.25)*(1i)*15*(k^2).*(rx.^2).*(ry.^2).*h3./(2*r.^6) - (0.25)*(1i)*(3*k^2)*(ry.^2).*h3./(2*r.^4) -... + ((0.25)*(1i)*(k^3).*(-rx.^2).*(ry.^2).*h1./(r.^5) + (0.25)*(1i)*(k^3).*(-rx.^2).*(ry.^2).*h4./(2*r.^5)) +... + (0.25)*(1i)*(3*k^3).*(rx.^2).*(ry.^2).*h1./(2*r.^5) + (0.25)*(1i)*(3*k^3).*(rx.^2).*(ry.^2).*h4./(4*r.^5) - ... + (0.25)*(1i)*(k^3).*(ry.^2).*h1./(2*r.^3) - (0.25)*(1i)*(k^3).*(ry.^2).*h4./(4*r.^3) +... + (0.25)*(1i)*(k^3).*(-rx.^2).*h1./(2*r.^3) + (0.25)*(1i)*(k^3).*(rx.^2).*(ry.^2).*h1./(2*r.^5)-... + (0.25)*(1i)*(k^4).*(rx.^2).*(ry.^2).*h3./(4*r.^4) - (0.25)*(1i)*(k^3).*(rx.^2).*h4./(4*r.^3) + ... + (0.25)*(1i)*(k^3).*(rx.^2).*(ry.^2).*h4./(4*r.^5) -(0.25)*(1i)*(k^4)*(rx.^2).*(ry.^2).*h3./(8*r.^4)+... + (0.25)*(1i)*(k^4).*(rx.^2).*(ry.^2).*h5./(8*r.^4); + + forth(:, :, 6) = (-0.25)*(1i)*k*((-9*rx.*ry./(r.^5) + 15*rx.*(ry.^3)./(r.^7)).*h1 + 2*k*rx.*(ry).*h3./(r.^4)-... + k*(-rx./(r.^3) + 3*rx.*(ry.^2)./(r.^5)).*ry.*h3./(2*r) - 4*k*rx.*(ry.^3).*h3./(r.^6) + ... + k*rx.*(ry).*h3./(2*r.^4) - k*rx.*(ry.^3).*h3./(2*r.^6) - (k^2)*rx.*(ry.^3).*h1./(2*r.^5) - (k^2)*rx.*(ry.^3).*h4./(4*r.^5)-... + (k^2)*rx.*(ry.^3).*h1./(r.^5) - (k^2)*rx.*(ry.^3).*h4./(2*r.^5) + 3*k*rx.*ry.*h3./(2*r.^4) -... + 3*k*rx.*(ry.^3).*h3./(2*r.^6) + (k^2)*rx.*ry.*h1./(2*r.^3) + (k^2)*rx.*ry.*h4./(4*r.^3) - (k^2)*rx.*(ry.^3).*h1./(2*r.^5)-... + (k^2)*rx.*(ry.^3).*h4./(4*r.^5) + (k^2)*rx.*(ry).*h1./(2*r.^3) +(k^2)*rx.*ry.*h4./(4*r.^3) - ... + (k^2)*rx.*(ry.^3).*h1./(2*r.^5) - (k^2)*rx.*(ry.^3).*h4./(4*r.^5) + (k^2)*rx.*ry.*h1./(2*r.^3)-... + (k^2)*rx.*(ry.^3).*h1./(2*r.^5) +(k^3)*rx.*(ry.^3).*h3./(4*r.^4) + (k^2)*rx.*(ry).*h4./(4*r.^3)-... + (k^2)*rx.*(ry.^3).*h4./(4*r.^5) - (k^3)*rx.*(-ry.^3).*h3./(8*r.^4) - (k^3)*rx.*(ry.^3).*h5./(8*r.^4)); + + forth(:, :, 7) = (-0.25)*(1i)*k*((-9*rx.*ry./(r.^5) + 15*rx.*(ry.^3)./(r.^7)).*h1 + 2*k*rx.*(ry).*h3./(r.^4)-... + k*(-rx./(r.^3) + 3*rx.*(ry.^2)./(r.^5)).*ry.*h3./(2*r) - 4*k*rx.*(ry.^3).*h3./(r.^6) + ... + k*rx.*(ry).*h3./(2*r.^4) - k*rx.*(ry.^3).*h3./(2*r.^6) - (k^2)*rx.*(ry.^3).*h1./(2*r.^5) - (k^2)*rx.*(ry.^3).*h4./(4*r.^5)-... + (k^2)*rx.*(ry.^3).*h1./(r.^5) - (k^2)*rx.*(ry.^3).*h4./(2*r.^5) + 3*k*rx.*ry.*h3./(2*r.^4) -... + 3*k*rx.*(ry.^3).*h3./(2*r.^6) + (k^2)*rx.*ry.*h1./(2*r.^3) + (k^2)*rx.*ry.*h4./(4*r.^3) - (k^2)*rx.*(ry.^3).*h1./(2*r.^5)-... + (k^2)*rx.*(ry.^3).*h4./(4*r.^5) + (k^2)*rx.*(ry).*h1./(2*r.^3) +(k^2)*rx.*ry.*h4./(4*r.^3) - ... + (k^2)*rx.*(ry.^3).*h1./(2*r.^5) - (k^2)*rx.*(ry.^3).*h4./(4*r.^5) + (k^2)*rx.*ry.*h1./(2*r.^3)-... + (k^2)*rx.*(ry.^3).*h1./(2*r.^5) +(k^3)*rx.*(ry.^3).*h3./(4*r.^4) + (k^2)*rx.*(ry).*h4./(4*r.^3)-... + (k^2)*rx.*(ry.^3).*h4./(4*r.^5) - (k^3)*rx.*(-ry.^3).*h3./(8*r.^4) - (k^3)*rx.*(ry.^3).*h5./(8*r.^4)); + + forth(:, :, 8) = (-0.25)*(1i)*(3*k)*h1./(r.^3) + 0.25*(1i)*(18*k)*(ry.^2).*h1./(r.^5) - (0.25)*(1i)*(15*k).*(ry.^4).*h1./(r.^7) +... + (0.25)*(1i)*(3*k^2).*(h3)./(2*r.^2) - (0.25)*(1i)*(9*k^2).*(ry.^2).*h3./(r.^4) + (0.25)*(1i)*(15*k^2).*(ry.^4).*h3./(2*r.^6) +... + (0.25)*(1i)*(k^3)*(-ry.^2).*h1./(r.^3) - (0.25)*(1i)*(k^3).*(ry.^2).*h4./(2*r.^3) - (0.25)*(1i)*(k^3).*(-ry.^4).*h1./(r.^5) + ... + (0.25)*(1i)*(k^3)*(ry.^4).*h4./(2*r.^5) - (0.25)*(1i)*(3*k^3).*(ry.^2).*h1./((2*r.^3)) - (0.25)*(1i)*(3*k^3).*(ry.^2).*h4./(4*r.^3) +... + (0.25)*(1i)*(3*k^3)*(ry.^4).*h1./(2*r.^5) + (0.25)*(1i)*(3*k^3).*(ry.^4).*h4./(4*r.^5) - (0.25)*(1i)*(k^3).*(ry.^2).*h1./(2*r.^3) + ... + (0.25)*(1i)*(k^3)*(ry.^4).*h1./(2*r.^5) - (0.25)*(1i)*(k^4).*(ry.^4).*h3./(4*r.^4) - (0.25)*(1i)*(k^3).*(ry.^2).*h4./(4*r.^3) + ... + (0.25)*(1i)*(k^3)*(ry.^4).*h4./(4*r.^5) + (0.25)*(1i)*(k^4).*(-ry.^4).*h3./(8*r.^4) + (0.25)*(1i)*(k^4)*(ry.^4).*h5./(8*r.^4); + +end +end + + + + + + + + + + + diff --git a/chunkie/+chnk/+lap2d/kern.m b/chunkie/+chnk/+lap2d/kern.m index 3913961..67d9fb8 100644 --- a/chunkie/+chnk/+lap2d/kern.m +++ b/chunkie/+chnk/+lap2d/kern.m @@ -33,6 +33,15 @@ submat = (-grad(:,:,1).*ny + grad(:,:,2).*nx); end +if strcmpi(type,'hilb') + srcnorm = chnk.normal2d(srcinfo); + [~,grad] = chnk.lap2d.green(src,targ,true); + nx = repmat((srcnorm(1,:)),nt,1); + ny = repmat((srcnorm(2,:)),nt,1); + + submat = 2*(grad(:,:,1).*ny - grad(:,:,2).*nx); +end + if strcmpi(type,'sgrad') [~,grad] = chnk.lap2d.green(src,targ,true); submat = reshape(permute(grad,[3,1,2]),2*nt,ns); diff --git a/chunkie/circle.m b/chunkie/circle.m new file mode 100644 index 0000000..4f0c42c --- /dev/null +++ b/chunkie/circle.m @@ -0,0 +1,23 @@ + + + + + +% define the geometry of the circle +function [r, d, d2]= circle(t) + +x = cos(t); +y = sin(t); + +dx = -sin(t); +dy = cos(t); + +dxx = cos(t); +dyy = sin(t); + +r = [(x(:)).';(y(:)).']; +d = [(dx(:)).'; (dy(:)).']; +d2 = [(dxx(:)).'; (dyy(:)).']; + + +end \ No newline at end of file diff --git a/chunkie/demo/flexural/test_free_plate_circle.m b/chunkie/demo/flexural/test_free_plate_circle.m new file mode 100644 index 0000000..cfd3b33 --- /dev/null +++ b/chunkie/demo/flexural/test_free_plate_circle.m @@ -0,0 +1,201 @@ +clear +clc + +kvec = [1;-1.5]; + +% + +zk = sqrt(norm(kvec)); % our k (wave number) +nu = 1/3; +cparams = []; + +cparams.eps = 1e-6; +%cparams.nover = 0; +cparams.maxchunklen = 4./zk; % setting a chunk length helps when the + % frequency is known' + + + +chnkr = chunkerfunc(@(t) starfish(t), cparams); + +figure(1) % plot the chunker-object (supposed to be a circle centered at 1 with radius 1) +clf +plot(chnkr, '-x') +title('Chunkr object') +hold on +quiver(chnkr) +axis equal +drawnow() + +coefs = [nu; 0]; +opts = []; +opts.sing = 'log'; +fkern = @(s,t) chnk.helm2d.kern(zk, s, t, 'free plate first part', coefs); % build the desired kernel + + +fkern1 = @(s,t) chnk.helm2d.kern(zk, s, t, 'free plate hilbert subtract', coefs); % hilbert subtraction kernels in K11 +fkern2 = @(s,t) chnk.helm2d.kern(zk, s, t, 'free plate coupled hilbert', coefs); +hilbert = @(s,t) chnk.lap2d.kern(s, t, 'hilb'); +double = @(s,t) chnk.lap2d.kern(s,t, 'd'); + +fkern3 = @(s,t) chnk.helm2d.kern(zk, s, t, 'free plate K21 first part', coefs); % singularity subtration kernel in K21 (including swapping its Asmyptotics expansions) + +fkern4 = @(s,t) chnk.helm2d.kern(zk, s, t, 'free plate K21 second part', coefs); % kernels in K21 needs to multiply by curvature + +fkern5 = @(s,t) chnk.helm2d.kern(zk, s, t, 'free plate K21 hilbert part', coefs); % kernels in K21 coupled with hilbert transforms and needs to multiply by curvature + +fkern6 = @(s,t) chnk.helm2d.kern(zk, s, t, 'free plate K22 second part', coefs); % kernels in K22 needs to multiply by curvature + + +start = tic; +sysmat = chunkermat(chnkr,fkern, opts); +sysmat1 = chunkermat(chnkr, fkern1, opts); +sysmat2 = chunkermat(chnkr, fkern2, opts); +K21 = chunkermat(chnkr, fkern3, opts); +K21second = chunkermat(chnkr, fkern4, opts); +K21hilbert = chunkermat(chnkr, fkern5, opts); + +K22second = chunkermat(chnkr, fkern6, opts); + +D = chunkermat(chnkr, double, opts); +t1 = toc(start); +fprintf('%5.2e s : time to assemble matrix\n',t1) + +opts2 = []; +opts2.sing = 'pv'; + +start = tic; +H = chunkermat(chnkr, hilbert, opts2); % Assemble hilbert transforms +t1 = toc(start); +fprintf('%5.2e s : time to assemble matrix\n',t1) + + +hilb = sysmat1*H - ((1+nu)/2).*(D*D)- ((1+nu)*nu/2).*(D*D); +hilb2 = sysmat2*H + K21hilbert*H; + +mat1 = sysmat(1:2:end, 1:2:end); +mat4 = sysmat(2:2:end, 2:2:end); + + + +sysmat(1:2:end, 1:2:end) = mat1 + hilb; +sysmat(2:2:end, 1:2:end) = K21 + hilb2 + K21second; +sysmat(2:2:end, 2:2:end) = mat4 + K22second; + + + + +A = [-1/2 + (1/8)*(1+nu).^2, 0; 0, 1/2]; % jump matrix (for interior problem) + +M = kron(eye(chnkr.npt), A); + +lhs = M + sysmat ; + +[hess, third, ~] = chnk.helm2d.thirdforth_derivatives(zk, [0;0], chnkr.r); +[hessK, thirdK, ~] = chnk.helm2d.thirdforth_derivatives(zk*(1i), [0;0], 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 = 1/(2*zk^2).*(hess(:, :, 1).*(nx.*nx) + hess(:, :, 2).*(2*nx.*ny) + hess(:, :, 3).*(ny.*ny))-... + 1/(2*zk^2).*(hessK(:, :, 1).*(nx.*nx) + hessK(:, :, 2).*(2*nx.*ny) + hessK(:, :, 3).*(ny.*ny))+... + coefs(1)/(2*zk^2).*(hess(:, :, 1).*(taux.*taux) + hess(:, :, 2).*(2*taux.*tauy) + hess(:, :, 3).*(tauy.*tauy))-... + coefs(1)/(2*zk^2).*(hessK(:, :, 1).*(taux.*taux) + hessK(:, :, 2).*(2*taux.*tauy) + ... + hessK(:, :, 3).*(tauy.*tauy)); + +secondbc = -1./(2*zk^2).*(third(:, :, 1).*(nx.*nx.*nx) + third(:, :, 2).*(3*nx.*nx.*ny) +... + third(:, :, 4).*(3*nx.*ny.*ny) + third(:, :, 6).*(ny.*ny.*ny)) + ... + 1./(2*zk^2).*(thirdK(:, :, 1).*(nx.*nx.*nx) + thirdK(:, :, 2).*(3*nx.*nx.*ny)+... + thirdK(:, :, 4).*(3*nx.*ny.*ny) + thirdK(:, :, 6).*(ny.*ny.*ny)) -... + (2-coefs(1))/(2*zk^2).*(third(:, :, 1).*(taux.*taux.*nx) + third(:, :, 2).*(taux.*taux.*ny + 2*taux.*tauy.*nx) +... + third(:, :, 4).*(2*taux.*tauy.*ny+ tauy.*tauy.*nx) +... + + third(:, :, 6).*(tauy.*tauy.*ny)) + ... + (2-coefs(1))/(2*zk^2).*(thirdK(:, :, 1).*(taux.*taux.*nx) + thirdK(:, :, 2).*(taux.*taux.*ny + 2*taux.*tauy.*nx) +... + thirdK(:, :, 4).*(2*taux.*tauy.*ny+ tauy.*tauy.*nx) +... + + thirdK(:, :, 6).*(tauy.*tauy.*ny)) +... + (1-coefs(1)).*(1/(2*zk^2).*(hess(:, :, 1).*taux.*taux + hess(:, :, 2).*(2*taux.*tauy) + hess(:, :, 3).*tauy.*tauy)-... + 1/(2*zk^2).*(hessK(:, :, 1).*taux.*taux + hessK(:, :, 2).*(2*taux.*tauy) + ... + hessK(:, :, 3).*tauy.*tauy)-... + (1/(2*zk^2).*(hess(:, :, 1).*nx.*nx + hess(:, :, 2).*(2*nx.*ny) + hess(:, :, 3).*ny.*ny)-... + 1/(2*zk^2).*(hessK(:, :, 1).*nx.*nx + hessK(:, :, 2).*(2*nx.*ny) + hessK(:, :, 3).*ny.*ny))); + +[nt, ~] = size(sysmat); + +rhs = zeros(nt, 1); +rhs(1:2:end) = firstbc ; +rhs(2:2:end) = secondbc; + +tic +%sol = gmres(lhs, rhs, [], 1e-13, 200); +sol = lhs\rhs; +toc; + +rho1 = sol(1:2:end); % first density +rho2 = sol(2:2:end); + + +xs = -3:0.01:3; % generate some targets +ys = -3:0.01:3; +[X,Y] = meshgrid(xs, ys); +targets = [X(:).'; Y(:).']; +[~,na] = size(targets); + +tic +in = chunkerinterior(chnkr, targets); +out = ~in; +trad = targets(1,:).^2+targets(2,:).^2; +%out = out.*(sqrt(trad.')>2); +out = find(out); +toc + + +ikern1 = @(s,t) chnk.helm2d.kern(zk, s, t, 'free plate eval first', coefs); % build the kernel of evaluation +ikern2 = @(s,t) chnk.helm2d.kern(zk, s, t, 'free plate eval second'); +ikern3 = @(s,t) chnk.helm2d.kern(zk, s, t, 'free plate eval first hilbert',coefs); + +opts_off = []; +opts_off.forcesmooth = true; +opts_off.accel = false; +%coupled = chunkerkernevalmat(chnkr, ikern3, targets(:, out), opts_off); + +start1 = tic; +Hrho = H*rho1; +Dsol1 = chunkerkerneval(chnkr, ikern1,rho1, targets(:, out),opts_off); +Dsol2 = chunkerkerneval(chnkr, ikern3, Hrho, targets(:, out),opts_off); +Dsol3 = chunkerkerneval(chnkr, ikern2, rho2, targets(:,out),opts_off); +Dsol = Dsol1 + Dsol2 + Dsol3; +t2 = toc(start1); +fprintf('%5.2e s : time for kernel eval (for plotting)\n',t2) + +true_sol = zeros(na, 1); +utarg = zeros(na, 1); + +[val,~] = chnk.helm2d.green(zk,[0;0],targets(:,out)); % Hankel part + +zkimag = (1i)*zk; +[valK,~] = chnk.helm2d.green(zkimag,[0;0], targets(:,out)); % modified bessel K part + +trueval = 1/(2*zk^2).*val - 1/(2*zk^2).*valK; + + +utarg(out) = Dsol; +true_sol(out) = trueval; + + +uerr = utarg - true_sol; +uerr = reshape(uerr,size(X)); +figure(2) +h = pcolor(X,Y,log10(abs(uerr))); +set(h,'EdgeColor','None'); hold on; +title("Absolute error (free plate kernel on a circle)", 'FontSize',16) +plot(chnkr,'w-','LineWidth',2); + +colorbar + diff --git a/chunkie/ellipse.m b/chunkie/ellipse.m new file mode 100644 index 0000000..15ad679 --- /dev/null +++ b/chunkie/ellipse.m @@ -0,0 +1,16 @@ +function [r, d, d2] = ellipse(t) + +x = 3.*cos(t); +y = sin(t); + +x1 = -3.*sin(t); +y1 = cos(t); + +x2 = -3.*cos(t); +y2 = -sin(t); + +r = [(x(:)).'; (y(:)).']; +d = [(x1(:)).'; (y1(:)).']; +d2 = [(x2(:)).'; (y2(:)).']; + +end \ No newline at end of file