From 1c26fae95c9b764fb31c6e5be78c59a8d2cdfbd8 Mon Sep 17 00:00:00 2001 From: Manas Rachh Date: Tue, 26 Sep 2023 16:43:34 -0400 Subject: [PATCH 1/2] beta test --- chunkie/+chnk/+helm1d/kern.m | 63 +++++++++++++++++++++++ chunkie/+chnk/+quadggq/getremovablequad.m | 37 +++++++++++++ chunkie/+chnk/+quadggq/setup.m | 7 ++- chunkie/@kernel/helm1d.m | 35 +++++++++++++ chunkie/@kernel/kernel.m | 4 ++ chunkie/chunkermat.m | 13 +++++ 6 files changed, 158 insertions(+), 1 deletion(-) create mode 100644 chunkie/+chnk/+helm1d/kern.m create mode 100644 chunkie/+chnk/+quadggq/getremovablequad.m create mode 100644 chunkie/@kernel/helm1d.m diff --git a/chunkie/+chnk/+helm1d/kern.m b/chunkie/+chnk/+helm1d/kern.m new file mode 100644 index 00000000..1671ba55 --- /dev/null +++ b/chunkie/+chnk/+helm1d/kern.m @@ -0,0 +1,63 @@ + +function submat= kern(zk,srcinfo,targinfo,type,varargin) +%CHNK.HELM1D.KERN standard Helmholtz layer potential kernels in 2D +% +% Syntax: submat = chnk.heml1d.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) = e^{i zk |t-t'|} where t parametrization +% coordinate +% +% D(x,y) = \nabla_{n_y} G(x,y) +% S(x,y) = G(x,y) +% S'(x,y) = \nabla_{n_x} 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 == 'c', combined layer kernel D + i eta S +% varargin{1} - eta in the combined layer formula, 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(1,:); +src = src(:).'; +targ = targinfo.r(1,:); +targ = targ(:); + +[~,ns] = size(src); +[~,nt] = size(targ); + + +if strcmpi(type,'s') + submat = exp(1j*zk*abs(targ-src))/2/1j/zk; +end +end + + diff --git a/chunkie/+chnk/+quadggq/getremovablequad.m b/chunkie/+chnk/+quadggq/getremovablequad.m new file mode 100644 index 00000000..186da344 --- /dev/null +++ b/chunkie/+chnk/+quadggq/getremovablequad.m @@ -0,0 +1,37 @@ +function [xs0, ws0] = getremovablequad(k, nfac) +%CHNK.QUADGGQ.GETREMOVABLEQUAD returns a quadrature +% for integrating piecewise smooth quadrature, where the rule +% on each panel is upsampled by a factor of nfac, i.e. +% an order nfac*k quadrature rule is used on both panels + +if nargin < 2 + nfac = 1; +end + + +kover = ceil(k*nfac); +[xleg, wleg, ~, ~] = lege.exps(k); +if kover == k + xover = xleg; + wover = wleg; +else + [xover, wover, ~, ~] = lege.exps(kover); +end + +xover01 = (xover+1)/2; +wover01 = wover/2; + +xs0 = cell(k,1); +ws0 = cell(k,1); +for i=1:k + xn = xleg(i); + xls = xover01*(xn+1)-1; + wls = wover01*(xn+1); + xrs = xover01*(1-xn)+xn; + wrs = wover01*(1-xn); + xs0{i} = [xls; xrs]; + ws0{i} = [wls; wrs]; +end + +end + diff --git a/chunkie/+chnk/+quadggq/setup.m b/chunkie/+chnk/+quadggq/setup.m index a47bbdc7..66584ac0 100644 --- a/chunkie/+chnk/+quadggq/setup.m +++ b/chunkie/+chnk/+quadggq/setup.m @@ -11,6 +11,9 @@ % [[[---------------------------]-------------------]-------------------] % type = 'log' 'pv' 'hs' % +% f_{0}(x) + abs|x-x_{j}|*f_{1}(x) +% type = 'removable' +% % where the f_i are polynomials of order 2*k % % input: @@ -39,6 +42,8 @@ [xs0,wts0] = chnk.quadggq.gethqsuppquad(k,1); elseif strcmpi(type,'hs') [xs0,wts0] = chnk.quadggq.gethqsuppquad(k,2); + elseif strcmpi(type,'removable') + [xs0,wts0] = chnk.quadggq.getremovablequad(k,1); end ainterp1 = lege.matrin(k,xs1); @@ -58,4 +63,4 @@ auxquad.wts0 = wts0; auxquad.ainterp1 = ainterp1; auxquad.ainterps0 = ainterps0; - \ No newline at end of file + diff --git a/chunkie/@kernel/helm1d.m b/chunkie/@kernel/helm1d.m new file mode 100644 index 00000000..eb89bc05 --- /dev/null +++ b/chunkie/@kernel/helm1d.m @@ -0,0 +1,35 @@ +function obj = helm1d(type, zk, coefs) +%KERNEL.HELM2D Construct the Helmholtz kernel. +% KERNEL.HELM2D('s', ZK) or KERNEL.HELM2D('single', ZK) constructs the +% single-layer Helmholtz kernel with wavenumber ZK. +% +% +% See also CHNK.HELM2D.KERN. + +if ( nargin < 1 ) + error('Missing Helmholtz kernel type.'); +end + +if ( nargin < 2 ) + error('Missing Helmholtz wavenumber.'); +end + +obj = kernel(); +obj.name = 'helmholtz1d'; +obj.params.zk = zk; +obj.opdims = [1 1]; + +switch lower(type) + + case {'s', 'single'} + obj.type = 's'; + obj.eval = @(s,t) chnk.helm1d.kern(zk, s, t, 's'); + obj.fmm = []; + obj.sing = 'removable'; + + otherwise + error('Unknown Helmholtz kernel type ''%s''.', type); + +end + +end diff --git a/chunkie/@kernel/kernel.m b/chunkie/@kernel/kernel.m index f8248115..9c69973c 100644 --- a/chunkie/@kernel/kernel.m +++ b/chunkie/@kernel/kernel.m @@ -9,6 +9,7 @@ % ---- ---- % 'laplace' ('lap', 'l') 's', 'd', 'sp', 'c' % 'helmholtz' ('helm', 'h') 's', 'd', 'sp', 'dp', 'c' +% 'helmholtz1d' ('helm1d', 'h1d') 's' % 'helmholtz difference' ('helmdiff', 'hdiff') 's', 'd', 'sp', 'dp' % 'elasticity' ('elast', 'e') 's', 'strac', 'd', 'dalt' % 'stokes' ('stok', 's') 'svel', 'spres', 'strac', @@ -70,6 +71,8 @@ obj = kernel.lap2d(varargin{:}); case {'helmholtz', 'helm', 'h'} obj = kernel.helm2d(varargin{:}); + case {'helmholtz1d', 'helm1d', 'h1d'} + obj = kernel.helm1d(varargin{:}); case {'helmholtz difference', 'helmdiff', 'hdiff'} obj = kernel.helm2ddiff(varargin{:}); case {'stokes', 'stok', 's'} @@ -105,6 +108,7 @@ obj = lap2d(varargin); obj = helm2d(varargin); + obj = helm1d(varargin); obj = helm2ddiff(varargin); obj = stok2d(varargin); obj = elast2d(varargin); diff --git a/chunkie/chunkermat.m b/chunkie/chunkermat.m index 5e1897fc..ec6acd36 100644 --- a/chunkie/chunkermat.m +++ b/chunkie/chunkermat.m @@ -30,6 +30,7 @@ % is not defined by the kernel object. Supported % types are: % smooth => smooth kernels +% removable => piecewise smooth kernels % log => logarithmically singular kernels or % smooth times log + smooth % pv => principal value singular kernels + log @@ -134,6 +135,9 @@ nsub = 40; adaptive_correction = false; sing = 'log'; +if ~isempty(kern.sing) + sing = kern.sing; +end % get opts from struct if available @@ -346,6 +350,15 @@ auxquads = chnk.quadggq.setup(k,type); opts.auxquads.ggqlog = auxquads; end + elseif strcmpi(sing, 'removable') + type = 'removable'; + if (isfield(opts,'auxquads') && isfield(opts.auxquads,'ggqremovable')) + auxquads = opts.auxquads.ggqremovable; + else + k = chnkr.k; + auxquads = chnk.quadggq.setup(k, type); + opts.auxquads.ggqremovable = auxquads; + end elseif strcmpi(singi,'log') type = 'log'; if (isfield(opts,'auxquads') &&isfield(opts.auxquads,'ggqlog')) From c981e1a7baf0bd338a6db60125ca5469ed08aaa4 Mon Sep 17 00:00:00 2001 From: Manas Rachh Date: Mon, 10 Feb 2025 22:17:33 -0500 Subject: [PATCH 2/2] fixing small bug in chunkermat --- chunkie/chunkermat.m | 3 --- 1 file changed, 3 deletions(-) diff --git a/chunkie/chunkermat.m b/chunkie/chunkermat.m index 41369fe3..9e800990 100644 --- a/chunkie/chunkermat.m +++ b/chunkie/chunkermat.m @@ -138,9 +138,6 @@ nsub = 40; adaptive_correction = false; sing = 'log'; -if ~isempty(kern.sing) - sing = kern.sing; -end % get opts from struct if available