more stuff in fixing chunker mat, adding non-legendre points option t…

…o proxyfun, and derivative of c in kernel in helmholtz, fixing small typo of abs vs sqrt(z2)

chunkie/+chnk/+flam/proxy_square_pts.m

@@ -1,4 +1,4 @@
-function [pr,ptau,pw,pin] = proxy_square_pts(porder)
+function [pr,ptau,pw,pin] = proxy_square_pts(porder, opts)
 %PROXY_SQUARE_PTS return the proxy points, unit tangents, weights, and 
 % function handle for determining if points are within proxy surface.
 % This function is for the proxy surface around a unit box centered at 
@@ -16,12 +16,31 @@
     porder = 64;
 end
 
+if nargin < 2
+    opts = [];
+end
+
+iflege = 0;
+if isfield(opts, 'iflege')
+    iflege = opts.iflege;
+end
+
+
 assert(mod(porder,4) == 0, ...
     'number of proxy points on square should be multiple of 4')
 
 po4 = porder/4;
 
-pts = -1.5+3*(0:(po4-1))*1.0/po4;
+if ~iflege
+    pts = -1.5+3*(0:(po4-1))*1.0/po4;
+    pw = ones(porder,1)*(4.0*3.0/porder);
+else
+    [xleg, wleg] = lege.exps(po4);
+    pts = xleg.'*1.5;
+    wleg = wleg*1.5;   
+    pw = [wleg; wleg; wleg; wleg];
+end
+
 one = ones(size(pts));
 
 pr = [pts, one*1.5, -pts, -1.5*one;
@@ -30,8 +49,6 @@
 ptau = [one, one*0, -one, one*0;
     one*0, one, one*0, -one];
 
-pw = ones(porder,1)*(4.0*3.0/porder);
-
 pin = @(x) max(abs(x),[],1) < 1.5;
 
 end

chunkie/+chnk/+helm2d/fmm.m

@@ -38,6 +38,8 @@
 %                   note that targinfo must contain targ.n
 %                type = 'dprime' normal derivative of double layer,
 %                   note that targinfo must contain targ.n
+%                type = 'cprime' normal derivative of combined layer,
+%                   note that targinfo must contain targ.n
 %   sigma - density
 %   varargin{1} - coef in the combined layer formula, otherwise
 %                does nothing
@@ -61,7 +63,7 @@
     case {'d', 'dprime'}
         srcuse.dipstr = sigma(:).';
         srcuse.dipvec = srcinfo.n(1:2,:);
-    case 'c'
+    case {'c', 'cprime'}
         coefs = varargin{1};
         srcuse.charges = coefs(2)*sigma(:).';
         srcuse.dipstr  = coefs(1)*sigma(:).';
@@ -83,7 +85,7 @@
     switch lower(type)
         case {'s', 'd', 'c'}
             varargout{1} = U.pottarg.';
-        case {'sprime', 'dprime'}
+        case {'sprime', 'dprime', 'cprime'}
             if ( ~isfield(targinfo, 'n') )
                 error('CHUNKIE:helm2d:fmm:normals', ...
                     'Targets require normal info when evaluating Helmholtz kernel ''%s''.', type);

chunkie/+chnk/+helm2d/kern.m

@@ -133,14 +133,14 @@
 % normal derivative of combined field
 if strcmpi(type,'cprime')
   coef = ones(2,1);
-  if(nargin == 5); coef = varargin{1}; end;
+  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);
+  [~,grad,hess] = chnk.helm2d.green(zk,src,targ);
 
   nxsrc = repmat(srcnorm(1,:),nt,1);
   nysrc = repmat(srcnorm(2,:),nt,1);
@@ -161,7 +161,7 @@
 % Dirichlet and neumann data corresponding to combined field
 if strcmpi(type,'c2trans') 
   coef = ones(2,1);
-  if(nargin == 5); coef = varargin{1}; end;
+  if(nargin == 5); coef = varargin{1}; end
   targnorm = targinfo.n;
   srcnorm = srcinfo.n;
 

