src/draw.h

    Drawing functions for Basilisk View

    #include "fractions.h"
    #include "gl/font.h"

    clear(): removes all objects previously drawn

    void clear()
    {
      bview * view = get_view();
      if (view->active)
        view->active = false;
      draw();
    }

    view(): sets up viewing parameters

    • tx, ty: shifts the camera center point.
    • fov: changes the field-of-view.
    • quat[]: the quaternion defining the camera angles.
    • sx, sy, sz: stretch factors for each axis.
    • width, height, samples: the width and height (in pixels) of the image to render (default is 800 x 800). The image can optionally be generated by first rendering an image with samples times more pixels in each direction followed by subsampling. This provides a form of antialiasing. Default is four samples.
    • bg[]: an array of red, green, blue values between 0 and 1 which defines the background color.
    • theta, phi, psi: Euler-like angles (in radians), used (instead of quat[]) to define the camera angle.
    • relative: whether the theta and phi angles are absolute or relative to the current position (i.e. increments of the current angles).
    • camera: predefined camera angles: “left”, “right”, “top”, “bottom”, “front”, “back” and “iso”.
    • map: an optional coordinate mapping function.
    struct _view_set {
      float tx, ty;
      float fov;
      float quat[4];
      float sx, sy, sz;
      unsigned width, height, samples;
      float bg[3];
      float theta, phi, psi;
      bool relative;
      float res;
      char * camera;
      void (* map) (coord *);
      float p1x, p1y, p2x, p2y; // for trackball
      bview * view;
    };
    
    void view (struct _view_set p)
    {
      bview * v = p.view ? p.view : get_view();
      if (p.fov) {
        if (p.relative)
          v->fov += (0.1 + 3.*v->fov)*p.fov;
        else
          v->fov = p.fov;
        v->fov = clamp(v->fov,0.01,100.);
      }
      for (int i = 0; i < 4; i++)
        if (p.quat[i]) {
          for (int j = 0; j < 4; j++)
    	v->quat[j] = p.quat[j];
          break;
        }
      if (p.tx) v->tx = p.relative ? v->tx + p.tx*0.02*(0.01 + 3.*v->fov) : p.tx;
      if (p.ty) v->ty = p.relative ? v->ty + p.ty*0.02*(0.01 + 3.*v->fov) : p.ty;
      if (p.sx) v->sx = p.sx;
      if (p.sy) v->sy = p.sy;
      if (p.sz) v->sz = p.sz;
      if (p.bg[0] || p.bg[1] || p.bg[2])
        for (int i = 0; i < 3; i++)
          v->bg[i] = p.bg[i];
      
      if (p.camera) {
        v->gfsview = false;
        if (strlen(p.camera) >= 4 &&
    	!strcmp (&p.camera[strlen(p.camera) - 4], ".gfv")) {
          FILE * fp = fopen (p.camera, "r");
          if (!fp) {
    	perror (p.camera);
    	exit (1);
          }
          char s[81];
          float q[4], fov;
          int nq = 0, nf = 0;
          while (fgets (s, 81, fp) && (!nq || !nf)) {
    	if (!nq)
    	  nq = sscanf (s, "  q0 = %f q1 = %f q2 = %f q3 = %f",
    		       &q[0], &q[1], &q[2], &q[3]);
    	if (!nf)
    	  nf = sscanf (s, "  fov = %f", &fov);
          }
          if (nq != 4 || nf != 1) {
    	fprintf (stderr, "%s: not a valid gfv file\n", p.camera);
    	exit (1);
          }
          for (int j = 0; j < 4; j++)
    	v->quat[j] = q[j];
          v->fov = fov;
          v->gfsview = true;
        }
        else if (!strcmp (p.camera, "left"))
          gl_axis_to_quat ((float[]){0,1,0}, - pi/2., v->quat);
        else if (!strcmp (p.camera, "right"))
          gl_axis_to_quat ((float[]){0,1,0}, pi/2., v->quat);
        else if (!strcmp (p.camera, "top"))
          gl_axis_to_quat ((float[]){1,0,0}, - pi/2., v->quat);
        else if (!strcmp (p.camera, "bottom"))
          gl_axis_to_quat ((float[]){1,0,0}, pi/2., v->quat);
        else if (!strcmp (p.camera, "front"))
          gl_axis_to_quat ((float[]){0,0,1}, 0., v->quat);
        else if (!strcmp (p.camera, "back"))
          gl_axis_to_quat ((float[]){0,1,0}, pi, v->quat);
        else if (!strcmp (p.camera, "iso")) {
          gl_axis_to_quat ((float[]){0,1,0}, pi/4., v->quat);
          float q[4];
          gl_axis_to_quat ((float[]){1,0,0}, - pi/4., q);
          gl_add_quats(q, v->quat, v->quat);
        }
        else {
          fprintf (stderr, "view(): unknown camera '%s'\n", p.camera);
          exit (1);
        }
      }
      else if (p.theta || p.phi || p.psi) {
        v->gfsview = false;
        float q[4];
        gl_axis_to_quat ((float[]){1,0,0}, - p.phi, q);
        if (p.relative) {
          float q1[4];
          gl_axis_to_quat ((float[]){0,1,0}, p.theta, q1);
          gl_add_quats(q, q1, q1);
          float q2[4];
          gl_axis_to_quat ((float[]){0,0,1}, p.psi, q2);
          gl_add_quats(q1, q2, q2);
          gl_add_quats(q2, v->quat, v->quat);
        }
        else {
          gl_axis_to_quat ((float[]){0,1,0}, p.theta, v->quat);
          gl_add_quats(q, v->quat, v->quat);
          gl_axis_to_quat ((float[]){0,0,1}, p.psi, q);
          gl_add_quats(q, v->quat, v->quat);
        }
      }
    
      if (p.map)
        v->map = p.map;
      
      if (p.p1x || p.p1y || p.p2x || p.p2y) { // trackball
        float q[4];
        gl_trackball(q, p.p1x, p.p1y, p.p2x, p.p2y);
        gl_add_quats (q, v->quat, v->quat);
      }
    
      if (p.res)
        v->res = p.res;
      
      if ((p.width && p.width != v->width) ||
          (p.height && p.height != v->height) ||
          (p.samples && p.samples != v->samples)) {
        v->width = v->width/v->samples;
        v->height = v->height/v->samples;
        if (p.width) v->width = p.width;
        if (p.height) v->height = p.height;
        if (p.samples) v->samples = p.samples;
        v->width *= v->samples;
        v->height *= v->samples;
        framebuffer_destroy (v->fb);
        v->fb = framebuffer_new (v->width, v->height);
        init_gl();
      }
    
      clear();
    }

    translate(): translates the origin.

    The block following this command will be drawn in a translated coordinate system.

    struct _translate {
      float x, y, z;
    };
    
    void begin_translate (struct _translate p)
    {
      bview * view = draw();
      glMatrixMode (GL_MODELVIEW);
      glPushMatrix();
      glTranslatef (p.x, p.y, p.z);
      gl_get_frustum (&view->frustum);
    }
    
    void end_translate()
    {
      bview * view = draw();
      glMatrixMode (GL_MODELVIEW);
      glPopMatrix();
      gl_get_frustum (&view->frustum);
    }

    mirror(): symmetry relative to a plane.

    The block following this command will be drawn in a coordinate system symmetric relative to the given plane. The plane is given by n and \alpha as explained in squares().

    struct _mirror {
      coord n;
      double alpha;
    };
    
    void begin_mirror (struct _mirror p)
    {
      bview * view = draw();
      glMatrixMode (GL_MODELVIEW);
      glPushMatrix();
      normalize (&p.n);
      GLfloat s[16], t[16];
      s[0] = 1. - 2.*p.n.x*p.n.x;
      s[1] = - 2.*p.n.x*p.n.y;  s[2] = - 2.*p.n.x*p.n.z;
      s[3] = 0.;
      s[4] = s[1];
      s[5] = 1. - 2.*p.n.y*p.n.y; s[6] = - 2.*p.n.y*p.n.z;
      s[7] = 0.;
      s[8] = s[2];   s[9] = s[6];  s[10] = 1. - 2.*p.n.z*p.n.z; 
      s[11] = 0.;
      s[12] = 0.;    s[13] = 0.;   s[14] = 0.;                    
      s[15] = 1.;
    
      t[0] = 1.;  t[1] = 0.;   t[2] = 0.;  t[3] = 0.;
      t[4] = 0.;  t[5] = 1.;   t[6] = 0.;  t[7] = 0.;
      t[8] = 0.;  t[9] = 0.;   t[10] = 1.; t[11] = 0.;
      t[12] = - 2.*p.n.x*p.alpha; 
      t[13] = - 2.*p.n.y*p.alpha;  
      t[14] = - 2.*p.n.z*p.alpha; 
      t[15] = 1.;
      matrix_multiply (s, t);
      glMultMatrixf (s);
      gl_get_frustum (&view->frustum);
      view->reversed = !view->reversed;
    }
    
    void end_mirror() {
      end_translate();
      bview * view = draw();
      view->reversed = !view->reversed;
    }

    Utility functions

    The tree structure is used to traverse only the cells which are within the field of view of the camera.

    static void mapped_position (bview * view, coord * p, double * r)
    {
      double x = p->x, y = p->y, z = p->z, rm = 0.;
      view->map (p);
      for (int i = -1; i <= 1; i += 2)
        for (int j = -1; j <= 1; j += 2)
          for (int k = -1; k <= 1; k += 2) {
    	coord q = {x + i**r, y + j**r, z + k**r};
    	view->map (&q);
    	double pq = sq(p->x - q.x) + sq(p->y - q.y) + sq(p->z - q.z);
    	if (pq > rm)
    	  rm = pq;
          }
      *r = sqrt (rm);
    }
    
    @def foreach_visible(view)
    foreach_cell() {
    #if dimension == 2
      double _r = Delta*0.71;
    #else // dimension == 3
      double _r = Delta*0.87;
    #endif
      coord _p = {x, y, z};
      if ((view)->map)
        mapped_position (view, &_p, &_r);
      if (!sphere_in_frustum (_p.x, _p.y, _p.z, _r, &(view)->frustum))
        continue;
      if (is_leaf(cell) ||
          sphere_diameter (_p.x, _p.y, _p.z, _r/L0, &(view)->frustum)
          < (view)->res) {
        if (is_active(cell) && is_local(cell)) {
    @
    @def end_foreach_visible()
        }
        continue;
      }
    }
    end_foreach_cell();
    @

    A similar technique can be used to traverse the cells which are both visible and intersected by a plane defined by \displaystyle n_x x + n_y y + n_z z = \alpha

    #if dimension == 3
    static void glnormal3d (bview * view, double x, double y, double z) {
      // fixme: mapping? (see glvertex3d)
      if (view->gfsview || view->reversed)
        glNormal3d (- x, - y, - z);
      else
        glNormal3d (x, y, z);
    }
    
    @def foreach_visible_plane(view, n1, alpha1)
    coord n = {(n1).x, (n1).y, (n1).z};
    double _alpha = 0.9999999*(alpha1);
    {
      double norm = sqrt(sq(n.x) + sq(n.y) + sq(n.z));
      if (!norm)
        n.z = 1.;
      else
        n.x /= norm, n.y /= norm, n.z /= norm, _alpha /= norm;
    }
    glnormal3d (view, n.x, n.y, n.z); // do not use normal inversion
    foreach_cell() {
      // fixme: coordinate mapping
      double _r = Delta*0.87, alpha = (_alpha - n.x*x - n.y*y - n.z*z)/Delta;
      if (fabs(alpha) > 0.87 || !sphere_in_frustum (x, y, z, _r, &(view)->frustum))
        continue;
      if (is_leaf(cell) ||
          sphere_diameter (x, y, z, _r/L0, &(view)->frustum) < (view)->res) {
        if (is_active(cell) && is_local(cell)) {
    @
    @def end_foreach_visible_plane()
        }
        continue;
      }
    }
    end_foreach_cell();
    @
    #endif // dimension == 3
    
    static scalar lookup_field (const char * name)
    {
      if (name)
        for (scalar s in all)
          if (!strcmp (s.name, name))
    	return s;
      return (scalar){-1};
    }
    
    static vector lookup_vector (const char * name)
    {
      if (name) {
        char component[strlen(name) + 3];
        strcpy (component, name);
        strcat (component, ".x");
        for (scalar s in all)
          if (!strcmp (s.name, component))
    	return s.v;
      }
      return (vector){{-1}};
    }
    
    static bool _reversed = false;
    
    static void begin_draw_lines (bview * view, float color[3], float lw)
    {
      glMatrixMode (GL_PROJECTION);
      glPushMatrix();
      glTranslatef (0., 0., view->lc*view->fov/24.);
      glColor3f (color[0], color[1], color[2]);
      glLineWidth (view->samples*(lw > 0. ? lw : 1.));
      _reversed = view->reversed;
      view->reversed = false;
    }
    
    static void end_draw_lines()
    {    
      glMatrixMode (GL_PROJECTION);
      glPopMatrix();
      bview * view = draw();
      view->reversed = _reversed;
    }
    
    static inline double interp (Point point, coord p, scalar col) {
      struct _interpolate _r = { col, x + p.x*Delta, y + p.y*Delta, z + p.z*Delta };
      return interpolate_linear (point, _r);
    }
    
    #define colorize_args(args)						\
      scalar col = {-1};							\
      if (args.color && strcmp (args.color, "level")) {			\
        col = lookup_field (args.color);					\
        if (col.i < 0) {							\
          fprintf (stderr, "colorize_args(): no field named '%s'\n", args.color); \
          return false;							\
        }									\
      }									\
    									\
      double cmap[NCMAP][3];						\
      if (args.color) {							\
        if (args.min == 0 && args.max == 0) {				\
          if (col.i < 0) /* level */					\
    	args.min = 0, args.max = depth();				\
          else {								\
    	stats s = statsf (col);						\
    	double avg = s.sum/s.volume;					\
    	if (args.spread < 0.)						\
    	  args.min = s.min, args.max = s.max;				\
    	else {								\
    	  double spread = (args.spread ? args.spread : 5.)*s.stddev;	\
    	  args.min = avg - spread; args.max = avg + spread;		\
    	}								\
          }									\
        }									\
        if (!args.map)							\
          args.map = jet;							\
        args.map (cmap);							\
      }									\
      									\
      if ((dimension > 2 || args.linear) &&					\
          !args.fc[0] && !args.fc[1] && !args.fc[2])			\
        args.fc[0] = args.fc[1] = args.fc[2] = 1.;
    
    #define color_facet(args)						\
      if (args.color && (!args.linear || col.i < 0)) {			\
        color b = colormap_color (cmap, col.i < 0 ?				\
    			      (double) level : val(col,0,0,0),		\
    			      args.min, args.max);			\
        glColor3f (b.r/255., b.g/255., b.b/255.);				\
      }
    
    #define color_vertex(args, val)						\
      if (args.color && args.linear && col.i >= 0) {			\
        if (view->vector) {							\
          color b = colormap_color (cmap, val, args.min, args.max);		\
          glColor3f (b.r/255., b.g/255., b.b/255.);				\
        }									\
        else {								\
          double _v = val;							\
          glTexCoord1d (clamp(((_v) - args.min)/(args.max - args.min), 0., 1.)); \
        }									\
      }
    
    static void begin_colorized (float fc[3],
    			     double cmap[NCMAP][3], bool use_texture)
    {
      // do not use textures for vector graphics (not supported by GL2PS)
      if (use_texture) {
        GLfloat texture[3*256];
        for (int i = 0; i < 256; i++) {
          color c = colormap_color (cmap, i/255., 0, 1);
          texture[3*i] = c.r/255.;
          texture[3*i + 1] = c.g/255.;
          texture[3*i + 2] = c.b/255.;
        }
        glTexImage1D (GL_TEXTURE_1D, 0, GL_RGB, 256,0, GL_RGB, GL_FLOAT, texture);
        glTexParameteri (GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
        glTexParameteri (GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
        glTexParameteri (GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
        glTexParameteri (GL_TEXTURE_1D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
        glEnable (GL_TEXTURE_1D);
      }
      glColor3f (fc[0], fc[1], fc[2]);
    }
    
    static void end_colorized() {
      glDisable (GL_TEXTURE_1D);
    }
    
    #define colorize() colorized (p.fc, cmap, !view->vector &&		\
    			      p.color && p.linear && col.i >= 0)

    draw_vof(): displays VOF-reconstructed interfaces

    • c: the name (as a string) of the Volume-Of-Fluid field.
    • s: the (optional) name of the face fraction field.
    • edges: whether to display the edges of the facets.
    • larger: makes each cell larger by this factor. This helps close the gaps in the VOF interface representation. Default is 1.1 in 3D and when edges are not displayed, otherwise it is 1.
    • filled: in 2D, whether to fill the inside (1) or outside (-1).
    • color: use this field to color each interface fragment.
    • min, max: the minimum and maximum values to use for color mapping.
    • spread: the “spread factor” to use if min and max are not defined. The maximum and minimum values will be taken as the average plus or minus spread times the standard deviation. Default is 5. If negative, the minimum and maximum values of the field are used.
    • linear: if true the color will be linearly interpolated for each vertex of the facet.
    • map: the colormap to use. Default is jet.
    • fc[]: an array of red, green, blue values between 0 and 1 which defines the facet color.
    • lc[]: an array of red, green, blue values between 0 and 1 which defines the line color.
    • lw: the line width.
    struct _draw_vof {
      char * c;
      char * s;
      bool edges;
      double larger;
      int filled;
    
      char * color;
      double min, max, spread;
      bool linear;
      colormap map;
      float fc[3], lc[3], lw;
    };

    The somewhat complicated function below checks whether an interface fragment is present within a given cell. The interface is defined by the volume fraction field c. cmin is the threshold below which a fragment is considered too small.

    static bool cfilter (Point point, scalar c, double cmin)
    {
      double cmin1 = 4.*cmin;
      if (c[] <= cmin) {
        foreach_dimension()
          if (c[1] >= 1. - cmin1 || c[-1] >= 1. - cmin1)
    	return true;
        return false;
      }
      if (c[] >= 1. - cmin) {
        foreach_dimension()
          if (c[1] <= cmin1 || c[-1] <= cmin1)
    	return true;
        return false;
      }
      int n = 0;
      double min = HUGE, max = - HUGE;
      foreach_neighbor(1) {
        if (c[] > cmin && c[] < 1. - cmin && ++n >= (1 << dimension))
          return true;
        if (c[] > max) max = c[];
        if (c[] < min) min = c[];
      }
      return max - min > 0.5;
    }
    
    #if dimension <= 2
    static void glvertex2d (bview * view, double x, double y) {
      if (view->map) {
        coord p = {x, y, 0.};
        view->map (&p);
        glVertex2d (p.x, p.y);
      }
      else
        glVertex2d (x, y);
    }
    #else // dimension > 2
    static void glvertex3d (bview * view, double x, double y, double z) {
      if (view->map) {
        coord p = {x, y, z};
        view->map (&p);
        glVertex3d (p.x, p.y, p.z);
      }
      else
        glVertex3d (x, y, z);
    }
    #endif // dimension > 2
    
    trace
    bool draw_vof (struct _draw_vof p)
    {
      scalar c = lookup_field (p.c);
      if (c.i < 0) {
        fprintf (stderr, "draw_vof(): no field named '%s'\n", p.c);
        return false;
      }
      face vector s = lookup_vector (p.s);
      
      colorize_args (p);
      
      double cmin = 1e-3; // do not reconstruct fragments smaller than this
    
    #if TREE
      // make sure we prolongate properly
      void (* prolongation) (Point, scalar) = c.prolongation;
      if (prolongation != fraction_refine) {
        c.prolongation = fraction_refine;
        boundary ({c});
      }
    #endif // TREE
        
      bview * view = draw();
    #if dimension == 2
      if (p.filled) {
        glColor3f (p.fc[0], p.fc[1], p.fc[2]);
        glNormal3d (0, 0, view->reversed ? -1 : 1);
        foreach_visible (view) {
          if ((p.filled > 0 && c[] >= 1.) || (p.filled < 0 && c[] <= 0.)) {
    	glBegin (GL_QUADS);
    	glvertex2d (view, x - Delta/2., y - Delta/2.);
    	glvertex2d (view, x + Delta/2., y - Delta/2.);
    	glvertex2d (view, x + Delta/2., y + Delta/2.);
    	glvertex2d (view, x - Delta/2., y + Delta/2.);
    	glEnd();
    	view->ni++;
          }
          else if (c[] > 0. && c[] < 1.) {
    	coord n = facet_normal (point, c, s), s = {1.,1.};
    	if (p.filled < 0)
    	  foreach_dimension()
    	    n.x = - n.x;
    	double alpha = plane_alpha (p.filled < 0. ? 1. - c[] : c[], n);
    	alpha += (n.x + n.y)/2.;
    	foreach_dimension()
    	  if (n.x < 0.) alpha -= n.x, n.x = - n.x, s.x = - 1.;
    	coord v[5];
    	int nv = 0;
    	if (alpha >= 0. && alpha <= n.x) {
    	  v[nv].x = alpha/n.x, v[nv++].y = 0.;
    	  if (alpha <= n.y)
    	    v[nv].x = 0., v[nv++].y = alpha/n.y;
    	  else if (alpha >= n.y && alpha - n.y <= n.x) {
    	    v[nv].x = (alpha - n.y)/n.x, v[nv++].y = 1.;
    	    v[nv].x = 0., v[nv++].y = 1.;
    	  }
    	  v[nv].x = 0., v[nv++].y = 0.;
    	}
    	else if (alpha >= n.x && alpha - n.x <= n.y) {
    	  v[nv].x = 1., v[nv++].y = (alpha - n.x)/n.y;
    	  if (alpha >= n.y && alpha - n.y <= n.x) {
    	    v[nv].x = (alpha - n.y)/n.x, v[nv++].y = 1.;
    	    v[nv].x = 0., v[nv++].y = 1.;
    	  }
    	  else if (alpha <= n.y)
    	    v[nv].x = 0., v[nv++].y = alpha/n.y;
    	  v[nv].x = 0., v[nv++].y = 0.;
    	  v[nv].x = 1., v[nv++].y = 0.;
    	}
    	glBegin (GL_POLYGON);
    	if (s.x*s.y < 0.)
    	  for (int i = nv - 1; i >= 0; i--)
    	    glvertex2d (view, x + s.x*(v[i].x - 0.5)*Delta,
    			y + s.y*(v[i].y - 0.5)*Delta);
    	else
    	  for (int i = 0; i < nv; i++)
    	    glvertex2d (view, x + s.x*(v[i].x - 0.5)*Delta,
    			y + s.y*(v[i].y - 0.5)*Delta);
    	glEnd ();
    	view->ni++;
          }
        }
      }
      if (!p.filled || p.edges)
        draw_lines (view, p.lc, p.lw) {
          glBegin (GL_LINES);
          foreach_visible (view)
    	if (cfilter (point, c, cmin)) {
    	  coord n = facet_normal (point, c, s);
    	  double alpha = plane_alpha (c[], n);
    	  coord segment[2];
    	  if (facets (n, alpha, segment) == 2) {
    	    glvertex2d (view, x + segment[0].x*Delta, y + segment[0].y*Delta);
    	    glvertex2d (view, x + segment[1].x*Delta, y + segment[1].y*Delta);
    	    view->ni++;
    	  }
    	}
          glEnd ();
        }
    #else // dimension == 3
      double larger =
        p.larger ? p.larger : p.edges || (p.color && !p.linear) ? 1. : 1.1;
      colorize() {
        foreach_visible (view)
          if (cfilter (point, c, cmin)) {
    	coord n = facet_normal (point, c, s);
    	double alpha = plane_alpha (c[], n);
    	coord v[12];
    	int m = facets (n, alpha, v, larger);
    	if (m > 2) {
    	  color_facet (p);
    	  glnormal3d (view, n.x, n.y, n.z);
    	  glBegin (GL_POLYGON);
    	  for (int i = 0; i < m; i++) {
    	    color_vertex (p, interp (point, v[i], col));
    	    glvertex3d (view,
    			x + v[i].x*Delta, y + v[i].y*Delta, z + v[i].z*Delta);
    	  }
    	  glEnd ();
    	  view->ni++;
    	}
          }
      }
      if (p.edges)
        draw_lines (view, p.lc, p.lw) {
          foreach_visible (view)
    	if (cfilter (point, c, cmin)) {
    	  coord n = facet_normal (point, c, s);
    	  double alpha = plane_alpha (c[], n);
    	  coord v[12];
    	  int m = facets (n, alpha, v, larger);
    	  if (m > 2) {
    	    glBegin (GL_LINE_LOOP);
    	    for (int i = 0; i < m; i++)
    	      glvertex3d (view,
    			  x + v[i].x*Delta, y + v[i].y*Delta, z + v[i].z*Delta);
    	    glEnd ();
    	    view->ni++;
    	  }
    	}
        }
    #endif // dimension == 3
    
    #if TREE
      // revert prolongation
      if (prolongation != fraction_refine) {
        c.prolongation = prolongation;
        boundary ({c});
      }
    #endif // TREE
    
      return true;
    }

    isoline(): displays isolines

    Draws a single isoline at val of field phi, or n isolines between min and max (included).

    Extra parameters are the same as for draw_vof().

    struct _isoline {
      char * phi;
      double val;
      int n;
      
      // all fields below must be identical to struct _draw_vof above
      char * c;
      char * s;
      bool edges;
      double larger;
      int filled;
    
      char * color;
      double min, max, spread;
      bool linear;
      colormap map;
      float fc[3], lc[3], lw;
    };
    
    trace
    bool isoline (struct _isoline p)
    {
    #if dimension == 2
      if (!p.color) p.color = p.phi;
      colorize_args (p);
      scalar phi = col, fiso[];
      face vector siso[];
      p.c = "fiso", p.s = "siso";
      struct _draw_vof a = *((struct _draw_vof *)&p.c);
      if (p.n < 2) {
        fractions (phi, fiso, siso, p.val);
        draw_vof (a);
      }
      else if (p.max > p.min) {
        double dv = (p.max - p.min)/(p.n - 1);
        for (p.val = p.min; p.val <= p.max; p.val += dv) {
          fractions (phi, fiso, siso, p.val);
          draw_vof (a);
        }
      }
    #else // dimension == 3
      assert (false);
    #endif // dimension == 3
      return true;
    }

    cells(): displays grid cells

    In 3D the intersections of the cells with a plane are displayed. The default plane is z=0. This can be changed by setting n and alpha which define the plane \displaystyle n_x x + n_y y + n_z z = \alpha

    struct _cells {
      coord n;
      double alpha;
      float lc[3], lw; // the line color and width
    };
    
    trace
    void cells (struct _cells p)
    {
      bview * view = draw();
      draw_lines (view, p.lc, p.lw) {
    #if dimension == 2
        foreach_visible (view) {
          glBegin (GL_LINE_LOOP);
          glvertex2d (view, x - Delta/2., y - Delta/2.);
          glvertex2d (view, x + Delta/2., y - Delta/2.);
          glvertex2d (view, x + Delta/2., y + Delta/2.);
          glvertex2d (view, x - Delta/2., y + Delta/2.);
          glEnd();
          view->ni++;
        }
    #else // dimension == 3
        foreach_visible_plane (view, p.n, p.alpha) {
          coord v[12];
          int m = facets (n, alpha, v, 1.);
          if (m > 2) {
    	glBegin (GL_LINE_LOOP);
    	for (int i = 0; i < m; i++)
    	  glvertex3d (view, x + v[i].x*Delta, y + v[i].y*Delta, z + v[i].z*Delta);
    	glEnd ();
    	view->ni++;
          }
        }
    #endif // dimension == 3
      }
    }

    vectors(): displays vector fields

    The vectors are scaled using the scale factor.

    struct _vectors {
      char * u;
      double scale;    // default is one
      float lc[3], lw; // the line color and width
    };
    
    trace
    void vectors (struct _vectors p)
    {
    #if dimension == 2
      vector u;
      struct { char x, y, z; } index = {'x', 'y', 'z'};
      foreach_dimension() {
        char name[80];
        sprintf (name, "%s.%c", p.u, index.x);
        u.x = lookup_field (name);
      }
      bview * view = draw();
      draw_lines (view, p.lc, p.lw) {
        double scale = p.scale ? p.scale : 1.;
        glBegin (GL_LINES);
        foreach_visible (view)
          if (u.x[] != nodata) {
    	glvertex2d (view, x, y);
    	glvertex2d (view, x + scale*u.x[], y + scale*u.y[]);
    	view->ni++;
          }
        glEnd();
      }
    #else // dimension == 3
      assert (false); // not implemented yet
    #endif // dimension == 3
    }

    squares(): displays colormapped fields

    The field name is given by color. The min, max, spread, map etc. arguments work as described in draw_vof().

    In 3D the intersections of the field with a plane are displayed. The default plane is z=0. This can be changed by setting n and alpha which define the plane \displaystyle n_x x + n_y y + n_z z = \alpha

    struct _squares {
      char * color;
      double min, max, spread;
      bool linear;
      colormap map;
      float fc[3], lc[3];
      
      coord n;
      double alpha;
    };
    
    trace
    bool squares (struct _squares p)
    {
      colorize_args (p);
      scalar f = col;
      
      bview * view = draw();
      glShadeModel (GL_SMOOTH);
      if (p.linear) {
        colorize() {
    #if dimension == 2
          glNormal3d (0, 0, view->reversed ? -1 : 1);
          foreach_visible (view)
            if (f[] != nodata) {
    	  glBegin (GL_TRIANGLE_FAN);
    	  color_vertex (p, (4.*f[] +
    			    2.*(f[1] + f[-1] + f[0,1] + f[0,-1]) +
    			    f[-1,-1] + f[1,1] + f[-1,1] + f[1,-1])/16.);
    	  glvertex2d (view, x, y);
    	  color_vertex (p, (f[] + f[-1] + f[-1,-1] + f[0,-1])/4.);
    	  glvertex2d (view, x - Delta/2., y - Delta/2.);
    	  color_vertex (p, (f[] + f[1] + f[1,-1] + f[0,-1])/4.);
    	  glvertex2d (view, x + Delta/2., y - Delta/2.);
    	  color_vertex (p, (f[] + f[1] + f[1,1] + f[0,1])/4.);
    	  glvertex2d (view, x + Delta/2., y + Delta/2.);
    	  color_vertex (p, (f[] + f[-1] + f[-1,1] + f[0,1])/4.);
    	  glvertex2d (view, x - Delta/2., y + Delta/2.);
    	  color_vertex (p, (f[] + f[-1] + f[-1,-1] + f[0,-1])/4.);
    	  glvertex2d (view, x - Delta/2., y - Delta/2.);
    	  glEnd();
    	  view->ni++;
    	}
    #else // dimension == 3
          foreach_visible_plane (view, p.n, p.alpha)
    	if (f[] != nodata) {
    	  coord v[12];
    	  int m = facets (n, alpha, v, 1.);
    	  if (m > 2) {
    	    coord c = {0,0,0};
    	    for (int i = 0; i < m; i++)
    	      foreach_dimension()
    		c.x += v[i].x/m;
    	    glBegin (GL_TRIANGLE_FAN);
    	    color_vertex (p, interp (point, c, f));
    	    glvertex3d (view, x + c.x*Delta, y + c.y*Delta, z + c.z*Delta);
    	    for (int i = 0; i < m; i++) {
    	      color_vertex (p, interp (point, v[i], f));
    	      glvertex3d (view,
    			  x + v[i].x*Delta, y + v[i].y*Delta, z + v[i].z*Delta);
    	    }
    	    color_vertex (p, interp (point, v[0], f));
    	    glvertex3d (view,
    			x + v[0].x*Delta, y + v[0].y*Delta, z + v[0].z*Delta);
    	    glEnd ();
    	    view->ni++;
    	  }
    	}
    #endif // dimension == 3
        }
      }
      else { // !p.linear
    #if dimension == 2
        glNormal3d (0, 0, view->reversed ? -1 : 1);
        glBegin (GL_QUADS);
        foreach_visible (view)
          if (f[] != nodata) {
    	color_facet (p);
    	glvertex2d (view, x - Delta/2., y - Delta/2.);
    	glvertex2d (view, x + Delta/2., y - Delta/2.);
    	glvertex2d (view, x + Delta/2., y + Delta/2.);
    	glvertex2d (view, x - Delta/2., y + Delta/2.);
    	view->ni++;
          }
        glEnd();
    #else // dimension == 3
        foreach_visible_plane (view, p.n, p.alpha)
          if (f[] != nodata) {
    	coord v[12];
    	int m = facets (n, alpha, v, 1.);
    	if (m > 2) {
    	  color_facet (p);
    	  glBegin (GL_POLYGON);
    	  for (int i = 0; i < m; i++)
    	    glvertex3d (view,
    			x + v[i].x*Delta, y + v[i].y*Delta, z + v[i].z*Delta);
    	  glEnd ();
    	  view->ni++;
    	}
          }
    #endif // dimension == 3
      }
      return true;
    }

    box(): displays box boundaries and axis coordinates

    • notics: do not draw tick marks (default is false).
    • lc[]: an array of red, green, blue values between 0 and 1 which defines the line color.
    • lw: the line width.
    struct _box {
      bool notics;
      float lc[3], lw;
    };
    	  
    trace
    bool box (struct _box p)
    {
      bview * view = draw();
      draw_lines (view, p.lc, p.lw) {
    
        float height = 0.5*gl_StrokeHeight();
        float width = gl_StrokeWidth ('1'), scale = L0/(60.*width), length;
        float Z1 = dimension == 2 ? 0. : Z0;
        char label[80];
      
        glMatrixMode (GL_MODELVIEW);
    
        if (!p.notics) {
          int nt = 8;
          for (int i = 0; i <= nt; i++) {
    	glPushMatrix();
    	glTranslatef (X0 + i*L0/nt - height/2.*scale, Y0 - width/3.*scale, Z1);
    	glRotatef (-90, 0, 0, 1);
    	glScalef (scale, scale, scale);
    	sprintf (label, "%g", X0 + i*L0/nt);
    	gl_StrokeString (label);
    	glPopMatrix();
    
    	glPushMatrix();
    	sprintf (label, "%g", Y0 + i*L0/nt);
    	length = gl_StrokeLength (label);
    	glTranslatef (X0 - (length + width/3.)*scale,
    		      Y0 + i*L0/nt - height/2.*scale, Z1);
    	glScalef (scale, scale, scale);
    	gl_StrokeString (label);
    	glPopMatrix();
    
    #if dimension > 2
    	glPushMatrix();
    	sprintf (label, "%g", Z0 + i*L0/nt);
    	length = gl_StrokeLength (label);
    	glTranslatef (X0 - (length + width/3.)*scale,
    		      Y0, Z0 + i*L0/nt + height/2.*scale);
    	glRotatef (-90, 1, 0, 0);
    	glScalef (scale, scale, scale);
    	gl_StrokeString (label);
    	glPopMatrix();
    #endif
          }
    
          glPushMatrix();
          sprintf (label, "%g", X0 + L0/2.);
          length = gl_StrokeLength (label);
          glTranslatef (X0 + L0/2 - height*scale, Y0 - (length + 4.*width)*scale, Z1);
          glScalef (2.*scale, 2.*scale, 2.*scale);
          gl_StrokeString ("X");
          glPopMatrix();
    
      
          glPushMatrix();
          sprintf (label, "%g", Y0 + L0/2.);
          length = gl_StrokeLength (label);
          glTranslatef (X0 - (length + 4.*width)*scale,
    		    Y0 + L0/2. - height*scale, Z1);
          glScalef (2.*scale, 2.*scale, 2.*scale);
          gl_StrokeString ("Y");
          glPopMatrix();
    
    #if dimension > 2
          glPushMatrix();
          sprintf (label, "%g", Z0 + L0/2.);
          length = gl_StrokeLength (label);
          glTranslatef (X0 - (length + 4.*width)*scale,
    		    Y0, Z0 + L0/2. + height*scale);
          glRotatef (-90, 1, 0, 0);
          glScalef (2.*scale, 2.*scale, 2.*scale);
          gl_StrokeString ("Z");
          glPopMatrix();
    #endif
        }
      
    #if dimension == 2
        foreach_level (0) {
          glBegin (GL_LINE_LOOP);
          glvertex2d (view, x - Delta/2., y - Delta/2.);
          glvertex2d (view, x + Delta/2., y - Delta/2.);
          glvertex2d (view, x + Delta/2., y + Delta/2.);
          glvertex2d (view, x - Delta/2., y + Delta/2.);
          glEnd ();
          view->ni++;
        }  
    #else // dimension != 2
        foreach_level (0) {
          for (int i = -1; i <= 1; i += 2) {
    	glBegin (GL_LINE_LOOP);
    	glvertex3d (view, x - Delta/2., y - Delta/2., z + i*Delta/2.);
    	glvertex3d (view, x + Delta/2., y - Delta/2., z + i*Delta/2.);
    	glvertex3d (view, x + Delta/2., y + Delta/2., z + i*Delta/2.);
    	glvertex3d (view, x - Delta/2., y + Delta/2., z + i*Delta/2.);
    	glEnd ();
    	view->ni++;
    	glBegin (GL_LINES);
    	for (int j = -1; j <= 1; j += 2) {
    	  glvertex3d (view, x + i*Delta/2., y + j*Delta/2., z - Delta/2.);
    	  glvertex3d (view, x + i*Delta/2., y + j*Delta/2., z + Delta/2.);
    	}
    	glEnd ();
    	view->ni++;
          }
        }
    #endif // dimension != 2
      }
      return true;
    }

    isosurface(): displays an isosurface of a field

    • f: the name (as a string) of the field.
    • v: the value of the isosurface.
    • edges: whether to draw the edges of isosurface facets.
    • color: use this field to color each interface fragment.

    The min, max, spread, map etc. arguments work as described in draw_vof().

    struct _isosurface {
      char * f;
      double v;
      bool edges;
      
      char * color;
      double min, max, spread;
      bool linear;
      colormap map;
      float fc[3], lc[3], lw;
    };
    
    trace
    bool isosurface (struct _isosurface p)
    {
    #if dimension > 2
      scalar f = lookup_field (p.f);
      if (f.i < 0) {
        fprintf (stderr, "isosurface(): no field named '%s'\n", p.f);
        return false;
      }
    
      colorize_args (p);
    
      vertex scalar v[];
      foreach_vertex()
        v[] = (f[] + f[-1] + f[0,-1] + f[-1,-1] +
    	   f[0,0,-1] + f[-1,0,-1] + f[0,-1,-1] + f[-1,-1,-1])/8.;
      boundary ({v});
      
      vector n[];
      foreach()
        foreach_dimension()
          n.x[] = center_gradient(f);
      boundary ((scalar *){n});
    
      bview * view = draw();
      glShadeModel (GL_SMOOTH);
      colorize() {
        foreach_visible (view) {
          double val[8] = {
    	v[0,0,0], v[1,0,0], v[1,0,1], v[0,0,1],
    	v[0,1,0], v[1,1,0], v[1,1,1], v[0,1,1]
          };
          double t[5][3][3];
          int nt = polygonize (val, p.v, t);
          for (int i = 0; i < nt; i++) {
    	color_facet (p);
    	glBegin (GL_POLYGON);
    	for (int j = 0; j < 3; j++) {
    	  coord v = {t[i][j][0], t[i][j][1], t[i][j][2]}, np;
    	  foreach_dimension()
    	    np.x = interp (point, v, n.x);
    	  glnormal3d (view, np.x, np.y, np.z);
    	  color_vertex (p, interp (point, v, col));
    	  glvertex3d (view, x + v.x*Delta, y + v.y*Delta, z + v.z*Delta);
    	}
    	glEnd ();
    	view->ni++;
          }
        }
      }
      if (p.edges)
        draw_lines (view, p.lc, p.lw) {
          foreach_visible (view) {
    	double val[8] = {
    	  v[0,0,0], v[1,0,0], v[1,0,1], v[0,0,1],
    	  v[0,1,0], v[1,1,0], v[1,1,1], v[0,1,1]
    	};
    	double t[5][3][3];
    	int nt = polygonize (val, p.v, t);
    	for (int i = 0; i < nt; i++) {
    	  glBegin (GL_LINE_LOOP);
    	  for (int j = 0; j < 3; j++) {
    	    coord v = {t[i][j][0], t[i][j][1], t[i][j][2]};
    	    glvertex3d (view, x + v.x*Delta, y + v.y*Delta, z + v.z*Delta);
    	  }	
    	  glEnd ();
    	  view->ni++;
    	}
          }
        }
    #endif // dimension > 2
      return true;
    }

    travelling(): moves the camera to a different viewpoint

    • start: starting time of the camera motion.
    • end: time at which the viewpoint should be reached.
    • tx, ty, quat, fov: definition of the target viewpoint.
    struct _travelling {
      double start, end;
      float tx, ty, quat[4], fov;
    };
    
    #define interpo(v)							\
      (!p.v ? v : ((t - p.start)*(p.v) + (p.end - t)*(v))/(p.end - p.start))
    
    void travelling (struct _travelling p)
    {
      static float tx, ty, quat[4], fov;
      static double told = -1.;
      if (told < p.start && t >= p.start) {
        bview * view = get_view();
        tx = view->tx, ty = view->ty, fov = view->fov;
        for (int i = 0; i < 4; i++)
          quat[i] = view->quat[i];
      }
      if (t >= p.start && t <= p.end)
        view (tx = interpo (tx), ty = interpo (ty),
    	  fov = interpo (fov),
    	  quat = {interpo(quat[0]), interpo(quat[1]),
    	          interpo(quat[2]), interpo(quat[3])});
      if (told < p.end && t >= p.end) {
        bview * view = get_view();
        tx = view->tx, ty = view->ty, fov = view->fov;
        for (int i = 0; i < 4; i++)
          quat[i] = view->quat[i];
      }  
      told = t;  
    }
    
    #undef interpo

    draw_string(): draws strings on a separate layer (for annotations)

    • str: string to display.
    • pos: position: “0” bottom left, “1” top left, “2” top right and “3” bottom right (default 0).
    • size: the size of the text, given as the number of characters which can fit within the width of the screen. Default is 40.
    • lc[]: an array of red, green, blue values between 0 and 1 which defines the text color.
    • lw: the line width.
    struct _draw_string {
      char * str;
      int pos;
      float size;
      float lc[3], lw;
    };
    
    trace
    bool draw_string (struct _draw_string p)
    {
      bview * view = draw();
      
      glMatrixMode (GL_PROJECTION);
      glPushMatrix();             
      glLoadIdentity();
    
      glMatrixMode (GL_MODELVIEW);
      glPushMatrix();
      glLoadIdentity();
        
      glColor3f (p.lc[0], p.lc[1], p.lc[2]);
      glLineWidth (view->samples*(p.lw > 0. ? p.lw : 1.));
    
      float width  = gl_StrokeWidth ('1'), height = gl_StrokeHeight();
      if (!p.size)
        p.size = 40;
      float hscale = 2./(p.size*width), vscale = hscale*view->width/view->height;
      float vmargin = width/2.*vscale;
      if (p.pos == 0)
        glTranslatef (-1., -1. + vmargin, 0.);
      else if (p.pos == 1)
        glTranslatef (-1., 1. - height*vscale, 0.);
      else if (p.pos == 2)
        glTranslatef (1. - strlen(p.str)*width*hscale, 1. - height*vscale, 0.);
      else
        glTranslatef (1. - strlen(p.str)*width*hscale, -1. + vmargin, 0.);    
      glScalef (hscale, vscale, 1.);
      gl_StrokeString (p.str); 
      
      glMatrixMode (GL_MODELVIEW);
      glPopMatrix();
      glMatrixMode (GL_PROJECTION);
      glPopMatrix();  
    
      return true;
    }

    Usage