/*
What's new:
C44 - [Martin Zoubek] first public release
C49 - [Martin Zoubek] * Solver_Init fixed for periodic boundary conditions,
      works correctly only if periodic sides are not rotated
C51 - [Martin Zoubek] * added rotations of velocity vector, solver should
      now work with any periodic boundaries
*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "Memory.h"
#include "FileIO.h"
#include "Mesh.h"
#include "Progress.h"
#include "Solver.h"
/*#include "Tetrahedron.h" */
#include "Triangle.h"
#include "Types.h"

tSolver_Boundary *Solver_BoundaryCreate(tMesh_Mesh *m){
  uint nbc,i;
  tMesh_SurfaceElement *se;
  tSolver_Boundary *b;

  if(!m)
    return 0;

  nbc=0;  /* nbc is number of surfaces */
  se=m->SurfaceElement;
  while(se){
    if(se->surfnr>nbc)
      nbc=se->surfnr;
    se=se->Next;
    }
  if(!Memory_GetM(b,tSolver_Boundary,nbc))
    return 0;
  for(i=1;i<=nbc;i++)  /* set all surfaces to 'slip' */
    Solver_BoundarySetSlip(b,i);
  return b;
  }

void Solver_BoundarySetInletOutlet(tSolver_Boundary *b,uint surfnr,double rho,
  double v1,double v2,double p,BOOL fixed,double kappa){

  //    printf("\n\n\n\n OPPOSITE INLET / OUTLET !!!!!!!!!!!!!!!!!!!\n\n\n\n");

  //   printf("\n RINGLEB FLOW INLET / OUTLET !!!!!!!!!!!!!!!!!!!\n"); 

  if(surfnr == 1) b[surfnr-1].type=cSolver_BoundaryType_Inlet;
  else if(surfnr == 2) b[surfnr-1].type=cSolver_BoundaryType_Outlet;
  else if(surfnr == 3) b[surfnr-1].type=cSolver_BoundaryType_Outlet;
  else if(surfnr == 4) b[surfnr-1].type=cSolver_BoundaryType_Outlet; 
  else if(surfnr == 5) b[surfnr-1].type=cSolver_BoundaryType_Outlet; 
  else if(surfnr == 6) b[surfnr-1].type=cSolver_BoundaryType_Outlet; 
  else if(surfnr == 7) b[surfnr-1].type=cSolver_BoundaryType_Outlet; 
  else {
    printf("I do not know if it is inlet or outlet !!!!!\n");
    return;
  }
  b[surfnr-1].w[0]=rho;
  b[surfnr-1].w[1]=rho*v1;
  b[surfnr-1].w[2]=rho*v2;
  b[surfnr-1].w[3]=p/(kappa-1)+0.5*rho*(v1*v1+v2*v2);
  b[surfnr-1].fixed=fixed;
  }

void Solver_BoundarySetSlip(tSolver_Boundary *b,uint surfnr){
  b[surfnr-1].type=cSolver_BoundaryType_Slip;
  }

void Solver_SetInitialCondition_Sod(tMesh_Mesh *m,
  double rhoL,double v1L,double v2L, double pL,
  double rhoR,double v1R,double v2R, double pR,
  double kappa){
  tMesh_VolumeElement *ve;
  uint j,i;
  m->kappa=kappa;
  ve=m->VolumeElement;
  i=0;
  while(ve){
    for(j=0;j<3;j++){
      if(ve->P[j]->Point[0] < 0){
	ve->w[j][0]=rhoL;
	ve->w[j][1]=v1L*rhoL;
	ve->w[j][2]=v2L*rhoL;
	ve->w[j][3]=pL/(kappa-1)+0.5*rhoL*(v1L*v1L+v2L*v2L);
      }
      else{
	ve->w[j][0]=rhoR;
	ve->w[j][1]=v1R*rhoR;
	ve->w[j][2]=v2R*rhoR;
	ve->w[j][3]=pR/(kappa-1)+0.5*rhoR*(v1R*v1R+v2R*v2R);
      }
    }
    ve=ve->Next;
    i++;
  }
 printf("End of Solver_SetInitialCondition_Sod\n");
}
void Solver_SetInitialCondition(tMesh_Mesh *m,double rho,double v1,double v2,
       double p,double kappa){
  tMesh_VolumeElement *ve;
  tMesh_Point *p1,*p2;
  double rho_l,v1_l, v2_l, p_l;
  double2 xp;
  double4 qj,wj;
  uint j,k,counter;
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};
 
  rho_l = rho;
  v1_l = v1;
  v2_l = v2;
  p_l = p;

  m->kappa=kappa;
  
  ve=m->VolumeElement;
  while(ve){
    double bar_rho=0, bar_v1=0, bar_v2=0, bar_p=0;
    double epsilon = 5, radius2,pi=3.1415927, f;
    for(j=0;j<3;j++){

      /* >Vortex moving */
      p1=ve->P[perm[j][0]];
      p2=ve->P[perm[j][1]];
      // computation of integration nodes 
      double2_weight_awer(xp,p1->Point,p2->Point,0.5,0.5);
      xp[0] = xp[0] -5;
      xp[1] = xp[1] -5;
      radius2 = double2_norm2(xp);
      bar_v1 = epsilon/2/pi*exp(0.5*(1-radius2))*(-xp[1]);
      bar_v2 = epsilon/2/pi*exp(0.5*(1-radius2))*xp[0];
      f = -(kappa-1)*epsilon*epsilon/8/kappa/pi/pi*exp(1-radius2);
      bar_rho = f/(kappa-1);
      bar_p = bar_rho*kappa;

      rho_l = rho + bar_rho;
      v1_l = v1 + bar_v1;
      v2_l = v2 + bar_v2;;
      p_l = p + bar_p;

      ve->w[j][0]=rho_l;
      ve->w[j][1]=v1_l*rho_l;
      ve->w[j][2]=v2_l*rho_l;
      ve->w[j][3]=p_l/(kappa-1)+0.5*rho_l*(v1_l*v1_l+v2_l*v2_l);
    }
    ve=ve->Next;
  }
  
/*  RINGLEB FLOW */
/* ve=m->VolumeElement;
  while(ve){
    for(j=0;j<3;j++){
      // points forming the boundary edge 
      p1=ve->P[perm[j][0]];
      p2=ve->P[perm[j][1]];
      // computation of integration nodes 
      double2_weight_awer(xp,p1->Point,p2->Point,0.5,0.5);
      // computation of the solution in integration nodes
      calc_ringleb(kappa,xp,qj);
      // transformation from physical quantities to the state vector 
      vector_q2w(qj,wj,kappa);
      for(k=0;k<4;k++) ve->w[j][k]=wj[k];
      if(counter == 17) 
	printf("%7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %d %d \n",
	       xp[0],xp[1],
	       ve->w[j][0],ve->w[j][1],ve->w[j][2],ve->w[j][3],counter,j);

    }
    counter++;
    ve=ve->Next;
  }
*/
}

  tMesh_VolumeElement **vSolver_Init_QuicksortBase;

static int Solver_Init_QuicksortCompareFunction(const void *a,const void *b){
  Progress_Update((tMesh_VolumeElement **)a-vSolver_Init_QuicksortBase);
  return (pointer)*(tMesh_VolumeElement **)a-
    (pointer)*(tMesh_VolumeElement **)b;
  }

int Curved_Element_Detecting(tMesh_Mesh *m){
  FILE *f0;
  int i,j,k,nbpc,ipoin;
  int ibpoin[10]; /*maximal number of curved part of boundary */
  uint closed[10];
  const int maxlen=72000; /* maximal number of nodes in PROFILES  */ 
  double xb[maxlen][2]; /* coordinates of curved boundary */
  static double zero=1E-2;  /*precision for position */
  uint counter;
  tMesh_SurfaceElement *se;

  #define Exit(s,s1){\
    printf(s,s1);\
    return 0;\
    }

  if((f0=fopen("profiles","rb"))==0)
    Exit("Error: Can't open file'%s'\n","profiles");
  nbpc=FileIO_ReadInteger(f0); 

  ibpoin[0] = 0;
  for(i=0;i<nbpc;i++){
    double dist;
    ipoin=FileIO_ReadInteger(f0); 
    ibpoin[i+1] = ibpoin[i] + ipoin;
    for(j=0;j<ipoin;j++){
      if(ibpoin[i]+j -1>maxlen) 
	Exit("Error: The declared array is too small in comparison with \n",
	     "profiles");
      xb[ibpoin[i]+j][0] = FileIO_ReadFloat(f0);
      xb[ibpoin[i]+j][1] = FileIO_ReadFloat(f0); 
    }

    /* checking if the loop of node on profiles is closed */
    if(Identical_Nodes(xb[ibpoin[i]],xb[ibpoin[i+1]-1])  )
      closed[i] = 1;
    else
      closed[i] = 0;
  }
  close(f0);

  se=m->SurfaceElement;
  m->ncurvedelements = 0;
 
  Progress_Init("Looking up curved elements",m->nbound); 
  counter=0;
  while(se){
    double mindist1 = 10.,mindist2 = 10., dist;
    int min_i1,min_i2,min_i,ipart1,ipart2;
    se->curveelement=0;
    if(!( se->surfnr == 1 ||  se->surfnr == 2)){  /* NOT inlet or outlet */
       /* break; */
      /*
       if(se->surfnr == 3) 
	printf("%7.5g %7.5g boundary \n%7.5g %7.5g\n%\n%",
	       se->P[0]->Point[0],se->P[0]->Point[1],
	       se->P[1]->Point[0],se->P[1]->Point[1]); 
      */
      for(i=0;i<nbpc;i++)   /* seeking the same nodes in xb[] field */
	for(k=ibpoin[i];k<ibpoin[i+1];k++){
	  dist = sqrt(double2_distance2(se->P[0]->Point,xb[k]));
	  if(mindist1 >dist){
	    mindist1 = dist;
	    min_i1 = k;
	    ipart1 = i;
	  }
	  dist = sqrt(double2_distance2(se->P[1]->Point,xb[k]));
	  if(mindist2 >dist){
	    mindist2 = dist;
	    min_i2 = k;
	    ipart2 = i;
	  }

	}
      /*
      printf("%7.5g %7.5g ||(%d, %d) %7.5g %d %7.5g %d Curved \n",
	     se->P[0]->Point[0],se->P[0]->Point[1],ibpoin[i],ibpoin[i+1],
	     mindist1, min_i1,mindist2, min_i2  ); 
      */

      if(mindist1 < zero && mindist2 < zero){
	if(ipart1!=ipart2)
	  printf("DEFERENT PART OF CURVED BOUNDARY DETECTED!!! %d %d \n",
		 ipart1, ipart2);
	
	if( closed[ipart1]){/* checking if the final point was found corectly*/
	/*                    some nodes are two times */
	  if(min_i1 == ibpoin[ipart1] && 
	     inabs(ibpoin[ipart1+1]-1 -min_i2) < inabs(min_i1-min_i2) )
	    min_i1 = ibpoin[ipart1+1]-1;
	  if(min_i1 == ibpoin[i+1]-1 && 
	     inabs(ibpoin[ipart1]-min_i2) < inabs(min_i1-min_i2) )
	    min_i1 = ibpoin[ipart1];
	  if(min_i2 == ibpoin[ipart1] && 
	     inabs(ibpoin[ipart1+1]-1 -min_i1) < inabs(min_i1-min_i2) )
	    min_i2 = ibpoin[ipart1+1]-1;
	  if(min_i2 == ibpoin[ipart1+1]-1 && 
	     inabs(ibpoin[ipart1]-min_i1) < inabs(min_i1-min_i2) )
	    min_i2 = ibpoin[ipart1];
	}

	if(inabs(min_i1-min_i2)<2){
	  printf("Any other node is not presented\n");
	  printf("!!!!!!!!  %d %d %d %d  \n",                           
		min_i1,min_i2,ibpoin[ipart1+1],ibpoin[ipart1]);
	}
	else if( inabs(min_i1-min_i2) > (ibpoin[ipart1+1]-ibpoin[ipart1])/2){
	  /*	  printf("!!!!!!!!  %d %d %d %d  \n",
		 min_i1,min_i2,ibpoin[ipart1+1],ibpoin[ipart1]);
	  printf("!!!!!!!!  %d  \n",
	  (min_i1 + min_i2 - (ibpoin[ipart1+1]-ibpoin[ipart1]) )/2); */
	  se->curveelement=1;
	  min_i = (min_i1 + min_i2 - (ibpoin[ipart1+1]-ibpoin[ipart1] -1) )/2;
	  if(min_i < 0) min_i = ibpoin[ipart1+1]-ibpoin[ipart1]-1 - min_i;
	  for(j=0;j<2;j++)
	    se->xb[j] = xb[min_i][j];

	  /*	  
	  printf("%7.5g %7.5g cur...\n",se->P[0]->Point[0],se->P[0]->Point[1]);
	  printf("%7.5g %7.5g \n",se->xb[0],se->xb[1]);
	  printf("%7.5g %7.5g \n\n",se->P[1]->Point[0],se->P[1]->Point[1]);
	  */
	  m->ncurvedelements++;
	}
	else{
	  se->curveelement=1;
	  min_i = (min_i1 + min_i2)/2;
	  for(j=0;j<2;j++)
	    se->xb[j] = xb[min_i][j];
	  /*
	  printf("%7.5g %7.5g curved\n",se->P[0]->Point[0],se->P[0]->Point[1]);
	  printf("%7.5g %7.5g \n",se->xb[0],se->xb[1]);
	  printf("%7.5g %7.5g \n\n",se->P[1]->Point[0],se->P[1]->Point[1]);
	  */
	/*	if(m->ncurvedelements<3 || m->ncurvedelements > 143)
	  printf("\n %7.5g  %7.5g  %7.5g  %7.5g %7.5g  %d  Distance %d  \n",
	 se->xb[0],se->xb[1],xb[min_i][0],xb[min_i][1],mindist1,min_i,
         m->ncurvedelements); */
	  m->ncurvedelements++;
	
	}
      }
    }
    Progress_Update(++counter); 
    se=se->Next;
  }
  Progress_Done();
  printf("The %d nodes detected on curved part of boundary\n",
	 m->ncurvedelements);

  return 1;
}

