

/*
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){

  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 {
    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] = p1->Point[0] -5;
      xp[1] = p1->Point[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;
    }
}

  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=25000; /* 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 || %7.5g %d %7.5g %d Curved \n",
	     se->P[0]->Point[0],se->P[0]->Point[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");
	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];
	  
	  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];

	/*	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 *), uint implic){
  tSolver_State *s;
  tMesh_VolumeElement *ve,**element;
  tMesh_SurfaceElement *se;
  tSolver_Curve_Element *cve;
  uint counter,i,j,k,l,ic,jpr;
  int itest = -355;
  #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;

  s->implicit_method=implic;

  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->integrate=integrateproc;
  s->t=0;
  s->nextt=1e20;  //infinity
  s->linfty=1e20;  //infinity
		       /*inicialization of derivation of test functions */
/*  
  s->dertestfunction[0][0] = 0;
  s->dertestfunction[0][1] = -2;
  s->dertestfunction[1][0] = 2;
  s->dertestfunction[1][1] = 2;
  s->dertestfunction[2][0] = -2;
  s->dertestfunction[2][1] = 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]); */
    }

    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;
	
	
        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];
	      /*	      if(i == 544) 
		printf("neigh-inn %d %d %d %i %i %d \n",
		       i,k,l,s->element[i].ve->P[k],
		       s->element[s->face[j].j].ve->P[l],
		       s->face[j].jj); */

	      s->face[j].boundary =0;
	    }
	  }
	  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 %i %i %i %d \n",
		       i,k,l,s->element[i].ve->P[k],
		       s->element[i].ve->P[k]->PeriodicCopy,
		       s->element[s->face[j].j].ve->P[l],
		       s->face[j].jj); */

	  }
	
	  /* 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);
  if(s->implicit_method) /* implicit method ... vertex base functions */
    Linear_Mapping_Computing_Vertex(s);
  else  /* explicit method  ... middle edge base function */
    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_Vertex(tSolver_State *s){
  tMesh_VolumeElement *ve;
  uint i,j,k,l,m;
  int ic,ib;
  static int itest =-1;
  double a[2];    /* VERTEX ASSOCIATED !!!!! */
  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_Vertex\n");
  for(i=0;i<s->nelem;i++){  /* in fact dertest times area(K_i)  */
    ve= s->element[i].ve;   /* VERTEX ASSOCIATED !!!!! */
    s->element[i].dertest[0][0] =-(ve->P[2]->Point[1] - ve->P[1]->Point[1])/2;
    s->element[i].dertest[0][1] = (ve->P[2]->Point[0] - ve->P[1]->Point[0])/2;
    s->element[i].dertest[1][0] =-(ve->P[0]->Point[1] - ve->P[2]->Point[1])/2;
    s->element[i].dertest[1][1] = (ve->P[0]->Point[0] - ve->P[2]->Point[0])/2;
    s->element[i].dertest[2][0] =-(ve->P[1]->Point[1] - ve->P[0]->Point[1])/2;
    s->element[i].dertest[2][1] = (ve->P[1]->Point[0] - ve->P[0]->Point[0])/2;

    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("Node %4d",i);
      for(k=0;k<3;k++)
	for(l=0;l<2;l++)
	  printf("%9.4g", s->element[i].dertest[k][l]);
      printf("\n");
    }
    
    /* 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]= -( 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])/2;
       e[j]= -(ve->P[perm[ib][1]]->Point[j] - ve->P[perm[ib][2]]->Point[j])/2;
    }
     /*
       printf("%7.5g  %7.5g curved \n",
       s->curvedelement[ic].xb[0],s->curvedelement[ic].xb[1]);
       printf("%7.5g  %7.5g ...... \n\n",
       0.5*(ve->P[perm[ib][0]]->Point[0]+ve->P[perm[ib][1]]->Point[0]),
       0.5*(ve->P[perm[ib][0]]->Point[1]+ve->P[perm[ib][1]]->Point[1])); 
     */

     /*     b[0] = 0;
	    b[1] = 0; */
                               /* VERTEX ASSOCIATED !!!!! */
    s->curvedelement[ic].dertest[perm[ib][0]][0][perm[ib][2]] = -b[1]-e[1];
    s->curvedelement[ic].dertest[perm[ib][0]][0][perm[ib][0]] = -e[1];
    s->curvedelement[ic].dertest[perm[ib][0]][0][perm[ib][1]] = -b[1]-e[1];
    s->curvedelement[ic].dertest[perm[ib][0]][1][perm[ib][2]] = b[0]+e[0];
    s->curvedelement[ic].dertest[perm[ib][0]][1][perm[ib][0]] = e[0];
    s->curvedelement[ic].dertest[perm[ib][0]][1][perm[ib][1]] = b[0]+e[0];

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

    s->curvedelement[ic].dertest[perm[ib][2]][0][perm[ib][2]] = e[1]-d[1];
    s->curvedelement[ic].dertest[perm[ib][2]][0][perm[ib][0]] = e[1]-d[1]-b[1];
    s->curvedelement[ic].dertest[perm[ib][2]][0][perm[ib][1]] = e[1]-d[1]+b[1];
    s->curvedelement[ic].dertest[perm[ib][2]][1][perm[ib][2]] = d[0]-e[0];
    s->curvedelement[ic].dertest[perm[ib][2]][1][perm[ib][0]] = d[0]-e[0]+b[0];
    s->curvedelement[ic].dertest[perm[ib][2]][1][perm[ib][1]] = d[0]-e[0]-b[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];

     if(i == itest)
       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("-----------------------");   */


    for(j=0;j<2;j++){ /* VERTEX ASSOCIATED !!!!! */
      b[j]*=-2;
      d[j]*=-2;
      e[j]*=-2;
    } 


    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(i == itest){ 
      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("Length: %d %d  %8.4g %8.4g %8.4g %8.4g\n",
	     i,ic,
	     s->curvedelement[ic].length[0],
	     s->curvedelement[ic].length[1],
	     s->curvedelement[ic].length[2],
	     s->curvedelement[ic].length[3]); 
      
      
    } 
    }
  }
  printf(" END of Linear_Mapping_Computing_Vertex\n");
}

