Merge branch 'master' of https://github.com/fastalgorithms/chunkie in…

…to rcip-interp

.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

README.md

@@ -76,7 +76,7 @@ and optionally on the
 the library will not function without FLAM and its subdirectories included in the matlab path.
 
 
-## chunkie developers and antecedents
+## chunkie team
 
 chunkie has benefitted from the contributions of several developers: Travis Askham, 
 Manas Rachh, Michael O'Neil, Jeremy Hoskins, Dan Fortunato, Shidong Jiang, 

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/kernbyindex.m

@@ -1,5 +1,5 @@
 % TODO: rewrite this documentation with support for multiple chunkers
-function mat = kernbyindex(i,j,chnkrs,whts,kern,opdims_mat,spmat,l2scale)
+function mat = kernbyindex(i,j,chnkobj,kern,opdims_mat,spmat,l2scale)
 %% 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.
@@ -26,7 +26,6 @@
 % j - array of col indices to compute
 % chnkr - chunker object describing boundary, 
 %         or chunkgraph.echunks, which is an array of chunker objects.
-% whts - smooth weights on chunker object (already repeated according to opdims)
 % kern - kernel function of the form kern(s,t,stau,ttau) where s and t 
 %    are source and target points and stau and ttau are the local unit
 %    tangents
@@ -42,14 +41,26 @@
     kern = kern.eval;
 end
 
-if nargin < 8
+if nargin < 7
     l2scale = false;
 end
 
-% find unique underlying points
+if class(chnkobj) == "chunker"
+    chnkrs = chnkobj;
+elseif class(chnkobj) == "chunkgraph"
+    chnkrs = chnkobj.echnks;
+else
+    msg = "In CHNK.FLAM.KERNBYINDEX, invalid chunk object";
+    error(msg);
+end
 
 
 nchunkers = length(chnkrs);
+if length(size(opdims_mat)) == 2
+    opdims_mat = opdims_mat(:);
+    opdims_mat = repmat(opdims_mat,[nchunkers nchunkers]);
+    opdims_mat = reshape(opdims_mat,[2 nchunkers nchunkers]);
+end
 
 
 irowlocs = zeros(nchunkers+1,1);
@@ -79,10 +90,31 @@
     mat_row_end   = irowlocs(itrg+1)-1;
     flag_i = i>=mat_row_start & i<=mat_row_end;
     mat_row_ind = i(flag_i);
+
+
 
-    if length(mat_row_ind) == 0
+    if isempty(mat_row_ind)
         continue
     end