tSolver_State *Solver_Init(tMesh_Mesh *m,tSolver_Boundary *b,double CFL,
  void (*integrateproc)(tSolver_State *)){
  tSolver_State *s;
  tMesh_VolumeElement *ve,**element;
  tMesh_SurfaceElement *se;
  tSolver_Curve_Element *cve;
  uint counter,i,j,k,l,ic,jpr;
  int itest = -118;
  #define FindElement(ve,i){\
    uint _a,_b;\
    _a=0;\
    _b=s->nelem-1;\
    while(element[i=(_a+_b)>>1]!=(ve))\
      if((pointer)element[i]<(pointer)(ve))\
        _a=i+1;\
      else\
        _b=i-1;\
    }
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};
  if(!m||!b||!Memory_GetM(s,tSolver_State,1))
    return 0;

  Mesh_FindNeighbours(m);

  s->m=m;
  s->DOF = 3;  /* degree of freedom, piecewise linear */
  s->boundary=b;
  s->CFL=CFL;
  s->kappa=m->kappa;
  s->kappa2=m->kappa*(m->kappa-1);
   s->flux=Solver_Flux_Batten;
  //s->flux=Solver_Flux_Roe;
  s->integrate=integrateproc;
  s->t=0;
  s->nextt=1e20;  /* infinity */
  s->linfty=1e20;  /* infinity */

  s->cL =0.;
  s->cD =0.;
  s->cDp = 0.;

  s->nelem=m->nelem;
  s->nbound=s->nbc=0;
  s->ncurvedelements=m->ncurvedelements;
  se=m->SurfaceElement;
  while(se){
    s->nbound++;
    if(se->surfnr>s->nbc)
      s->nbc=se->surfnr;
    se=se->Next;
    }
  s->nface=(s->nelem<<1)+(s->nbound>>1);  /* number of faces */
  if(!Memory_GetM(s->element,tSolver_Element,s->nelem)||
     !Memory_GetM(s->face,tSolver_Face,s->nface)||
     !Memory_GetM(element,tMesh_VolumeElement *,s->nelem) )
    return 0;
  if(s->ncurvedelements)
    if(!Memory_GetM(s->curvedelement,tSolver_Curve_Element,s->ncurvedelements))
      return 0;
     
  ve=m->VolumeElement;  /* assign elements */
  i=0;
  while(ve){
    element[i++]=ve;
    ve=ve->Next;
    }

  Progress_Init("Sorting elements",s->nelem);
  vSolver_Init_QuicksortBase=element;
  qsort((void *)element,s->nelem,sizeof(tMesh_VolumeElement *),
    Solver_Init_QuicksortCompareFunction);
  Progress_Done();


  //printf("Before !!!!! LOOKING P BC,  %d  %d %d\n",m->nper,s->nelem, s->DOF);

  /*  Looking up periodic boundary faces */
  for(jpr=0;jpr<m->nper;jpr++){
    double2 xper;
    xper[0] = m->rper[2*jpr];
    xper[1] = m->rper[2*jpr+1];
    for(i=0;i<s->nelem;i++){
      for(k=0;k<s->DOF;k++){
	if(!element[i]->Neighbour[k]){ /* no neighbour, candidate for per BC*/
	  double2 x1,x2,xp1,xp2,xd1,xd2;
	  int ii,kk;
	  for(l=0;l<2; l++){
	    x1[l]= element[i]->P[perm[k][0]]->Point[l];
	    x2[l]= element[i]->P[perm[k][1]]->Point[l];
	  }
	  for(ii=0;ii<s->nelem;ii++){
	    for(kk=0;kk<s->DOF;kk++){
	      if(!element[ii]->Neighbour[kk]){ /* candidate for per BC*/
		for(l=0;l<2; l++){
		  xp1[l]= element[ii]->P[perm[kk][0]]->Point[l];
		  xp2[l]= element[ii]->P[perm[kk][1]]->Point[l];
		}
		double2_sub(xd1,x1,xp2);
		double2_sub(xd2,x2,xp1);
		if(double2_distance2(xd1,xper) < 1E-05 && 
		   double2_distance2(xd2,xper) < 1E-05){ 
		  /*  periodic pair found */
		  element[ii]->Neighbour[kk] = element[i];
		  element[i]->Neighbour[k] = element[ii];

		  /*printf("%7.5g %7.5g \n",
			 (x1[0]+x2[0])/2,(x1[1]+x2[1])/2);
		  printf("%7.5g %7.5g \n\n",
			 (xp1[0]+xp2[0])/2,(xp1[1]+xp2[1])/2);
		  */
		  
		}
		// second possibility
		double2_sub(xd1,xp2,x1);
		double2_sub(xd2,xp1,x2);
		if(double2_distance2(xd1,xper) < 1E-05 && 
		   double2_distance2(xd2,xper) < 1E-05){ 
		  /*  periodic pair found */
		  element[ii]->Neighbour[kk] = element[i];
		  element[i]->Neighbour[k] = element[ii];

		  /*printf("%7.5g %7.5g \n",
			 (x1[0]+x2[0])/2,(x1[1]+x2[1])/2);
		  printf("%7.5g %7.5g \n\n",
			 (xp1[0]+xp2[0])/2,(xp1[1]+xp2[1])/2);
		  */
		  
		}
	      }
	    }
	  }
	}  /* end of if */
      }  /* end of k */
    } /* end of i */
  }  /* end of jper */

  //  printf("Before FACEs\n"); 


  Progress_Init("Looking up faces",s->nelem); 
  counter=0;
  ic=0;  /* index of a curved element */
  for(i=j=0;i<s->nelem;i++){
    for(k=0;k<s->DOF;k++){
      double x2,p,rho,v1,theta,theta0,e;
      s->element[i].ve=element[i];
      s->element[i].w[k][0]=element[i]->w[k][0];
      s->element[i].w[k][1]=element[i]->w[k][1];
      s->element[i].w[k][2]=element[i]->w[k][2];
      s->element[i].w[k][3]=element[i]->w[k][3];

      /* NSexact solution for NS */
      /*      x2=(s->element[i].ve->P[perm[k][0]]->Point[1]+
          s->element[i].ve->P[perm[k][1]]->Point[1])/2;
      p=1.6;
      theta0=4.257142857;
      v1=x2;
      theta=theta0-0.257142857*x2*x2;\\ ^^^^^^^^^^^ = 0.5*Pr/kappa  
      rho=p/theta/(s->kappa-1);
      e=p/(s->kappa-1)+0.5*rho*v1*v1;
      s->element[i].w[k][0]=rho;
      s->element[i].w[k][1]=rho*v1;
      s->element[i].w[k][2]=0;
      s->element[i].w[k][3]=e;
      //if(i==itest)
        printf("\n IC %d %17.15g  %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g\n",
        k,rho,v1,e,p,rho*v1,theta,s->element[i].w[k][1]);
      */
      /*
      printf("IC %d %d %7.5g %7.5g %7.5g %7.5g \n",
	     i,k,s->element[i].w[k][0],s->element[i].w[k][1],
	     s->element[i].w[k][2],s->element[i].w[k][3]);
      */
    }


    s->element[i].oneovervolume=1/Triangle_Volume(element[i]->P[0]->Point,
      element[i]->P[1]->Point,element[i]->P[2]->Point);
    s->element[i].maxovervolume=0;
    s->element[i].circumference=0;
    s->element[i].icve=-1; /*non curved element */
    for(k=0;k<3;k++){
      double volume;
      volume=Edge_Length(element[i]->P[perm[k][0]]->Point,
        element[i]->P[perm[k][1]]->Point);
      s->element[i].circumference+=volume;

      if(s->element[i].maxovervolume<volume)
        s->element[i].maxovervolume=volume;
      if((pointer)element[i]->Neighbour[k]<(pointer)element[i]){
        Triangle_NormalVector(element[i]->P[perm[k][0]]->Point,
          element[i]->P[perm[k][1]]->Point,element[i]->P[perm[k][2]]->Point,
          s->face[j].n);
        s->face[j].volume=volume;
        s->face[j].i=i;
        s->face[j].ii=k;
	//s->element[i].idxface[k] = j;
	

        if(element[i]->Neighbour[k]){
          FindElement(element[i]->Neighbour[k],s->face[j].j);
	  if(!element[i]->Surface[k]){  /* inner edge */
	    for(l=0;l<3;l++){
	      if(s->element[i].ve->P[k] == 
		 s->element[s->face[j].j].ve->P[l])s->face[j].jj=perm[l][2];
	      s->face[j].boundary =0;
	      /*	      if(i == 544) */
	      /*
	      printf("neigh-inn %d: %d %d %d  %d \n",
		     j,s->face[j].j,s->face[j].jj,s->face[j].i,s->face[j].ii);
	      */
	    }
	  }
	  else{ /* periodic boundary edge */
	    /*	    for(l=0;l<3;l++){
	      if(s->element[i].ve->P[k]->PeriodicCopy == 
		 s->element[s->face[j].j].ve->P[l])s->face[j].jj=perm[l][2];
	    }
	    */
	    for(jpr=0;jpr<m->nper;jpr++){
	      int ii,kk;
	      double2 xper;
	      double2 x1,x2,xp1,xp2,xd1,xd2;
	      xper[0] = m->rper[2*jpr];
	      xper[1] = m->rper[2*jpr+1];
	      for(l=0;l<2; l++){
		x1[l]= element[i]->P[perm[k][0]]->Point[l];
		x2[l]= element[i]->P[perm[k][1]]->Point[l];
	      }
	      ii = s->face[j].j;
	      for(kk=0; kk<s->DOF; kk++){
		for(l=0;l<2; l++){
		  xp1[l]= element[ii]->P[perm[kk][0]]->Point[l];
		  xp2[l]= element[ii]->P[perm[kk][1]]->Point[l];
		}
		double2_sub(xd1,x1,xp2);
		double2_sub(xd2,x2,xp1);
		if(double2_distance2(xd1,xper) < 1E-05 && 
		   double2_distance2(xd2,xper) < 1E-05){ 
		  /*  periodic pair found */
		  s->face[j].jj=kk;
		}	
		double2_sub(xd1,xp2,x1);
		double2_sub(xd2,xp1,x2);
		if(double2_distance2(xd1,xper) < 1E-05 && 
		   double2_distance2(xd2,xper) < 1E-05){ 
		  /*  periodic pair found */
		  s->face[j].jj=kk;
		}	
	      }
	    }
	    s->face[j].boundary =-1;
	    
	    /*	      if(i == 544) */
	    
	    /*
	      printf("neigh-per %d: %d %d %d  %d \n",
		     j,s->face[j].j,s->face[j].jj,s->face[j].i,s->face[j].ii);
	    */
	  }
	  /* s->face[j].boundary=element[i]->Surface[k]?-1:0;*/
	  /*periodic-inside*/
	  /*	  if(s->face[j].boundary == -1){
		  printf("######### %d %d %d %d \n",
		  i,k,s->face[j].j,s->face[j].jj); 
		  }*/
	}
        else{
          se=element[i]->Surface[k];
	  switch(b[se->surfnr-1].type){
            case cSolver_BoundaryType_Inlet:
              s->face[j].boundary=se->surfnr;
              break;
            case cSolver_BoundaryType_Outlet:
              s->face[j].boundary=se->surfnr;
              break;
            case cSolver_BoundaryType_Slip:
              s->face[j].boundary=se->surfnr;
              break;
	     
	  }
	  if(se->curveelement == 1){
	    s->curvedelement[ic].xb[0] = se->xb[0];
	    s->curvedelement[ic].xb[1] = se->xb[1];
	    s->curvedelement[ic].ib = k;
	    s->element[i].icve=ic;
	    ic++;
	  }
	}

	/*
	printf("Face: %d %d %d %d %d %d \n",
	       j,s->face[j].boundary,s->face[j].i,s->face[j].j,
	       s->face[j].ii,s->face[j].jj);  
	*/
        j++;
      }
    }
    s->element[i].maxovervolume*=s->element[i].oneovervolume; 
    /* maximal length of edge over area of triangle */
    //    Progress_Update(++counter);
  }
  
  /*  for(i=0;i<s->nbc;i++)
      printf( "21 Boundary %d %d  \n",i,b[i].type); */
     

  /*  for(i=0;i<s->nelem;i++){
    if( s->element[i].icve>=0){
      ic=s->element[i].icve;
      printf( "%7.5g %7.5g Curved edge %d  %d  \n",
	     s->curvedelement[ic].xb[0], s->curvedelement[ic].xb[1],i,ic);
    }
    }  */

  /*    for(i=0;i<s->nelem;i++){
    printf(" %10.5g %10.5g \n" ,s->element[i].F[2][0],s->element[i].F[2][1]);
    printf(" %10.5g %10.5g \n" ,s->element[i].F[0][0]+s->element[i].F[2][0],
	   s->element[i].F[0][1]+s->element[i].F[2][1]);
    printf(" %10.5g %10.5g \n" ,s->element[i].F[1][0]+s->element[i].F[2][0],
	   s->element[i].F[1][1]+s->element[i].F[2][1]);
    printf("\n");
    } */

  s->nface=j;  /* for periodic boundary conditions, do not remove! */
  Progress_Done();
  Seek_Index_Neighbours(s);
  Linear_Mapping_Computing(s);

  Memory_Free(element);
  return s;
}
void  Seek_Index_Neighbours(tSolver_State *s){
  int i;
  
  for(i=0;i<s->nface;i++){
    uint _i,_j,_ii,_jj;
    _i=s->face[i].i;
    _ii= s->face[i].ii;
    if(s->face[i].boundary<=0){  /* inside or periodic */
      _j=s->face[i].j;
      _jj=s->face[i].jj;
      s->element[_i].neigh[_ii]=_j;
      s->element[_i].neighind[_ii]=_jj;
      s->element[_j].neigh[_jj]=_i;
      s->element[_j].neighind[_jj]=_ii;
    }
    else{
      s->element[_i].iface=i;
      if(s->boundary[s->face[i].boundary-1].type==
	    cSolver_BoundaryType_Inlet)
	s->element[_i].neigh[_ii]=-1;
      else if(s->boundary[s->face[i].boundary-1].type==
	      cSolver_BoundaryType_Outlet)
	s->element[_i].neigh[_ii]=-2;
      else if(s->boundary[s->face[i].boundary-1].type==
	      cSolver_BoundaryType_Slip)
	s->element[_i].neigh[_ii]=-3;
      else
	printf("The boundary edge DOES NOT DETECTED !!!!!!!\n");
    }
  }
  /*
  for(i=0;i<s->nelem;i++){
    printf("Neighbours %5d: %5d(%d)  %5d(%d)  %5d(%d)\n",i,
	   s->element[i].neigh[0], s->element[i].neighind[0],
	   s->element[i].neigh[1], s->element[i].neighind[1],
	   s->element[i].neigh[2], s->element[i].neighind[2]);
  } 
  */ 
}

void Linear_Mapping_Computing(tSolver_State *s){
  tMesh_VolumeElement *ve;
  uint i,j,k,l,m;
  int ic,ib;
  static int itest = -52;
  double a[2];
  static double xhat[3][2]={{0.5,0.5},{0,0.5},{0.5,0}};
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};

  a[0] = (1/sqrt(3) + 1)/2;
  a[1] = (-1/sqrt(3) + 1)/2;

  printf(" Begin of Linear_Mapping_Computing\n");
  for(i=0;i<s->nelem;i++){  /* in fact dertest times area(K_i)  */
    ve= s->element[i].ve;
    s->element[i].dertest[0][0] =  ve->P[1]->Point[1] - ve->P[0]->Point[1];
    s->element[i].dertest[0][1] =-(ve->P[1]->Point[0] - ve->P[0]->Point[0]);
    s->element[i].dertest[1][0] =  ve->P[2]->Point[1] - ve->P[1]->Point[1];
    s->element[i].dertest[1][1] =-(ve->P[2]->Point[0] - ve->P[1]->Point[0]);
    s->element[i].dertest[2][0] =  ve->P[0]->Point[1] - ve->P[2]->Point[1];
    s->element[i].dertest[2][1] =-(ve->P[0]->Point[0] - ve->P[2]->Point[0]);

    
    /* for curved elements, curved edge is P[1] P[2]  */
    if( s->element[i].icve>=0){
      double2 b,d,e;
      ic=s->element[i].icve;
      ib=s->curvedelement[ic].ib;
      
     for(j=0;j<2;j++){   /* b[j] is in fact b[j]/2 !!!!! */
       b[j]= 2*( s->curvedelement[ic].xb[j] -0.5* 
	  (ve->P[perm[ib][0]]->Point[j]+ve->P[perm[ib][1]]->Point[j]) );
       d[j]= ve->P[perm[ib][0]]->Point[j] - ve->P[perm[ib][2]]->Point[j];
       e[j]= ve->P[perm[ib][1]]->Point[j] - ve->P[perm[ib][2]]->Point[j];
    }
     /*      b[0] = 0;
	b[1] = 0; */

    s->curvedelement[ic].dertest[perm[ib][0]][0][perm[ib][0]] = e[1]-d[1];
    s->curvedelement[ic].dertest[perm[ib][0]][0][perm[ib][1]] = e[1]-d[1]-b[1];
    s->curvedelement[ic].dertest[perm[ib][0]][0][perm[ib][2]] = e[1]-d[1]+b[1];
    s->curvedelement[ic].dertest[perm[ib][0]][1][perm[ib][0]] = d[0]-e[0];
    s->curvedelement[ic].dertest[perm[ib][0]][1][perm[ib][1]] = d[0]-e[0]+b[0];
    s->curvedelement[ic].dertest[perm[ib][0]][1][perm[ib][2]] = d[0]-e[0]-b[0];

    s->curvedelement[ic].dertest[perm[ib][1]][0][perm[ib][0]] = -b[1]-e[1];
    s->curvedelement[ic].dertest[perm[ib][1]][0][perm[ib][1]] = -e[1];
    s->curvedelement[ic].dertest[perm[ib][1]][0][perm[ib][2]] = -b[1]-e[1];
    s->curvedelement[ic].dertest[perm[ib][1]][1][perm[ib][0]] = b[0]+e[0];
    s->curvedelement[ic].dertest[perm[ib][1]][1][perm[ib][1]] = e[0];
    s->curvedelement[ic].dertest[perm[ib][1]][1][perm[ib][2]] = b[0]+e[0];

    s->curvedelement[ic].dertest[perm[ib][2]][0][perm[ib][0]] = b[1]+d[1];
    s->curvedelement[ic].dertest[perm[ib][2]][0][perm[ib][1]] = b[1]+d[1];
    s->curvedelement[ic].dertest[perm[ib][2]][0][perm[ib][2]] = d[1];
    s->curvedelement[ic].dertest[perm[ib][2]][1][perm[ib][0]] = -b[0]-d[0];
    s->curvedelement[ic].dertest[perm[ib][2]][1][perm[ib][1]] = -b[0]-d[0];
    s->curvedelement[ic].dertest[perm[ib][2]][1][perm[ib][2]] = -d[0];
      
    /*    if(i == itest) { 
      printf("%7.5g  %7.5g  Element\n",ve->P[0]->Point[0],ve->P[0]->Point[1]);
      printf("%7.5g  %7.5g  Element\n",ve->P[1]->Point[0],ve->P[1]->Point[1]);
      printf("%7.5g  %7.5g  Element\n",ve->P[2]->Point[0],ve->P[2]->Point[1]);
      printf("%7.5g  %7.5g  Node\n",
             s->curvedelement[ic].xb[0],s->curvedelement[ic].xb[1]);
      printf("b[0] = %7.5g, b[1]= %7.5g\n",b[0],b[1]);
      printf("Node %4d %4d %4d",i,ic,ib);
      for(k=0;k<3;k++)
	for(l=0;l<2;l++)
	  printf("%9.4g", s->element[i].dertest[k][l]);
      printf("\n");
      for(m=0;m<3;m++){
	printf("Curv %4d %4d %4d",i,ib,m);
	for(k=0;k<3;k++)
	  for(l=0;l<2;l++)
	    printf("%9.4g", s->curvedelement[ic].dertest[k][l][m] );
	printf("\n");
      }
      
      } 
      */  
 /*    printf("---------------------------------------------------\n"); */
    
    for(j=0;j<3;j++)
      s->curvedelement[ic].detJF[perm[ib][j]]=
	2*xhat[j][0]*(b[1]*d[0]-b[0]*d[1])+
	2*xhat[j][1]*(b[0]*e[1]-b[1]*e[0]) + d[0]*e[1] - d[1]*e[0];

        /* for(j=0;j<3;j++)
	   printf("Weights: %d  %d %d  %7.5g %7.5g  %7.5g \n",i,ic,ib,
	   s->curvedelement[ic].detJF[0],
	   s->curvedelement[ic].detJF[1],s->curvedelement[ic].detJF[2]); 
	   printf("b,c,d %7.5g %7.5g  %7.5g \n",2*b[0],d[0],e[0]);
	   printf("b,c,d %7.5g %7.5g  %7.5g \n",2*b[1],d[1],e[1]);
	   printf("-----------------------");   */
    
    s->curvedelement[ic].normal[0][0] =  -2*b[1]*(a[1]-a[0])-d[1]+e[1];
    s->curvedelement[ic].normal[0][1] =   2*b[0]*(a[1]-a[0])+d[0]-e[0];
    s->curvedelement[ic].normal[1][0] =  -2*b[1]*(a[0]-a[1])-d[1]+e[1];
    s->curvedelement[ic].normal[1][1] =   2*b[0]*(a[0]-a[1])+d[0]-e[0];
    s->curvedelement[ic].normal[2][0] =   2*b[1]-d[1]+e[1];
    s->curvedelement[ic].normal[2][1] =  -2*b[0]+d[0]-e[0];
    s->curvedelement[ic].normal[3][0] =  -2*b[1]-d[1]+e[1];
    s->curvedelement[ic].normal[3][1] =   2*b[0]+d[0]-e[0];

    for(k=0;k<4;k++){
	s->curvedelement[ic].length[k]=
	  sqrt(double2_norm2(s->curvedelement[ic].normal[k]));
	for(l=0;l<2;l++)
	  s->curvedelement[ic].normal[k][l]*=1/s->curvedelement[ic].length[k];
    }
      if(ic == -31){ 
          printf("Normals: %d %d  %8.4g %8.4g %8.4g %8.4g \n",
                 i,ic,
		 s->curvedelement[ic].normal[0][0],
		 s->curvedelement[ic].normal[1][0],
		 s->curvedelement[ic].normal[2][0],
		 s->curvedelement[ic].normal[3][0]); 
          printf("Normals: %d %d  %8.4g %8.4g %8.4g %8.4g\n",
                 i,ic,
		 s->curvedelement[ic].normal[0][1],
		 s->curvedelement[ic].normal[1][1],
		 s->curvedelement[ic].normal[2][1],
		 s->curvedelement[ic].normal[3][1]); 


	  } 
    }


  }
  printf(" END of Linear_Mapping_Computing\n");
}

void Solver_Flux_Roe(const tSolver_State *s,const double4 wi,
		     const double4 wj,const double2 n,double4 H){
  
  double wi_rot[4],wj_rot[4],qL[4],qR[4],
    f_roe[4],rlam[4],r[4][4],dw[5],gamma[4],f1[4],
    qqL,qqR,qqq,rhat,uhat,vhat,ehat,vhat2,ahat,pres;
  uint j,k;
  /* rotation */
  wi_rot[0] = wi[0];
  wi_rot[1] = n[0]*wi[1] + n[1]*wi[2];
  wi_rot[2] = -n[1]*wi[1] + n[0]*wi[2];
  wi_rot[3] = wi[3];

  wj_rot[0] = wj[0];
  wj_rot[1] = n[0]*wj[1] + n[1]*wj[2];
  wj_rot[2] = -n[1]*wj[1] + n[0]*wj[2];
  wj_rot[3] = wj[3];
  /*   physical variables rhu,u,v,H (enthalpy)  */
  qL[0] = wi_rot[0];
  qL[1] = wi_rot[1]/wi_rot[0];
  qL[2] = wi_rot[2]/wi_rot[0];
  qL[3] = (wi_rot[3]*s->kappa - 0.5*(s->kappa-1)*
	   (wi_rot[1]*wi_rot[1]+wi_rot[2]*wi_rot[2])/wi_rot[0])/wi_rot[0];

  qR[0] = wj_rot[0];
  qR[1] = wj_rot[1]/wj_rot[0];
  qR[2] = wj_rot[2]/wj_rot[0];
  qR[3] = (wj_rot[3]*s->kappa - 0.5*(s->kappa-1)*
	   (wj_rot[1]*wj_rot[1]+wj_rot[2]*wj_rot[2])/wj_rot[0])/wj_rot[0];

  /*  printf("qL %7.5g %7.5g %7.5g %7.5g \n",
	      qL[0],qL[1],qL[2],qL[3]);

  printf("qR %7.5g %7.5g %7.5g %7.5g \n",
	      qR[0],qR[1],qR[2],qR[3]);
  */
  qqL = pow_fast(qL[0],0.5);
  qqR = pow_fast(qR[0],0.5);
  qqq = qqL+qqR;

  rhat = qqq/2;
  uhat = (qqL*qL[1] + qqR*qR[1] )/qqq;
  vhat = (qqL*qL[2] + qqR*qR[2] )/qqq;
  ehat = (qqL*qL[3] + qqR*qR[3] )/qqq;
  vhat2 = uhat*uhat+vhat*vhat;
  ahat = pow_fast((s->kappa-1)*(ehat-vhat2/2),0.5);

  /*  printf("q^ %7.5g %7.5g %7.5g %7.5g %7.5g \n",      
	 rhat,uhat,vhat,ehat,ahat);
  */
  rlam[0] = uhat-ahat;
  rlam[1] = uhat;
  rlam[2] = uhat;
  rlam[3] = uhat+ahat;

  /*
  printf("Eigenvals: %7.5g  %7.5g  %7.5g  %7.5g \n",
	 rlam[0],rlam[1],rlam[2],rlam[3]);
  */

  r[0][0] = 1.;
  r[0][1] = uhat-ahat;
  r[0][2] = vhat;
  r[0][3] = ehat-uhat*ahat;

  r[1][0] = 1.;
  r[1][1] = uhat;
  r[1][2] = vhat;
  r[1][3] = vhat2/2;

  r[2][0] = 0.;
  r[2][1] = 0.;
  r[2][2] = 1.;
  r[2][3] = vhat;

  r[3][0] = 1.;
  r[3][1] = uhat+ahat;
  r[3][2] = vhat;
  r[3][3] = ehat+uhat*ahat;

  for(k=0;k<4;k++)
    dw[k] = wj_rot[k] - wi_rot[k];
  
  dw[4] = dw[3] -(dw[2]-vhat*dw[0])*vhat;

  gamma[1] = (s->kappa-1)/ahat/ahat*
    (dw[0]*(ehat-uhat*uhat) + uhat*dw[1]-dw[4]);
  gamma[0] = 0.5/ahat*(dw[0]*(uhat+ahat) - dw[1] - ahat*gamma[1]);

  gamma[3] = dw[0] - gamma[0] - gamma[1];
  gamma[2] = dw[2] - vhat*dw[0];

  /*
  printf("Gammas: %7.5g  %7.5g  %7.5g  %7.5g \n",
	 gamma[0],gamma[1],gamma[2],gamma[3]);
  */

  //  compute_pressure(pres,wi_rot,s->kappa);
  pres = (wi[3]-0.5*(wi[1]*wi[1]+wi[2]*wi[2])/wi[0])*(s->kappa-1);

  /*
    printf("Pressure= %7.5g: %7.5g, %7.5g, %7.5g, %7.5g\n",
	 pres,wi[0],wi[1],wi[2],wi[3]);
  */
  f1[0] = qL[1];
  f1[1] = qL[1]*qL[1]*qL[0]+pres;
  f1[2] = qL[1]*qL[2]*qL[0];
  f1[3] = (pres + wi[3])*qL[1]/wi[0];

  /*
  printf("Flux_ex = %7.5g: %7.5g, %7.5g, %7.5g, %7.5g\n",
	 pres,f1[0],f1[1],f1[2],f1[3]);
  */

  for(k=0;k<4;k++){
    f_roe[k] = f1[k];
    /*    if(k == 0)
      printf("0th %7.5g %7.5g \n",
	     f_roe[k],pres);
    */
    for(j=0;j<4;j++){
      if(rlam[j] <= 0.){
	f_roe[k]+=  gamma[j]*rlam[j]*r[j][k];
	/*
	if(k == 0)
	  printf("%dth %7.5g %7.5g %7.5g %7.5g %7.5g  \n",j+1,
	       f_roe[k],gamma[j]*rlam[j]*r[j][k],gamma[j],rlam[j],r[j][k]);
	*/
      }
    }
  }

  H[0] = f_roe[0];
  H[1] = f_roe[1]*n[0] - f_roe[2]*n[1];
  H[2] = f_roe[1]*n[1] + f_roe[2]*n[0];
  H[3] = f_roe[3];

}

