sandbox/yonghui/smalltest/ImportMesh

    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
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    
    /**
    It's a code to import mesh with no other utilities. So it can be fast to pre-check if we have set the origin & size & ect correctly. * /
    
      /**
    ## Code to import mesh */
    
      #include "grid/octree.h"
      #include "navier-stokes/centered.h"
      #include "utils.h"
      #include "distance.h"
      #include "fractions.h"
      #include "view.h"
      #define MIN_LEVEL 5
      #define LEVEL 7
      #define MAX_LEVEL 9 
      void fraction_from_stl (scalar f, FILE * fp, int maxlevel)
    {
      coord * p = input_stl (fp);
      coord min, max;
      bounding_box (p, &min, &max);
      double maxl = -HUGE;
      foreach_dimension()
        if (max.x - min.x > maxl)
        maxl = max.x - min.x;
        scalar d[];
      distance (d, p);
      while (adapt_wavelet ({d}, (double[]){1e-3*L0}, MAX_LEVEL).nf);
      vertex scalar phi[];
      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.;
      fractions (phi, f); 
    }
    scalar f0[];
    
    int main()
    {
      init_grid (64);
      //1 mm = 1
      size (10);
      origin (-L0/2 ,-5., -L0/2);
      //chaneg filename
      FILE * fp = fopen ("FILENAME.stl", "r");
      fraction_from_stl (f0, fp, LEVEL);
      f0.refine = f0.prolongation = fraction_refine;
      fclose (fp); 
      clear();
      draw_vof ("f0", edges = true, lw = 0.5);
      save ("vof.png");
      dump (file = "dump_init");
    } 
    
    /**
    ## calculate equivelent radius*/
    //calculate volume
    double volume = statsf(f0).sum;
    fprintf(stderr, "total volume : %g \n", volume);
    //calculate equivalent surface
    double eq_s =volume /L0 ;
    fprintf(stderr, "equivalent surface: %g \n",eq_s);
    //calculate equivalent radius
    double eq_r = sqrt(eq_s / M_PI) ;
    fprintf(stderr, "equivalent radius :%g \n",eq_r);
    
    /**
    ## view.bv useful commande */
    box();
    squares("u.x",n={0,1,0},alpha=0.);
    draw_vof("f");