#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"


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;
  }

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

  for(i=0;i<s->nface;i++){
    uint _i,_j,_ii,_jj;
    double4 wi1, wj1, H1, wi2, wj2, H2, wpom;  /*two integration nodes */
/*    double4 H; */
    double2 Hper;  //keeps rotated H[1]..H[2] for periodic b.c.
    
    /*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);

    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);
         
       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]);

       }

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

       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{  //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("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]); */

/*        Triangle_VectorRotation(p1->PeriodicCopy->Point,
          p2->PeriodicCopy->Point,
          Mesh_EdgeOpposite(vej,p1->PeriodicCopy,p2->PeriodicCopy)->Point,
          p1->Point,p2->Point,p3->Point,v1, w+1); //rotate velocity vector v_j
   */
        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];
          }
        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);
        }
      else{  //slip
        double p;
        /*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;
	  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;
	}
	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 ){ 
      printf("\nFlux %d %d %6.4e %6.4e %6.4e %6.4e %6.4e %6.4e \n",i,k, 
          s->element[i].H[0][0],
          s->element[i].H[1][0],
          s->element[i].H[2][0],0,0,0);
      for(k=1;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-1],-s->element[i].HV[1][k-1],
	       -s->element[i].HV[2][k-1]); 
        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<3;j++)
        s->element[i].H[k][j+1]-=s->element[i].HV[k][j]; 
}

}


void Solver_Compute_Inviscid_Volumes(tSolver_State *s){
  uint i,j,k;
  int ic;
  const int itest = -355,itest0=-366;
  double rho,p,e;
  double2 v,n;
  double4 VI[s->DOF];
  
  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 coordiante %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_Compute_Viscous_Edges(tSolver_State *s){
  int i,j;
  uint k,l;
  int ic;
  const int lastcomponent=3;   /* =3  or =2, in  HV[.][lastcomponent-1]  */
  const int itest = -355, itest0 =-366;
  double3 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<3;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<3;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<3;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;
    double4 wl,wr,wDl,wDr;
    double2 noutL,noutR;
    double3 theta,detJF;
    lintransf local_dertest[3];
    double3 Dwnablaphix[s->DOF][s->DOF][3],Dwnablaphiy[s->DOF][s->DOF][3];  
    /*                     1. integ. node, 2. phi, 3. compon. phi*/

    /* computation of <D(w)nablaphi> in all integration nodes */

    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];
          }
      if(ic>=0)   /* curved element */
        detJF[j]=2/s->curvedelement[ic].detJF[j];
      else
        detJF[j]=s->element[i].oneovervolume;
    }
    Compute_Dw_nablaPhi(s->DOF,s->element[i].w,s->element[i].nablarho,
                        detJF,Re, kappaRePr,
                        local_dertest, Dwnablaphix, Dwnablaphiy);
    if(i == itest|| i == itest0){
      for(j=0;j<1;j++)
        for(k=0;k<3;k++)
	  for(m=0;m<3;m++){
	    printf("Dwnablaphix[%d][%d][%d]=(%7.5g,%7.5g,%7.5g)\n",
		   j,k,m,Dwnablaphix[j][k][m][0],Dwnablaphix[j][k][m][1],
		   Dwnablaphix[j][k][m][2]);
	  }
      for(j=0;j<1;j++)
        for(k=0;k<3;k++)
	  for(m=0;m<3;m++){
	    printf("Dwnablaphiy[%d][%d][%d]=(%7.5g,%7.5g,%7.5g)\n",
		   j,k,m,Dwnablaphiy[j][k][m][0],Dwnablaphiy[j][k][m][1],
		   Dwnablaphiy[j][k][m][2]);
	  }
    } 

    for(j=0;j<3;j++){ /* j goes over edges  */
      _i=s->element[i].neigh[j];
      ic=s->element[i].icve;
      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 */
        area = 0.5/Re;   /* 0.5 is the weight */
        for(l=0;l<3;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<3;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]]);
            if(i == itest || i == itest0) 
             printf("Temperaures %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;
            }
            wDl[3] = wDl[0]*(theta[0]+theta[2]-theta[1]); 
            wDr[3] = wDr[0]*(theta[0]+theta[1]-theta[2]); 
