///
/// 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_lagrangian2.icc The Mossolov problem by the augmented Lagrangian method -- solver class body
mosolov_augmented_lagrangian::mosolov_augmented_lagrangian()
: Bi(0), n(1), r(1), tol(1e-10), max_iter(1000000), 
  Xh(), Th(), lh(), a(), b(), inv_mt(), pa()
{}
void mosolov_augmented_lagrangian::reset (geo omega, string approx) {
  Xh = space (omega, approx);
  Xh.block ("boundary");
  string grad_approx = "P" + itos(Xh.degree()-1) + "d";
  Th = space (omega, grad_approx, "vector");
  trial u (Xh), sigma(Th);
  test  v (Xh), tau  (Th);
  lh = integrate(2*v);
  a = integrate (dot(grad(u),grad(v)));
  b = integrate (dot(grad(u),tau));
  integrate_option iopt;
  iopt.invert = true;
  inv_mt = integrate(dot(sigma,tau), iopt);
  pa = problem (a);
}
void
mosolov_augmented_lagrangian::initial (field& sigma_h, field& uh) const {
  uh = field(Xh);
  uh ["boundary"] = 0;
  pa.solve (lh, uh);
  test tau (Th);
  field mt_grad_uh = integrate(dot(grad(uh),tau)); 
  sigma_h = inv_mt*mt_grad_uh;
}
void mosolov_augmented_lagrangian::put (odiststream& out, 
    field& sigma_h, field& uh) const
{
  out << catchmark("Bi")    << Bi << endl
      << catchmark("n")     << n  << endl
      << catchmark("r")     << r  << endl
      << catchmark("sigma") << sigma_h
      << catchmark("u")     << uh;
}
