% Solve a PDE on a surface mesh with surfacefun
% (https://github.com/danfortunato/surfacefun).
%
% Generated by mesh-pde-solver: the parameter values below come from the UI,
% and this exact script is what the Solve button runs in the browser (via
% numbl). To run it in desktop MATLAB, keep the mesh file (the app's
% converted .msh download) next to this script; without the mip package
% manager, comment out the `mip load` line and put surfacefun and its
% dependencies on the path yourself.

mip load --install surfacefun;

% ---- parameters -----------------------------------------------------------
meshfile = '{{MESHFILE}}';
pde      = '{{PDE}}';              % 'poisson' or 'helmholtz'
frhs     = @(x, y, z) {{F_EXPR}};  % right-hand side f
crhs     = @(x, y, z) {{C_EXPR}};  % zeroth-order coefficient c (helmholtz only)
p        = {{ORDER}};              % polynomial order per patch
% With the parameters baked in as constants, one branch of each pde/closed
% `if` below is statically dead — suppress the resulting lint file-wide.
closed   = {{CLOSED}};             %#ok<*BDLGI,*BDSCI,*UNRCH> true if the mesh is closed
% ----------------------------------------------------------------------------

dom = surfacemesh.import(meshfile, 'gmsh');
dom = resample(dom, p + 1);

f = surfacefun(@(x, y, z) frhs(x, y, z) + 0*x, dom);

pdo = [];
pdo.lap = 1;
isPoisson = strcmp(pde, 'poisson');
if ~isPoisson
    pdo.c = @(x, y, z) crhs(x, y, z) + 0*x;
end

if isPoisson && closed
    % The closed-surface Laplace-Beltrami problem is rank-deficient by one
    % and only solvable for mean-zero data; project the RHS accordingly.
    f = f - mean(f);
end

L = surfaceop(dom, pdo, f);
if closed
    if isPoisson
        L.rankdef = true;
    end
    u = L.solve();
else
    u = L.solve(0);   % zero Dirichlet boundary data on open surfaces
end

% In desktop MATLAB, visualize the solution with:  plot(u), colorbar

% ---- pack the solution as result.json for the host app --------------------
% One flat x/y/z/u array per patch: a column-major n-by-n grid for quad
% patches, or the n*(n+1)/2-point trianglepts(n) set for triangle patches.
% result.n is the number of points per patch edge in both cases.
np = length(dom);
px = cell(1, np);
py = cell(1, np);
pz = cell(1, np);
pu = cell(1, np);
umin = inf;
umax = -inf;
for k = 1:np
    px{k} = real(dom.x{k}(:).');
    py{k} = real(dom.y{k}(:).');
    pz{k} = real(dom.z{k}(:).');
    vals = real(u.vals{k}(:).');
    pu{k} = vals;
    umin = min(umin, min(vals));
    umax = max(umax, max(vals));
end
result = struct();
result.type = 'solution';
if ( dom.ptype(1) == surfacemesh.patchtype.tri )
    npts = length(dom.x{1});
    result.n = round((sqrt(8*npts + 1) - 1) / 2);
    result.ptype = 'tri';
else
    result.n = size(dom.x{1}, 1);
    result.ptype = 'quad';
end
result.npatches = np;
result.x = px;
result.y = py;
result.z = pz;
result.u = pu;
result.umin = umin;
result.umax = umax;
result.pde = pde;

fid = fopen('result.json', 'w');
fprintf(fid, '%s', jsonencode(result));
fclose(fid);