terrain.h

Go to the documentation of this file.

/** @file terrain.h
 */
#include <stdarg.h>
#include <kdt/kdt.h>
#pragma autolink -L$BASILISK/kdt -lkdt
 
#if _OPENMP
# define NPROC omp_get_max_threads()
# define PROC tid()
#else
# define NPROC 1
# define PROC 0
#endif
 
attribute {
  void ** kdt;
  scalar nt, dmin, dmax;
}
 
static int includes_point (KdtRect rect, Point point)
{
  Delta_x /= 2.; Delta_y /= 2.;
  return (rect[0].l >= x - Delta_x && rect[0].h <= x + Delta_x &&
    rect[1].l >= y - Delta_y && rect[1].h <= y + Delta_y);
}
  
static int includes (KdtRect rect, Point * p)
{
  return includes_point (rect, *p);
}
 
static int intersects_point (KdtRect rect, Point point)
{
  Delta_x /= 2.; Delta_y /= 2.;
  return (rect[0].l <= x + Delta_x && rect[0].h >= x - Delta_x &&
    rect[1].l <= y + Delta_y && rect[1].h >= y - Delta_y);
}
 
static int intersects (KdtRect rect, Point * p)
{
  return intersects_point (rect, *p);
}
 
#if MULTIGRID
static void reconstruct_terrain (Point point, scalar zb)
{
  KdtSum s;
  kdt_sum_init (&s);
  Delta_x /= 2.; Delta_y /= 2.;
  KdtRect rect = {{x - Delta_x, x + Delta_x},
      {y - Delta_y, y + Delta_y}};
  for (Kdt ** kdt = (Kdt **) zb.kdt[PROC]; *kdt; kdt++)
    kdt_query_sum (*kdt,
       (KdtCheck) includes,
       (KdtCheck) intersects, &point,
       rect, &s);
  scalar n = zb.nt, dmin = zb.dmin, dmax = zb.dmax;
  n[] = s.n;
  if (s.w > 0.) {
    zb[] = s.H0/s.w;
    dmin[] = s.Hmin;
    dmax[] = s.Hmax;
  }
  else {
    /* not enough points in database, use bilinear interpolation
       from coarser level instead */
    if (level > 0)
      zb[] = (9.*coarse(zb,0,0) + 
        3.*(coarse(zb,child.x,0) + coarse(zb,0,child.y)) + 
        coarse(zb,child.x,child.y))/16.;
    else
      zb[] = 0.; // no points at level 0!
    dmin[] = nodata;
    dmax[] = nodata;
  }
}
 
void refine_terrain (Point point, scalar zb)
{
  for (int _c = 0; _c < 4; _c++) /* foreach_child */
    reconstruct_terrain (point, zb);
}
#endif // MULTIGRID
 
void delete_terrain (scalar zb)
{
  if (!zb.kdt)
    return;
  for (int i = 0; i < NPROC; i++) {
    for (Kdt ** kdt = (Kdt **) zb.kdt[i]; *kdt; kdt++)
      kdt_destroy (*kdt);
    free (zb.kdt[i]);
  }
  free (zb.kdt);
  zb.kdt = NULL;
}
 
@define CHARP char * // fixme: workaround for va_arg macro
 
trace
void terrain (scalar zb, ...)
{
  zb.kdt = qcalloc (NPROC, void *);
  zb.delete = delete_terrain;
  
  int nt = 0;
  va_list ap;
  va_start (ap, zb);
  char * name;
  while ((name = va_arg (ap, CHARP))) {
    for (int i = 0; i < NPROC; i++) {
      Kdt ** kdt = (Kdt **) zb.kdt[i];
      zb.kdt[i] = qrealloc (kdt, nt + 2, Kdt *);
      kdt[nt] = kdt_new();
      kdt[nt + 1] = NULL;
      char * fname = name;
      if (name[0] == '~') {
  char * home = getenv ("HOME");
  if (home != NULL) {
    fname = malloc(sizeof(char)*(strlen(home) + strlen(name)));
    strcpy (fname, home);
    strcat (fname, name + 1);
  }
      }
      if (kdt_open (kdt[nt], fname)) {  
  fprintf (stderr, "terrain: could not open terrain database '%s'\n", 
     fname);
  exit (1);
      }
      if (fname != name)
  free (fname);
    }
    nt++;
  }
  va_end (ap);
 
  scalar n = {0} /* new scalar */;
  scalar dmin = {0} /* new scalar */;
  scalar dmax = {0} /* new scalar */;
  zb.nt = n;
  zb.dmin = dmin;
  zb.dmax = dmax;
 
#if TREE
  zb.refine = refine_terrain;
  n.refine = no_restriction;
  n.restriction = no_restriction;
  dmin.refine = no_restriction;
  dmin.restriction = no_restriction;
  dmax.refine = no_restriction;
  dmax.restriction = no_restriction;
#endif
 
  /* trash: zb */;
#if MULTIGRID && !_GPU
  for (int l = 0; l <= depth(); l++) {
    for (int _l = 0; _l < l; _l++)
      reconstruct_terrain (point, zb);
    boundary_level ({zb}, l);
  }
#else
  foreach (cpu) {
    KdtSum s;
    int niter = 8;
    do {
      kdt_sum_init (&s);
      KdtRect rect = {{x - Delta_x/2., x + Delta_x/2.},
          {y - Delta_y/2., y + Delta_y/2.}};
      for (Kdt ** kdt = (Kdt **) zb.kdt[PROC]; *kdt; kdt++)
  kdt_query_sum (*kdt,
           (KdtCheck) includes,
           (KdtCheck) intersects, &point,
           rect, &s);
      Delta_x *= 2., Delta_y *= 2.;
    } while (!s.w && niter--);
    scalar n = zb.nt, dmin = zb.dmin, dmax = zb.dmax;
    n[] = s.n;
    if (s.w > 0.) {
      zb[] = s.H0/s.w;
      dmin[] = s.Hmin;
      dmax[] = s.Hmax;
    }
    else {
      zb[] = 0.;
      dmin[] = nodata;
      dmax[] = nodata;
    }
  }
#endif
#if !TREE
  delete_terrain (zb);
#endif
}