/*
What's new:
C44 - [Martin Zoubek] first public release
*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "FileIO.h"
#include "Memory.h"
#include "Mesh.h"
#include "Solver.h"
#include "Time.h"
#include "Types.h"

int main(int argc,char **argv){
  FILE *f0,*f1,*f2,*f3;
  char *quantity,*gnufile,*stype;
  tMesh_VolumeElement *ve;
  char sf[7];
  tMesh_Mesh *m;
  double qmin,qmax,CFL;
  uint nisol,ignu;
  uint numwall;  /* number of par of bundary to draw */
  uint i,j;
  BOOL imaxmin;  /* =1  => the min and max values for plot are given by user */
  tSolver_Boundary *b;
  void (*integrateproc)(tSolver_State *)=0;
  tSolver_State *s;

  #define WriteTimeStep	300
  #define Exit0(s,c){\
    printf(s);\
    return c;\
    }
  #define Exit1(s,s1,c){\
    printf(s,s1);\
    return c;\
    }
  #define eq(a,b) (fabs((a)-(b))<1e-20)


  if(argc<3)
    Exit0("Syntax: Solve inifile\n",1);

  if((f1=fopen(argv[1],"rb"))==0)
    Exit1("Error: Can't open file '%s'\n",argv[1],2);
  if((f2=fopen(argv[2],"rb"))==0)
    Exit1("Error: Can't open file '%s'\n",argv[2],2);
  m=Mesh_Read(f1,f2);  
  fclose(f1);
  fclose(f2);

  m->kappa = 1.4;
  Mesh_FindNeighbours(m);

  /*  printf("Before  Curved_Element_Detecting \n"); */
  if(!Curved_Element_Detecting(m))
         Exit0("Error: No file with curved boundary describtion\n",3);
  /*  printf("After Curved_Element_Detecting \n"); */

  s=Solver_Init(m,b,CFL,integrateproc);
  printf("Solver_Init finished\n");


  if((f0=fopen("plot.dgm","rb"))==0)
    Exit1("Error: Can't open file '%s'\n","plot.dgm",2);

  /*  if((f4=fopen("gnu-files","rb"))==0)
    Exit1("Error: Can't open file '%s'\n","gnu-files",2);*/

 
  sf[0] = 'g';
  sf[1] = 'n';
  sf[2] = 'u';
  sf[3] = '.';
  sf[6] = '\0';

  quantity=FileIO_ReadString(f0);  /* quantity to draw */
  ignu = 1;

  while(strcmp(quantity,"EOD")){
    if(ignu < 10){
      sf[4] = '0';
      sf[5] = ignu+'0';
    }
    else if (ignu >=10 && ignu <100){
      sf[4] = (ignu - ignu%10)/10 +'0';
      sf[5] = ignu%10 + '0';
    }

    /*    gnufile=FileIO_ReadString(f4);*/
    gnufile = sf;
    f3=fopen(gnufile,"wb");

    if(!strcmp(quantity,"E2V")){ 
      printf("%s: Edge2Volume ",gnufile);
      Edge2Volume(f3,m);
      printf("\n");
    }
    else{
      stype=FileIO_ReadString(f0);
      if(!strcmp(stype,"ISO")) {
	printf("%s: Isolines: ",gnufile);
	nisol = FileIO_ReadInteger(f0);
	imaxmin = FileIO_ReadInteger(f0);
	if(imaxmin){
	  qmin= FileIO_ReadFloat(f0);
	  qmax= FileIO_ReadFloat(f0);
	}
	else{
	  qmin = 1E+10;
	  qmax = 1E-10;
	}
	/* Relimiting(m); */
	PlotIsolines(f3,m,s,quantity,nisol,imaxmin,qmin,qmax);
      }
      else if(!strcmp(stype,"WALL")){ 
	printf("%s: Wall plot: ",gnufile);
	numwall = FileIO_ReadInteger(f0);
	for (i=0;i<numwall;i++)
	  PlotWall(f3,m,s,quantity,FileIO_ReadInteger(f0));
	printf("\n");
      }
      else if(!strcmp(stype,"CUT")){ 
	double2 x1,x2;
	printf("%s: Diag cut  ",gnufile);
	x1[0] = FileIO_ReadFloat(f0);
	x1[1] = FileIO_ReadFloat(f0);
	x2[0] = FileIO_ReadFloat(f0);
	x2[1] = FileIO_ReadFloat(f0);
	PlotCut(f3,m,s,quantity,x1,x2);
	printf("\n");
      }
      else
	printf("Any other subroutine than ISO\n");
    }
    fclose(f3);
    quantity=FileIO_ReadString(f0); /* quantity to draw, EOD => end of data */
    ignu++;
  }
  fclose(f0);
}

