///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================
//! @examplefile mosolov_augmented_lagrangian1.icc The Mossolov problem by the augmented Lagrangian method -- solver class body
#include "vector_projection.h"
int mosolov_augmented_lagrangian::solve (field& sigma_h, field& uh) const {
  test v(Xh);
  derr << "# k residue" << endl;
  for (size_t k = 0; true; ++k) {
    field grad_uh = inv_mt*(b*uh);
    auto c = compose(vector_projection(Bi,n,1,r), norm(sigma_h+r*grad_uh));
    field gamma_h = interpolate(Th, c*(sigma_h + r*grad_uh));
    field delta_sigma_h = r*(grad_uh - gamma_h);
    sigma_h += delta_sigma_h;
    Float residue = delta_sigma_h.max_abs();
    derr << k << " " << residue << endl;
    if (residue <= tol || k >= max_iter) {
      derr << endl << endl;
      return (pow(residue,3) <= tol) ? 0 : 1;
    }
    field rhs = (1/r)*(lh - integrate(dot(sigma_h - r*gamma_h, grad(v))));
    pa.solve (rhs, uh);
  }
}
