src/draw.h

This file contains the definitions for the drawing primitives of Basilisk View.

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

clear(): removes all objects previously drawn

void clear()
{
  bview * view = get_view();
  if (view->active) {
    if (view->list >= 0)
      glEndList();
    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: 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”.
struct _view_set {
  float tx, ty;
  float fov;
  float quat[4];
  float sx, sy, sz;
  unsigned width, height, samples;
  float bg[3];
  float θ, φ;
  bool relative;
  float res;
  char * camera;
  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}, - π/2., v->quat);
    else if (!strcmp (p.camera, "right"))
      gl_axis_to_quat ((float[]){0,1,0}, π/2., v->quat);
    else if (!strcmp (p.camera, "top"))
      gl_axis_to_quat ((float[]){1,0,0}, - π/2., v->quat);
    else if (!strcmp (p.camera, "bottom"))
      gl_axis_to_quat ((float[]){1,0,0}, π/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}, π, v->quat);
    else if (!strcmp (p.camera, "iso")) {
      gl_axis_to_quat ((float[]){0,1,0}, π/4., v->quat);
      float q[4];
      gl_axis_to_quat ((float[]){1,0,0}, - π/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.θ || p.φ) {
    v->gfsview = false;
    float q[4];
    gl_axis_to_quat ((float[]){1,0,0}, - p.φ, q);
    if (p.relative) {
      float q1[4];
      gl_axis_to_quat ((float[]){0,1,0}, p.θ, q1);
      gl_add_quats(q, q1, q1);
      gl_add_quats(q1, v->quat, v->quat);
    }
    else {
      gl_axis_to_quat ((float[]){0,1,0}, p.θ, v->quat);
      gl_add_quats(q, v->quat, v->quat);
    }
  }

  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();
}

Utility functions

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