void Solver_Flux_Batten(const tSolver_State *s,const double4 wi,
  const double4 wj,const double2 n,double4 H){
  #define eps1 1e-6
  #define eps2 1e-6
  double alpha=1,rrl,rrr,ql,qr,pl,pr,pl1,pr1,pl2,pr2,cl,cr,ps,lconst,rconst,
    c1l,c1r,ml,mr,pn,newml,newmr,t1,t2,sm,ut,vt,wt,sv,ri,qi,pi,eexp,rstar,
    sexp,uu,vv,ww;
  char conv=0;
  uint i;
  rrl=1/wi[0];
  rrr=1/wj[0];
  ql=rrl*(n[0]*wi[1]+n[1]*wi[2]);
  qr=rrr*(n[0]*wj[1]+n[1]*wj[2]);
  pl=(s->kappa-1)*(wi[3]-0.5*rrl*(wi[1]*wi[1]+wi[2]*wi[2]));
  pr=(s->kappa-1)*(wj[3]-0.5*rrr*(wj[1]*wj[1]+wj[2]*wj[2]));
  if(pl <=0 || pr <=0 ){ 
    printf("ZERO pressure %7.5g %7.5g\n ",pl,pr);
    printf(" %7.5g %7.5g  %7.5g %7.5g  %7.5g %7.5g %7.5g %7.5g  \n",
	   wj[0],wj[1],wj[2],wj[3],
	   0.5*rrr*(wj[1]*wj[1]+wj[2]*wj[2]),
	   wj[3]-0.5*rrr*(wj[1]*wj[1]+wj[2]*wj[2]),rrr,(s->kappa-1) );
  }
  if(pl <=0) pl=0.001;
  if(pr <=0) pr=0.001;
  pl1=1/pl;
  pr1=1/pr;
  cl=sqrt(s->kappa*pl*rrl);
  cr=sqrt(s->kappa*pr*rrr);
  ps=0.5*(pl+pr);
  lconst=sqrt(wi[0]*pl);  /* start of exact Riemann solver */
  rconst=sqrt(wj[0]*pr);
  c1l=(0.4/2.3664319)*lconst;
  c1r=(0.4/2.3664319)*rconst;
  pl2=ps*pl1;
  pr2=ps*pr1;
  if(pl2<0.999)
    ml=c1l*(1-pl2)/(1-pow_fast(pl2,0.14285714));
  else
    ml=lconst*sqrt(1.2*pl2+0.2);
  if(pr2<0.999)
    mr=c1r*(1-pr2)/(1-pow_fast(pr2,0.14285714));
  else
    mr=rconst*sqrt(1.2*pr2+0.2);
  do{  /* Godunov's iteration */
    for(i=0;i<20;i++){
      pn=(pl*mr+pr*ml+(ql-qr)*ml*mr)/(ml+mr);
      if(pn<eps1)
        pn=eps1;
      ps*=1-alpha;
      ps+=alpha*pn;
      pl2=ps*pl1;
      pr2=ps*pr1;
      if(pl2<0.999)
        newml=c1l*(1-pl2)/(1-pow_fast(pl2,0.14285714));
      else
        newml=lconst*sqrt(1.2*pl2+0.2);
      if(pr2<0.999)
        newmr=c1r*(1-pr2)/(1-pow_fast(pr2,0.14285714));
      else
        newmr=rconst*sqrt(1.2*pr2+0.2);
      t1=newml-ml;
      if(t1<0)
        t1=-t1;
      t2=newmr-mr;
      if(t2<0)
        t2=-t2;
      ml=newml;
      mr=newmr;
      if((t1<eps1&&t2<eps1)||alpha<eps2){
        conv=1;
        break;
        }
      }
    alpha*=0.8;
    }while(!conv);
  sm=(pl-pr+ml*ql+mr*qr)/(ml+mr);  /* now find what is happening at the interface */
  if(sm<0){  /* right of the contact */
    ut=wj[1]*rrr-qr*n[0];
    vt=wj[2]*rrr-qr*n[1];
    if(ps>pr){	/* right travelling shock */
      sv=qr+mr*rrr;
      if(sv<0){	 /* supersonic from right to left */
        ri=wj[0];
        qi=qr;
        pi=pr;
        }
      else{  /* between contact and right moving shock */
        ri=mr/(sv-sm);
        qi=sm;
        pi=ps;
        }
      }
    else{  /* right travelling expansion */
      eexp=qr+cr;
      if(eexp<0){  /* supersonic from right to left */
        ri=wj[0];
        qi=qr;
        pi=pr;
        }
      else{
	rstar=wj[0]*exp(log(ps*pr1)/s->kappa);
        if(s->kappa*ps>rstar*sm*sm){  /* between contact and start of expansion */
          ri=rstar;
          qi=sm;
          pi=ps;
          }
        else{  /* in the expansion wave */
          qi=((s->kappa-1)*qr-2*cr)/(s->kappa+1);
          ri=wj[0]*pow_fast(1+0.5*(s->kappa-1)*(qi-qr)/cr,2/(s->kappa-1));
          pi=pr*pow_fast(ri*rrr,s->kappa);
          }
        }
      }
    }
  else{  /* left of the contact */
    ut=wi[1]*rrl-ql*n[0];
    vt=wi[2]*rrl-ql*n[1];
    if(ps>pl){  /* left travelling shock */
      sv=ql-ml*rrl;
      if(sv>0){  /* supersonic from left to right */
        ri=wi[0];
        qi=ql;
        pi=pl;
        }
      else{  /* between shock and the contact */
        ri=ml/(sm-sv);
        qi=sm;
        pi=ps;
        }
      }
    else{  /* left travelling expansion */
      sexp=ql-cl;
      if(sexp>0){  /* supersonic from left to right */
        ri=wi[0];
        qi=ql;
        pi=pl;
        }
      else{
	rstar=wi[0]*exp(log(ps*pl1)/s->kappa);
        if(rstar*sm*sm>s->kappa*ps){  /* in the expansion wave */
          qi=((s->kappa-1)*ql+2*cl)/(s->kappa+1);
          ri=wi[0]*pow_fast(1-0.5*(s->kappa-1)*(qi-ql)/cl,2/(s->kappa-1));
          pi=pl*pow_fast(ri*rrl,s->kappa);
          }
        else{  /* between expansion and contact */
          ri=rstar;
          qi=sm;
          pi=ps;
          }
        }
      }
    }
  uu=n[0]*qi+ut;
  vv=n[1]*qi+vt;
  H[0]=ri*qi;
  H[1]=(H[0]*uu+pi*n[0]);
  H[2]=(H[0]*vv+pi*n[1]);
  H[3]=(pi/(s->kappa-1)+0.5*ri*(uu*uu+vv*vv)+pi)*qi;
  }

double Compute_Outlet_Pressure(tSolver_State *s){
  uint i,_i,_ii;
  double ptot,length;

  /*   computing of the mean value of pressure at outlet */
  ptot = 0;
  length=0;
  for(i=0;i<s->nface;i++){
    _i=s->face[i].i;
    _ii= s->face[i].ii;
    if(( s->face[i].boundary>0) &&
       (s->boundary[s->face[i].boundary-1].type==
        cSolver_BoundaryType_Inlet ||
	s->boundary[s->face[i].boundary-1].type==
        cSolver_BoundaryType_Outlet )){
      double oneoverrho,norm,product,p;
      product=s->face[i].n[0]*s->element[_i].w[_ii][1]+
	s->face[i].n[1]*s->element[_i].w[_ii][2];
      if(0<product){
	oneoverrho=1/s->element[_i].w[_ii][0];
	norm=s->element[_i].w[_ii][1]*s->element[_i].w[_ii][1]+
	s->element[_i].w[_ii][2]*s->element[_i].w[_ii][2];
	p = (s->kappa-1)*(s->element[_i].w[_ii][3]-0.5*norm*oneoverrho);
	ptot+=p*s->face[i].volume;
	length+=s->face[i].volume;
      }
    }
  }
  return ptot/length;
}

