sandbox/Antoonvh/frac-dist.h

    Distance to a fraction field

    We use a slightly modified distance() function to find the distance field to a surface that is represented by volume fractions.

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

    Therefore, we need triangles; Using Alexis Berny’s method for triangularization of a fraction field.

    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;
    }

    distance to surface

    The user interface is provided via distance_to_surface()

    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= {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});
      }
    }