From d0d16d29cbe913cfaafc0909c8deac69464756b5 Mon Sep 17 00:00:00 2001
From: Manas Rachh <mrachh@flatironinstitute.org>
Date: Fri, 18 Oct 2024 15:19:30 -0700
Subject: [PATCH 01/15] smoother

---
 chunkie/+chnk/+smoother/get_pseudo_normals.m |  6 ++++++
 chunkie/+chnk/+smoother/smooth_curve.m       | 12 ++++++++++++
 2 files changed, 18 insertions(+)
 create mode 100644 chunkie/+chnk/+smoother/get_pseudo_normals.m
 create mode 100644 chunkie/+chnk/+smoother/smooth_curve.m

diff --git a/chunkie/+chnk/+smoother/get_pseudo_normals.m b/chunkie/+chnk/+smoother/get_pseudo_normals.m
new file mode 100644
index 00000000..e1110d61
--- /dev/null
+++ b/chunkie/+chnk/+smoother/get_pseudo_normals.m
@@ -0,0 +1,6 @@
+function pseudo_normals = get_pseudo_normals(verts)
+    [~, nv] = size(verts);
+    verts_ext = [verts, verts(:,1)];
+    d = verts_ext(:,2:end) - verts_ext(:,1:nv);
+    
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+smoother/smooth_curve.m b/chunkie/+chnk/+smoother/smooth_curve.m
new file mode 100644
index 00000000..05827387
--- /dev/null
+++ b/chunkie/+chnk/+smoother/smooth_curve.m
@@ -0,0 +1,12 @@
+nv = 10;
+z = exp(1j*2*pi*(1:nv-1)/nv);
+verts = [real(z); imag(z)];
+nchs = randi([1,10], 1, nv);
+
+umesh = [];
+umesh.verts = verts;
+umesh.pseudo_normals = get_pseudo_normal(verts);
+
+qmesh = [];
+qmesh
+dmesh = [];

From 2e6cd6d47bdc189389abde7449ccbd8f65b685b4 Mon Sep 17 00:00:00 2001
From: jghoskins <32715900+jghoskins@users.noreply.github.com>
Date: Fri, 18 Oct 2024 15:35:02 -0700
Subject: [PATCH 02/15] kernels!

---
 chunkie/+chnk/+smoother/expeval.m  | 76 ++++++++++++++++++++++++++++++
 chunkie/+chnk/+smoother/gpsi_all.m | 65 +++++++++++++++++++++++++
 chunkie/+chnk/+smoother/gpsi_loc.m | 48 +++++++++++++++++++
 chunkie/+chnk/+smoother/green.m    | 61 ++++++++++++++++++++++++
 chunkie/+chnk/+smoother/psi_eval.m |  6 +++
 5 files changed, 256 insertions(+)
 create mode 100644 chunkie/+chnk/+smoother/expeval.m
 create mode 100644 chunkie/+chnk/+smoother/gpsi_all.m
 create mode 100644 chunkie/+chnk/+smoother/gpsi_loc.m
 create mode 100644 chunkie/+chnk/+smoother/green.m
 create mode 100644 chunkie/+chnk/+smoother/psi_eval.m