void Solver_ComputeFluxes(tSolver_State *s){
  uint i,j,k;
  const int itest = -52, itest1 = -118;
  const double oneoversqrt3 = 0.57735026918962584;
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};
  int ic;
  double p_out;

  /* p_out = Compute_Outlet_Pressure(s); */

  for(i=0;i<s->nface;i++){
    uint _i,_j,_ii,_jj;
    double4 wi1, wj1, H1, wi2, wj2, H2, wpom;  /*two integration nodes */
    double4 H3,H4; 
    double rez1,rez2;
    double2 Hper;  /* keeps rotated H[1]..H[2] for periodic b.c. */
    double4 q0,q1,q2,qi1,qi2;
    
    /*printf("\n Solver_ComputeFluxes: %d %d %d \n",i,s->nface, 
       s->face[i].boundary); */

    _i=s->face[i].i;
    _ii= s->face[i].ii;

/*     if(_i == itest1) 
       for(j=0;j<4;j++)
       printf("Fluxes initial %d %d %d %7.5g %7.5g %7.5g \n",
          _i,perm[_ii][0],j,s->element[_i].H[0][j],
        s->element[_i].H[1][j],s->element[_i].H[2][j]); */

    double4_sub(wpom,s->element[_i].w[perm[_ii][2]],
              s->element[_i].w[perm[_ii][1]]);
    double4_multiply(wpom, oneoversqrt3, wpom);
    double4_add(wi1,s->element[_i].w[perm[_ii][0]],wpom);
    double4_sub(wi2,s->element[_i].w[perm[_ii][0]],wpom);


    /* recomputing in a new way  */
    /*
    vector_1_w2q(s->element[_i].w[perm[_ii][0]],q0,s->kappa);
    vector_1_w2q(s->element[_i].w[perm[_ii][1]],q1,s->kappa);
    vector_1_w2q(s->element[_i].w[perm[_ii][2]],q2,s->kappa);

    double4_sub(wpom,q2,q1);
    double4_multiply(wpom, oneoversqrt3, wpom);
    double4_add(qi1,q0,wpom);
    double4_sub(qi2,q0,wpom);
    vector_1_q2w(qi1,wi1,s->kappa);
    vector_1_q2w(qi2,wi2,s->kappa);
    */

    if(s->face[i].boundary<=0){  /* inside or periodic */

      _j=s->face[i].j;
      _jj=s->face[i].jj;
      if(!s->face[i].boundary ||s->face[i].boundary ){  /* inside */
	/* inner as well as periodic */
       /*printf("Inner \n");*/

       double4_sub(wpom,s->element[_j].w[perm[_jj][2]],
              s->element[_j].w[perm[_jj][1]]);
       double4_multiply(wpom, oneoversqrt3, wpom);
       double4_add(wj1,s->element[_j].w[perm[_jj][0]],wpom);
       double4_sub(wj2,s->element[_j].w[perm[_jj][0]],wpom);

	/* recomputing in a new way  */
       /*
	vector_1_w2q(s->element[_j].w[perm[_jj][0]],q0,s->kappa);
	vector_1_w2q(s->element[_j].w[perm[_jj][1]],q1,s->kappa);
	vector_1_w2q(s->element[_j].w[perm[_jj][2]],q2,s->kappa);
	double4_sub(wpom,q2,q1);
	double4_multiply(wpom, oneoversqrt3, wpom);
	double4_add(qi1,q0,wpom);
	double4_sub(qi2,q0,wpom);
	vector_1_q2w(qi1,wj1,s->kappa);
	vector_1_q2w(qi2,wj2,s->kappa);
       */
         
       /*
       if(s->face[i].boundary<0){  // periodic 
	  printf("i=%d, ii=%d, j=%d, jj=%d, nn=(%7.5g, %7.5g), ibc=%d\n", 
	   _i,_ii,_j,_jj,
		 s->face[i].n[0]*s->face[i].volume,
		 s->face[i].n[1]*s->face[i].volume, s->face[i].boundary); 
	  printf("Per1: i=%d, ii=%d, x=(%10.5g, %10.5g)\n",
		 _i,_ii, 
		 s->element[_i].ve->P[perm[_ii][0]]->Point[0],
		 s->element[_i].ve->P[perm[_ii][0]]->Point[1]);
	  printf("Per2: j=%d, jj=%d, x=(%10.5g, %10.5g)\n",
		 _j,_jj, 
		 s->element[_j].ve->P[perm[_jj][1]]->Point[0],
		 s->element[_j].ve->P[perm[_jj][1]]->Point[1]);
       }
       */
       if(_i == itest || _j == itest) {
       /*       if( _i == 2153 && _j == 1999){ */
	 printf("i ii, j jj,, n_ij %d %d, %d %d, %7.5g %7.5g %d\n", 
	   _i,_ii,_j,_jj,
	s->face[i].n[0]*s->face[i].volume ,s->face[i].n[1]*s->face[i].volume,
           s->face[i].boundary); 

       printf("wi1 wi2 %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",
         wi1[0],wi1[1],wi1[2],wi1[3], wi2[0],wi2[1],wi2[2],wi2[3]);
       printf("wj1 wj2 %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",
         wj1[0],wj1[1],wj1[2],wj1[3], wj2[0],wj2[1],wj2[2],wj2[3]);

       }

       /*       wi1[0] = 1.;
       wi1[1] = 2.;
       wi1[2] = 0.;
       wi1[3] = 1.5;
       wj2[0] = 1.1;
       wj2[1] = 1.1;
       wj2[2] = -0.1;
       wj2[3] = 1.05;
       
       s->face[i].n[0] = 0.6;
       s->face[i].n[1] = 0.8;
       
       printf("wi1 wi2 %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",
	      wi1[0],wi1[1],wi1[2],wi1[3], wi2[0],wi2[1],wi2[2],wi2[3]);
       printf("wj1 wj2 %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",
	      wj1[0],wj1[1],wj1[2],wj1[3], wj2[0],wj2[1],wj2[2],wj2[3]);

       s->flux(s,wi1,wj2,s->face[i].n,H1); 
       printf("FLUXES: %7.5g %7.5g %7.5g %7.5g  \n",
	        H1[0],H1[1],H1[2],H1[3]); 
       break;
       printf("FLUXES: %4d %4d %7.5g %7.5g %7.5g %7.5g %7.5g  \n",
	       _i,_j, H1[0],H1[1],H1[2],H1[3],s->face[i].volume); 
       */

       s->flux(s,wi1,wj2,s->face[i].n,H1); 
       s->flux(s,wi2,wj1,s->face[i].n,H2); 

       /*
	 s->flux1(s,wi1,wj2,s->face[i].n,H3); 
	 s->flux1(s,wi2,wj1,s->face[i].n,H4); 
	 
	 double4_norm(rez1,H1,H3);
	 double4_norm(rez2,H2,H4);
       */
       /*
	 if(rez1 > 1E-06 || rez2 > 1E-06){
	 printf("FLUXES: %4d %4d %7.5g %7.5g %7.5g %7.5g %7.5g  \n",
	 _i,_j, H1[0],H1[1],H1[2],H1[3],rez1); 
	 printf("FLUXES: %4d %4d %7.5g %7.5g  %7.5g %7.5g  \n",
	 _i,_j, H3[0],H3[1],H3[2],H3[3]); 
	 printf("FLUXES: %4d %4d %7.5g %7.5g %7.5g %7.5g %7.5g  \n",
	 _i,_j, H2[0],H2[1],H2[2],H2[3],rez2); 
	 printf("FLUXES: %4d %4d %7.5g %7.5g  %7.5g %7.5g  \n",
	 _i,_j, H4[0],H4[1],H4[2],H4[3]); 
	 printf("_____________________________________\n");
	 
	 }
       */
       if(_i == itest || _j == itest){ 
	 /* if( _i == 2153 && _j == 1999){  */
	 /*      if( _i <2005 && _i > 1995 || ( _j <2005 && _j > 1995)){ */
	printf("FLUXES: %4d %4d %7.5g %7.5g %7.5g %7.5g %7.5g  \n",
	       _i,_j, H1[0],H1[1],H1[2],H1[3],s->face[i].volume); 
	printf("FLUXES: %4d %4d %7.5g %7.5g %7.5g %7.5g %7.5g  \n",
	       _i,_j, H2[0],H2[1],H2[2],H2[3],s->face[i].volume); 
        printf("\n");
       }

  }
      else{  /* old periodic */
        double4 w;
        double2 v1;
        tMesh_VolumeElement *vei,*vej;
        tMesh_Point *p1,*p2,*p3;
        uint k;
        w[0]=s->element[_j].w[0][0];
        w[3]=s->element[_j].w[0][3];
        v1[0]=s->element[_j].w[0][1];
        v1[1]=s->element[_j].w[0][2];
        vei=s->element[_i].ve;
        vej=s->element[_j].ve;

	 printf("\n\n\n #############################\n\n\n\n");
/*        printf("Per: %10.5g %10.5g %10.5g %10.5g\n",
        s->element[_i].w[0][0],s->element[_i].w[0][1],s->element[_i].w[0][2],
           s->element[_i].w[0][3]);
        printf("Per: %10.5g %10.5g %10.5g %10.5g\n",
          w[0],v1[0],v1[1],w[3]);  */

        for(k=0;vei->Neighbour[k]!=vej;k++);
        p1=vei->P[perm[k][0]];
        p2=vei->P[perm[k][1]];

/*        printf("Per: %10.5g %10.5g %10.5g %10.5g\n",
         p1->Point[0],p1->Point[1],
          p1->PeriodicCopy->Point[0],p1->PeriodicCopy->Point[1]); */

        w[1] = v1[0]; w[2] = v1[1];

/*        printf("Vel: %10.5g %10.5g %10.5g %10.5g\n",
        v1[0],v1[1],w[1],w[2]); */

        s->flux(s,s->element[_i].w[0],w,s->face[i].n,H1);

/*        printf("Flu: %10.5g %10.5g %10.5g %10.5g\n",
           H[0],H[1],H[2],H[3]); */
/*
        Triangle_VectorRotation(p1->Point,p2->Point,p3->Point,
          p1->PeriodicCopy->Point,p2->PeriodicCopy->Point,
          Mesh_EdgeOpposite(vej,p1->PeriodicCopy,p2->PeriodicCopy)->Point,
          H+1,Hper);  rotate result back  */
    } 
    }
    else{          /* boundary */
       if(s->boundary[s->face[i].boundary-1].type==
        cSolver_BoundaryType_Inlet ||
	 s->boundary[s->face[i].boundary-1].type==
	  cSolver_BoundaryType_Outlet ){  /* inlet/outlet */
        /*double4 w; */
        /*printf("Inlet/outlet \n"); */
	 if(s->boundary[s->face[i].boundary-1].fixed){  /* fixed */
	   /*
	     for(j=0;j<4;j++){
	     wj1[j]=s->boundary[s->face[i].boundary-1].w[j];
	     wj2[j]=s->boundary[s->face[i].boundary-1].w[j];
	   }
	   */
	   /* RINGLEB FLOW */
	   tMesh_VolumeElement *ve;
	   tMesh_Point *p1,*p2;
	   double a[2];
	   double2 xp[2];
	   double4 qj1,qj2;

	   a[0] = (1/sqrt(3) + 1)/2;
	   a[1] = (-1/sqrt(3) + 1)/2;
    
	   _i=s->face[i].i;
	   _ii= s->face[i].ii;
	   ve=s->element[_i].ve;

	   /* points forming the boundary edge */
	   p1=ve->P[perm[_ii][0]];
	   p2=ve->P[perm[_ii][1]];
	   /*
	     printf("%7.5g %7.5g \n",p1->Point[0],p1->Point[1]);
	     printf("%7.5g %7.5g \n",p2->Point[0],p2->Point[1]);
	   */
	   
	   /* computation of integration nodes */
	   double2_weight_awer(xp[0],p1->Point,p2->Point,a[0],a[1]);
	   double2_weight_awer(xp[1],p1->Point,p2->Point,a[1],a[0]);
	   

	   /* computation of the solution in integration nodes */
	   calc_ringleb(s->kappa,xp[0],qj1);
	   calc_ringleb(s->kappa,xp[1],qj2);
	   /* transform form physival quantities to the state vector */
	   vector_q2w(qj1,wj1,s->kappa);
	   vector_q2w(qj2,wj2,s->kappa);

	   /*
	     if(s->face[i].boundary-1 == 3){
	     printf("%7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",
		    xp[0][0],xp[0][1],qj1[1],qj1[2],
		    s->face[i].n[0],s->face[i].n[1],
		    qj1[1]*s->face[i].n[0]+qj1[2]*s->face[i].n[1]);
	     printf("%7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n\n",
		    xp[1][0],xp[1][1],qj2[1],qj2[2],
		    s->face[i].n[0],s->face[i].n[1],
		    qj2[1]*s->face[i].n[0]+qj2[2]*s->face[i].n[1]);
		    }
	   */


	   /*
	   printf("%9.7g %9.7g %9.7g %9.7g \n",qj1[0],qj1[1],qj1[2],qj1[3]);
	   printf("%9.7g %9.7g %9.7g %9.7g\n---\n",qj2[0],qj2[1],qj2[2],qj2[3]);
	   printf("%9.7g %9.7g %9.7g %9.7g \n",wj1[0],wj1[1],wj1[2],wj1[3]);
	   printf("%9.7g %9.7g %9.7g %9.7g\n---\n",wj2[0],wj2[1],wj2[2],wj2[3]);
	   */
	 }
	 else{  /* physical */
          double oneoverrho,norm,c2,product,product2;
          oneoverrho=1/s->element[_i].w[_ii][0];
          norm=s->element[_i].w[_ii][1]*s->element[_i].w[_ii][1]+
            s->element[_i].w[_ii][2]*s->element[_i].w[_ii][2];
          c2=s->kappa2*oneoverrho*(s->element[_i].w[_ii][3]-
             0.5*oneoverrho*norm);
          product=oneoverrho*(s->face[i].n[0]*s->element[_i].w[_ii][1]+
            s->face[i].n[1]*s->element[_i].w[_ii][2]);
          product2=product*product;
          if(0<product&&c2<=product2)  /* supersonic outlet */
            for(j=0;j<4;j++){
              /*w[j]=s->element[_i].w[_ii][j];*/
                wj1[j] = wi1[j];
                wj2[j] = wi2[j];
            }
          else
            if(0<product){  /* subsonic outlet */
              for(j=0;j<4;j++){
                /* w[j]=s->element[_i].w[_ii][j]; */
                wj1[j] = wi1[j];
                wj2[j] = wi2[j];
              }

              wj1[3]=s->boundary[s->face[i].boundary-1].w[3]+
                0.5*((wj1[1]*wj1[1]+wj1[2]*wj1[2])/wj1[0]-
                (s->boundary[s->face[i].boundary-1].w[1]*
                s->boundary[s->face[i].boundary-1].w[1]+
                s->boundary[s->face[i].boundary-1].w[2]*
                s->boundary[s->face[i].boundary-1].w[2])/
                s->boundary[s->face[i].boundary-1].w[0]);

              wj2[3]=s->boundary[s->face[i].boundary-1].w[3]+
                0.5*((wj2[1]*wj2[1]+wj2[2]*wj2[2])/wj2[0]-
                (s->boundary[s->face[i].boundary-1].w[1]*
                s->boundary[s->face[i].boundary-1].w[1]+
                s->boundary[s->face[i].boundary-1].w[2]*
                s->boundary[s->face[i].boundary-1].w[2])/
                s->boundary[s->face[i].boundary-1].w[0]);
              }
            else
              if(product2<c2){  /* subsonic inlet */
                for(j=0;j<3;j++){
                  wj1[j]=s->boundary[s->face[i].boundary-1].w[j];
                  wj2[j]=s->boundary[s->face[i].boundary-1].w[j];
                }

                wj1[3]=wi1[3]+ 0.5*((wj1[1]*wj1[1]+wj1[2]*wj1[2])/wj1[0]-
                  (wi1[1]*wi1[1] + wi1[2]*wi1[2])/wi1[0]);
                wj2[3]=wi2[3]+ 0.5*((wj2[1]*wj2[1]+wj2[2]*wj2[2])/wj2[0]-
                  (wi2[1]*wi2[1] + wi2[2]*wi2[2])/wi2[0]);

/*                printf("BC: %d %7.5g %7.5g %7.5g %7.5g \n",
                   i,wj1[0],wj1[1],wj1[2],wj1[3]);
                printf("wi: %d %7.5g %7.5g %7.5g %7.5g \n",
                   i,wi1[0],wi1[1],wi1[2],wi1[3]); */
                }
              else  /* supersonic inlet */
                for(j=0;j<4;j++){
                  wj1[j]=s->boundary[s->face[i].boundary-1].w[j];
                  wj2[j]=s->boundary[s->face[i].boundary-1].w[j];
                }
          }
        s->flux(s,wi1,wj1,s->face[i].n,H1);
        s->flux(s,wi2,wj2,s->face[i].n,H2);
	/*
	 s->flux1(s,wi1,wj2,s->face[i].n,H3); 
	 s->flux1(s,wi2,wj1,s->face[i].n,H4); 
	 
	 double4_norm(rez1,H1,H3);
	 double4_norm(rez2,H2,H4);
	*/
	/*
	 if(rez1 > 1E-08 || rez2 > 1E-08){
	 printf("FLUXES: %4d %4d %7.5g %7.5g %7.5g %7.5g %7.5g  \n",
	 _i,_j, H1[0],H1[1],H1[2],H1[3],rez1); 
	 printf("FLUXES: %4d %4d %7.5g %7.5g  %7.5g %7.5g  \n",
	 _i,_j, H3[0],H3[1],H3[2],H3[3]); 
	 printf("FLUXES: %4d %4d %7.5g %7.5g %7.5g %7.5g %7.5g  \n",
	 _i,_j, H2[0],H2[1],H2[2],H2[3],rez2); 
	 printf("FLUXES: %4d %4d %7.5g %7.5g  %7.5g %7.5g  \n",
	 _i,_j, H4[0],H4[1],H4[2],H4[3]); 
	 printf("_____________________________________\n");
	 
	 }
	*/
        }
       else{  /* slip */
        double p;

	   /*	   if(s->face[i].boundary-1 == 3 && ic >= 0){
	     printf("%7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",
		    xp[0][0],xp[0][1],qj1[1],qj1[2],
		    s->curvedelement[ic].normal[0][0],
		    s->curvedelement[ic].normal[0][1],
		    qj1[1]*s->curvedelement[ic].normal[0][0]
		    +qj1[2]*s->curvedelement[ic].normal[0][1]);
	     printf("%7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n\n",
		    xp[1][0],xp[1][1],qj2[1],qj2[2],
		    s->curvedelement[ic].normal[1][0],
		    s->curvedelement[ic].normal[1][1],
		    qj2[1]*s->curvedelement[ic].normal[1][0]
		    +qj2[2]*s->curvedelement[ic].normal[1][1]);
	   }
	   */

	   /*
	   printf("%9.7g %9.7g %9.7g %9.7g \n",qj1[0],qj1[1],qj1[2],qj1[3]);
	   printf("%9.7g %9.7g %9.7g %9.7g\n---\n",qj2[0],qj2[1],qj2[2],qj2[3]);
	   printf("%9.7g %9.7g %9.7g %9.7g \n",wj1[0],wj1[1],wj1[2],wj1[3]);
	   printf("%9.7g %9.7g %9.7g %9.7g\n---\n",wj2[0],wj2[1],wj2[2],wj2[3]);
	   */

        /*printf("Wall \n"); */
	/* printf("%7.5g %7.5g Wall \n",
	      s->element[_i].ve->P[perm[_ii][2]]->Point[0],
	      s->element[_i].ve->P[perm[_ii][2]]->Point[1] ); */
	ic=s->element[_i].icve;	
	if(ic >= 0){  /*curved elements  */
         /*	  printf("Normals %d %d %d \n",i,_i,ic); */
         /* printf("Normals: %7.5g  %7.5g  %7.5g  %7.5g %7.5g  %7.5g  \n",
		 s->face[i].n[0],s->face[i].n[1], 
		 s->curvedelement[ic].normal[0][0],
		 s->curvedelement[ic].normal[0][1],
		 s->curvedelement[ic].normal[1][0],
		 s->curvedelement[ic].normal[1][1]); */

        /*   printf("%7.5g  %7.5g %7.5g  %7.5g  \n",
           s->curvedelement[ic].xb[0],s->curvedelement[ic].xb[1],
             s->face[i].n[0],s->face[i].n[1]);
          printf("%7.5g  %7.5g  \n \n",
           s->curvedelement[ic].xb[0]+s->face[i].n[0]*s->face[i].volume,
           s->curvedelement[ic].xb[1]+s->face[i].n[1]*s->face[i].volume);
          printf("%7.5g  %7.5g  %7.5g  %7.5g  \n",
           s->curvedelement[ic].xb[0]-0.01,s->curvedelement[ic].xb[1],
          s->curvedelement[ic].normal[0][0],s->curvedelement[ic].normal[0][1]);
          printf("%7.5g  %7.5g  \n\n",
           s->curvedelement[ic].xb[0]-0.01+
            s->curvedelement[ic].normal[0][0]*s->curvedelement[ic].length[0],
           s->curvedelement[ic].xb[1]+
            s->curvedelement[ic].normal[0][1]*s->curvedelement[ic].length[0]);
          printf("%7.5g  %7.5g  %7.5g  %7.5g  \n",
           s->curvedelement[ic].xb[0]+0.01,s->curvedelement[ic].xb[1],
          s->curvedelement[ic].normal[1][0],s->curvedelement[ic].normal[1][1]);
          printf("%7.5g  %7.5g  \n",
           s->curvedelement[ic].xb[0]+0.01+
             s->curvedelement[ic].normal[1][0]*s->curvedelement[ic].length[1],
           s->curvedelement[ic].xb[1]+
             s->curvedelement[ic].normal[1][1]*s->curvedelement[ic].length[1]);
          printf("-------------- \n");
          printf(" %7.5g %7.5g \n",
          s->element[_i].ve->P[0]->Point[0],s->element[_i].ve->P[0]->Point[1]);
          printf(" %7.5g %7.5g \n",
          s->element[_i].ve->P[1]->Point[0],s->element[_i].ve->P[1]->Point[1]);
          printf(" %7.5g %7.5g \n",
          s->element[_i].ve->P[2]->Point[0],s->element[_i].ve->P[2]->Point[1]);
          printf("-------------- \n");  */


	  p=(s->kappa-1)*(wi1[3]-0.5*(wi1[1]*wi1[1]+wi1[2]*wi1[2])/wi1[0]);
	  H1[0]=0;
	  H1[1]=p*s->curvedelement[ic].normal[0][0];
	  H1[2]=p*s->curvedelement[ic].normal[0][1];
	  H1[3]=0;
	 if(_i == itest) 
	    printf("WALLS1: %d %7.5g   %7.5g  %7.5g  %7.5g  %7.5g  \n",
		   _i,p,s->curvedelement[ic].normal[0][0],
		   s->curvedelement[ic].normal[0][1],H1[1],H1[2]);

	  p=(s->kappa-1)*(wi2[3]-0.5*(wi2[1]*wi2[1]+wi2[2]*wi2[2])/wi2[0]);
	  H2[0]=0;
	  H2[1]=p*s->curvedelement[ic].normal[1][0];
	  H2[2]=p*s->curvedelement[ic].normal[1][1];
	  H2[3]=0;
	  
	   if(_i == itest) 
	    printf("WALLS2: %d %7.5g   %7.5g  %7.5g  %7.5g  %7.5g  \n",
		   _i,p,s->curvedelement[ic].normal[1][0],
		   s->curvedelement[ic].normal[1][1],H2[1],H2[2]);
	}
	else{   /*non curved elements */
	  p=(s->kappa-1)*(wi1[3]-0.5*(wi1[1]*wi1[1]+wi1[2]*wi1[2])/wi1[0]);
	  H1[0]=0;
	  H1[1]=p*s->face[i].n[0];
	  H1[2]=p*s->face[i].n[1];
	  H1[3]=0;
	  p=(s->kappa-1)*(wi2[3]-0.5*(wi2[1]*wi2[1]+wi2[2]*wi2[2])/wi2[0]);
	  H2[0]=0;
	  H2[1]=p*s->face[i].n[0];
	  H2[2]=p*s->face[i].n[1];
	  H2[3]=0;
	}
      }
    }

    /* multiplication by the edge length and the weight =0.5 */
    ic = s->element[_i].icve;
    if(ic >=0 &&s->boundary[s->face[i].boundary-1].type==
       cSolver_BoundaryType_Slip )  /*curved edge */
      for(j=0;j<4;j++){
	H1[j]*=s->curvedelement[ic].length[0]*0.5;
	H2[j]*=s->curvedelement[ic].length[1]*0.5;
      }
    else  /* noncurved edges */
      for(j=0;j<4;j++){
	H1[j]*=s->face[i].volume*0.5;
	H2[j]*=s->face[i].volume*0.5;
      }

    /* combination of fluxes H with the values of the test functions */
    for(j=0;j<4;j++){
      s->element[_i].H[perm[_ii][0]][j]+=H1[j]+H2[j];
      s->element[_i].H[perm[_ii][1]][j]+=(-H1[j]+H2[j])*oneoversqrt3;
      s->element[_i].H[perm[_ii][2]][j]+=( H1[j]-H2[j])*oneoversqrt3;
      /*      if(_i == itest1) */
      if( _i <2005 && _i > 2995)
	printf("Fluxes adding1 %d %d %d %7.5g %7.5g %7.5g \n",
	       _i,perm[_ii][0],j,
       s->element[_i].H[0][j],s->element[_i].H[1][j],s->element[_i].H[2][j]);
    }

    if(s->face[i].boundary<=0){
      /*      if(s->face[i].boundary==-1){  // periodic fix
        for(j=0;j<2;j++)
          H[j+1]=Hper[j]*s->face[i].volume;
      } */
    /* combination of fluxes H with the values of the test functions */
      for(j=0;j<4;j++){  /* opisite sign of outer vector */
        s->element[_j].H[perm[_jj][0]][j]-=H1[j]+H2[j];
        s->element[_j].H[perm[_jj][1]][j]-=( H1[j]-H2[j])*oneoversqrt3;
        s->element[_j].H[perm[_jj][2]][j]-=(-H1[j]+H2[j])*oneoversqrt3;

	if(_j == itest1)
	/*	if( _j <2005 && _j > 1995) */
	  printf("Fluxes adding2 %d %d %d %7.5g %7.5g %7.5g \n",
         _j,perm[_jj][0],j,
	 s->element[_j].H[0][j],s->element[_j].H[1][j],s->element[_j].H[2][j]);

      }
     } 
    }
  /* division by the area of element, from the time derivative term */
  /* multiplication by s->DOF, the weights from the time derivative term */
  for(i=0;i<s->nelem;i++){
    int ic;
    ic=s->element[i].icve;  
    if(ic >= 0){  /* curved elements */
      for(k=0;k<s->DOF;k++)
	for(j=0;j<4;j++)
	  s->element[i].H[k][j]*=2/s->curvedelement[ic].detJF[k]*s->DOF; 
         /*   printf("Volumes %d %7.5g  %7.5g %7.5g %7.5g \n",i,
	     s->element[i].oneovervolume,2/s->curvedelement[ic].detJF[0],
         2/s->curvedelement[ic].detJF[1],2/s->curvedelement[ic].detJF[2]); */
      if(i == itest1) 
	/*      if( i <2005 && i > 1995) */
	printf("Fluxes dividing1 %d  %d %7.5g %7.5g %7.5g \n",
           i,perm[i][0],
	   s->element[i].H[0][0],s->element[i].H[1][0],s->element[i].H[2][0]);
    }
    else{
      for(k=0;k<s->DOF;k++) /* non curved elements */
	for(j=0;j<4;j++)
	  s->element[i].H[k][j]*=s->element[i].oneovervolume*s->DOF; 
      if(i == itest1)
	/*  if( i <2005 && i > 1995) */
          printf("Fluxes dividing2  %d %d %7.5g %7.5g %7.5g %7.5g %d\n",
          i,perm[i][0],
         s->element[i].H[0][0],s->element[i].H[1][0],s->element[i].H[2][0],
          s->element[i].oneovervolume,s->DOF );
    }
  }

  for(i=0;i<s->nelem;i++){
    if(i == itest1 ){ 
      for(k=0;k<4;k++)
         printf("Flux %d %d %6.4e %6.4e %6.4e %6.4e %6.4e %6.4e \n",i,k, 
	       s->element[i].H[0][k],
               s->element[i].H[1][k],
	       s->element[i].H[2][k], 
	       -s->element[i].HV[0][k],-s->element[i].HV[1][k],
	       -s->element[i].HV[2][k]); 
        printf("\n");
    }
  }

if(s->Re >0){
/* adding viscous terms to inviscid ones */
  for(i=0;i<s->nelem;i++)
    for(k=0;k<s->DOF;k++)
      for(j=0;j<4;j++)
        s->element[i].H[k][j]-=s->element[i].HV[k][j]; 
}

}

void Solver_Limiting(tSolver_State *s){
  uint i;
  double jumps[s->nelem];
  const uint k=0;   /* control of jumps of density */
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};
  int itest =-404;


  s->number_of_limiting=0;
  s->number_of_nolimiting=0;

  for(i=0;i<s->nelem;i++){ 
    jumps[i]=0;
    s->element[i].limited =0;
  }
  for(i=0;i<s->nface;i++){
    if(s->face[i].boundary<=0){  /* inside or periodic */
      double jump,local;
      uint _i,_j,_ii,_jj;

      _i = s->face[i].i;
      _ii= s->face[i].ii;
      _j = s->face[i].j;
      _jj= s->face[i].jj;
      local = s->element[_i].w[perm[_ii][0]][k] - 
                s->element[_j].w[perm[_jj][0]][k];
      jump = pow2(local);
      local = s->element[_i].w[perm[_ii][2]][k] -
              s->element[_i].w[perm[_ii][1]][k]+
              s->element[_j].w[perm[_jj][2]][k] -
              s->element[_j].w[perm[_jj][1]][k]; 
      jump+= pow2(local)/3;
      jump*= s->face[i].volume;
      jumps[_i] +=jump;
      jumps[_j] +=jump;
    }
  }

  for(i=0;i<s->nelem;i++){
    double indicator;

    indicator = jumps[i]/s->element[i].circumference*
                pow_fast(s->element[i].oneovervolume,0.75);

       if(i == itest)
           printf("Limiting: %d %7.5g %7.5g %7.5g %7.5g %7.5g \n",
               i,jumps[i],s->element[i].circumference,
               s->element[i].oneovervolume,
               pow_fast(s->element[i].oneovervolume,0.75),indicator); 

    if(indicator > s->limiting){
      int neigh,j,i1,j1;
      double gi,gj,wiL,wiR,wjL,wjR,aver;
      BOOL limit;
      uint count=0;

      if(i == -21){
	aver = (s->element[i].w[0][0]+s->element[i].w[1][0]+
		s->element[i].w[2][0])/3;
	printf("Lim: %d %4.2g %7.5g %7.2e %7.5g %7.2e %7.5g %7.2e %7.5g\n",
	       i,indicator,
	       s->element[i].w[0][0],s->element[i].w[0][0]-aver,
	       s->element[i].w[1][0],s->element[i].w[1][0]-aver,
	       s->element[i].w[2][0],s->element[i].w[2][0]-aver,
	       aver);
      }
      limit = TRUE;

      // in some cases we do not perform the limiting
      for(i1=0;i1<3;i1++){
	gi = 2*s->element[i].w[perm[i1][0]][0] -
	  s->element[i].w[perm[i1][1]][0] - s->element[i].w[perm[i1][2]][0];
	j = s->element[i].neigh[i1];

	if(j < 0) {
	  count ++;
	}
	else{
	  j1 = s->element[i].neighind[i1];
	  gj = 2*s->element[j].w[perm[j1][0]][0] -
	    s->element[j].w[perm[j1][1]][0]- s->element[j].w[perm[j1][2]][0];

	  if(gi*gj<0){
	    wiL = s->element[i].w[perm[i1][0]][0] +
	      s->element[i].w[perm[i1][1]][0]
	      - s->element[i].w[perm[i1][2]][0];

	    wiR = s->element[i].w[perm[i1][0]][0] +
	      s->element[i].w[perm[i1][2]][0]
	      - s->element[i].w[perm[i1][1]][0];

	    wjR = s->element[j].w[perm[j1][0]][0] +
	      s->element[j].w[perm[j1][1]][0]
	      - s->element[j].w[perm[j1][2]][0];

	    wjL = s->element[j].w[perm[j1][0]][0] +
	      s->element[j].w[perm[j1][2]][0]
	      - s->element[j].w[perm[j1][1]][0];

	    if(gi>0 && wiL<=wjL && wiR <=wjR ) count++;
	    if(gi<0 && wiL>=wjL && wiR >=wjR ) count++;
	  }
	}
      }

      if(count == 3){
	limit = FALSE;
	s->number_of_nolimiting++;
      }


      /* if(!limit){
	     	printf("L: %d %d  %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g\n",
	       i,j,
	       s->element[i].w[perm[i1][0]][0],
	       s->element[i].w[perm[i1][1]][0],
	       s->element[i].w[perm[i1][2]][0],
	       s->element[j].w[perm[j1][0]][0],
	       s->element[j].w[perm[j1][1]][0],
	       s->element[j].w[perm[j1][2]][0]); 
		printf("L: %d %d  %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g\n\n",
		i,j,gi,gj,wiL,wiR,wjL,wjR); 
		} 
      */

      if(limit){
	Solver_Limiting_Element(s,i,indicator);


	s->element[i].limited =1;
	/*      for(j=0;j<3;j++)
		if((neigh=s->element[i].neigh[j]) >=0){
		Solver_Limiting_Element(s,neigh);
		} */
	s->number_of_limiting++;

	/*	printf("\n %7.5g %7.5g %7.5g  %7.5g %7.5g  %7.5g \n",
	       s->element[i].ve->P[0]->Point[0],
	       s->element[i].ve->P[0]->Point[1],
	       s->element[i].w[0][0],s->element[i].w[0][1],
	       s->element[i].w[0][2],s->element[i].w[0][3]);
	printf("%7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",
	       s->element[i].ve->P[1]->Point[0],
	       s->element[i].ve->P[1]->Point[1],
	       s->element[i].w[1][0],s->element[i].w[1][1],
	       s->element[i].w[1][2],s->element[i].w[1][3]);
	printf("%7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",
	       s->element[i].ve->P[2]->Point[0],
	       s->element[i].ve->P[2]->Point[1],
	       s->element[i].w[2][0],s->element[i].w[2][1],
	       s->element[i].w[2][2],s->element[i].w[2][3]);
	*/
      }
    }
  }

}

