#include "eta.icc"
#include "nu.icc"
#include "dirichlet.icc"
p_laplacian::p_laplacian (Float p1, const geo& omega, string approx)
 : p(p1), Xh(), lh(), m(), sm(), qopt(), a1(), sa1() {
  Xh = space (omega, approx);
  Xh.block ("boundary");
  qopt.set_family(quadrature_option_type::gauss);
  qopt.set_order(2*Xh.degree()-1);
  trial u (Xh); test v (Xh);
  lh = integrate (v);
  m = integrate (u*v);
  sm = solver (m.uu());
}
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)),
                      qopt);
  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)),
                  qopt);
  sa1 = ldlt (a1.uu());
}
field p_laplacian::derivative_solve (const field& rh) const {
  field delta_uh (Xh,0);
  delta_uh.set_u() = sa1.solve(rh.u());
  return delta_uh;
}
