navier_stokes_dg::navier_stokes_dg (
 Float Re1, const geo& omega, string approx)
 : Re(Re1), Xh(), Qh(), qopt(), a0(), b(), c(), mu(), mp(), lh0(), lh(), kh(), 
   smu(), smp(), a1(), stokes1()
{
  Xh = space (omega, approx, "vector");
  Qh = space (omega, approx);
  qopt.set_family(quadrature_option_type::gauss);
  qopt.set_order(2*Xh.degree()+1);
  stokes_dirichlet_dg (Xh, Qh, a0, b, c, mp, lh0, kh, qopt);
  trial u (Xh); test v (Xh);
  lh = lh0 + Re*inertia_fix_rhs (v, qopt);
  mu = integrate (dot(u,v), qopt);
  smu = solver(mu.uu());
  smp = solver(mp.uu());
}
navier_stokes_dg::value_type
navier_stokes_dg::initial (string restart) const {
  value_type xh(2);
  xh[0] = field (Xh, 0);
  xh[1] = field (Qh, 0);
  Float Re0 = 0;
  if (restart == "") {
    solver_abtb stokes0 (a0.uu(), b.uu(), c.uu(), mp.uu());
    stokes0.solve (lh0.u(), kh.u(), xh[0].set_u(), xh[1].set_u());
  } else {
    idiststream in (restart);
    in >> catchmark("Re") >> Re0
       >> catchmark("u")  >> xh[0]
       >> catchmark("p")  >> xh[1];
    check_macro (xh[1].get_space() == Qh, "unexpected " << xh[0].get_space().stamp()
	<< " approximation in file \"" << restart << "\" (" << Xh.stamp() << " expected)");
  }
  derr << "# continuation: from Re=" << Re0 << " to " << Re << endl;
  return xh;
}
navier_stokes_dg::value_type
navier_stokes_dg::residue (const value_type& xh) const {
  trial u (Xh); test v (Xh);
  form a = a0 + Re*inertia(xh[0], u, v, qopt);
  value_type mrh(2);
  mrh[0] = a*xh[0] + b.trans_mult(xh[1]) - lh;
  mrh[1] = b*xh[0] -            c*xh[1]  - kh;
  return mrh;
}
void navier_stokes_dg::update_derivative (const value_type& xh) const {
  trial u (Xh); test v (Xh);
  a1 = a0 + Re*(inertia(xh[0], u, v, qopt) + inertia(u, xh[0], v, qopt));
  stokes1 = solver_abtb (a1.uu(), b.uu(), c.uu(), mp.uu());
}
navier_stokes_dg::value_type
navier_stokes_dg::derivative_solve (const value_type& mrh) const {
  value_type delta_xh(2);
  delta_xh[0] = field (Xh, 0);
  delta_xh[1] = field (Qh, 0);
  stokes1.solve (mrh[0].u(), mrh[1].u(), 
	         delta_xh[0].set_u(), delta_xh[1].set_u());
  return delta_xh;
}
navier_stokes_dg::value_type
navier_stokes_dg::derivative_trans_mult (const value_type& mrh) const {
  value_type rh(2);
  rh[0] = field (Xh);
  rh[1] = field (Qh);
  rh[0].set_u() = smu.solve(mrh[0].u());
  rh[1].set_u() = smp.solve(mrh[1].u());
  value_type mgh(2);
  mgh[0] = a1.trans_mult(rh[0]) + b.trans_mult(rh[1]);
  mgh[1] =             b*rh[0]  -            c*rh[1];
  return mgh;
}
