// FullNeumannPoisson-CG-v0.edp
// modified by Masashi Katsurada on 2026/6/23.
// 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

int withprojection = 0;

func f = -5.0 * pi * pi * cos(pi * x) * cos(2.0 * pi * y);
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);

// the projection to Image A
real[int] one(Vh.ndof);
one = 1.0;
one = one / sqrt(one' * one);

func real[int] Proj(real[int] &uu){
  real utmp = one' * uu;
  //cout << "(1,u)=" << utmp << endl;
  uu = uu - utmp * one;
  return uu;
}

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

Vh u;
u[] = 0.0;
if (withprojection) ff = Proj(ff);
LinearCG(OpA, u[], ff, eps=1.e-12, nbiter=1000, verbosity=10);
if (withprojection)
  plot(u, cmm="with projection", wait = true);
else
  plot(u, cmm="without projection", wait = true);

real tmp = one' * u[];
cout << "solution satisfies orthogonality as " << tmp << endl;