#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "Memory.h"
#include "FileIO.h"
#include "Matrix.h"
#include "Mesh.h"
#include "Progress.h"
#include "Solver.h"
#include "Triangle.h"
#include "Types.h"
#include "Implicit.h"
#include "Time.h"

int Solver_Implicit_Init(tSolver_State *s){
  tSolver_Matrix_Element *me;
  uint i,j,k,p,q,counter;
  int ic;
  const int itest=  -10;
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};
  double a1,a2;
  
  s->rows=s->nelem*s->DOF*4;  /* #rows of matrix = #unknowns */
  s->stencil=28;   /* =4*7,  stencil for vertex associated test functions */

  if(!Memory_GetM(s->matrix_row,tSolver_Matrix_Element*,s->nelem*s->DOF*4)||
     !Memory_GetM(s->precond_row,tSolver_Matrix_Element*,s->nelem*s->DOF*4)||
     !Memory_GetM(s->matrix,double,s->rows*s->stencil) ||
     !Memory_GetM(s->matrix_index,int,s->rows*s->stencil) ||
     !Memory_GetM(s->matrix_RHS,double,s->rows) )
    return 0;


  Clear_Matrix_Init(s);

  /* old version of matrix representation */
  /*
    for(i=0;i<s->nelem;i++){
    for(j=0;j<s->DOF;j++){
      for(k=0;k<4;k++){
	int i1;
	uint j1,k1,l,l1;
	index2element(p,i,j,k);
	for(l=0;l<s->DOF;l+=2){  // two neighboring triangles, first&last 
	  i1=s->element[i].neigh[perm[j][l]];
	  if(i1>=0){  // inner element 
	    uint j2,k2;
	    l1=s->element[i].neighind[perm[j][l]];
	    for(j1=0;j1<2;j1++)
	      for(k1=0;k1<4;k1++){
		k2=3-k1;
		index2element(q,i1,perm[l1][j1],k2); 
		Matrix_ElementAdd(s,0,p,q);  
	      } // end of k1,j1 
	      } // end of if(i1) 
	}

	for(j1=0;j1<s->DOF;j1++){
	  uint j2,k2;
	  j2=s->DOF-1-j1;
	  for(k1=0;k1<4;k1++){
	    k2=3-k1;
	    index2element(q,i,j2,k2);  // the same triangle 
	    Matrix_ElementAdd(s,0,p,q);  
	  }  // end of k1 
	}  // end of j1 

      }  // end of k 
    }  // end of j 
    if(i == itest) Show_matrix(s,i);
  }
  */

  /*  NEW ARRAYS !!!!  */
  for(i=0;i<s->nelem;i++){
    for(j=0;j<s->DOF;j++){
      for(k=0;k<4;k++){
	int i1;
	int j1,k1,l,l1;
	counter = 0;
	index2element(p,i,j,k);

	for(j1=0;j1<s->DOF;j1++){   /* the same triangle */
	  uint j2,k2;
	  //	  j2=perm[j][j1];
	  for(k1=0;k1<4;k1++){
	    index2element(q,i,j1,k1); 
	    s->matrix_index[p*s->stencil+counter] = q;
	    counter++;
	  }  /* end of k1 */
	}  /* end of j1 */


	for(l=2;l>=0;l-=2){      /* two neighboring triangles, first&last */
                                 /* have nonzero "trace" */
	  i1=s->element[i].neigh[perm[j][l]];
	  if(i1>=0){  /* inner element */
	    uint j2,k2;
	    l1=s->element[i].neighind[perm[j][l]];
	    for(j1=1;j1>=0;j1--) /* only two adjacent vertices are in stencil*/
	      for(k1=0;k1<4;k1++){
		index2element(q,i1,perm[l1][j1],k1); 
		s->matrix_index[p*s->stencil+counter] = q;
		counter++;
	      } /* end of k1,j1 */
	  }
	  else { /* else of if(i1) */
	    counter+=8;  /* i is neighbouring=> we skip the space */
		    
	  } /* end of if(i1) */
	}
	/*
	  printf("Line %d\n",p);
	  me=s->matrix_row[p];
	  while(me){
	  printf("%d, ",me->q);
	  me=me->Next;
	  }
	  printf("\n");
	  for(l=0;l<s->stencil;l++){
	  if(s->matrix_index[p*s->stencil+l] >=0)
	  printf("%d, ",s->matrix_index[p*s->stencil+l]);
	  }
	  printf("\n");
	*/

      }  /* end of k */
    }  /* end of j */
    if(i == itest) Show_matrix(s,i);
  }  /* end of i */
  Compute_A_Init(s->A1,s->A2,s->kappa);
  /*  Solver_Implicit_Output_Init(s,s->nelem*s->DOF*4,28); */

  if(s->implicit_method == 2){
    /* implict, second order in time */
    /* wold = w^{k}, K1 = w^{k-1} */
    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->element[i].K1[j][k] = s->element[i].w[j][k];
	}
      }
    }
  }

}
void Solver_Check_matrix(tSolver_State *s){
  tSolver_Matrix_Element *me;
  uint i,j,k,l,p,q,i1,j1,k1;
  static double w_ex[4]={1,1,0,4.5};
  double sum, tol = 1E-06;

  printf("\nChecking the matrix elements, Tol=%7.5g\n",tol);
  for(i=0;i<s->nelem;i++)
    for(j=0;j<s->DOF;j++)
      for(k=0;k<4;k++){
	sum = 0;
	index2element(p,i,j,k);
	for(q=0;q<s->stencil;q++){
	  l=p*s->stencil+q;
	  if(s->matrix_index[l] >=0){
	    element2index(s->matrix_index[l],i1,j1,k1);
	    sum+=s->matrix[l]*s->element[i1].w[j1][k1];
	    if(p == 1){
	      printf("## %d %d %d %d %10.7g %10.7g %10.7g %10.7g %10.7g \n",
		     p,q,l,s->matrix_index[l],s->matrix[l],
		     s->element[i1].w[j1][k1],
		     s->matrix[l]*s->element[i1].w[j1][k1],sum,
		     s->matrix_RHS[p]);
	    }
	    
	  }
	}
	if(fabs(sum-s->matrix_RHS[p]) > tol){
	    printf("Results: %d(%d,%d,%d): %7.5g =? %7.5g, reziduum =%7.5g \n",
	       p,i,j,k,sum,s->matrix_RHS[p],sum-s->matrix_RHS[p]);
	  /*	  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 \n",
		 s->element[i].ve->P[2]->Point[0],
		 s->element[i].ve->P[2]->Point[1]);
	*/
	}
      }
}

void Show_matrix(tSolver_State *s, uint i){
  tSolver_Matrix_Element *me;
  uint j,k,p,q,i1,j1,k1;
  uint count;

/* old version of matrix representation */
  /*
     
  for(j=0;j<s->DOF;j++)
    for(k=0;k<4;k++){
      index2element(p,i,j,k);
      
      me=s->matrix_row[p];
      printf("A[%d(%d,%d,%d)][(. . .)]=\n",p,i,j,k);
      count=0;
      while(me){
	element2index(me->q,i1,j1,k1);

	printf("%7.5g(%d:%d,%d,%d), ",me->a,me->q,i1,j1,k1);
	count++;
	if(count==4){
	  printf("\n");
	  count=0;
	}
	me=me->Next;
      }
      printf("\n");
    }
      printf("\n*********** O L D *************\n\n");
  */

  for(j=0;j<s->DOF;j++){
    for(k=0;k<4;k++){
      index2element(p,i,j,k);
      count=0;

      printf("A[%d(%d,%d,%d)][(. . .)]=\n",p,i,j,k);
      for(q=0;q<s->stencil;q++){
	matrix2index(q,i1,j1,k1,s->stencil);
	printf("%7.5g(%d:%d,%d,%d), ",
	       s->matrix[p*s->stencil+q],p*s->stencil+q,q,j1,k1);
	
	/* printf("%d(%d:%d,%d,%d), ",
	   s->matrix_index[p*s->stencil+q],p*s->stencil+q,q,j1,k1);*/
	count++;
	if(count==4){
	  printf("\n");
	  count=0;
	}
      }
      printf("\n");
    }
  }
    
  printf("b[(%d, . .)]=\n",i);
  for(j=0;j<s->DOF;j++){
    for(k=0;k<4;k++){
      index2element(p,i,j,k);
      printf("%7.5g(%d,%d), ",s->matrix_RHS[p],j,k);
    }
    printf("\n");
  }

}
void Show_matrix_precond(tSolver_State *s, uint i){
  tSolver_Matrix_Element *me;
  uint j,k,p,q,i1,j1,k1;
  uint count;

  for(j=0;j<s->DOF;j++)
    for(k=0;k<4;k++){
      index2element(p,i,j,k);
      
      me=s->precond_row[p];
      printf("A*[(%d,%d,%d)][(. . .)]=\n",i,j,k);
      count=0;
      while(me){
	element2index(me->q,i1,j1,k1);

	printf("%7.5g(%d,%d,%d), ",me->a,i1,j1,k1);
	count++;
	if(count==4){
	  printf("\n");
	  count=0;
	}
	me=me->Next;
      }
      printf("\n");
    }
}

