///
/// 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 navier_stokes_dg1.icc The Navier-Stokes equations with the discontinuous Galerkin method -- class body
navier_stokes_dg::navier_stokes_dg (
 Float Re1, const geo& omega, string approx)
 : Re(Re1), Xh(), Qh(), iopt(), a0(), b(), c(), mu(), mp(), lh0(), lh(), kh(), 
   pmu(), pmp(), a1(), stokes1()
{
  Xh = space (omega, approx, "vector");
  Qh = space (omega, approx);
  iopt.set_family(integrate_option::gauss);
  iopt.set_order(2*Xh.degree()+1);
  stokes_dirichlet_dg (Xh, Qh, a0, b, c, mp, lh0, kh, iopt);
  trial u (Xh); test v (Xh);
  lh = lh0 + Re*inertia_fix_rhs (v, iopt);
  mu = integrate (dot(u,v), iopt);
  pmu = problem (mu);
  pmp = problem (mp);
}
navier_stokes_dg::value_type
navier_stokes_dg::initial (string restart) const {
  value_type xh = { field(Xh, 0), field(Qh, 0) };
  Float Re0 = 0;
  if (restart == "") {
    problem_mixed stokes0 (a0, b, c);
    stokes0.set_metric (mp);
    stokes0.solve (lh0, kh, xh[0], xh[1]);
  } 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().name() << " approximation in file \""
       << restart << "\" (" << Xh.name() << " 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, iopt);
  value_type mrh = { a*xh[0] + b.trans_mult(xh[1]) - lh,
                     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, iopt) + inertia(u, xh[0], v, iopt));
  stokes1 = problem_mixed (a1, b, c);
  stokes1.set_metric (mp);
}
navier_stokes_dg::value_type
navier_stokes_dg::derivative_solve (const value_type& mrh) const {
  value_type delta_xh = { field(Xh, 0), field(Qh, 0) };
  stokes1.solve (mrh[0], mrh[1], delta_xh[0], delta_xh[1]);
  return delta_xh;
}
navier_stokes_dg::value_type
navier_stokes_dg::derivative_trans_mult (const value_type& mrh) const {
  value_type rh = { field (Xh), field (Qh) };
  pmu.solve (mrh[0], rh[0]);
  pmp.solve (mrh[1], rh[1]);
  value_type mgh = { a1.trans_mult(rh[0]) + b.trans_mult(rh[1]),
                     b*rh[0]  -            c*rh[1]             };
  return mgh;
}
