Intro to surfacefun: interactive numbl example project
26 changed files+657−0
.github/workflows/deploy.ymladded+40−0View file
@@ -0,0 +1,40 @@
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+
32+ deploy:
33+ needs: build
34+ runs-on: ubuntu-latest
35+ environment:
36+ name: github-pages
37+ url: ${{ steps.deployment.outputs.page_url }}
38+ steps:
39+ - id: deployment
40+ 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+88−0View file
@@ -0,0 +1,88 @@
1+# Intro to surfacefun
2+
3+An interactive introduction to
4+[surfacefun](https://github.com/danfortunato/surfacefun) — a package for
5+numerically computing with functions on surfaces, and solving PDEs on them, with
6+high-order accuracy. Everything here runs in your browser through
7+[numbl](https://numbl.org): pick a script on the left and click **Run**. Each
8+script is self-contained and draws one or more figures (use the figure tabs).
9+
10+Every script begins with
11+
12+```matlab
13+mip load --install flatironinstitute/flatironinstitute/surfacefun
14+```
15+
16+which fetches surfacefun (and its dependency, Chebfun) on first use.
17+
18+## 1. Surface meshes
19+
20+A `surfacemesh` is a collection of high-order quadrilateral patches whose union
21+defines a surface in 3-D.
22+
23+- [`mesh_single_patch.m`](mesh_single_patch.m) — build a surface from one
24+ high-order Chebyshev patch.
25+- [`mesh_sphere.m`](mesh_sphere.m) — the built-in "cubed sphere", plus mesh
26+ queries (patches, order, area, volume, bounding box).
27+- [`mesh_blob.m`](mesh_blob.m) — a smooth random "blob".
28+- [`mesh_torus.m`](mesh_torus.m) — a genus-1 torus.
29+- [`mesh_twisted_torus.m`](mesh_twisted_torus.m) — patches need not join
30+ smoothly.
31+- [`mesh_mobius.m`](mesh_mobius.m) — surfaces need not be orientable.
32+- [`mesh_refine.m`](mesh_refine.m) — change the order (`resample`) or split
33+ patches (`refine`).
34+
35+## 2. Scalar functions
36+
37+A `surfacefun` is a scalar field on a `surfacemesh`, built from a handle of the
38+Cartesian coordinates `(x,y,z)`.
39+
40+- [`fun_plot.m`](fun_plot.m) — define and plot a function on the sphere.
41+- [`fun_arithmetic.m`](fun_arithmetic.m) — overloaded `+`, `*`, `abs`, ...
42+- [`fun_contour.m`](fun_contour.m) — contour (level-set) plot on the surface.
43+- [`fun_gradient.m`](fun_gradient.m) — the surface gradient `grad(f)`.
44+- [`fun_laplacian.m`](fun_laplacian.m) — the surface Laplacian `lap(f)`.
45+- [`fun_integral_norm.m`](fun_integral_norm.m) — `integral`, `mean`, and norms.
46+
47+## 3. Vector fields
48+
49+A `surfacefunv` is a three-component vector field over a `surfacemesh`.
50+
51+- [`vec_normals.m`](vec_normals.m) — the unit normal field, as a quiver plot.
52+- [`vec_gradient.m`](vec_gradient.m) — `grad(f)` is tangent to the surface.
53+- [`vec_mean_curvature.m`](vec_mean_curvature.m) — mean curvature from
54+ `div(n)/2`.
55+- [`vec_curl.m`](vec_curl.m) — the surface curl of a vector field.
56+
57+## 4. Solving PDEs
58+
59+`surfaceop` is a fast direct solver for second-order linear elliptic PDEs on
60+surfaces.
61+
62+- [`pde_laplace_beltrami.m`](pde_laplace_beltrami.m) — a Laplace-Beltrami
63+ problem on the sphere, checked against an exact spherical-harmonic solution.
64+- [`pde_helmholtz.m`](pde_helmholtz.m) — a variable-coefficient Helmholtz
65+ problem on a stellarator.
66+- [`pde_open_surface.m`](pde_open_surface.m) — an open surface with Dirichlet
67+ boundary data, reusing one factorization for several right-hand sides.
68+
69+## 5. Worked examples
70+
71+- [`ex_hodge.m`](ex_hodge.m) — the Hodge decomposition of a tangential vector
72+ field into curl-free, divergence-free, and harmonic parts.
73+- [`ex_reaction_diffusion.m`](ex_reaction_diffusion.m) — the complex
74+ Ginzburg-Landau equation, integrated in time with a reused implicit solver.
75+
76+## Notes
77+
78+The scripts keep the polynomial order and refinement modest so they run quickly
79+in the browser. Raise `p` and `nref` for higher accuracy and smoother pictures.
80+The numerics follow the surfacefun documentation at
81+<https://surfacefun.readthedocs.io>.
82+
83+## How it's deployed
84+
85+On every push to `main`, the workflow in
86+[`.github/workflows/deploy.yml`](.github/workflows/deploy.yml) bundles these
87+files with the numbl browser IDE and publishes the result to GitHub Pages.
88+Edit `numbl-project.json` to change the title or which file opens first.
ex_hodge.madded+50−0View file
@@ -0,0 +1,50 @@
1+% ex_hodge.m — Hodge decomposition of a tangential vector field.
2+%
3+% Any tangential vector field f on a surface splits into curl-free,
4+% divergence-free, and harmonic parts:
5+%
6+% f = grad(u) + n x grad(v) + w.
7+%
8+% hodge() computes this by solving two Laplace-Beltrami problems. We use a
9+% modest torus here so the demo stays quick (the docs use a finer mesh).
10+
11+mip load --install flatironinstitute/flatironinstitute/surfacefun
12+
13+p = 8;
14+nu = 8;
15+nv = 16;
16+dom = surfacemesh.torus(p + 1, nu, nv);
17+
18+% A smooth random tangential vector field.
19+rng(0)
20+bb = boundingbox(dom);
21+gx = randnfun3(10, bb);
22+gy = randnfun3(10, bb);
23+gz = randnfun3(10, bb);
24+g = cross([0 1 1], surfacefunv(@(x, y, z) gx(x, y, z), ...
25+ @(x, y, z) gy(x, y, z), ...
26+ @(x, y, z) gz(x, y, z), dom));
27+vn = normal(dom);
28+f = -cross(vn, vn, g); % project g onto the tangent plane
29+
30+% Compute the decomposition.
31+[u, v, w] = hodge(f);
32+
33+figure(1);
34+quiver(f, 0.3, 4)
35+title('Tangential field f');
36+
37+figure(2);
38+quiver(grad(u), 0.3, 4)
39+title('Curl-free part grad(u)');
40+
41+figure(3);
42+quiver(cross(vn, grad(v)), 0.3, 4)
43+title('Divergence-free part n x grad(v)');
44+
45+figure(4);
46+quiver(w, 0.3, 4)
47+title('Harmonic part w');
48+
49+% The harmonic remainder should be nearly divergence-free.
50+fprintf('|| div(w) || = %.3e\n', norm(div(w)));
ex_reaction_diffusion.madded+51−0View file
@@ -0,0 +1,51 @@
1+% ex_reaction_diffusion.m — A time-dependent PDE on a surface.
2+%
3+% The complex Ginzburg-Landau equation
4+%
5+% du/dt = delta (1 + i*alpha) Delta_Gamma u + u - (1 + i*beta) u |u|^2
6+%
7+% is integrated with a semi-implicit (backward-Euler) scheme: the
8+% diffusion term is treated implicitly, the reaction term explicitly. The
9+% diffusion operator is time-independent, so one surfaceop factorization
10+% is reused at every step. We use a sphere and a short run to stay quick.
11+
12+mip load --install flatironinstitute/flatironinstitute/surfacefun
13+
14+% Parameters.
15+dt = 0.1;
16+nsteps = 60;
17+alpha = 0;
18+beta = 1.5;
19+delta = 5e-3;
20+
21+% Nonlinear (reaction) term.
22+N = @(u) u - (1 + beta*1i)*u .* (abs(u).^2);
23+
24+% Mesh and the reusable implicit solver for (I - dt*delta*(1+i*alpha) Lap).
25+p = 8;
26+nref = 1;
27+dom = surfacemesh.sphere(p + 1, nref);
28+pdo = struct('lap', -dt*delta*(1 + alpha*1i), 'c', 1);
29+L = surfaceop(dom, pdo);
30+L.build();
31+
32+% Random initial condition.
33+rng(1)
34+f = randnfun3(0.5, boundingbox(dom));
35+u = surfacefun(@(x, y, z) f(x, y, z), dom);
36+
37+figure(1);
38+plot(real(u)), colorbar
39+title('Initial condition');
40+
41+% Time stepping: each step is one solve with a fresh right-hand side.
42+for k = 1:nsteps
43+ L.rhs = u + dt*N(u);
44+ u = L.solve();
45+end
46+
47+figure(2);
48+plot(real(u)), colorbar
49+title(sprintf('Re(u) after %d steps (t = %.1f)', nsteps, nsteps*dt));
50+
51+fprintf('done: %d steps, ||Re(u)|| = %.4f\n', nsteps, norm(real(u)));
fun_arithmetic.madded+19−0View file
@@ -0,0 +1,19 @@
1+% fun_arithmetic.m — Overloaded arithmetic on surfacefuns.
2+%
3+% Standard MATLAB arithmetic is overloaded for surfacefuns, so you can
4+% combine fields with +, -, .*, abs, etc.
5+
6+mip load --install flatironinstitute/flatironinstitute/surfacefun
7+
8+p = 16;
9+nref = 2;
10+dom = surfacemesh.sphere(p + 1, nref);
11+
12+f = surfacefun(@(x, y, z) cos(6*x) .* y + exp(z), dom);
13+x = surfacefun(@(x, y, z) x, dom);
14+
15+g = abs(f + 2*x);
16+
17+figure(1);
18+plot(g), colorbar
19+title('g = |f + 2x|');
fun_contour.madded+18−0View file
@@ -0,0 +1,18 @@
1+% fun_contour.m — Contour plot of a surfacefun.
2+%
3+% contour draws level sets of a scalar field on the surface. Contouring
4+% is more expensive than a plain plot, so this script uses a coarse mesh
5+% to stay quick; raise p / nref for smoother curves.
6+
7+mip load --install flatironinstitute/flatironinstitute/surfacefun
8+
9+p = 8;
10+nref = 0;
11+dom = surfacemesh.sphere(p + 1, nref);
12+
13+f = surfacefun(@(x, y, z) cos(6*x) .* y + exp(z), dom);
14+
15+figure(1);
16+contour(f, linewidth=2)
17+axis off
18+title('Contours of f');
fun_gradient.madded+20−0View file
@@ -0,0 +1,20 @@
1+% fun_gradient.m — Surface gradient of a scalar function.
2+%
3+% grad differentiates a surfacefun, automatically accounting for the
4+% on-surface metric. It returns the three Cartesian components of the
5+% surface gradient.
6+
7+mip load --install flatironinstitute/flatironinstitute/surfacefun
8+
9+p = 16;
10+nref = 2;
11+dom = surfacemesh.sphere(p + 1, nref);
12+
13+f = surfacefun(@(x, y, z) cos(6*x) .* y + exp(z), dom);
14+
15+[fx, fy, fz] = grad(f);
16+
17+figure(1);
18+subplot(131), plot(fx), title('\partial_x f')
19+subplot(132), plot(fy), title('\partial_y f')
20+subplot(133), plot(fz), title('\partial_z f')
fun_integral_norm.madded+17−0View file
@@ -0,0 +1,17 @@
1+% fun_integral_norm.m — Integrals, means and norms of a surfacefun.
2+%
3+% These reductions integrate over the surface using the on-surface
4+% measure. No plot here — just printed quantities.
5+
6+mip load --install flatironinstitute/flatironinstitute/surfacefun
7+
8+p = 16;
9+nref = 2;
10+dom = surfacemesh.sphere(p + 1, nref);
11+
12+f = surfacefun(@(x, y, z) cos(6*x) .* y + exp(z), dom);
13+
14+fprintf('integral(f) = %.6f\n', integral(f)); % definite integral over the surface
15+fprintf('mean(f) = %.6f\n', mean(f)); % integral / surface area
16+fprintf('L2 norm = %.6f\n', norm(f)); % sqrt(integral of f^2)
17+fprintf('Linf norm = %.6f\n', norm(f, inf)); % max |f|
fun_laplacian.madded+16−0View file
@@ -0,0 +1,16 @@
1+% fun_laplacian.m — The surface Laplacian (Laplace-Beltrami operator).
2+%
3+% Higher-order operators are built by composing the basic ones. lap(f)
4+% is the surface Laplacian of f.
5+
6+mip load --install flatironinstitute/flatironinstitute/surfacefun
7+
8+p = 16;
9+nref = 2;
10+dom = surfacemesh.sphere(p + 1, nref);
11+
12+f = surfacefun(@(x, y, z) cos(6*x) .* y + exp(z), dom);
13+
14+figure(1);
15+plot(lap(f)), colorbar
16+title('\Delta_\Gamma f');
fun_plot.madded+16−0View file
@@ -0,0 +1,16 @@
1+% fun_plot.m — A scalar function on a surface (surfacefun).
2+%
3+% A surfacefun represents a scalar field on a surfacemesh, built from a
4+% function handle of the Cartesian coordinates (x,y,z).
5+
6+mip load --install flatironinstitute/flatironinstitute/surfacefun
7+
8+p = 16;
9+nref = 2;
10+dom = surfacemesh.sphere(p + 1, nref);
11+
12+f = surfacefun(@(x, y, z) cos(6*x) .* y + exp(z), dom);
13+
14+figure(1);
15+plot(f), hold on, plot(dom), colorbar
16+title('f = cos(6x) y + e^z');
mesh_blob.madded+18−0View file
@@ -0,0 +1,18 @@
1+% mesh_blob.m — A random "blob" surface.
2+%
3+% surfacemesh.blob radially perturbs a spherical mesh by a smooth random
4+% field, so each run with a fresh seed produces a different shape.
5+
6+mip load --install flatironinstitute/flatironinstitute/surfacefun
7+
8+p = 16;
9+nref = 2;
10+
11+rng(0)
12+dom = surfacemesh.blob(p + 1, nref);
13+
14+figure(1);
15+plot(dom), camlight
16+title('Random blob');
17+
18+fprintf('patches: %d\n', length(dom));
mesh_mobius.madded+17−0View file
@@ -0,0 +1,17 @@
1+% mesh_mobius.m — A Mobius band.
2+%
3+% Surfaces need not even be orientable. surfacemesh.mobius builds the
4+% classic one-sided band.
5+
6+mip load --install flatironinstitute/flatironinstitute/surfacefun
7+
8+p = 16;
9+nu = 30;
10+nv = 7;
11+dom = surfacemesh.mobius(p + 1, nu, nv);
12+
13+figure(1);
14+plot(dom), camlight
15+title('Mobius band');
16+
17+fprintf('patches: %d\n', length(dom));
mesh_refine.madded+28−0View file
@@ -0,0 +1,28 @@
1+% mesh_refine.m — Changing a mesh's resolution.
2+%
3+% A mesh can be modified by changing the polynomial degree per patch
4+% ("p-refinement", via resample) or by splitting patches ("h-refinement",
5+% via refine).
6+
7+mip load --install flatironinstitute/flatironinstitute/surfacefun
8+
9+p = 16;
10+dom = surfacemesh.sphere(p + 1, 1);
11+
12+figure(1);
13+plot(dom, surface='off')
14+title(sprintf('Original: %d patches, order %d', length(dom), order(dom)));
15+
16+% p-refinement: resample each patch down to 4 nodes per dimension.
17+dom2 = resample(dom, 4);
18+
19+figure(2);
20+plot(dom2, surface='off')
21+title(sprintf('Resampled to order %d', order(dom2)));
22+
23+% h-refinement: split each patch into four.
24+dom3 = refine(dom2);
25+
26+figure(3);
27+plot(dom3, surface='off')
28+title(sprintf('Refined: %d patches', length(dom3)));
mesh_single_patch.madded+22−0View file
@@ -0,0 +1,22 @@
1+% mesh_single_patch.m — A surface from a single high-order patch.
2+%
3+% Every patch of a surfacemesh is represented on tensor-product Chebyshev
4+% nodes in [-1,1]^2. Here we build one patch and map it to the graph of
5+% z = cos(u)cos(v). The constructor takes cell arrays of node coordinates
6+% (one cell per patch).
7+
8+mip load --install flatironinstitute/flatironinstitute/surfacefun
9+
10+n = 16;
11+[u, v] = chebpts2(n);
12+
13+x = u;
14+y = v;
15+z = cos(u) .* cos(v);
16+dom = surfacemesh({x}, {y}, {z});
17+
18+figure(1);
19+plot(dom), camlight
20+title('Single high-order patch');
21+
22+fprintf('number of patches: %d\n', length(dom));
mesh_sphere.madded+23−0View file
@@ -0,0 +1,23 @@
1+% mesh_sphere.m — The built-in "cubed sphere".
2+%
3+% surfacemesh.sphere inflates a cube mesh onto the sphere. The first
4+% argument is the number of nodes per dimension (order p plus one); the
5+% second is the number of refinement levels.
6+
7+mip load --install flatironinstitute/flatironinstitute/surfacefun
8+
9+p = 16;
10+nref = 2;
11+dom = surfacemesh.sphere(p + 1, nref);
12+
13+figure(1);
14+plot(dom), camlight
15+title('Cubed sphere');
16+
17+% A few queries about the mesh:
18+fprintf('patches: %d\n', length(dom));
19+fprintf('order per patch: %d\n', order(dom));
20+fprintf('degrees of freedom:%d\n', numel(dom));
21+fprintf('surface area err: %.2e (vs 4*pi)\n', surfacearea(dom) - 4*pi);
22+fprintf('enclosed vol err: %.2e (vs 4/3*pi)\n', volume(dom) - 4/3*pi);
23+fprintf('bounding box: [%g %g %g %g %g %g]\n', boundingbox(dom));
mesh_torus.madded+18−0View file
@@ -0,0 +1,18 @@
1+% mesh_torus.m — A torus (genus-1 surface).
2+%
3+% surfacemesh.torus(p+1, nu, nv) builds a smooth torus with nu patches
4+% around the tube and nv patches around the ring. Surfaces of any genus
5+% are supported.
6+
7+mip load --install flatironinstitute/flatironinstitute/surfacefun
8+
9+p = 16;
10+nu = 8;
11+nv = 24;
12+dom = surfacemesh.torus(p + 1, nu, nv);
13+
14+figure(1);
15+plot(dom), camlight
16+title('Torus');
17+
18+fprintf('patches: %d\n', length(dom));
mesh_twisted_torus.madded+17−0View file
@@ -0,0 +1,17 @@
1+% mesh_twisted_torus.m — A twisted torus.
2+%
3+% The mesh need not be smooth between patches. surfacemesh.twisted_torus
4+% shears the patches as it wraps them around the ring.
5+
6+mip load --install flatironinstitute/flatironinstitute/surfacefun
7+
8+p = 16;
9+nu = 4;
10+nv = 32;
11+dom = surfacemesh.twisted_torus(p + 1, nu, nv);
12+
13+figure(1);
14+plot(dom), camlight
15+title('Twisted torus');
16+
17+fprintf('patches: %d\n', length(dom));
numbl-project.jsonadded+4−0View file
@@ -0,0 +1,4 @@
1+{
2+ "title": "Intro to surfacefun",
3+ "entry": "README.md"
4+}
pde_helmholtz.madded+25−0View file
@@ -0,0 +1,25 @@
1+% pde_helmholtz.m — A variable-coefficient surface Helmholtz problem.
2+%
3+% surfaceop handles general elliptic operators on surfaces of any genus.
4+% Here we solve (Delta_Gamma + c(x)) u = f on a genus-1 stellarator, with
5+% a spatially varying coefficient c.
6+
7+mip load --install flatironinstitute/flatironinstitute/surfacefun
8+
9+p = 16;
10+nu = 8;
11+nv = 24;
12+dom = surfacemesh.stellarator(p + 1, nu, nv);
13+
14+% Operator: Laplacian plus a zeroth-order term c(x,y,z).
15+pdo = [];
16+pdo.lap = 1;
17+pdo.c = @(x, y, z) 300*(1 - z);
18+
19+f = -1;
20+L = surfaceop(dom, pdo, f);
21+u = L.solve();
22+
23+figure(1);
24+plot(u), colorbar
25+title('Variable-coefficient Helmholtz on a stellarator');
pde_laplace_beltrami.madded+37−0View file
@@ -0,0 +1,37 @@
1+% pde_laplace_beltrami.m — Solving a PDE on a surface (surfaceop).
2+%
3+% surfaceop is a fast direct solver for second-order linear elliptic PDEs
4+% on surfaces. Here we solve a Laplace-Beltrami problem on the sphere,
5+% where the spherical harmonics give an exact solution to check against.
6+
7+mip load --install flatironinstitute/flatironinstitute/surfacefun
8+
9+p = 16;
10+nref = 2;
11+dom = surfacemesh.sphere(p + 1, nref);
12+
13+% Exact solution: a spherical harmonic Y_lm is an eigenfunction of the
14+% Laplace-Beltrami operator with eigenvalue -l(l+1).
15+l = 3;
16+m = 2;
17+sol = spherefun.sphharm(l, m);
18+sol = surfacefun(@(x, y, z) sol(x, y, z), dom);
19+f = -l*(l + 1)*sol;
20+
21+% Define the operator (pdo.lap = 1 is the Laplace-Beltrami operator).
22+pdo = [];
23+pdo.lap = 1;
24+
25+L = surfaceop(dom, pdo, f);
26+
27+% The closed-surface problem is rank-deficient by one; impose the
28+% mean-zero condition.
29+L.rankdef = true;
30+
31+u = L.solve();
32+
33+figure(1);
34+plot(u), colorbar
35+title('Laplace-Beltrami solution on the sphere');
36+
37+fprintf('error || u - sol || = %.3e\n', norm(u - sol));
pde_open_surface.madded+36−0View file
@@ -0,0 +1,36 @@
1+% pde_open_surface.m — A PDE on an open surface, with boundary data.
2+%
3+% On an open surface, surfaceop.solve takes Dirichlet boundary data. And
4+% because surfaceop is a direct solver, one factorization can be reused
5+% for many right-hand sides — much cheaper than rebuilding it each time.
6+
7+mip load --install flatironinstitute/flatironinstitute/surfacefun
8+
9+% Build an open surface by keeping a subset of a blob's patches.
10+rng(0)
11+p = 16;
12+nref = 2;
13+dom = surfacemesh.blob(p + 1, nref);
14+dom = surfacemesh(dom.x(1:16), dom.y(1:16), dom.z(1:16));
15+
16+figure(1);
17+plot(dom), view(-110, 30), camlight
18+title('Open surface (subset of a blob)');
19+
20+% Solve Delta_Gamma u = -1 with zero Dirichlet data on the boundary.
21+pdo = [];
22+pdo.lap = 1;
23+L = surfaceop(dom, pdo, -1);
24+u = L.solve(0);
25+
26+figure(2);
27+plot(u), view(-110, 30), colorbar
28+title('Solution, rhs = -1, u = 0 on boundary');
29+
30+% Reuse the same factorization with a new right-hand side.
31+L.rhs = @(x, y, z) sin(x .* y);
32+u = L.solve(0);
33+
34+figure(3);
35+plot(u), view(-110, 30), colorbar
36+title('Reused solver, rhs = sin(xy)');
vec_curl.madded+18−0View file
@@ -0,0 +1,18 @@
1+% vec_curl.m — The surface curl of a vector field.
2+%
3+% curl of a surfacefunv is again a surfacefunv.
4+
5+mip load --install flatironinstitute/flatironinstitute/surfacefun
6+
7+rng(0)
8+p = 16;
9+nref = 2;
10+dom = surfacemesh.blob(p + 1, nref);
11+
12+v = surfacefunv(@(x, y, z) cos(2*x), ...
13+ @(x, y, z) sin(4*y), ...
14+ @(x, y, z) sin(3*z), dom);
15+
16+figure(1);
17+quiver(curl(v), 0.1, 6)
18+title('Surface curl of a vector field');
vec_gradient.madded+20−0View file
@@ -0,0 +1,20 @@
1+% vec_gradient.m — The surface gradient as a vector field.
2+%
3+% grad(f) of a surfacefun is a surfacefunv. The gradient is tangent to
4+% the surface, as the quiver plot shows.
5+
6+mip load --install flatironinstitute/flatironinstitute/surfacefun
7+
8+rng(0)
9+p = 16;
10+nref = 2;
11+dom = surfacemesh.blob(p + 1, nref);
12+
13+f = surfacefun(@(x, y, z) cos(6*x) .* y + exp(z), dom);
14+
15+figure(1);
16+quiver(grad(f), 0.05, 6)
17+title('Surface gradient of f (tangent to the surface)');
18+
19+% Identity: the divergence of the gradient is the Laplacian.
20+fprintf('|| div(grad f) - lap(f) || = %.3e\n', norm(div(grad(f)) - lap(f)));
vec_mean_curvature.madded+17−0View file
@@ -0,0 +1,17 @@
1+% vec_mean_curvature.m — Mean curvature from the normal field.
2+%
3+% The mean curvature of a surface is half the surface divergence of its
4+% unit normal field: H = div(n)/2.
5+
6+mip load --install flatironinstitute/flatironinstitute/surfacefun
7+
8+rng(0)
9+p = 16;
10+nref = 2;
11+dom = surfacemesh.blob(p + 1, nref);
12+
13+v = normal(dom);
14+
15+figure(1);
16+plot(div(v)/2), colorbar
17+title('Mean curvature div(n)/2');
vec_normals.madded+18−0View file
@@ -0,0 +1,18 @@
1+% vec_normals.m — The surface normal field (surfacefunv).
2+%
3+% A surfacefunv is a three-component vector field on a surface. normal(dom)
4+% returns the unit outward normal. quiver(v, scale, n) draws n vectors per
5+% patch, scaled by `scale`.
6+
7+mip load --install flatironinstitute/flatironinstitute/surfacefun
8+
9+rng(0)
10+p = 16;
11+nref = 2;
12+dom = surfacemesh.blob(p + 1, nref);
13+
14+v = normal(dom);
15+
16+figure(1);
17+quiver(v, 0.2, 6)
18+title('Surface normals');