// Atsushi Suzuki (c) 5 Nov. 2025
// direct solver for enhanced matirx with Lagrange multiplier 
// external force satisfies int2d(Th)(f) = 0.0 but with perturbation
real delta = 1e-3;
func f = -5.0 * pi * pi * cos(pi * x) * cos(2.0 * pi * y) + delta;
int nn = 100;
mesh Th = square(nn, nn);

fespace Vh(Th, P1);

varf Poisson(u, v) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v)) + int2d(Th)(f * v);

matrix A = Poisson(Vh, Vh);
real[int] one(Vh.ndof);
one = 1.0;
matrix K = [[A, one], [one', 0.0]];
set (K, solver = sparsesolver);
real[int] ff = Poisson(0, Vh);
real tmp = one' * ff;
cout << "f belongs to image space? (1, f)="<< tmp << endl;
real[int] bb(Vh.ndof+1), uu(Vh.ndof+1);
bb = [ff, 0.0];
uu = K^-1 * bb;
Vh u;
[u[], tmp] = uu;
cout << "Lagrange multiplier = " << tmp << endl;
plot(Th, u, cmm="without projection", wait = true);
// RHS is forced in the image space
tmp = one' * ff;
tmp = tmp / (one' * one); // normalize without changing one vector
ff = ff - tmp * one; // projection onto image space = orthogonal to span[1]
tmp = one' * ff;
bb = [ff, 0.0];
uu = K^-1 * bb;
[u[], tmp] = uu;
cout << "Lagrange multiplier = " << tmp << endl;
plot(Th, u, cmm="with projection", wait = true);