% pde_open_surface.m — A PDE on an open surface, with boundary data.
%
% On an open surface, surfaceop.solve takes Dirichlet boundary data. And
% because surfaceop is a direct solver, one factorization can be reused
% for many right-hand sides — much cheaper than rebuilding it each time.

mip load --install flatironinstitute/flatironinstitute/surfacefun

% Build an open surface by keeping a subset of a blob's patches.
rng(0)
p = 16;
nref = 2;
dom = surfacemesh.blob(p + 1, nref);
dom = surfacemesh(dom.x(1:16), dom.y(1:16), dom.z(1:16));

figure(1);
plot(dom), view(-110, 30), camlight
title('Open surface (subset of a blob)');

% Solve Delta_Gamma u = -1 with zero Dirichlet data on the boundary.
pdo = [];
pdo.lap = 1;
L = surfaceop(dom, pdo, -1);
u = L.solve(0);

figure(2);
plot(u), view(-110, 30), colorbar
title('Solution, rhs = -1, u = 0 on boundary');

% Reuse the same factorization with a new right-hand side.
L.rhs = @(x, y, z) sin(x .* y);
u = L.solve(0);

figure(3);
plot(u), view(-110, 30), colorbar
title('Reused solver, rhs = sin(xy)');