#include "fractions.h"
#include "distance2.h"

struct cff {
  scalar c;
  face vector f; // optional
};

coord * coords_from_fractions (struct cff input) {
  coord * p;
  scalar c = input.c;
  if (!input.f.x.i)
    input.f.x.i = -1;
  //We first count how many triangles there will be:
  long int nr_tri = 0, nr_p = 0;
  foreach(reduction(+:nr_tri)) {
    if (c[] > 1e-6 && c[] < 1. - 1e-6) {
      coord n = facet_normal (point, c, input.f);
      double alpha = plane_alpha (c[], n);
      coord v[12];
      nr_tri += facets (n, alpha, v, 1.) - 2;
    }
  }
  //Each triangle has 3 coordinates and we need a `nodata` coordinate:
  p = (coord*) malloc ((nr_tri*3 + 1)*sizeof(coord));
  // Fill `p` with triplets
  foreach(){
    if (c[] > 1e-6 && c[] < 1. - 1e-6) {
      coord n = facet_normal (point, c, input.f);
      double alpha = plane_alpha (c[], n);
      coord v[12];
      int m = facets (n, alpha, v, 1.); //m is the number of coordinates
      for (int j = 0; j < m - 2; j++) {
	coord temp1 = {x + v[0].x*Delta  , y + v[0].y*Delta  , z + v[0].z*Delta}; 
	coord temp2 = {x + v[j+1].x*Delta, y + v[j+1].y*Delta, z + v[j+1].z*Delta};
	coord temp3 = {x + v[j+2].x*Delta, y + v[j+2].y*Delta, z + v[j+2].z*Delta};
	p[nr_p++] = temp1;
	p[nr_p++] = temp2;
	p[nr_p++] = temp3;
      }
    }
  }
  coord w = {nodata, nodata, nodata};
  p[nr_p] = w;
  return p;
}

struct dts {
  scalar c;           // Volume fraction field
  face vector f;      // Optional face fraction field
  scalar d;           // Optional output: centered distance field
  vertex scalar phi;  // Optional output: distance at vertices
};

void distance_to_surface (struct dts dff) {
  assert (dimension == 3); //fixme: there should be a 2D analog
  struct cff nja = (struct cff){dff.c, dff.f};
  coord * p = coords_from_fractions (nja);
  
  scalar d = automatic(dff.d);

  distance (d, p);

  if (dff.phi.i) { //Distance at vertices
    vertex scalar phi = dff.phi;
    boundary ({d});
    foreach_vertex()
      phi[] = (d[] + d[-1] + d[0,-1] + d[-1,-1] +
	       d[0,0,-1] + d[-1,0,-1] + d[0,-1,-1] + d[-1,-1,-1])/8.;
    boundary ({phi});
  }
}