int Matrix_ElementAdd(tSolver_State *s, double a, uint p, uint q){
tSolver_Matrix_Element *me, *me1;


/* matrix A  */
 if(!s || !Memory_GetM(me,tSolver_Matrix_Element,1) )
    return 0;
  me->a = a;
  me->p = p;
  me->q = q;

  me->Last_row=FALSE;
  me->Prev=0;
  me->Next=s->matrix_row[p]; 
  if(me->Next)
    me->Next->Prev=me;
  else
    me->Last_row=TRUE;
  s->matrix_row[p]=me;
  s->nummatrix++;

  return 1;

} 
void Solver_Implicit(tSolver_State *s){
  uint l,i,j,k;
  static uint impl_level_iter=1;

  //  if(s->limiting >0) Solver_Limiting_Vertex(s);
  if(s->limiting >0) Solver_Limiting_Vertex(s);
  
  for(l=0;l<impl_level_iter;l++){

    if(s->implicit_method == 2){
      /* implict, second order in time */
      
    }
    Clear_Matrix(s);
    /*    printf("$$$ %7.5g %7.5g %7.5g %7.5g \n",
	   s->element[0].wold[0][0],s->element[0].wold[0][1],
	   s->element[0].wold[0][2],s->element[0].wold[0][3]);
    printf("... %7.5g %7.5g %7.5g %7.5g \n-----------------\n",
	   s->element[0].w[0][0],s->element[0].w[0][1],
	   s->element[0].w[0][2],s->element[0].w[0][3]);
    */
    /*for(i=0;i<100;i++)
      printf("%7.5g(%d,%d)",s->matrix[i],i,s->matrix_index[i]);*/
    

  if(s->implicit_method == 2){
    /* implict, second order in time */
    /* wold = w^{k}, K1 = w^{k-1} */
    for(i=0; i<s->nelem; i++){
      for(j=0; j<s->DOF; j++){
	for(k=0; k<4; k++){
	  s->element[i].w[j][k] = 4*s->element[i].wold[j][k] 
	    -s->element[i].K1[j][k];
	}
	/*	printf("(%d %d): %7.5g %7.5g %7.5g %7.5g\n",i,j,
	       s->element[i].w[j][0],s->element[i].wold[j][1],
	       s->element[i].K1[j][0],s->element[i].w[j][3]);
	*/
      }
    }
  }


    Solver_Implicit_Time(s);
    /*     printf("\n\nSolver_Implicit_Time - done\n");   */

    /*Show_matrix(s,0); */

  if(s->implicit_method == 2){
    /* implict, second order in time */
    /* wold = w^{k}, K1 = w^{k-1} */
    for(i=0; i<s->nelem; i++){
      for(j=0; j<s->DOF; j++){
	for(k=0; k<4; k++){
	  s->element[i].w[j][k] = 2*s->element[i].wold[j][k] 
	    -s->element[i].K1[j][k];
	}
	/*
	printf("(%d %d): %7.5g %7.5g %7.5g %7.5g\n",i,j,
	       s->element[i].w[j][0],s->element[i].wold[j][1],
	       s->element[i].K1[j][0],s->element[i].w[j][3]);
	*/
      }
    }
  }

    Solver_Implicit_Volume(s);
    /*    printf("\n\nSolver_Implicit_Volume - done\n");  */

    Solver_Implicit_Edge(s);
    /*     printf("\n\nSolver_Implicit_Edge - done\n");  */

    /*Solver_Check_matrix(s); */

    /*  Show_matrix(s,0); */

  /* Solver_Implicit_Preconditioning(s,s->nelem*s->DOF*4,28); */
    Solver_Implicit_Output(s,s->rows,s->stencil);
  /*       printf("Solver_Implicit_Output - done\n");  */

  if(s->implicit_method == 2){
    /* implict, second order in time */
    /* wold = w^{k}, K1 = w^{k-1} */
    for(i=0; i<s->nelem; i++){
      for(j=0; j<s->DOF; j++){
	for(k=0; k<4; k++){
	  s->element[i].K1[j][k] = s->element[i].wold[j][k];
	}

      }
    }
  }

  }
}

