///
/// 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 burgers_diffusion_operators.icc The diffusive Burgers equation -- operators
field lh (Float epsilon, Float t, const test& v) {
#ifdef  NEUMANN
  return field (v.get_vf_space(), 0);
#else //  NEUMANN
  size_t d = v.get_vf_space().get_geo().dimension();
  size_t k = v.get_vf_space().degree();
  Float beta = (k+1)*(k+d)/Float(d);
  return epsilon*integrate ("boundary", 
         beta*penalty()*g(epsilon,t)*v
       - g(epsilon,t)*dot(grad_h(v),normal()));
#endif //  NEUMANN
}
field gh (Float epsilon, Float t, const field& uh, const test& v) {
  return - integrate (dot(compose(f,uh),grad_h(v)))
         + integrate ("internal_sides",
             compose (phi, normal(), inner(uh), outer(uh))*jump(v))
         + integrate ("boundary",
             compose (phi, normal(), uh, g(epsilon,t))*v);
}
