#include "nsf4t.h"
scalar * tracers = NULL;
#define funcux (-y*exp(-sq(x) - sq(y)))
#define funcuy ( x*exp(-sq(x) - sq(y)))

double funu_x (double x, double y) {
  return funcux;
}

double funu_y (double x, double y) {
  return funcuy;
}

double dfunux_x (double x, double y) {
  return 2.*x*y*exp(-sq(x) - sq(y));
}
double dfunux_y (double x, double y) {
  return (2.*sq(y) - 1)*exp(-sq(x) - sq(y));
}

double dfunuy_x (double x, double y) {
  return (1 - 2.*sq(x))*exp(-sq(x) - sq(y));
}

double dfunuy_y (double x, double y) {
  return -2*x*y*exp(-sq(x) - sq(y));
}

double d2funux_x (double x, double y) {
  return (2*y - 4.*sq(x)*y)*exp(-sq(x) - sq(y));
}

double d2funux_y (double x, double y) {
  return (6*y - 4.*cube(y))*exp(-sq(x) - sq(y));
}

double d2funuy_x (double x, double y) {
  return (4.*cube(x) - 6*x)*exp(-sq(x) - sq(y));
}

double d2funuy_y (double x, double y) {
  return (4.*sq(y)*x - 2*x)*exp(-sq(x) - sq(y));
}

double pressure (double x, double y) {
  return exp(-2*(sq(x) + sq(y)))/4.;
}

double dp_x (double x, double y) {
  return -x*exp(-2*(sq(x) + sq(y)));
}

double dp_y (double x, double y) {
  return -y*exp(-2*(sq(x) + sq(y)));
}

double ue, expo = 1.;
int main() {
  L0 = 15.;
  X0 = -L0/2. + pi/25.;
  Y0 = -L0/2. - exp(1.)/25.;
  N = 8;
  run();
}

event check_errors (i = 0) {
  for (ue = 1e-1; ue > 1e-9; ue /= 4) {
    do {
      foreach_face() 
	u.x[] = Gauss6_x (x, y, Delta, funu_x);
      boundary ((scalar*){u});
    } while (adapt_flow (ue, 99, expo).nf > 1);
    double e1 = 0, einf = -1;
    vector v[];
    foreach_dimension() {
      v.x.prolongation = refine_vert5;
      v.x.restriction = restriction_vert;
    }
    // Face to vertex
    foreach() {
      foreach_dimension() {
	v.x[] = FACE_TO_VERTEX_4(u.x);
	double el = fabs(funu_x(x - Delta/2, y - Delta/2) - v.x[]);
	e1 += el*sq(Delta);
	if (einf < el)
	  einf = el;
	v.x[] = funu_x(x - Delta/2, y - Delta/2); //...
      }
    }
    boundary ((scalar*){v});
    vector * grads = NULL;
    for (scalar s in {v}) {
    vector dsd = new_vector ("gradient");
    grads = vectors_add (grads, dsd);
    }
    for (vector grad in grads) {
      foreach_dimension() {
	grad.x.prolongation = refine_vert5;
	grad.x.restriction = restriction_vert;
      }
    }
    compact_upwind ((scalar*){v}, grads, v);
    double u1d = 0, uinfd = -1;
    // Gradients
    foreach() {
      vector gradux = grads[0]; 
      vector graduy = grads[1]; 
      foreach_dimension() {
	double el = fabs (gradux.y[] -
			  dfunux_y (x - Delta/2., y - Delta/2));
	
	u1d += sq(Delta)*el;
	if (uinfd < el) {
	  uinfd = el;
	}
	el = fabs (graduy.y[] -
		   dfunuy_y (x - Delta/2., y - Delta/2));
	u1d += sq(Delta)*el;
	if (uinfd < el) {
	  uinfd = el;
	}
      }
    }
    double u1d2 = 0, uinfd2 = -1;
    // double derivatives
    foreach() {
      scalar s = v.x;
      foreach_dimension() {
	double el = fabs (D2SDX2/sq(Delta) -
			  d2funux_x (x - Delta/2., y - Delta/2.));
	u1d2 += sq(Delta)*el;
	if (uinfd2 < el)
	  uinfd2 = el;
      }
      s = v.y;
      foreach_dimension() {
	double el = fabs (D2SDX2/sq(Delta) -
			  d2funuy_x (x - Delta/2., y - Delta/2.));
	u1d2 += sq(Delta)*el;
	if (uinfd2 < el)
	  uinfd2 = el;
      }
    }
    double u1r = 0, uinfr = -1;
    // Vertex to face
    foreach_face() {
      double el =  fabs (Gauss6_x(x, y, Delta, funu_x)
			 - VERTEX_TO_FACE_4(v.x));
      u1r += sq(Delta)*el;
      if (uinfr < el)
	uinfr = el;
    }
    // p.prolongation = refine_5th; for 4th order L_{inf}
    foreach()
      p[] = Gauss6 (x, y , Delta, pressure);
    boundary ({p});
    double uf1 = 0, ufinf = -1;
    // Face gradient
    foreach_face() {
      double el = fabs ((p[-2,0] - 15.*p[-1,0] +
			 15.*p[0,0] - p[1,0])/(12*Delta) -
			Gauss6_x (x, y, Delta, dp_x));
      uf1 += sq(Delta)*el;
      if (ufinf < el)
	ufinf = el;
    }
    free (grads); grads = NULL;
    fprintf (stderr,"%g %g %g %g %g %g %g %g %g %g %g\n",sqrt(grid->tn),
	     e1, einf, u1d, uinfd, u1d2, uinfd2, u1r, uinfr, uf1, ufinf);
  }
  scalar lev[];
  foreach()
    lev[] = level;
  output_ppm (lev, file = "lev.png", n = 350, max = depth() + 0.5);
  return 1;
}