merging axissym kernels and master, green_laptest fails

.readthedocs.yaml

@@ -0,0 +1,32 @@
+# .readthedocs.yaml
+# Read the Docs configuration file
+# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
+
+# Required
+version: 2
+
+# Set the OS, Python version and other tools you might need
+build:
+  os: ubuntu-22.04
+  tools:
+    python: "3.12"
+    # You can also specify other tool versions:
+    # nodejs: "19"
+    # rust: "1.64"
+    # golang: "1.19"
+
+# Build documentation in the "docs/" directory with Sphinx
+sphinx:
+   configuration: docs/conf.py
+
+# Optionally build your docs in additional formats such as PDF and ePub
+# formats:
+#    - pdf
+#    - epub
+
+# Optional but recommended, declare the Python requirements required
+# to build your documentation
+# See https://docs.readthedocs.io/en/stable/guides/reproducible-builds.html
+python:
+  install:
+    - requirements: docs/requirements.txt

chunkie/+chnk/+elast2d/kern.m

@@ -33,12 +33,17 @@
 % - type: string 
 %   * type == 'S', single layer kernel
 %   * type == 'Strac', traction of single layer kernel
+%   * type == 'Sgrad', gradient of single layer kernel
+%                kernel is 4 x 2, where entries are organized as
+%                          Sgrad = [grad S(1,1) grad S(1,2) 
+%                                   grad S(2,1) grad S(2,2)] 
 %   * type == 'D', double layer kernel
 %   * type == 'Dalt', alternative (smooth) double layer kernel  
 %   * type == 'Daltgrad', gradient of alternative (smooth) double layer 
 %                kernel is 4 x 2, where entries are organized as
 %                          Daltgrad = [grad Dalt(1,1) grad Dalt(1,2) 
 %                                      grad Dalt(2,1) grad Dalt(2,2)] 
+%   * type == 'Dalttrac', traction of alternative (smooth) double layer 
 %               
 % outpts 
 % - mat: (2 nt) x (2 ns) matrix 
@@ -53,8 +58,8 @@
 eta = mu/(2*pi*(lam+2*mu));
 zeta = (lam+mu)/(pi*(lam+2*mu));
 
-n = size(s.r,2);
-m = size(t.r,2);
+n = size(s.r(:,:),2);
+m = size(t.r(:,:),2);
 
 x = (t.r(1,:)).' - s.r(1,:);
 y = (t.r(2,:)).' - s.r(2,:);
@@ -81,6 +86,37 @@
         kron(rdotv./r2,[1 0; 0 1]));
     mat = term1+term2;
 end
+if (strcmpi(type,'sgrad'))
+    mat = zeros(4*m,2*n);
+    tmp = beta*x./r2;
+    mat(1:4:end,1:2:end) = tmp;
+    mat(3:4:end,2:2:end) = tmp;
+    tmp = beta*y./r2;
+    mat(2:4:end,1:2:end) = tmp;
+    mat(4:4:end,2:2:end) = tmp;
+
+    % x der of x^2/ r^2
+    tmp = gamma*(2*r2.*x - x.^2.*(2*x))./r4;
+    mat(1:4:end,1:2:end) = mat(1:4:end,1:2:end) + tmp;
+    % y der of x^2/ r^2
+    tmp = gamma*(-x.^2.*(2*y))./r4;
+    mat(2:4:end,1:2:end) = mat(2:4:end,1:2:end) + tmp;
+    % x der of xy/ r^2
+    tmp = gamma*(r2.*y-x.*y.*(2*x))./r4;
+    mat(3:4:end,1:2:end) = mat(3:4:end,1:2:end) + tmp;
+    mat(1:4:end,2:2:end) = mat(1:4:end,2:2:end) + tmp;
+    % y der of xy/ r^2
+    tmp = gamma*(r2.*x-x.*y.*(2*y))./r4;
+    mat(4:4:end,1:2:end) = mat(4:4:end,1:2:end) + tmp;
+    mat(2:4:end,2:2:end) = mat(2:4:end,2:2:end) + tmp;
+    % x der of y^2/ r^2
+    tmp = gamma*(- y.^2.*(2*x))./r4;
+    mat(3:4:end,2:2:end) = mat(3:4:end,2:2:end) + tmp;
+    % y der of y^2/ r^2
+    tmp = gamma*(2*r2.*y-y.^2.*(2*y))./r4;
+    mat(4:4:end,2:2:end) = mat(4:4:end,2:2:end) + tmp;
+
+end
 if (strcmpi(type,'d'))
     dirx = s.n(1,:); dirx = dirx(:).';
     diry = s.n(2,:); diry = diry(:).';