void Solver_Implicit_Time(tSolver_State *s){
  uint i,j,k,row,p,q,q1,q2;
  int ic;
  const int itest=  -10;
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};
  double a1,a2,aa1,aa2;
  tSolver_Matrix_Element *me;


  /* better for linear system */
  /*  for(i=0;i<s->nelem;i++) s->element[i].tauloc=1; */
  
  for(i=0;i<s->nelem;i++){ /*time derivative terms */

    ic = s->element[i].icve;
    for(j=0;j<s->DOF;j++){
      if( ic>=0){
	a1=s->curvedelement[ic].detJF[perm[j][2]]/12;
	a2=s->curvedelement[ic].detJF[perm[j][1]]/12;
      }
      else{
	a1 = 1/s->element[i].oneovervolume/12; /* 12= 3  *  2  *  2    */
	a2 = 1/s->element[i].oneovervolume/12; /* weight,aver(i,j), phi(x_l) */
      }
      a1/=s->element[i].tauloc;
      a2/=s->element[i].tauloc;

      aa1 = a1;
      aa2 = a2;

      if(s->implicit_method == 2){
	/* implict, second order in time */
	a1=a1*3/2;   /* ?????? a1*=3/2 does not work ??????????  */
	a2=a2*3/2;
	aa1=aa1/2;
	aa2=aa2/2;
      }

      for(k=0;k<4;k++){
	double rhs;

	/* old version of matrix representation */
	/*
     	index2element(p,i,perm[j][0],k);
	index2element(q1,i,perm[j][1],k);
	index2element(q2,i,perm[j][2],k);
	Update_Matrix_Element(s,p,p,a1+a2);
	Update_Matrix_Element(s,p,q1,a1);
	Update_Matrix_Element(s,p,q2,a2);
	*/

	index2element(row,i,perm[j][0],k); 
	index2matrix(p,row,perm[j][0],k,s->stencil);
	s->matrix[p]+=a1+a2;
	index2matrix(p,row,perm[j][1],k,s->stencil);
	s->matrix[p]+=a1;
	index2matrix(p,row,perm[j][2],k,s->stencil);
	s->matrix[p]+=a2;

	rhs = (aa1+aa2)*s->element[i].w[perm[j][0]][k] 
	  + aa1*s->element[i].w[perm[j][1]][k]
	  + aa2*s->element[i].w[perm[j][2]][k];
	//	Update_Matrix_RHS(s,p,rhs);
	s->matrix_RHS[row]+=rhs;
	/*	printf("!!! (%d,%d) %d %d %d %7.5g %7.5g \n",
		       i,j,p,q1,q2,a1,a2); */
      }
     
    } 
    if(i == itest){
      printf("Time step %7.5g\n",s->element[i].tauloc);
      Show_matrix(s,i); 
            printf("----------- END OF TIME DERIVATIVES ---------\n");

    }
  }
  /*  i=0;
  me=s->matrix_row[i];
  printf("Element %d %d %d %7.5g\n",i,me->p,me->q,me->a);
  me=me->Next;
  printf("Element %d %d %d %7.5g\n",i,me->p,me->q,me->a);
  */
}
void Solver_Implicit_Volume(tSolver_State *s){
  uint i,j,k,l,row,p,q1,q2;
  double a;
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};
  matrix4 A1,A2; 
  const int itest= -99;
  
  /*  Compute_A_Init(A1,A2,s->kappa); */


  for(i=0;i<s->nelem;i++){ /* volume Integral terms */
    int ic;
    double4 wi,wqi,wq1,wq2;
    uint m,n;
    double nablaphi[2];
    ic=s->element[i].icve;  

    for(l=0;l<3;l++){  /* l goes over integration nodes */
      /* conservative variables */
      /*      
       double4_aver(wi,s->element[i].w[perm[l][1]],s->element[i].w[perm[l][2]]);
      */

      /* physical variables */
      /*      printf("Conser: %d %d %7.5g %7.5g %7.5g %7.5g \n",i,l,
	     s->element[i].w[perm[l][1]][0],s->element[i].w[perm[l][1]][1],
	     s->element[i].w[perm[l][1]][2],s->element[i].w[perm[l][1]][3]); 
      */     
     
      /*      printf("Physic: %d %d %7.5g %7.5g %7.5g %7.5g \n",i,l,
	     wq1[0],wq1[1],wq1[2],wq1[3]); 
      vector_q2w(wq1,wi,s->kappa);
      printf("Bacwar: %d %d %7.5g %7.5g %7.5g %7.5g \n--------------\n",i,l,
	     wi[0],wi[1],wi[2],wi[3]); 
      */
      vector_w2q(s->element[i].w[perm[l][1]],wq1,s->kappa);
      vector_w2q(s->element[i].w[perm[l][2]],wq2,s->kappa);
      double4_aver(wqi,wq1,wq2);
      vector_q2w(wqi,wi,s->kappa); 

      if(Compute_A(wi,s->A1,s->A2,s->kappa) )  printf("Problem!!\n"); 
      /*if(Compute_A_Linear(wi,s->A1,s->A2,s->kappa) )printf("Problem!!\n"); */

      for(j=0;j<s->DOF;j++){ /* j goes over test functions */
	if(ic>=0){  /* curved elements */
	  for(m=0;m<2;m++)
	    nablaphi[m]=s->curvedelement[ic].dertest[j][m][l];
	  /*	  *2/  s->curvedelement[ic].detJF[l];  */
	}
	else{  /* noncurved elements */
	  for(m=0;m<2;m++)
	    nablaphi[m]=s->element[i].dertest[j][m];
	  /*  /s->element[i].oneovervolume; */
	}

	for(m=0;m<4;m++)
	  for(n=0;n<4;n++){
	    /* old version of matrix representation */
	    /*
     	    index2element(p,i,j,m);
	    index2element(q1,i,perm[l][1],n);
	    index2element(q2,i,perm[l][2],n);
	    a=-(s->A1[m][n]*nablaphi[0] + s->A2[m][n]*nablaphi[1])/6;// 2*3 
	    Update_Matrix_Element(s,p,q1,a);
	    Update_Matrix_Element(s,p,q2,a);
	    */

	    index2element(row,i,j,m);
	    index2matrix(q1,row,perm[l][1],n,s->stencil);
	    index2matrix(q2,row,perm[l][2],n,s->stencil);
	    a=-(s->A1[m][n]*nablaphi[0] + s->A2[m][n]*nablaphi[1])/6;/* 2*3 */ 
	    s->matrix[q1]+=a;
	    s->matrix[q2]+=a;
	  }

      }  /* end of j */
    } /* end of l */
    /*if(i == itest){
      Show_matrix(s,i); 
      printf("----------- END OF VOLUME -------------\n");
      }
    */
  } /* end of i */
}
void Solver_Implicit_Edge(tSolver_State *s){
  uint i,j,k,l,p,row,q1,q2,q3,q4,qq1,qq2,qq3,qq4;
  double a1,a2,a3,a4,slip_pressure;
  static uint perm[3][3]={{0,1,2},{1,2,0},{2,0,1}};
  static uint perm2[4][2]={{0,1},{1,0},{0,1},{1,0}};
  matrix4 Pp,Pm;
  double length,phi;
  double2 n_unit;
  const double t[2]={0.7886751345948129,0.2113248654051871};
  const int itest=  -10;

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

  for(k=0;k<s->nface;k++){ /* edge integral terms */
    uint ii,jj,z;
    double4 w,wi,wj,wq1,wq2,wqi,wqj,wq;
    double area;
    int ic;
    uint m,n;
    BOOL slip;

    i=s->face[k].i;
    ii=s->face[k].ii;
    area = 1/s->element[i].oneovervolume;

    for(l=0;l<2;l++) n_unit[l]=s->face[k].n[l];
    length=s->face[k].volume;

    for(z=0;z<2;z++){  /* integration nodes on one edge */

      /* conservative variables */
      /*      double4_lin_comb(wi,s->element[i].w[perm[ii][0]],
		       s->element[i].w[perm[ii][1]],
		       t[perm2[z][0]],t[perm2[z][1]]);
      */

      /* physical variables */
      vector_w2q(s->element[i].w[perm[ii][0]],wq1,s->kappa);
      vector_w2q(s->element[i].w[perm[ii][1]],wq2,s->kappa);
      double4_lin_comb(wqi,wq1,wq2, t[perm2[z][0]],t[perm2[z][1]]);

      slip=FALSE;
      if(s->face[k].boundary<=0){/* inner or periodic edge */
	j=s->face[k].j;
	jj=s->face[k].jj;

	/* conservative variables */
	/*        double4_lin_comb(wj,s->element[j].w[perm[jj][0]],
			 s->element[j].w[perm[jj][1]],
			 t[perm2[z][1]],t[perm2[z][0]]);
	*/
	/* physical variables */
	vector_w2q(s->element[j].w[perm[jj][0]],wq1,s->kappa);
	vector_w2q(s->element[j].w[perm[jj][1]],wq2,s->kappa);
	double4_lin_comb(wqj,wq1,wq2, t[perm2[z][1]],t[perm2[z][0]]);
	
      } /* end of inner edge */
      else if(s->boundary[s->face[k].boundary-1].type==
	      cSolver_BoundaryType_Slip){

	/* conservative variables */
	/*	double4_lin_comb(wi,s->element[i].w[perm[ii][0]],
			 s->element[i].w[perm[ii][1]],
			 t[perm2[z][0]],t[perm2[z][1]]); */

	/*	printf("Conser: %d %d %7.5g %7.5g %7.5g %7.5g \n",i,z,
		wi[0],wi[1],wi[2],wi[3]); */

	/* physical variables rho, v1, v2, p */
	vector_w2q(s->element[i].w[perm[ii][0]],wq1,s->kappa);
	vector_w2q(s->element[i].w[perm[ii][1]],wq2,s->kappa);
	double4_lin_comb(wqi,wq1,wq2, t[perm2[z][0]],t[perm2[z][1]]);
	vector_q2w(wqi,wi,s->kappa);

	/**  !!!!  CHANGE WHEN vector_w2q and vector_q2w are CHANGED */
	slip_pressure=wqi[3]; 
	/* slip_pressure=wqi[3]*wqi[3]*wqi[0]/s->kappa; */

	/*	printf("Physic: %d %d %7.5g %7.5g %7.5g %7.5g \n------\n",i,z,
		wi[0],wi[1],wi[2],wi[3]); */

	slip=TRUE;

      }
      else{ /* boundary edge */  /* only Inlet & Outlet */
	  tMesh_VolumeElement *ve; 
	  tMesh_Point *p1,*p2;
	  double2 xp;
	  double4 qj1;

	  /*  for evaluation of Pp and Pm  */
	  vector_w2q(s->element[i].w[perm[ii][0]],wq1,s->kappa);
	  vector_w2q(s->element[i].w[perm[ii][1]],wq2,s->kappa);
	  /*	printf("Physic: %d %d %7.5g %7.5g %7.5g %7.5g \n------\n",i,z,
		wi[0],wi[1],wi[2],wi[3]); */
	  double4_lin_comb(wqi,wq1,wq2,t[perm2[z][0]],t[perm2[z][1]]);
	  vector_q2w(wqi,wi,s->kappa);

	  if(s->init !=4 && s->init !=5) { /* Original  BC */
	    Set_Boundary_State(s,k,i,wi,wj); /* only Inlet & Outlet */
	    vector_w2q(wj,wqj,s->kappa);
	  }
	  else{   /* 	  RINGLEB FLOW */
	    /* points forming the boundary edge */
	    ve=s->element[i].ve;
	    p1=ve->P[perm[ii][0]];
	    p2=ve->P[perm[ii][1]];
	    /*
	      printf("%7.5g %7.5g \n",p1->Point[0],p1->Point[1]);
	      printf("%7.5g %7.5g \n",p2->Point[0],p2->Point[1]);
	    */
	    
	    /* computation of integration nodes */
	    double2_weight_awer(xp,p1->Point,p2->Point,
				t[perm2[z][0]],t[perm2[z][1]]);
	    /*
	      printf("%7.5g %7.5g \n",xp[0][0],xp[0][1]);
	      printf("%7.5g %7.5g \n\n",xp[1][0],xp[1][1]);
	    */
	    
	    /* computation of the solution in integration nodes */
	    calc_ringleb(s->kappa,xp,qj1);
	    /* transform form physical quantities to the state vector */
	    vector_or_q2w(qj1,wj,s->kappa);
	    vector_w2q(wj,wqj,s->kappa);
	  } /* end of Ringleb */

      }   /* end of Inlet/outlet */
   
      /* LINEAR SYSTEM !!!!! */
      /*     if(slip) 
	     for(m=0;m<4;m++) wj[m]=wi[m];
	     slip = FALSE;
	     } */
      if(!slip){

	double4_aver(wq,wqi,wqj);

	/**  !!!!  CHANGE inside Matrixpm WHEN vector_w2q(q2w) are CHANGED */

	Matrixpm(wq,area,n_unit,length,Pp,Pm,s->element[i].tauloc,
		s->kappa,s->CFL);  
	/* Matrixpm_Linear(w,area,n_unit,length,Pp,Pm,s->element[i].tauloc,
	   s->kappa,s->CFL); */
      }
      /*
      else{
	int ic;
	uint k1,l1,m1;

	ic=s->element[i].icve;	
	if(ic >= 0){  // curved elements  
	  Compute_A_slip(wi, s->curvedelement[ic].normal[z],
			 s->curvedelement[ic].length[z],Pp, s->kappa);
	}
	else{ //*non curved elements 
	  Compute_A_slip(wi, n_unit, length, Pp, s->kappa);
	}
	slip=0;  
      }
      */
	/*	
	  printf("As(w=(%7.5g,%7.5g,%7.5g,%7.5g))        n=(%7.5g,%7.5g)),:%d\n",
	  wi[0],wi[1],wi[2],wi[3],
	  s->curvedelement[ic].normal[z][0]*s->curvedelement[ic].length[z],
	  s->curvedelement[ic].normal[z][1]*s->curvedelement[ic].length[z],
	  ic);
	  
	  for(k1=0;k1<4;k1++){
	  for(l1=0;l1<4;l1++)
	  printf("%7.5g,",Pp[k1][l1]);
	  printf("\n");
	  }
	  printf("\n\n");
	*/
      
      
      /*      if((i== itest || j==itest) && !slip){ 
	uint k1,l1,m1;
	printf("Pp(w=(%7.5g,%7.5g,%7.5g,%7.5g))        Pm(n=(%7.5g,%7.5g))\n",
	       w[0],w[1],w[2],w[3],n_unit[0]*length,n_unit[1]*length);
	for(k1=0;k1<4;k1++){
	  for(l1=0;l1<4;l1++)
	    printf("%7.5g,",Pp[k1][l1]);
	  printf("       ");
	  for(l1=0;l1<4;l1++)
	    printf("%7.5g,",Pm[k1][l1]);
	  printf("\n");
	}
	printf("\n");
      } 
      */  /* end of itest */

      for(l=0;l<2;l++){ /* over test functions in the stencil, only 2!!!! */
	if(!slip){
	  phi=t[perm2[z+l][0]]/2;  /*  0.5 is a weight */
	  for(m=0;m<4;m++){
	    for(n=0;n<4;n++){
	      /* test functions on K_i */
	      /* old version of matrix representation */
	      /*
     	      index2element(p,i,perm[ii][l],m);
	      index2element(q1,i,perm[ii][0],n);
	      index2element(q2,i,perm[ii][1],n);
	      a1=Pp[m][n]*phi*t[perm2[z][0]];
	      a2=Pp[m][n]*phi*t[perm2[z][1]];
	      Update_Matrix_Element(s,p,q1,a1);
	      Update_Matrix_Element(s,p,q2,a2);
	      */

	      index2element(row,i,perm[ii][l],m);
	      index2matrix(qq1,row,perm[ii][0],n,s->stencil);
	      index2matrix(qq2,row,perm[ii][1],n,s->stencil);
	      a1=Pp[m][n]*phi*t[perm2[z][0]];
	      a2=Pp[m][n]*phi*t[perm2[z][1]];
	      s->matrix[qq1]+=a1;
	      s->matrix[qq2]+=a2;

	      /*
		printf("Test A: %d %d %d %d %d %d\n",
		p,q1,q2,row,s->matrix_index[qq1],s->matrix_index[qq2]);
	      */

	      if(s->face[k].boundary<=0){/* inner edge on K_i  */
		/* old version of matrix representation */
		/*
     		index2element(q3,j,perm[jj][1],n);
		index2element(q4,j,perm[jj][0],n);
		a3=Pm[m][n]*phi*t[perm2[z][0]];
		a4=Pm[m][n]*phi*t[perm2[z][1]];
		Update_Matrix_Element(s,p,q3,a3);
		Update_Matrix_Element(s,p,q4,a4);
		*/
		index2matrix(qq3,row,5-2*l,n,s->stencil);
		index2matrix(qq4,row,5-2*l+1,n,s->stencil);
		a3=Pm[m][n]*phi*t[perm2[z][0]];
		a4=Pm[m][n]*phi*t[perm2[z][1]];
		s->matrix[qq3]+=a3;
		s->matrix[qq4]+=a4;

		/*
		  printf("Test B: %d %d %d %d %d %d\n",
		  p,q3,q4,row,s->matrix_index[qq3],s->matrix_index[qq4]);
		*/
	      } 
	    } /* end of n */
	  } /* end of m */

	  if(s->face[k].boundary<=0){ /* test functions on K_j */
	    phi=t[perm2[z+l][1]]/2;  /*  0.5 is a weight */
	    for(m=0;m<4;m++){
	      for(n=0;n<4;n++){
		/* old version of matrix representation */
		/*
     
		index2element(p,j,perm[jj][l],m);
		index2element(q1,i,perm[ii][0],n);
		index2element(q2,i,perm[ii][1],n);
		index2element(q3,j,perm[jj][1],n);
		index2element(q4,j,perm[jj][0],n);
		
		a1=-Pp[m][n]*phi*t[perm2[z][0]];
		a2=-Pp[m][n]*phi*t[perm2[z][1]];
		a3=-Pm[m][n]*phi*t[perm2[z][0]];
		a4=-Pm[m][n]*phi*t[perm2[z][1]];
		
		Update_Matrix_Element(s,p,q1,a1);
		Update_Matrix_Element(s,p,q2,a2);
		Update_Matrix_Element(s,p,q3,a3);
		Update_Matrix_Element(s,p,q4,a4);
		*/

		index2element(row,j,perm[jj][l],m);
		index2matrix(qq4,row,perm[jj][0],n,s->stencil);
		index2matrix(qq3,row,perm[jj][1],n,s->stencil);
		index2matrix(qq2,row,5-2*l,n,s->stencil);
		index2matrix(qq1,row,5-2*l+1,n,s->stencil);

		a1=-Pp[m][n]*phi*t[perm2[z][0]];
		a2=-Pp[m][n]*phi*t[perm2[z][1]];
		a3=-Pm[m][n]*phi*t[perm2[z][0]];
		a4=-Pm[m][n]*phi*t[perm2[z][1]];

		s->matrix[qq1]+=a1;
		s->matrix[qq2]+=a2;
		s->matrix[qq3]+=a3;
		s->matrix[qq4]+=a4;

		/*
		  printf("Test C: %d %d %d %d %d\n",
			 p,q1,q2,q3,q4);
		  printf("Test  : %d %d %d %d %d\n\n",
			 row,s->matrix_index[qq1],s->matrix_index[qq2],
			 s->matrix_index[qq3],s->matrix_index[qq4]);
		*/
		
	      } /* end of n */
	    } /* end of m */
	  }  /* end of test functions on K_j */
	}  /* end of !slip  */

	
	if(s->face[k].boundary>0){ /* boundary edge ==> right hand side */
	  phi=t[perm2[z+l][0]]/2;  /*  0.5 is a weight */

	  if(slip){  /* impermeable walls  */
	    double24 Fb;
	    int m1;
	    double ah;

	    Set_Slip(s,k,i,z,slip_pressure ,wq, w, Fb);
	    
	    for(m=1;m<3;m++){  /* linear part of linearization */
	      index2element(row,i,perm[ii][l],m);
	      //printf("L.H.S. [%d %d %d %d %d]\n",i,ii,l,m,row);
	      for(n=0;n<4;n++){ 
		index2matrix(qq1,row,perm[ii][0],n,s->stencil);
		index2matrix(qq2,row,perm[ii][1],n,s->stencil);
		a1=Fb[m-1][n]*phi*t[perm2[z][0]];
		a2=Fb[m-1][n]*phi*t[perm2[z][1]];
		s->matrix[qq1]+=a1;
		s->matrix[qq2]+=a2; 
		//if(row <2)
		//  printf("%d (%d) %7.5g %d (%d) %7.5g [%7.5g]\n",
		//	 qq1,perm[ii][0],a1,qq2,perm[ii][1],a2,t[perm2[z][0]]);
	      }
	      //	      printf("\n");
	    }
	    /* EW */
	    for(m=1;m<3;m++){ /* RIGHT_HAND_SIDE */
	      //printf("R.H.S. [%d %d]",ii,l);
	      index2element(row,i,perm[ii][l],m);
	      a3=-w[m]*phi;       /* absolute part */
	      for(m1=0;m1<4;m1++){  /* linear part */
		a3+=Fb[m-1][m1]*wi[m1]*phi; 
		/*
		  (s->element[i].w[perm[ii][0]][m1]*t[perm2[z][0]]+
		  s->element[i].w[perm[ii][1]][m1]*t[perm2[z][1]]); */
		//if(row <2)
		//printf("%7.5g %7.5g (%d,%d)\n",-Fb[m-1][m1]*phi,wi[m1],m1,row);
	      }
	      //printf("\n");
	      s->matrix_RHS[row]+=a3;
	    }	  
	    //printf("\n---------------------\n");


	    /* OLD  */
	    /*
	      for(m=1;m<3;m++){ // RIGHT_HAND_SIDE 
	      index2element(row,i,perm[ii][l],m);
	      a3=-w[m]*phi;
	      s->matrix_RHS[row]+=a3;
	    }	  
	    */
	  }
	  else{ /* inlet outlet */
	    Matrix_Times_Vector4(Pm,wj,w);
	
	    for(m=0;m<4;m++){ /* RIGHT_HAND_SIDE */
	      /*
		index2element(p,i,perm[ii][l],m);
		a3=-w[m]*phi;
		Update_Matrix_RHS(s,p,a3);
	      */
	      
	      index2element(row,i,perm[ii][l],m);
	      a3=-w[m]*phi;
	      s->matrix_RHS[row]+=a3;
	    }
	  }
	  /*	  if((i== itest || j==itest) && slip){
	    printf("bound state %7.5g %7.5g %7.5g %7.5g\n",
		   w[0],w[1],w[2],w[3]);
	  } */
	}
      } /* end of l */
    } /*end of z */
  } /* end of k */
  if(itest>=0) Show_matrix(s,itest); 
}