void Linear_Mapping_Computing(tSolver_State *s){
  tMesh_VolumeElement *ve;
  uint i,j,k,l,m;
  int ic,ib;
  static int itest =-32;
  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_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;

  for(i=0;i<s->nelem;i++) 
    jumps[i]=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){
      uint j,k;
      double aver;
       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->number_of_limiting++;
    }
  }
}


void Solver_Limiting_Vertex(tSolver_State *s){
  uint i;
  double *jumps;
  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 = -673;

  jumps=calloc(s->nelem,sizeof(double));
  s->number_of_limiting=0;

  for(i=0;i<s->nelem;i++){ 
    jumps[i]=0; 
    s->element[i].limited =0; /* non limited */
  }
  
  for(i=0;i<s->nface;i++){
    if(s->face[i].boundary<=0){ /* inside or periodic */
      double jump,local1,local2;
      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;
      local1 = s->element[_i].w[perm[_ii][0]][k] - 
	s->element[_j].w[perm[_jj][1]][k];
      local2 = s->element[_i].w[perm[_ii][1]][k] - 
                s->element[_j].w[perm[_jj][0]][k];
      jump = pow2(local1) + local1*local2 + pow2(local2);
      jump*= s->face[i].volume/3;
      jumps[_i] +=jump;
      jumps[_j] +=jump;
    }
  }
  for(i=0;i<s->nelem;i++){
    double indicator;
    uint j,l;

    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); 
      printf("%d %7.5g %7.5g %7.5g \n",i,s->element[i].w[0][0],
	     s->element[i].w[1][0],s->element[i].w[2][0]);
      for(j=0;j<3;j++)
	if((l=s->element[i].neigh[j]) >=0)
	  printf("%d %7.5g %7.5g %7.5g \n",l,s->element[l].w[0][0],
		 s->element[l].w[1][0],s->element[l].w[2][0]);
    }

    if(indicator > s->limiting){
      int neigh,j;

      Solver_Limiting_Element(s,i);
      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->element[neigh].limited =1;
	  s->number_of_limiting++;
	}
      */
      s->number_of_limiting++;
    }
  }
  free(jumps);
}

