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)
11msh = load_gmsh_quads(mshfile);
12dom = surfacemesh_from_quads(msh);
13dom = resample(dom, params.p + 1);
15fh = eval(['@(x, y, z) ', params.f]);
16f = surfacefun(@(x, y, z) fh(x, y, z) + 0*x, dom);
18pdo = [];
19pdo.lap = 1;
20isPoisson = strcmp(params.pde, 'poisson');
21if ~isPoisson
22 ch = eval(['@(x, y, z) ', params.c]);
23 pdo.c = @(x, y, z) ch(x, y, z) + 0*x;
24end
26closed = params.closed;
27if isPoisson && closed
28 % The closed-surface Laplace-Beltrami problem is rank-deficient by one
29 % and only solvable for mean-zero data; project the RHS accordingly.
30 f = f - mean(f);
31end
33L = surfaceop(dom, pdo, f);
34if closed
35 if isPoisson
36 L.rankdef = true;
37 end
38 u = L.solve();
39else
40 u = L.solve(0); % zero Dirichlet boundary data on open surfaces
41end
43result = pack_solution(dom, u);
44result.pde = params.pde;
45end
47function data = pack_solution(dom, u)
48% One flat (column-major) x/y/z/u array per patch, each an n-by-n grid —
49% the same layout surfacefun-interactive's figure app uses.
50np = length(dom);
51px = cell(1, np);
52py = cell(1, np);
53pz = cell(1, np);
54pu = cell(1, np);
55umin = inf;
56umax = -inf;
57for k = 1:np
58 px{k} = real(dom.x{k}(:).');
59 py{k} = real(dom.y{k}(:).');
60 pz{k} = real(dom.z{k}(:).');
61 vals = real(u.vals{k}(:).');
62 pu{k} = vals;
63 umin = min(umin, min(vals));
64 umax = max(umax, max(vals));
65end
66data = struct();
67data.type = 'solution';
68data.n = size(dom.x{1}, 1);
69data.npatches = np;
70data.x = px;
71data.y = py;
72data.z = pz;
73data.u = pu;
74data.umin = umin;
75data.umax = umax;
76end