@@ -151,4 +187,59 @@
 
 end
 
+if (strcmpi(type,'dalttrac'))
+    dirx = s.n(1,:); dirx = dirx(:).';
+    diry = s.n(2,:); diry = diry(:).';
+    rdotv = x.*dirx + y.*diry;
+    r6 = r4.*r2;
+
+    matg = zeros(4*m,2*n);
+
+    % i = 1, j = 1, l = 1
+    aij_xl = -zeta*(-4*x.^3.*rdotv./r6+(2*x.*rdotv+x.^2.*dirx)./r4) ...
+        -2*eta*(dirx./r2 - 2*rdotv.*x./r4);
+    matg(1:4:end,1:2:end) = aij_xl;
+
+    % i = 1, j = 1, l = 2
+    aij_xl = -zeta*(-4*x.^2.*y.*rdotv./r6+(x.^2.*diry)./r4) ...
+        -2*eta*(diry./r2 - 2*rdotv.*y./r4);
+    matg(2:4:end,1:2:end) = aij_xl;
+
+    % i = 2, j = 1, l = 1
+    aij_xl = -zeta*(-4*x.^2.*y.*rdotv./r6+(y.*rdotv+x.*y.*dirx)./r4);
+    matg(3:4:end,1:2:end) = aij_xl;
+
+    % i = 2, j = 1, l = 2
+    aij_xl = -zeta*(-4*x.*y.^2.*rdotv./r6+(x.*rdotv+x.*y.*diry)./r4);
+    matg(4:4:end,1:2:end) = aij_xl;
+
+    % i = 1, j = 2, l = 1
+    aij_xl = -zeta*(-4*x.^2.*y.*rdotv./r6+(y.*rdotv+x.*y.*dirx)./r4);
+    matg(1:4:end,2:2:end) = aij_xl;
+
+    % i = 1, j = 2, l = 2
+    aij_xl = -zeta*(-4*x.*y.^2.*rdotv./r6+(x.*rdotv + x.*y.*diry)./r4);
+    matg(2:4:end,2:2:end) = aij_xl;
+
+    % i = 2, j = 2, l = 1
+    aij_xl = -zeta*(-4*y.^2.*x.*rdotv./r6+(y.^2.*dirx)./r4) ...
+        -2*eta*(dirx./r2 - 2*rdotv.*x./r4);
+    matg(3:4:end,2:2:end) = aij_xl;
+
+    % i = 2, j = 2, l = 2
+    aij_xl = -zeta*(-4*y.^3.*rdotv./r6+(2*y.*rdotv+y.^2.*diry)./r4) ...
+        -2*eta*(diry./r2 - 2*rdotv.*y./r4);
+    matg(4:4:end,2:2:end) = aij_xl;
+
+    mat = zeros(2*m,2*n);
+    n1 = t.n(1,:); n1 = n1(:); n2 = t.n(2,:); n2 = n2(:);
+    tmp = matg(1:4:end,:)+matg(4:4:end,:);
+    mat(1:2:end,:) = lam*n1.*tmp;
+    mat(2:2:end,:) = lam*n2.*tmp;
+    tmp = mu*(matg(2:4:end,:) + matg(3:4:end,:));
+    mat(1:2:end,:) = mat(1:2:end,:) + tmp.*n2 + 2*mu*matg(1:4:end,:).*n1;
+    mat(2:2:end,:) = mat(2:2:end,:) + tmp.*n1 + 2*mu*matg(4:4:end,:).*n2;
+
+end
+
 end

chunkie/+chnk/+flam/kernbyindexr.m

@@ -1,4 +1,4 @@
-function mat = kernbyindexr(i,j,targs,chnkr,kern,opdims,spmat)
+function mat = kernbyindexr(i,j,targobj,chnkr,kern,opdims,spmat)
 %% evaluate system matrix by entry index utility function for 
 % general kernels, with replacement for specific entries and the
 % ability to add a low-rank modification, rectangular version
@@ -17,6 +17,10 @@
 %
 % i - array of row indices to compute
 % j - array of col indices to compute
+% targobj - object describing the target points, can be specified as
+%       * array of points
+%       * chunker object
+%       * chunkgraph object
 % chnkr - chunker object describing boundary
 % whts - smooth integration weights on chnkr
 % kern - kernel function of the form kern(s,t,stau,ttau) where s and t 
