#include "embed.h"
#include "navier-stokes/centered.h"
#include "fractions.h"
#include "PointTriangle.h"
#include "view.h"

u.n[embed] = dirichlet (0.);
u.t[embed] = dirichlet (0.);

double muv = 0.0001;

face vector muc[];

double angle = 20; //channel Angle
double L = 2.8;    //Length of a channel
double R = 0.3;    //Radius of cirle
double sect = 25;  //sections in circle
 

int get_coords_Tesla (coord * p) {
  double dydx = tan(pi*angle/180.);
  double O = R/(dydx); 
  double y1 = dydx*L/2, y2 = dydx * (L/2 - O);
  int c = 0;
  double start = -2;
  for (int j = 0; j < 5; j++) {
    for (int l = 1; l > -2; l -= 2) {
      p[c++] = (coord){start, 0.};
      p[c++] = (coord){start + L/2, l*y1};
      start += L/2 - O;
      p[c++] = (coord){start, l*y2};
      p[c++] = (coord){start + y2/dydx, 0.};
      for (double theta = 0; theta < pi; theta += pi/sect + 1e-3) {
	double xr = start + O, yr = l*(y1 - R);
	p[c++] = (coord){xr + R*sin(theta), yr + R*cos(theta)};
	p[c++] = (coord){xr + R*sin(theta + pi/sect), yr + R*cos(theta + pi/sect)};
      }
      start += L/2 - O;
    }
  }
  p[c].x = nodata;
  return c;
}

int main() {
  periodic (left);
  double dydx = tan(pi*angle/180.);
  double O = R/(dydx); 
  L0 = 4*(L - 2*O);
  printf ("%g %g\n", L0, O);
  N = 256;
  Y0 = -2;
  const face vector av[] = {0.01, 0};
  a = av;
  mu = muc;
  run();
  const face vector av2[] = {-0.01, 0};
  a = av2;
  run();
}

event properties (i++) {
  foreach_face()
    muc.x[] = muv*fs.x[];
  boundary((scalar*){muc});
}

event init (t = 0) {
  coord p[999];
  get_coords_Tesla(p);
  vertex scalar d[];
  foreach_vertex() {
    double min = HUGE;
    coord * a = p;
    coord vert = {x, y};
    while (a->x != nodata) {
      coord r;
      double d, d2 = PointSegmentDistance (&vert, a, a + 1, &r, &d);
      if (d2 < min)
	min = d2;
      a += 2;
    }
    d[] = min > 0 ? -sqrt (min) : -100;
  }  fractions (d, cs, fs, -0.1);
}

event adapt (i++) {
  adapt_wavelet({cs, u}, (double[]){1e-3, 1e-2, 1e-2}, 8);
}

event movie (t += 0.1) {
  scalar U[];
  foreach() 
    U[] = sqrt(sq(u.x[]) + sq(u.y[]));
  boundary ({U});
  view (fov = 10, tx = -0.5,
	height = 300, width = 600);
  draw_vof ("cs", "fs", filled = -1, fc = {0.6, 0.6, 0.6});
  draw_vof ("cs", "fs");
  squares ("U", min = 0, max = 3e-1);
  if (constant(a.x) > 0)
    draw_string ("  Flow ->", 1, lw = 4, size = 30);
  else
    draw_string ("  <- Flow", 1, lw = 4, size = 30);
  save ("tesla.mp4");
  if (fabs(t - 49.9) < 1e-5 && constant(a.x) > 0) 
    save ("tesla.png");
}

event flow_rate (i++) {
  double flx = 0;
  foreach_face(x) {
    if (fabs(x - X0) < 1e-5) 
      flx += uf.x[]*fm.x[]*Delta;
  }
  if (constant(a.x) > 0)
    printf ("%g %g %g\n", t, flx, statsf(u.x).sum);
  else
    fprintf (stderr, "%g %g %g\n", t, flx, statsf(u.x).sum);
}

event stop (t = 50);