int Clear_Matrix_Init(tSolver_State *s){
  uint i;
  for(i=0;i<s->nelem*s->DOF*4;i++){
    s->matrix_row[i]=0;
    s->matrix_RHS[i]=0;
  }
  for(i=0;i<s->rows;i++)
    s->matrix_RHS[i]=0;
  for(i=0;i<s->rows*s->stencil;i++)
    s->matrix_index[i]=-1;
  
  return 0;
}

int Clear_Matrix(tSolver_State *s){
  tSolver_Matrix_Element *me;
  uint i,j,p;

  for(i=0;i<s->nelem*s->DOF*4;i++){
    s->matrix_RHS[i]=0;
    me=s->matrix_row[i];
    while(me){
      me->a=0;
      me=me->Next;
    }
  }

  for(i=0;i<s->rows;i++)
    s->matrix_RHS[i] = 0;
  for(i=0;i<s->rows*s->stencil;i++)
    s->matrix[i]=0;
  
  return 0;
}

int Update_Matrix_RHS(tSolver_State *s, uint p,double a){
  s->matrix_RHS[p]+=a;
  return 0;
}

int Update_Matrix_Element(tSolver_State *s, uint p,uint q, double a){
  tSolver_Matrix_Element *me;
  me=s->matrix_row[p];
  while(me){
    if(me->q == q){
      me->a+=a;
      return 1;
    }
    me=me->Next;
  }
  printf("*******  E L E M E N T   D O E S   N O T   F O U N D *****\n");
  return 0;
}
tSolver_Matrix_Element *Find_Precond_Element(tSolver_State *s, uint p, uint q){
  tSolver_Matrix_Element *me;
  me = s->precond_row[p];
  while(me){
    if(me->q == q)
      return me;
    me=me->Next;
  }
  return 0;
}