void  Solver_Limiting_Element(tSolver_State *s, uint i, double indicator){
  uint j,k;
  double aver,alpha;

  //  alpha=(indicator-1)/(1+indicator);
  //alpha = 0.8;
  for(k=0;k<4;k++){
    aver = 0;
    for(j=0;j<s->DOF;j++)
      aver+= s->element[i].w[j][k];
    aver/=3;
    for(j=0;j<s->DOF;j++)
        s->element[i].w[j][k]=aver; 
    /*s->element[i].w[j][k]=alpha*aver + (1-alpha)*s->element[i].w[j][k];*/
  }

}

void Solver_Compute_Inviscid_Volumes(tSolver_State *s){
  uint i,j,k;
  int ic;
  const int itest = -883,itest0= -2399;
  double rho,p,e;
  double2 v,n;
  double VI[s->DOF][4];
  
  for(i=0;i<s->nelem;i++){
    double dertestfunction[s->DOF][2];
    /* derivatives of test functions after transformation 
       see page DGM-INT-16a   */
    for(j=0;j<3;j++)
      for(k=0;k<2;k++)
	dertestfunction[j][k] = s->element[i].dertest[j][k];

    if(i == itest || i == itest0){
      printf("%7.5g %7.5g triangle coordinates %d \n",
	 s->element[i].ve->P[0]->Point[0],s->element[i].ve->P[0]->Point[1],i);
      printf("%7.5g %7.5g\n",
         s->element[i].ve->P[1]->Point[0],s->element[i].ve->P[1]->Point[1]);
      printf("%7.5g %7.5g\n",
         s->element[i].ve->P[2]->Point[0],s->element[i].ve->P[2]->Point[1]);
   }

    if(i == itest || i == itest0){
       printf("Test functions %d %10.5g %10.5g %10.5g \n",i,
	  dertestfunction[0][0],dertestfunction[1][0],dertestfunction[2][0]);
       printf("Test functions %d %10.5g %10.5g %10.5g \n",i,
	   dertestfunction[0][1],dertestfunction[1][1],dertestfunction[2][1]);}

    for(j=0;j<4;j++)
      for(k=0;k<s->DOF;k++)
	VI[k][j] = 0;   /* clearing of arrays */
    
    for(j=0;j<s->DOF;j++){   /* j goes over integration nodes */
      rho=s->element[i].w[j][0];
      v[0]=s->element[i].w[j][1]/rho;
      v[1]=s->element[i].w[j][2]/rho;
      e=s->element[i].w[j][3];
      p=(s->kappa-1)*(e-0.5*rho*double2_norm2(v));
      
      ic=s->element[i].icve;
      if(ic >= 0){ /*  curved elements */
	for(k=0;k<s->DOF;k++){   /* k goes over test functions */
	  VI[k][0]+= rho*v[0]*s->curvedelement[ic].dertest[k][0][j] +
	    rho*v[1]*s->curvedelement[ic].dertest[k][1][j];
	  VI[k][1]+= (rho*v[0]*v[0]+p)*s->curvedelement[ic].dertest[k][0][j] +
	    rho*v[0]*v[1]*s->curvedelement[ic].dertest[k][1][j];
	  VI[k][2]+= rho*v[0]*v[1]*s->curvedelement[ic].dertest[k][0][j] +
	    (rho*v[1]*v[1]+p)*s->curvedelement[ic].dertest[k][1][j];
	  VI[k][3]+= (e+p)*v[0]*s->curvedelement[ic].dertest[k][0][j] +
	    (e+p)*v[1]*s->curvedelement[ic].dertest[k][1][j];
	  }

	for(k=0;k<s->DOF;k++){
          if((i == itest || i == itest0) && j == 0){
	  printf("Der test functions %d %d x %7.5g  %7.5g %7.5g %7.5g \n",i, k,
		 dertestfunction[k][0],s->curvedelement[ic].dertest[k][0][0],
		 s->curvedelement[ic].dertest[k][0][1],
		 s->curvedelement[ic].dertest[k][0][2] ); 
	  printf("Der test functions %d %d y %7.5g  %7.5g %7.5g %7.5g \n",i, k,
		 dertestfunction[k][1],s->curvedelement[ic].dertest[k][1][0],
		 s->curvedelement[ic].dertest[k][1][1],
		 s->curvedelement[ic].dertest[k][1][2] ); 
       	  }
        }
	if(i == itest || i == itest0) 
            for(k=0;k<s->DOF;k++)
              printf("VI[k] %d %d %d %6.4g %6.4g %6.4g %6.4g  \n",
	       i,j,k,VI[k][0],VI[k][1],VI[k][2],VI[k][3]);  

      }
      else{ /* non curved elements */
	for(k=0;k<s->DOF;k++){   /* k goes over test functions */
	  VI[k][0]+= rho*v[0]*dertestfunction[k][0] +
	    rho*v[1]*dertestfunction[k][1];
	  VI[k][1]+= (rho*v[0]*v[0]+p)*dertestfunction[k][0] +
	    rho*v[0]*v[1]*dertestfunction[k][1];
	  VI[k][2]+= rho*v[0]*v[1]*dertestfunction[k][0] +
	    (rho*v[1]*v[1]+p)*dertestfunction[k][1];
	  VI[k][3]+= (e+p)*v[0]*dertestfunction[k][0] +
	    (e+p)*v[1]*dertestfunction[k][1];
	if(i == itest || i == itest0) 
            printf("VI[k] %d %d %d %6.4g %6.4g %6.4g %6.4g  \n",
	       i,j,k,VI[k][0],VI[k][1],VI[k][2],VI[k][3]);  

	       
	
	} 
      }
    }
    for(k=0;k<s->DOF;k++)
      for(j=0;j<4;j++)
	s->element[i].H[k][j] = -VI[k][j]/s->DOF; 
    
    for(j=0;j<4;j++)
      if(i == itest || i == itest0) 
	printf("Volume Integr %d %d  %7.5g %7.5g %7.5g \n",
	       i,j,s->element[i].H[0][j],
	       s->element[i].H[1][j],s->element[i].H[2][j]);
    
      
  }
}
void  Write_Triangle(tSolver_State *s, int i){
uint j,k;
for(j=0;j<4;j++){
  k=j;
  if(j == 3) k=0;
  printf("%9.6g %9.6g triangle %d \n",
   s->element[i].ve->P[k]->Point[0],s->element[i].ve->P[k]->Point[1],i);
}
printf("\n");
}
void  Write_Edge(tSolver_State *s, int k){
uint j,i;
i = s->face[k].i;
j = s->face[k].ii;
printf("%9.6g %9.6g edge %d \n",
   s->element[i].ve->P[j]->Point[0],s->element[i].ve->P[j]->Point[1],k);
j++;
if(j ==3) j=0;
printf("%9.6g %9.6g edge %d \n\n",
   s->element[i].ve->P[j]->Point[0],s->element[i].ve->P[j]->Point[1],k);

}
double Edge_Integral(double aL, double aR, double bL, double bR){
  return (aL*(2*bL + bR ) + aR*(bL + 2*bR ))/6.;
  /* Simpson integration rule   */
}

double Temperature(double4 w){
return (w[3]-0.5*(w[1]*w[1]+w[2]*w[2])/w[0])/w[0];
}
double Pressure(double4 w,double kappa){
return (w[3]-0.5*(w[1]*w[1]+w[2]*w[2])/w[0])*(kappa-1);
}

void Solver_Inviscid_Stabilization(tSolver_State *s){
  int i,j,l,_i,_j;
  double4  wjumpL, wjumpR;
  double area;
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};

  for(i=0;i<s->nelem;i++){  /* computation of [phi] * [w] */
    for(j=0;j<s->DOF;j++){
      _i=s->element[i].neigh[j];
      if(_i>=0){   /*inner edge */
	_j=s->element[i].neighind[j];
	/* computation of [w] over edge j */
	Compute_Wjump(s->element[i].w,s->element[_i].w,j,_j,wjumpL,wjumpR);
	/* computation of [w] [phi] /Re */
	area = 1E-01;   /* 0.5 is the weight */
	for(l=0;l<4;l++){  /* for K_i test functions */
	  s->element[i].H[perm[j][0]][l]-=( wjumpL[l] + wjumpR[l])*area;
	  s->element[i].H[perm[j][1]][l]-=(-wjumpL[l] + wjumpR[l])*area;
	  s->element[i].H[perm[j][2]][l]-=( wjumpL[l] - wjumpR[l])*area;
	} 
	if( fabs( (wjumpL[l] + wjumpR[l]))*area >= 1E-01 ||
	    fabs( (wjumpL[l] - wjumpR[l]))*area >= 1E-01){
	  printf("JUM-L %d %d %d %7.5g %7.5g %7.5g %7.5g \n",
		 i,j,_i,wjumpL[0],wjumpL[1],wjumpL[2],wjumpL[3]);
	  printf("JUM-R          %7.5g %7.5g %7.5g %7.5g \n",
		 wjumpR[0],wjumpR[1],wjumpR[2],wjumpR[3]);
	}
      }  /* end of (_i>=0)  */
    }  /* end of j */
  }  /* end of i */
}