/*            wDl[3] = wDl[0]*(3*theta[1]-theta[0])/2; 
            wDr[3] = wDr[0]*(3*theta[0]-theta[1])/2; */
/*            wDl[3] = wDl[0]*(theta[1]+theta[0])/2;
            wDr[3] = wDr[0]*(theta[0]+theta[1])/2; */
         }
         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  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  */
        double3_vertex4(wl, s->element[i].w[perm[j][0]], 
		    s->element[i].w[perm[j][2]],s->element[i].w[perm[j][1]]);
        double3_vertex4(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<3;n++){
          wjumpL[n]=wl[n]-wDl[n+1]; 
          wjumpR[n]=wr[n]-wDr[n+1]; 
        }
       /* computation of w phi */
        area = 0.5/Re;   /* 0.5 is the weight */
       /* computation for test functions having support on edge k  */
        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){
	  printf("\n wl,    wR     : %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",
                  wl[0],wl[1],wl[2],wr[0],wr[1],wr[2]);
	  printf("\n wDl,   wDR    : %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",
                  wDl[1],wDl[2],wDl[3],wDr[1],wDr[2],wDr[3]);
	  printf("\n wjumpL, wjumpR: %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",
                  wjumpL[0],wjumpL[1],wjumpL[2],wjumpR[0],wjumpR[1],wjumpR[2]);
	  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,ibd);
      }

      if(ibd == 0) printf("BAD CYCLING!!!!!! \n");
      for(k=0;k<s->DOF;k++){  /* k goes over test function */
        for(l=0;l<3;l++){  
	  /*              l goes over component of test function phi*/
	  /* computation of <Dw nabla phi> over edge j */
	  for(m=0;m<3;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<3;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<3;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))
            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 == itest ){
       printf("\n PART OF  <D(w)nabla phi>[w]  -  <2> or (2)\n");
       Show_Viscous_Term(s,itest);
      }
    } /* end of j */
  } /* end of i */

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

  // 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<3;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<3;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 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<3;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,
			 double3 oneoverarea, double Re,  double kappaRePr, 
			 lintransf dertestfunction[3],
			 double3 Dwnablaphix[3][3][3],
                         double3 Dwnablaphiy[3][3][3]){ 
/*                              originaly [DOF][DOF][3]  */
  uint j,k,l,m;
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};
  double3 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<3;l++){  /* l goes over component of test function */
      for(m=0;m<3;m++){
	nablaphix[m]=0;
	nablaphiy[m]=0;
        phi[m+1] = 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] = 1;  /* only hear is nonzero */
	ViscousTerm(w[j],phi,nablarho,nablaphix,nablaphiy,Re,kappaRePr,
		    wlocx[j],wlocy[j]);
      }
      for(j=0;j<3;j++){ 
	double3_vertex(Dwnablaphix[j][k][l],wlocx[perm[j][0]],
		       wlocx[perm[2][0]],wlocx[perm[j][1]]);
	double3_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,
		   double3 wjumpL, double3 wjumpR){
  uint l;
  double3 wloc;
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};

  /* element K_i  */
  double3_vertex4(wjumpL, wi[perm[i][0]], wi[perm[i][2]], wi[perm[i][1]]);
  double3_vertex4(wjumpR, wi[perm[i][0]], wi[perm[i][1]], wi[perm[i][2]]);
  /* element K_j, neighboor of K_i  */
  double3_vertex4(wloc, wj[perm[j][0]], wj[perm[j][1]], wj[perm[j][2]]);
  for(l=0;l<3;l++) wjumpL[l]-=wloc[l];
  double3_vertex4(wloc, wj[perm[j][0]], wj[perm[j][2]], wj[perm[j][1]]);
  for(l=0;l<3;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, double3 Dleft, double3 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= -366;

  if(i == itest || j == itest || i == itest0 || j == itest0){
    for(l=0;l<3;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<3;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<3;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 \n", 
	   i,Dleft[0],Dleft[1],Dleft[2],nL[0],nR[0]);
    printf("***  Dright %d %7.5g %7.5g %7.5g %7.5g %7.5g \n", 
	   i,Dright[0],Dright[1],Dright[2],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<3;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<3;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<3;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<3;l++){
      Dleft[l]*=0.5;
      Dright[l]*=0.5;
    }

  }   /* end of boundary edge  */
  for(l=0;l<3;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 \n", 
	   i,Dleft[0],Dleft[1],Dleft[2],nL[0],nR[0]);
    printf("***  Dright %d %7.5g %7.5g %7.5g %7.5g %7.5g \n", 
	   i,Dright[0],Dright[1],Dright[2],nL[1],nR[1]);
  }

}