void Relimiting(tMesh_Mesh *m){
  tMesh_VolumeElement *ve, *ve1;
  double jumps[m->nelem][3];
  uint indexjump[m->nelem][3];
  tMesh_VolumeElement *vejump[m->nelem][3];
  double local,jump;
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};
  const uint k=0;  /* limiting in density */
  uint i,j,l;

  for(i=0;i<m->nelem;i++)
    for(j=0;j<3;j++)
      jumps[i][j]=0;

  /* computing of jumps */
  i=0;
  ve=m->VolumeElement;
  while(ve){
    for(j=0;j<3;j++){
      double local,jump,lenght;
      double2 Ploc;
      if(ve1 = ve->Neighbour[j])
	for(l=0;l<3;l++)
	    if( (ve->P[perm[j][0]] == ve1->P[perm[l][1]]) &&
		(ve->P[perm[j][1]] == ve1->P[perm[l][0]]) ){
	      vejump[i][j] = ve1;
	      indexjump[i][j]=l;
	      local = ve->w[perm[j][0]][k] - ve1->w[perm[l][0]][k];
	      jump = pow2(local);
	      local = ve->w[perm[j][2]][k] - ve->w[perm[l][1]][k] +
		ve1->w[perm[j][2]][k] - ve1->w[perm[l][1]][k];
	      jump+= pow2(local)/3;
	      jumps[i][j] =jump;
	      break;
	    }
    } 
    ve=ve->Next;
    i++;
  }

/* recomputing the values of w */
  i=0;
  ve=m->VolumeElement;
  while(ve){
    double area,indicator;
    area = Triangle_Volume(ve->P[0]->Point,ve->P[1]->Point,ve->P[2]->Point);
    area = pow_fast(area,0.75);
    for(j=0;j<3;j++){
      indicator = jumps[i][j]/pow_fast(area,0.75);
      if(indicator > 0.5E-01) {
	/* printf("%7.5g %7.5g %7.5g Neighbours: %d %d  %7.5g \n",
	       ve->P[perm[j][0]]->Point[0],ve->P[perm[j][0]]->Point[1],
	       ve->w[j][0],	       i,j,indicator);
	printf("%7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n \n",
	       ve->P[perm[j][1]]->Point[0],ve->P[perm[j][1]]->Point[1],
	       vejump[i][j]->w[indexjump[i][j]][0],
	       vejump[i][j]->w[0][0],vejump[i][j]->w[1][0],
	       vejump[i][j]->w[2][0]);  */

	for(l=0;l<4;l++){
	  ve->w[j][l]= (ve->w[j][l]+vejump[i][j]->w[indexjump[i][j]][l])/2;
	  vejump[i][j]->w[indexjump[i][j]][l] = ve->w[j][l];
	}
      }
      
    } 
    ve=ve->Next;
    i++;
  }
} 

void PlotCut(FILE *gnu, tMesh_Mesh *m, tSolver_State *s, 
	     char *quantity, double2 x1, double2 x2){
  tMesh_VolumeElement *ve;
  uint i,j,num_inter;
  double2 x0, x[3], xint[3];
  double3 q;
  uint inters;
  double t0, t1;
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};

  CalculateQuantity( m, s, quantity);

  ve=m->VolumeElement; 
  while(ve){ 

    /* courent element */

    for(j=0;j<3;j++){
      x[j][0] = ve->P[j]->Point[0]; 
      x[j][1] = ve->P[j]->Point[1];
      q[j] = ve->q[perm[j][0]] + ve->q[perm[j][2]] - ve->q[perm[j][1]];
    }

    /* we seek intersection of triangle with the cut line */
    num_inter = 0;
    for(j=0;j<3;j++){
      inters = Intersection( x1, x2, x[perm[j][0]], x[perm[j][1]], x0);
      // printf("After %7.5g %7.5g     inters=%d\n", x0[0],x0[1],inters);
      if(inters==1){
	xint[num_inter][0] = x0[0];
	xint[num_inter][1] = x0[1];
	num_inter++;
      }
    }

    if(num_inter==3){//intersection goes through the triangle into the vertex
      if(double2_distance2(xint[1],xint[2])>1E-10){
      }
      else{
	xint[0][0] =xint[2][0];
	xint[0][1] =xint[2][1];
      }
    }
    else if(num_inter==1){//intersection goes out of the 
      //                    triangle into the vertex
      xint[1][0] =xint[0][0];
      xint[1][1] =xint[0][1];
    }
    if(num_inter>0){ //plot the intersection between  xint[0] and  xint[1]
      double q0, q1;

      
      q0=Interpolate(x[0],x[1],x[2],xint[0],q);
      q1=Interpolate(x[0],x[1],x[2],xint[1],q);
 
      fprintf(gnu," %7.5g %7.5g %7.5g  \n",xint[0][0],q0,xint[0][1]);
      fprintf(gnu," %7.5g %7.5g %7.5g  \n\n",xint[1][0],q1,xint[1][1]);
    }

    ve = ve->Next;
  }
}
double Interpolate(double2 x0,double2 x1, double2 x2, double2 xint,
		   double3 q){
  
  double det[4], q0=0;
  int j;
  
  det[3] = x0[0]*(x1[1]-x2[1])+ x1[0]*(x2[1]-x0[1])+ x2[0]*(x0[1]-x1[1]);
  det[0] = xint[0]*(x1[1]-x2[1])+ x1[0]*(x2[1]-xint[1])+x2[0]*(xint[1]-x1[1]);

  det[1] = x0[0]*(xint[1]-x2[1])+ xint[0]*(x2[1]-x0[1])+x2[0]*(x0[1]-xint[1]);
  det[2] = x0[0]*(x1[1]-xint[1])+ x1[0]*(xint[1]-x0[1])+xint[0]*(x0[1]-x1[1]);

  for(j=0;j<3;j++)
    q0 += q[j]*det[j]/det[3];

  return q0;
}


