#ifndef eulerPROBLEM_IMPL_H_
#define eulerPROBLEM_IMPL_H_

#include <core/mfilename.h>
#include <matrices/tnlMatrixSetter.h>
#include <solvers/pde/tnlExplicitUpdater.h>
#include <solvers/pde/tnlLinearSystemAssembler.h>
#include <solvers/pde/tnlBackwardTimeDiscretisation.h>

template< typename Mesh,
          typename BoundaryCondition,
          typename RightHandSide,
          typename DifferentialOperator >
tnlString
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
getTypeStatic()
{
   return tnlString( "eulerProblem< " ) + Mesh :: getTypeStatic() + " >";
}

template< typename Mesh,
          typename BoundaryCondition,
          typename RightHandSide,
          typename DifferentialOperator >
tnlString
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
getPrologHeader() const
{
   return tnlString( "euler" );
}

template< typename Mesh,
          typename BoundaryCondition,
          typename RightHandSide,
          typename DifferentialOperator >
void
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
writeProlog( tnlLogger& logger, const tnlParameterContainer& parameters ) const
{
   /****
    * Add data you want to have in the computation report (log) as follows:
    * logger.writeParameter< double >( "Parameter description", parameter );
    */
}

template< typename Mesh,
          typename BoundaryCondition,
          typename RightHandSide,
          typename DifferentialOperator >
bool
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
setup( const tnlParameterContainer& parameters )
{
   if( ! this->boundaryCondition.setup( parameters, "boundary-conditions-" ) ||
       ! this->rightHandSide.setup( parameters, "right-hand-side-" ) )
      return false;
   return true;
}

template< typename Mesh,
          typename BoundaryCondition,
          typename RightHandSide,
          typename DifferentialOperator >
typename eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::IndexType
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
getDofs( const MeshType& mesh ) const
{
   /****
    * Return number of  DOFs (degrees of freedom) i.e. number
    * of unknowns to be resolved by the main solver.
    */
   return 3*mesh.template getEntitiesCount< typename MeshType::Cell >();
}

template< typename Mesh,
          typename BoundaryCondition,
          typename RightHandSide,
          typename DifferentialOperator >
void
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
bindDofs( const MeshType& mesh,
          DofVectorType& dofVector )
{
}

template< typename Mesh,
          typename BoundaryCondition,
          typename RightHandSide,
          typename DifferentialOperator >
bool
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
setInitialCondition( const tnlParameterContainer& parameters,
                     const MeshType& mesh,
                     DofVectorType& dofs,
                     MeshDependentDataType& meshDependentData )
{
   typedef typename MeshType::Cell Cell;
   double gamma = parameters.getParameter< double >( "gamma" );
   double rhoL = parameters.getParameter< double >( "left-density" );
   double velL = parameters.getParameter< double >( "left-velocity" );
   double preL = parameters.getParameter< double >( "left-pressure" );
   double eL = ( preL / (gamma - 1) ) + 0.5 * rhoL * velL * velL;
   double rhoR = parameters.getParameter< double >( "right-density" );
   double velR = parameters.getParameter< double >( "right-velocity" );
   double preR = parameters.getParameter< double >( "right-pressure" );
   double eR = ( preR / (gamma - 1) ) + 0.5 * rhoR * velR * velR;
   double x0 = parameters.getParameter< double >( "riemann-border" );
   cout << gamma << " " << rhoL << " " << velL << " " << preL << " " << eL << " " << rhoR << " " << velR << " " << preR << " " << eR << " " << x0 << " " << gamma << endl;
   int count = mesh.template getEntitiesCount< Cell >()/3;
   this->rho.bind(dofs,0,count);
   this->rhoVel.bind(dofs,count,count);
   this->energy.bind(dofs,2 * count,count);
   this->data.setSize(2*count);
   this->pressure.bind(this->data,0,count);
   this->velocity.bind(this->data,count,count);
   for(long int i = 0; i < count; i++)
      if (i < x0 * count )
         {
            this->rho[i] = rhoL;
            this->rhoVel[i] = rhoL * velL;
            this->energy[i] = eL;
            this->velocity[i] = velL;
            this->pressure[i] = preL;
         }
      else
         {
            this->rho[i] = rhoR;
            this->rhoVel[i] = rhoR * velR;
            this->energy[i] = eR;
            this->velocity[i] = velR;
            this->pressure[i] = preR;
         };
   this->gamma = gamma;
   cout << "dofs = " << dofs << endl;
   getchar();
  
   
   /*
   const tnlString& initialConditionFile = parameters.getParameter< tnlString >( "initial-condition" );
   if( ! dofs.load( initialConditionFile ) )
   {
      cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." << endl;
      return false;
   }
   */
   return true; 
}
template< typename Mesh,
          typename BoundaryCondition,
          typename RightHandSide,
          typename DifferentialOperator >
   template< typename Matrix >
