/**
# Basic geometric functions

These basic geometric functions are mostly related to Volume-Of-Fluid
computations.

We consider a square cell of size unity centered on the origin, cut by
a straight line.

![Cell and interface](/src/figures/square.svg)

The line can be described by the equation
$$
n_xx+n_yy=\alpha
$$
where $\mathbf{n}$ is a vector normal to the interface and $\alpha$ is
the intercept. We note $c$ the volume of the part of the square cell
which lies "inside" the interface, where "inside" is defined by
convention as the opposite direction to the normal vector $\mathbf{n}$
(i.e. the normal vector is pointing "outside").

With these definitions, the interface is uniquely defined by providing
$\mathbf{n}$ and either $\alpha$ or $c$ i.e. there is a unique
function which computes $\alpha$ given $c$ and $\mathbf{n}$. We call
this function `line_alpha()` and define it as: */

#if dimension == 1
double line_alpha (double c, coord n)
{
  return clamp (c, 0., 1.) - 0.5;
}
#endif // dimension == 1

#if dimension >= 2
double line_alpha (double c, coord n)
{
  double alpha, n1, n2;
  
  n1 = fabs (n.x); n2 = fabs (n.y);
  if (n1 > n2)
    swap (double, n1, n2);

  c = clamp (c, 0., 1.);
  double v1 = n1/2.;
  if (c <= v1/n2)
    alpha = sqrt (2.*c*n1*n2);
  else if (c <= 1. - v1/n2)
    alpha = c*n2 + v1;
  else
    alpha = n1 + n2 - sqrt (2.*n1*n2*(1. - c));

  if (n.x < 0.)
    alpha += n.x;
  if (n.y < 0.)
    alpha += n.y;

  return alpha - (n.x + n.y)/2.;
}
#endif // dimension >= 2

#if dimension >= 3
double plane_alpha (double c, coord n)
{
  double alpha;
  coord n1;
  
  n1.x = fabs (n.x); n1.y = fabs (n.y); n1.z = fabs (n.z);

  double m1, m2, m3;
  m1 = min(n1.x, n1.y);
  m3 = max(n1.x, n1.y);
  m2 = n1.z;
  if (m2 < m1) {
    double tmp = m1;
    m1 = m2;
    m2 = tmp;
  }
  else if (m2 > m3) {
    double tmp = m3;
    m3 = m2;
    m2 = tmp;
  }
  double m12 = m1 + m2;
  double pr = max(6.*m1*m2*m3, 1e-50);
  double V1 = m1*m1*m1/pr;
  double V2 = V1 + (m2 - m1)/(2.*m3), V3;
  double mm;
  if (m3 < m12) {
    mm = m3;
    V3 = (m3*m3*(3.*m12 - m3) + m1*m1*(m1 - 3.*m3) + m2*m2*(m2 - 3.*m3))/pr;
  }
  else {
    mm = m12;
    V3 = mm/(2.*m3);
  }

  c = clamp (c, 0., 1.);
  double ch = min(c, 1. - c);
  if (ch < V1)
    alpha = pow (pr*ch, 1./3.);
  else if (ch < V2)
    alpha = (m1 + sqrt(m1*m1 + 8.*m2*m3*(ch - V1)))/2.;
  else if (ch < V3) {
    double p12 = sqrt (2.*m1*m2);
    double q = 3.*(m12 - 2.*m3*ch)/(4.*p12);
    double teta = acos(clamp(q,-1.,1.))/3.;
    double cs = cos(teta);
    alpha = p12*(sqrt(3.*(1. - cs*cs)) - cs) + m12;
  }
  else if (m12 <= m3)
    alpha = m3*ch + mm/2.;
  else {
    double p = m1*(m2 + m3) + m2*m3 - 1./4., p12 = sqrt(p);
    double q = 3.*m1*m2*m3*(1./2. - ch)/(2.*p*p12);
    double teta = acos(clamp(q,-1.,1.))/3.;
    double cs = cos(teta);
    alpha = p12*(sqrt(3.*(1. - cs*cs)) - cs) + 1./2.;
  }
  if (c > 1./2.) alpha = 1. - alpha;

  if (n.x < 0.)
    alpha += n.x;
  if (n.y < 0.)
    alpha += n.y;
  if (n.z < 0.)
    alpha += n.z;

  return alpha - (n.x + n.y + n.z)/2.;;
}
#else // dimension < 3
# define plane_alpha line_alpha
#endif

/**
Conversely there is a unique function computing $c$ as a function of
$\mathbf{n}$ and $\alpha$. We call this function `line_area()` and
define it as: */

#if dimension == 1
double line_area (double nx, double alpha)
{
  double area;

  alpha += nx/2.;
  if (nx < 0.) {
    alpha -= nx;
    nx = -nx;
  }

  if (alpha <= 0.)
    return 0.;

  if (alpha >= nx)
    return 1.;

  area = alpha/nx;
  return clamp (area, 0., 1.);
}
#define plane_volume(n, alpha) line_area(n.x, alpha)
#endif // dimension == 1

