% pde_laplace_beltrami.m — Solving a PDE on a surface (surfaceop).
%
% surfaceop is a fast direct solver for second-order linear elliptic PDEs
% on surfaces. Here we solve a Laplace-Beltrami problem on the sphere,
% where the spherical harmonics give an exact solution to check against.

mip load --install flatironinstitute/flatironinstitute/surfacefun

p = 16;
nref = 2;
dom = surfacemesh.sphere(p + 1, nref);

% Exact solution: a spherical harmonic Y_lm is an eigenfunction of the
% Laplace-Beltrami operator with eigenvalue -l(l+1).
l = 3;
m = 2;
sol = spherefun.sphharm(l, m);
sol = surfacefun(@(x, y, z) sol(x, y, z), dom);
f = -l*(l + 1)*sol;

% Define the operator (pdo.lap = 1 is the Laplace-Beltrami operator).
pdo = [];
pdo.lap = 1;

L = surfaceop(dom, pdo, f);

% The closed-surface problem is rank-deficient by one; impose the
% mean-zero condition.
L.rankdef = true;

u = L.solve();

figure(1);
plot(u), colorbar
title('Laplace-Beltrami solution on the sphere');

fprintf('error || u - sol || = %.3e\n', norm(u - sol));