int Solver_Implicit_Preconditioning(tSolver_State *s, uint nmax, uint nzeromax){
  tSolver_Matrix_Element *me, *me1, *me_kk, *me_ij, *me_ik, *me_kj;
  double spau[nmax*nzeromax],spal[nmax*nzeromax];
  uint i,j,k,m,p,q,count;
  int itest = -10;

  if( itest >=0) Show_matrix(s,itest);

  for(i=0;i<nmax;i++){
    me = s->matrix_row[i];
    me1 = s->precond_row[i];
    while(me){
      me1->a = me->a;
      me=me->Next;
      me1=me1->Next;
    }
  }
    
  if( itest >=0) Show_matrix_precond(s,itest);

  for(k=0;k<nmax;k++){
    printf("prec iter %d\n",k);

    me_kk = Find_Precond_Element(s,k,k);
    for(i=k+1;i<nmax;i++){
      if( (me_ik = Find_Precond_Element(s,i,k)) ){
	me_ik->a/=me_kk->a;
	for(j=k+1;j<nmax;j++){
	  if( (me_ij = Find_Precond_Element(s,i,j)) &&
	      (me_kj = Find_Precond_Element(s,k,j)) )
	    me_ij->a -= me_ik->a*me_kj->a;
	}  /* end of j */
      }  /*end of if me_ik */      
    }  /* end of i  */
  }  /* end of k  */

  if( itest >=0) Show_matrix_precond(s,itest);

  for(i=0;i<nmax*nzeromax;i++){
    spau[i]=0;
    spal[i]=0;
  }
  count=0;
  for(p=0;p<nmax;p++){
    me=s->precond_row[p];
    while(me){
      q=me->q;
      if(p <=q){
	spau[count]=me->a;
	me1=s->matrix_row[q];
	while(me1){
	  if(me1->q == p){
	    spal[count]=me1->a;
	    break;
	  }
	  me1=me1->Next;
	}
	count++;
      }
      me=me->Next;
    }
  } 

  /* calling f77 subroutine */
  precondmatrix_(&nmax, &nzeromax, spau, spal);

}


int Solver_Implicit_Output_Init(tSolver_State *s, uint nmax, uint nzeromax){
  tSolver_Matrix_Element *me, *me1;
  uint i,j,k,m,p,q,count;
  int *index, *irwst;
  
  index=calloc(nmax*nzeromax,sizeof(int));
  irwst=calloc(nmax+1,sizeof(int));

  count=0;
  irwst[0]=count;
  for(p=0;p<nmax;p++){
    me=s->matrix_row[p];
    while(me){
      q=me->q;
      index[count]=q;
      count++;
      if(count > nmax*nzeromax) {
	printf("count > nmax*nzeromax in Solver_Implicit_Output_Init\n");
	return -1;
      }

      /*      if(p <=q){
	index[count]=q;
	count++;
	}*/
      me=me->Next;
    }
    irwst[p+1]=count;
  } 

  /* calling f77 subroutine */
  initmatrix_(&nmax, &nzeromax, irwst, index);

  free(index);
  free(irwst);
}