#if dimension >= 2
double line_area (double nx, double ny, double alpha)
{
  double a, v, area;

  alpha += (nx + ny)/2.;
  if (nx < 0.) {
    alpha -= nx;
    nx = - nx;
  }
  if (ny < 0.) {
    alpha -= ny;
    ny = - ny;
  }

  if (alpha <= 0.)
    return 0.;

  if (alpha >= nx + ny)
    return 1.;

  if (nx < 1e-10)
    area = alpha/ny;
  else if (ny < 1e-10)
    area = alpha/nx;
  else {
    v = sq(alpha);

    a = alpha - nx;
    if (a > 0.)
      v -= a*a;
    
    a = alpha - ny;
    if (a > 0.)
      v -= a*a;

    area = v/(2.*nx*ny);
  }

  return clamp (area, 0., 1.);
}
#endif // dimension >= 2

#if dimension >= 3
double plane_volume (coord n, double alpha)
{
  double al = alpha + (n.x + n.y + n.z)/2. +
    max(0., -n.x) + max(0., -n.y) + max(0., -n.z);
  if (al <= 0.)
    return 0.;
  double tmp = fabs(n.x) + fabs(n.y) + fabs(n.z);
  if (al >= tmp)
    return 1.;
  if (tmp < 1e-10)
    return 0.;
  double n1 = fabs(n.x)/tmp;
  double n2 = fabs(n.y)/tmp;
  double n3 = fabs(n.z)/tmp;
  al = max(0., min(1., al/tmp));
  double al0 = min(al, 1. - al);
  double b1 = min(n1, n2);
  double b3 = max(n1, n2);
  double b2 = n3;
  if (b2 < b1) {
    tmp = b1;
    b1 = b2;
    b2 = tmp;
  }
  else if (b2 > b3) {
    tmp = b3;
    b3 = b2;
    b2 = tmp;
  }
  double b12 = b1 + b2;
  double bm = min(b12, b3);
  double pr = max(6.*b1*b2*b3, 1e-50);
  if (al0 < b1)
    tmp = al0*al0*al0/pr;
  else if (al0 < b2)
    tmp = 0.5*al0*(al0 - b1)/(b2*b3) +  b1*b1*b1/pr;
  else if (al0 < bm)
    tmp = (al0*al0*(3.*b12 - al0) + b1*b1*(b1 - 3.*al0) +
	   b2*b2*(b2 - 3.*al0))/pr;
  else if (b12 < b3)
    tmp = (al0 - 0.5*bm)/b3;
  else
    tmp = (al0*al0*(3. - 2.*al0) + b1*b1*(b1 - 3.*al0) + 
	   b2*b2*(b2 - 3.*al0) + b3*b3*(b3 - 3.*al0))/pr;

  double volume = al <= 0.5 ? tmp : 1. - tmp;
  return clamp (volume, 0., 1.);
}
#elif dimension == 2
# define plane_volume(n, alpha) line_area(n.x, n.y, alpha)
#endif

/**
VOF algorithms require the computation of volume fractions on
(rectangular) parts of the initial square cell.

We first define a function which takes an interface definition
($\mathbf{n}$, $\alpha$), the coordinates of the lower-left `a`
and upper-right `b` corners of a rectangle and returns the
fraction of this rectangle which lies inside the interface. */

double rectangle_fraction (coord n, double alpha, coord a, coord b)
{
  coord n1;
  foreach_dimension() {
    alpha -= n.x*(b.x + a.x)/2.;
    n1.x = n.x*(b.x - a.x);
  }
  return plane_volume (n1, alpha);
}

/**
From the interface definition, it is also possible to compute the
coordinates of the segment in 2D, or facet in 3D, representing the
interface in the unit cell.

In two dimensions, the function below returns the 0,1 or 2 coordinates
(stored in the `p` array provided by the user) of the corresponding
interface segments. The case where only 1 coordinate is returned
corresponds to the degenerate case where the interface intersects the
cell exactly on a vertex. 

In three dimensions, the function returns up to 12 coordinates of the
planar fragment. */