void Solver_Compute_Viscous_Edges(tSolver_State *s){
  int i,j;
  uint k,l;
  int ic;
  const int lastcomponent=4;   /* =4  or =3, in  HV[.][lastcomponent]  */
  const int itest = -883, itest0 = -883, itest1 = -58, itest2 = -2339;
  double4 Dleft, Dright,  wjumpL, wjumpR,DwnL,DwnR, wloc;
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};
  double Re,kappaRePr;
  float ibd;
  /*  printf("Step 1\n"); */
  Re = s->Re;
  kappaRePr = s->kappaRePr;
  
  for(k=0;k<s->nface;k++){  /* computation of <D(w)nabla w> * [phi] */
    int ii,jj;
    double area;
      if(s->face[k].boundary<=0) { /* inner edge */
      i = s->face[k].i;
      ii= s->face[k].ii;
      j = s->face[k].j;
      jj= s->face[k].jj;
      
      Compute_DwnablaW(s,i,ii,j,jj,s->face[k].n,s->face[k].n,
                       s->face[k].volume, s->face[k].volume, Dleft,Dright);

      /* computation for test functions having support on edge k  */
      /*area = s->element[i].oneovervolume*0.5; */  /* 0.5 is the weight */
      area = 0.5;   /* 0.5 is the weight */
      for(l=0;l<4;l++){ /* for K_i test functions */
	s->element[i].HV[perm[ii][0]][l]+=( Dleft[l] + Dright[l])*area;
	s->element[i].HV[perm[ii][1]][l]+=(-Dleft[l] + Dright[l])*area;
	s->element[i].HV[perm[ii][2]][l]+=( Dleft[l] - Dright[l])*area;
      }

      /*area = s->element[j].oneovervolume*0.5; */  /* 0.5 is the weight */
      area = 0.5;   /* 0.5 is the weight */
      for(l=0;l<4;l++){ /* for K_j test functions, n_j = -n_i */
	s->element[j].HV[perm[jj][0]][l]-=( Dleft[l] + Dright[l])*area;
	s->element[j].HV[perm[jj][1]][l]-=( Dleft[l] - Dright[l])*area;
	s->element[j].HV[perm[jj][2]][l]-=(-Dleft[l] + Dright[l])*area;
      }
    } /* end of inner edge */
    else if(s->boundary[s->face[k].boundary-1].type==
	    cSolver_BoundaryType_Inlet ||
	    s->boundary[s->face[k].boundary-1].type==
	    cSolver_BoundaryType_Slip){  /* Dirichlet boundary condition */

      i = s->face[k].i;
      ii= s->face[k].ii;

/*          Write_Triangle(s,i); */
/*      Write_Edge(s,k); */

      ic=s->element[i].icve;	
      if(ic >= 0){  /*curved elements  */
/*          printf("Normals: %d %d %d %8.2e %8.2e %8.2e %8.2e%8.2e %8.2e\n",
                 i,ii,ic,
		 s->face[k].n[0],s->face[k].n[1], 
		 s->curvedelement[ic].normal[2][0],
		 s->curvedelement[ic].normal[2][1],
		 s->curvedelement[ic].normal[3][0],
		 s->curvedelement[ic].normal[3][1]); */

        Compute_DwnablaW(s,i,ii,-1,0,
                         s->curvedelement[ic].normal[2],
                         s->curvedelement[ic].normal[3],
                         s->curvedelement[ic].length[2],
                         s->curvedelement[ic].length[3], Dleft,Dright);
      }
      else /*noncurved elements */
        Compute_DwnablaW(s,i,ii,-1,0,s->face[k].n,s->face[k].n,
                         s->face[k].volume,s->face[k].volume, Dleft,Dright);


      /* computation for test functions having support on edge k  */
      area = 0.5;   /* 0.5 is the weight */
/*      for(l=0;l<lastcomponent;l++){  */
      for(l=0;l<4;l++){  
	s->element[i].HV[perm[ii][0]][l]+=( Dleft[l] + Dright[l])*area;
	s->element[i].HV[perm[ii][1]][l]+=(-Dleft[l] + Dright[l])*area;
	s->element[i].HV[perm[ii][2]][l]+=( Dleft[l] - Dright[l])*area;
      }
    }   /* end of noundary edges */
    
  }  /* end  of  for(k=0;k<s->nface;k++) */

  if(itest >=0){
    printf("\n after <D(w)nabla w>[phi], D(w)nabla w* phi  -  <1>, (1)\n");
    Show_Viscous_Term(s,itest);
    }
  if(itest0 >=0){
    printf("\n after <D(w)nabla w>[phi], D(w)nabla w* phi  -  <1>, (1)\n");
    Show_Viscous_Term(s,itest0);
    }


    /* printf("Step 2\n"); */

  for(i=0;i<s->nelem;i++){  /* computation of <D(w)nablaphi> * [w] */
    int _i,_j,m,n,edge,ic;
    double scalprod,rv2,area,theta_wall;
    double4 wl,wr,wDl,wDr;
    double2 noutL,noutR;
    double3 theta,detJF,pressure, density;
    double2 local_nablarho[3];
    double local_nablarho1[3][2];
    //    lintransf local_dertest1[3];
    bilintransf local_dertest;
    double Dwnablaphix[s->DOF][s->DOF][4][4],Dwnablaphiy[s->DOF][s->DOF][4][4];
    /*                     1. integ. node, 2. phi, 3. compon. phi*/

    /* computation of <D(w)nablaphi> in all integration nodes */
    area = 0.5/Re;   /* stabilization term, sigma*0.5 is the weight */
    //        area = 0.5/Re*16; // A3
    area = 0.5/Re*1;

    theta_wall = s->thetainfty*(1+(s->kappa-1)/2*s->Minfty*s->Minfty);

    ic=s->element[i].icve;
    for(j=0;j<s->DOF;j++){
      for(m=0;m<s->DOF;m++){
        for(n=0;n<2;n++){
          if(ic>=0) /* curved element */
	    local_dertest[m][n][j]=s->curvedelement[ic].dertest[m][n][j];
          else
	    local_dertest[m][n][j]=s->element[i].dertest[m][n];
	}
      }
      for(n=0;n<2;n++){
      	if(ic>=0){ /* curved element */
      	  local_nablarho[j][n]=s->curvedelement[ic].nablarho[j][n];
      	  local_nablarho1[j][n]=s->curvedelement[ic].nablarho[j][n];
	}
	else{
	  local_nablarho[j][n]=s->element[i].nablarho[n];
	  local_nablarho1[j][n]=s->element[i].nablarho[n];
	}
      }
      if(ic>=0)   /* curved element */
	detJF[j]=2/s->curvedelement[ic].detJF[j];
      else
        detJF[j]=s->element[i].oneovervolume;
      
    }

    for(j=0;j<3;j++)
      for(n=0;n<2;n++){
	if(fabs(local_nablarho[j][n] - local_nablarho1[j][n]) > 1E-06)
	  printf("INCOMPATIBALE VALUES!!! \n\n\n\n");
      }
    
    Compute_Dw_nablaPhi(s->DOF,s->element[i].w,local_nablarho,
                        detJF,Re, kappaRePr,
                        local_dertest, Dwnablaphix, Dwnablaphiy);
    if(i == itest1|| i == itest2){
      for(j=0;j<3;j++)
        printf("w   (%6.3g,%6.3g,%6.3g,%6.3g) %d \n", 
	       s->element[i].w[j][0],s->element[i].w[j][1],
	       s->element[i].w[j][2],s->element[i].w[j][3],ic);
      for(j=0;j<3;j++)
	printf("&phi (%6.3g,%6.3g), (%6.3g,%6.3g), (%6.3g,%6.3g) \n", 
	       local_dertest[0][0][j],local_dertest[0][1][j],
	       local_dertest[1][0][j],local_dertest[1][1][j],
	       local_dertest[2][0][j],local_dertest[2][1][j]);
	
      for(j=0;j<3;j++)
        for(k=0;k<3;k++)
	  for(m=0;m<4;m++){
	    printf("Dw&phi[%d][%d][%d]=(%6.3g,%6.3g,%6.3g,%6.3g),(%6.3g,%6.3g,%6.3g,%6.3g)\n",
		   j,k,m,Dwnablaphix[j][k][m][0],Dwnablaphix[j][k][m][1],
		   Dwnablaphix[j][k][m][2],Dwnablaphix[j][k][m][3],
		   Dwnablaphiy[j][k][m][0],Dwnablaphiy[j][k][m][1],
		   Dwnablaphiy[j][k][m][2],Dwnablaphiy[j][k][m][3] );
	  }  
      printf("-------------------------------------------------\n");
    } 

    for(j=0;j<3;j++){ /* j goes over edges  */
      _i=s->element[i].neigh[j];
      ic=s->element[i].icve;
      if(i == itest ||i == itest0 )
	printf("EDGE %d %d %d %d  \n", i,j,_i,ic);

      if(_i < 0 && ic>=0){ /* curved element */
        noutL[0]= s->curvedelement[ic].normal[2][0]
                 *s->curvedelement[ic].length[2];
        noutL[1]= s->curvedelement[ic].normal[2][1]
                 *s->curvedelement[ic].length[2];
        noutR[0]= s->curvedelement[ic].normal[3][0]
                 *s->curvedelement[ic].length[3];
        noutR[1]= s->curvedelement[ic].normal[3][1]
                 *s->curvedelement[ic].length[3];
      }
      else{  /*noncurved element */
        noutL[0]=s->element[i].dertest[j][0];  /* outer normal * length */
        noutL[1]=s->element[i].dertest[j][1];
        noutR[0]=s->element[i].dertest[j][0];  
        noutR[1]=s->element[i].dertest[j][1];
      }
      if(_i>=0){   /*inner edge */
	_j=s->element[i].neighind[j];
	/* computation of [w] over edge j */
	Compute_Wjump(s->element[i].w,s->element[_i].w,j,_j,wjumpL,wjumpR);
        /* computation of [w] [phi] /Re */
/*        for(l=0;l<4;l++){ !!!!!!*/ /* for K_i test functions */
	for(l=0;l<4;l++){  /* for K_i test functions */
	  s->element[i].HV[perm[j][0]][l]-=( wjumpL[l] + wjumpR[l])*area;
          s->element[i].HV[perm[j][1]][l]-=(-wjumpL[l] + wjumpR[l])*area;
          s->element[i].HV[perm[j][2]][l]-=( wjumpL[l] - wjumpR[l])*area;
	} 

/*        for(l=0;l<4;l++){ 
	  s->element[_i].HV[perm[_j][0]][l]+=( wjumpL[l] + wjumpR[l])*area;
          s->element[_i].HV[perm[_j][1]][l]+=( wjumpL[l] - wjumpR[l])*area;
	  s->element[_i].HV[perm[_j][2]][l]+=(-wjumpL[l] + wjumpR[l])*area;
        } */
        if(i == itest ||i == itest0 ){
	    printf("\n after  [w][phi]  -  <3> , (elemts: %d,%d) \n",
	        i,_i);
          Show_Viscous_Term(s,i);
        }
        ibd = 2;
       }
      else{  /* boundary edge with the index "edge" */
        edge = s->element[i].iface;
        if(s->boundary[s->face[edge].boundary-1].type==
	    cSolver_BoundaryType_Inlet ||
          s->boundary[s->face[edge].boundary-1].type==
	    cSolver_BoundaryType_Slip){  /* Dirichlet boundary condition */

      /*    setting of boundary state wD  */
      /* computation of absolute temperature in two others integration nodes*/
          for (l=0;l<3;l++){
	    theta[l] = Temperature(s->element[i].w[perm[j][l]]);
	    pressure[l] = Pressure(s->element[i].w[perm[j][l]],s->kappa);
	    density[l] = s->element[i].w[perm[j][l]][0];
	  }
	  if(i == itest || i == itest0) 
	    printf("Temperatures %d %d %7.5g %7.5g \n",i,j,theta[0],theta[1]); 

          if(s->boundary[s->face[edge].boundary-1].type==
	     cSolver_BoundaryType_Slip ){ 
       /* NSexact   */
           /* && s->face[edge].boundary != 4  ){ */

            wDl[0]=s->element[i].w[perm[j][0]][0]
                  +s->element[i].w[perm[j][2]][0]
                  -s->element[i].w[perm[j][1]][0]; 
            wDr[0]=s->element[i].w[perm[j][0]][0]
                  +s->element[i].w[perm[j][1]][0]
                  -s->element[i].w[perm[j][2]][0]; 
            for(n=0;n<2;n++){
	      wDl[n+1] = 0;     /* noslip condition */
              wDr[n+1] = 0;
            }
	    /* extrapolation of temperature */
	    /* adiabatic wall */
	    //wDl[3] = wDl[0]*(theta[0]+theta[2]-theta[1]); 
	    //wDr[3] = wDr[0]*(theta[0]+theta[1]-theta[2]); 
	    /*	    printf("Temp: %7.5g %7.5g %7.5g %7.5g \n",
		    theta[0],theta[1],theta[2],theta_wall); */
	    /* prescribed temperature */
	    wDl[3] = wDl[0]*theta_wall;
	    wDr[3] = wDr[0]*theta_wall; 

	    /* extrapolation of pressure */
	    /* wDl[3] = (pressure[0]+pressure[2]-pressure[1])/(s->kappa-1); 
	       wDr[3] = (pressure[0]+pressure[1]-pressure[2])/(s->kappa-1); */
	    /* extrapolation of density */
	     /*   wDl[3] = (pressure[0]+pressure[2]-pressure[1])/(s->kappa-1); 
	       wDr[3] = (pressure[0]+pressure[1]-pressure[2])/(s->kappa-1);
	       wDl[0] = wDl[3]/(theta[0]+theta[2]-theta[1]); 
	       wDr[0] = wDr[3]/(theta[0]+theta[1]-theta[2]);  */

	    /* extrapolation of energy */
	    /*            wDl[3]=s->element[i].w[perm[j][0]][3]
                  +s->element[i].w[perm[j][2]][3]
                  -s->element[i].w[perm[j][1]][3]; 
            wDr[3]=s->element[i].w[perm[j][0]][3]
                  +s->element[i].w[perm[j][1]][3]
                  -s->element[i].w[perm[j][2]][3];  */
         }
         else{ /* inlet  */
          /* NSexact */
         /* if( s->face[edge].boundary != 4){  /*  */
             for(l=0;l<3;l++){
	       wDl[l]=s->boundary[s->face[edge].boundary-1].w[l];
               wDr[l]=s->boundary[s->face[edge].boundary-1].w[l];
             }
            /* NSexact */
            /*   NSexact(s->element[i].ve->P[perm[j][0]]->Point[0],
                            s->element[i].ve->P[perm[j][0]]->Point[1],wDl);
               NSexact(s->element[i].ve->P[perm[j][1]]->Point[0],
                            s->element[i].ve->P[perm[j][1]]->Point[1],wDr);
            }*/
         /* NSexact */
         /* else{
             wDl[0]=1;
             wDr[0]=1;
             wDl[1]=1*wDl[0];
             wDr[1]=1*wDr[0];
             wDl[2]=0;
             wDr[2]=0;
           }  */
           rv2 = 0.5*(wDl[1]*wDl[1] + wDl[2]*wDl[2])/wDl[0];
            wDl[3] = wDl[0]*(theta[0]+theta[2]-theta[1])+rv2; 
            wDr[3] = wDr[0]*(theta[0]+theta[1]-theta[2])+rv2; 
/*           wDl[3] = wDl[0]*(3*theta[1]-theta[0])/2 + rv2;
           wDr[3] = wDr[0]*(3*theta[0]-theta[1])/2 + rv2; */
/*           wDl[3] = wDl[0]*(theta[1]+theta[0])/2 + rv2;
           wDr[3] = wDr[0]*(theta[0]+theta[1])/2 + rv2;*/
	 }
        if(i == itest || i == itest0){
             printf("wDl: %7.5g %7.5g %7.5g %7.5g \n",
                    wDl[0],wDl[1],wDl[2],wDl[3]);
             printf("wDr: %7.5g %7.5g %7.5g %7.5g \n",
                    wDr[0],wDr[1],wDr[2],wDr[3]);
       }
       /* computation of   w-wD over boundary edge j  */
        double4_vertex(wl, s->element[i].w[perm[j][0]], 
		    s->element[i].w[perm[j][2]],s->element[i].w[perm[j][1]]);
        double4_vertex(wr, s->element[i].w[perm[j][0]], 
		    s->element[i].w[perm[j][1]],s->element[i].w[perm[j][2]]);
/*        if(i == itest || i == itest0){
             printf("wl %d %7.5g %7.5g %7.5g %7.5g \n",0,
             -1,s->element[i].w[perm[j][0]][0], 
             s->element[i].w[perm[j][2]][0],s->element[i].w[perm[j][1]][0]);
             printf("wl %d %7.5g %7.5g %7.5g %7.5g \n",1,
             wl[0],s->element[i].w[perm[j][0]][1], 
             s->element[i].w[perm[j][2]][1],s->element[i].w[perm[j][1]][1]);
             printf("wl %d %7.5g %7.5g %7.5g %7.5g \n",2,
             wl[1],s->element[i].w[perm[j][0]][2], 
             s->element[i].w[perm[j][2]][2],s->element[i].w[perm[j][1]][2]);
             printf("wl %d %7.5g %7.5g %7.5g %7.5g \n--------\n",3,
             wl[2],s->element[i].w[perm[j][0]][3], 
             s->element[i].w[perm[j][2]][3],s->element[i].w[perm[j][1]][3]);
        } */

        for(n=0;n<4;n++){
          wjumpL[n]=wl[n]-wDl[n]; 
          wjumpR[n]=wr[n]-wDr[n]; 
        }
       /* computation of w phi */
       /* computation for test functions having support on edge k  */
       /* for(l=0;l<lastcomponent;l++){!!!!!! */ /* for K_i test functions */
        for(l=0;l<lastcomponent;l++){ /* for K_i test functions */
	  s->element[i].HV[perm[j][0]][l]-=( wjumpL[l] + wjumpR[l])*area;
          s->element[i].HV[perm[j][1]][l]-=(-wjumpL[l] + wjumpR[l])*area;
          s->element[i].HV[perm[j][2]][l]-=( wjumpL[l] - wjumpR[l])*area;
        }
	//	 if(i == itest || i == itest0){
	/*	if( s->boundary[s->face[edge].boundary-1].type==
	    cSolver_BoundaryType_Slip) {
	    printf("wl,    wR     : %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",
                  wl[0],wl[1],wl[2],wl[3],wr[0],wr[1],wr[2],wr[3]);
	  printf("wDl,   wDR    : %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",
	  wDl[0],wDl[1],wDl[2],wDl[3],wDr[0],wDr[1],wDr[2],wDr[3]); 
	  printf("wjumpL, wjumpR: %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",
                  wjumpL[0],wjumpL[1],wjumpL[2],wjumpL[3],wjumpR[0],wjumpR[1],wjumpR[2], wjumpR[3]);
		  } */
	/*  printf("\n after  w*phi boundary - (4)+(5) area=%7.5g\n",area);
	  Show_Viscous_Term(s,i);
	  } */ 

        ibd = 1;
       }  /* end of inlet or outlet */
       else{
          ibd = 0;
         break; /* outlet */
       }
      } /* end of inner or boundary */
      if(i == itest || i == itest0){ 
         printf("ibd=%g\n",ibd);
         printf(".terms <2> and (2)+(3) %d %d %d:",i,j,k);
      }

      if(ibd == 0) printf("BAD CYCLING!!!!! \n");
      for(k=0;k<s->DOF;k++){  /* k goes over test function */
        for(l=0;l<4;l++){  
	  /*              l goes over component of test function phi*/
	  /* computation of <Dw nabla phi> over edge j */
	  for(m=0;m<4;m++){
	    DwnL[m] = Dwnablaphix[perm[j][0]][k][l][m]*noutL[0] +
	      Dwnablaphiy[perm[j][0]][k][l][m]*noutL[1];
	    DwnR[m] = Dwnablaphix[perm[j][1]][k][l][m]*noutR[0] +
	      Dwnablaphiy[perm[j][1]][k][l][m]*noutR[1];
	  }
	  scalprod=0;

/*          if((itest == i  || i == itest0) && k==1 && l==0){
            printf("\n normals:%7.5g %7.5g %7.5g %7.5g\n",
              noutL[0],noutL[1],noutR[0],noutR[1]);
            for(m=0;m<4;m++) printf("SS %d %d %d %d %7.5g %7.5g %7.5g %7.5g \n",
               i,k,l,m,DwnL[m],DwnR[m],wjumpL[m], wjumpR[m]);
          }  */
          
	  for(m=0;m<4;m++)  
	    /* computing of int_{\Gamma} <Dw nabla phi> \CDOT [w] dS */
	     scalprod+= Edge_Integral(DwnL[m],DwnR[m],wjumpL[m], wjumpR[m]);
          if(!(l > lastcomponent-1 && ibd == 1))
	    /* NONSYMMETRIC */
	    s->element[i].HV[k][l]-=scalprod/ibd; /*ibd inner:<.>=(.+.)/2*/
	    /* SYMMETRIC */
	    /*s->element[i].HV[k][l]+=scalprod/ibd;*/ /* ibd inner: <.>=(.+.)/2 */

	  if(i == itest || i == itest0) printf("%10.3e",-scalprod/ibd);
        }   /* end of l */

	if(itest == i  || i == itest0) 
	  if(k == 2) printf("\n-------------------------------\n");
	  else printf("%7.4g %7.4g\n terms <2> and (2)+(3) %d %d %d:",
		      s->element[i].dertest[j][0],s->element[i].dertest[j][1],
		      i,j,k+1);
      
       if( i == itest0) Show_Viscous_Term(s,itest0);
       if( i == itest) Show_Viscous_Term(s,itest);
      }
      if(i == itest ){
	printf("\n PART OF  <D(w)nabla phi>[w]  -  <2> or (2)\n");
	Show_Viscous_Term(s,itest);
      }
      if(i == itest0 ){
	printf("\n PART OF  <D(w)nabla phi>[w]  -  <2> or (2)\n");
	Show_Viscous_Term(s,itest0);
      }
    } /* end of j */
  } /* end of i */

  if(itest0 >=0){
    printf("$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n");
    printf("\n after <D(w)nabla phi>[w]  -  <2> or (2)\n");
    Show_Viscous_Term(s,itest0);
    }
  if(itest >=0){
    printf("\n after <D(w)nabla phi>[w]  -  <2> or (2)\n");
    Show_Viscous_Term(s,itest);
    }

  /* printf("Step 3\n"); */


  /* division by the area of element, from the time derivative term */
  /* multiplication by s->DOF, the weights from the time derivative term */
  for(i=0;i<s->nelem;i++){
    int ic;
    ic=s->element[i].icve;  
    if(ic >= 0){  /* curved elements */
      for(k=0;k<s->DOF;k++)
	for(j=0;j<4;j++)
	  s->element[i].HV[k][j]*=2/s->curvedelement[ic].detJF[k]*s->DOF; 
         /*   printf("Volumes %d %7.5g  %7.5g %7.5g %7.5g \n",i,
	     s->element[i].oneovervolume,2/s->curvedelement[ic].detJF[0],
         2/s->curvedelement[ic].detJF[1],2/s->curvedelement[ic].detJF[2]); */
      if(i == itest || i == itest0) 
	/*      if( i <2005 && i > 1995) */
	printf("Fluxes dividing1 %d %7.5g %7.5g %7.5g \n",
           i,s->element[i].HV[0][0],s->element[i].HV[1][0],
             s->element[i].HV[2][0]);
    }
    else{
      for(k=0;k<s->DOF;k++) /* non curved elements */
	for(j=0;j<4;j++)
	  s->element[i].HV[k][j]*=s->element[i].oneovervolume*s->DOF; 
      if(i == itest || i == itest0)
	/*  if( i <2005 && i > 1995) */
          printf("Fluxes dividing2  %d %7.5g %7.5g %7.5g %7.5g %d\n",
          i,s->element[i].HV[0][0],s->element[i].HV[1][0],
          s->element[i].HV[2][0], s->element[i].oneovervolume,s->DOF );
    }
  }

 

}

void NSexact(double x1, double x2, double4 w){
double p, theta0, v1, theta, rho, kappa, e;

kappa = 1.4;

      p=1.6;
      theta0=4.257142857;
      v1=x2;
      theta=theta0-0.257142857*x2*x2;/*   ^^^^^^^^^^^ = 0.5*Pr/kappa   */
      rho=p/theta/(kappa-1);
      e=p/(kappa-1)+0.5*rho*v1*v1;
      w[0]=rho;
      w[1]=rho*v1;
      w[2]=0;
      w[3]=e;

/*  printf("NSEXACT: %7.5g  %7.5g %7.5g  %7.5g  %7.5g  %7.5g \n",x1,x2,
          w[0],w[1],w[2],w[3]);*/
}

void Show_Viscous_Term(tSolver_State*s, uint i){
  uint j,k;
  printf("Viscous term, element %d \n",i);
  printf("           HV[0][k]     HV[1][k]     HV[2][k]\n"); 
  for(k=0;k<4;k++)
    printf("HV[j][%d]: %9.5g  %9.5g  %9.5g\n",k,
	   s->element[i].HV[0][k],s->element[i].HV[1][k],
	   s->element[i].HV[2][k]);

}
void Compute_Dw_nablaPhi(int DOF,   statevec w, double2 nablarho[3],
			 double3 oneoverarea, double Re,  double kappaRePr, 
			 bilintransf dertestfunction,
			 double Dwnablaphix[3][3][4][4],
                         double Dwnablaphiy[3][3][4][4]){ 
/*                              originaly [DOF][DOF][4][4]  */
  uint j,k,l,m;
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};
  double4 nablaphix, nablaphiy, wlocx[3], wlocy[3];
  double4 phi;

  /* computation of <D(w)\nabla phi>, supp(phi)=K_i ==> only one K_i */
  phi[0] = 0; /* always zero */
  for(k=0;k<DOF;k++){  /* k goes over test function */
    for(l=0;l<4;l++){  /* l goes over component of test function */
      for(m=0;m<4;m++){
	nablaphix[m]=0;
	nablaphiy[m]=0;
        phi[m] = 0;
      }
      for(j=0;j<3;j++){   /* j goes over edges */
        nablaphix[l] = dertestfunction[k][0][j]*oneoverarea[j];
        nablaphiy[l] = dertestfunction[k][1][j]*oneoverarea[j];
	/* dertestfunction = outer normal times length of edge */
        if(k==j) phi[l] = 1;  /* only this component is nonzero */
	ViscousTerm(w[j],phi,nablarho[j],nablaphix,nablaphiy,Re,kappaRePr,
		    wlocx[j],wlocy[j]);
      }
      for(j=0;j<3;j++){ 
	double4_vertex(Dwnablaphix[j][k][l],wlocx[perm[j][0]],
		       wlocx[perm[2][0]],wlocx[perm[j][1]]);
	double4_vertex(Dwnablaphiy[j][k][l],wlocy[perm[j][0]],
		       wlocy[perm[2][0]],wlocy[perm[j][1]]);
      }
	/*	if(k==1 && l==2 && j==2) 
	  printf("ViscousTerm: %6.4g %6.4g %6.4g %6.4g %6.4g %6.4g %6.4g \n",
	  w[j][0],w[j][1],w[j][2],w[j][3],wloc1[0],wloc1[1],wloc1[2]); */
	/*	if(k==1 && l==1 && j==2) wloc[0] = -10; */
	/*	if(k==1 && l==1 && j==2){ 
	  printf("nabla phi    : %6.4g %6.4g %6.4g %6.4g %6.4g %6.4g\n",
		 nablaphix[0],nablaphix[1],nablaphix[2],
		 nablaphiy[0],nablaphiy[1],nablaphiy[2]);
	  printf("ViscousTerm: %6.4g %6.4g %6.4g %6.4g %6.4g %6.4g %6.4g \n",
		 w[j][0],w[j][1],w[j][2],w[j][3],wloc2[0],wloc2[1],wloc2[2]);
		 } */
    }   /* end  of l  */
  }   /* end  of k  */
}

void Compute_Wjump(statevec wi, statevec wj,uint i, uint j,
		   double4 wjumpL, double4 wjumpR){
  uint l;
  double4 wloc;
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};

  /* element K_i  */
  double4_vertex(wjumpL, wi[perm[i][0]], wi[perm[i][2]], wi[perm[i][1]]);
  double4_vertex(wjumpR, wi[perm[i][0]], wi[perm[i][1]], wi[perm[i][2]]);
  /* element K_j, neighboor of K_i  */
  double4_vertex(wloc, wj[perm[j][0]], wj[perm[j][1]], wj[perm[j][2]]);
  for(l=0;l<4;l++) wjumpL[l]-=wloc[l];
  double4_vertex(wloc, wj[perm[j][0]], wj[perm[j][2]], wj[perm[j][1]]);
  for(l=0;l<4;l++) wjumpR[l]-=wloc[l];
}


void Compute_DwnablaW(tSolver_State*s, int i, uint ii, int j, uint jj, 
                      double2 nL, double2 nR, double edgelengthL, 
                      double edgelengthR, double4 Dleft, double4 Dright){
  uint l;
  /*  double3 Dxleft,Dyleft, Dxright, Dyright; */
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};
  int itest = -355, itest0= -118;

  if(i == itest || j == itest || i == itest0 || j == itest0){
    for(l=0;l<4;l++){
      printf("1 Dsxnablaw %d %d %d %7.5g %7.5g %7.5g \n", i,j,l,
	     s->element[i].Dsxnablaw[perm[ii][0]][l],
	     s->element[i].Dsxnablaw[perm[ii][1]][l],
	     s->element[i].Dsxnablaw[perm[ii][2]][l]);
    }
    for(l=0;l<4;l++){
      printf("1 Dsynablaw %d %d %d %7.5g %7.5g %7.5g \n", i,j,l,
	     s->element[i].Dsynablaw[perm[ii][0]][l],
	     s->element[i].Dsynablaw[perm[ii][1]][l],
	     s->element[i].Dsynablaw[perm[ii][2]][l]);
    }
  }

  /*  triangle K_i  */

  for(l=0;l<4;l++){
    Dleft[l] = s->element[i].Dsxnablaw[perm[ii][0]][l]*nL[0] +
      s->element[i].Dsynablaw[perm[ii][0]][l]*nL[1];
    Dright[l] = s->element[i].Dsxnablaw[perm[ii][1]][l]*nR[0] +
      s->element[i].Dsynablaw[perm[ii][1]][l]*nR[1];
  }

  if(i == itest || j == itest || i == itest0 || j == itest0){
    printf("***  Dleft %d %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n", 
	   i,Dleft[0],Dleft[1],Dleft[2],Dleft[3],nL[0],nR[0]);
    printf("***  Dright %d %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n", 
	   i,Dright[0],Dright[1],Dright[2],Dright[3],nL[1],nR[1]);
  }

  if(j>=0){  /* inner edge, triangle K_j  */
    if(i == itest || j == itest || i == itest0 || j == itest0){
      for(l=0;l<4;l++){
	printf("2 Dsxnablaw %d %d %d %7.5g %7.5g %7.5g \n", i,j,l,
	       s->element[j].Dsxnablaw[perm[jj][0]][l],
	       s->element[j].Dsxnablaw[perm[jj][1]][l],
	       s->element[j].Dsxnablaw[perm[jj][2]][l]);
      }
      for(l=0;l<4;l++){
	printf("2 Dsynablaw %d %d %d %7.5g %7.5g %7.5g \n", i,j,l,
	       s->element[j].Dsynablaw[perm[jj][0]][l],
	       s->element[j].Dsynablaw[perm[jj][1]][l],
	       s->element[j].Dsynablaw[perm[jj][2]][l]);
      }
    }
    for(l=0;l<4;l++){
      Dleft[l] += s->element[j].Dsxnablaw[perm[jj][1]][l]*nL[0] +
	s->element[j].Dsynablaw[perm[jj][1]][l]*nL[1];
      Dright[l] += s->element[j].Dsxnablaw[perm[jj][0]][l]*nR[0] +
	s->element[j].Dsynablaw[perm[jj][0]][l]*nR[1];
    }
    for(l=0;l<4;l++){
      Dleft[l]*=0.5;
      Dright[l]*=0.5;
    }

  }   /* end of boundary edge  */
  for(l=0;l<4;l++){
    Dleft[l]*=edgelengthL;
    Dright[l]*=edgelengthR;
    }
  if(i == itest || j == itest || i == itest0 || j == itest0){
    printf("***  Dleft %d %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n", 
	   i,Dleft[0],Dleft[1],Dleft[2],Dleft[3],nL[0],nR[0]);
    printf("***  Dright %d %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n", 
	   i,Dright[0],Dright[1],Dright[2],Dright[3],nL[1],nR[1]);
  }

}

