1% Solve a PDE on a surface mesh with surfacefun
2% (https://github.com/danfortunato/surfacefun).
3%
4% Generated by mesh-pde-solver: the parameter values below come from the UI,
5% and this exact script is what the Solve button runs in the browser (via
6% numbl). To run it in desktop MATLAB, keep the mesh file (the app's
7% converted .msh download) next to this script; without the mip package
8% manager, comment out the `mip load` line and put surfacefun and its
9% dependencies on the path yourself.
11mip load --install surfacefun;
13% ---- parameters -----------------------------------------------------------
14meshfile = '{{MESHFILE}}';
15pde = '{{PDE}}'; % 'poisson' or 'helmholtz'
16frhs = @(x, y, z) {{F_EXPR}}; % right-hand side f
17crhs = @(x, y, z) {{C_EXPR}}; % zeroth-order coefficient c (helmholtz only)
18p = {{ORDER}}; % polynomial order per patch
19% With the parameters baked in as constants, one branch of each pde/closed
20% `if` below is statically dead — suppress the resulting lint file-wide.
21closed = {{CLOSED}}; %#ok<*BDLGI,*BDSCI,*UNRCH> true if the mesh is closed
22% ----------------------------------------------------------------------------
24dom = surfacemesh.import(meshfile, 'gmsh');
25dom = resample(dom, p + 1);
27f = surfacefun(@(x, y, z) frhs(x, y, z) + 0*x, dom);
29pdo = [];
30pdo.lap = 1;
31isPoisson = strcmp(pde, 'poisson');
32if ~isPoisson
33 pdo.c = @(x, y, z) crhs(x, y, z) + 0*x;
34end
36if isPoisson && closed
37 % The closed-surface Laplace-Beltrami problem is rank-deficient by one
38 % and only solvable for mean-zero data; project the RHS accordingly.
39 f = f - mean(f);
40end
42L = surfaceop(dom, pdo, f);
43if closed
44 if isPoisson
45 L.rankdef = true;
46 end
47 u = L.solve();
48else
49 u = L.solve(0); % zero Dirichlet boundary data on open surfaces
50end
52% In desktop MATLAB, visualize the solution with: plot(u), colorbar
54% ---- pack the solution as result.json for the host app --------------------
55% One flat x/y/z/u array per patch: a column-major n-by-n grid for quad
56% patches, or the n*(n+1)/2-point trianglepts(n) set for triangle patches.
57% result.n is the number of points per patch edge in both cases.
58np = length(dom);
59px = cell(1, np);
60py = cell(1, np);
61pz = cell(1, np);
62pu = cell(1, np);
63umin = inf;
64umax = -inf;
65for k = 1:np
66 px{k} = real(dom.x{k}(:).');
67 py{k} = real(dom.y{k}(:).');
68 pz{k} = real(dom.z{k}(:).');
69 vals = real(u.vals{k}(:).');
70 pu{k} = vals;
71 umin = min(umin, min(vals));
72 umax = max(umax, max(vals));
73end
74result = struct();
75result.type = 'solution';
76if ( dom.ptype(1) == surfacemesh.patchtype.tri )
77 npts = length(dom.x{1});
78 result.n = round((sqrt(8*npts + 1) - 1) / 2);
79 result.ptype = 'tri';
80else
81 result.n = size(dom.x{1}, 1);
82 result.ptype = 'quad';
83end
84result.npatches = np;
85result.x = px;
86result.y = py;
87result.z = pz;
88result.u = pu;
89result.umin = umin;
90result.umax = umax;
91result.pde = pde;
93fid = fopen('result.json', 'w');
94fprintf(fid, '%s', jsonencode(result));
95fclose(fid);