#ifndef _Solver_h
  #define _Solver_h

  #include "Mesh.h"
  #include "Types.h"


  typedef struct{
    double2 xb;
    bilintransf dertest;  /* derivat. of test function in 3 integ. nodes */
    double detJF[3];  /* value of Jacobian in 3 integration nodes */
    uint ib;  /* ib=0,1,2  ... index of curved edge within the triangle */
    double normal[4][2];   /* outer normal for curved edge, vertices, Gauss */
    double length[4];  /* length of normal = increase of \Gamma */
    double2 nablarho[3];  /* for viscous terms computation */
    }tSolver_Curve_Element;

  typedef struct{
    tMesh_VolumeElement *ve;
    int icve;
    statevec w,wold,H,K1,K2,K3, HV,Dsxnablaw,Dsynablaw;
    lintransf dertest;
    double2 nablarho;
    double oneovervolume;
    double maxovervolume;  //for timestep
    double circumference; /* the lenght of \partial K_i */
    double tauloc;  /* for local time stepping */
    int neigh[3];
    uint neighind[3];
    uint limited;
    int iface;
    int idxface[3];
    }tSolver_Element;

  typedef struct{
    double2 n;  //normal vector
    double volume;
    uint i,j;  //j is defined only if boundary==0
    uint ii,jj;  /* indexes 0,1,2 of edges for triangles i,j */
    int boundary;  //0=inside, -1=periodic, >0 slip/i/o
    }tSolver_Face;

  typedef enum{
    cSolver_BoundaryType_Inlet,cSolver_BoundaryType_Outlet,
      cSolver_BoundaryType_Slip
    }tSolver_BoundaryType;

  typedef struct{
    double4 w;  //inlet, outlet
    tSolver_BoundaryType type;
    BOOL fixed;  //inlet, outlet
    }tSolver_Boundary;

  typedef struct _tSolver_State{
    tMesh_Mesh *m;
    uint nelem,nbound,nbc,nface,local,ncurvedelements;
    uint number_of_limiting,number_of_nolimiting;
    uint DOF;   /* degree of freedon per one element */
    /*   DOF <=3 */
    double CFL, cL, cD, cDp;
    double kappa,kappa2;  //kappa2=kappa*(kappa-1)
    tSolver_Boundary *boundary;
    tSolver_Element *element;
    tSolver_Face *face;
    tSolver_Curve_Element *curvedelement;
    void (*flux)(const struct _tSolver_State *,const double4,const double4,
      const double2,double4);
    void (*flux1)(const struct _tSolver_State *,const double4,const double4,
      const double2,double4);
    void (*integrate)(struct _tSolver_State *);
    double t,nextt,tau,linfty,Machmax, Machmin,limiting, Re, kappaRePr;
    double Minfty,thetainfty;
    }tSolver_State;

  tSolver_Boundary *Solver_BoundaryCreate(tMesh_Mesh *);
  void Solver_BoundarySetInletOutlet(tSolver_Boundary *,uint,double,double,
    double,double,BOOL,double);
  void Solver_BoundarySetSlip(tSolver_Boundary *,uint);

  void Solver_SetInitialCondition(tMesh_Mesh *,double,double,double,
    double,double);

  void Solver_SetInitialCondition_Sod(tMesh_Mesh *,double,double,double,
    double,double,double,double,double,double);

  void Solver_Flux_Batten(const tSolver_State *,const double4,const double4,
    const double2,double4);
  void Solver_Flux_Roe(const tSolver_State *,const double4,const double4,
    const double2,double4);

  void Solver_Integrate_Euler(tSolver_State *);
  void Solver_Integrate_EulerCauchy(tSolver_State *);
  void Solver_Integrate_RK4(tSolver_State *);

  tSolver_State *Solver_Init(tMesh_Mesh *,tSolver_Boundary *,double,
    void (*)(tSolver_State *));
  void Solver_Solve(tSolver_State *);
  void Solver_MeshUpdate(tSolver_State *);

  void  Seek_Index_Neighbours(tSolver_State *);
  void  Solver_Compute_Inviscid_Volumes(tSolver_State *);
  void  Solver_Compute_Viscous_Volumes(tSolver_State *);
  void  Solver_Compute_Viscous_Edges(tSolver_State *);

  void Solver_Limiting(tSolver_State *);
  void Solver_Limiting_Element(tSolver_State *, uint, double );

  void Linear_Mapping_Computing(tSolver_State *);
  int Curved_Element_Detecting(tMesh_Mesh *);

  int ViscousTerm(double4, double4, double2, double4, double4, double, double, 
		  double4, double4);
  double Edge_Integral(double, double, double, double);
  void Compute_DwnablaW(tSolver_State*, int , uint, int , uint , 
			double2,double2, double, double, double4 , double4 );
  void Compute_Dw_nablaPhi(int , statevec, double2[3], double3, double,
			    double , bilintransf,
			    double [3][3][4][4],double [3][3][4][4]);
  void Show_Viscous_Term(tSolver_State*, uint );
  void Compute_Wjump(statevec , statevec ,uint , uint ,
		     double4, double4);
  BOOL Identical_Nodes(double2, double2);
  double Temperature(double4);
  double Pressure(double4, double);
  void Coeficients(tSolver_State *, double2, double2);
  void Sod_Error(tSolver_State *, double, double, double, double, 
                 double, double, double, double );
  void NSexact(double , double , double4 );
  double Compute_Outlet_Pressure(tSolver_State *);
  void Solver_Inviscid_Stabilization(tSolver_State *);
  void calc_ringleb(double ,double2 , double4 );
  double cspeed_ringleb(double ,double ,double );
  double fun_ringleb(double ,double ,double ,double );
  void state_ringleb(double ,double ,double ,double4 );
 void  Ringleb_Error(tSolver_State *);

void PlotIsolines(FILE *,tMesh_Mesh *, tSolver_State *, char *, uint, BOOL, 
		  double, double);
void PlotWall(FILE *,tMesh_Mesh *, tSolver_State *, char *, int );
void PlotCut(FILE *,tMesh_Mesh *, tSolver_State *, char *, double2, double2);
void CalculateQuantity(tMesh_Mesh *, tSolver_State *, char *);
void Relimiting(tMesh_Mesh *);
void Edge2Volume(FILE *,tMesh_Mesh *);
void Solver_Compute_Velocity_Deformation(tSolver_State *);
int Intersection( double2, double2, double2, double2,  double2 );
double Interpolate(double2 ,double2 , double2 , double2, double3);

#endif