int Intersection( double2 x1,double2 x2,double2 x3,double2 x4,double2 x0){
  double  eps = 0.001, st, rlen1, rlen2, rtol,t1,t2;
  uint inters =0;

  st = (x2[1]-x1[1])*(x4[0]-x3[0])+(x1[0]-x2[0])*(x4[1]-x3[1]);
  if(st != 0){
    /* nonparalel edges */
    
    t1 = (x2[0]*(x3[1]-x1[1])+ x1[0]*(x2[1]-x3[1])+ x3[0]*(x1[1]-x2[1]))/st;
    t2 = (x1[1]*(x3[0]-x4[0])+ x3[1]*(x4[0]-x1[0])+ x4[1]*(x1[0]-x3[0]))/st;
    
    /*
    rlen1 = pow_fast((x1[0]-x2[0])*(x1[0]-x2[0]) + 
		     (x1[1]-x2[1])*(x1[1]-x2[1]),0.5);
    rlen2 = pow_fast((x3[0]-x4[0])*(x3[0]-x4[0]) + 
		     (x3[1]-x4[1])*(x3[1]-x4[1]),0.5);
    
    //   prescribed tolerance 
    rtol = eps*min(rlen1, rlen2);
    */
    
      /* intersection */
    if( t1 >= 0 &&  t1 <= 1 && t2 >= 0 && t2 <= 1 ) inters = 1;
    x0[0] = x1[0] + t2*(x2[0]-x1[0]);
    x0[1] = x1[1] + t2*(x2[1]-x1[1]);
 
    /*
    printf("\n---------\nINSC  %7.5g %7.5g  %7.5g %7.5g inters=%d \n\n", 
	   x0[0],x0[1],t1,t2,inters);
    */    
	
    /*     x0[0] = x3[0] + t1*(x4[0]-x3[0]);
	   y0[0] = x3[1] + t1*(x4[1]y-x3[1]);
    */
  }
  return inters;
}


void PlotWall(FILE *gnu,tMesh_Mesh *m,  tSolver_State *s,
	      char *quantity, int iwall){
  tMesh_SurfaceElement *se;
  tMesh_VolumeElement *ve;
  uint i,j;
  uint ic;  /* ic =0, horizontal,  ic=1  vertical */
  double totallength,progreslength;
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};

  CalculateQuantity( m, s, quantity);

  totallength = 0;
  se=m->SurfaceElement; 
  while(se){
    double2 PQ;
    if(se->surfnr == inabs(iwall)){
      double2_sub(PQ,se->P[0]->Point,se->P[1]->Point);
      totallength += sqrt(double2_norm2(PQ));
    }
    se=se->Next;
  }
  printf( "%s (nbc=%d) ,",quantity, iwall);

  if(iwall == 1 || iwall == 2)
    ic = 1;
  else
    ic = 0;

  progreslength=0;
  ve=m->VolumeElement; 
  while(ve){
    for(i=0;i<3;i++){
      if(!(ve->Neighbour[i]) &&  ve->Surface[i]->surfnr==inabs(iwall)){
	double2 PQ;
	double length,q;
	double2_sub(PQ,ve->P[perm[i][0]]->Point,ve->P[perm[i][1]]->Point);
	length = sqrt(double2_norm2(PQ));

	q=ve->q[perm[i][0]] + ve->q[perm[i][2]] - ve->q[perm[i][1]];
	fprintf(gnu," %7.5g %7.5g  \n",ve->P[perm[i][0]]->Point[ic],q);
	/*		progreslength/totallength,q); */
	q=ve->q[perm[i][0]] + ve->q[perm[i][1]] - ve->q[perm[i][2]];
	fprintf(gnu," %7.5g %7.5g  \n\n",ve->P[perm[i][1]]->Point[ic],q);
	/*(progreslength+length)/totallength,q);*/

	progreslength+=length;

      }
    }
    ve=ve->Next;
  }

}


