colonizedtree.c

Potential flow past a digitally-generated tree,
- Initialization

Potential flow past a digitally-generated tree,

on a digitally-generated surface:

Do not mind the blue sky

#include "grid/octree.h"
#include "embed.h"
#include "poisson.h"
#include "run.h"
#include "colonization.h"
#include "../ABL/perlin.h"
#include "../tracer-particles.h"
#include "view.h"
#include "../scatter2.h"

scalar pot[], laplace[]; 
Particles P;

pot[front] = dirichlet (1);
pot[back] = dirichlet (0);
pot[embed] = neumann (0.);
vector u[];
u.n[front] = neumann(0);
u.n[back] = neumann(0);
int ja = 1000; double R1 = 3, H = 1.5, Hg = 1;

double sy = 1.4;
int main () {
  L0 = 10;
  X0 = Z0 = -L0/2;
  N = 64;
  run();
}

Initialization

In the initialization stage we first find the tree nodes via the space colonization algorithm and convert those to Branches. Next, the tree interface can be reconstructed in terms of tree-cell-volume fractions. A periodic Perlin topography is added before an equation for the flow potential (pot) is solved. This is then converted into the velocity field. Finally this field is visualized with tracer particles.

face vector fp[];
scalar cp[];
event init (t = 0) {
  //Seed atraction points in a prolate semi spheriod
  srand (0);
  coord atr[ja];
  int j = 0;
  while (j < ja) {
    double R = R1*pow(fabs(noise()), 1./3.1);
    coord vec = {noise(), noise(), noise()};
    normalize (&vec);
    if (vec.y > 0) {
      foreach_dimension()
	atr[j].x = R*vec.x;
      atr[j].y *= sy;
      atr[j].y += H + Hg;
      j++;
    }
  }

  //Get tree nodes
  int tn = 1;
  Tnode * nodes = calloc (tn, sizeof(Tnode));
  nodes[0].x = 0;
  nodes[0].y = Hg;
  nodes[0].z = 0;
  tn = colonize (atr, ja, &nodes, 1, 2, 0.2, 0.3);

 // Convert nodes to branches
  Branch * trees = calloc (499, sizeof(Branch)); //?
  printf ("# nodes: %d\n", tn);
  nb = nodestotree (nodes, tn, &trees);
  printf ("# branches: %d\n", nb);
  fflush (stdout);

  //Convert branches to volume fractions
  int it = 0;
  do {
    tree_interface (trees, cs, fs);
    boundary({cs});
    printf ("# Tree Reconstruction Iteration %d\n", ++it);
  } while (adapt_wavelet ({cs}, (double[]){0.1}, 9, 4).nf > (grid->tn/100));
  tree_interface (trees, cs, fs);
  boundary({cs, fs});
  fractions_cleanup (cs, fs);

  //Obtain surface volume and face fractions
  cp.prolongation = cs.prolongation;
  vertex scalar phi[];
  int nx = 5;
  int ny = 3;
  init_perlin (nx, ny);
  it = 0;
  do {
    foreach_vertex()
      phi[] =  y - Hg - (perlin (x, z, nx, ny));
    boundary ({phi});
    fractions (phi, cp, fp);
    boundary ({cp});
    printf ("# Surface Reconstruction Iteration %d\n", ++it);
  }  while (adapt_wavelet ({cs, cp}, (double[]){0.1, 0.01}, 8, 4).nf > (grid->tn/1000));
  foreach_vertex()
    phi[] =  y - Hg - (perlin (x, z, nx, ny));
  boundary ({phi});
  fractions (phi, cp, fp);
 
  //Combine tree and surface
  foreach()
    cs[] = min (cs[], cp[]);
  foreach_face()
    fs.x[] = min(fs.x[], fp.x[]);
  boundary ({cs, fs});
  fractions_cleanup (cs, fs);
  
  //Obain potential
  printf ("# Start solving for the Potential\n");
  NITERMIN = 10;
  mgstats mg = poisson (pot, laplace, alpha = fs,

unknown function parameter ‘embed_flux’

			embed_flux = embed_flux, tolerance = HUGE);
  printf ("# Solved with a max. residual %g\n", mg.resa);

  //Get potential flow field
  foreach() {
    foreach_dimension()
      u.x[] = center_gradient(pot) * cs[];
  }
  boundary ((scalar*){u});

  //Place tracer particles
  P = init_tp_square (n = 30, ym = H + Hg, l = 2*R1, zp = -L0/4);
  DT = 0.1;

  //free stuff!
  free (trees);
  free (nodes);
  free (gradp);
}

event set_dtmax (i++) {
  dt = dtnext(DT);
}

double th = 0;
event movie (t += 1.) {
  view (bg = {0.3, 0.3, 0.9}, theta = th, phi = 0.2);
  view (ty = -.3);
  draw_vof ("cs", "fs", fc = {0.73,0.4,0.2});
  translate (y = 0.01)
    draw_vof ("cp", "fp", fc = {0.2,0.9,0.2});
  scatter (P);
  box();
  save ("tree.mp4");
  th += 0.006;
}

event stop (t = 500);