void Solver_Compute_Viscous_Volumes(tSolver_State *s){
  uint i,j,k,l;
  int ic;
  const int itest = -355, itest0=-366;
  double3 VI[s->DOF];  /* first component is always =0 */
  double3 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<3;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<3;j++)
         wx[j]= wy[j]=0;
      /* computing of \partial rho/x, rho/y   */
      for(j=0;j<2;j++) 
         nablarho[j] =0;
      for(j=0;j<s->DOF;j++){
        for(k=0;k<3;k++){
	  wx[k]+=dertestfunction[j][0]*s->element[i].w[j][k+1];
          wy[k]+=dertestfunction[j][1]*s->element[i].w[j][k+1];
        }
        for(k=0;k<2;k++) 
          nablarho[k]+=dertestfunction[j][k]*s->element[i].w[j][0];
      }
      for(k=0;k<3;k++){  /* dertestfunction is in fact \partial\phi *area  */
        wx[k]*=s->element[i].oneovervolume;
        wy[k]*=s->element[i].oneovervolume;
      }
      for(k=0;k<2;k++) nablarho[k]*=s->element[i].oneovervolume;
    }  /* 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<3;l++)
       wx[l]= wy[l]=0;
      /* computing of \partial rho/x, rho/y   */
      for(l=0;l<2;l++) 
        nablarho[l] =0;
        for(l=0;l<s->DOF;l++){
          for(k=0;k<3;k++){
            wx[k]+=s->curvedelement[ic].dertest[l][0][j]*
                   s->element[i].w[l][k+1];
            wy[k]+=s->curvedelement[ic].dertest[l][1][j]*
                   s->element[i].w[l][k+1];
          }
          for(k=0;k<2;k++) 
            nablarho[k]+=s->curvedelement[ic].dertest[l][k][j]*
                         s->element[i].w[l][0];
        }
        for(k=0;k<3;k++){  /* dertestfunction is in fact \partial\phi *area  */
          wx[k]*=2/s->curvedelement[ic].detJF[j];
          wy[k]*=2/s->curvedelement[ic].detJF[j];
        }
        for(k=0;k<2;k++) nablarho[k]*=2/s->curvedelement[ic].detJF[j];
      }
      ViscousTerm(s->element[i].w[j],s->element[i].w[j],
                  nablarho,wx,wy,s->Re,s->kappaRePr, Dx[j],Dy[j]);
    }
    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<3;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<3;l++)
	    VI[k][l]+=Dx[j][l]*dertestfunction[k][0] + 
	      Dy[j][l]*dertestfunction[k][1];
