///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2018 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 transport_tensor_exact.icc The tensorial transport benchmark -- right-hand-side and exact solution
point u(const point& x) { return point (-x[1], x[0]); }
tensor grad_u = {{0,-1},{1,0}};

struct sigma_exact {  
 Float f (const point& x) const {
    return 0.5*exp( - nu*t - norm2(x-xt)/sqr(r0));
 }
 Float df_dt (const point& x) const {
    return ( - nu + 2*dot(vt,x-xt)/sqr(r0))*f(x);
 }
 tensor operator() (const point& x) const {
    Float c0 = f(x);
    tensor s;
    s(0,0) = c0*(1+cos(2*t));
    s(1,1) = c0*(1-cos(2*t));
    s(0,1) = 
    s(1,0) = c0*sin(2*t);
    return s;
 }
 tensor time_derivative (const point& x) const {
    Float c0 = f(x);
    Float c1 = df_dt(x);
    tensor s;
    s(0,0) = c1*(1+cos(2*t)) - 2*c0*sin(2*t);
    s(1,1) = c1*(1-cos(2*t)) + 2*c0*sin(2*t);
    s(0,1) =
    s(1,0) = c1*sin(2*t)     + 2*c0*cos(2*t);
    return s;
 }
 sigma_exact (Float nu1, Float t1 = 0)
  : nu(nu1), t(t1), r0(0.1), x0(0.25,0), xt(), vt()
 {
    xt = point( x0[0]*cos(t) - x0[1]*sin(t), 
                x0[0]*sin(t) + x0[1]*cos(t));
    vt = point(-x0[0]*sin(t) - x0[1]*cos(t), 
                x0[0]*cos(t) - x0[1]*sin(t));
 }
 Float nu, t, r0;
 point x0, xt, vt;
};
struct chi {
 tensor operator() (const point& x) const {
    return - _s.time_derivative(x);
 }
 chi (Float nu, Float t = 0): _s(nu, t) {}
 sigma_exact _s;
};
typedef sigma_exact sigma_g;
