1function result = solve_pde(mshfile, params)
2%SOLVE_PDE Load the quad mesh and solve the selected PDE on it.
3% params fields (from the host UI):
4% pde - 'poisson' (lap u = f) or 'helmholtz' ((lap + c) u = f)
5% f - right-hand side, a MATLAB expression in x, y, z
6% c - zeroth-order coefficient expression (helmholtz only)
7% p - polynomial order per patch
8% closed - true if every mesh edge is shared by exactly two quads
9% (determined host-side from the connectivity)
11dom = surfacemesh.import(mshfile, 'gmsh');
12dom = resample(dom, params.p + 1);
14fh = eval(['@(x, y, z) ', params.f]);
15f = surfacefun(@(x, y, z) fh(x, y, z) + 0*x, dom);
17pdo = [];
18pdo.lap = 1;
19isPoisson = strcmp(params.pde, 'poisson');
20if ~isPoisson
21 ch = eval(['@(x, y, z) ', params.c]);
22 pdo.c = @(x, y, z) ch(x, y, z) + 0*x;
23end
25closed = params.closed;
26if isPoisson && closed
27 % The closed-surface Laplace-Beltrami problem is rank-deficient by one
28 % and only solvable for mean-zero data; project the RHS accordingly.
29 f = f - mean(f);
30end
32L = surfaceop(dom, pdo, f);
33if closed
34 if isPoisson
35 L.rankdef = true;
36 end
37 u = L.solve();
38else
39 u = L.solve(0); % zero Dirichlet boundary data on open surfaces
40end
42result = pack_solution(dom, u);
43result.pde = params.pde;
44end
46function data = pack_solution(dom, u)
47% One flat (column-major) x/y/z/u array per patch, each an n-by-n grid —
48% the same layout surfacefun-interactive's figure app uses.
49np = length(dom);
50px = cell(1, np);
51py = cell(1, np);
52pz = cell(1, np);
53pu = cell(1, np);
54umin = inf;
55umax = -inf;
56for k = 1:np
57 px{k} = real(dom.x{k}(:).');
58 py{k} = real(dom.y{k}(:).');
59 pz{k} = real(dom.z{k}(:).');
60 vals = real(u.vals{k}(:).');
61 pu{k} = vals;
62 umin = min(umin, min(vals));
63 umax = max(umax, max(vals));
64end
65data = struct();
66data.type = 'solution';
67data.n = size(dom.x{1}, 1);
68data.npatches = np;
69data.x = px;
70data.y = py;
71data.z = pz;
72data.u = pu;
73data.umin = umin;
74data.umax = umax;
75end