@@ -42,15 +46,32 @@
 [juni,~,ijuni] = unique(jpts);
 
 % matrix-valued entries of kernel for unique points
-
-ri = targs(:,iuni); rj = chnkr.r(:,juni);
+rj = chnkr.r(:,juni);
 dj = chnkr.d(:,juni); d2j = chnkr.d2(:,juni);
 nj = chnkr.n(:,juni);
-%di = bsxfun(@rdivide,di,sqrt(sum(di.^2,1)));
-%dj = bsxfun(@rdivide,dj,sqrt(sum(dj.^2,1)));
 srcinfo = []; srcinfo.r = rj; srcinfo.d = dj; srcinfo.d2 = d2j;
 srcinfo.n = nj;
-targinfo = []; targinfo.r = ri;
+% Assign appropriate object to targinfo
+targinfo = [];
+if isa(targobj, "chunker") || isa(targobj, "chunkgraph")
+    targinfo.r = targobj.r(:,iuni);
+    targinfo.d = targobj.d(:,iuni);
+    targinfo.d2 = targobj.d2(:,iuni);
+    targinfo.n = targobj.n(:,iuni);
+elseif isstruct(targobj)
+    if isfield(targobj,'d')
+        targinfo.d = targobj.d(:,iuni);
+    end
+    if isfield(targobj,'d2')
+        targinfo.d = targobj.d(:,iuni);
+    end
+    if isfield(targobj,'n')
+        targinfo.n = targobj.n(:,iuni);
+    end
+    targinfo.r = targobj.r(:,iuni);
+else
+    targinfo.r = targobj(:,iuni);
+end
 
 matuni = kern(srcinfo,targinfo);
 

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/+helm2d/transmission_helper.m

@@ -151,7 +151,7 @@
                 c1 = coefs(d1);
                 c2 = coefs(d2);
                 ind1 = sum(nchs(1:i-1))*ngl*2+(1:2:2*nchs(i)*ngl);
-                ind2 = sum(nchs(1:i-1))*ngl*2+(1:2:2*nchs(i)*ngl);
+                ind2 = sum(nchs(1:i-1))*ngl*2+(2:2:2*nchs(i)*ngl);
                 targnorm = chnkrs(i).n;
                 nx = targnorm(1,:); nx = nx(:);
                 ny = targnorm(2,:); ny = ny(:);

chunkie/+chnk/+quadadap/buildmat.m

@@ -22,7 +22,6 @@
 adj = chnkr.adj;
 d = chnkr.d;
 d2 = chnkr.d2;
-h = chnkr.h;
 n = chnkr.n;
 
 [t,wts,u] = lege.exps(k);
@@ -77,7 +76,7 @@
         dt = d(:,:,ibefore);
         d2t = d2(:,:,ibefore);
         nt = n(:,:,ibefore);
-        submat = chnk.adapgausswts(r,d,n,d2,h,t,bw,i,rt,dt,nt,d2t, ...
+        submat = chnk.adapgausswts(r,d,n,d2,t,bw,i,rt,dt,nt,d2t, ...
             kern,opdims,t2,w2);
 
         imat = 1 + (ibefore-1)*k*opdims(1);
@@ -90,7 +89,7 @@
         dt = d(:,:,iafter);
         nt = n(:,:,iafter);
         d2t = d2(:,:,iafter);
-        submat = chnk.adapgausswts(r,d,n,d2,h,t,bw,i,rt,dt,nt,d2t, ...
+        submat = chnk.adapgausswts(r,d,n,d2,t,bw,i,rt,dt,nt,d2t, ...
             kern,opdims,t2,w2);
 
         imat = 1 + (iafter-1)*k*opdims(1);
@@ -101,7 +100,7 @@
 
     % self
 
-    submat = chnk.quadggq.diagbuildmat(r,d,n,d2,h,[],i,kern,opdims,...
+    submat = chnk.quadggq.diagbuildmat(r,d,n,d2,[],i,kern,opdims,...
         xs0,wts0,ainterps0kron,ainterps0);
 
     imat = 1 + (i-1)*k*opdims(1);
@@ -163,7 +162,7 @@
         d2t = d2(:,targfix);
         nt = n(:,targfix);
 
-        submat = chnk.adapgausswts(r,d,n,d2,h,t,bw,i,rt,dt,nt,d2t, ...
+        submat = chnk.adapgausswts(r,d,n,d2,t,bw,i,rt,dt,nt,d2t, ...
                 kern,opdims,t2,w2);
 
         imats = bsxfun(@plus,(1:opdims(1)).',opdims(1)*(targfix(:)-1).');