Initial commit: chunkie examples for numbl
15 changed files+560−0
.github/workflows/deploy.ymladded+43−0View file
@@ -0,0 +1,43 @@
1+name: Deploy numbl project to GitHub Pages
2+
3+# Bundles this numbl project with the browser IDE and publishes it to GitHub
4+# Pages on every push to main.
5+#
6+# One-time setup: Settings → Pages → "Build and deployment" → Source →
7+# "GitHub Actions".
8+
9+on:
10+ push:
11+ branches: [main]
12+ workflow_dispatch:
13+
14+permissions:
15+ contents: read
16+ pages: write
17+ id-token: write
18+
19+concurrency:
20+ group: pages
21+ cancel-in-progress: false
22+
23+jobs:
24+ build:
25+ runs-on: ubuntu-latest
26+ steps:
27+ - uses: actions/checkout@v4
28+ - uses: flatironinstitute/numbl/.github/actions/build-site@main
29+ with:
30+ project-dir: .
31+ # Build numbl from the main branch (development version) rather than
32+ # the published npm release
33+ numbl-ref: main
34+
35+ deploy:
36+ needs: build
37+ runs-on: ubuntu-latest
38+ environment:
39+ name: github-pages
40+ url: ${{ steps.deployment.outputs.page_url }}
41+ steps:
42+ - id: deployment
43+ uses: actions/deploy-pages@v4
.gitignoreadded+4−0View file
@@ -0,0 +1,4 @@
1+# Local output from `numbl build-site`
2+/_site/
3+/dist/
4+node_modules/
README.mdadded+56−0View file
@@ -0,0 +1,56 @@
1+# chunkie examples
2+
3+An interactive companion to
4+[chunkie](https://github.com/fastalgorithms/chunkie) — a MATLAB package for
5+boundary integral equations in two dimensions. Everything here runs in your
6+browser through [numbl](https://numbl.org): pick a script on the left and click
7+**Run**.
8+
9+Every script begins with
10+
11+```matlab
12+mip load --install magland/magland/chunkie
13+```
14+
15+which fetches chunkie on first use.
16+
17+## Geometry
18+
19+chunkie represents a curve as a `chunker`: a panelization into high-order
20+Gauss-Legendre chunks.
21+
22+- [`chunkie_ex01_circle.m`](chunkie_ex01_circle.m) — chunk up a circle; plot
23+ points and normals.
24+- [`chunkie_ex02_bymode.m`](chunkie_ex02_bymode.m) — a star-shaped curve from a
25+ Fourier series for the radius.
26+- [`chunkie_ex03_starfish.m`](chunkie_ex03_starfish.m) — the classic starfish.
27+- [`chunkie_ex04_chunkerfit.m`](chunkie_ex04_chunkerfit.m) — fit a smooth curve
28+ to sample points (`chunkerfit`).
29+- [`chunkie_ex05_barbell.m`](chunkie_ex05_barbell.m) — a rounded polygon
30+ (`chunkerpoly`).
31+- [`chunkie_ex06_combined.m`](chunkie_ex06_combined.m) — transform and merge
32+ several curves into one multi-component domain.
33+- [`chunkie_ex07_interior.m`](chunkie_ex07_interior.m) — find which points lie
34+ inside a domain (`chunkerinterior`).
35+
36+## Solving integral equations
37+
38+Each PDE is recast as a boundary integral equation, discretized with
39+`chunkermat`, solved with GMRES, and evaluated on a grid with
40+`chunkerkerneval`.
41+
42+- [`chunkie_ex00_starfish.m`](chunkie_ex00_starfish.m) — Helmholtz scattering
43+ off a starfish (combined-field representation).
44+- [`chunkie_ex08_laplace.m`](chunkie_ex08_laplace.m) — a Laplace Neumann
45+ problem (single-layer with nullspace regularization).
46+- [`chunkie_ex09_helmholtz.m`](chunkie_ex09_helmholtz.m) — sound-soft Helmholtz
47+ scattering off a peanut.
48+- [`chunkie_ex10_stokes.m`](chunkie_ex10_stokes.m) — Stokes flow through a
49+ multiply-connected pipe.
50+
51+## How it's deployed
52+
53+On every push to `main`, the workflow in
54+[`.github/workflows/deploy.yml`](.github/workflows/deploy.yml) bundles these
55+files with the numbl browser IDE and publishes the result to GitHub Pages. Edit
56+`numbl-project.json` to change the title or which file opens first.
chunkie_ex00_starfish.madded+53−0View file
@@ -0,0 +1,53 @@
1+mip load --install magland/magland/chunkie;
2+tic;
3+
4+% planewave definitions
5+
6+kvec = 20*[1;-1.5];
7+zk = norm(kvec);
8+planewave = @(kvec,r) exp(1i*sum(bsxfun(@times,kvec(:),r(:,:)))).';
9+
10+% discretize domain
11+
12+narms = 5;
13+amp = 0.5;
14+chnkr = chunkerfunc(@(t) starfish(t,narms,amp),struct('maxchunklen',4/zk));
15+
16+% build CFIE and solve
17+
18+fkern = kernel('helm','c',zk,[1,-zk*1i]);
19+sysmat = chunkermat(chnkr,fkern);
20+sysmat = 0.5*eye(chnkr.k*chnkr.nch) + sysmat;
21+
22+rhs = -planewave(kvec(:),chnkr.r(:,:));
23+sol = gmres(sysmat,rhs,[],1e-13,100);
24+
25+% evaluate at targets
26+
27+x1 = linspace(-3,3,400);
28+[xxtarg,yytarg] = meshgrid(x1,x1);
29+targets = [xxtarg(:).';yytarg(:).'];
30+
31+in = chunkerinterior(chnkr,targets);
32+out = ~in;
33+
34+uscat = chunkerkerneval(chnkr,fkern,sol,targets(:,out));
35+
36+uin = planewave(kvec,targets(:,out));
37+utot = uscat(:)+planewave(kvec,targets(:,out));
38+
39+% plot
40+
41+maxu = max(abs(utot(:)));
42+figure()
43+zztarg = nan(size(xxtarg));
44+zztarg(out) = utot;
45+h=pcolor(xxtarg,yytarg,imag(zztarg));
46+set(h,'EdgeColor','none')
47+hold on
48+plot(chnkr,'LineWidth',2)
49+axis equal tight
50+colormap(redblue)
51+caxis([-maxu,maxu])
52+
53+toc;
chunkie_ex01_circle.madded+19−0View file
@@ -0,0 +1,19 @@
1+% Example 1: Chunk up a circle and plot geometry + normals.
2+
3+mip load --install magland/magland/chunkie;
4+
5+tic;
6+
7+rad = 2; ctr = [1.0;-0.5];
8+circfun = @(t) ctr + rad*[cos(t(:).');sin(t(:).')];
9+
10+chnkr1 = chunkerfunc(circfun);
11+
12+toc;
13+
14+figure(1)
15+clf
16+plot(chnkr1,'b-x')
17+hold on
18+quiver(chnkr1,'r')
19+axis equal tight
chunkie_ex02_bymode.madded+19−0View file
@@ -0,0 +1,19 @@
1+% Example 2: Star-shaped domain defined by a cosine/sine series for the radius.
2+
3+mip load --install magland/magland/chunkie;
4+
5+tic;
6+
7+rng(0)
8+modes = randn(11,1); modes(1) = 1.1*sum(abs(modes(2:end)));
9+ctr = [1.0;-0.5];
10+chnkr2 = chunkerfunc(@(t) chnk.curves.bymode(t,modes,ctr));
11+
12+toc;
13+
14+figure(2)
15+clf
16+plot(chnkr2,'b-x')
17+hold on
18+quiver(chnkr2,'r')
19+axis equal tight
chunkie_ex03_starfish.madded+18−0View file
@@ -0,0 +1,18 @@
1+% Example 3: Classic starfish domain.
2+
3+mip load --install magland/magland/chunkie;
4+
5+tic;
6+
7+narms = 5;
8+amp = 0.5;
9+chnkr3 = chunkerfunc(@(t) starfish(t,narms,amp));
10+
11+toc;
12+
13+figure(3)
14+clf
15+plot(chnkr3,'b-x')
16+hold on
17+quiver(chnkr3,'r')
18+axis equal tight
chunkie_ex04_chunkerfit.madded+43−0View file
@@ -0,0 +1,43 @@
1+% Example 4: Fit a spline curve to sample points using chunkerfit, with
2+% varying tolerance and an option to split between sample points.
3+
4+mip load --install magland/magland/chunkie;
5+
6+tic;
7+
8+rng(0)
9+n = 20;
10+tt = sort(2*pi*rand(n,1));
11+r = chnk.curves.bymode(tt, [2 0.5 0.2 0.7]);
12+
13+opts = [];
14+opts.ifclosed = true;
15+opts.cparams = [];
16+opts.cparams.eps = 1e-3;
17+opts.pref = [];
18+opts.pref.k = 16;
19+chnkr = chunkerfit(r, opts);
20+
21+opts.splitatpoints = true;
22+chnkrsp = chunkerfit(r, opts);
23+
24+opts.cparams.eps = 1e-9;
25+opts.splitatpoints = false;
26+chnkr2 = chunkerfit(r, opts);
27+
28+opts.splitatpoints = true;
29+chnkrsp2 = chunkerfit(r, opts);
30+
31+toc;
32+
33+figure(7)
34+clf
35+tiledlayout(2,2,"TileSpacing","compact");
36+nexttile; plot(chnkr,'k-x'); title(sprintf("nch = %d, eps=1e-3",chnkr.nch));
37+hold on; plot(r(1,:),r(2,:),'bd');
38+nexttile; plot(chnkrsp,'k-x'); title(sprintf("nch = %d, eps=1e-3\n (splitting)",chnkrsp.nch));
39+hold on; plot(r(1,:),r(2,:),'bd');
40+nexttile; plot(chnkr2,'k-x'); title(sprintf("nch = %d, eps=1e-9",chnkr2.nch));
41+hold on; plot(r(1,:),r(2,:),'bd');
42+nexttile; plot(chnkrsp2,'k-x'); title(sprintf("nch = %d, eps=1e-9\n (splitting)",chnkrsp2.nch));
43+hold on; plot(r(1,:),r(2,:),'bd');
chunkie_ex05_barbell.madded+17−0View file
@@ -0,0 +1,17 @@
1+% Example 5: Rounded polygon (barbell) via chunkerpoly.
2+
3+mip load --install magland/magland/chunkie;
4+
5+tic;
6+
7+verts = chnk.demo.barbell(2.0,2.0,1.0,1.0);
8+chnkr4 = chunkerpoly(verts);
9+
10+toc;
11+
12+figure(4)
13+clf
14+plot(chnkr4,'b-x')
15+hold on
16+quiver(chnkr4,'r')
17+axis equal tight
chunkie_ex06_combined.madded+53−0View file
@@ -0,0 +1,53 @@
1+% Example 6: Combine a transformed copy of a star-shaped domain and a
2+% transformed copy of a circle into a single multi-component chunker, with
3+% consistent outward-pointing normals.
4+
5+mip load --install magland/magland/chunkie;
6+
7+tic;
8+
9+% recreate the circle (chnkr1)
10+rad = 2; ctr = [1.0;-0.5];
11+circfun = @(t) ctr + rad*[cos(t(:).');sin(t(:).')];
12+chnkr1 = chunkerfunc(circfun);
13+
14+% recreate the random-mode star (chnkr2)
15+rng(0)
16+modes = randn(11,1); modes(1) = 1.1*sum(abs(modes(2:end)));
17+ctr2 = [1.0;-0.5];
18+chnkr2 = chunkerfunc(@(t) chnk.curves.bymode(t,modes,ctr2));
19+
20+% rotate, reflect, and reverse the star
21+chnkr5 = chnkr2;
22+theta = pi/4;
23+chnkr5 = chnkr5.rotate(theta);
24+chnkr5 = chnkr5.reflect(pi/2);
25+chnkr5 = chnkr5.reverse();
26+
27+% affine-transform the circle and reverse its orientation
28+chnkr6 = chnkr1;
29+A = 0.5*[2 -1; 1 1];
30+r1 = [-1;0.5];
31+chnkr6 = r1 + A*chnkr6;
32+chnkr6 = chnkr6.reverse();
33+
34+% recenter the circle inside the star and scale it down so it fits
35+mins5 = min(chnkr5); maxs5 = max(chnkr5);
36+mins6 = min(chnkr6); maxs6 = max(chnkr6);
37+ctr5 = 0.5*(mins5+maxs5);
38+ctr6 = 0.5*(mins6+maxs6);
39+ext5 = min(maxs5 - mins5);
40+ext6 = max(maxs6 - mins6);
41+scale = 0.3*ext5/ext6;
42+chnkr6 = chnkr6 + (-ctr6);
43+chnkr6 = scale*chnkr6;
44+chnkr6 = chnkr6 + ctr5;
45+
46+% merge into one multi-component domain
47+chnkr7 = merge([chnkr5,chnkr6]);
48+
49+toc;
50+
51+figure(5); clf
52+plot(chnkr7,'b-x'); hold on; quiver(chnkr7,'r')
53+axis equal tight
chunkie_ex07_interior.madded+61−0View file
@@ -0,0 +1,61 @@
1+% Example 7: Find the points on the interior of a (multi-component)
2+% chunker domain using chunkerinterior, and visualize them.
3+
4+mip load --install magland/magland/chunkie;
5+
6+tic;
7+
8+% recreate the circle (chnkr1)
9+rad = 2; ctr = [1.0;-0.5];
10+circfun = @(t) ctr + rad*[cos(t(:).');sin(t(:).')];
11+chnkr1 = chunkerfunc(circfun);
12+
13+% recreate the random-mode star (chnkr2)
14+rng(0)
15+modes = randn(11,1); modes(1) = 1.1*sum(abs(modes(2:end)));
16+ctr2 = [1.0;-0.5];
17+chnkr2 = chunkerfunc(@(t) chnk.curves.bymode(t,modes,ctr2));
18+
19+% rebuild the merged domain (chnkr7)
20+chnkr5 = chnkr2;
21+chnkr5 = chnkr5.rotate(pi/4);
22+chnkr5 = chnkr5.reflect(pi/2);
23+chnkr5 = chnkr5.reverse();
24+
25+chnkr6 = chnkr1;
26+A = 0.5*[2 -1; 1 1];
27+r1 = [-1;0.5];
28+chnkr6 = r1 + A*chnkr6;
29+chnkr6 = chnkr6.reverse();
30+
31+% recenter the circle inside the star and scale it down so it fits
32+mins5 = min(chnkr5); maxs5 = max(chnkr5);
33+mins6 = min(chnkr6); maxs6 = max(chnkr6);
34+ctr5 = 0.5*(mins5+maxs5);
35+ctr6 = 0.5*(mins6+maxs6);
36+ext5 = min(maxs5 - mins5);
37+ext6 = max(maxs6 - mins6);
38+scale = 0.3*ext5/ext6;
39+chnkr6 = chnkr6 + (-ctr6);
40+chnkr6 = scale*chnkr6;
41+chnkr6 = chnkr6 + ctr5;
42+
43+chnkr7 = merge([chnkr5,chnkr6]);
44+
45+% grid of test points
46+mins = min(chnkr7); maxs = max(chnkr7);
47+x1 = linspace(mins(1),maxs(1)); y1 = linspace(mins(2),maxs(2));
48+[xx,yy] = meshgrid(x1,y1);
49+pts = [xx(:).'; yy(:).'];
50+
51+% interior test
52+in = chunkerinterior(chnkr7,pts);
53+
54+zz = nan(size(xx));
55+zz(in) = 1;
56+
57+toc;
58+
59+figure(6)
60+h = pcolor(xx,yy,zz); set(h,'EdgeColor','none');
61+axis equal tight
chunkie_ex08_laplace.madded+41−0View file
@@ -0,0 +1,41 @@
1+% Example 8: Laplace Neumann BVP on a starfish domain via single-layer
2+% representation, with the W (onesmat) regularization added so that the
3+% second-kind integral equation is invertible.
4+
5+mip load --install magland/magland/chunkie;
6+
7+tic;
8+
9+% chunker discretization of a starfish domain
10+chnkr = chunkerfunc(@(t) starfish(t));
11+
12+% boundary condition. by symmetry of the starfish this integrates to zero
13+pwfun = @(r) r(1,:).^2.*r(2,:);
14+rhs = pwfun(chnkr.r); rhs = rhs(:);
15+
16+% kernel: normal derivative of the Laplace single-layer
17+kernsp = kernel('lap','sprime');
18+
19+% matrix discretization of the BIE
20+sysmat = chunkermat(chnkr,kernsp);
21+sysmat = sysmat + 0.5*eye(chnkr.npt) + onesmat(chnkr);
22+
23+% solve
24+sigma = gmres(sysmat,rhs);
25+
26+% grid for plotting
27+x1 = linspace(-2,2,300);
28+[xx,yy] = meshgrid(x1,x1);
29+targs = [xx(:).'; yy(:).'];
30+in = chunkerinterior(chnkr,targs);
31+uu = nan(size(xx));
32+
33+% evaluate solution using the single layer (not its derivative)
34+kerns = kernel('lap','s');
35+uu(in) = chunkerkerneval(chnkr,kerns,sigma,targs(:,in));
36+
37+toc;
38+
39+figure(1); clf
40+h = pcolor(xx,yy,uu); set(h,'EdgeColor','none'); colorbar
41+hold on; plot(chnkr,'k'); colormap(redblue)
chunkie_ex09_helmholtz.madded+47−0View file
@@ -0,0 +1,47 @@
1+% Example 9: Helmholtz sound-soft scattering by a peanut shape, using the
2+% combined-field (CFIE) representation D - i*k*S.
3+
4+mip load --install magland/magland/chunkie;
5+
6+tic;
7+
8+% chunker discretization of a peanut-shaped domain
9+modes = [1.25,-0.25,0,0.5];
10+ctr = [0;0];
11+chnkr = chunkerfunc(@(t) chnk.curves.bymode(t,modes,ctr));
12+
13+% incident plane wave defines the boundary condition
14+kwav = 3*[-1,-5];
15+pwfun = @(r) exp(1i*kwav*r(:,:));
16+rhs = -pwfun(chnkr.r); rhs = rhs(:);
17+
18+% CFIE kernel D - i*k*S
19+zk = norm(kwav);
20+coefs = [1,-1i*zk];
21+kerncfie = kernel('h','c',zk,coefs);
22+
23+% matrix discretization of the BIE (with 0.5*I from exterior jump)
24+sysmat = chunkermat(chnkr,kerncfie);
25+sysmat = sysmat + 0.5*eye(chnkr.npt);
26+
27+% solve
28+sigma = gmres(sysmat,rhs,[],1e-10,100);
29+
30+% grid for plotting (exterior)
31+x1 = linspace(-5,5,300);
32+[xx,yy] = meshgrid(x1,x1);
33+targs = [xx(:).'; yy(:).'];
34+in = chunkerinterior(chnkr,{x1,x1});
35+uu = nan(size(xx));
36+
37+% evaluate scattered field outside, then add incident field
38+uu(~in) = chunkerkerneval(chnkr,kerncfie,sigma,targs(:,~in));
39+uu(~in) = uu(~in) + pwfun(targs(:,~in)).';
40+
41+toc;
42+
43+figure(1); clf
44+h = pcolor(xx,yy,real(uu)); set(h,'EdgeColor','none'); colorbar
45+colormap(redblue); umax = max(abs(uu(:))); clim([-umax,umax]);
46+hold on
47+plot(chnkr,'k')
chunkie_ex10_stokes.madded+82−0View file
@@ -0,0 +1,82 @@
1+% Example 10: Stokes flow in a multiply-connected pipe-like domain, using
2+% the combined-layer Stokes representation u = (D - S)[sigma] with the
3+% nullspace correction W = normonesmat.
4+
5+mip load --install magland/magland/chunkie;
6+
7+tic;
8+
9+% outer boundary: peanut shape
10+modes = [2.5,0,0,1];
11+ctr = [0;0];
12+chnkrouter = chunkerfunc(@(t) chnk.curves.bymode(t,modes,ctr));
13+
14+% inner boundaries: small circles, reversed so orientations are consistent
15+chnkrcirc = chunkerfunc(@(t) chnk.curves.bymode(t,0.25,[0;0]));
16+chnkrcirc = reverse(chnkrcirc);
17+centers = [ [-2:2, -2:2]; [(0.7 + 0.2*(-1).^(-2:2)) , ...
18+ (-0.7 + 0.2*(-1).^(-2:2))]];
19+centers = centers + 0.1*randn(size(centers));
20+
21+% assemble outer + shifted circles into one chunker
22+chnkrlist = [chnkrouter];
23+for j = 1:size(centers,2)
24+ chnkr1 = chnkrcirc;
25+ chnkr1 = chnkr1.move([0;0],centers(:,j));
26+ chnkrlist = [chnkrlist chnkr1];
27+end
28+chnkr = merge(chnkrlist);
29+
30+% boundary velocity (Gaussian profile on the outer boundary, zero on holes)
31+wid = 0.3;
32+f = @(r) [exp(-r(2,:).^2/(2*wid^2)); zeros(size(r(2,:)))];
33+rhsout = f(chnkrouter.r(:,:)); rhsout = rhsout(:);
34+rhs = [rhsout; zeros(2*10*chnkrcirc.npt,1)];
35+
36+% combined-layer Stokes kernel D - S
37+c = -1;
38+mu = 1;
39+kerncvel = kernel('stok','dvel',mu) + c*kernel('stok','svel',mu);
40+
41+% matrix discretization of the BIE
42+cmat = chunkermat(chnkr,kerncvel);
43+
44+% identity term and nullspace correction
45+W = normonesmat(chnkr);
46+sysmat = cmat - 0.5*eye(2*chnkr.npt) + W;
47+
48+% solve
49+sigma = gmres(sysmat,rhs,[],1e-10,100);
50+
51+% grid for plotting (interior)
52+x1 = linspace(-3.75,3.75,200);
53+y1 = linspace(-2,2,100);
54+[xx,yy] = meshgrid(x1,y1);
55+targs = [xx(:).'; yy(:).'];
56+in = chunkerinterior(chnkr,{x1,y1});
57+uu = nan([2,size(xx)]);
58+pres = nan(size(xx));
59+
60+% velocity from layer potential
61+uu(:,in) = reshape(chunkerkerneval(chnkr,kerncvel,sigma,targs(:,in)),2,nnz(in));
62+
63+% pressure from the corresponding pressure kernel
64+kerncpres = kernel('stok','dpres',mu) + c*kernel('stok','spres',mu);
65+opts = []; opts.eps = 1e-3;
66+pres(in) = chunkerkerneval(chnkr,kerncpres,sigma,targs(:,in),opts);
67+
68+toc;
69+
70+figure(1); clf
71+h = pcolor(xx,yy,pres); set(h,'EdgeColor','none'); colorbar
72+colormap("parula");
73+hold on
74+plot(chnkr,'k','LineWidth',2); axis equal
75+
76+u = reshape(uu(1,:,:),size(xx)); v = reshape(uu(2,:,:),size(xx));
77+startt = linspace(pi-pi/12,pi+pi/12,20);
78+startr = 0.99*chnk.curves.bymode(startt,modes,[0;0]);
79+startx = startr(1,:); starty = startr(2,:);
80+
81+sl = streamline(xx,yy,u,v,startx,starty);
82+set(sl,'LineWidth',2,'Color','w')
numbl-project.jsonadded+4−0View file
@@ -0,0 +1,4 @@
1+{
2+ "title": "chunkie examples",
3+ "entry": "README.md"
4+}