diff --git a/chunkie/+chnk/+smoother/expeval.m b/chunkie/+chnk/+smoother/expeval.m
new file mode 100644
index 00000000..13c19f69
--- /dev/null
+++ b/chunkie/+chnk/+smoother/expeval.m
@@ -0,0 +1,76 @@
+function [cvals] = expeval(czs,ifexps)
+
+    sz = size(czs);
+    czs = czs(:);
+    cvals= zeros(size(czs));
+
+    ifloc = find(abs(czs) <= 2);
+    iffar = find(abs(czs) >2);
+    
+    if (numel(ifloc))
+    cloc = expeval_tayl(czs(ifloc));
+    cvals(ifloc) = cloc;
+    end
+    if (numel(iffar)>0)
+    cfar = expeval_lentz(czs(iffar));
+    cvals(iffar) = cfar;
+    end
+
+    if ((nargin>1) && ifexps)
+        cvals = cvals.*exp(-czs);
+    end
+
+    cvals = reshape(cvals,sz);
+
+end
+
+
+function [cvals] = expeval_tayl(czs)
+
+        gamma = 0.577215664901532860606512090082;
+
+        cvals = -gamma-log(czs);
+        cpow = czs;
+        rfac = 1;
+
+        eps = 1E-24;
+
+        for ii=1:31
+            cvals = cvals+cpow.*rfac/ii;
+            cpow = -czs.*cpow;
+            rfac = rfac/(ii+1);
+            % if (max(abs(cpow.*rfac)) < eps)
+            %     break;
+            % end
+        end
+
+         cvals = cvals.*exp(czs);
+
+end
+
+
+function [cvals] = expeval_lentz(czs)
+
+        nterms = 40;
+        cf = 1./czs;
+        cd = cf;
+        cv = 1/1E-30;
+
+        for ii=1:80
+        can = ii;
+        cd  = 1./(1+ii*cd);
+        cv  = 1 + ii./cv;
+        del = cd.*cv;
+        cf  = cf.*del;
+        cd  = 1./(czs+ii*cd);
+        cv  = czs + ii./cv;
+        del = cd.*cv;
+        cf  = cf.*del;
+        % if (max(abs(del-1)) < 1E-14)
+        %     break;
+        % end
+        end
+
+        cvals = cf;
+
+end
diff --git a/chunkie/+chnk/+smoother/gpsi_all.m b/chunkie/+chnk/+smoother/gpsi_all.m
new file mode 100644
index 00000000..37e69076
--- /dev/null
+++ b/chunkie/+chnk/+smoother/gpsi_all.m
@@ -0,0 +1,65 @@
+function [val,grad,hess,hess_sig] = gpsi_all(dx,dy,dr,dsigt)
+
+    [val,vald] = chnk.smoother.psi_eval(dr,dsigt);
+
+    da = 1/(4*pi);
+    z  = dr.^2./(2*dsigt.*dsigt);
+    val = val + log(z)*da;
+
+    if (nargout==1) 
+        return
+    end
+
+    grad = zeros(numel(dx),2);
+    grad(:,1) = vald.*dx./dr;
+    grad(:,2) = vald.*dy./dr;
+
+    grad(:,1) = grad(:,1)+2*dx./(dr.*dr)*da;
+    grad(:,2) = grad(:,2)+2*dy./(dr.*dr)*da;
+
+    if (nargout==2)
+        return
+    end
+
+    hess = zeros(numel(dx),2,2);
+
+    hess(:,1,1) = 0;
+    hess(:,1,2) = 0;
+    hess(:,2,1) = 0;
+    hess(:,2,2) = 0;
+
+    eterm = exp(-dr.*dr./(2*dsigt.^2));
+
+    hess(:,1,1) = eterm./(dr.*dr);
+    hess(:,2,2) = eterm./(dr.*dr);
+
+    fact = -(2*dsigt.^2+dr.^2)./(dsigt.^2*dr.^2).*eterm;
+    hess(:,1,1) = hess(:,1,1) + (dx.*dx./dr.^2).*fact;
+    hess(:,1,2) = hess(:,1,2) + (dx.*dy./dr.^2).*fact;
+    hess(:,2,1) = hess(:,2,1) + (dx.*dy./dr.^2).*fact;
+    hess(:,2,2) = hess(:,2,2) + (dy.*dy./dr.^2).*fact;
+
+    hess(:,1,1) = hess(:,1,1)/(2*pi);
+    hess(:,2,1) = hess(:,2,1)/(2*pi);
+    hess(:,1,2) = hess(:,1,2)/(2*pi);
+    hess(:,2,2) = hess(:,2,2)/(2*pi);
+
+    hess(:,1,1) = hess(:,1,1) - 2./(dr.*dr)*da;
+    hess(:,2,2) = hess(:,2,2) - 2./(dr.*dr)*da;
+
+    fact = da*4./dr.^4;
+    hess(:,1,1) = hess(:,1,1) + dx.*dx.*fact;
+    hess(:,1,2) = hess(:,1,2) + dx.*dy.*fact;
+    hess(:,2,1) = hess(:,2,1) + dx.*dy.*fact;
+    hess(:,2,2) = hess(:,2,2) + dy.*dy.*fact;
+
+    if (nargout==3)
+        return
+    end
+
+
+    fact = eterm./(dsigt.^3)/(2*pi);
+    hess_sig(:,1) = -dx*fact;
+    hess_sig(:,2) = -dy*fact;
+
+end
diff --git a/chunkie/+chnk/+smoother/gpsi_loc.m b/chunkie/+chnk/+smoother/gpsi_loc.m
new file mode 100644
index 00000000..7a3ea8d7
--- /dev/null
+++ b/chunkie/+chnk/+smoother/gpsi_loc.m
@@ -0,0 +1,48 @@
+function [val,grad,hess,hess_sig] = gpsi_loc(dx,dy,dr,dsigt)
+
+    degam = 0.57721566490153286060651209008240243104;
+    da = sqrt(2)*dsigt;
+    z = dr./da;
+    z2 = z.*z;
+    val = (z2-z2.^2/4+z2.^3/18-z2.^4/96+z2.^5/600-z2.^6/4320);
+    val = val./(4*pi);
+    val = val - degam/(4*pi); 
+
+    if (nargout==1) 
+        return
+    end
+
+    grad = zeros(numel(dx),2);
+    vvd = (1-z2/2+z2.^2/6-z2.^3/24+z2.^4/120-z2.^5/720);
+    grad(:,1) = dx.*vvd./(2*pi*da.^2);
+    grad(:,2) = dy.*vvd./(2*pi*da.^2);
+
+    if (nargout==2) 
+        return
+    end
+
+    hess = zeros(numel(dx),2,2);
+    hess(:,1,1) = -vvd./(2*pi*da.^2);
+    hess(:,2,2) = -vvd./(2*pi*da.^2);
+    hess(:,1,2) = 0;
+    hess(:,2,1) = 0;
+
+    vvvd = (-1+4*z2/6-6*z2.^2/24+8*z2.^3/120-10*z2.^4/720);
+    vvvd = -vvvd./(2*pi*da.^2)./da./da;
+    hess(:,1,1) = hess(:,1,1) + dx.*dx.*vvvd;
+    hess(:,2,2) = hess(:,2,2) + dy.*dy.*vvvd;
+    hess(:,1,2) = hess(:,1,2) + dx.*dy.*vvvd;
+    hess(:,2,1) = hess(:,2,1) + dx.*dy.*vvvd;
+
+    if (nargout==3)
+        return
+    end
+
+    hess_sig = zeros(numel(dx),2);
+
+    eterm = exp(-dr.*dr./(2*dsigt.^2));
+    fact = eterm./(dsigt.^3)/(2*pi);
+    hess_sig(:,1) = -dx.*fact;
+    hess_sig(:,2) = -dy.*fact;
+
+end
diff --git a/chunkie/+chnk/+smoother/green.m b/chunkie/+chnk/+smoother/green.m
new file mode 100644
index 00000000..66a44c5c
--- /dev/null
+++ b/chunkie/+chnk/+smoother/green.m
@@ -0,0 +1,61 @@
+function [val,grad,hess,hess_sig] = green(rs,rt,dsigt)
+
+[~,ns] = size(rs);
+[~,nt] = size(rt);
+
+xs = repmat(rs(1,:),nt,1);
+ys = repmat(rs(2,:),nt,1);
+
+xt = repmat(rt(1,:).',1,ns);
+yt = repmat(rt(2,:).',1,ns);
+dst = repmat(dsigt.',1,ns);
+
+dx = xt-xs;
+dy = yt-ys;
+dr = sqrt(dx.^2+dy.^2);
+
+sz = size(dx);
+dx = dx(:);
+dy = dy(:);
+dr = dr(:);
+dst = dst(:);
+
+da = sqrt(2.0d0)*dsigt;
+z = dr./da;
+
+iloc = find(z < 0.1);
+ifar = find(z >= 0.1);
+
+vs = zeros(ns*nt,1);
+gs = zeros(ns*nt,2);
+hs = zeros(ns*nt,2,2);
+hsigs = zeros(ns*nt,2);
+
+if (numel(iloc) ~=0)
+    [vl,gl,hl,hsl] = chnk.smoother.gpsi_loc(dx,dy,dr,dst);
+    vs(iloc) = vl;
+    gs(iloc,:) = gl;
+    hs(iloc,:,:) = hl;
+    hsigs(iloc,:) = hsl;
+end
+if (numel(ifar) ~=0)
+    [vf,gf,hf,hsf] = chnk.smoother.gpsi_all(dx,dy,dr,dsigt);
+    vs(ifar) = vf;
+    gs(ifar,:) = gf;
+    hs(ifar,:,:) = hf;
+    hsigs(ifar,:) = hsf;
+end
+
+val = reshape(vs,sz);
+if (nargout > 1)
+    grad = squeeze(reshape(gs,[sz,2]));
+end
+if (nargout > 2)
+    hess = squeeze(reshape(hs,[sz,2,2]));
+end
+if (nargout > 3)
+    hess_sig = squeeze(reshape(hsigs,[sz,2]));
+
+end
+
+end
diff --git a/chunkie/+chnk/+smoother/psi_eval.m b/chunkie/+chnk/+smoother/psi_eval.m
new file mode 100644
index 00000000..2148005d
--- /dev/null
+++ b/chunkie/+chnk/+smoother/psi_eval.m
@@ -0,0 +1,6 @@
+function [val,vald] = psi_eval(x,dsig)
+    z = x.*x./(2*dsig.*dsig);
+    vexp = chnk.smoother.expeval(z,1);
+    val  = vexp/(4*pi);
+    vald = -exp(-z)./(2*pi*x);
+end

From 8ea18a58fd45e5fb4ad2006d20fb10ace02ec4d0 Mon Sep 17 00:00:00 2001
From: Manas Rachh <mrachh@flatironinstitute.org>
Date: Fri, 18 Oct 2024 16:15:17 -0700
Subject: [PATCH 03/15] curve discretizer

---
 chunkie/+chnk/+smoother/get_pseudo_normals.m |  6 ----
 chunkie/+chnk/+smoother/get_umesh.m          | 23 +++++++++++++++
 chunkie/+chnk/+smoother/smooth_curve.m       | 31 +++++++++++++++-----
 3 files changed, 46 insertions(+), 14 deletions(-)
 delete mode 100644 chunkie/+chnk/+smoother/get_pseudo_normals.m
 create mode 100644 chunkie/+chnk/+smoother/get_umesh.m

diff --git a/chunkie/+chnk/+smoother/get_pseudo_normals.m b/chunkie/+chnk/+smoother/get_pseudo_normals.m
deleted file mode 100644
index e1110d61..00000000
--- a/chunkie/+chnk/+smoother/get_pseudo_normals.m
+++ /dev/null
@@ -1,6 +0,0 @@
-function pseudo_normals = get_pseudo_normals(verts)
-    [~, nv] = size(verts);
-    verts_ext = [verts, verts(:,1)];
-    d = verts_ext(:,2:end) - verts_ext(:,1:nv);
-    
-end
\ No newline at end of file
diff --git a/chunkie/+chnk/+smoother/get_umesh.m b/chunkie/+chnk/+smoother/get_umesh.m
new file mode 100644
index 00000000..3808cc6d
--- /dev/null
+++ b/chunkie/+chnk/+smoother/get_umesh.m
@@ -0,0 +1,23 @@
+function [umesh] = get_umesh(verts)
+    umesh = [];
+    umesh.verts = verts;
+    [~, nv] = size(verts);
+    verts_ext = [verts, verts(:,1)];
+    d = verts_ext(:,2:end) - verts_ext(:,1:nv);
+    umesh.centroids = 0.5*(verts_ext(:,2:end) + verts_ext(:,1:nv));
+    dd = sqrt(sum(d.^2, 1));
+    umesh.lengths = dd;
+    face_normals = zeros(2,nv);
+    face_normals(1,:) = d(2,:)./dd;
+    face_normals(2,:) = -d(1,:)./dd;
+
+    umesh.face_normals = face_normals;
+    face_normals_ext = [face_normals(:,end), face_normals];
+    pseudo_normals = 0.5*(face_normals_ext(:,2:end) + ...
+               face_normals_ext(:,1:nv));
+    dd = sqrt(sum(pseudo_normals.^2, 1));
+    pseudo_normals = pseudo_normals./dd;
+
+    umesh.pseudo_normals = pseudo_normals;
+    
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+smoother/smooth_curve.m b/chunkie/+chnk/+smoother/smooth_curve.m
index 05827387..fe314f2e 100644
--- a/chunkie/+chnk/+smoother/smooth_curve.m
+++ b/chunkie/+chnk/+smoother/smooth_curve.m
@@ -1,12 +1,27 @@
-nv = 10;
-z = exp(1j*2*pi*(1:nv-1)/nv);
+nv = 3;
+z = exp(1j*2*pi*(1:nv)/nv);
 verts = [real(z); imag(z)];
 nchs = randi([1,10], 1, nv);
+k = 16;
 
-umesh = [];
-umesh.verts = verts;
-umesh.pseudo_normals = get_pseudo_normal(verts);
+kquad = 32;
 
-qmesh = [];
-qmesh
-dmesh = [];
+
+
+umesh = chnk.smoother.get_umesh(verts);
+dmesh = chnk.smoother.get_mesh(umesh, nchs, k);
+qmesh = get_mesh(umesh, nchs, kquad);
+
+figure(1)
+clf
+plot(verts(1,:), verts(2,:), 'k.','MarkerSize',20); hold on;
+plot(umesh.centroids(1,:), umesh.centroids(2,:), 'r.', 'MarkerSize',20);
+quiver(verts(1,:), verts(2,:), umesh.pseudo_normals(1,:), umesh.pseudo_normals(2,:),'k')
+quiver(umesh.centroids(1,:), umesh.centroids(2,:), umesh.face_normals(1,:), umesh.face_normals(2,:),'r')
+
+figure(2)
+clf
+plot(dmesh.r(1,:), dmesh.r(2,:),'k.'); hold on;
+axis equal;
+% quiver(dmesh.r(1,:), dmesh.r(2,:), dmesh.n(1,:), dmesh.n(2,:),'k');
+quiver(dmesh.r(1,:), dmesh.r(2,:), dmesh.pseudo_normals(1,:), dmesh.pseudo_normals(2,:),'r');
\ No newline at end of file

From 839b7efbc8c9b64423d3364ceea4ae90f245c19f Mon Sep 17 00:00:00 2001
From: jghoskins <32715900+jghoskins@users.noreply.github.com>
Date: Fri, 18 Oct 2024 16:19:38 -0700
Subject: [PATCH 04/15] Create get_sigs.m

---
 chunkie/+chnk/+smoother/get_sigs.m | 40 ++++++++++++++++++++++++++++++
 1 file changed, 40 insertions(+)
 create mode 100644 chunkie/+chnk/+smoother/get_sigs.m

diff --git a/chunkie/+chnk/+smoother/get_sigs.m b/chunkie/+chnk/+smoother/get_sigs.m
new file mode 100644
index 00000000..99672d93
--- /dev/null
+++ b/chunkie/+chnk/+smoother/get_sigs.m
@@ -0,0 +1,40 @@
+function [val,grad] = get_sigs(src_input,rt,sig0,dlam)
+
+        rc = src_input.centroids;
+        dlens = src_input.lengths;
+
+        [~,ns] = size(rc);
+        [~,nt] = size(rt);
+
+        xc = repmat(rc(1,:),nt,1);
+        yc = repmat(rc(2,:),nt,1);
+        dl = repmat(dlens,nt,1);
+
+        xt = repmat(rt(1,:).',1,ns);
+        yt = repmat(rt(2,:).',1,ns);
+
+        dx = xt-xc;
+        dy = yt-yc;
+        dr = sqrt(dx.^2+dy.^2);
+        dmin = min(dr,[],2);
+
+        dsigj = dl/dlam;
+        dexpn = exp(-(dr.^2-dmin.^2)./(2*sig0^2));
+        dsignum = sum(dsigj.*dexpn,2);
+        dsigden = sum(dexpn,2);
+        dsignumx = -dsigj.*dexpn.*(dx./sig0^2);
+        dsignumy = -dsigj.*dexpn.*(dy./sig0^2);
+        dsig = dsignum./dsigden;
+        dsigx = dsignumx./dsigden;
+        dsigy = dsignumy./dsigden;
+
+        val = dsig;
+
+        sz = size(dsig);
+        grad = zeros([sz,2]);
+        grad(:,:,1) = dsigx;
+        grad(:,:,2) = dsigy;
+        
+        grad = squeeze(grad);
+
+end
\ No newline at end of file

From a161d987f0094b2a7696b4d90179be2bba0becb9 Mon Sep 17 00:00:00 2001
From: jghoskins <32715900+jghoskins@users.noreply.github.com>
Date: Fri, 18 Oct 2024 16:25:25 -0700
Subject: [PATCH 05/15] Update get_sigs.m

---
 chunkie/+chnk/+smoother/get_sigs.m | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/chunkie/+chnk/+smoother/get_sigs.m b/chunkie/+chnk/+smoother/get_sigs.m
index 99672d93..a0493a13 100644
--- a/chunkie/+chnk/+smoother/get_sigs.m
+++ b/chunkie/+chnk/+smoother/get_sigs.m
@@ -22,8 +22,8 @@
         dexpn = exp(-(dr.^2-dmin.^2)./(2*sig0^2));
         dsignum = sum(dsigj.*dexpn,2);
         dsigden = sum(dexpn,2);
-        dsignumx = -dsigj.*dexpn.*(dx./sig0^2);
-        dsignumy = -dsigj.*dexpn.*(dy./sig0^2);
+        dsignumx = -sum(dsigj.*dexpn.*(dx./sig0^2),2);
+        dsignumy = -sum(dsigj.*dexpn.*(dy./sig0^2),2);
         dsig = dsignum./dsigden;
         dsigx = dsignumx./dsigden;
         dsigy = dsignumy./dsigden;

From 5df96f25003d1db76d84a77e6ee40d03d245f9e3 Mon Sep 17 00:00:00 2001
From: Manas Rachh <mrachh@flatironinstitute.org>
Date: Fri, 18 Oct 2024 22:28:30 -0700
Subject: [PATCH 06/15] updating smoother

---
 chunkie/+chnk/+smoother/get_sigs.m     |   6 +-
 chunkie/+chnk/+smoother/gpsi_all.m     |   6 +-
 chunkie/+chnk/+smoother/green.m        | 117 +++++++++++++------------
 chunkie/+chnk/+smoother/smooth_curve.m |  65 +++++++++++++-
 4 files changed, 126 insertions(+), 68 deletions(-)

diff --git a/chunkie/+chnk/+smoother/get_sigs.m b/chunkie/+chnk/+smoother/get_sigs.m
index a0493a13..dcdb13fb 100644
--- a/chunkie/+chnk/+smoother/get_sigs.m
+++ b/chunkie/+chnk/+smoother/get_sigs.m
@@ -1,7 +1,7 @@
-function [val,grad] = get_sigs(src_input,rt,sig0,dlam)
+function [val,grad] = get_sigs(umesh,rt,sig0,dlam)
 
-        rc = src_input.centroids;
-        dlens = src_input.lengths;
+        rc = umesh.centroids;
+        dlens = umesh.lengths;
 
         [~,ns] = size(rc);
         [~,nt] = size(rt);
diff --git a/chunkie/+chnk/+smoother/gpsi_all.m b/chunkie/+chnk/+smoother/gpsi_all.m
index 37e69076..9b9e4eb3 100644
--- a/chunkie/+chnk/+smoother/gpsi_all.m
+++ b/chunkie/+chnk/+smoother/gpsi_all.m
@@ -33,7 +33,7 @@
     hess(:,1,1) = eterm./(dr.*dr);
     hess(:,2,2) = eterm./(dr.*dr);
 
-    fact = -(2*dsigt.^2+dr.^2)./(dsigt.^2*dr.^2).*eterm;
+    fact = -(2*dsigt.^2+dr.^2)./(dsigt.^2.*dr.^2).*eterm;
     hess(:,1,1) = hess(:,1,1) + (dx.*dx./dr.^2).*fact;
     hess(:,1,2) = hess(:,1,2) + (dx.*dy./dr.^2).*fact;
     hess(:,2,1) = hess(:,2,1) + (dx.*dy./dr.^2).*fact;
@@ -59,7 +59,7 @@
 
 
     fact = eterm./(dsigt.^3)/(2*pi);
-    hess_sig(:,1) = -dx*fact;
-    hess_sig(:,2) = -dy*fact;
+    hess_sig(:,1) = -dx.*fact;
+    hess_sig(:,2) = -dy.*fact;
 
 end
diff --git a/chunkie/+chnk/+smoother/green.m b/chunkie/+chnk/+smoother/green.m
index 66a44c5c..83736972 100644
--- a/chunkie/+chnk/+smoother/green.m
+++ b/chunkie/+chnk/+smoother/green.m
@@ -1,61 +1,62 @@
 function [val,grad,hess,hess_sig] = green(rs,rt,dsigt)
-
-[~,ns] = size(rs);
-[~,nt] = size(rt);
-
-xs = repmat(rs(1,:),nt,1);
-ys = repmat(rs(2,:),nt,1);
-
-xt = repmat(rt(1,:).',1,ns);
-yt = repmat(rt(2,:).',1,ns);
-dst = repmat(dsigt.',1,ns);
-
-dx = xt-xs;
-dy = yt-ys;
-dr = sqrt(dx.^2+dy.^2);
-
-sz = size(dx);
-dx = dx(:);
-dy = dy(:);
-dr = dr(:);
-dst = dst(:);
-
-da = sqrt(2.0d0)*dsigt;
-z = dr./da;
-
-iloc = find(z < 0.1);
-ifar = find(z >= 0.1);
-
-vs = zeros(ns*nt,1);
-gs = zeros(ns*nt,2);
-hs = zeros(ns*nt,2,2);
-hsigs = zeros(ns*nt,2);
-
-if (numel(iloc) ~=0)
-    [vl,gl,hl,hsl] = chnk.smoother.gpsi_loc(dx,dy,dr,dst);
-    vs(iloc) = vl;
-    gs(iloc,:) = gl;
-    hs(iloc,:,:) = hl;
-    hsigs(iloc,:) = hsl;
-end
-if (numel(ifar) ~=0)
-    [vf,gf,hf,hsf] = chnk.smoother.gpsi_all(dx,dy,dr,dsigt);
-    vs(ifar) = vf;
-    gs(ifar,:) = gf;
-    hs(ifar,:,:) = hf;
-    hsigs(ifar,:) = hsf;
-end
-
-val = reshape(vs,sz);
-if (nargout > 1)
-    grad = squeeze(reshape(gs,[sz,2]));
-end
-if (nargout > 2)
-    hess = squeeze(reshape(hs,[sz,2,2]));
-end
-if (nargout > 3)
-    hess_sig = squeeze(reshape(hsigs,[sz,2]));
-
-end
+    
+    [~,ns] = size(rs);
+    [~,nt] = size(rt);
+    
+    xs = repmat(rs(1,:),nt,1);
+    ys = repmat(rs(2,:),nt,1);
+    
+    xt = repmat(rt(1,:).',1,ns);
+    yt = repmat(rt(2,:).',1,ns);
+    dst = repmat(dsigt.',1,ns);
+    
+    dx = xt-xs;
+    dy = yt-ys;
+    dr = sqrt(dx.^2+dy.^2);
+    
+    sz = size(dx);
+    dx = dx(:);
+    dy = dy(:);
+    dr = dr(:);
+    dst = dst(:);
+    
+    da = sqrt(2.0d0)*dst;
+    z = dr./da;
+    
+    iloc = find(z < 0.1);
+    ifar = find(z >= 0.1);
+    
+    vs = zeros(ns*nt,1);
+    gs = zeros(ns*nt,2);
+    hs = zeros(ns*nt,2,2);
+    hsigs = zeros(ns*nt,2);
+    
+    if (numel(iloc) ~=0)
+        [vl,gl,hl,hsl] = chnk.smoother.gpsi_loc(dx(iloc), dy(iloc), ...
+                dr(iloc), dst(iloc));
+        vs(iloc) = vl;
+        gs(iloc,:) = gl;
+        hs(iloc,:,:) = hl;
+        hsigs(iloc,:) = hsl;
+    end
+    if (numel(ifar) ~=0)
+        [vf,gf,hf,hsf] = chnk.smoother.gpsi_all(dx(ifar), dy(ifar), ...
+              dr(ifar), dst(ifar));
+        vs(ifar) = vf;
+        gs(ifar,:) = gf;
+        hs(ifar,:,:) = hf;
+        hsigs(ifar,:) = hsf;
+    end
+    
+    val = reshape(vs,sz);
+    if (nargout > 1)
+        grad = squeeze(reshape(gs,[sz,2]));
+    end
+    if (nargout > 2)
+        hess = squeeze(reshape(hs,[sz,2,2]));
+    end
+    if (nargout > 3)
+        hess_sig = squeeze(reshape(hsigs,[sz,2]));
+    end
 
 end
diff --git a/chunkie/+chnk/+smoother/smooth_curve.m b/chunkie/+chnk/+smoother/smooth_curve.m
index fe314f2e..d92d740a 100644
--- a/chunkie/+chnk/+smoother/smooth_curve.m
+++ b/chunkie/+chnk/+smoother/smooth_curve.m
@@ -1,7 +1,7 @@
 nv = 3;
 z = exp(1j*2*pi*(1:nv)/nv);
 verts = [real(z); imag(z)];
-nchs = randi([1,10], 1, nv);
+nchs = 3*ones(nv,1);
 k = 16;
 
 kquad = 32;
@@ -10,7 +10,7 @@
 
 umesh = chnk.smoother.get_umesh(verts);
 dmesh = chnk.smoother.get_mesh(umesh, nchs, k);
-qmesh = get_mesh(umesh, nchs, kquad);
+qmesh = chnk.smoother.get_mesh(umesh, nchs, kquad);
 
 figure(1)
 clf
@@ -23,5 +23,62 @@
 clf
 plot(dmesh.r(1,:), dmesh.r(2,:),'k.'); hold on;
 axis equal;
-% quiver(dmesh.r(1,:), dmesh.r(2,:), dmesh.n(1,:), dmesh.n(2,:),'k');
-quiver(dmesh.r(1,:), dmesh.r(2,:), dmesh.pseudo_normals(1,:), dmesh.pseudo_normals(2,:),'r');
\ No newline at end of file
+quiver(dmesh.r(1,:), dmesh.r(2,:), dmesh.pseudo_normals(1,:), dmesh.pseudo_normals(2,:),'r');
+
+
+[~, nd] = size(dmesh.r);
+h = zeros(nd,1);
+n_newton = 5;
+
+sig0 = sqrt(5)*max(umesh.lengths);
+lam = 10;
+
+
+ww = qmesh.wts(:).';
+dpx = dmesh.pseudo_normals(1,:).';
+dpy = dmesh.pseudo_normals(2,:).';
+
+rnx = qmesh.n(1,:).';
+rny = qmesh.n(2,:).';
+
+for i=1:n_newton
+    rt = dmesh.r;
+    rt(1,:) = rt(1,:) + (h.*dpx).';
+    rt(2,:) = rt(2,:) + (h.*dpy).';
+
+    [sig, sig_grad] = chnk.smoother.get_sigs(umesh, rt, sig0, lam);
+    [val, grad, hess, hess_sig] = chnk.smoother.green(qmesh.r, rt, sig);
+    gx = grad(:,:,1); 
+    gy = grad(:,:,2);
+    
+    phi = -(gx*(rnx.*qmesh.wts) + gy*(rny.*qmesh.wts)) - 0.5;
+    
+    
+    h11 = hess(:,:,1,1).*ww;
+    h12 = hess(:,:,1,2).*ww;
+    h21 = hess(:,:,2,1).*ww;
+    h22 = hess(:,:,2,2).*ww;
+    
+    h1sig = hess_sig(:,:,1).*ww;
+    h2sig = hess_sig(:,:,2).*ww;
+
+    dx1 = h11*rnx + h12*rny;
+    dy1 = h21*rnx + h22*rny;
+    
+    dx2 = h1sig*rnx + h2sig*rny;
+    dy2 = dx2;
+    
+    dsigx = sig_grad(:,1);
+    dsigy = sig_grad(:,2);
+
+    dx2 = dx2.*dsigx;
+    dy2 = dy2.*dsigy;
+    
+    
+    dphidh = (dx1 + dx2).*dpx + (dy1 + dy2).*dpy;
+    h = h - phi./dphidh;
+
+    figure(i)
+    plot(rt(1,:), rt(2,:), 'k.')
+ 
+end
\ No newline at end of file

From a53812ab12651a7e981c14c2dd391fd54d37d45f Mon Sep 17 00:00:00 2001
From: Manas Rachh <mrachh@flatironinstitute.org>
Date: Sat, 19 Oct 2024 09:14:50 -0700
Subject: [PATCH 07/15] fixing get mesh depedency

---
 chunkie/+chnk/+smoother/get_mesh.m | 40 ++++++++++++++++++++++++++++++
 1 file changed, 40 insertions(+)
 create mode 100644 chunkie/+chnk/+smoother/get_mesh.m

diff --git a/chunkie/+chnk/+smoother/get_mesh.m b/chunkie/+chnk/+smoother/get_mesh.m
new file mode 100644
index 00000000..44346fbe
--- /dev/null
+++ b/chunkie/+chnk/+smoother/get_mesh.m
@@ -0,0 +1,40 @@
+function qmesh = get_mesh(umesh, nchs, k)
+    [x, w] = lege.exps(k);
+    
+    nchtot = sum(nchs(:));
+    n = nchtot*k;
+    qmesh.r = zeros(2,n);
+    qmesh.n = zeros(2,n);
+    qmesh.pseudo_normals = zeros(2,n);
+    qmesh.wts = zeros(n,1);
+    
+    [~, ne] = size(umesh.verts);
+    verts_ext = [umesh.verts, umesh.verts(:,1)];
+    pseudo_normals_ext = [umesh.pseudo_normals, umesh.pseudo_normals(:,1)];
+    istart = 1;
+    for i=1:ne
+        [xext, wext] = get_nodes(x, w, nchs(i));
+        ll = length(xext);
+        iend = istart+ll-1;
+        qmesh.r(:,istart:iend) = verts_ext(:,i) + ...
+             (verts_ext(:,i+1) - verts_ext(:,i)).*xext.';
+        qmesh.n(:,istart:iend) = repmat(umesh.face_normals(:,i),[1,ll]);
+        qmesh.pseudo_normals(:,istart:iend) = ...
+             repmat(umesh.pseudo_normals(:,i), [1,ll]) + ...
+             (pseudo_normals_ext(:,i+1) - pseudo_normals_ext(:,i)).*xext.';
+        qmesh.wts(istart:iend) = wext.*umesh.lengths(i);
+        istart = istart + ll;
+    end
+
+end
+
+
+function [xext, wext] = get_nodes(x, w, nch)
+    tabs = 0:1/nch:1;
+    tas = tabs(1:nch);
+    tbs = tabs(2:end);
+    xext = (x(:)+1)/2.*(tbs - tas) + tas;
+    xext = xext(:);
+    wext = w(:)/2.*(tbs-tas);
+    wext = wext(:);    
+end
\ No newline at end of file

From d5ff3cee26145262ac0b7c2c77f93fa77e4ed27c Mon Sep 17 00:00:00 2001
From: jghoskins <32715900+jghoskins@users.noreply.github.com>
Date: Sat, 19 Oct 2024 11:15:17 -0700
Subject: [PATCH 08/15] smoothed smoother

Added newt_step and some support for multiple regions.
---
 chunkie/+chnk/+smoother/merge.m         | 102 ++++++++++++++++++++++++
 chunkie/+chnk/+smoother/newt_step.m     |  63 +++++++++++++++
 chunkie/+chnk/+smoother/smooth_curve.m  |   4 +-
 chunkie/+chnk/+smoother/smooth_curve2.m |  88 ++++++++++++++++++++
 chunkie/+chnk/+smoother/smooth_curve3.m |  87 ++++++++++++++++++++
 5 files changed, 342 insertions(+), 2 deletions(-)
 create mode 100644 chunkie/+chnk/+smoother/merge.m
 create mode 100644 chunkie/+chnk/+smoother/newt_step.m
 create mode 100644 chunkie/+chnk/+smoother/smooth_curve2.m
 create mode 100644 chunkie/+chnk/+smoother/smooth_curve3.m

diff --git a/chunkie/+chnk/+smoother/merge.m b/chunkie/+chnk/+smoother/merge.m
new file mode 100644
index 00000000..5c217cd0
--- /dev/null
+++ b/chunkie/+chnk/+smoother/merge.m
@@ -0,0 +1,102 @@
+function [umesh,dmesh,qmesh,scales,levels] = merge(umeshes,dmeshes,...
+    qmeshes,src_codes,targ_codes)
+
+    nv = 0;
+    for ii=1:numel(umeshes)
+        nv = nv + size(umeshes{ii}.verts,2);
+    end
+
+    verts = zeros(2,nv);
+    centroids = zeros(2,nv);
+    lengths = zeros(1,nv);
+    face_normals = zeros(2,nv);
+    pseudo_normals = zeros(2,nv);
+
+    ind = 1;
+    for ii=1:numel(umeshes)
+        nv_loc = size(umeshes{ii}.verts,2);
+        ind2 = ind + nv_loc - 1;
+        verts(:,ind:ind2)  = umeshes{ii}.verts;
+        centroids(:,ind:ind2)  = umeshes{ii}.centroids;
+        lengths(:,ind:ind2) = umeshes{ii}.lengths;
+        face_normals(:,ind:ind2)  = umeshes{ii}.face_normals;
+        pseudo_normals(:,ind:ind2)  = umeshes{ii}.pseudo_normals;
+        ind = ind2 + 1;
+    end
+
+    umesh = [];
+    umesh.verts = verts;
+    umesh.centroids = centroids;
+    umesh.lengths = lengths;
+    umesh.face_normals = face_normals;
+    umesh.pseudo_normals = pseudo_normals;
+
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    npd = 0;
+    for ii=1:numel(dmeshes)
+        npd = npd + size(dmeshes{ii}.r,2);
+    end
+
+    r = zeros(2,npd);
+    n = zeros(2,npd);
+    pseudo_normals = zeros(2,npd);
+    wts = zeros(npd,1);
+    levels = ones(npd,1);
+    
+    ind = 1;
+    for ii=1:numel(dmeshes)
+        np_loc = size(dmeshes{ii}.r,2);
+        ind2 = ind + np_loc - 1;
+        r(:,ind:ind2) = dmeshes{ii}.r;
+        n(:,ind:ind2) = dmeshes{ii}.n;
+        pseudo_normals(:,ind:ind2) = dmeshes{ii}.pseudo_normals;
+        wts(ind:ind2,1) = dmeshes{ii}.wts;
+        if (nargin>4)
+            levels(ind:ind2) = targ_codes(ii);
+        end
+        ind = ind2 + 1;
+    end
+
+    dmesh = [];
+    dmesh.r = r;
+    dmesh.n = n;
+    dmesh.pseudo_normals = pseudo_normals;
+    dmesh.wts = wts;
+
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    nqd = 0;
+    for ii=1:numel(qmeshes)
+        nqd = nqd + size(qmeshes{ii}.r,2);
+    end
+
+    r = zeros(2,nqd);
+    n = zeros(2,nqd);
+    pseudo_normals = zeros(2,nqd);
+    wts = zeros(nqd,1);
+    scales = ones(nqd,1);
+    
+    ind = 1;
+    for ii=1:numel(qmeshes)
+        np_loc = size(qmeshes{ii}.r,2);
+        ind2 = ind + np_loc - 1;
+        r(:,ind:ind2) = qmeshes{ii}.r;
+        n(:,ind:ind2) = qmeshes{ii}.n;
+        pseudo_normals(:,ind:ind2) = qmeshes{ii}.pseudo_normals;
+        wts(ind:ind2,1) = qmeshes{ii}.wts;
+        if (nargin>3)
+            scales(ind:ind2) = src_codes(ii);
+        end
+        ind = ind2 + 1;
+    end
+
+    qmesh = [];
+    qmesh.r = r;
+    qmesh.n = n;
+    qmesh.pseudo_normals = pseudo_normals;
+    qmesh.wts = wts;
+
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+smoother/newt_step.m b/chunkie/+chnk/+smoother/newt_step.m
new file mode 100644
index 00000000..368abe11
--- /dev/null
+++ b/chunkie/+chnk/+smoother/newt_step.m
@@ -0,0 +1,63 @@
+function [h] = newt_step(h,umesh,dmesh,qmesh,sig0,lam,opts)
+
+level = 0.5;
+step_fact = 1;
+ww = qmesh.wts(:).';
+
+if (nargin > 6)
+if (isfield(opts,"levels"))
+    level = opts.levels;
+end
+
+if (isfield(opts,"scales"))
+    ww = ww.*opts.scales.';
+end
+
+if (isfield(opts,"step_fact"))
+    step_fact = opts.step_fact;
+end
+end
+
+dpx = dmesh.pseudo_normals(1,:).';
+dpy = dmesh.pseudo_normals(2,:).';
+
+rnx = qmesh.n(1,:).';
+rny = qmesh.n(2,:).';
+
+    rt = dmesh.r;
+    rt(1,:) = rt(1,:) + (h.*dpx).';
+    rt(2,:) = rt(2,:) + (h.*dpy).';
+
+    [sig, sig_grad] = chnk.smoother.get_sigs(umesh, rt, sig0, lam);
+    [val, grad, hess, hess_sig] = chnk.smoother.green(qmesh.r, rt, sig);
+    gx = grad(:,:,1); 
+    gy = grad(:,:,2);
+    
+    phi = -(gx*(rnx.*ww.') + gy*(rny.*ww.')) - level;
+    
+    
+    h11 = hess(:,:,1,1).*ww;
+    h12 = hess(:,:,1,2).*ww;
+    h21 = hess(:,:,2,1).*ww;
+    h22 = hess(:,:,2,2).*ww;
+    
+    h1sig = hess_sig(:,:,1).*ww;
+    h2sig = hess_sig(:,:,2).*ww;
+
+    dx1 = h11*rnx + h12*rny;
+    dy1 = h21*rnx + h22*rny;
+    
+    dx2 = h1sig*rnx + h2sig*rny;
+    dy2 = dx2;
+    
+    dsigx = sig_grad(:,1);
+    dsigy = sig_grad(:,2);
+
+    dx2 = dx2.*dsigx;
+    dy2 = dy2.*dsigy;
+    
+    
+    dphidh = (dx1 + dx2).*dpx + (dy1 + dy2).*dpy;
+    h = h - step_fact*phi./dphidh;
+
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+smoother/smooth_curve.m b/chunkie/+chnk/+smoother/smooth_curve.m
index d92d740a..2fc92322 100644
--- a/chunkie/+chnk/+smoother/smooth_curve.m
+++ b/chunkie/+chnk/+smoother/smooth_curve.m
@@ -28,10 +28,10 @@
 
 [~, nd] = size(dmesh.r);
 h = zeros(nd,1);
-n_newton = 5;
+n_newton = 200;
 
 sig0 = sqrt(5)*max(umesh.lengths);
-lam = 10;
+lam = 2.5;
 
 
 ww = qmesh.wts(:).';
diff --git a/chunkie/+chnk/+smoother/smooth_curve2.m b/chunkie/+chnk/+smoother/smooth_curve2.m
new file mode 100644
index 00000000..62795253
--- /dev/null
+++ b/chunkie/+chnk/+smoother/smooth_curve2.m
@@ -0,0 +1,88 @@
+nv = 3;
+z = exp(1j*2*pi*(1:nv)/nv);
+verts = [real(z); imag(z)];
+nchs = 3*ones(nv,1);
+k = 16;
+
+kquad = 32;
+
+
+
+umesh = chnk.smoother.get_umesh(verts);
+dmesh = chnk.smoother.get_mesh(umesh, nchs, k);
+qmesh = chnk.smoother.get_mesh(umesh, nchs, kquad);
+
+figure(1)
+clf
+plot(verts(1,:), verts(2,:), 'k.','MarkerSize',20); hold on;
+plot(umesh.centroids(1,:), umesh.centroids(2,:), 'r.', 'MarkerSize',20);
+quiver(verts(1,:), verts(2,:), umesh.pseudo_normals(1,:), umesh.pseudo_normals(2,:),'k')
+quiver(umesh.centroids(1,:), umesh.centroids(2,:), umesh.face_normals(1,:), umesh.face_normals(2,:),'r')
+
+figure(2)
+clf
+plot(dmesh.r(1,:), dmesh.r(2,:),'k.'); hold on;
+axis equal;
+quiver(dmesh.r(1,:), dmesh.r(2,:), dmesh.pseudo_normals(1,:), dmesh.pseudo_normals(2,:),'r');
+
+
+[~, nd] = size(dmesh.r);
+h = zeros(nd,1);
+n_newton = 3;
+
+sig0 = sqrt(5)*max(umesh.lengths);
+lam = 10;
+
+
+ww = qmesh.wts(:).';
+dpx = dmesh.pseudo_normals(1,:).';
+dpy = dmesh.pseudo_normals(2,:).';
+
+rnx = qmesh.n(1,:).';
+rny = qmesh.n(2,:).';
+
+for i=1:n_newton
+    [h] = chnk.smoother.newt_step(h,umesh,dmesh,qmesh,sig0,lam);
+    % rt = dmesh.r;
+    % rt(1,:) = rt(1,:) + (h.*dpx).';
+    % rt(2,:) = rt(2,:) + (h.*dpy).';
+    % 
+    % [sig, sig_grad] = chnk.smoother.get_sigs(umesh, rt, sig0, lam);
+    % [val, grad, hess, hess_sig] = chnk.smoother.green(qmesh.r, rt, sig);
+    % gx = grad(:,:,1); 
+    % gy = grad(:,:,2);
+    % 
+    % phi = -(gx*(rnx.*qmesh.wts) + gy*(rny.*qmesh.wts)) - 0.5;
+    % 
+    % 
+    % h11 = hess(:,:,1,1).*ww;
+    % h12 = hess(:,:,1,2).*ww;
+    % h21 = hess(:,:,2,1).*ww;
+    % h22 = hess(:,:,2,2).*ww;
+    % 
+    % h1sig = hess_sig(:,:,1).*ww;
+    % h2sig = hess_sig(:,:,2).*ww;
+    % 
+    % dx1 = h11*rnx + h12*rny;
+    % dy1 = h21*rnx + h22*rny;
+    % 
+    % dx2 = h1sig*rnx + h2sig*rny;
+    % dy2 = dx2;
+    % 
+    % dsigx = sig_grad(:,1);
+    % dsigy = sig_grad(:,2);
+    % 
+    % dx2 = dx2.*dsigx;
+    % dy2 = dy2.*dsigy;
+    % 
+    % 
+    % dphidh = (dx1 + dx2).*dpx + (dy1 + dy2).*dpy;
+    % h = h - phi./dphidh;
+
+    figure(i)    
+    rt = dmesh.r;
+    rt(1,:) = rt(1,:) + (h.*dpx).';
+    rt(2,:) = rt(2,:) + (h.*dpy).';
+    plot(rt(1,:), rt(2,:), 'k.')
+ 
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+smoother/smooth_curve3.m b/chunkie/+chnk/+smoother/smooth_curve3.m
new file mode 100644
index 00000000..e1435051
--- /dev/null
+++ b/chunkie/+chnk/+smoother/smooth_curve3.m
@@ -0,0 +1,87 @@
+nv = 3;
+z = exp(1j*2*pi*(1:nv)/nv);
+verts = [real(z); imag(z)];
+nchs = 3*ones(nv,1);
+k = 16;
+
+kquad = 32;
+
+
+
+umesh = chnk.smoother.get_umesh(verts);
+dmesh = chnk.smoother.get_mesh(umesh, nchs, k);
+qmesh = chnk.smoother.get_mesh(umesh, nchs, kquad);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+nv = 3;
+z = exp(1j*2*pi*(1:nv)/nv)+(0.9+1*1i);
+verts2 = [real(z); imag(z)];
+nchs = 3*ones(nv,1);
+k = 16;
+
+kquad = 32;
+
+
+
+umesh2 = chnk.smoother.get_umesh(verts2);
+dmesh2 = chnk.smoother.get_mesh(umesh2, nchs, k);
+qmesh2 = chnk.smoother.get_mesh(umesh2, nchs, kquad);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+umeshes = {umesh,umesh2};
+dmeshes = {dmesh,dmesh2};
+qmeshes = {qmesh,qmesh2};
+src_codes = [1,-1];
+targ_codes = [0.5,0.5]*0;
+[umesh,dmesh,qmesh,scales,levels] = chnk.smoother.merge(umeshes,dmeshes,...
+    qmeshes,src_codes,targ_codes);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+figure(1)
+clf
+verts = umesh.verts;
+plot(verts(1,:), verts(2,:), 'k.','MarkerSize',20); hold on;
+plot(umesh.centroids(1,:), umesh.centroids(2,:), 'r.', 'MarkerSize',20);
+quiver(verts(1,:), verts(2,:), umesh.pseudo_normals(1,:), umesh.pseudo_normals(2,:),'k')
+quiver(umesh.centroids(1,:), umesh.centroids(2,:), umesh.face_normals(1,:), umesh.face_normals(2,:),'r')
+
+figure(2)
+clf
+plot(dmesh.r(1,:), dmesh.r(2,:),'k.'); hold on;
+axis equal;
+quiver(dmesh.r(1,:), dmesh.r(2,:), dmesh.pseudo_normals(1,:), dmesh.pseudo_normals(2,:),'r');
+
+
+[~, nd] = size(dmesh.r);
+h = zeros(nd,1);
+n_newton = 20;
+
+sig0 = sqrt(5)*max(umesh.lengths);
+lam = 10;
+
+
+ww = qmesh.wts(:).';
+dpx = dmesh.pseudo_normals(1,:).';
+dpy = dmesh.pseudo_normals(2,:).';
+
+rnx = qmesh.n(1,:).';
+rny = qmesh.n(2,:).';
+opts = [];
+opts.scales = scales;
+opts.levels  = levels;
+opts.step_fact = 0.9;
+for i=1:n_newton
+    [h] = chnk.smoother.newt_step(h,umesh,dmesh,qmesh,sig0,lam,opts);
+    figure(i)    
+    rt = dmesh.r;
+    rt(1,:) = rt(1,:) + (h.*dpx).';
+    rt(2,:) = rt(2,:) + (h.*dpy).';
+    plot(rt(1,:), rt(2,:), 'k.')
+ 
+end
\ No newline at end of file

From dfd428477ade88ecc1d28a70b742af5c8eb9e33c Mon Sep 17 00:00:00 2001
From: jghoskins <32715900+jghoskins@users.noreply.github.com>
Date: Sat, 19 Apr 2025 18:55:10 -0500
Subject: [PATCH 09/15] Create chnks2chunker.m

Needs polishing, but takes panel coords to a chunker object.
---
 chunkie/chnks2chunker.m | 54 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 54 insertions(+)
 create mode 100644 chunkie/chnks2chunker.m

diff --git a/chunkie/chnks2chunker.m b/chunkie/chnks2chunker.m
new file mode 100644
index 00000000..1f8fe589
--- /dev/null
+++ b/chunkie/chnks2chunker.m
@@ -0,0 +1,54 @@
+function chnkr = chnks2chunker(r)
+
+pref = [];
+
+ifclosed=true;
+dim = size(r,1);
+k   = size(r,2);
+nch = size(r,3); 
+ 
+pref.dim = dim;
+pref.k   = k;
+
+[xs,~,us,vs] = lege.exps(k);
+dermat = (vs*[lege.derpol(us); zeros(1,k)]).';
+
+
+%       . . . start chunking
+
+adjs = zeros(2,nch);
+
+adjs(1,:) = (1:nch) - 1;
+adjs(2,:) = (1:nch) + 1;
+
+adjs(1,1)=nch;
+adjs(2,nch)=1;
+
+
+%       up to here, everything has been done in parameter space, [ta,tb]
+%       . . . finally evaluate the k nodes on each chunk, along with 
+%       derivatives and chunk lengths
+
+chnkr = chunker(pref); % empty chunker
+chnkr = chnkr.addchunk(nch);
+
+
+for i = 1:nch
+    
+    rtmp  = squeeze(r(:,:,i));
+    dtmp  = rtmp*dermat;
+    d2tmp = dtmp*dermat;
+    chnkr.rstor(:,:,i) = rtmp;
+    chnkr.dstor(:,:,i) = dtmp;
+    chnkr.d2stor(:,:,i) = d2tmp;
+end
+
+chnkr.adjstor(:,1:nch) = adjs(:,1:nch);
+
+% Set normals
+chnkr.nstor(:,:,1:nch) = normals(chnkr);
+
+% Set weights
+chnkr.wtsstor(:,1:nch) = weights(chnkr);
+
+end

From 9e1eae26197618fb40aad4baf9fe262a637cb93f Mon Sep 17 00:00:00 2001
From: Manas Rachh <mrachh@flatironinstitute.org>
Date: Mon, 28 Apr 2025 13:53:35 -0400
Subject: [PATCH 10/15] plotting outside for smoother test

---
 chunkie/+chnk/+smoother/smooth_curve2.m | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/chunkie/+chnk/+smoother/smooth_curve2.m b/chunkie/+chnk/+smoother/smooth_curve2.m
index 62795253..47ea7e7d 100644
--- a/chunkie/+chnk/+smoother/smooth_curve2.m
+++ b/chunkie/+chnk/+smoother/smooth_curve2.m
@@ -79,10 +79,13 @@
     % dphidh = (dx1 + dx2).*dpx + (dy1 + dy2).*dpy;
     % h = h - phi./dphidh;
 
-    figure(i)    
+    
     rt = dmesh.r;
     rt(1,:) = rt(1,:) + (h.*dpx).';
     rt(2,:) = rt(2,:) + (h.*dpy).';
-    plot(rt(1,:), rt(2,:), 'k.')
+    
  
-end
\ No newline at end of file
+end
+
+figure(3)    
+plot(rt(1,:), rt(2,:), 'k.')
\ No newline at end of file

From ae01c250a76b3fc1e64122b9b4afd062d18db2ac Mon Sep 17 00:00:00 2001
From: jghoskins <32715900+jghoskins@users.noreply.github.com>
Date: Mon, 28 Apr 2025 14:56:48 -0400
Subject: [PATCH 11/15] wrapping smoother

wrapped smooth routine and added test
---
 chunkie/+chnk/+smoother/newt_step.m           |  8 +-
 chunkie/+chnk/+smoother/smooth.m              | 88 +++++++++++++++++++
 chunkie/+chnk/+smoother/smooth_curve3.m       | 11 ++-
 chunkie/@chunker/chunker.m                    |  1 +
 .../chunkerpoints.m}                          |  2 +-
 devtools/test/smootherTest.m                  |  9 ++
 6 files changed, 115 insertions(+), 4 deletions(-)
 create mode 100644 chunkie/+chnk/+smoother/smooth.m
 rename chunkie/{chnks2chunker.m => @chunker/chunkerpoints.m} (96%)
 create mode 100644 devtools/test/smootherTest.m

diff --git a/chunkie/+chnk/+smoother/newt_step.m b/chunkie/+chnk/+smoother/newt_step.m
index 368abe11..b87a362b 100644
--- a/chunkie/+chnk/+smoother/newt_step.m
+++ b/chunkie/+chnk/+smoother/newt_step.m
@@ -1,4 +1,4 @@
-function [h] = newt_step(h,umesh,dmesh,qmesh,sig0,lam,opts)
+function [h,varargout] = newt_step(h,umesh,dmesh,qmesh,sig0,lam,opts)
 
 level = 0.5;
 step_fact = 1;
@@ -60,4 +60,10 @@
     dphidh = (dx1 + dx2).*dpx + (dy1 + dy2).*dpy;
     h = h - step_fact*phi./dphidh;
 
+    if (nargout > 1)
+        varargout{1} = max(abs(phi(:)));
+    end
+    if (nargout > 2)
+        varargout{2} = abs(phi(:));
+    end
 end
\ No newline at end of file
diff --git a/chunkie/+chnk/+smoother/smooth.m b/chunkie/+chnk/+smoother/smooth.m
new file mode 100644
index 00000000..ce8182f7
--- /dev/null
+++ b/chunkie/+chnk/+smoother/smooth.m
@@ -0,0 +1,88 @@
+function [chnkr,varargout] = smooth(verts,opts)
+
+    nv = size(verts,2);
+
+    k = 16;
+    if isfield(opts,'k')
+       k = opts.k; 
+    end
+    
+    kquad = 32;
+    if isfield(opts,'kquad')
+        kquad = 32;
+    end
+    
+    nch  = 3;
+    nchs = nch*ones(nv,1);
+
+    if isfield(opts,'nchs')
+        if (numel(opts.nchs) == 1)
+            nchs = opts.nchs*ones(nv,1);
+        elseif(size(opts.nchs,2) == nchs)
+            nchs = opts.nchs;
+        end
+    end
+
+    n_newton = 200;
+    if isfield(opts,'n_newton')
+        n_newton = opts.n_newton;
+    end
+
+    etol = 1E-6;
+    if isfield(opts,'etol')
+        etol = opts.etol;
+    end
+
+    lam = 2.5;
+    if isfield(opts,'lam')
+        lam = opts.lam;
+    end
+
+    umesh = chnk.smoother.get_umesh(verts);
+    dmesh = chnk.smoother.get_mesh(umesh, nchs, k);
+    qmesh = chnk.smoother.get_mesh(umesh, nchs, kquad);
+
+    sig0 = sqrt(5)*max(umesh.lengths);
+    if isfield(opts,'sig0')
+        sig0 = opt.sig0;
+    end
+
+    dpx = dmesh.pseudo_normals(1,:).';
+    dpy = dmesh.pseudo_normals(2,:).';
+
+    opts_newt = [];
+    opts_newt.scales = 1;
+    opts_newt.levels = 0.5;
+
+    opts_newt.step_fact = 1;
+    if isfield(opts,'step_fact')
+        opts_newt.step_fact = opts.step_fact;
+    end
+
+    [~, nd] = size(dmesh.r);
+    h = zeros(nd,1);
+
+    for i=1:n_newton
+        [h,err,err_by_pt] = chnk.smoother.newt_step(h,umesh,dmesh,...
+            qmesh,sig0,lam,opts_newt);
+        rt = dmesh.r;
+        rt(1,:) = rt(1,:) + (h.*dpx).';
+        rt(2,:) = rt(2,:) + (h.*dpy).';
+        if (err < etol)
+            break
+        end
+    end
+
+    if (err >= etol)
+        error('ERROR: surface smoother failed to converge');
+        return
+    end
+
+    chnkr = chunker.chunkerpoints(rt);
+    if (nargout > 1)
+        varargout{1} = err;
+    end
+    if (nargout > 2)
+        varargout{2} = err_by_pt;
+    end
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+smoother/smooth_curve3.m b/chunkie/+chnk/+smoother/smooth_curve3.m
index e1435051..1e76359c 100644
--- a/chunkie/+chnk/+smoother/smooth_curve3.m
+++ b/chunkie/+chnk/+smoother/smooth_curve3.m
@@ -76,12 +76,19 @@
 opts.scales = scales;
 opts.levels  = levels;
 opts.step_fact = 0.9;
+
+%%
+
+n_newton = 100;
+etol     = 1E-10;
 for i=1:n_newton
-    [h] = chnk.smoother.newt_step(h,umesh,dmesh,qmesh,sig0,lam,opts);
+    [h,err,err_p_pt] = chnk.smoother.newt_step(h,umesh,dmesh,qmesh,sig0,lam,opts);
     figure(i)    
     rt = dmesh.r;
     rt(1,:) = rt(1,:) + (h.*dpx).';
     rt(2,:) = rt(2,:) + (h.*dpy).';
     plot(rt(1,:), rt(2,:), 'k.')
- 
+    if (err < etol)
+        break
+    end
 end
\ No newline at end of file
diff --git a/chunkie/@chunker/chunker.m b/chunkie/@chunker/chunker.m
index e0e2b35e..3ce79cd4 100644
--- a/chunkie/@chunker/chunker.m
+++ b/chunkie/@chunker/chunker.m
@@ -369,6 +369,7 @@
     methods(Static)
         obj = chunkerfunc(fcurve,varargin)
         obj = chunkerpoly(verts,varargin)
+        obj = chunkerpoints(r,varargin)
         function lvlrfac = lvlrfacdefault()
             lvlrfac = 2.1;
         end
diff --git a/chunkie/chnks2chunker.m b/chunkie/@chunker/chunkerpoints.m
similarity index 96%
rename from chunkie/chnks2chunker.m
rename to chunkie/@chunker/chunkerpoints.m
index 1f8fe589..3b41f0df 100644
--- a/chunkie/chnks2chunker.m
+++ b/chunkie/@chunker/chunkerpoints.m
@@ -1,4 +1,4 @@
-function chnkr = chnks2chunker(r)
+function chnkr = chunkerpoints(r)
 
 pref = [];
 
diff --git a/devtools/test/smootherTest.m b/devtools/test/smootherTest.m
new file mode 100644
index 00000000..2f439e0b
--- /dev/null
+++ b/devtools/test/smootherTest.m
@@ -0,0 +1,9 @@
+nv = 3;
+z = exp(1j*2*pi*(1:nv)/nv);
+verts = [real(z); imag(z)];
+
+opts = [];
+opts.lam = 10;
+[chnkr,err,err_by_pt] = chnk.smoother.smooth(verts,opts);
+
+assert(err<1E-6)

From d0a013ab856764b3a68bfb1f1acf2bdd4889548e Mon Sep 17 00:00:00 2001
From: jghoskins <32715900+jghoskins@users.noreply.github.com>
Date: Mon, 28 Apr 2025 15:24:08 -0400
Subject: [PATCH 12/15] Update chunkerpoints.m

---
 chunkie/@chunker/chunkerpoints.m | 52 ++++++++++++++++++++++++++++----
 1 file changed, 46 insertions(+), 6 deletions(-)

diff --git a/chunkie/@chunker/chunkerpoints.m b/chunkie/@chunker/chunkerpoints.m
index 3b41f0df..c3ff91ed 100644
--- a/chunkie/@chunker/chunkerpoints.m
+++ b/chunkie/@chunker/chunkerpoints.m
@@ -1,8 +1,34 @@
-function chnkr = chunkerpoints(r)
+function chnkr = chunkerpoints(src)
+
+d  = [];
+d2 = [];
+if (strcmp(class(src),'double'))
+    r = src;
+elseif (strcmp(class(src),'struct'))
+    if (isfield(src,'r'))
+        r = src.r;
+    else
+        error('ERROR: missing field r in chunkerpoints');
+    end
+    if (isfield(src,'d') && isequal(size(r),size(src.d)))
+        d = src.d;
+    end
+    if (isfield(src,'d2') && isequal(size(r),size(src.d)))
+        d2 = src.d2;
+    end
+end
+
+ifclosed = true;
+
+if (nargin >1)
+    if (isfield(opts,'ifclosed'))
+        ifclosed = opts.ifclosed;
+    end
+
+end
 
 pref = [];
 
-ifclosed=true;
 dim = size(r,1);
 k   = size(r,2);
 nch = size(r,3); 
@@ -21,8 +47,14 @@
 adjs(1,:) = (1:nch) - 1;
 adjs(2,:) = (1:nch) + 1;
 
-adjs(1,1)=nch;
-adjs(2,nch)=1;
+if ifclosed
+    adjs(1,1)=nch;
+    adjs(2,nch)=1;
+else
+    adjs(1,1)=-1;
+    adjs(2,nch)=-1;
+end
+
 
 
 %       up to here, everything has been done in parameter space, [ta,tb]
@@ -36,8 +68,16 @@
 for i = 1:nch
     
     rtmp  = squeeze(r(:,:,i));
-    dtmp  = rtmp*dermat;
-    d2tmp = dtmp*dermat;
+    if (~isempty(d))
+        dtmp = squeeze(d(:,:,i));
+    else
+        dtmp  = rtmp*dermat;
+    end
+    if (~isempty(d2))
+        d2tmp = squeeze(d2(:,:,i));
+    else
+        d2tmp = dtmp*dermat;    
+    end
     chnkr.rstor(:,:,i) = rtmp;
     chnkr.dstor(:,:,i) = dtmp;
     chnkr.d2stor(:,:,i) = d2tmp;

From 485902cca56cb77305363b1c2fc8becb160766ec Mon Sep 17 00:00:00 2001
From: jghoskins <32715900+jghoskins@users.noreply.github.com>
Date: Mon, 28 Apr 2025 15:25:14 -0400
Subject: [PATCH 13/15] Update chunkerpoints.m

---
 chunkie/@chunker/chunkerpoints.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/chunkie/@chunker/chunkerpoints.m b/chunkie/@chunker/chunkerpoints.m
index c3ff91ed..ccc145a0 100644
--- a/chunkie/@chunker/chunkerpoints.m
+++ b/chunkie/@chunker/chunkerpoints.m
@@ -1,4 +1,4 @@
-function chnkr = chunkerpoints(src)
+function chnkr = chunkerpoints(src,opts)
 
 d  = [];
 d2 = [];

From 7713ed40836ee53a9b643dd12b5170ab56eb21ce Mon Sep 17 00:00:00 2001
From: jghoskins <32715900+jghoskins@users.noreply.github.com>
Date: Mon, 28 Apr 2025 16:42:43 -0400
Subject: [PATCH 14/15] some documentation

For travis
---
 chunkie/+chnk/+smoother/smooth.m | 58 ++++++++++++++++++++++++++++++--
 chunkie/@chunker/chunkerpoints.m | 20 +++++++++++
 2 files changed, 76 insertions(+), 2 deletions(-)

diff --git a/chunkie/+chnk/+smoother/smooth.m b/chunkie/+chnk/+smoother/smooth.m
index ce8182f7..1f62e6ad 100644
--- a/chunkie/+chnk/+smoother/smooth.m
+++ b/chunkie/+chnk/+smoother/smooth.m
@@ -1,4 +1,58 @@
 function [chnkr,varargout] = smooth(verts,opts)
+%SMOOTH attempts to generate a smoothed curve from a polygon defined 
+%   by a collection of vertices, using the method in:
+%
+%   "A fast boundary integral method for high-order 
+%           multiscale mesh generation"
+%
+%   Felipe Vico, Leslie Greengard, Michael O'Neil, Manas Rachh
+%
+%   Though it should be noted that nothing about this code is 
+%   expected to be fast.
+%
+% Syntax: chnkr = chnk.smoother.smooth(verts,opts)
+%
+% Input:
+%   verts - a 2 by n matrix containing the x and y locations of the 
+%           vertices of the original polygon
+%   opts  - options structure. available options (default settings)
+%           opts.k  - number of Legendre nodes per panel on the 
+%                     output chunker object                     (16)
+%           opts.kquad - number of quadrature nodes used by
+%                     surface smoother in computing integrals. Must 
+%                     be greater than k.                        (32)
+%           opts.nchs  - the number of panels per edge in the final
+%                     chunker. Can either by a single positive integer, 
+%                     in which case each edge will have the same number 
+%                     of panels, or a vector of length (# of vertices), 
+%                     in which case the ith entry will specify the number 
+%                     of panels on the edge joining the ith vertex to the
+%                     (i+1)st vertex.                            (3)
+%           opts.n_newton - the maximum number of steps of newton to 
+%                     take when trying to smooth the curve.     (200)
+%           opts.etol   - the error tolerance to be used when running 
+%                     newton.                                   (1E-10)
+%           opts.lam   - the smoothing parameter. Larger lambdas will 
+%                     make the outputted curve less smooth. Small lambdas 
+%                     increase the risk of the smoother failing to 
+%                     converge. If the smoother fails, try running with
+%                     a smaller lam.                            (10)
+%           opts.step_fact - a parameter to multiply the step lengths 
+%                     by when doing newton steps. If the smoother fails 
+%                     to converge, then sometimes decreasing this can 
+%                     help. Should always be <=1                (1)
+
+% Output:
+%   chnkr - the requested chunker containing the smoothed geometry
+%   err (optional) - an approximate error of the curve smoother error
+%   err_by_pt (optional) - a vector containing point by point errors
+%
+% The method works by convolving the polygon with a normalized Gaussian 
+% (with a spatially dependent variance). The outputted curve ideally 
+% is chosen to lie on the 1/2 level set of the resulting function. Errors
+% here represent the difference between the values of this function at 
+% the final points on the outputted curve and 1/2.
+%
 
     nv = size(verts,2);
 
@@ -28,12 +82,12 @@
         n_newton = opts.n_newton;
     end
 
-    etol = 1E-6;
+    etol = 1E-10;
     if isfield(opts,'etol')
         etol = opts.etol;
     end
 
-    lam = 2.5;
+    lam = 10;
     if isfield(opts,'lam')
         lam = opts.lam;
     end
diff --git a/chunkie/@chunker/chunkerpoints.m b/chunkie/@chunker/chunkerpoints.m
index ccc145a0..ade8b0bf 100644
--- a/chunkie/@chunker/chunkerpoints.m
+++ b/chunkie/@chunker/chunkerpoints.m
@@ -1,4 +1,24 @@
 function chnkr = chunkerpoints(src,opts)
+%CHUNKERPOINTS create a chunker object given a set of Gauss-Legendre 
+% panels discretizing a curve.
+%
+% Syntax: chnkr = chunkerpoints(src,opts)
+%
+% Input: 
+%   src - If src is of class 'double', then it should be a 
+%         dim x k x nch array. dim is the dimension in which the curve 
+%         lives (typically 2), k is the number of Gauss-Legendre nodes 
+%         per panel (typically 16), and nch is the total number of panels.
+%   opts - a structure containing optional argumensts (defaults)
+%       opts.ifclosedtrue - a boolean. True if the desired curve 
+%           is closed, and False if not. Note that no checking is done. 
+%           All this changes is the corresponding chunker adjacency info.
+%
+% Output:
+%   chnkr - a chunker object containing the discretization of the domain
+%   
+% see also CHUNKERPOLY, CHUNKERPREF, CHUNKER
+%
 
 d  = [];
 d2 = [];

From 6e0cfe5bd48c1094286571b637ce3d5bde94129e Mon Sep 17 00:00:00 2001
From: Travis Askham <askhamwhat@gmail.com>
Date: Mon, 28 Apr 2025 16:50:40 -0400
Subject: [PATCH 15/15] small tweak

---
 chunkie/+chnk/+smoother/smooth.m | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/chunkie/+chnk/+smoother/smooth.m b/chunkie/+chnk/+smoother/smooth.m
index 1f62e6ad..a29b830f 100644
--- a/chunkie/+chnk/+smoother/smooth.m
+++ b/chunkie/+chnk/+smoother/smooth.m
@@ -41,7 +41,7 @@
 %                     by when doing newton steps. If the smoother fails 
 %                     to converge, then sometimes decreasing this can 
 %                     help. Should always be <=1                (1)
-
+%
 % Output:
 %   chnkr - the requested chunker containing the smoothed geometry
 %   err (optional) - an approximate error of the curve smoother error
@@ -139,4 +139,4 @@
     if (nargout > 2)
         varargout{2} = err_by_pt;
     end
-end
\ No newline at end of file
+end