///
/// 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 p_laplacian1.icc The p-Laplacian problem by the Newton method -- class body
#include "eta.h"
#include "nu.h"
#include "dirichlet.icc"
p_laplacian::p_laplacian (Float p1, const geo& omega, string approx)
 : p(p1), Xh(), lh(), m(), pm(), a1(), pa1() {
  Xh = space (omega, approx);
  Xh.block ("boundary");
  trial u (Xh); test v (Xh);
  lh = integrate (v);
  m = integrate (u*v);
  pm = problem (m);
}
field p_laplacian::initial() const {
  field uh (Xh, 0);
  dirichlet (lh, uh);
  return uh;
}
field p_laplacian::residue (const field& uh) const {
  trial u (Xh); test v (Xh);
  form a = integrate (compose(eta(p), norm2(grad(uh)))*dot(grad(u),grad(v)));
  field mrh = a*uh - lh;
  mrh.set_b() = 0;
  return mrh;
}
void p_laplacian::update_derivative (const field& uh) const {
  size_t d = Xh.get_geo().dimension();
  trial u (Xh); test v (Xh);
  a1 = integrate (dot(compose(nu<eta>(eta(p),d), grad(uh))*grad(u),grad(v)));
  pa1 = problem (a1);
}
field p_laplacian::derivative_solve (const field& rh) const {
  field delta_uh (Xh,0);
  pa1.solve (rh, delta_uh);
  return delta_uh;
}
