// FullNeumannPoisson-LM-v1.edp
// modified by Masashi Katsurada on 2026/6/23.
// Atsushi Suzuki (c) 5 Nov. 2025
// direct solver for enhanced matirx with Lagrange multiplier 
// external force satisfies int2d(Th)(f) = 0.0

real perturbation = 0; // 1e-3;
func f = -5.0 * pi * pi * cos(pi * x) * cos(2.0 * pi * y) + perturbation;
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] ff = Poisson(0, Vh);

real[int] one(Vh.ndof);
one = 1.0;
cout << "Does f belong to image space? (1, f)="<< (one' * ff) << endl;

matrix K = [[A, one], [one', 0.0]];
set (K, solver = sparsesolver);

real[int] bb(Vh.ndof+1), uu(Vh.ndof+1);
bb = [ff, 0.0];
uu = K^-1 * bb;

Vh u;
real lambda;
[u[], lambda] = uu;
cout << "without projection" << endl;
cout << "Lagrange multiplier (should be equal to 0) = " << lambda << endl;

plot(u, cmm="without projection", wait = true);

// RHS is forced in the image space
real 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 << endl << "with projection" << endl;
cout << "Lagrange multiplier (should be equal to 0) = " << tmp << endl;
plot(u, cmm="with projection", wait = true);