@def foreach_visible(view)
foreach_cell() {
#if dimension == 2
  double r = Δ*0.71;
#else // dimension == 3
  double r = Δ*0.87;
#endif
  if (!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()
    }
    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 nxx+nyy+nzz=α

#if dimension == 3
@def foreach_visible_plane(view, n1, alpha1)
coord n = {(n1).x, (n1).y, (n1).z};
double _alpha = (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 (n.x, n.y, n.z);
foreach_cell() {
  double r = Δ*0.87, α = (_alpha - n.x*x - n.y*y - n.z*z)/Δ;
  if (fabs(α) > 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)
{
  for (scalar s in all)
    if (!strcmp (s.name, name))
      return s;
  return (scalar){-1};
}

static void draw_lines (bview * view)
{
  glMatrixMode (GL_PROJECTION);
  glPushMatrix();
  glTranslatef (0., 0., view->lc);
  glColor3f (0., 0., 0.); // default color
  glLineWidth (view->samples);
}

static inline double interp (Point point, coord p, scalar col) {
  struct _interpolate r = { col, x + p.x*Δ, y + p.y*Δ, z + p.z*Δ };
  return interpolate_linear (point, r);
}

draw_vof(): displays VOF-reconstructed interfaces

  • c: the name (as a string) of the Volume-Of-Fluid 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.
  • 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.
  • linear: if true the color will be linearly interpolated for each vertex of the facet.
  • map: the colormap to use. Default is jet.
struct _draw_vof {
  char * c;
  bool edges;
  double larger;
  
  char * color;
  double min, max, spread;
  bool linear;
  colormap map;
};

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;
}

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;
  }

  scalar col = {-1};
  if (p.color && strcmp (p.color, "level")) {
    col = lookup_field (p.color);
    if (col.i < 0) {
      fprintf (stderr, "draw_vof(): no field named '%s'\n", p.color);
      return false;
    }
  }
    
  double cmap[NCMAP][3];
  if (p.color) {
    if (p.min == 0 && p.max == 0) {
      if (col.i < 0) // level
	p.min = 0, p.max = depth();
      else {
	stats s = statsf (col);
	double avg = s.sum/s.volume,
	  spread = (p.spread ? p.spread : 5.)*s.stddev;
	p.min = avg - spread; p.max = avg + spread;
      }
    }
    if (!p.map)
      p.map = jet;
    p.map (cmap);
  }
  
  double cmin = 1e-3; // do not reconstruct fragments smaller than this

#if TREE
  // make sure we prolongate properly
  if (c.prolongation != fraction_refine) {
    c.prolongation = c.refine = fraction_refine;
    boundary ({c});
  }
#endif // TREE
    
  bview * view = draw();
#if dimension == 2
  draw_lines (view);
  glBegin (GL_LINES);
  foreach_visible (view)
    if (cfilter (point, c, cmin)) {
      coord n = mycs (point, c);
      double α = plane_alpha (c[], n);
      coord segment[2];
      if (facets (n, α, segment) == 2) {
	glVertex2d (x + segment[0].x*Δ, y + segment[0].y*Δ);
	glVertex2d (x + segment[1].x*Δ, y + segment[1].y*Δ);
	view->ni++;
      }
    }
  glEnd ();
  glPopMatrix ();
#else // dimension == 3
  double larger = p.larger ? p.larger : p.edges ? 1 : 1.1;
  glColor3f (1., 1., 1.); // default color
  foreach_visible (view)
    if (cfilter (point, c, cmin)) {
      coord n = mycs (point, c);
      double α = plane_alpha (c[], n);
      coord v[12];
      int m = facets (n, α, v, larger);
      if (m > 2) {
	if (p.color && (!p.linear || col.i < 0)) {
	  color b = colormap_color (cmap, col.i < 0 ? level : col[],
				    p.min, p.max);
	  glColor3f (b.r/255., b.g/255., b.b/255.);
	}
	normalize (&n);
	if (view->gfsview)
	  glNormal3d (- n.x, - n.y, - n.z);
	else
	  glNormal3d (n.x, n.y, n.z);
	glBegin (GL_POLYGON);
	for (int i = 0; i < m; i++) {
	  if (p.color && p.linear && col.i >= 0) {
	    color b = colormap_color (cmap, interp (point, v[i], col),
				      p.min, p.max);
	    glColor3f (b.r/255., b.g/255., b.b/255.);
	  }
	  glVertex3d (x + v[i].x*Δ, y + v[i].y*Δ, z + v[i].z*Δ);
	}
	glEnd ();
	view->ni++;
      }
    }
  if (p.edges) {
    draw_lines (view);
    foreach_visible (view)
      if (cfilter (point, c, cmin)) {
	coord n = mycs (point, c);
	double α = plane_alpha (c[], n);
	coord v[12];
	int m = facets (n, α, v, larger);
	if (m > 2) {
	  glBegin (GL_LINE_LOOP);
	  for (int i = 0; i < m; i++)
	    glVertex3d (x + v[i].x*Δ, y + v[i].y*Δ, z + v[i].z*Δ);
	  glEnd ();
	  view->ni++;
	}
      }
    glPopMatrix ();
  }
#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 nxx+nyy+nzz=α

struct _cells {
  coord n;
  double α;
};

trace
void cells (struct _cells p)
{
  bview * view = draw();
  draw_lines (view);
#if dimension == 2
  foreach_visible (view) {
    glBegin (GL_LINE_LOOP);
    glVertex2d (x - Δ/2., y - Δ/2.);
    glVertex2d (x + Δ/2., y - Δ/2.);
    glVertex2d (x + Δ/2., y + Δ/2.);
    glVertex2d (x - Δ/2., y + Δ/2.);
    glEnd();
    view->ni++;
  }
#else // dimension == 3
  foreach_visible_plane (view, p.n, p.α) {
    coord v[12];
    int m = facets (n, α, v, 1.);
    if (m > 2) {
      glBegin (GL_LINE_LOOP);
      for (int i = 0; i < m; i++)
	glVertex3d (x + v[i].x*Δ, y + v[i].y*Δ, z + v[i].z*Δ);
      glEnd ();
      view->ni++;
    }
  }
#endif // dimension == 3
  glPopMatrix ();  
}

squares(): displays colormapped fields

The field name is given by f. The min, max, spread and map 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 nxx+nyy+nzz=α

struct _squares {
  char * f;
  double min, max, spread;
  bool linear;
  colormap map;
  coord n;
  double α;
};

#define gltexcoord(v) glTexCoord1d (clamp(((v) - p.min)/(p.max - p.min), 0, 1))

trace
bool squares (struct _squares p)
{
  scalar f = lookup_field (p.f);
  if (f.i < 0) {
    fprintf (stderr, "squares(): no field named '%s'\n", p.f);
    return false;
  }
  
  if (p.min == 0 && p.max == 0) {
    stats s = statsf (f);
    double avg = s.sum/s.volume, spread = (p.spread ? p.spread : 5.)*s.stddev;
    p.min = avg - spread; p.max = avg + spread;
  }
  if (!p.map)
    p.map = jet;
  double cmap[NCMAP][3];
  p.map (cmap);

  bview * view = draw();

  glShadeModel (GL_SMOOTH);
  if (p.linear) {
    if (view->vector) { // does not use textures for vector graphics
			// (not supported by GL2PS)
#if dimension == 2
      glBegin (GL_QUADS);
      foreach_visible (view)
	if (f[] != nodata) {
	  color c;
	  c = colormap_color (cmap, (f[] + f[-1] + f[-1,-1] + f[0,-1])/4.,
			      p.min, p.max);
	  glColor3f (c.r/255., c.g/255., c.b/255.);
	  glVertex2d (x - Δ/2., y - Δ/2.);
	  c = colormap_color (cmap, (f[] + f[1] + f[1,-1] + f[0,-1])/4.,
			      p.min, p.max);
	  glColor3f (c.r/255., c.g/255., c.b/255.);
	  glVertex2d (x + Δ/2., y - Δ/2.);
	  c = colormap_color (cmap, (f[] + f[1] + f[1,1] + f[0,1])/4.,
			      p.min, p.max);
	  glColor3f (c.r/255., c.g/255., c.b/255.);
	  glVertex2d (x + Δ/2., y + Δ/2.);
	  c = colormap_color (cmap, (f[] + f[-1] + f[-1,1] + f[0,1])/4.,
			      p.min, p.max);
	  glColor3f (c.r/255., c.g/255., c.b/255.);
	  glVertex2d (x - Δ/2., y + Δ/2.);
	  view->ni++;
	}
      glEnd();
#else // dimension == 3
      foreach_visible_plane (view, p.n, p.α)
	if (f[] != nodata) {
	  coord v[12];
	  int m = facets (n, α, v, 1.);
	  if (m > 2) {
	    glBegin (GL_POLYGON);
	    for (int i = 0; i < m; i++) {
	      color c = colormap_color (cmap, interp (point, v[i], f),
					p.min, p.max);
	      glColor3f (c.r/255., c.g/255., c.b/255.);
	      glVertex3d (x + v[i].x*Δ, y + v[i].y*Δ, z + v[i].z*Δ);
	    }
	    glEnd ();
	    view->ni++;
	  }
	}
#endif // dimension == 3
    }
    else { // not vector graphics, better results with textures
      glEnable (GL_TEXTURE_1D);
      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);
      glColor3f (1., 1., 1.);
#if dimension == 2
      foreach_visible (view)
	if (f[] != nodata) {
	  glBegin (GL_TRIANGLE_FAN);
	  gltexcoord ((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 (x, y);
	  gltexcoord ((f[] + f[-1] + f[-1,-1] + f[0,-1])/4.);
	  glVertex2d (x - Δ/2., y - Δ/2.);
	  gltexcoord ((f[] + f[1] + f[1,-1] + f[0,-1])/4.);
	  glVertex2d (x + Δ/2., y - Δ/2.);
	  gltexcoord ((f[] + f[1] + f[1,1] + f[0,1])/4.);
	  glVertex2d (x + Δ/2., y + Δ/2.);
	  gltexcoord ((f[] + f[-1] + f[-1,1] + f[0,1])/4.);
	  glVertex2d (x - Δ/2., y + Δ/2.);
	  gltexcoord ((f[] + f[-1] + f[-1,-1] + f[0,-1])/4.);
	  glVertex2d (x - Δ/2., y - Δ/2.);
	  glEnd();
	  view->ni++;
	}
#else // dimension == 3
      foreach_visible_plane (view, p.n, p.α)
	if (f[] != nodata) {
	  coord v[12];
	  int m = facets (n, α, 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);
	    gltexcoord (interp (point, c, f));
	    glVertex3d (x + c.x*Δ, y + c.y*Δ, z + c.z*Δ);
	    for (int i = 0; i < m; i++) {
	      gltexcoord (interp (point, v[i], f));
	      glVertex3d (x + v[i].x*Δ, y + v[i].y*Δ, z + v[i].z*Δ);
	    }
	    gltexcoord (interp (point, v[0], f));
	    glVertex3d (x + v[0].x*Δ, y + v[0].y*Δ, z + v[0].z*Δ);
	    glEnd ();
	    view->ni++;
	  }
	}
#endif // dimension == 3
      glDisable (GL_TEXTURE_1D);
    }
  }
  else { // !p.linear
#if dimension == 2
    glBegin (GL_QUADS);
    foreach_visible (view)
      if (f[] != nodata) {	
	color c = colormap_color (cmap, f[], p.min, p.max);
	glColor3f (c.r/255., c.g/255., c.b/255.);
	glVertex2d (x - Δ/2., y - Δ/2.);
	glVertex2d (x + Δ/2., y - Δ/2.);
	glVertex2d (x + Δ/2., y + Δ/2.);
	glVertex2d (x - Δ/2., y + Δ/2.);
	view->ni++;
      }
    glEnd();
#else // dimension == 3
    foreach_visible_plane (view, p.n, p.α)
      if (f[] != nodata) {
	coord v[12];
	int m = facets (n, α, v, 1.);
	if (m > 2) {
	  color c = colormap_color (cmap, f[], p.min, p.max);
	  glColor3f (c.r/255., c.g/255., c.b/255.);
	  glBegin (GL_POLYGON);
	  for (int i = 0; i < m; i++)
	    glVertex3d (x + v[i].x*Δ, y + v[i].y*Δ, z + v[i].z*Δ);
	  glEnd ();
	  view->ni++;
	}
      }
#endif // dimension == 3
  }
  return true;
}

box(): displays box boundaries and axis coordinates

trace
bool box()
{
  bview * view = draw();
  draw_lines (view);

  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);
  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 (x - Δ/2., y - Δ/2.);
    glVertex2d (x + Δ/2., y - Δ/2.);
    glVertex2d (x + Δ/2., y + Δ/2.);
    glVertex2d (x - Δ/2., y + Δ/2.);
    glEnd ();
    view->ni++;
  }  
