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){
      int neigh,j,i1,j1;
      double gi,gj,wiL,wiR,wjL,wjR,aver;
      BOOL limit;
      uint count=0;

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

      for(i1=0;i1<3;i1++){
	gi = 2*s->element[i].w[perm[i1][0]][0] -
	  s->element[i].w[perm[i1][1]][0] - s->element[i].w[perm[i1][2]][0];
	j = s->element[i].neigh[i1];
	j1 = s->element[i].neighind[i1];
	gj = 2*s->element[j].w[perm[j1][0]][0] -
	  s->element[j].w[perm[j1][1]][0] - s->element[j].w[perm[j1][2]][0];
	if(gi*gj<0){
	  wiL = s->element[i].w[perm[i1][0]][0] +
	  s->element[i].w[perm[i1][1]][0] - s->element[i].w[perm[i1][2]][0];
	  wiR = s->element[i].w[perm[i1][0]][0] +
	  s->element[i].w[perm[i1][2]][0] - s->element[i].w[perm[i1][1]][0];
	  wjR = s->element[j].w[perm[j1][0]][0] +
	  s->element[j].w[perm[j1][1]][0] - s->element[j].w[perm[j1][2]][0];
	  wjL = s->element[j].w[perm[j1][0]][0] +
	  s->element[j].w[perm[j1][2]][0] - s->element[j].w[perm[j1][1]][0];
	  if(gi>0 && wiL<=wjL && wiR <=wjR ) count++;
	  if(gi<0 && wiL>=wjL && wiR >=wjR ) count++;
	}
      }
      if(count == 3) limit = FALSE;

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

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