#ifndef FINITE_VOLUMES_EULER_UNKNOWNS_HPP
#define FINITE_VOLUMES_EULER_UNKNOWNS_HPP

template <typename TMeshData>
class FiniteVolumesEulerUnknowns
{
public:
  typedef TMeshData MeshDataType;
  typedef typename MeshDataType::MeshType MeshType;

  static constexpr size_t dimension = MeshType::dimension;
  typedef TinyVector<dimension> Rd;

private:
  const MeshDataType& m_mesh_data;
  const MeshType& m_mesh;

  Kokkos::View<double*> m_rhoj;
  Kokkos::View<Rd*> m_uj;
  Kokkos::View<double*> m_Ej;

  Kokkos::View<double*> m_ej;
  Kokkos::View<double*> m_pj;
  Kokkos::View<double*> m_gammaj;
  Kokkos::View<double*> m_cj;
  Kokkos::View<double*> m_mj;
  Kokkos::View<double*> m_inv_mj;
  Kokkos::View<double*> m_kj;
  Kokkos::View<Rd*> m_Sj;
  Kokkos::View<Rd*> m_uL;
  Kokkos::View<Rd*> m_uR;
  Kokkos::View<double*> m_kL;
  Kokkos::View<double*> m_kR;

public:
  Kokkos::View<double*> rhoj()
  {
    return m_rhoj;
  }

  const Kokkos::View<const double*> rhoj() const
  {
    return m_rhoj;
  }

  Kokkos::View<Rd*> uj()
  {
    return m_uj;
  }

  const Kokkos::View<const Rd*> uj() const
  {
    return m_uj;
  }

  Kokkos::View<double*> Ej()
  {
    return m_Ej;
  }

  const Kokkos::View<const double*> Ej() const
  {
    return m_Ej;
  }

  Kokkos::View<double*> ej()
  {
    return m_ej;
  }

  const Kokkos::View<const double*> ej() const
  {
    return m_ej;
  }

  Kokkos::View<double*> pj()
  {
    return m_pj;
  }

  const Kokkos::View<const double*> pj() const
  {
    return m_pj;
  }

  Kokkos::View<double*> gammaj()
  {
    return m_gammaj;
  }

  const Kokkos::View<const double*> gammaj() const
  {
    return m_gammaj;
  }

  Kokkos::View<double*> cj()
  {
    return m_cj;
  }

  const Kokkos::View<const double*> cj() const
  {
    return m_cj;
  }

  Kokkos::View<double*> mj()
  {
    return m_mj;
  }

  const Kokkos::View<const double*> mj() const
  {
    return m_mj;
  }

  Kokkos::View<double*> invMj()
  {
    return m_inv_mj;
  }

  const Kokkos::View<const double*> invMj() const
  {
    return m_inv_mj;
  }

  Kokkos::View<double*> kj()
  {
    return m_kj;
  }

  const Kokkos::View<const double*> kj() const
  {
    return m_kj;
  }

   Kokkos::View<Rd*> Sj()
  {
   return m_Sj;
  }

  const Kokkos::View<const Rd*> Sj() const
  {
    return m_Sj;
  }

  Kokkos::View<Rd*> uL()
  {
    return m_uL;
  }

  const Kokkos::View<const Rd*> uL() const
  {
    return m_uL;
  }

  Kokkos::View<Rd*> uR()
  {
    return m_uR;
  }

  const Kokkos::View<const Rd*> uR() const
  {
    return m_uR;
  }

  Kokkos::View<double*> kL()
  {
    return m_kL;
  }

  const Kokkos::View<const double*> kL() const
  {
    return m_kL;
  }

  Kokkos::View<double*> kR()
  {
    return m_kR;
  }

  const Kokkos::View<const double*> kR() const
  {
    return m_kR;
  }

  
  // --- Acoustic Solver ---
  
  void initializeSod()
  {
    const Kokkos::View<const Rd*> xj = m_mesh_data.xj();

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
  	if (xj[j][0]<0.5) {
  	  m_rhoj[j]=1;
  	} else {
  	  m_rhoj[j]=0.125;
  	}
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
  	if (xj[j][0]<0.5) {
  	  m_pj[j]=1;
  	} else {
  	  m_pj[j]=0.1;
  	}
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
  	m_uj[j] = zero;
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
  	m_gammaj[j] = 1.4;
      });

    BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj);
    block_eos.updateEandCFromRhoP();

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
  	m_Ej[j] = m_ej[j]+0.5*(m_uj[j],m_uj[j]);
	//m_Ej[j] = 2.;
      });

    const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
  	m_mj[j] = m_rhoj[j] * Vj[j];
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
  	m_inv_mj[j] = 1./m_mj[j];
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
  	m_kj[j] =  0.00001*0.5; 
      });

     // Conditions aux bords de Dirichlet sur u et k
    
    m_uL[0] = zero;
    m_uR[0] = zero;
    m_kL[0] = 0.00001*0.5;
    m_kR[0] = 0.00001*0.5;
  }
  
  
  /* 

  // DIFFUSION PURE

    void initializeSod()
    {
    const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
    double pi = 4.*std::atan(1.);

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
	m_rhoj[j]=1.;
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
	  m_pj[j]=1.;
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
	m_uj[j] = std::sin(pi*xj[j][0]);
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
	m_gammaj[j] = 2.;
      });

    BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj);
    block_eos.updateEandCFromRhoP();

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
	//m_Ej[j] = m_ej[j]+0.5*(m_uj[j],m_uj[j]);
	m_Ej[j] = 2.;
      });

    const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
	m_mj[j] = m_rhoj[j] * Vj[j];
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
	m_inv_mj[j] = 1./m_mj[j];
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
	m_kj[j] = 2.; // k constant
	//m_kj[j] = xj[j][0]; // k non constant, k = x
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
    	m_Sj[j] = std::sin(pi*xj[j][0])*(pi*pi*xj[j][0]-1.) - std::cos(xj[j][0]*pi)*pi;
      });

    // Conditions aux bords de Dirichlet sur u et k
    
    m_uL[0] = zero;
    m_uR[0] = zero;
    m_kL[0] = 2.;
    m_kR[0] = 2.;

  }
*/
  


  FiniteVolumesEulerUnknowns(const MeshDataType& mesh_data)
    : m_mesh_data(mesh_data),
      m_mesh(m_mesh_data.mesh()),
      m_rhoj("rhoj",m_mesh.numberOfCells()),
      m_uj("uj",m_mesh.numberOfCells()),
      m_Ej("Ej",m_mesh.numberOfCells()),
      m_ej("ej",m_mesh.numberOfCells()),
      m_pj("pj",m_mesh.numberOfCells()),
      m_gammaj("gammaj",m_mesh.numberOfCells()),
      m_cj("cj",m_mesh.numberOfCells()),
      m_mj("mj",m_mesh.numberOfCells()),
      m_inv_mj("inv_mj",m_mesh.numberOfCells()),
      m_kj("kj",m_mesh.numberOfCells()),
      m_Sj("Sj",m_mesh.numberOfCells()),
  m_uL("uL", 1),
  m_uR("uR", 1),
  m_kL("kL", 1),
  m_kR("kR", 1)
  
  {
    ;
  }

};

#endif // FINITE_VOLUMES_EULER_UNKNOWNS_HPP