diff --git a/chunkie/+chnk/+quadggq/getremovablequad.m b/chunkie/+chnk/+quadggq/getremovablequad.m new file mode 100644 index 0000000..186da34 --- /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 a47bbdc..66584ac 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 0000000..eb89bc0 --- /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 a5c9e2b..39efb84 100644 --- a/chunkie/@kernel/kernel.m +++ b/chunkie/@kernel/kernel.m @@ -10,6 +10,7 @@ % 'laplace' ('lap', 'l') 's', 'd', 'sp', 'c' % 'helmholtz' ('helm', 'h') 's', 'd', 'sp', 'dp', 'c' % 'cp' +% 'helmholtz1d' ('helm1d', 'h1d') 's' % 'helmholtz difference' 's', 'd', 'sp', 'dp' % ('helmdiff', 'hdiff', 'helm_diff') % 'elasticity' ('elast', 'e') 's', 'strac', 'd', 'dalt' @@ -100,6 +101,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', 'helm_diff'} obj = kernel.helm2ddiff(varargin{:}); case {'stokes', 'stok', 's'} @@ -157,6 +160,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 1c48c6f..4b3b0b3 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 @@ -379,6 +380,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'))