Merge pull request #152 from fastalgorithms/removable-quad

Removable quad

Author: jghoskins

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
+    

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;
-  
+  

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

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);

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'))