#if dimension <= 2
int facets (coord n, double alpha, coord p[2])
{
  int i = 0;
  for (double s = -0.5; s <= 0.5; s += 1.)
    foreach_dimension()
      if (fabs (n.y) > 1e-4 && i < 2) {
	double a = (alpha - s*n.x)/n.y;
	if (a >= -0.5 && a <= 0.5) {
	  p[i].x   = s;
	  p[i++].y = a;
	}
      }
  return i;
}
#else
static coord cube_edge[12][2] = {
  {{0.,0.,0.},{1.,0.,0.}},{{0.,0.,1.},{1.,0.,1.}},
  {{0.,1.,1.},{1.,1.,1.}},{{0.,1.,0.},{1.,1.,0.}},
  {{0.,0.,0.},{0.,1.,0.}},{{0.,0.,1.},{0.,1.,1.}},
  {{1.,0.,1.},{1.,1.,1.}},{{1.,0.,0.},{1.,1.,0.}},
  {{0.,0.,0.},{0.,0.,1.}},{{1.,0.,0.},{1.,0.,1.}},
  {{1.,1.,0.},{1.,1.,1.}},{{0.,1.,0.},{0.,1.,1.}}
};

/* first index is the edge number, second index is the edge
   orientation (0 or 1), third index are the edges which this edge may
   connect to in order */
static int cube_connect[12][2][4] = {
  {{9, 1, 8}, {4, 3, 7}},   /* 0 */
  {{6, 2, 5},  {8, 0, 9}}, /* 1 */
  {{10, 3, 11},  {5, 1, 6}},  /* 2 */
  {{7, 0, 4},   {11, 2, 10}},  /* 3 */
  {{3, 7, 0},   {8, 5, 11}},  /* 4 */
  {{11, 4, 8},  {1, 6, 2}},  /* 5 */
  {{2, 5, 1},  {9, 7, 10}}, /* 6 */
  {{10, 6, 9}, {0, 4, 3}},   /* 7 */
  {{5, 11, 4},  {0, 9, 1}}, /* 8 */
  {{1, 8, 0}, {7, 10, 6}}, /* 9 */
  {{6, 9, 7},  {3, 11, 2}},   /* 10 */
  {{2, 10, 3},   {4, 8, 5}}    /* 11 */
};

int facets (coord n, double alpha, coord v[12], double h)
{
  coord a[12];
  int orient[12];

  for (int i = 0; i < 12; i++) {
    coord e, d;
    double den = 0., t = alpha;
    foreach_dimension() {
      d.x = h*(cube_edge[i][0].x - 0.5);
      e.x = h*(cube_edge[i][1].x - 0.5);
      den += n.x*(e.x - d.x);
      t -= n.x*d.x;
    }
    orient[i] = -1;
    if (fabs (den) > 1e-10) {
      t /= den;
      if (t >= 0. && t < 1.) {
	double s = - alpha;
	foreach_dimension() {
	  a[i].x = d.x + t*(e.x - d.x);
	  s += n.x*e.x;
	}
	orient[i] = (s > 0.);
      }
    }
  }

  for (int i = 0; i < 12; i++) {
    int nv = 0, e = i;
    while (orient[e] >= 0) {
      int m = 0, * ne = cube_connect[e][orient[e]];
      v[nv++] = a[e];
      orient[e] = -1;
      while (m < 3 && orient[e] < 0)
	e = ne[m++];
    }
    if (nv > 2)
      return nv;
  }
  return 0;
}
#endif // dimension == 3

/**
This function fills the coordinates *p* of the centroid of the
interface fragment and returns the length/area of the fragment. */

double line_length_center (coord m, double alpha, coord * p)
{
  alpha += (m.x + m.y)/2.;

  coord n = m;
  foreach_dimension(2)
    if (n.x < 0.) {
      alpha -= n.x;
      n.x = - n.x;
    }

  p->x = p->y = p->z = 0.;

  if (alpha <= 0. || alpha >= n.x + n.y)
    return 0.;
  
  foreach_dimension(2)
    if (n.x < 1e-4) {
      p->x = 0.;
      p->y = (m.y < 0. ? 1. - alpha : alpha) - 0.5;
      return 1.;
    }
  
  if (alpha >= n.x) {
    p->x += 1.;
    p->y += (alpha - n.x)/n.y;
  }
  else
    p->x += alpha/n.x;

  double ax = p->x, ay = p->y;
  if (alpha >= n.y) {
    p->y += 1.;
    ay -= 1.;
    p->x += (alpha - n.y)/n.x;
    ax -= (alpha - n.y)/n.x;
  }
  else {
    p->y += alpha/n.y;
    ay -= alpha/n.y;
  }

  foreach_dimension(2) {
    p->x /= 2.;
    p->x = clamp (p->x, 0., 1.);
    if (m.x < 0.)
      p->x = 1. - p->x;
    p->x -= 0.5;
  }

  return sqrt (ax*ax + ay*ay);
}

