struct Dissipation {
  scalar dis;
  vector u;
  face vector mu;
#if COMPRESSIBLE
  face vector lambdav;
#endif
};

#if AXI
# define gamma (1./y)
#else
# define gamma (1.)
#endif

void dissipation (struct Dissipation p)
{
  vector u = p.u; 
  scalar dis = p.dis;
  (const) face vector mu = p.mu;

#if COMPRESSIBLE
  (const) face vector lambdav = p.lambdav;
#endif

#if AXI
  scalar ur = u.y;
#endif
#if TREE
  /* conservative coarse/fine discretisation (2nd order) */

  foreach () {
#if AXI
    dis[] = ((mu.x[] + mu.x[1] + mu.y[] + mu.y[0,1])/2.*gamma*sq(ur[]/y)
#if COMPRESSIBLE
	     + (lambdav.x[] + lambdav.x[1] +lambdav.y[]+lambdav.y[0,1])/4.*
	     (ur[]/y + (u.x[1] -u.x[-1] + u.y[0,1] -u.y[0,-1])/2./Delta)*ur[]/y
#endif
	     );
#else
    dis[] = 0.;
#endif
  }

  foreach_dimension() {
    face vector taux[];
#if COMPRESSIBLE
    face vector tauc[];
#if AXI
    face vector axic[];
#endif
#endif
    foreach_face(x) {
#if COMPRESSIBLE
#if AXI
      axic.x[] = lambdav.x[]*(ur[]+ur[-1])
	*(u.x[] - u.x[-1])/Delta/2.;
#endif
      tauc.x[] = lambdav.x[]*(u.x[] - u.x[-1]
     #if dimension > 1
			      + (u.y[0,1] + u.y[-1,1])/4
			      - (u.y[0,-1] + u.y[-1,-1])/4.
    #endif
    #if dimension > 2
			      + (u.z[0,0,1] + u.z[-1,0,1])/4
			      - (u.z[0,0,-1] + u.z[-1,0,-1])/4.
    #endif
			      )*(u.x[] - u.x[-1])/sq(Delta);
#endif
      taux.x[] =  2.*mu.x[]*sq((u.x[] - u.x[-1])/Delta);
    }
    #if dimension > 1
      foreach_face(y)
	taux.y[] = mu.y[]*(u.x[] - u.x[0,-1] +
			   (u.y[1,-1] + u.y[1,0])/4. -
			   (u.y[-1,-1] + u.y[-1,0])/4.)
	*(u.x[] - u.x[0,-1])/sq(Delta);
    #endif
    #if dimension > 2
      foreach_face(z)
	taux.z[] = mu.z[]*(u.x[] - u.x[0,0,-1] +
			   (u.z[1,0,-1] + u.z[1,0,0])/4. -
			   (u.z[-1,0,-1] + u.z[-1,0,0])/4.)
	*(u.x[] - u.x[0,0,-1])/sq(Delta);
    #endif  
      foreach () {
        double d = 0;
	foreach_dimension()
	  d += (taux.x[1] + taux.x[])*gamma;
	dis[] += (d
#if COMPRESSIBLE 
               + tauc.x[1] + tauc.x[]
#if AXI
               + (axic.x[1] + axic.x[])/y
#endif
#endif
	       )/2.;
      }
  }
#else
/*   /\* "naive" discretisation (only 1st order on trees) *\/ */
  foreach () {
#if AXI
    dis[] = ((mu.x[] + mu.x[1] + mu.y[] + mu.y[0,1])/2.*gamma*sq(ur[]/y)
#if COMPRESSIBLE
	     + (lambdav.x[] + lambdav.x[1] +lambdav.y[]+lambdav.y[0,1])/4.*
	     (ur[]/y + (u.x[1] -u.x[-1] + u.y[0,1] -u.y[0,-1])/2./Delta)*ur[]/y
#endif
	     );
#else
    dis[] = 0.;
#endif

    foreach_dimension() 
	     dis[] += ((mu.x[1,0]*sq(u.x[1] - u.x[])
		+ mu.x[]*sq(u.x[] - u.x[-1])
      #if dimension > 1
		+ mu.y[0,1]*(u.x[0,1] - u.x[] +
		(u.y[1,0] + u.y[1,1])/4. -
		(u.y[-1,0] + u.y[-1,1])/4.)*(u.x[0,1] - u.x[])/2.
		+ mu.y[]*(u.x[] - u.x[0,-1] +
		(u.y[1,-1] + u.y[1,0])/4. -
		(u.y[-1,-1] + u.y[-1,0])/4.)*(u.x[] - u.x[0,-1])/2.
      #endif
      #if dimension > 2
		+ mu.z[0,0,1]*(u.x[0,0,1] - u.x[] +
			       (u.z[1,0,0] + u.z[1,0,1])/4. -
			       (u.z[-1,0,0] + u.z[-1,0,1])/4.)*(u.x[0,0,1] - u.x[])/2.
		- mu.z[]*(u.x[] - u.x[0,0,-1] +
			  (u.z[1,0,-1] + u.z[1,0,0])/4. -
			  (u.z[-1,0,-1] + u.z[-1,0,0])/4.)*(u.x[] - u.x[0,0,-1])/2.
       #endif
		)*gamma
#if COMPRESSIBLE
	     + lambdav.x[1]*sq(u.x[1] - u.x[])/2.
	     + lambdav.x[]*sq(u.x[] - u.x[-1])/2.
               #if dimension > 1
	     + lambdav.x[1]*((u.y[1,1] + u.y[0,1])/4 -
	                     (u.y[1,-1] + u.y[0,-1])/4.)*(u.x[1] - u.x[])/2.
	     + lambdav.x[]*((u.y[0,1] + u.y[-1,1])/4 -
			    (u.y[0,-1] + u.y[-1,-1])/4.)*(u.x[] - u.x[-1])/2.
#if AXI
	     + lambdav.x[1]*(ur[1] + ur[])*(u.x[1] - u.x[])/4./y*Delta
	     + lambdav.x[]*(ur[-1] + ur[])*(u.x[] - u.x[-1])/4./y*Delta
#endif
               #endif
               #if dimension > 2
		  + lambdav.x[1]*((u.z[1,0,1] + u.z[0,0,1])/4 -
				  (u.z[1,0,-1] + u.z[0,0,-1])/4.)*(u.x[1] - u.x[])/2.
		  + lambdav.x[]*((u.z[0,0,1] + u.z[-1,0,1])/4 -
				    (u.z[0,0,-1] + u.z[-1,0,-1])/4.)*(u.x[] - u.x[-1])/2.
              #endif
#endif
	     )/sq(Delta);
}
#endif
}	     
  
#undef gamma