#else // dimension != 2
  foreach_level (0) {
    for (int i = -1; i <= 1; i += 2) {
      glBegin (GL_LINE_LOOP);
      glVertex3d (x - Δ/2., y - Δ/2., z + i*Δ/2.);
      glVertex3d (x + Δ/2., y - Δ/2., z + i*Δ/2.);
      glVertex3d (x + Δ/2., y + Δ/2., z + i*Δ/2.);
      glVertex3d (x - Δ/2., y + Δ/2., z + i*Δ/2.);
      glEnd ();
      view->ni++;
      glBegin (GL_LINES);
      for (int j = -1; j <= 1; j += 2) {
	glVertex3d (x + i*Δ/2., y + j*Δ/2., z - Δ/2.);
	glVertex3d (x + i*Δ/2., y + j*Δ/2., z + Δ/2.);
      }
      glEnd ();
      view->ni++;
    }
  }
#endif // dimension != 2

  glMatrixMode (GL_PROJECTION);
  glPopMatrix();
  return true;
}

isosurface(): displays an isosurface of a field

Work in progress.

struct _isosurface {
  char * f;
  double v;
};

trace
bool isosurface (struct _isosurface p)
{
  scalar f = lookup_field (p.f);
  if (f.i < 0) {
    fprintf (stderr, "isosurface(): no field named '%s'\n", p.f);
    return false;
  }

  bview * view = draw();
  foreach_visible (view) {
    
  }
  return true;
}

Usage