#if dimension == 2
#define plane_area_center(m,a,p) line_length_center(m,a,p)
#else // dimension == 3
double plane_area_center (coord m, double alpha, coord * p)
{
  foreach_dimension()
    if (fabs (m.x) < 1e-4) {
      coord n, q;
      ((double *)&n)[0] = m.y;
      ((double *)&n)[1] = m.z;
      double length = line_length_center (n, alpha, &q);
      p->x = 0.;
      p->y = ((double *)&q)[0];
      p->z = ((double *)&q)[1];
      return length;
    }

  alpha += (m.x + m.y + m.z)/2.;
  coord n = m;
  foreach_dimension()
    if (n.x < 0.) {
      alpha -= n.x;
      n.x = - n.x;
    }

  double amax = n.x + n.y + n.z;
  if (alpha < 0. || alpha > amax) {
    p->x = p->y = p->z = 0.;
    return 0.;
  }

  double area = sq(alpha);
  p->x = p->y = p->z = area*alpha;

  foreach_dimension() {
    double b = alpha - n.x;
    if (b > 0.) {
      area -= b*b;
      p->x -= b*b*(2.*n.x + alpha);
      p->y -= b*b*b;
      p->z -= b*b*b;
    }
  }

  amax = alpha - amax;
  foreach_dimension() {
    double b = amax + n.x;
    if (b > 0.) {
      area += b*b;
      p->y += b*b*(2.*n.y + alpha - n.z);
      p->z += b*b*(2.*n.z + alpha - n.y);
      p->x += b*b*b;
    }
  }

  area *= 3.;
  foreach_dimension() {
    if (area) {
      p->x /= area*n.x;
      p->x = clamp (p->x, 0., 1.);
    }
    else
      p->x = 0.;
    if (m.x < 0.) p->x = 1. - p->x;
    p->x -= 0.5;
  }

  return area*sqrt (1./(sq(n.x)*sq(n.y)) +
		    1./(sq(n.x)*sq(n.z)) +
		    1./(sq(n.z)*sq(n.y)))/6.;
}
#endif // dimension == 3

/**
This function fills the coordinates *p* of the centroid of the
fraction *a* of a square cell lying under the line $(m,alpha)$. */

void line_center (coord m, double alpha, double a, coord * p)
{
  alpha += (m.x + m.y)/2.;
  
  coord n = m;
  foreach_dimension(2)
    if (n.x < 0.) {
      alpha -= n.x;
      n.x = - n.x;
    }

  p->z = 0.;
  if (alpha <= 0.) {
    p->x = p->y = -0.5;
    return;
  }

  if (alpha >= n.x + n.y) {
    p->x = p->y = 0.;
    return;
  }

  foreach_dimension(2)
    if (n.x < 1e-4) {
      p->x = 0.;
      p->y = sign(m.y)*(a/2. - 0.5);
      return;
    }

  p->x = p->y = cube(alpha);
  
  foreach_dimension(2) {
    double b = alpha - n.x;
    if (b > 0.) {
      p->x -= sq(b)*(alpha + 2.*n.x);
      p->y -= cube(b);
    }
  }

  foreach_dimension(2) {
    p->x /= 6.*sq(n.x)*n.y*a;
    p->x = sign(m.x)*(p->x - 0.5);
  }
}

/**
This function fills the coordinates *p* of the centroid of the
fraction *a* of a cubic cell lying under the plane $(m,alpha)$. */

#if dimension == 2
#define plane_center(m,alpha,a,p) line_center(m,alpha,a,p)
#else // dimension == 3
void plane_center (coord m, double alpha, double a, coord * p)
{
  foreach_dimension()
    if (fabs (m.x) < 1e-4) {
      coord n, q;
      ((double *)&n)[0] = m.y;
      ((double *)&n)[1] = m.z;
      line_center (n, alpha, a, &q);
      p->x = 0.;
      p->y = ((double *)&q)[0];
      p->z = ((double *)&q)[1];
      return;
    }

  alpha += (m.x + m.y + m.z)/2.;
  coord n = m;
  foreach_dimension()
    if (n.x < 0.) {
      alpha -= n.x;
      n.x = - n.x;
    }

  if (alpha <= 0. || a == 0.) {
    p->x = p->y = p->z = -0.5;
    return;
  }

  if (alpha >= n.x + n.y + n.z || a == 1.) {
    p->x = p->y = p->z = 0.;
    return;
  }

  p->x = p->y = p->z = sq(sq(alpha));
  foreach_dimension() {
    double b = alpha - n.x;
    if (b > 0.) {
      p->x -= cube(b)*(3.*n.x + alpha);
      p->y -= sq(sq(b));
      p->z -= sq(sq(b));
    }
  }

  double amax = alpha - (n.x + n.y + n.z);
  foreach_dimension() {
    double b = amax + n.z;
    if (b > 0.) {
      p->x += cube(b)*(3.*n.x + alpha - n.y);
      p->y += cube(b)*(3.*n.y + alpha - n.x);
      p->z += sq(sq(b));
    }
  }

  double b = 24.*n.x*n.y*n.z*a;
  foreach_dimension() {
    p->x /= b*n.x;
    p->x = sign(m.x)*(p->x - 0.5);
  }
}
#endif // dimension == 3