int Solver_Implicit_Output(tSolver_State *s, uint nmax, uint nzeromax){
  FILE *f1, *f2;
  tSolver_Matrix_Element *me, *me1;
  uint count,iterations,count1;
  uint i,j,k,m,p,q;
  double reziduum;
  /*  double *sparsu,*sparsl,*b,*x_init;*/
  double *spars,*b,*x_init;
  int *index, *irwst;
  double *spars1;
  int *index1, *irwst1;
  static double time0=0;
  double  time1=0, time2=0,time3=0,time4=0, timea=0, timeb=0;
  int *bcnods; 
  const int itest = -16;

  time1 = Time_CPUGet();

  /*  printf("F1A\n"); */

	       
  spars=calloc(nmax*nzeromax,sizeof(double));
  index=calloc(nmax*nzeromax,sizeof(int));
  irwst=calloc(nmax+1,sizeof(int));

  spars1=calloc(nmax*nzeromax,sizeof(double));
  index1=calloc(nmax*nzeromax,sizeof(int));
  irwst1=calloc(nmax+1,sizeof(int));

  b=calloc(nmax,sizeof(double));
  x_init=calloc(nmax,sizeof(double));
  bcnods = calloc(nmax,sizeof(int));

  for(i=0;i<nmax*nzeromax;i++){
    spars[i]=0;
    /*    sparsu[i]=0;
	  sparsl[i]=0; */
  }
  count=0;
  count1=0;
  irwst[0]=count+1; 
  irwst1[0]=count1+1; 
  timea = Time_CPUGet();
  for(p=0;p<nmax;p++){
    uint i,j,k;
    uint iyes = 0;
    element2index(p, i, j, k);

    me=s->matrix_row[p];

    if( s->element[i].limited == 1){/*on i-th element limitation was done */
      //if( s->element[i].limited == 1 && i > 50 && i<10){  /* any limitation*/
      uint l,l1;
      uint index_local[36];
      double spars_local[36];
      uint count_local=0;

      /*
      for(l=0;l<s->stencil;l++)
	printf("%d ",l);
      printf("\n");
      for(l=0;l<s->stencil;l++)
	printf("%d ",s->matrix_index[p*s->stencil+l]);
      printf("\n");
      for(l=0;l<s->stencil;l++)
	printf("%d ",s->matrix_index[(p+1)*s->stencil+l]);
      printf("\n");
      for(l=0;l<s->stencil;l++)
	printf("%d ",s->matrix_index[(p+2)*s->stencil+l]);
      printf("\n");
      */

      for(l=0;l<12;l++){ /* i-th element, first part of 3 rows */
	index_local[count_local]=s->matrix_index[p*s->stencil+l];
	spars_local[count_local]=s->matrix[p*s->stencil+l];
	for(l1=1;l1<3;l1++)
	  spars_local[count_local]+=s->matrix[(p+l1)*s->stencil+l];
	/*
	  printf("*** %d %d %d %7.5g \n",l,count_local,index_local[count_local],
	       spars_local[count_local]);
	*/
	count_local++;
      }

      for(l=0;l<8;l++){ /* first neighbour */
	index_local[count_local]=s->matrix_index[p*s->stencil+12+l];
	spars_local[count_local]=s->matrix[p*s->stencil+12+l];
	spars_local[count_local]+=s->matrix[(p+2)*s->stencil+20+l];
	/*
	printf("* * %d %d %d %d %7.5g \n",
	       l,count_local,index_local[count_local],
	       s->matrix_index[(p+2)*s->stencil+20+l],
	       spars_local[count_local]);
	*/
	count_local++;
      }

      for(l=0;l<8;l++){ /* second neighbour */
	index_local[count_local]=s->matrix_index[p*s->stencil+20+l];
	spars_local[count_local]=s->matrix[p*s->stencil+20+l];
	spars_local[count_local]+=s->matrix[(p+1)*s->stencil+12+l];
	/*
	printf("*.* %d %d %d %d %7.5g \n",
	       l,count_local,index_local[count_local],
	       s->matrix_index[(p+1)*s->stencil+12+l],
	       spars_local[count_local]);
	*/
	count_local++;
      }

      for(l=0;l<8;l++){ /* third neighbour */
	index_local[count_local]=s->matrix_index[(p+1)*s->stencil+20+l];
	spars_local[count_local]=s->matrix[(p+1)*s->stencil+20+l];
	spars_local[count_local]+=s->matrix[(p+2)*s->stencil+12+l];
	/*	printf("*!* %d %d %d %d %7.5g %7.5g %7.5g \n",
	       l,count_local,index_local[count_local],
	       s->matrix_index[(p+2)*s->stencil+12+l],
	       spars_local[count_local],
	       s->matrix[(p+1)*s->stencil+20+l],
	       s->matrix[(p+2)*s->stencil+12+l]);
	*/
	count_local++;
      }

      //      while(me){  /* first row of i-th element (for j=0)  */
      //	index_local[count_local]=me->q+1;
      //	spars_local[count_local]=me->a;
	/*printf("LOCAL %d %d %7.5g\n",
		count_local,index_local[count_local],
		spars_local[count_local]); */
      //	count_local++;
      //	me=me->Next;
      //}
      //for(l=1;l<=2;l++){ /* we add the second and third rows */
                         /* of i-th element (for j=1,2) to the first one */
      //	me=s->matrix_row[p+l];
      //	while(me){  
      //	  uint elem_find=0;
      //	  for(l1=0;l1<count_local;l1++){
      //    if(index_local[l1] == me->q+1){
      //      spars_local[l1]+=me->a;
	      /*	      printf("L_OLD %d %d %7.5g %7.5g\n",
			      l1,index_local[l1],spars_local[l1],me->a); */
      //      elem_find=1;
      //      break;
      //    }
      //	  }
      //  if(!elem_find){  /* the column q doen't find, we create a new */
      //    index_local[count_local]=me->q+1;
      //    spars_local[count_local]=me->a;
	    /*	    printf("L_NEW %d %d %7.5g\n",
		count_local,index_local[count_local],
		spars_local[count_local]); */
      //    count_local++;
      //  }
      //me=me->Next;
      //}
      //}  /* end of l  */
      
      for(l=0;l<count_local;l++){          /* i-th row is complete,  */
	index1[count1] = index_local[l]+1;  /* we fill in the global arrays  */
	spars1[count1] = spars_local[l];
	count1++;
	/*printf("## %d %d %d %7.5g\n",l,count1,index1[count1],spars1[count1]);*/
      }
      irwst1[p+1]=count1+1;
      b[p]=s->matrix_RHS[p]+s->matrix_RHS[p+1]+s->matrix_RHS[p+2];
      s->matrix_RHS[p]=0;

    /* we create two new rows, that wi1 = wi0  and wi2=wi0  */
      for(l=1;l<=2;l++){
	p++;
	spars1[count1] = +1;       /* actual element, j=l */
	index1[count1] = p+1;
	count1++;
	spars1[count1] = -1;        /* original element, j=0 */
	index1[count1] = p+1-l;
	count1++;
	irwst1[p+1]=count1+1;
	b[p] = 0;               /* right hand side */
	s->matrix_RHS[p]=0;
      }
    }
    else{  /* element without limiting */
      /* old version of matrix representation */
      /*
      while(me){
	spars[count]=me->a;
	index[count]=me->q+1;
	count++;
	me=me->Next;
	
	//	if(fabs(spars[count-1]-s->matrix[count1]) > 1E-05) {
	//  printf("row %d: (%d %d )%9.5e;  (%d %d)%9.5e;  %9.5e\n",
	// p,count-1,index[count-1],spars[count-1],
	// count1,s->matrix_index[count1]+1,s->matrix[count1],
	// spars[count-1]-s->matrix[count1]);
	//}
	//count1++;
	//while(s->matrix_index[count1]==-1)
	//  count1++;

      }
      irwst[p+1]=count+1;
      */

      /*  NEW ARRAYS  */
      for(j=0;j<s->stencil;j++){
	i=p*s->stencil+j;
	if(s->matrix_index[i] >=0){
	  spars1[count1]=s->matrix[i];
	  index1[count1]=s->matrix_index[i]+1;
	  count1++;
	}
      }

      irwst1[p+1]=count1+1;
      /* right hand side */
      b[p]=s->matrix_RHS[p];
      s->matrix_RHS[p]=0;
    }
  }

  /*
  for(p=0;p<s->rows;p++){
    printf("$$ %d %d %d\n",irwst[p+1],irwst1[p+1],irwst[p+1]-irwst1[p+1]);
    for(i=irwst1[p];i<irwst1[p+1];i++)
      printf("$$ %d:  %d %d %d %7.5g %7.5g %7.5g \n",
	     i,index[i],index1[i],index[i]-index1[i],
	     spars[i], spars1[i], spars[i]-spars1[i]);
  }
  */

  timeb = Time_CPUGet();


  /*  for(p=0;p<nmax;p++) bcnods[p]=1;*/  /* we do not use it */
  for(i=0;i<s->nelem;i++)
    for(j=0;j<s->DOF;j++){
      for(m=0;m<4;m++){
	index2element(p,i,j,m);
	x_init[p]=s->element[i].w[j][m];
      }
    }

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

  /* checking of matrixes */
  /*
     printf("CHECKING OF MATRIXES\n");
  for(p=0;p<nmax;p++){
    int itest=-168;
    double eps=1E-03;
    double sum=0;

    element2index(p,i,j,m);
    
    if(i == itest) printf("... %d %d %d: ",p,irwst1[p],irwst1[p+1]);
    for(q=irwst1[p];q<irwst1[p+1];q++){
      sum+=spars1[q-1]*x_init[index1[q-1]-1];
      if(i == itest) 
	printf("%7.5g,%7.5g(%d,%d) ",
	       spars1[q-1],x_init[index1[q-1]-1],q,index1[q-1]-1);
    }
    if(i == itest)  printf("***\n");
    if(fabs(sum - b[p]) > eps || i==itest){
      element2index(p,i,j,m);
       printf("REZ %d(%d,%d,%d): %7.5g =? %7.5g, reziduum =%7.5g \n",
	       p,i,j,m,sum,b[p],fabs(sum-b[p]) );
           	  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 \n",
		 s->element[i].ve->P[2]->Point[0],
		 s->element[i].ve->P[2]->Point[1]);
      
    }
    } */  /* end of checking of matrices */

  /* calling f77 subroutine */

  time2 = Time_CPUGet();
 
  /*
  for(p=0;p<nmax;p++){
    element2index(p, i, j, m); 
    if( s->element[i].limited == 1){
      printf("... %d %d %d %d %7.5g \n",p,i,j,m,s->element[i].w[j][m]);
    }
  }
  */

  interface_(&nmax, &nzeromax, spars1, irwst1, index1, b,x_init,bcnods,
	     &iterations,&reziduum); 

  time3 = Time_CPUGet();

  s->iterations = iterations;
  s->reziduum = reziduum;

  for(p=0;p<nmax;p++){
    element2index(p, i, j, m); 
    /*
    if(fabs(s->element[i].w[j][m]-x_init[p]) > 1E-06){
      printf("$$ %d %d %d %d %7.5g  %7.5g  %7.5g \n",
	     p,i,j,m,s->element[i].w[j][m],x_init[p],
	     fabs(s->element[i].w[j][m]-x_init[p]));
    }
    */


      s->element[i].w[j][m]=x_init[p]; 


    /*
    if( s->element[i].limited == 1){
      printf("### %d %d %d %d %7.5g \n",p,i,j,m,s->element[i].w[j][m]);
	}
    */
  }


  free(spars);
  free(index);
  free(irwst);
  free(spars1);
  free(index1);
  free(irwst1);
  free(b);
  free(x_init);
  free(bcnods);

  time4 = Time_CPUGet();

  /*    printf("Time costs %5.3g %5.3g %5.3g %5.3g %5.3g | %5.3g %5.3g %5.3g %5.3g %5.3g  \n", */
  /*   printf("Time costs: %5.3g (%5.3g), %5.3g (%5.3g), %5.3g (%5.3g) | %5.3g (%5.3g) \n",
	(time1-time0),(time1-time0)/(time4-time0),
	(time2-time1),(time2-time1)/(time4-time0),
	(time3-time2),(time3-time2)/(time4-time0),
	//           (time4-time3),(time4-time3)/(time4-time0),
	(timeb-timea),(timeb-timea)/(time4-time0) );
	time0 =  Time_CPUGet();
  */
  
  return 0;
}