bool
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
setupLinearSystem( const MeshType& mesh,
                   Matrix& matrix )
{
   const IndexType dofs = this->getDofs( mesh );
   typedef typename Matrix::CompressedRowsLengthsVector CompressedRowsLengthsVectorType;
   CompressedRowsLengthsVectorType rowLengths;
   if( ! rowLengths.setSize( dofs ) )
      return false;
   tnlMatrixSetter< MeshType, DifferentialOperator, BoundaryCondition, CompressedRowsLengthsVectorType > matrixSetter;
   matrixSetter.template getCompressedRowsLengths< typename Mesh::Cell >( mesh,
                                                                          differentialOperator,
                                                                          boundaryCondition,
                                                                          rowLengths );
   matrix.setDimensions( dofs, dofs );
   if( ! matrix.setCompressedRowsLengths( rowLengths ) )
      return false;
   return true;
}

template< typename Mesh,
          typename BoundaryCondition,
          typename RightHandSide,
          typename DifferentialOperator >
bool
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
makeSnapshot( const RealType& time,
              const IndexType& step,
              const MeshType& mesh,
              DofVectorType& dofs,
              MeshDependentDataType& meshDependentData )
{
   cout << endl << "Writing output at time " << time << " step " << step << "." << endl;
   this->bindDofs( mesh, dofs );
   tnlString fileName;
   ofstream vysledek;
   cout << "pressure:" << endl;
   for (int i = 0; i<100; i++) cout << this->pressure[i] << " " ;
   vysledek.open("pressure" + to_string(step) + ".txt");
      for (int i = 0; i<101; i++)
      vysledek << 0.01*i << " " << pressure[i] << endl;
   vysledek.close();
   cout << " " << endl;
   cout << "velocity:" << endl;
   for (int i = 0; i<100; i++) cout << this->velocity[i] << " " ;
   vysledek.open("velocity" + to_string(step) + ".txt");
      for (int i = 0; i<101; i++)
      vysledek << 0.01*i << " " << pressure[i] << endl;
   vysledek.close();
   cout << "energy:" << endl;
   for (int i = 0; i<100; i++) cout << this->energy[i] << " " ;
   vysledek.open("energy" + to_string(step) + ".txt");
      for (int i = 0; i<101; i++)
      vysledek << 0.01*i << " " << energy[i] << endl;
   vysledek.close();
   cout << " " << endl;
   cout << "density:" << endl;
   for (int i = 0; i<100; i++) cout << this->rho[i] << " " ;
   vysledek.open("density" + to_string(step) + ".txt");
      for (int i = 0; i<101; i++)
      vysledek << 0.01*i << " " << rho[i] << endl;
   vysledek.close();
   getchar();

   FileNameBaseNumberEnding( "rho-", step, 5, ".tnl", fileName );
   if( ! rho.save( fileName ) )
      return false;
   FileNameBaseNumberEnding( "rhoVel-", step, 5, ".tnl", fileName );
   if( ! rhoVel.save( fileName ) )
      return false;
   FileNameBaseNumberEnding( "energy-", step, 5, ".tnl", fileName );
   if( ! energy.save( fileName ) )
      return false;
   return true;
}

template< typename Mesh,
          typename BoundaryCondition,
          typename RightHandSide,
          typename DifferentialOperator >
