1function result = solve_pde(mshfile, params)
2%SOLVE_PDE Load the surface 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 cells
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 x/y/z/u array per patch: a column-major n-by-n grid for quad
48% patches (the layout surfacefun-interactive's figure app uses), or the
49% n*(n+1)/2-point trianglepts(n) set for triangle patches. data.n is the
50% number of points per patch edge in both cases.
51np = length(dom);
52px = cell(1, np);
53py = cell(1, np);
54pz = cell(1, np);
55pu = cell(1, np);
56umin = inf;
57umax = -inf;
58for k = 1:np
59 px{k} = real(dom.x{k}(:).');
60 py{k} = real(dom.y{k}(:).');
61 pz{k} = real(dom.z{k}(:).');
62 vals = real(u.vals{k}(:).');
63 pu{k} = vals;
64 umin = min(umin, min(vals));
65 umax = max(umax, max(vals));
66end
67data = struct();
68data.type = 'solution';
69if ( dom.ptype(1) == surfacemesh.patchtype.tri )
70 npts = length(dom.x{1});
71 data.n = round((sqrt(8*npts + 1) - 1) / 2);
72 data.ptype = 'tri';
73else
74 data.n = size(dom.x{1}, 1);
75 data.ptype = 'quad';
76end
77data.npatches = np;
78data.x = px;
79data.y = py;
80data.z = pz;
81data.u = pu;
82data.umin = umin;
83data.umax = umax;
84end