chunkie/+chnk/+quadnative/buildmat.m

@@ -36,12 +36,7 @@
 targinfo = []; targinfo.r = rt; targinfo.d = dt; targinfo.d2 = d2t; targinfo.n = nt;
 hs = h(j); ht = h(i);
 
-dsnrms = sqrt(sum(abs(ds).^2,1));
-%taus = bsxfun(@rdivide,ds,dsnrms);
-
-%dtnrms = sqrt(sum(abs(dt).^2,1));
-%taut = bsxfun(@rdivide,dt,dtnrms);
-
+dsnrms = sqrt(sum(ds.^2,1));
 ws = kron(hs(:),wts(:));
 
 dsdt = dsnrms(:).*ws;

chunkie/@kernel/helm2d.m

@@ -14,6 +14,10 @@
 %   parameter ETA, i.e., COEFS(1)*KERNEL.HELM2D('d', ZK) + 
 %   COEFS(2)*KERNEL.HELM2D('s', ZK).
 %
+%   KERNEL.HELM2D('cp', ZK, COEFS) or KERNEL.HELM2D('combined', ZK, COEFS)
+%   constructs the derivative of the combined-layer Helmholtz kernel with 
+%   wavenumber ZK and  parameter ETA, i.e., COEFS(1)*KERNEL.HELM2D('dp', ZK) + 
+%   COEFS(2)*KERNEL.HELM2D('sp', ZK).
 % See also CHNK.HELM2D.KERN.
 
 if ( nargin < 1 )
@@ -66,6 +70,17 @@
         obj.fmm  = @(eps,s,t,sigma) chnk.helm2d.fmm(eps, zk, s, t, 'c', sigma, coefs);
         obj.sing = 'log';
 
+    case {'cp', 'cprime'}
+        if ( nargin < 3 )
+            warning('Missing combined layer coefficients. Defaulting to [1,1i].');
+            coefs = [1,1i];
+        end
+        obj.type = 'cprime';
+        obj.params.coefs = coefs;
+        obj.eval = @(s,t) chnk.helm2d.kern(zk, s, t, 'cprime', coefs);
+        obj.fmm  = @(eps,s,t,sigma) chnk.helm2d.fmm(eps, zk, s, t, 'cprime', sigma, coefs);
+        obj.sing = 'hs';
+
     otherwise
         error('Unknown Helmholtz kernel type ''%s''.', type);
 

chunkie/@kernel/kernel.m

@@ -9,6 +9,7 @@
 %      ----                                         ----
 %      'laplace'    ('lap', 'l')                    's', 'd', 'sp', 'c'
 %      'helmholtz'  ('helm', 'h')                   's', 'd', 'sp', 'dp', 'c'
+%                                                   'cp'
 %      'helmholtz difference' ('helmdiff', 'hdiff') 's', 'd', 'sp', 'dp'
 %      'elasticity' ('elast', 'e')                  's', 'strac', 'd', 'dalt'
 %      'stokes'     ('stok', 's')                   'svel', 'spres', 'strac',

chunkie/chunkerkernevalmat.m

@@ -6,8 +6,17 @@
 % Syntax: mat = chunkerkernevalmat(chnkr,kern,targs,opts)
 %
 % Input:
-%   chnkr - chunker object description of curve
-%   kern - integral kernel taking inputs kern(srcinfo,targinfo) 
+%   chnkr - chunker object describing boundary, currently
+%              only supports chunkers, and not chunkgraphs
+%   kern  - kernel function. By default, this should be a function handle
+%           accepting input of the form kern(srcinfo,targinfo), where srcinfo
+%           and targinfo are in the ptinfo struct format, i.e.
+%                ptinfo.r - positions (2,:) array
+%                ptinfo.d - first derivative in underlying
+%                     parameterization (2,:)
+%                ptinfo.n - unit normals (2,:)
+%                ptinfo.d2 - second derivative in underlying
+%                     parameterization (2,:)
 %   targobj - object describing the target points, can be specified as
 %       * array of points
 %       * chunker object