void CalculateQuantity(tMesh_Mesh *m,  tSolver_State *s, char *quantity){
  tMesh_VolumeElement *ve;
  uint j;
  uint i,k;
  double rho_inf=0,p_inf=0;
  double2 v_inf;
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};
  static uint perm2[2][2]={{0,1},{1,0}};
  double a[2];

  // Blasius problem  
  /*
  double Re=1.E+04;
  double xinput = -0.45;
  */

  // NACA 0012
  
  //double Re=73;
  //double Re=5000;
    double Re=500;
  double xinput = -2;
  

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


  if(!strcmp(quantity,"PC") || !strcmp(quantity,"CF") ){
    double weight=0,cosangle,rho,rv2,p;
    double2 v,xc;
     v_inf[0]=v_inf[1] = 0;
     ve=m->VolumeElement;
     while(ve){
       for(j=0;j<3;j++){
	 if(!ve->Neighbour[j]){
	   xc[0]=(ve->P[0]->Point[0]+ve->P[1]->Point[0]+ve->P[2]->Point[0])/3;
	   xc[1]=(ve->P[0]->Point[1]+ve->P[1]->Point[1]+ve->P[2]->Point[1])/3;
	   if(xc[0] < xinput){
	     v[0]=v[1]=rho=p=0;
	     for(k=0;k<3;k++){
	       for(j=0;j<2;j++)
		 v[j]+=ve->w[k][j+1]/ve->w[k][0];
	       rho+=ve->w[k][0];
	       rv2 = (ve->w[k][1]*ve->w[k][1]+ve->w[k][2]*ve->w[k][2])/
		 ve->w[k][0];
	       p+= (m->kappa-1)*(ve->w[k][3]-0.5*rv2);
	     }
	     for(j=0;j<2;j++) v[j]/=3;
	     rho/=3;
	     p/=3;
	     rho_inf+=rho;
	     for(j=0;j<2;j++) v_inf[j]+=v[j];
	     p_inf+=p;
	     weight+=1;
	     /*	     printf("%7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",
		    ve->P[1]->Point[0],ve->P[1]->Point[1],
		    rho,v[0],v[1],p); */
	   }
	 }
       }
       ve=ve->Next;
     }
     rho_inf/=weight;
     for(j=0;j<2;j++) v_inf[j]/=weight;
     p_inf/=weight;
     printf("\n INFINITY rho=%7.5g, v=(%7.5g, %7.5g),|v|=%7.5g, p=%7.5g \n",
	    rho_inf,v_inf[0],v_inf[1],pow_fast(double2_norm2(v_inf),0.5), p);
  }

  if(!strcmp(quantity,"CF")) { /* skin friction coefficient */
   double4 Dx,Dy;

    s->Re = Re;
    s->kappa = m->kappa;
    s->kappaRePr = m->kappa/Re/0.72;

    Solver_Compute_Velocity_Deformation(s);   /* in order to compute Dx,Dy */

    //    i=m->nelem;
    //ve=m->VolumeElement; 
    //while(ve){
    for(i=0;i<s->nelem;i++){
      double p,rv2, rv2a, det;
      double3 x,y,u,v;
      double2 n,t;
      double dij[2][2], tau[2][2];
      uint k,l;
      int itest = -19860;

      ve=s->element[i].ve; 
	
      for(k=0;k<3;k++){
	x[k] = (s->element[i].ve->P[perm[k][0]]->Point[0]+s->element[i].ve->P[perm[k][1]]->Point[0])/2;
	y[k] = (s->element[i].ve->P[perm[k][0]]->Point[1]+s->element[i].ve->P[perm[k][1]]->Point[1])/2;
	u[k] = s->element[i].w[k][1]/s->element[i].w[k][0];
	v[k] = s->element[i].w[k][2]/s->element[i].w[k][0];
      }

      // deformation of velocity 
      det = (x[1]-x[0])*(y[2]-y[0]) - (x[2]-x[0])*(y[1]-y[0]);
      dij[0][0] = ((u[1]-u[0])*(y[2]-y[0]) - (u[2]-u[0])*(y[1]-y[0]))/det;
      dij[0][1] = ((x[1]-x[0])*(u[2]-u[0]) - (x[2]-x[0])*(u[1]-u[0]))/det;
      dij[1][0] = ((v[1]-v[0])*(y[2]-y[0]) - (v[2]-v[0])*(y[1]-y[0]))/det;
      dij[1][1] = ((x[1]-x[0])*(v[2]-v[0]) - (x[2]-x[0])*(v[1]-v[0]))/det;
	
      det = 2./3.*(dij[0][0] + dij[1][1] );  // div v
      for(k=0;k<2;k++){
	for(l=0;l<2;l++){
	  tau[k][l] = dij[k][l] + dij[l][k];
	}
	tau[k][k] -= det;
      }

      if(i == itest){
	//if(s->element[i].icve >= 0){

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


	printf("\n%7.5g %7.5g %7.5g %7.5g %d \n",x[0],y[0],u[0],v[0],i);
	printf("%7.5g %7.5g %7.5g %7.5g \n",x[1],y[1],u[1],v[1]);
	printf("%7.5g %7.5g %7.5g %7.5g \n",x[2],y[2],u[2],v[2]);
	printf("%7.5g %7.5g \n\n",x[0],y[0]);
	printf("i: %i,dij: %7.5g %7.5g %7.5g %7.5g, tau: %7.5g %7.5g %7.5g  %7.5g\n",i,
	       dij[0][0],dij[0][1],dij[1][0],dij[1][1],
	       tau[0][0],tau[0][1],tau[1][0],tau[1][1]);
      }
      
      rv2 = 0;  // default for nonboundary elements 
      rv2a = 0;  // default for nonboundary elements 
      
      for(k=0;k<3;k++){
	//if(s->element[i].neigh[k]<0){  // boundary element
	  if(!ve->Neighbour[k]){  // boundary element
	  double rlen,rlen1;
	  int ic, ii;

	  ii = k;
	  for(j=0;j<2;j++){
	    ic=s->element[i].icve;	
	    //if(ic >= 0 && ic < 0){  /* curved elements */
	    if(ic >= 0){  /* curved elements */

	      double4_weight_awer(Dx,s->element[i].Dsxnablaw[perm[ii][0]],
				  s->element[i].Dsxnablaw[perm[ii][1]],
				  a[perm2[0][j]],a[perm2[1][j]]);
	      double4_weight_awer(Dy,s->element[i].Dsynablaw[perm[ii][0]],
				  s->element[i].Dsynablaw[perm[ii][1]],
				  a[perm2[0][j]],a[perm2[1][j]]);

	      if(i == itest){
		printf("%7.5g %7.5g %d curved \n",
		       s->element[i].ve->P[perm[k][1]]->Point[0],
		       s->element[i].ve->P[perm[k][1]]->Point[1],k);
		printf("%7.5g %7.5g %d \n",
		       s->element[i].ve->P[perm[k][0]]->Point[0],
		       s->element[i].ve->P[perm[k][0]]->Point[1],k);
	      }

	      n[0] = s->curvedelement[ic].normal[2+j][0];
	      n[1] = s->curvedelement[ic].normal[2+j][1];
	      t[0] = -n[1];
	      t[1] = n[0];
	    }
	    else{ /* nonucurved elements */

	      if(i == itest){
		printf("%7.5g %7.5g %d NONcurved \n",
		       s->element[i].ve->P[perm[k][1]]->Point[0],
		       s->element[i].ve->P[perm[k][1]]->Point[1],k);
		printf("%7.5g %7.5g %d \n",
		       s->element[i].ve->P[perm[k][0]]->Point[0],
		       s->element[i].ve->P[perm[k][0]]->Point[1],k);
	      }


	      double4_average(Dx,s->element[i].Dsxnablaw[perm[ii][0]],
			      s->element[i].Dsxnablaw[perm[ii][1]]);
	      double4_average(Dy,s->element[i].Dsynablaw[perm[ii][0]],
			      s->element[i].Dsynablaw[perm[ii][1]]); 
       	  
	      t[0] = s->element[i].ve->P[perm[k][1]]->Point[0] - s->element[i].ve->P[perm[k][0]]->Point[0];
	      t[1] = s->element[i].ve->P[perm[k][1]]->Point[1] - s->element[i].ve->P[perm[k][0]]->Point[1];

	      rlen = pow_fast(double2_norm2(t),0.5);
	      t[0]/=rlen;
	      t[1]/=rlen;
	      n[0] = t[1];
	      n[1] = -t[0];
	    }


	    rv2 = t[0]*(tau[0][0]*n[0] + tau[0][1]*n[1]) +
	      t[1]*(tau[1][0]*n[0] + tau[1][1]*n[1]);
	    
	    rv2a = t[0]*(Dx[1]*n[0] + Dy[1]*n[1]) +
	      t[1]*(Dx[2]*n[0] + Dy[2]*n[1]);

	    if(i == itest){
	      printf("Old:      %d %7.5g %7.5g %7.5g %7.5g \n",
		     i,tau[0][0],tau[0][1],tau[1][0],tau[1][1]);
	      printf("New: ver1 %d %7.5g %7.5g %7.5g %7.5g \n",
		     i,s->element[i].Dsxnablaw[perm[ii][0]][1],
		     s->element[i].Dsynablaw[perm[ii][0]][1],
		     s->element[i].Dsxnablaw[perm[ii][0]][2],
		     s->element[i].Dsynablaw[perm[ii][0]][2]);
	      printf("New: ver2 %d %7.5g %7.5g %7.5g %7.5g \n",
		     i,s->element[i].Dsxnablaw[perm[ii][1]][1],
		     s->element[i].Dsynablaw[perm[ii][1]][1],
		     s->element[i].Dsxnablaw[perm[ii][1]][2],
		     s->element[i].Dsynablaw[perm[ii][1]][2]);

	      printf("   : nod %d %7.5g %7.5g %7.5g %7.5g \n---------------\n",
		     i,Dx[1],Dy[1],Dx[2],Dy[2]);
	      printf("   a( , ) = %7.5g %7.5g \n----------------\n",a[0],a[1]);
	      printf("n( , )= %7.5g %7.5g, t( , )= %7.5g %7.5g \n",
		    n[0],n[1],t[0],t[1] );

	    }


    /*	    printf("%7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g\n",
		    ve->P[perm[k][0]]->Point[0],ve->P[perm[k][0]]->Point[1],
		   t[0],t[1],rlen, rho_inf, double2_norm2(v_inf),Re);
	    printf("%7.5g %7.5g %7.5g %7.5g \n",
		   ve->P[perm[k][1]]->Point[0],ve->P[perm[k][1]]->Point[1],
		   n[0],n[1]);
	    printf("%7.5g %7.5g,  tau=%7.5g\n\n",
		   ve->P[perm[k][2]]->Point[0],ve->P[perm[k][2]]->Point[1],
		   rv2);
	    */
	    //ve->q[perm[k][j]] = -rv2*2/rho_inf*double2_norm2(v_inf)/Re;
	    ve->q[perm[k][j]] = -rv2a*2/rho_inf*double2_norm2(v_inf)/Re;
	  }
	  ve->q[perm[k][2]] = ( ve->q[perm[k][0]] +  ve->q[perm[k][1]])/2;
	  

	  /* averaging for smoothing */
	  ve->q[perm[k][0]] = ve->q[perm[k][2]];
	  ve->q[perm[k][1]] = ve->q[perm[k][2]];

	}
      }
      //      i--;
      // ve=ve->Next;
    }
    return;
  }
  /* other quantities then CF */
  /* setting of values to compute */
  i=m->nelem;
  ve=m->VolumeElement; 
  while(ve){
    double p,rv2, det;
    //    double3 x,y,u,v;
    //double2 n,t;
    //double dij[2][2], tau[2][2];
    uint k,l;

    for(j=0;j<3;j++){
      if(!strcmp(quantity,"RO")){ /*density*/
	ve->q[j]=ve->w[j][0];  
      }
      else if(!strcmp(quantity,"RV")) { /* absolute value of velocity */
 	rv2 = (ve->w[j][1]*ve->w[j][1] + ve->w[j][2]*ve->w[j][2]);
	ve->q[j] = sqrt(rv2);
      }
      else if(!strcmp(quantity,"V")) { /* absolute value of velocity */
 	rv2 = (ve->w[j][1]*ve->w[j][1] + ve->w[j][2]*ve->w[j][2]);
	ve->q[j] = sqrt(rv2)/ve->w[j][0];
      }
      else if(!strcmp(quantity,"V1")) { /* first component of velocity */
	ve->q[j] = ve->w[j][1]/ve->w[j][0];
      }
      else if(!strcmp(quantity,"V2")) { /*  second component of velocity */
	ve->q[j] = ve->w[j][2]/ve->w[j][0];
      }
      else if(!strcmp(quantity,"P")) { /* pressure */
 	rv2 = (ve->w[j][1]*ve->w[j][1] + ve->w[j][2]*ve->w[j][2])/ve->w[j][0];
	ve->q[j] = (m->kappa-1)*(ve->w[j][3]-0.5*rv2);
      }
      else if(!strcmp(quantity,"PC")) { /* pressure coefficient */
 	rv2 = (ve->w[j][1]*ve->w[j][1] + ve->w[j][2]*ve->w[j][2])/ve->w[j][0];
	p = (m->kappa-1)*(ve->w[j][3]-0.5*rv2);
	ve->q[j] = (p-p_inf)*2/rho_inf*double2_norm2(v_inf);
      }
      else if(!strcmp(quantity,"M")) {   /* Mach number */
	rv2 = (ve->w[j][1]*ve->w[j][1] + ve->w[j][2]*ve->w[j][2])/ve->w[j][0];
	p = (m->kappa-1)*(ve->w[j][3]-0.5*rv2);
	ve->q[j]=sqrt(rv2/p/m->kappa);
	if(ve->q[j] < 0)
	  printf("Mach negative %i %d %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",
		 i,j,ve->w[j][0],ve->w[j][1],ve->w[j][2],rv2,p,ve->q[j]);
	/*	if(ve->q[j] > 2)
	  printf("Mach too large %i %d %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",
		 i,j,ve->w[j][0],ve->w[j][1],ve->w[j][2],rv2,p,ve->q[j]);
	*/
      }
      else if(!strcmp(quantity,"E")) {   /* total energy */
	ve->q[j]=ve->w[j][3];  
      }
      else if(!strcmp(quantity,"T")) {   /* absolute temperature */
	rv2 = (ve->w[j][1]*ve->w[j][1] + ve->w[j][2]*ve->w[j][2])/ve->w[j][0];
	ve->q[j]= (ve->w[j][3]-0.5*rv2)/ve->w[j][0];
      }
      else
	printf("Bad quantity to draw: %s\n",quantity);
    }
    i--;
    ve=ve->Next;
  }

}

