// Atsushi Suzuki (c) 5 Nov. 2025
// CG method with orthogonal projection onto Im A
// 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;
real tmp = sqrt(one' * one);
one = one / tmp;

int count = 0; 
func real[int] Proj(real[int] &uu){
  real utmp = one' * uu;
  uu = uu - utmp * one;
  return uu;
}

func real[int] OpA(real[int] &uu) {
  real[int] vv(uu.n);
  uu = Proj(uu);
  vv = A * uu;
  uu = Proj(vv);
  return uu;
}

real[int] ff = Poisson(0, Vh);
Vh u;
u[] = 0.0;
ff = Proj(ff);
LinearCG(OpA, u[], ff, eps=1.e-12, nbiter=1000, verbosity=10);
plot(Th, u, cmm="with projection", wait = true);
tmp = one' * u[];
cout << "solution satisfies orthogonality as " << tmp << endl;