/*          if(i == itest ){
            l= 0;
            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 the value for edge integral terms */
    for(k=0;k<2;k++) s->element[i].nablarho[k] = nablarho[k]; 
    /* storing values for  boundary integral terms in VERTEXES of triangles*/
    for(j=0;j<s->DOF;j++){  
      double3_vertex(s->element[i].Dsxnablaw[j],Dx[perm[j][0]],
		     Dx[perm[j][2]],Dx[perm[j][1]]);
      double3_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<3;k++)
	  printf("Dsynablaw: %7.5g  %7.5g %7.5g %7.5g %7.5g %7.5g\n",
		 Dy[0][k],Dy[1][k],Dy[2][k],s->element[i].Dsynablaw[0][k],
		s->element[i].Dsynablaw[1][k],s->element[i].Dsynablaw[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<3;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 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){
 
          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[1], wL,
                 s->element[i].ve->P[perm[ii][1]]->Point[1], wR); 
      }
    Progress_Update(i+1);
  }
  Progress_Done();
  close(f1);
  }

int ViscousTerm(double4 w, double4 phi, 
                double2 nablarho, double3 wx, double3 wy, 
	        double Re, double kappaRePr, double3 D1, double3 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] = (2*wx[0]-wy[1] - 2*phi[1]*derrho[0]+phi[2]*derrho[1])*2*mirho/3;
  D1[1] = (wx[1]+wy[0] - phi[2]*derrho[0]-phi[1]*derrho[1])*mirho;
  D1[2] = (D1[0]*w[1]+D1[1]*w[2] + 
           kappaRePr*(wx[2]-(w[1]*wx[0]+w[2]*wx[1])/w[0]
           +( (phi[1]*w[1]+phi[2]*w[2])/w[0] -phi[3])*derrho[0]) )/w[0];
  D2[0] = D1[1];
  D2[1] = (2*wy[1]-wx[0] - 2*phi[2]*derrho[1]+phi[1]*derrho[0])*2*mirho/3;
  D2[2] = (D2[0]*w[1]+D2[1]*w[2] + 
           kappaRePr*(wy[2] -(w[1]*wy[0]+w[2]*wy[1])/w[0]
           +( (phi[1]*w[1]+phi[2]*w[2])/w[0]-phi[3])*derrho[1]))/w[0];
  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){
uint i,ii,j,k,l;
double const maxcosangle=-0.8;  /* angle between velocity nad normal vector */
double rho_inf=0,weight=0,cosangle;
double2 v_inf,v;
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;
for(k=0;k<s->nface;k++){
double3 wx,wy,Dx,Dy;
double2 nablarho;

  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;
    /* computing of \partial w/x, w/y   */
    for(j=0;j<3;j++) wx[j]=wy[j]=0;
    /* computing of \partial rho/x, rho/y   */
    for(j=0;j<2;j++) nablarho[j] =0;
    for(j=0;j<s->DOF;j++){
      for(l=0;l<3;l++){
	wx[l]+=s->element[i].dertest[j][0]*s->element[i].w[j][l+1];
	wy[l]+=s->element[i].dertest[j][1]*s->element[i].w[j][l+1];
      }
      for(l=0;l<2;l++) 
        nablarho[l]+=s->element[i].dertest[j][l]*s->element[i].w[j][0];
    }
    for(l=0;l<3;l++){  /* dertestfunction is in fact \partial\phi *area  */
      wx[l]*=s->element[i].oneovervolume;
      wy[l]*=s->element[i].oneovervolume;
    }
    for(l=0;l<2;l++) nablarho[l]*=s->element[i].oneovervolume;
    ViscousTerm(s->element[i].w[ii],s->element[i].w[ii],
                nablarho,wx,wy,s->Re,s->kappaRePr, Dx,Dy);
    coef[0]+= ((Dx[0]-Pressure(s->element[i].w[ii],s->kappa))*s->face[k].n[0]+
               Dy[0]*s->face[k].n[1])*s->face[k].volume;
    coef[1]+= (Dx[1]*s->face[k].n[0]+
              (Dy[1]-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);

}

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);

}
