///
/// 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 yield_slip_augmented_lagrangian.icc The yield slip problem by the augmented Lagrangian method -- solver function
#include "projection.h"
int yield_slip_augmented_lagrangian (Float S, Float n, Float Cf,
    geo boundary, field lh, field& lambda_h, field& uh,
    Float tol, size_t max_iter, Float r)
{
  const space& Xh = uh.get_space();
  const space& Wh = lambda_h.get_space();
  trial u(Xh), lambda(Wh);
  test  v(Xh),     mu(Wh);
  form m  = integrate (u*v),
       a0 = integrate (dot(grad(u),grad(v))),
       a  = a0 + integrate (boundary, r*u*v),
       mb = integrate (lambda*mu),
       b  = integrate (boundary, u*mu);
  problem pa (a);
  derr << "# k residue" << endl;
  Float residue0 = 0;
  for (size_t k = 0; true; ++k) {
    field gamma_h = interpolate(Wh,
      compose(projection(S,n,Cf,r), lambda_h + r*uh["boundary"]));
    field delta_lambda_h = r*(uh["boundary"] - gamma_h);
    lambda_h += delta_lambda_h;
    Float residue = delta_lambda_h.max_abs();
    derr << k << " " << residue << endl;
    if (residue <= tol || k >= max_iter) return (residue <= tol) ? 0 : 1;
    field rhs = lh + b.trans_mult(r*gamma_h - lambda_h); 
    pa.solve (rhs, uh);
  }
}