void Solver_Compute_Viscous_Volumes(tSolver_State *s){
  uint i,j,k,l;
  int ic;
  const int itest = -19927, itest0= -19927;
  double VI[s->DOF][4];  
  double4 wx,wy,Dx[3],Dy[3];
  double wtest[3][4];
  double dertestfunction[s->DOF][2];
  double2 nablarho;
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};

  for(i=0;i<s->nelem;i++){
    double dertestfunction[s->DOF][2];
    /* derivatives of test functions after transformation 
       see page DGM-INT-16a   */
    for(j=0;j<3;j++)
      for(k=0;k<2;k++)
	dertestfunction[j][k] = s->element[i].dertest[j][k];

    for(j=0;j<4;j++)
      for(k=0;k<s->DOF;k++)
	VI[k][j] = 0;   /* clearing of arrays */
    
    ic=s->element[i].icve;
    if(ic < 0) {  /* noncurved elements, \nabla w is constant over K_i */
      /* computing of \partial w/x, w/y   */
      for(j=0;j<4;j++)
         wx[j]= wy[j]=0;
      for(j=0;j<s->DOF;j++){
        for(k=0;k<4;k++){
	  wx[k]+=dertestfunction[j][0]*s->element[i].w[j][k];
          wy[k]+=dertestfunction[j][1]*s->element[i].w[j][k];
        }
      }
      for(k=0;k<4;k++){  /* dertestfunction is in fact \partial\phi *area  */
        wx[k]*=s->element[i].oneovervolume;
        wy[k]*=s->element[i].oneovervolume;
      }
    /*  storing the value for edge integral terms */
     s->element[i].nablarho[0] = wx[0]; 
     s->element[i].nablarho[1] = wy[0]; 

    }  /* end of noncurved elements */
    /* computing of the term Ds(w)\nabla w */
    for(j=0;j<s->DOF;j++){   /* j goes over integration nodes */
      if(ic >= 0){ /*  curved elements another values of wx,wy,nablarho*/
      for(l=0;l<4;l++)
       wx[l]= wy[l]=0;
      for(l=0;l<s->DOF;l++){
          for(k=0;k<4;k++){
            wx[k]+=s->curvedelement[ic].dertest[l][0][j]*
                   s->element[i].w[l][k];
            wy[k]+=s->curvedelement[ic].dertest[l][1][j]*
                   s->element[i].w[l][k];
          }
        }
        for(k=0;k<4;k++){  /* dertestfunction is in fact \partial\phi *area  */
          wx[k]*=2/s->curvedelement[ic].detJF[j];
          wy[k]*=2/s->curvedelement[ic].detJF[j];
        }
    /*  storing the value for edge integral terms */
	s->curvedelement[ic].nablarho[j][0] = wx[0]; 
	s->curvedelement[ic].nablarho[j][1] = wy[0]; 
      }
      nablarho[0] = wx[0];
      nablarho[1] = wy[0];
      ViscousTerm(s->element[i].w[j],s->element[i].w[j],
                  nablarho,wx,wy,s->Re,s->kappaRePr, Dx[j],Dy[j]);
          if(i == itest ){
            for(l= 0;l<4;l++)
               printf("Check %d %d %d %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",i,j,l,
               s->element[i].w[j][0],s->element[i].w[j][1],
               s->element[i].w[j][2],s->element[i].w[j][3],
               wx[1],wx[2],Dx[j][l],Dy[j][l]);
          } 
    }
    for(j=0;j<s->DOF;j++){   /* j goes over integration nodes */
      ic=s->element[i].icve;
      if(ic >= 0){ /*  curved elements */
	for(k=0;k<s->DOF;k++)   /* k goes over test functions */
          for(l=0;l<4;l++)
	    VI[k][l]+=Dx[j][l]*s->curvedelement[ic].dertest[k][0][j]+
	      Dy[j][l]*s->curvedelement[ic].dertest[k][1][j];
      }
      else{ /* non curved elements */
	for(k=0;k<s->DOF;k++){   /* k goes over test functions */
	  for(l=0;l<4;l++)
	    VI[k][l]+=Dx[j][l]*dertestfunction[k][0] + 
	      Dy[j][l]*dertestfunction[k][1];
/*          if(i == itest ){
            l= 3;
            printf("Volume %d funct=%d %7.5g %7.5g %7.5g %7.5g %7.5g \n",j,k,
               VI[k][l],Dx[j][l],dertestfunction[k][0],
               Dy[j][l],dertestfunction[k][1]); 
          } */
        }
      }
    } /* end of j goes over integration nodes  */

    /* storing values for  boundary integral terms in VERTEXES of triangles*/
    for(j=0;j<s->DOF;j++){  
      double4_vertex(s->element[i].Dsxnablaw[j],Dx[perm[j][0]],
		     Dx[perm[j][2]],Dx[perm[j][1]]);
      double4_vertex(s->element[i].Dsynablaw[j],Dy[perm[j][0]],
		     Dy[perm[j][2]],Dy[perm[j][1]]);
    }
      if(i == itest)
	for(k=0;k<4;k++)
	  printf("##########Dsxnablaw: %7.5g  %7.5g %7.5g %7.5g %7.5g %7.5g\n",
		 Dx[0][k],Dx[1][k],Dx[2][k],s->element[i].Dsxnablaw[0][k],
		s->element[i].Dsxnablaw[1][k],s->element[i].Dsxnablaw[2][k]);
    /* computing of volume integral int_{K_i} D_s(w,nabla w) \phi_s dx */
    for(k=0;k<s->DOF;k++)
      for(j=0;j<4;j++)
	s->element[i].HV[k][j] = -VI[k][j]/s->DOF;
    
    if(i == itest || i == itest0) Show_Viscous_Term(s,i);
  }   /* end of i */

}

void Show_Inviscid_Term(tSolver_State*s, uint i){
  uint j,k;
     for(k=0;k<4;k++)
       printf("w[ ] %d %d %6.4e %6.4e %6.4e H[ ] %d %d  %6.4e %6.4e %6.4e\n",
            i,k, s->element[i].w[0][k],s->element[i].w[1][k],
		 s->element[i].w[2][k],i,k, 
		 s->element[i].H[0][k],s->element[i].H[1][k],
		 s->element[i].H[2][k]); 
        printf("\n");
  }
void Show_Inviscid_Term1(tSolver_State*s, uint i){
  uint j,k;
     for(k=0;k<4;k++)
       printf("w[ ] %d %d %6.4e %6.4e %6.4e K1[ ] %d %d  %6.4e %6.4e %6.4e\n",
            i,k, s->element[i].w[0][k],s->element[i].w[1][k],
		 s->element[i].w[2][k],i,k, 
		 s->element[i].K1[0][k],s->element[i].K1[1][k],
		 s->element[i].K1[2][k]); 
        printf(" --K1--\n");
  }


void Solver_Integrate_Euler(tSolver_State *s){
  int const itest= -100;
  if(s->Re>0){
    Solver_Compute_Viscous_Volumes(s);
    Solver_Compute_Viscous_Edges(s);
  }

  if(s->limiting >0) Solver_Limiting(s);
  Solver_Compute_Inviscid_Volumes(s);
  Solver_ComputeFluxes(s);
  /*  if(s->Re==0) Solver_Inviscid_Stabilization(s); */

  if(itest > 0) Show_Inviscid_Term(s,itest);
  if(itest > 0) Show_Inviscid_Term1(s,itest);
}

void Solver_Integrate_EulerCauchy(tSolver_State *s){
  int const itest= -100;
  uint i,j,k;
  Solver_Integrate_Euler(s);
  for(i=0;i<s->nelem;i++)
   for(k=0;k<s->DOF;k++)
    for(j=0;j<4;j++){
      s->element[i].K1[k][j]=s->element[i].H[k][j];
      s->element[i].H[k][j]=0;
      s->element[i].w[k][j]-=s->element[i].tauloc*s->element[i].K1[k][j];
      }
  if(itest > 0) Show_Inviscid_Term(s,itest);
  if(itest > 0) Show_Inviscid_Term1(s,itest);
  Solver_Integrate_Euler(s);
  for(i=0;i<s->nelem;i++)
   for(k=0;k<s->DOF;k++)
    for(j=0;j<4;j++){
      s->element[i].H[k][j]=0.5*(s->element[i].K1[k][j]+s->element[i].H[k][j]);
      s->element[i].w[k][j]+= s->element[i].tauloc*s->element[i].K1[k][j];
      }
  if(itest > 0) Show_Inviscid_Term(s,itest);
  }

void Solver_Integrate_RK4(tSolver_State *s){
  uint i,j,k;
  Solver_Integrate_Euler(s);
  for(i=0;i<s->nelem;i++)
   for(k=0;k<s->DOF;k++)
    for(j=0;j<4;j++){
      s->element[i].K1[k][j]=s->element[i].H[k][j];
      s->element[i].H[k][j]=0;
      s->element[i].w[k][j]-=0.5*s->element[i].tauloc*s->element[i].K1[k][j];
      }
  Solver_Integrate_Euler(s);
  for(i=0;i<s->nelem;i++)
   for(k=0;k<s->DOF;k++)
    for(j=0;j<4;j++){
      s->element[i].K2[k][j]=s->element[i].H[k][j];
      s->element[i].H[k][j]=0;
      s->element[i].w[k][j]-=0.5*s->element[i].tauloc*(s->element[i].K2[k][j]-
        s->element[i].K1[k][j]);
      }
  Solver_Integrate_Euler(s);
  for(i=0;i<s->nelem;i++)
   for(k=0;k<s->DOF;k++)
    for(j=0;j<4;j++){
      s->element[i].K3[k][j]=s->element[i].H[k][j];
      s->element[i].H[k][j]=0;
      s->element[i].w[k][j]-=s->element[i].tauloc*(s->element[i].K3[k][j]-
        0.5*s->element[i].K2[k][j]);
      }
  Solver_Integrate_Euler(s);
  for(i=0;i<s->nelem;i++)
   for(k=0;k<s->DOF;k++)
    for(j=0;j<4;j++){
      s->element[i].H[k][j]=(s->element[i].K1[k][j]+2*(s->element[i].K2[k][j]+
        s->element[i].K3[k][j])+s->element[i].H[k][j])/6;
      s->element[i].w[k][j]-=-s->element[i].tauloc*s->element[i].K3[k][j];
      }
  }

void Solver_Solve(tSolver_State *s){
  uint i,j,k;
  double err,Mach,rro,rv2,pp;
  int itest = -168, itest0= -118;

  s->tau=1E+02;
  s->Machmax = 0;
  s->Machmin = 100;
 
//  printf("Begin of  Solver_Solve\n"); 

  /* setting of time step  */
  for(i=0;i<s->nelem;i++){
    double oneoverrho,norm,dt,lambda;
    s->element[i].tauloc=1E+02;
    for(k=0;k<s->DOF;k++){
      for(j=0;j<4;j++)  /* set flux to zero */
        s->element[i].H[k][j]=0;
      if(s->Re > 0)
	for(j=0;j<3;j++)  /* set flux to zero */
	  s->element[i].HV[k][j]=0;


      oneoverrho=1/s->element[i].w[k][0];
      norm=s->element[i].w[k][1]*s->element[i].w[k][1]+
           s->element[i].w[k][2]*s->element[i].w[k][2];
      lambda=sqrt(s->kappa2*oneoverrho*
		 (s->element[i].w[k][3]- 0.5*oneoverrho*norm))+ 
	    oneoverrho*sqrt(norm); /* maximal eigenvalue  v+c */
      
      dt=1/(s->element[i].maxovervolume*2*lambda)/s->DOF; /* inviscid part */  
      if(s->Re > 0)
        dt = min(dt,s->Re/s->element[i].oneovervolume/8/s->DOF); /* viscous */
      /*       2, 8 number of inviscid, viscous terms in the DGM method */
      /*       DOF from time derivatives   */
      if(s->tau>dt)
	s->tau=dt;
      if(s->element[i].tauloc>dt)
        s->element[i].tauloc=dt;
    }
    s->element[i].tauloc = s->CFL*s->element[i].tauloc;

  }
  s->tau=s->CFL*s->tau;
  
  //printf("Time step: %7.5g\n",s->tau);

  if(s->t+s->tau>s->nextt)
    s->tau=s->nextt-s->t;

  if(!s->local)
    for(i=0;i<s->nelem;i++)
      s->element[i].tauloc = s->tau;

  fflush(stdout); 

  /* time integration: \partial w/\partial t = F(w) */
  /*    printf("Before integrating\n");  */
  for(i=0;i<s->nelem;i++)
    for(j=0;j<s->DOF;j++)
      for(k=0;k<4;k++)
        s->element[i].wold[j][k] = s->element[i].w[j][k];

  //printf("Before integrate \n");
  s->integrate(s);
  //printf("After integrate \n");
  /* upgrading the unknowns */
  for(i=0;i<s->nelem;i++){

    if(i == itest || i == itest0)
        printf("Time steps %d %10.5g %10.5g %10.5g %10.5g %10.5g \n",
        i,s->element[i].tauloc,s->tau,s->Re/s->element[i].oneovervolume/3,
        1/s->element[i].oneovervolume,1/s->element[i].maxovervolume );

    if(i == itest )Show_Inviscid_Term(s,i);

    for(k=0;k<s->DOF;k++){
      for(j=0;j<4;j++)
	s->element[i].w[k][j]-=s->element[i].tauloc*s->element[i].H[k][j];
    }

     if(s->element[i].limited)  /* post limiting */
       Solver_Limiting_Element(s,i,1); 
        
    for(k=0;k<s->DOF;k++){
      rro = 1/ s->element[i].w[k][0];
      rv2 = rro*(s->element[i].w[k][1]*s->element[i].w[k][1] + 
          s->element[i].w[k][2]*s->element[i].w[k][2]);
      pp = (s->kappa-1)*(s->element[i].w[k][3]-0.5*rv2);

      if(pp <=0){
         pp = 0.001;
         s->element[i].w[k][3] = pp/(s->kappa-1)+0.5*rv2;
/*      }

      if(pp <=-0E+20) { */
	printf("%7.5g %7.5g triangle with nonpositive pressure %d  %d %7.5g  %7.5g \n",
	       s->element[i].ve->P[0]->Point[0],
	       s->element[i].ve->P[0]->Point[1],i,k,rv2,pp);
	printf("%7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",
	       s->element[i].ve->P[1]->Point[0],
	       s->element[i].ve->P[1]->Point[1],s->element[i].w[k][0],
	       s->element[i].w[k][1],s->element[i].w[k][2],
	       s->element[i].w[k][3]);
	printf("%7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n\n",
	       s->element[i].ve->P[2]->Point[0],
	       s->element[i].ve->P[2]->Point[1],s->element[i].H[k][0],
	       s->element[i].H[k][1],s->element[i].H[k][2],
	       s->element[i].H[k][3],s->element[i].HV[k][0],
	       s->element[i].HV[k][1],s->element[i].HV[k][2]);
      }

      Mach=  sqrt(rv2/(s->kappa*pp));
      s->Machmin = min(s->Machmin,Mach);
      s->Machmax = max(s->Machmax,Mach);

    }  /* end of k */
    /* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  */
    /* Solver_Limiting_Element(s,i);  */

      if(i == itest )Show_Inviscid_Term(s,i);

  }  /* end of i */
  s->t+=s->tau;


  /* computing of reziduum  */
  s->linfty=0;
  for(i=0;i<s->nelem;i++){
 
    err = 0;
    for(k=0;k<s->DOF;k++){
    /* L_infinity norm  */
    /* err= fabs((s->element[i].w[k][0]-s->element[i].wold[k][0])/
           s->element[i].tauloc);
      if(s->linfty < err)
	s->linfty= err; */

    /* L^1 norm  */
     err+= fabs((s->element[i].w[k][0]-s->element[i].wold[k][0])/
           s->element[i].tauloc);
    }
    s->linfty+= err/ s->element[i].oneovervolume;

    /* if(err/ s->element[i].oneovervolume > 1E-03){
      //printf("Err: %d %7.5g %7.5g %7.5g\n",i,err,
      //     err/ s->element[i].oneovervolume,s->linfty/3); 
      printf("\n %7.5g %7.5g %7.5g  %7.5g %7.5g  %7.5g \n",
	     s->element[i].ve->P[0]->Point[0],
	     s->element[i].ve->P[0]->Point[1],
	     s->element[i].w[0][0],s->element[i].w[0][1],
	     s->element[i].w[0][2],s->element[i].w[0][3]);
      printf("%7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",
	     s->element[i].ve->P[1]->Point[0],
	     s->element[i].ve->P[1]->Point[1],
	     s->element[i].w[1][0],s->element[i].w[1][1],
	     s->element[i].w[1][2],s->element[i].w[1][3]);
      printf("%7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",
	     s->element[i].ve->P[2]->Point[0],
	     s->element[i].ve->P[2]->Point[1],
	     s->element[i].w[2][0],s->element[i].w[2][1],
	     s->element[i].w[2][2],s->element[i].w[2][3]);
	     } */
  }
  s->linfty/=s->DOF;

}

void Solver_MeshUpdate(tSolver_State *s){
  uint i,j,k;
  double2 coef, coef1;
  Progress_Init("Updating mesh",s->nelem);
  for(i=0;i<s->nelem;i++){
    for(k=0;k<s->DOF;k++)
      for(j=0;j<4;j++)
        s->element[i].ve->w[k][j]=s->element[i].w[k][j];
    Progress_Update(i+1);
    }
  Progress_Done();
  if(s->Re>0){
    Solver_Compute_Viscous_Volumes(s);   /* in order to compute Dx,Dy */
    Coeficients(s,coef, coef1);
    printf("Drag coeficient c_D=%7.5g, lift coeficient c_L=%7.5g \n",
    coef[0],coef[1]); 
    printf("Drag coef-pres c_Dp=%7.5g, drag coef-visc c_Dv=%7.5g \n",
    coef1[0],coef1[1]); 
    s->cL = coef[1];
    s->cD = coef[0];
    s->cDp = coef1[0];
  }
}

void Plot_Walls(tSolver_State *s){
  FILE *f1;
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};
  uint j;

  f1=fopen("walls","wb");

  Progress_Init("Plotting walls",s->nelem);
  for(j=0;j<s->nface;j++){
    uint i,ii;
    const uint k=0;   /*component of state vector */
    double wL, wR;
    i=s->face[j].i;
    if(s->face[j].boundary>0)
/*      if(s->boundary[s->face[j].boundary-1].type==cSolver_BoundaryType_Inlet ||
   s->boundary[s->face[j].boundary-1].type==cSolver_BoundaryType_Outlet){*/
      if(s->boundary[s->face[j].boundary-1].type==cSolver_BoundaryType_Slip){
 
          ii= s->face[j].ii;   /*computation of values of w[k] in nodes of T */
          wL= s->element[i].w[perm[ii][0]][k] +
              s->element[i].w[perm[ii][2]][k] -
              s->element[i].w[perm[ii][1]][k]; 
          wR= s->element[i].w[perm[ii][0]][k] +
              s->element[i].w[perm[ii][1]][k] -
              s->element[i].w[perm[ii][2]][k]; 
         fprintf(f1, " %10.8g %10.8g \n %10.8g %10.8g \n \n",
                 s->element[i].ve->P[perm[ii][0]]->Point[0], wL,
                 s->element[i].ve->P[perm[ii][1]]->Point[0], wR); 
      }
    Progress_Update(i+1);
  }
  Progress_Done();
  close(f1);
  }

int ViscousTerm(double4 w, double4 phi, 
                double2 nablarho, double4 phix, double4 phiy, 
	        double Re, double kappaRePr, double4 D1, double4 D2 ){
  int ier=-1;
  double mirho, derrho[2];

  mirho = 1/Re/w[0];
  derrho[0] = nablarho[0]/w[0];
  derrho[1] = nablarho[1]/w[0];
/*  derrho[0] = 0;
  derrho[1] = 0; */

  D1[0] = 0;
  D1[1] = (2*phix[1]-phiy[2] - 2*phi[1]*derrho[0]+phi[2]*derrho[1])*2*mirho/3;
  D1[2] = (phix[2]+phiy[1] - phi[2]*derrho[0]-phi[1]*derrho[1])*mirho;
  D1[3] = (D1[1]*w[1]+D1[2]*w[2] + 
           kappaRePr*(phix[3]-(w[1]*phix[1]+w[2]*phix[2])/w[0]
           +( (phi[1]*w[1]+phi[2]*w[2])/w[0] -phi[3])*derrho[0]) )/w[0];
  D2[0] = 0;
  D2[1] = D1[2];
  D2[2] = (2*phiy[2]-phix[1] - 2*phi[2]*derrho[1]+phi[1]*derrho[0])*2*mirho/3;
  D2[3] = (D2[1]*w[1]+D2[2]*w[2] + 
           kappaRePr*(phiy[3] -(w[1]*phiy[1]+w[2]*phiy[2])/w[0]
           +( (phi[1]*w[1]+phi[2]*w[2])/w[0]-phi[3])*derrho[1]))/w[0];
  return ier;
}