int Compute_A_Linear(double4 w, matrix4 A1, matrix4 A2, double kappa){
  uint k,l;
  for(k=0;k<4;k++)
    for(l=0;l<4;l++){
      A1[k][l] = 0;
      A2[k][l] = 0;
    }
  A1[0][0] = 1;
  A1[1][1] = 1;
  A1[2][2] = 1;
  A1[3][3] = 1;
  A2[0][0] = 1;
  A2[1][1] = 1;
  A2[2][2] = 1;
  A2[3][3] = 1;
  return 0;
}
int Matrixpm_Linear(double4 w,double volume,double2 n, double length, 
	     matrix4 pp,matrix4 pm, double deltat, double g, double cfl ){
  uint i,k,l;
  for(k=0;k<4;k++)
    for(l=0;l<4;l++){
      pp[k][l] = 0;
      pm[k][l] = 0;
    }

  pp[0][0] = ( 1*n[0]  +  1*n[1])*length;
  pp[1][1] = ( 1*n[0]  +  1*n[1])*length;
  pp[2][2] = ( 1*n[0]  +  1*n[1])*length;
  pp[3][3] = ( 1*n[0]  +  1*n[1])*length;

  for(i=0;i<4;i++){
    pm[i][i]=0;                                                        
    if(pp[i][i] < 0){
      pm[i][i]=pp[i][i];
      pp[i][i]=0;                                 
    }
  }
  return 0;
}

int Compute_A_Init(matrix4 A1, matrix4 A2, double kappa){
  double k1;
  k1=kappa-1;
  

  A1[0][0] = 0;
  A1[0][1] = 1;
  A1[0][2] = 0;
  A1[0][3] = 0;
  A1[1][3] = k1;
  A1[2][3] = 0;

  A2[0][0] = 0;
  A2[0][1] = 0;
  A2[0][2] = 1;
  A2[0][3] = 0;
  A2[1][3] = 0;
  A2[2][3] = k1;
  return 0;
}

int Compute_A_slip(double4 w,double2 n,double length, matrix4 A, double kappa){
  double u,v,rho,e,k3,k1,ker,vvv;
  double2 rn;

  rn[0] = n[0]*length;
  rn[1] = n[1]*length;

  rho = w[0];
  u = w[1]/rho;
  v = w[2]/rho;
  e = w[3];
  k3=kappa-3;
  k1=kappa-1;
  vvv = k1/2*(u*u+v*v);
  ker = kappa*e/rho-vvv;
  
  A[0][0] = 0;
  A[0][1] = rn[0];
  A[0][2] = rn[1];
  A[0][3] = 0;

  A[1][0] = vvv*rn[0];
  A[1][1] = -k3*u*rn[0]+v*rn[1];
  A[1][2] = -k1*v*rn[0]+u*rn[1];
  A[1][3] = rn[0]*k1;

  A[2][0] = vvv*rn[1];
  A[2][1] = v*rn[0]-k1*u*rn[1];
  A[2][2] = u*rn[0]-k3*v*rn[1];
  A[2][3] = rn[1]*k1;

  A[3][0] = 0;
  A[3][1] = ker*rn[0];
  A[3][2] = ker*rn[1];
  A[3][3] = 0;
}

int Compute_A(double4 w, matrix4 A1, matrix4 A2, double kappa){
  double u,v,rho,e,k3,k1,uv,ke,u2,v2;
  rho = w[0];
  u = w[1]/rho;
  v = w[2]/rho;
  e = w[3];
  uv= u*v;
  u2= u*u;
  v2=v*v;
  k3=kappa-3;
  k1=kappa-1;
  ke = kappa*e;
  
  A1[1][0] = k3*u2/2+k1*v2/2;
  A1[1][1] = -k3*u;
  A1[1][2] = -k1*v;

  A1[2][0] = -uv;
  A1[2][1] = v;
  A1[2][2] = u;

  A1[3][0] = -ke*u/rho+k1*u*(u2+v2);
  A1[3][1] = ke/rho-k1*(3*u2+v2)/2;
  A1[3][2] = -k1*uv;
  A1[3][3] = kappa*u;

  A2[1][0] = -uv;
  A2[1][1] = v;
  A2[1][2] = u;

  A2[2][0] = k3*v2/2+k1*u2/2;
  A2[2][1] = -k1*u;
  A2[2][2] = -k3*v;

  A2[3][0] = -ke*v/rho+k1*v*(u2+v2);
  A2[3][1] = -k1*uv;
  A2[3][2] = ke/rho-k1*(u2+3*v2)/2;
  A2[3][3] = kappa*v;
  return 0;
}

int Matrixpm(double4 w,double volume,double2 n, double length, 
	     matrix4 pp,matrix4 pm, double deltat, double g, double cfl ){
  /*-----------------------------------------------------------------------
  calculation of the matrices p-plus (pp) and p-minus(pm)              
  of matrix p, where:                                                  
                  p=n1*a(w)+n2*b(w)                                    
  a,b being jacobian matrices of f,g, respectively.                    
  -----------------------------------------------------------------------*/
  double n1,n2,nv1,nv2,c,h,r,u,v,p,c2,c22,g1,vvv,vvvv,vmax,dtx,del,dt;
  matrix4 t,t1;
  uint i,k,l;
  double4 dp,dm;

  /*  printf("TTT. %7.5g %7.5g %7.5g %7.5g \n",w[0],w[1],w[2],w[3]); 
  printf("TTT. %7.5g %7.5g %7.5g %7.5g \n",length,n[0],n[1],volume); */

  /*      common /const/ g, dt, cfl*/
  /*r=w[0];
  u=w[1]/w[0];
  v=w[2]/w[0];
  p=(g-1.)*(w[3]-0.5*r*(u*u+v*v)); */
  
  r=w[0];
  u=w[1];
  v=w[2];
  p=w[3]; 
  //  p=w[3]*w[3]*w[0]/g;
  
                                                                   
  /* r-density, u,v-velocity, p-pressure, g-poisson constant              
     c-the local speed of sound, h-enthalpy */
                                                                       
  if( g*p/r <0 ) return 1;
  c=sqrt(g*p/r);                                    
  h=c*c/(g-1.)+(u*u+v*v)/2.;                                         
  c2=c*c;                                            
  c22=2.*c2;
  g1=g-1.;
  vvv= length;
  nv1=n[0];
  nv2=n[1];
  n1 = nv1*length;
  n2 = nv2*length;

  /*---t and t1 - transformation t--------  */

  t[0][0]=1.;                                                        
  t[0][1]=0.;                                                         
  t[0][2]=1./c22;                                                     
  t[0][3]=1./c22;                                                     
  t[1][0]=u;                                                      
  t[1][1]=nv2;                                                        
  t[1][2]=(u+c*nv1)/c22;                                              
  t[1][3]=(u-c*nv1)/c22;                                              
  t[2][0]=v;                                                          
  t[2][1]=-nv1;                                                       
  t[2][2]=(v+c*nv2)/c22;
  t[2][3]=(v-c*nv2)/c22;                                              
  t[3][0]=(u*u+v*v)/2.;                                               
  t[3][1]=nv2*u-nv1*v;                                                
  t[3][2]=(h+c*(nv1*u+nv2*v))/c22;
  t[3][3]=(h-c*(nv1*u+nv2*v))/c22;                                    
                                                                       
  t1[0][0]=1-g1/c2*(u*u/2+v*v/2);                                     
  t1[0][1]=g1/c2*u;                                          
  t1[0][2]=g1/c2*v;                                                   
  t1[0][3]=-g1/c2;
  t1[1][0]=nv1*v-nv2*u;                                               
  t1[1][1]=nv2;                                                       
  t1[1][2]=-nv1;                                                      
  t1[1][3]=0.;                                                        
  t1[2][0]=-c*(nv1*u+nv2*v)+g1*(u*u/2+v*v/2);
  t1[2][1]=c*nv1-g1*u;                            
  t1[2][2]=c*nv2-g1*v;                                                
  t1[2][3]=g1;                                                        
  t1[3][0]=c*(nv1*u+nv2*v)+g1*(u*u/2+v*v/2);
  t1[3][1]=-c*nv1-g1*u;                             
  t1[3][2]=-c*nv2-g1*v;                                               
  t1[3][3]=g1;

  /*  printf("......\n");
   for(k=0;k<4;k++){
    for(l=0;l<4;l++)
      printf("%7.5g,",t[k][l]);
    printf("       ");
    for(l=0;l<4;l++)
      printf("%7.5g,",t1[k][l]);
    printf("\n");

  } */

  /*---dp ... eigenvalues of the matrix n1*a(w)+n2*b(w)---------- */

  dp[0]=n1*u+n2*v;    
  dp[1]=dp[0];                                                      
  dp[2]=dp[0]+c*vvv;
  dp[3]=dp[0]-c*vvv;


  /* --changing dt according to cfl-condition----  */

  /*  vmax=0;
  for(i=0;i<4;i++)
    if(fabs(dp[i])>vmax) vmax=fabs(dp[i]);

  vvv=vmax*dt/volume;
  if(vvv>cfl){                                                
    dtx = volume/vmax * cfl;
    dt = min( dt, dtx );                                           
  }
  vvvv = vmax*deltat/volume;
  if( vvvv > cfl ){
    del = volume / vmax * cfl;
    deltat = min( deltat, del );
  }                                  */                         

  /*---calculating the negative (dm) and the positive(dp)                  
       parts of eigenvalues-----------------------------------  */

  for(i=0;i<4;i++){
    dm[i]=0;                                                        
    if(dp[i] < 0){
      dm[i]=dp[i];
      dp[i]=0;                                 
    }
  }
  /*  printf("DP %7.5g %7.5g %7.5g %7.5g \n",dp[0],dp[1],dp[2],dp[3]);
  printf("DM %7.5g %7.5g %7.5g %7.5g \n",dm[0],dm[1],dm[2],dm[3]); 
   for(k=0;k<4;k++){
    for(l=0;l<4;l++)
      printf("%7.5g,",pp[k][l]);
    printf("       ");
    for(l=0;l<4;l++)
      printf("%7.5g,",pm[k][l]);
    printf("\n");
  } */

  /*---multiplication of matrices: pp=t*dp*t1, where dp is diagonal...-- */
  Threem(t,dp,t1,pp); 
  Threem(t,dm,t1,pm);

  /*   for(k=0;k<4;k++){
    for(l=0;l<4;l++)
      printf("%7.5g,",pp[k][l]);
    printf("       ");
    for(l=0;l<4;l++)
      printf("%7.5g,",pm[k][l]);
    printf("\n");
   } 
   printf("-----------------------\n"); */
   return 0;
}