void Edge2Volume(FILE *gnu,tMesh_Mesh *m){
  tMesh_VolumeElement *ve;
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};
  uint j,k;
  double q;

  ve=m->VolumeElement; 
  while(ve->Next){ 
    ve=ve->Next;
  }
  while(ve != m->VolumeElement ){ /*reverse order */
    for(j=0;j<3;j++){
      for(k=0;k<4;k++){
	q=ve->w[perm[j][0]][k] + ve->w[perm[j][2]][k] -
	  ve->w[perm[j][1]][k];
	fprintf(gnu,"%12.10g ",q);
      }
      fprintf(gnu,"\n");
    }
    ve=ve->Prev;
  }

  ve = m->VolumeElement;
  for(j=0;j<3;j++){
    for(k=0;k<4;k++){
      q=ve->w[perm[j][0]][k] + ve->w[perm[j][2]][k] -
	ve->w[perm[j][1]][k];
      fprintf(gnu,"%12.10g ",q);
    }
    fprintf(gnu,"\n");
  }

}

void PlotIsolines(FILE *gnu,tMesh_Mesh *m,  tSolver_State *s,
		  char *quantity, uint nisol, 
		  BOOL imaxmin, double qmin, double qmax){
  tMesh_VolumeElement *ve;
  tMesh_SurfaceElement *se;
  const double eps=1E-05;
  uint i,j,j1,k;
  int level,lev_per,istart,iend;
  double qisol[nisol];
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};
  uint counter=0;

  CalculateQuantity( m, s, quantity);

  /* seeking of qmin and qmax */
  if(!imaxmin){
    ve=m->VolumeElement; 
    while(ve){ 
      double q;
      for(j=0;j<3;j++){
	q=ve->q[perm[j][0]] + ve->q[perm[j][2]] - ve->q[perm[j][1]];
	qmin=min(qmin,q);
	qmax=max(qmax,q);
      }
      ve=ve->Next;
    }
  }

  printf(" %s, value interval ( %7.5g, %7.5g) \n",quantity,qmin,qmax);


  for(k=0;k<nisol;k++)
    qisol[k] = qmin+ (qmax-qmin)*k/(nisol-1);

  /* usualy only one times, two times for square */
  /*  for(lev_per=0;lev_per<m->nper;lev_per++){  */
  
  lev_per = 0;

  if(m->nper > 0) {
    istart = -1;
    iend = 1;
  }
  else{
    istart = 0;
    iend =0;
  }
   /* drawing period with -1,0,-1 times */
     for(level=istart;level<=iend;level++){ 
      ve=m->VolumeElement; 
      while(ve){ 
	double q[3],x[3][2],xp[2];
	BOOL iplot;
	
	counter++;
	
	for(j=0;j<3;j++){
	  x[j][0] = ve->P[j]->Point[0]+level*m->rper[2*lev_per];
	  x[j][1] = ve->P[j]->Point[1]+level*m->rper[2*lev_per+1];
	  q[j] = ve->q[perm[j][0]] + ve->q[perm[j][2]] - ve->q[perm[j][1]];
	}
	
	for(k=0;k<nisol;k++){
	  iplot = 0;
	  for(j=0;j<3;j++){
	    j1 = perm[j][1];
	    
	    if(fabs(q[j1] - q[j]) < eps){
	      if(fabs(q[j]-qisol[k]) < eps){
		fprintf(gnu," %7.5g %7.5g %6.4g \n",x[j][0],x[j][1],qisol[k]);
		fprintf(gnu," %7.5g %7.5g %6.4g \n",x[j1][0],x[j1][1],qisol[k]);
		iplot = 1;
	      }
	    }
	    else{
	      if(( q[j] < qisol[k] && qisol[k] <= q[j1]) ||
		 ( q[j] >= qisol[k] && qisol[k] > q[j1]) ){
		for(i=0;i<2;i++)
		  xp[i] = x[j][i] + (x[j1][i]-x[j][i])*
		    (qisol[k] - q[j])/(q[j1] - q[j]);
		fprintf(gnu," %7.5g %7.5g %6.4g \n",xp[0],xp[1],qisol[k]);
		iplot = 1;
	      }
	    }
	  }
	  if(iplot) fprintf(gnu,"\n\n");
	}
	ve=ve->Next;
      }
      

      /*  plottig of boundary */
      se=m->SurfaceElement; 
      while(se){ 
	double xp[2];
	if(se->surfnr != m->iper[2*lev_per] && 
	   se->surfnr != m->iper[2*lev_per+1] ){
	  xp[0] = se->P[0]->Point[0]+level*m->rper[2*lev_per];
	  xp[1] = se->P[0]->Point[1]+level*m->rper[2*lev_per+1];
	  fprintf(gnu," %7.5g %7.5g %4.2g \n",xp[0],xp[1],-1.);
	  
	  xp[0] = se->P[1]->Point[0]+level*m->rper[2*lev_per];
	  xp[1] = se->P[1]->Point[1]+level*m->rper[2*lev_per+1];
	  fprintf(gnu," %7.5g %7.5g %4.2g \n\n\n",xp[0],xp[1],-1.);
	}
	se=se->Next;
      }
    }
    /*  } */
}

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

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

    for(j=0;j<4;j++)
      for(k=0;k<s->DOF;k++){
	VI[k][j] = 0;   /* clearing of arrays */
	Dx[k][j] = 0;
	Dy[k][j] = 0;
      }

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

    }  /* end of noncurved elements */
    /* computing of the term Ds(w)\nabla w */
    for(j=0;j<s->DOF;j++){   /* j goes over integration nodes */
      if(ic >= 0){ /*  curved elements another values of wx,wy,nablarho*/
      for(l=0;l<4;l++)
       wx[l]= wy[l]=0;
      for(l=0;l<s->DOF;l++){
          for(k=0;k<4;k++){
            wx[k]+=s->curvedelement[ic].dertest[l][0][j]*
                   s->element[i].w[l][k];
            wy[k]+=s->curvedelement[ic].dertest[l][1][j]*
                   s->element[i].w[l][k];
          }
        }
        for(k=0;k<4;k++){  /* dertestfunction is in fact \partial\phi *area  */
          wx[k]*=2/s->curvedelement[ic].detJF[j];
          wy[k]*=2/s->curvedelement[ic].detJF[j];
        }
    /*  storing the value for edge integral terms */
	s->curvedelement[ic].nablarho[j][0] = wx[0]; 
	s->curvedelement[ic].nablarho[j][1] = wy[0]; 
      }
      nablarho[0] = wx[0];
      nablarho[1] = wy[0];
      /*      ViscousTerm(s->element[i].w[j],s->element[i].w[j],
                  nablarho,wx,wy,s->Re,s->kappaRePr, Dx[j],Dy[j]);
      */

      wx[1] /=s->element[i].w[j][0];
      wx[2] /=s->element[i].w[j][0];
      wy[1] /=s->element[i].w[j][0];
      wy[2] /=s->element[i].w[j][0];
      
      det = 2./3.*(wx[1] + wy[2] );

      Dx[j][1] = wx[1] + wx[1] - det;
      Dx[j][2] = wy[1] + wx[2];

      Dy[j][1] = wy[1] + wx[2];
      Dy[j][2] = wy[2] + wy[2] - det;

          if(i == itest ){
            for(l= 0;l<4;l++)
               printf("C: %d %d %d %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g %7.5g \n",i,j,l,
               s->element[i].w[j][0],s->element[i].w[j][1],
		      s->element[i].w[j][2],wx[1],wy[1],
               wx[2],wy[2],Dx[j][l],Dy[j][l]);
          } 
    }

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

}