int ViscousTerm_NEW(double4 w, double4 phi, 
                double2 nablarho, double4 phix, double4 phiy, 
	        double Re, double kappaRePr, double4 D1, double4 D2 ){
  int ier=-1;
  double mirho, vx,vy,vz, v2; 

  mirho = 1/Re/w[0];
  vx = w[1]/w[0];
  vy = w[2]/w[0];
  vz = w[3]/w[0];
  v2 = vx*vx + vy*vy;

  D1[0] = 0;
  D1[1] = (2*phix[1]-phiy[2] - 2*phix[0]*vx + phiy[0]*vy)*2*mirho/3;
  D1[2] = (phix[2]+phiy[1] - phix[0]*vy - phiy[0]*vx)*mirho;
  D1[3] = D1[1]*vx + D1[2]*vy + 
           kappaRePr*(phix[3]-phix[0]*vz - phix[1]*vx - phix[2]*vy
           - v2*phix[0] ) /w[0];
  D2[0] = 0;
  D2[1] = D1[2];
  D2[2] = (2*phiy[2]-phix[1] - 2*phiy[0]*vy + phix[0]*vx)*2*mirho/3;
  D2[3] = D2[1]*vx + D2[2]*vy + 
           kappaRePr*(phiy[3]-phiy[0]*vz - phiy[1]*vx - phiy[2]*vy
           - v2*phiy[0] ) /w[0];
/*  printf("DD: %7.5g  %7.5g  %7.5g %7.5g  %7.5g \n",
          D1[1],D1[2],D2[2],D1[3],D2[3]); */
  return ier;
}

BOOL Identical_Nodes( double2 a, double2 b){
  static double zero=1E-05;
  double norm;
  norm = sqrt(double2_distance2(a,b)/double2_norm2(a)/double2_norm2(b));
  /*  printf(" norm [%7.5g, %7.5g],  [%7.5g, %7.5g], %7.5g \n",
      a[0],a[1],b[0],b[1],norm); */
  if(norm < zero) 
    return TRUE;
  else
    return FALSE;
}
void Coeficients(tSolver_State *s, double2 coef, double2 coef1){
  uint i,ii,j,k,l;
  int ic;
  double const maxcosangle=-0.8; /* angle between velocity and normal vector */
  double rho_inf=0,weight=0,cosangle;
  double2 v_inf,v;
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};
  static uint perm2[2][2]={{0,1},{1,0}};
  const double oneoversqrt3 = 0.57735026918962584;
  double a[2];

  a[0] = (1/sqrt(3) + 1)/2;
  a[1] = (-1/sqrt(3) + 1)/2;

 v_inf[0]=v_inf[1] = 0;
/* computing the values at infinity */
 for(k=0;k<s->nface;k++){
   if(s->face[k].boundary > 0 && 
      s->boundary[s->face[k].boundary-1].type==cSolver_BoundaryType_Inlet){
     i = s->face[k].i;
     ii= s->face[k].ii;
     for(j=0;j<2;j++) v[j]=s->element[i].w[ii][j+1]/s->element[i].w[ii][0];
     cosangle=double2_dotprod(v,s->face[k].n)/sqrt(double2_norm2(v));
     if(cosangle < maxcosangle){
       rho_inf+=s->element[i].w[ii][0]*s->face[k].volume;
       for(j=0;j<2;j++) v_inf[j]+=v[j]*s->face[k].volume;
       weight+=s->face[k].volume;
       /*Write_Edge(s,k); */
     }
   }
 }
 rho_inf/=weight;
 for(j=0;j<2;j++) v_inf[j]/=weight;

 /*printf("At infty %7.5g %7.5g %7.5g \n",rho_inf,v_inf[0],v_inf[1]); */

 coef[0] = 0;
 coef[1] = 0;
 coef1[0] = 0;
 coef1[1] = 0;
 for(k=0;k<s->nface;k++){
   double4 Dx,Dy;
   if(s->face[k].boundary > 0 && 
      s->boundary[s->face[k].boundary-1].type==cSolver_BoundaryType_Slip){
     i = s->face[k].i;
     ii= s->face[k].ii;

     ic=s->element[i].icve;	
     if(ic >= 0){
       double pressure;
       double4 wi[2],wpom;

       //       printf("Curved: %d %d %d \n",k,i,ic);
       double4_sub(wpom,s->element[i].w[perm[ii][2]],
		   s->element[i].w[perm[ii][1]]);
       double4_multiply(wpom, oneoversqrt3, wpom);
       double4_add(wi[0],s->element[i].w[perm[ii][0]],wpom);
       double4_sub(wi[1],s->element[i].w[perm[ii][0]],wpom);

       for(j=0;j<2;j++){
	 double2_weight_awer(Dx,s->element[i].Dsxnablaw[perm[ii][0]],
			     s->element[i].Dsxnablaw[perm[ii][1]],
			     a[perm2[0][j]],a[perm2[1][j]]);
	 double2_weight_awer(Dy,s->element[i].Dsynablaw[perm[ii][0]],
			     s->element[i].Dsynablaw[perm[ii][1]],
			     a[perm2[0][j]],a[perm2[1][j]]);

	 // complete part of coefficients

	 coef[0]-= ((Dx[1]-Pressure(wi[j],s->kappa))
		    *s->curvedelement[ic].normal[j][0]+
		    Dy[1]*s->curvedelement[ic].normal[j][1])*
	   s->curvedelement[ic].length[j]/2;
	 coef[1]-= (Dx[2]*s->curvedelement[ic].normal[j][0]+
		    (Dy[2]-Pressure(wi[j],s->kappa))*
		    s->curvedelement[ic].normal[j][1])*
	   s->curvedelement[ic].length[j]/2;
	 /*       if(k == 1354) {
	  printf("Curved: %d %d %d \n",k,i,ic);
	  printf("Dx: %7.5g %7.5g %7.5g %7.5g \n",Dx[0],Dx[1],Dx[2],Dx[3]);
	  printf("Dy: %7.5g %7.5g %7.5g %7.5g \n",Dy[0],Dy[1],Dy[2],Dy[3]);
	  printf("pres, nn: %7.5g %7.5g %7.5g %7.5g \n",
		 Pressure(wi[j],s->kappa),
		 s->curvedelement[ic].normal[j][0],
		 s->curvedelement[ic].normal[j][1],
		 s->curvedelement[ic].length[j]);
	  printf("coeff:  %12.9g %12.9g \n",
		 ((Dx[1]-Pressure(wi[j],s->kappa))
		  *s->curvedelement[ic].normal[j][0]+
		  Dy[1]*s->curvedelement[ic].normal[j][1])*
		 s->curvedelement[ic].length[j],
		 (Dx[2]*s->curvedelement[ic].normal[j][0]+
		  (Dy[2]-Pressure(wi[j],s->kappa))*
		  s->curvedelement[ic].normal[j][1])*
		 s->curvedelement[ic].length[j]);
	  printf("*************************************** \n");
       }
	 */
     // inviscid part of coefficient
	 Dx[1]=Dx[2]=Dy[1]=Dy[2] = 0;
	 coef1[0]-= ((Dx[1]-Pressure(wi[j],s->kappa))
		    *s->curvedelement[ic].normal[j][0]+
		    Dy[1]*s->curvedelement[ic].normal[j][1])*
	   s->curvedelement[ic].length[j]/2;
       }
     }
     else{
       double4_average(Dx,s->element[i].Dsxnablaw[perm[ii][0]],
		       s->element[i].Dsxnablaw[perm[ii][1]]);
       double4_average(Dy,s->element[i].Dsynablaw[perm[ii][0]],
		       s->element[i].Dsynablaw[perm[ii][1]]); 
       
     // complete part of coefficients
       coef[0]-= ((Dx[1]-Pressure(s->element[i].w[ii],s->kappa))*
		  s->face[k].n[0]+
		  Dy[1]*s->face[k].n[1])*s->face[k].volume;
       coef[1]-= (Dx[2]*s->face[k].n[0]+
		  (Dy[2]-Pressure(s->element[i].w[ii],s->kappa))*
		  s->face[k].n[1])
	 *s->face[k].volume;
       // inviscid part of coefficient
       Dx[1]=Dx[2]=Dy[1]=Dy[2] = 0;
       coef1[0]-=((Dx[1]-Pressure(s->element[i].w[ii],s->kappa))*
		  s->face[k].n[0]+
		  Dy[1]*s->face[k].n[1])*s->face[k].volume;
       /*     coe1f[1]-=(Dx[2]*s->face[k].n[0]+
	      (Dy[2]-Pressure(s->element[i].w[ii],s->kappa))*s->face[k].n[1])
	      *s->face[k].volume;
	      */
     }
   }
 }
 for(j=0;j<2;j++) coef[j]*=2/rho_inf/double2_norm2(v_inf);
 for(j=0;j<1;j++) coef1[j]*=2/rho_inf/double2_norm2(v_inf);
 coef1[1] = coef[0] - coef1[0];
}

void Sod_Error(tSolver_State *s,double rhoL,double v1L,double v2L,double pL,
               double rhoR,double v1R,double v2R,double pR){
double error_L2=0,error_L2_smooth=0, error_local, hmax;
double rho,v,p,e,x,y, x1,time, kappa,x0, p2,cl,cr, rleft,vc,vs, rhstep;
uint i,j,ihstep;

  x0=0;
  time = s->t;
  kappa = s->kappa;
  rhstep = 0.01;
  hmax = 0;

  cl = pow_fast(kappa*pL/rhoL,0.5);
  cr = pow_fast(kappa*pR/rhoR,0.5);
  p2 = 0.3031302;

  rleft = pow_fast(2./kappa/(kappa-1),0.5) 
          *(p2/pR -1)/pow_fast(1+(kappa+1)/(kappa-1)*p2/pR,0.5);
  vc = cr*rleft;
  vs = cr*cr*(p2/pR-1)/kappa/vc;
  x1 = x0 + ((kappa +1)/2*vc-cl)*time;

  for(i=0;i<s->nelem;i++){
    error_local = 0;
    ihstep = 1;
    hmax = max(hmax,s->element[i].maxovervolume/s->element[i].oneovervolume);
    for(j=0;j<3;j++){
      x=s->element[i].ve->P[j]->Point[0];
      if( fabs(x0+vc*time - x) < rhstep || fabs(x0+vs*time-x)<rhstep ) 
        ihstep = 0;  /* near discontinuity  */

      if(x < x0-cl*time){ 
        p = pL;
        v = 0;
        rho = rhoL;
      }
      else if(x > x0-cl*time &&  x < x1) {
        p = pL*pow_fast(1-(kappa-1)/(kappa+1)*(1+(x-x0)/cl/time),
            2*kappa/(kappa-1));
        v = 2/(kappa+1)*((x-x0)/time + cl);
        rho =rhoL*pow_fast(1-(kappa-1)/(kappa+1)*(1+(x-x0)/cl/time),
              2/(kappa-1));
      }
      else if((x>x1)&&( x < x0 + vc*time)){ 
        p = p2;
        v = vc;
        rho = rhoL*pow_fast(p2/pL,1./kappa);
      }
      else if((x>x0+vc*time)&&(x < x0 + vs*time)) {
        p = p2;
        v = vc;
        rho = rhoR*(1+(kappa+1)/(kappa-1)*p2/pR)/
             ((kappa+1)/(kappa-1) + p2/pR);
      }
      else if(x > x0+vs*time){
        p = pR;
        v = 0;
        rho = rhoR;
      }
      else{
        printf ("discontinuites on  element %d %d\n", i,i);
        break;
      }

/*      if(ihstep == 0) printf("%7.5g %7.5g \n",x,rho); */

      error_local+=(rho-s->element[i].w[j][0])*(rho-s->element[i].w[j][0]);

    }  /* end of j */
    error_L2+=error_local/s->element[i].oneovervolume/3;
    if(ihstep == 1)
       error_L2_smooth+=error_local/s->element[i].oneovervolume/3;
  }
  error_L2=pow_fast(error_L2,0.5);
  error_L2_smooth=pow_fast(error_L2_smooth,0.5);

  printf("||u-u_h||_L^2 = %10.5e(whole), %12.5e(smooth), h_max=%7.5g(%d el)\n",
     error_L2,error_L2_smooth,hmax,s->nelem);

}

void calc_ringleb(double gamma,double2 coords, double4 soln){
  double cspeed;
  
  cspeed = cspeed_ringleb(gamma,coords[0],coords[1] );
  state_ringleb(gamma,coords[0],cspeed,soln);
}

double cspeed_ringleb(double gamma,double x,double y){
  double tol, clow, chigh, cnew, err;
  int nit;

  tol = 1.e-15;
  clow  = .5000000;
  chigh = .9999999;

  nit = 0;
  err = 1.0;
  while(nit < 100 && err > tol){
    cnew = (clow+chigh)/2.;
    if( fun_ringleb(gamma,x,y,cnew)*fun_ringleb(gamma,x,y,chigh) > 0.)
      chigh = cnew;
    else
      clow  = cnew;

    nit++;
    err  = fabs(chigh-clow);
  }

  if(err > tol){
    printf("Error tolerance was not achieved in bisection iteration: %7.5g > %7.5g\n",err,tol);
    exit;
  }
  return cnew;
}
 
double fun_ringleb(double gamma,double x,double y,double c){
  double jval, gb, rho, q2;

  gb = gamma - 1;
  rho = pow_fast(c,(2/gb));
  jval = 1./c + 1./(3.*pow_fast(c,3)) + 1./(5.*pow_fast(c,5)) - 
    0.5*log((1+c)/(1.-c) );
  q2 = 2.*(1.-c*c)/gb;
  return pow2(x-jval/2.) + y*y - 1./(4.*pow2(rho)*pow2(q2));
}

void state_ringleb(double gamma,double x,double c,double4 soln){
  double gb, jval, rho, qval, kval, u,v,p,en,q2;

  gb    =  gamma - 1.;

  rho = pow_fast(c,(2/gb));
  jval = 1./c + 1./(3.*pow_fast(c,3)) + 1./(5.*pow_fast(c,5)) - 
    0.5*log((1+c)/(1.-c) );
  q2   = 2.*(1.-c*c)/gb;
  qval = sqrt(q2);
  kval = sqrt( 2./(1./q2 -2.*rho*(x - jval/2.)));

  v     = q2/kval;
  u     = sqrt( fabs(q2 - v*v) );
  p     = c*c*rho/gamma;
  en    = p/gb + 0.5*rho*(u*u + v*v);
      
  soln[0] = rho;
  soln[1] = u;
  soln[2] = v;
  soln[3] = p;
}

void  Ringleb_Error(tSolver_State *s){
  double error_L2=0, norm_u =0, error_local, norm_local, hmax;
  double4 q_exact, w_exact,w_local;
  statevec w_vertices;
  double weight[10], baryc[3][10], x[3][2], x_local[2];
  uint i,j,k,l, order, ng;
  const iw=3;  /* first component (density) of w is considered */
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};

  /* weights and barycentric coordinates for quadrature rules */
  order = 5; /* maximum is order = 5 */

  if(order == 1) {
    ng = 1;
    weight[0] = 1.;
    baryc[0][0] = 1./3;
    baryc[1][0] = 1./3;
    baryc[2][0] = 1./3;
  }
  else if(order == 2) {
    ng = 3;
    weight[0] = 1./3.;
    baryc[0][0] = 4./6;
    baryc[1][0] = 1./6;
    baryc[2][0] = 1./6;
    weight[1] = 1./3.;
    baryc[0][1] = 1./6;
    baryc[1][1] = 4./6;
    baryc[2][1] = 1./6;
    weight[2] = 1./3.;
    baryc[0][2] = 1./6;
    baryc[1][2] = 1./6;
    baryc[2][2] = 4./6;
  }
  else if(order == 3) {
    ng = 4;
    weight[0] = -0.5625;
    baryc[0][0] = 1./3.;
    baryc[1][0] = 1./3.;
    baryc[2][0] = 1./3.;
    weight[1] = 5207.8125/9999.;
    baryc[0][1] = 3./5.;
    baryc[1][1] = 1./5.;
    baryc[2][1] = 1./5.;
    weight[2] = weight[1];
    baryc[0][2] = 1./5.;
    baryc[1][2] = 3./5.;
    baryc[2][2] = 1./5.;
    weight[3] = weight[1];
    baryc[0][3] = 1./5.;
    baryc[1][3] = 1./5.;
    baryc[2][3] = 3./5.;
  }
  else if(order == 4) {
    ng = 6;
    weight[0] = 0.223381589678011;
    baryc[0][0] = 0.108103018168070;
    baryc[1][0] = 0.445948490915965;
    baryc[2][0] = 0.445948490915965;
    weight[1] = 0.223381589678011;
    baryc[1][1] = 0.108103018168070;
    baryc[2][1] = 0.445948490915965;
    baryc[0][1] = 0.445948490915965;
    weight[2] = 0.223381589678011;
    baryc[2][2] = 0.108103018168070;
    baryc[0][2] = 0.445948490915965;
    baryc[1][2] = 0.445948490915965;
    weight[3] = 0.109951743655322;
    baryc[0][3] = 0.816847572980459;
    baryc[1][3] = 0.091576213509771;
    baryc[2][3] = 0.091576213509771;
    weight[4] = 0.109951743655322;
    baryc[1][4] = 0.816847572980459;
    baryc[2][4] = 0.091576213509771;
    baryc[0][4] = 0.091576213509771;
    weight[5] = 0.109951743655322;
    baryc[2][5] = 0.816847572980459;
    baryc[0][5] = 0.091576213509771;
    baryc[1][5] = 0.091576213509771;
    }
  else if(order == 5) {
    ng = 7;
    weight[0] = 0.225;
    baryc[0][0] = 1./3.;
    baryc[1][0] = 1./3.;
    baryc[2][0] = 1./3.;
    weight[1] = 0.132394152788506;
    baryc[0][1] = 0.059715871789770;
    baryc[1][1] = 0.470142064105115;
    baryc[2][1] = 0.470142064105115;
    weight[2] = 0.132394152788506;
    baryc[1][2] = 0.059715871789770;
    baryc[2][2] = 0.470142064105115;
    baryc[0][2] = 0.470142064105115;
    weight[3] = 0.132394152788506;
    baryc[2][3] = 0.059715871789770;
    baryc[0][3] = 0.470142064105115;
    baryc[1][3] = 0.470142064105115;
    weight[4] = 0.125939180544827;
    baryc[0][4] = 0.79742698535308;
    baryc[1][4] = 0.101286507323456;
    baryc[2][4] = 0.101286507323456;
    weight[5] = 0.125939180544827;
    baryc[1][5] = 0.79742698535308;
    baryc[2][5] = 0.101286507323456;
    baryc[0][5] = 0.101286507323456;
    weight[6] = 0.125939180544827;
    baryc[2][6] = 0.79742698535308;
    baryc[0][6] = 0.101286507323456;
    baryc[1][6] = 0.101286507323456;
  }
  else
    printf("Bad integration role, order=%d\n",order);

  /*  for(i=0;i<ng;i++)
    printf("&& %d %d %7.5g %7.5g %7.5g %7.5g \n",
	   i,ng,weight[i],baryc[0][i],baryc[1][i],baryc[2][i]);
  */  
  for(i=0;i<s->nelem;i++){ 
    //  for(i=0;i<1;i++){
    for(j=0;j<3;j++){
      x[j][0]=s->element[i].ve->P[j]->Point[0]; /* coordinates of vertices */
      x[j][1]=s->element[i].ve->P[j]->Point[1];
      /*      printf("%7.5g %7.5g !!!\n",x[j][0],x[j][1]); */
      for(k=0;k<4;k++){  /* values of w in vertices */
	w_vertices[j][k] = s->element[i].w[perm[j][0]][k] + 
	  s->element[i].w[perm[j][2]][k] - s->element[i].w[perm[j][1]][k];
      }
    }
    error_local = 0;
    norm_local = 0;
    for(j=0;j<ng;j++){ /* we do through integration nodes */
      for(k=0;k<4;k++){
	w_local[k] = 0;
	for(l=0;l<3;l++)
	  w_local[k] += w_vertices[l][k]*baryc[l][j];
      }
      x_local[0]= x_local[1]=0;
      for(l=0;l<3;l++){
	x_local[0] += x[l][0]*baryc[l][j];
	x_local[1] += x[l][1]*baryc[l][j];
      }
      calc_ringleb(s->kappa,x_local,q_exact);
      vector_q2w(q_exact,w_exact,s->kappa);
      for(k=0;k<iw+1;k++){
	error_local += weight[j]*pow2(w_local[k] - w_exact[k]);
	norm_local += weight[j]*pow2(w_exact[k]);
      }
      /*      printf("%7.5g %7.5g %7.5g %7.5g %7.5g %d %d %d %7.5g %7.5g\n",
	     x_local[0],x_local[1],w_local[iw],w_exact[iw],error_local,
	     i,j,ng,weight[j],1./s->element[i].oneovervolume); */
    }
    error_L2+=error_local/s->element[i].oneovervolume;
    norm_u += norm_local/s->element[i].oneovervolume;
  }
  error_L2=sqrt(error_L2);
  norm_u = sqrt(norm_u);
  printf("Error in L^2 norm (absolut, relative): %7.5e %7.5e\n",
	 error_L2, error_L2/norm_u);
}