void
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
getExplicitRHS( const RealType& time,
                const RealType& tau,
                const MeshType& mesh,
                DofVectorType& _u,
                DofVectorType& _fu,
                MeshDependentDataType& meshDependentData )
{
    typedef typename MeshType::Cell Cell;
    int count = mesh.template getEntitiesCount< Cell >()/3;
	//bind _u
    this->_uRho.bind(_u,0,count);
    this->_uRhoVelocity.bind(_u,count,count);
    this->_uEnergy.bind(_u,2 * count,count);
		
	//bind _fu
    this->_fuRho.bind(_u,0,count);
    this->_fuRhoVelocity.bind(_u,count,count);
    this->_fuEnergy.bind(_u,2 * count,count);

   MeshFunctionType velocity( mesh, this->velocity );
   MeshFunctionType pressure( mesh, this->pressure );
   MeshFunctionType uRho( mesh, _uRho ); 
   MeshFunctionType fuRho( mesh, _fuRho );
   MeshFunctionType uRhoVelocity( mesh, _uRhoVelocity ); 
   MeshFunctionType fuRhoVelocity( mesh, _fuRhoVelocity );
   MeshFunctionType uEnergy( mesh, _uEnergy ); 
   MeshFunctionType fuEnergy( mesh, _fuEnergy );

/*
   //velocity
   EulerVelGetter.setRhoVel(_urhoVel);
   EulerVelGetter.setRho(_urho); 
   tnlOperatorFunction< EulerVelGetter, MeshFunction, void, time > OFVel;
   velocity = OFVel;

   //pressure
   EulerPressureGetter.setRhoVel(_urhoVel);
   EulerPressureGetter.setVel(velocity);
   EulerPressureGetter.setGamma(gamma);
   EulerPressureGetter.setEnergy(energy); 
   tnlOperatorFunction<EulerPressureGetter,MeshFunction,void,time>OFPressure;
   pressure = OFPressure;

   //rho
   this->bindDofs( mesh, _u );
   LaxFridrichsContinuity.setTau(tau);
   LaxFridrichsContinuity.setVelocity(velocity);
   tnlExplicitUpdater< Mesh, MeshFunctionType, DifferentialOperator, BoundaryCondition, RightHandSide > explicitUpdater;
   explicitUpdater.template update< typename Mesh::Cell >( time,
                                                           mesh,
                                                           this->LaxFridrichsContinuity,
                                                           this->boundaryCondition,
                                                           this->rightHandSide,
                                                           uRho,
                                                           fuRho );

   //rhoVelocity
   LaxFridrichsMomentum.setTau(tau);
   LaxFridrichsMomentum.setVelocity(velocity);
   LaxFridrichsMomentum.setPressure(pressure);
   tnlExplicitUpdater< Mesh, MeshFunctionType, DifferentialOperator, BoundaryCondition, RightHandSide > explicitUpdater;
   explicitUpdater.template update< typename Mesh::Cell >( time,
                                                           mesh,
                                                           this->LaxFridrichsMomentum,
                                                           this->boundaryCondition,
                                                           this->rightHandSide,
                                                           uRhoVelocity,
                                                           fuRhoVelocity );
   

   //energy
   LaxFridrichsEnergy.setTau(tau);
   LaxFridrichsEnergy.setPressure(pressure);
   LaxFridrichsEnergy.setVelocity(velocity);
   tnlExplicitUpdater< Mesh, MeshFunctionType, DifferentialOperator, BoundaryCondition, RightHandSide > explicitUpdater;
   explicitUpdater.template update< typename Mesh::Cell >( time,
                                                           mesh,
                                                           this->LaxFridrichsEnergy,
                                                           this->boundaryCondition,
                                                           this->rightHandSide,
                                                           uEnergy,
                                                           fuEnergy );   
   
   tnlBoundaryConditionsSetter< MeshFunctionType, BoundaryCondition > boundaryConditionsSetter; 
   boundaryConditionsSetter.template apply< typename Mesh::Cell >( 
      this->boundaryCondition, 
      time + tau, 
       u );*/
 }

template< typename Mesh,
          typename BoundaryCondition,
          typename RightHandSide,
          typename DifferentialOperator >
   template< typename Matrix >
void
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
assemblyLinearSystem( const RealType& time,
                      const RealType& tau,
                      const MeshType& mesh,
                      DofVectorType& _u,
                      Matrix& matrix,
                      DofVectorType& b,
                      MeshDependentDataType& meshDependentData )
{
   tnlLinearSystemAssembler< Mesh,
                             MeshFunctionType,
                             DifferentialOperator,
                             BoundaryCondition,
                             RightHandSide,
                             tnlBackwardTimeDiscretisation,
                             Matrix,
                             DofVectorType > systemAssembler;

   tnlMeshFunction< Mesh > u( mesh, _u );
   systemAssembler.template assembly< typename Mesh::Cell >( time,
                                                             tau,
                                                             mesh,
                                                             this->differentialOperator,
                                                             this->boundaryCondition,
                                                             this->rightHandSide,
                                                             u,
                                                             matrix,
                                                             b );
}

#endif /* eulerPROBLEM_IMPL_H_ */