+    % Note using opdims from first column, opdims(1) has to be
+    % constant in this row.
+    opdims = opdims_mat(:,itrg,1);
+    ipts = idivide(int64(mat_row_ind(:)-mat_row_start),int64(opdims(1)))+1;
+    [iuni,~,iiuni] = unique(ipts);
+
+    iiuni2 = (iiuni-1)*opdims(1) + mod(mat_row_ind(:)-mat_row_start,opdims(1))+1;
+
+    ri = chnkr_trg.r(:,iuni);
+    di = chnkr_trg.d(:,iuni); 
+    ni = chnkr_trg.n(:,iuni);
+    d2i = chnkr_trg.d2(:,iuni);
+    targinfo = []; targinfo.r = ri; targinfo.d = di; targinfo.d2 = d2i;
+        targinfo.n = ni;
+
+    wtarg = chnkr_trg.wts(iuni); 
+    wtarg = repmat(wtarg(:).',opdims(1),1);
+    wtarg = wtarg(:);
+
 
     for isrc=1:nchunkers
         chnkr_src = chnkrs(isrc);
@@ -91,52 +123,50 @@
         flag_j = j>=mat_col_start & j<=mat_col_end;
         mat_col_ind   = j(flag_j);
 
-        if length(mat_col_ind) == 0
+        if isempty(mat_col_ind)
             continue
         end
 
         opdims = opdims_mat(:,itrg,isrc);
-        ipts = idivide(int64(mat_row_ind(:)-mat_row_start),int64(opdims(1)))+1;
+
         jpts = idivide(int64(mat_col_ind(:)-mat_col_start),int64(opdims(2)))+1;
-
-        [iuni,~,iiuni] = unique(ipts);
         [juni,~,ijuni] = unique(jpts);
 
-        ri = chnkr_trg.r(:,iuni); rj = chnkr_src.r(:,juni);
-        di = chnkr_trg.d(:,iuni); dj = chnkr_src.d(:,juni);
-        nj = chnkr_src.n(:,juni); ni = chnkr_trg.n(:,iuni);
-        d2i = chnkr_trg.d2(:,iuni); d2j = chnkr_src.d2(:,juni);
+        rj = chnkr_src.r(:,juni);
+        dj = chnkr_src.d(:,juni);
+        nj = chnkr_src.n(:,juni); 
+        d2j = chnkr_src.d2(:,juni);
         srcinfo = []; srcinfo.r = rj; srcinfo.d = dj; srcinfo.d2 = d2j;
         srcinfo.n = nj;
-        targinfo = []; targinfo.r = ri; targinfo.d = di; targinfo.d2 = d2i;
-        targinfo.n = ni;
+
+        wsrc = chnkr_src.wts(juni); 
+        wsrc = repmat( (wsrc(:)).', opdims(2), 1);
+        wsrc = ( wsrc(:) ).';
+
 
-        if length(kern) == 1
+        if size(kern) == 1
             matuni = kern(srcinfo,targinfo);
         else
             matuni = kern{itrg,isrc}(srcinfo,targinfo);
         end
-        iiuni2 = (iiuni-1)*opdims(1) + mod(mat_row_ind(:)-mat_row_start,opdims(1))+1;
+
+        % scale matrix by weights
+        if(l2scale)
+            matuni = sqrt(wtarg) .* matuni .* sqrt(wsrc);
+        else
+            matuni = matuni .* wsrc;
+        end
+
+
         ijuni2 = (ijuni-1)*opdims(2) + mod(mat_col_ind(:)-mat_col_start,opdims(2))+1;
         mat(flag_i,flag_j) = matuni(iiuni2,ijuni2);
     end
 end
 
-% scale columns by weights
-
-whts = whts(:);
-
-if l2scale    
-
-    mat = sqrt(whts(i)) .* mat .* sqrt(whts(j).');
-
-else
-    mat = mat .* whts(j).' ;
-end
 
 % overwrite any entries given as nonzeros in sparse matrix
 
-if nargin > 6
+if nargin > 5
     [isp,jsp,vsp] = find(spmat(i,j));
     linsp = isp + (jsp-1)*length(i(:));
     mat(linsp) = vsp;
@@ -145,51 +175,3 @@
 return;
 
 
-% the original single chunker code
-
-% ipts = idivide(int64(i(:)-1),int64(opdims(1)))+1;
-% jpts = idivide(int64(j(:)-1),int64(opdims(2)))+1;
-
-% [iuni,~,iiuni] = unique(ipts);
-% [juni,~,ijuni] = unique(jpts);
-
-% % matrix-valued entries of kernel for unique points
-
-% ri = chnkr.r(:,iuni); rj = chnkr.r(:,juni);
-% di = chnkr.d(:,iuni); dj = chnkr.d(:,juni);
-% nj = chnkr.n(:,juni); ni = chnkr.n(:,iuni);
-% d2i = chnkr.d2(:,iuni); d2j = chnkr.d2(:,juni);
-% srcinfo = []; srcinfo.r = rj; srcinfo.d = dj; srcinfo.d2 = d2j;
-% srcinfo.n = nj;
-% targinfo = []; targinfo.r = ri; targinfo.d = di; targinfo.d2 = d2i;
-% targinfo.n = ni;
-% %di = bsxfun(@rdivide,di,sqrt(sum(di.^2,1)));
-% %dj = bsxfun(@rdivide,dj,sqrt(sum(dj.^2,1)));
-
-% matuni = kern(srcinfo,targinfo);
-
-% % relevant rows and columns in matuni
-
-% iiuni2 = (iiuni-1)*opdims(1) + mod(i(:)-1,opdims(1))+1;
-% ijuni2 = (ijuni-1)*opdims(2) + mod(j(:)-1,opdims(2))+1;
-
-% mat = matuni(iiuni2,ijuni2);
-
-% whts = whts(:);
-
-% % scale columns by weights
-
-% if l2scale
-%     mat = sqrt(whts(ipts(:))) .* mat .* sqrt(whts(jpts(:)).');
-% else
-%     mat = mat .* whts(jpts(:)).';
-% end
-% % overwrite any entries given as nonzeros in sparse matrix
-
-% if nargin > 6
-%     [isp,jsp,vsp] = find(spmat(i,j));
-%     linsp = isp + (jsp-1)*length(i(:));
-%     mat(linsp) = vsp;
-% end
-
-% end

chunkie/+chnk/+flam/kernbyindexr.m

@@ -1,4 +1,4 @@
-function mat = kernbyindexr(i,j,targs,chnkr,whts,kern,opdims,spmat)
+function mat = kernbyindexr(i,j,targs,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
@@ -63,12 +63,13 @@
 
 % scale columns by weights
 
+whts = chnkr.wts;
 wj = whts(jpts(:));
 mat = bsxfun(@times,mat,wj.');
 
 % overwrite any entries given as nonzeros in sparse matrix
 
-if nargin > 7
+if nargin > 6
     [isp,jsp,vsp] = find(spmat(i,j));
     linsp = isp + (jsp-1)*length(i(:));
     mat(linsp) = vsp;