int Threem(matrix4 a,double4 b,matrix4 c,matrix4 d){
  /*-----------------------------------------------------------------------
    matrices multiplication d=a*b*c, where b is diagonal, in fact        
    b is a vector consisting of diagonal elements                        
    ----------------------------------------------------------------------- */
  uint i,j,k;
  
  for(i=0;i<4;i++)
  for(j=0;j<4;j++){
    d[i][j]=0;
    for(k=0;k<4;k++)
      d[i][j]+=a[i][k]*b[k]*c[k][j];
  }
}

int Set_Slip(tSolver_State *s,uint i, uint _i, uint z,
	      double p,double4 wq, double4 wj, double24 Fb ){
  /* i ..index of edge, _i .. index of element */
  uint j;
  int ic;
  double n1,n2, v2;

  v2 = wq[1]*wq[1]+wq[2]*wq[2];

  ic=s->element[_i].icve;	
  if(ic >= 0){  /*curved elements  */
    n1 = s->curvedelement[ic].normal[z][0]*s->curvedelement[ic].length[z];
    n2 = s->curvedelement[ic].normal[z][1]*s->curvedelement[ic].length[z];
  }
  else{   /*non curved elements */
    n1 = s->face[i].n[0]*s->face[i].volume;
    n2 = s->face[i].n[1]*s->face[i].volume;
  }

  wj[0]=0;
  wj[1]=p*n1; 
  wj[2]=p*n2;
  wj[3]=0;
  n1*=s->kappa-1;
  n2*=s->kappa-1;
  Fb[0][0] = v2*n1/2;
  Fb[0][1] = -wq[1]*n1;
  Fb[0][2] = -wq[2]*n1;
  Fb[0][3] = n1;
  Fb[1][0] = v2*n2/2;
  Fb[1][1] = -wq[1]*n2;
  Fb[1][2] = -wq[2]*n2;
  Fb[1][3] = n2;

  return 1;
}

void Set_Boundary_State(tSolver_State *s,uint i, uint _i,
			double4 wi,double4 wj){
  /* i ..index of edge, _i .. index of element */
  uint j;
  int ic;
  double press;
  double amplitude=0.02, omega = 6.283185307e-00/2;

  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 */
    if(s->boundary[s->face[i].boundary-1].fixed){  /* fixed */
      for(j=0;j<4;j++)
	wj[j]=s->boundary[s->face[i].boundary-1].w[j];
    }
    else{  /* physical  */
      double oneoverrho,norm,c2,product,product2;
      oneoverrho=1/wi[0];
      norm=wi[1]*wi[1]+wi[2]*wi[2];
      c2=s->kappa2*oneoverrho*(wi[3]-0.5*oneoverrho*norm);
      product=oneoverrho*(s->face[i].n[0]*wi[1]+s->face[i].n[1]*wi[2]);
      product2=product*product;
      if(0<product&&c2<=product2)  /* supersonic outlet */
	for(j=0;j<4;j++)
	  wj[j] = wi[j];
      else
	if(0<product){  /* subsonic outlet */
	  for(j=0;j<4;j++)
	    wj[j] = wi[j];

	  /* ORIGINAL */
	  /*
	  wj[3]=s->boundary[s->face[i].boundary-1].w[3]+
	    0.5*((wj[1]*wj[1]+wj[2]*wj[2])/wj[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]);
	  */
	  press=(s->boundary[s->face[i].boundary-1].w[3]-
		 0.5*((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]))*(s->kappa-1);
	  
	  /* GAMM channel, time dependent output !!!!!!!!!!!!!!!  */

	  press*=(1+amplitude*sin(omega*s->t));

	  /*  !!!!!!!!!!!!!!! */
	  wj[3]=press/(s->kappa-1)+0.5*(wj[1]*wj[1]+wj[2]*wj[2])/wj[0];
	  /*printf("Outlet: time=%7.5g,  press=%7.5g,  ener=%7.5g\n", 
		 s->t,press,wj[3]);
	  */
	}
	else
	  if(product2<c2){  /* subsonic inlet */
	    for(j=0;j<3;j++)
	      wj[j]=s->boundary[s->face[i].boundary-1].w[j];
	    
	    wj[3]=wi[3]+ 0.5*((wj[1]*wj[1] + wj[2]*wj[2])/wj[0]-
			       (wi[1]*wi[1] + wi[2]*wi[2])/wi[0]);

	  }
	  else  /* supersonic inlet */
	    for(j=0;j<4;j++)
	      wj[j]=s->boundary[s->face[i].boundary-1].w[j];

    }
  }
  else{  /* noslip condition  */
    for(j=0;j<4;j++)
      wj[j] = wi[j];
  }
}

  
int Matrix_Times_Vector44(matrix4 A,double4 x, double4 b){
uint i,j;

 for(i=0;i<4;i++) b[i]=0;
 for(i=0;i<4;i++)
   for(j=0;j<4;j++)
     b[i]+=A[i][j]*x[j];
 return 0;
}


int Transform_w2q(double4 wi, double2 n, double kappa, double4 w, double4 q){
  double u,v,c,uv2,np,nm,upc,umc,vpc,vmc,uc,vc,kappa1,ckapp;

  /* transform w->q, to a system, which is diagonal */
  /*  tranform to local coordinate system */
  u = (wi[1]*n[0] + wi[2]*n[1])/wi[0];
  v = (wi[2]*n[0] - wi[1]*n[1])/wi[0];
  uv2 = (u*u+v*v)/2;

  kappa1=kappa-1;
  c=sqrt(kappa*kappa1*(wi[3]/wi[0]-uv2));

  uc = u*c;
  vc = v*c;
  np = n[0]+n[1];
  nm = n[0]-n[1];
  upc= u+c;
  umc= u-c;
  vpc= v+c;
  vmc= v-c;
  ckapp = c*c/kappa1+uv2;

  q[0] = w[0] + nm*w[1] + np*w[2] + w[3];
  q[1] = umc*w[0] + u*nm*w[1] + u*np*w[2] + upc*w[3];
  q[2] = v*w[0] + (vmc*n[0]-vpc*n[1])*w[1] + (vmc*n[1]+vpc*n[0])*w[2] + v*w[3];
  q[3] = (ckapp-uc)*w[0] + (uv2*nm-vc*np)*w[1] + (uv2*np+vc*nm)*w[2] +
    (ckapp+uc)*w[3];
}

int Transform_q2w(double4 wi, double2 n, double kappa, double4 q, double4 w){
  double u,v,c,c2,uv2,np,nm,uc,vc,kappa1,ckapp;
  int i;
  /* reverse transformation */
  /*  tranform to local coordinate system */
  u = (wi[1]*n[0] + wi[2]*n[1])/wi[0];
  v = (wi[2]*n[0] - wi[1]*n[1])/wi[0];
  uv2 = (u*u+v*v)/2;

  kappa1=kappa-1;
  c=sqrt(kappa*kappa1*(wi[3]/wi[0]-uv2));

  uc = u*c;
  vc = v*c;
  c2 = c*c;
  np = n[0]+n[1];
  nm = n[0]-n[1];
  ckapp = c*c/kappa1+uv2;

  w[0] = (uv2*kappa1+uc)*q[0] - (c+u*kappa1)*q[1] - v*kappa1*q[2] + kappa1*q[3];
  w[1] = (vc*np-(kappa1*uv2-c2)*nm)*q[0] + u*nm*kappa1*q[1] + 
    (v*kappa1*nm-np*c)*q[2] - kappa1*nm*q[3];
  w[2] = -(np*(kappa1*uv2-c2)+vc*nm)*q[0] + u*kappa1*np*q[1] + 
    (np*v*kappa1+nm*c)*q[2] - np*kappa1*q[3];
  w[3] = (uv2*kappa1-uc)*q[0] + (c-u*kappa1)*q[1] - v*kappa1*q[2] + kappa1*q[3];
  
  for(i=0;i<4;i++)
    w[i]*=0.5/c2;
}

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

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

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

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

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

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

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

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

  gb    =  gamma - 1.;

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

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

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

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

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

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