#define W 5.0
#define A 0.25

#define NX 1.
#define NY 1.

#include "embed.h"
#include "src/embed_utils.h"
#include "run.h"

scalar phi[], rho = unity;
face vector a[];
(const) face vector alpha = unityf, epsilon = unityf;

#include "src/stress.h"

static double exact (double x, double y)
{ 
  return A*sin(2.*NX*pi/W*x)*sin(2.*NY*pi/W*y);
}

#define FX(x,y) -4.*sq(A)*NX*(sq(NX)+sq(NY))*sq(sin(2.*NY*pi/W*y))*	\
			 sin(4.*NX*pi/W*x)*pow(pi/W,3.)
#define FY(x,y) -4.*sq(A)*NY*(sq(NX)+sq(NY))*sq(sin(2.*NX*pi/W*x))*	\
			 sin(4.*NY*pi/W*y)*pow(pi/W,3.)


phi[embed] = dirichlet(exact(x,y));
phi[left] = exact(x-Delta/2.,y);
phi[right] = exact(x+Delta/2.,y);
phi[top] = exact(x,y + Delta/2.);
phi[bottom] = exact(x,y - Delta/2.);


int main()
{
  size(W);
  origin(-W/2., -W/2.);
  for (int lvl = 5; lvl <= 8; lvl++) {
    N = 1 << lvl;
    run();
  }
}  

event defaults (i = 0) {
  vertex scalar phiv[];
  foreach_vertex()
    phiv[] = sq(x) + sq(y) - sq(1.1);
  
  fractions (phiv, cs, fs);

  epsilon = fs;
  alpha = fs;
  rho = cs;
  foreach() {
    double xc = x, yc = y;
#if 0    
      if (cs[] > 0. && cs[] < 1.) {
	coord n = facet_normal (point, cs, fs), p;
	double alpha = plane_alpha (cs[], n);
	line_center (n, alpha, cs[], &p);
	xc += p.x*Delta, yc += p.y*Delta;
      }
#endif      
    phi[] = cs[] > 0. ? exact(xc, yc) : nodata;
  }

  //  output_cells(stdout);
  //  output_facets(cs,stdout,fs);
}

event dibujar (i = 0) {
  face vector err[];

  foreach_face (x) {
    err.x[] = a.x[] - FX(x,y);
    printf ("x: %g %g %g %g %g  \n", x, y, a.x[], FX(x,y), err.x[]);
  }
     
  foreach_face (y) {
    err.y[] = a.y[] - FY(x,y);
    printf ("y: %g %g %g %g %g \n", x, y, a.y[], FY(x,y),err.y[]);
  }
  
  norm nx = normf (err.x), ny = normf (err.y);
     
  fprintf (stderr, "%d %g %g %g %g %g %g\n", N,
	   nx.avg, nx.rms, nx.max,
	   ny.avg, ny.rms, ny.max);

}