@@ -41,6 +50,37 @@
 
 % author: Travis Askham ([email protected])
 
+% convert kernel to kernel object, put in singularity info 
+% opts.sing provides a default value for singularities if not 
+% defined for kernels
+
+if isa(kern,'function_handle')
+    kern2 = kernel(kern);
+    kern = kern2;
+elseif isa(kern,'cell')
+    sz = size(kern);
+    kern2(sz(1),sz(2)) = kernel();
+    for j = 1:sz(2)
+        for i = 1:sz(1)
+            if isa(kern{i,j},'function_handle')
+                kern2(i,j) = kernel(kern{i,j});
+            elseif isa(kern{i,j},'kernel')
+                kern2(i,j) = kern{i,j};
+            else
+                msg = "Second input is not a kernel object, function handle, " ...
+                    + "or cell array";
+                error(msg);
+            end
+        end
+    end
+    kern = kern2;
+
+elseif ~isa(kern,'kernel')
+    msg = "Second input is not a kernel object, function handle, " ...
+                + "or cell array";
+    error(msg);
+end
+
 % determine operator dimensions using first two points
 
 
@@ -51,7 +91,9 @@
 targinfo.r = chnkr.r(:,2); targinfo.d = chnkr.d(:,2); 
 targinfo.d2 = chnkr.d2(:,2); targinfo.n = chnkr.n(:,2);
 
-ftemp = kern(srcinfo,targinfo);
+ftmp = kern.eval;
+
+ftemp = ftmp(srcinfo,targinfo);
 opdims = size(ftemp);
 
 if nargin < 4
@@ -96,24 +138,26 @@
 optsadap = []; 
 optsadap.eps = eps;
 
+
+
 if forcesmooth
-    mat = chunkerkernevalmat_smooth(chnkr,kern,opdims,targinfo, ...
+    mat = chunkerkernevalmat_smooth(chnkr,ftmp,opdims,targinfo, ...
         [],optssmooth);
     return
 end
 
 if forceadap
-    mat = chunkerkernevalmat_adap(chnkr,kern,opdims, ...
+    mat = chunkerkernevalmat_adap(chnkr,ftmp,opdims, ...
         targinfo,[],optsadap);
     return
 end
 
 if forcepqud
     optsflag = []; optsflag.fac = fac;
     flag = flagnear(chnkr,targinfo.r,optsflag);
-    spmat = chunkerkernevalmat_ho(chnkr,kern,opdims, ...
+    spmat = chunkerkernevalmat_ho(chnkr,ftmp,opdims, ...
         targinfo,flag,optsadap);
-    mat = chunkerkernevalmat_smooth(chnkr,kern,opdims,targinfo, ...
+    mat = chunkerkernevalmat_smooth(chnkr,ftmp,opdims,targinfo, ...
         flag,opts);
     mat = mat + spmat;
     return
@@ -123,15 +167,15 @@
 
 optsflag = []; optsflag.fac = fac;
 flag = flagnear(chnkr,targinfo.r,optsflag);
-spmat = chunkerkernevalmat_adap(chnkr,kern,opdims, ...
+spmat = chunkerkernevalmat_adap(chnkr,ftmp,opdims, ...
         targinfo,flag,optsadap);
 
 if nonsmoothonly
     mat = spmat;
     return;
 end
 
-mat = chunkerkernevalmat_smooth(chnkr,kern,opdims,targinfo, ...
+mat = chunkerkernevalmat_smooth(chnkr,ftmp,opdims,targinfo, ...
     flag,opts);
 
 mat = mat + spmat;

chunkie/chunkermat.m

@@ -7,7 +7,7 @@
 % Syntax: sysmat = chunkermat(chnkr,kern,opts)
 %
 % Input:
-%   chnkr - chunker object describing boundary
+%   chnkobj - chunker object describing boundary
 %   kern  - kernel function. By default, this should be a function handle
 %           accepting input of the form kern(srcinfo,targinfo), where srcinfo
 %           and targinfo are in the ptinfo struct format, i.e.