src/test/hf3D.c

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    
    #include "grid/octree.h"
    #include "fractions.h"
    #include "curvature.h"
    #include "utils.h"
    
    scalar f[];
    
    int main()
    {
      init_grid (16);
      origin (-0.5 - 0.125, -0.5 - 0.125, -0.5 - 0.125);
    #if TREE
      refine (level == 4 && fabs (x) < 0.375 && fabs (y) < 0.375);
      refine (level <= 5 && fabs (x) < 0.31 && fabs (y) < 0.31);
      f.refine = f.prolongation = fraction_refine;
    #endif
      fraction (f, - (0.2 - sqrt(sq(x+0.2) + sq(y+0.2) + sq(z))));
      vector h[];
      heights (f, h);
      scalar kappa[];
      curvature (f, kappa);
    #if 0
      foreach() {
        if (h.x[] != nodata)
          fprintf (stderr, "%g %g %g %g hx\n", x, y, height(h.x[]), Delta);
        if (h.y[] != nodata)
          fprintf (stderr, "%g %g %g %g hy\n", x, y, height(h.y[]), Delta);
        if (kappa[] != nodata)
          fprintf (stderr, "%g %g %g kappa\n", x, y, kappa[]);
      }
    #endif
      stats s = statsf (kappa);
      fprintf (stderr, "kappa min: %g avg: %g stddev: %g max: %g\n",
    	   s.min, s.sum/s.volume, s.stddev, s.max);
    #if 0
      FILE * fp = popen ("gfsview3D -s hf.gfv", "w");
      output_gfs (fp);
    #endif
    }