void  Solver_Limiting_Element(tSolver_State *s, uint i){
  uint j,k;
  double aver;
  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].wold[j][k]=aver;
    }
  }

}
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(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=-96;
  double CFL;
  double alpha=0.2;  /* limitation of CFL numer */

  /*  printf("Begin of  Solver_Solve\n"); */

  CFL = s->CFL;
  s->tau=1E+02;
  s->Machmax = 0;
  s->Machmin = 100;

  //  if(s->init == 1 || s->init ==4){  /* s->init ==4 means Ringleb */
  if(s->init == 1 || s->init ==4){
    /*  initial computation, we restrict CFL  */
    CFL = CFL - (s->CFL-1)*exp(-alpha*s->t);
    /*printf("## %7.5g %7.5g %7.5g %d\n",s->CFL,CFL,s->t,s->init);  */
  }


  /* 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;
    } /* end of k */
    s->element[i].tauloc = CFL*s->element[i].tauloc;

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

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


  fflush(stdout); 

  /* time integration: \partial w/\partial t = F(w) */
  /* printf("Before integrating\n");  */

  /*  if(itest >=0) 
    for(k=0;k<3;k++)
      printf("### W^{k+1}:%7.5g %7.5g %7.5g %7.5g \n",
	     s->element[itest].w[k][0], s->element[itest].w[k][1],
	     s->element[itest].w[k][2], s->element[itest].w[k][3]);
  */

  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];
  
  s->integrate(s);


  /* 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(!s->implicit_method){ /* Only explicit methods */
      for(k=0;k<s->DOF;k++){
	/* local time stepping & global time stepping the same procedure */
	for(j=0;j<4;j++)
	  s->element[i].w[k][j]-=s->element[i].tauloc*s->element[i].H[k][j];
      }
    }
    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(i == itest){ 
	printf("@@@ W^{k+1}:%7.5g %7.5g %7.5g %7.5g %7.5g\n",
	       s->element[i].w[k][0], s->element[i].w[k][1],
	       s->element[i].w[k][2], s->element[i].w[k][3],pp);
	printf("... H^{k+1}:%7.5g %7.5g %7.5g %7.5g %7.5g\n",
	       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].tauloc);
      }
      if(pp <=0){
	/*	double4 w_old,H_old;
	double tau_new,p_old,pp_old;

	pp_old = pp;
	for(j=0;j<4;j++){ 
	  H_old[j] = s->element[i].H[k][j];
	  w_old[j] = s->element[i].w[k][j] + H_old[j]*s->element[i].tauloc;
	}
	p_old=(s->kappa-1)*(w_old[3]-0.5*
			    (w_old[1]*w_old[1]+w_old[2]*w_old[2])/w_old[0]);
	tau_new = (p_old-0.01)/15 /((s->kappa-1)*(H_old[3]-0.5*
			    (w_old[1]*H_old[1]+w_old[2]*H_old[2])/w_old[0]));

	for(j=0;j<4;j++)
	  s->element[i].w[k][j]=w_old[j]-tau_new*H_old[j];

	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);
	
	printf("!!!!  Time step on triangle %d(%d) is %7.5g(old) %7.5g(new) \n",
	       i,k,s->element[i].tauloc,tau_new);
	printf("%7.5g %7.5g w^{k}  :%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],
	       w_old[0],w_old[1],w_old[2],w_old[3],p_old,pp_old);
	printf("%7.5g %7.5g H^{k}  :%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],
	       H_old[0],H_old[1],H_old[2],H_old[3]);
	printf("%7.5g %7.5g w^{k+1}:%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].w[k][0], s->element[i].w[k][1],
	       s->element[i].w[k][2], s->element[i].w[k][3],pp);
	*/

	printf("%7.5g %7.5g p=%7.5g \n",
	       s->element[i].ve->P[0]->Point[0],
	       s->element[i].ve->P[0]->Point[1],pp);
	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);

      /*      if(Mach > 2.0 && k==1){
	double xbar[2],v,c;
	barycentre(s,i,xbar);
	printf("%7.5g %7.5g Mach (%d,%d) %7.5g %7.5g %7.5g %7.5g %7.5g \n",
	       xbar[0],xbar[1],i,k,Mach,s->element[i].tauloc,
	       sqrt(rv2/s->element[i].w[k][0]),
	       pp,sqrt(s->kappa*pp/s->element[i].w[k][0]));
	       }*/
    }  /* end of k */
    /*      if(i == itest || i == itest0)Show_Inviscid_Term(s,i); */

  }  /* end of i */

  s->t+=s->tau;

  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 > 1E-05)  */
      /*      err/ s->element[i].oneovervolume,s->linfty); */

  }
  s->linfty/=s->DOF;
}

void Solver_MeshUpdate(tSolver_State *s){
  uint i,j,k;
  double2 coef;
  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){
    Coeficients(s,coef);
    printf("Drag coeficient c_D=%7.5g, lift coeficient c_L=%7.5g \n",
    coef[0],coef[1]); 
  }
}

