src/view.h

Basilisk View

This module defines functions which compute various graphical representations (“drawings”) of Basilisk fields, including reconstructed Volume-Of-Fluid facets and colorscale representations of cross-sections of scalar fields. These representations are rendered using OpenGL and can be saved in various formats (PPM, Gnuplot, OBJ, KML, PDF, SVG etc.).

Installation

In contrast with other Basilisk modules, this module relies on additional libraries which needs to be installed and linked with the Basilisk program. See gl/INSTALL for instructions.

Usage

A simple example would look like:

...
#include view.h
...
event image (t = 1) {
  clear();
  draw_vof ("f");
  box();
  save ("image.ppm");
}

The clear() function resets the image, draw_vof() and box() add two graphical representations and save() saves the resulting image in PPM format. See User functions for a detailed documentation.

The resulting program needs to be linked with the appropriate libraries using e.g.:

qcc -Wall -O2 program.c -o program \
    -L$BASILISK/gl -lglutils -lfb_osmesa -lGLU -lOSMesa -lm

or

qcc -Wall -O2 program.c -o program \
    -L$BASILISK/gl -lglutils -lfb_glx -lGLU -lGLEW -lGL -lX11

depending on which version of OpenGL should be used.

Implementation

We include the various helper functions defined either by the system or by the Basilisk libraries in gl/.

#if defined(__APPLE__)
#  include <OpenGL/gl.h>
#  include <OpenGL/glu.h>
#else
#  include <GL/gl.h>
#  include <GL/glu.h>
#endif

#include "gl/framebuffer.h"
#include "gl/gl2ps/gl2ps.h"
#include "gl/trackball.h"
#include "gl/utils.h"
#include "utils.h"
#include "input.h"

The bview class

Contains the definition of the current view.

struct _bview {
  float tx, ty, sx, sy, sz;
  float quat[4];
  float fov;

  bool gfsview;   // rotate axis to match gfsview
  bool vector;    // are we doing vector graphics?
  
  float bg[3];
  float lc;
  float res;

  unsigned width, height, samples;

  framebuffer * fb;
  Frustum frustum; // the view frustum
  
  int ni; // number of items drawn

  int list; // display list, for symmetries, not used yet
  bool active;
};

typedef struct _bview bview;

The allocator method.

bview * bview_new() {
  bview * p = calloc (1, sizeof(bview));

  p->tx = p->ty = 0;
  p->sx = p->sy = p->sz = 1.;
  p->quat[0] = p->quat[1] = p->quat[2] = 0; p->quat[3] = 1;
  p->fov = 24.;
  gl_trackball (p->quat, 0.0, 0.0, 0.0, 0.0);

#if dimension <= 2
  p->bg[0] = 1; p->bg[1] = 1; p->bg[2] = 1;
#else
  p->bg[0] = 0.3; p->bg[1] = 0.4; p->bg[2] = 0.6;
#endif
  p->res = 1.;
  p->lc = 0.001;

  p->vector = false;
  
  p->samples = 4;
  p->width = 800*p->samples, p->height = 800*p->samples;

  p->fb = framebuffer_new (p->width, p->height);
  
  init_gl();
  p->list = -1; // glGenLists (1);
  p->active = false;
  
  return p;
}

The destructor method.

void bview_destroy (bview * p)
{
  if (p->list >= 0)
    glDeleteLists (p->list, 1);
  framebuffer_destroy (p->fb);
  free (p);
}

For the moment there is a single (static) current view.

static bview * _view = NULL;

The current view needs to be destroyed when we exit Basilisk. This is done by adding this callback to the free_solver() lists of destructors.

static void destroy_view()
{
  assert (_view);
  bview_destroy (_view);
}

bview * get_view() {
  if (!_view) {
    _view = bview_new();
    free_solver_func_add (destroy_view);
  }
  return _view;
}

The main drawing function.

static void redraw() {
  bview * view = get_view();
    
  /* OpenGL somehow generates floating-point exceptions... turn them off */
  disable_fpe (FE_DIVBYZERO|FE_INVALID);

  //  glViewport (100, 100, view->width, view->height);
  glMatrixMode (GL_PROJECTION);
  glLoadIdentity ();
  double max = 2.;    
#if 0  
  GList * symmetries = get_symmetries (list);
  max = gfs_gl_domain_extent (domain, symmetries);
#endif

  gluPerspective (view->fov, view->width/(float)view->height, 1., 1. + 2.*max);
  glMatrixMode (GL_MODELVIEW);
	    
  glLoadIdentity ();
  glTranslatef (view->tx, view->ty, - (1. + max));
  //  glTranslatef (0, 0, - (1. + max));
    
  GLfloat m[4][4];
  gl_build_rotmatrix (m, view->quat);
  glMultMatrixf (&m[0][0]);
  
  if (view->gfsview) { // rotate to match gfsview parameters
    m[0][0] = 0., m[0][1] =  0., m[0][2] =  -1.;
    m[1][0] = 0., m[1][1] = -1., m[1][2] =   0.;
    m[2][0] = 1., m[2][1] =  0., m[2][2] =   0.;
    glMultMatrixf (&m[0][0]);
  }
  
  glScalef (view->sx/L0, view->sy/L0, view->sz/L0);

  glClearColor (view->bg[0], view->bg[1], view->bg[2], 0.);
  glClear (GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);

  gl_get_frustum (&view->frustum);
  
  view->active = true;
  view->ni = 0;
}

This is called by graphics primitives before drawing.

bview * draw() {
  bview * view = get_view();
  if (!view->active) {
    redraw();
    if (view->list >= 0)
      glNewList (view->list, GL_COMPILE);
  }
  return view;
}

Helper function for parallel image composition

compose_image() returns an image buffer made by composition of the framebuffer images on each of the MPI processes.

typedef void * pointer; // fixme: trace is confused by pointers

#if !_MPI
trace
static pointer compose_image (bview * view) {
  return framebuffer_image((view)->fb);
}
#else // _MPI
#if dimension <= 2
typedef struct {
  GLubyte a[4];
} RGBA;

static void compose_image_op (void * pin, void * pout, int * len,
			      MPI_Datatype * dptr)
{
  RGBA * in = pin, * out = pout;
  for (int i = 0; i < *len; i++,in++,out++)
    if (out->a[3] == 0)
      *out = *in;
}

trace
static pointer compose_image (bview * view)
{
  unsigned char * image = framebuffer_image (view->fb);
  if (npe() > 1) {
    MPI_Op op;
    MPI_Op_create (compose_image_op, true, &op);    
    MPI_Datatype rgba;
    MPI_Type_contiguous (4, MPI_BYTE, &rgba);
    MPI_Type_commit (&rgba);
    int size = view->width*view->height;
    if (pid() == 0)
      MPI_Reduce (MPI_IN_PLACE, image, size, rgba, op, 0, MPI_COMM_WORLD);
    else
      MPI_Reduce (image, image, size, rgba, op, 0, MPI_COMM_WORLD);
    MPI_Op_free (&op);
    MPI_Type_free (&rgba);
  }
  return image;
}
#else /* 3D */
typedef struct {
  GLubyte a[4];
  float depth;
} RGBA;

static void compose_image_op (void * pin, void * pout, int * len,
			      MPI_Datatype * dptr)
{
  RGBA * in = pin, * out = pout;
  for (int i = 0; i < *len; i++,in++,out++)
    if (out->depth > in->depth)
      *out = *in;
}

trace
static pointer compose_image (bview * view)
{
  unsigned char * image = framebuffer_image (view->fb);
  if (npe() > 1) {
    MPI_Op op;
    MPI_Op_create (compose_image_op, true, &op);
    MPI_Datatype rgba;
    MPI_Type_create_struct (2,
			    (int[]){4,1},
			    (MPI_Aint[]){0,4},
			    (MPI_Datatype[]){MPI_BYTE, MPI_FLOAT},
			    &rgba);
    MPI_Type_commit (&rgba);
    fbdepth_t * depth = framebuffer_depth (view->fb);
    int size = view->width*view->height;
    RGBA * buf = malloc (size*sizeof(RGBA));
    unsigned char * ptr = image;
    fbdepth_t * dptr = depth;
    for (int i = 0; i < size; i++) {
      for (int j = 0; j < 4; j++)
	buf[i].a[j] = *ptr++;
      buf[i].depth = *dptr++;
    }
    if (pid() == 0) {
      MPI_Reduce (MPI_IN_PLACE, buf, size, rgba, op, 0, MPI_COMM_WORLD);
      unsigned char * ptr = image;
      for (int i = 0; i < size; i++)
	for (int j = 0; j < 4; j++)
	  *ptr++ = buf[i].a[j];
    }
    else
      MPI_Reduce (buf, buf, size, rgba, op, 0, MPI_COMM_WORLD);
    free (buf);
    MPI_Op_free (&op);
    MPI_Type_free (&rgba);
  }
  return image;
}
#endif /* 3D */
#endif /* _MPI */

User functions

Drawing user functions are defined in draw.h.

#include "draw.h"

load(): read drawing commands from a file or buffer

The commands are calls of user functions. They can be read from a file defined by fp or file, or from the memory buffer buf.

If history is not NULL, the commands are appended to this buffer.

Besides the load(), save() and drawing functions defined in draw.h, valid drawing commands also include:

  • restore()
  • dump()
  • input_gfs()
  • display(): forces redrawing.
  • show(): writes on standard error the content of the current command history.
  • quit(): stops parsing commands.
struct _load {
  FILE * fp;   // read commands from this file
  char * file; // read commands from this file
  Array * buf; // read commands from this buffer
  Array * history; // append history to this one
};

bool load (struct _load p);

save(): saves drawing to a file in various formats

The file to write to is given either using its name file or the file pointer fp. If neither is specified, the default is stdout.

The format to use is given either explicitly using format, or, if a file name is given, using the file name extension (i.e. .ppm, .bv, etc.). If neither is specified, the default is “ppm”.

For vector graphics, the base line width can be specified using lw. The default is one.

The recognised file formats are:

  • “ppm”: Portable PixMap format. A basic uncompressed image format.
  • “png”, “jpg”: Compressed image formats. Will only work if the convert command from ImageMagick is installed on the system.
  • “bv”: Basilisk View format. Saves all Basilisk function calls necessary to reproduce the figure. Use together with load().
  • “gnu”: Gnuplot format. Saves a vector graphics (3D) representation of the objects.
  • “obj”: Wavefront 3D object format. Can be read by a number of 3D visualisation tools.
  • “kml”: Keyhole Markup Language. Can be used with Google Earth and other GIS.
  • “ps”, “eps”, “tex”, “pdf”, “svg”, “pgf”: the various vector graphics formats supported by gl2ps.

Note that MPI-parallel output is only implemented for the “ppm” format at the moment.

struct _save {
  char * file, * format;
  FILE * fp;
  float lw; /* base line width for vector drawings */
  int sort, options;

  Array * history; // command history
  bview * view;
};

static void bview_draw (bview * view)
{
  if (!view->active)
    return;
  
  if (view->list >= 0)
    glEndList();
  view->active = false;
    
#if 0  
  GfsFrustum frustum;
  gfs_gl_get_frustum (view, symmetries, &frustum);
  
  GList * i = list;
  while (i) {
    GfsGl * gl = i->data;

    if (GFS_IS_GL_CLIP_PLANE (gl)) {
      gl->format = p->format;
      gfs_gl_clip_plane_disable (GFS_GL_CLIP_PLANE (gl));
    }
    i = i->next;
  }
  i = list;
  while (i) {
    if (GFS_IS_GL_CUT_PLANE (i->data))
      GFS_GL_CUT_PLANE (i->data)->list = list;
    i = i->next;
  }

  GSList * clip = NULL;
  gboolean firstclip = TRUE;
  i = list;
  while (i) {
    GfsGl * gl = i->data;
    gl->format = p->format;
    if (GFS_IS_GL_CLIP_PLANE (gl)) {
      if (firstclip) {
	g_slist_foreach (clip, (GFunc) gfs_gl_clip_plane_disable, NULL);
	g_slist_free (clip); clip = NULL;
	firstclip = FALSE;	  
      }
      gfs_gl_draw (gl, &frustum);
      clip = g_slist_prepend (clip, gl);
    }
    else {
      gfs_gl_draw (gl, &frustum);
      firstclip = TRUE;
    }
    i = i->next;
  }
  g_slist_free (clip);
  glEndList();
  
  gfs_gl_symmetry_apply (symmetries, display_list);
  gfs_gl_frustum_free (&frustum);
  g_list_free (symmetries);
#else
  if (view->list >= 0)
    glCallList (view->list);
#endif
  glFinish ();

  enable_fpe (FE_DIVBYZERO|FE_INVALID);
}

static void redraw_feedback (struct _save * p)
{
  bview * view = p->view ? p->view : get_view();
  assert (p->history);
  if (p->history->len) {
    float res = view->res;
    view->res = 0.;
    view->vector = true; // vector graphics
    redraw();
    // we use a display list, just as a workaround for some buggy
    // feedback buffer implementations
    view->list = glGenLists (1);
    glNewList (view->list, GL_COMPILE);
    load (buf = p->history);
    bview_draw (view);
    glDeleteLists (view->list, 1);
    view->list = -1;
    view->vector = false;
    view->res = res;
  }
}

#define MAXBUFFSIZE (1 << 28) // 1 GB of feedback buffer

trace
bool save (struct _save p)
{
  if (p.file && (p.fp = fopen (p.file, "w")) == NULL) {
    perror (p.file);
    return false;
  }
  if (!p.fp)
    p.fp = stdout;

  char ppm[] = "ppm";
  if (!p.format) {
    p.format = ppm;
    if (p.file) {
      char * s = strchr (p.file, '.'), * dot = s;
      while (s) {
	dot = s;
	s = strchr (s + 1, '.');
      }
      if (dot)
	p.format = dot + 1;
    }
  }

  bview * view = p.view ? p.view : get_view();

  if (!strcmp (p.format, "ppm")) {
    bview_draw (view);
    unsigned char * image = compose_image (view);
    if (pid() == 0)
      gl_write_image (p.fp, image, view->width, view->height, view->samples);    
  }

  else if (!strcmp (p.format, "png") ||
	   !strcmp (p.format, "jpg")) {
    bview_draw (view);
    unsigned char * image = compose_image (view);
    if (pid() == 0) {
      char command[strlen ("convert ppm:- ") + strlen (p.file) + 1];
      strcpy (command, "convert ppm:- ");
      strcat (command, p.file);
      FILE * fp = popen (command, "w");
      gl_write_image (fp, image, view->width, view->height, view->samples);
      pclose (fp);
    }
  }

  else if (!strcmp (p.format, "bv")) {
    assert (p.history);
    fprintf (p.fp,
	     "view (fov = %g, quat = {%g,%g,%g,%g}, "
	     "tx = %g, ty = %g, "
	     "bg = {%g,%g,%g}, "
	     "width = %d, height = %d, samples = %d"
	     ");\n",
	     view->fov,
	     view->quat[0], view->quat[1], view->quat[2], view->quat[3],
	     view->tx, view->ty,
	     view->bg[0], view->bg[1], view->bg[2],
	     view->width/view->samples, view->height/view->samples,
	     view->samples);
    fwrite (p.history->p, 1, p.history->len, p.fp);
  }
  
  else if (!strcmp (p.format, "gnu") ||
	   !strcmp (p.format, "obj") ||
	   !strcmp (p.format, "kml")) {
    int format = (!strcmp (p.format, "gnu") ? FEEDBACK_GNU :
		  !strcmp (p.format, "obj") ? FEEDBACK_OBJ :
		  !strcmp (p.format, "kml") ? FEEDBACK_KML :
		  -1);
    unsigned buffsize = 1 << 24;
    bool done = false;
    while (!done && buffsize <= MAXBUFFSIZE) {
      float * f = gl_feedback_begin (buffsize);
      redraw_feedback (&p);
      done = gl_feedback_end (f, p.fp, format);
      buffsize *= 2;
    }
    if (!done)
      fprintf (stderr, "save(): error: exceeded maximum feedback buffer size\n");
  }
  
  else if (!strcmp (p.format, "ps") ||
	   !strcmp (p.format, "eps") ||
	   !strcmp (p.format, "tex") ||
	   !strcmp (p.format, "pdf") ||
	   !strcmp (p.format, "svg") ||
	   !strcmp (p.format, "pgf")) {
    GLint format = (!strcmp (p.format, "ps") ? GL2PS_PS :
		    !strcmp (p.format, "eps") ? GL2PS_EPS :
		    !strcmp (p.format, "tex") ? GL2PS_TEX :
		    !strcmp (p.format, "pdf") ? GL2PS_PDF :
		    !strcmp (p.format, "svg") ? GL2PS_SVG :
		    !strcmp (p.format, "pgf") ? GL2PS_PGF :
		    -1);
    GLint state = GL2PS_OVERFLOW;
    GLint sort = p.sort ? p.sort : GL2PS_SIMPLE_SORT;
    GLint options = p.options ? p.options : (GL2PS_SIMPLE_LINE_OFFSET |
					     GL2PS_SILENT |
					     GL2PS_BEST_ROOT |
					     GL2PS_OCCLUSION_CULL |
					     GL2PS_USE_CURRENT_VIEWPORT |
					     GL2PS_TIGHT_BOUNDING_BOX);
    unsigned buffsize = 1 << 24;
    while (state == GL2PS_OVERFLOW && buffsize <= MAXBUFFSIZE) {
      gl2psBeginPage ("", "bview",
		      NULL,
		      format, sort, options, 
		      GL_RGBA, 0, NULL, 
		      0, 0, 0,
		      buffsize, p.fp, "");
      redraw_feedback (&p);
      disable_fpe (FE_DIVBYZERO|FE_INVALID);
      state = gl2psEndPage();
      enable_fpe (FE_DIVBYZERO|FE_INVALID);
      buffsize *= 2;
    }
    if (state == GL2PS_OVERFLOW)
      fprintf (stderr, "save(): error: exceeded maximum feedback buffer size\n");
  }

  else {
    fprintf (stderr, "save(): unknown format '%s'\n", p.format);
    if (p.file) {
      fclose (p.fp);
      remove (p.file);
    }
    return false;
  }

  fflush (p.fp);
  if (p.file)
    fclose (p.fp);

  return true;
}

Implementation of the load() function.

The functions below parse a text file and perform the corresponding function calls.

static char * remove_blanks (char * line)
{
  while (strchr (" \t", *line)) line++;
  char * s = line, * cur = line;
  bool instring = false;
  while (*s != '\0' && *s != '#') {
    if (*s == '"')
      instring = !instring;
    if (instring || !strchr (" \t", *s))
      *cur++ = *s;
    s++;
  }
  *cur = '\0';
  return line;
}

static void fields_stats()
{
  fprintf (stderr, "# t = %g, fields = {", t);
  for (scalar s in all)
    fprintf (stderr, " %s", s.name);
  fputs (" }\n", stderr);
}

static void draw_append (char * buf, Array * history, FILE * interactive)
{
  if (interactive) {
    if (history->len)
      load (buf = history);
    save (fp = interactive);
  }
  array_append (history, buf, strlen(buf)*sizeof(char));
}

The draw_get.h file is generated automatically by params.awk and contains parsing commands for the functions defined in draw.h.

#include "draw_get.h"

static bool process_line (char * line, Array * history, FILE * interactive)
{
  if (line[0] == '\0')
    return true;
  char * buf = strdup (line);
  char * s = strtok (remove_blanks (line), "(");
  if (!s) {
    free (buf);
    return true;
  }

  if (!strcmp (s, "restore")) {
    char * file = NULL;
    parse_params ((Params[]){{"file", pstring, &file}, {NULL}});
    if (file) {
      if (!restore (file = file, list = all))
	fprintf (stderr, "could not restore from '%s'\n", file);
      else {
	fields_stats();
	clear();
	// rebuild display list using history
	if (history->len && load (buf = history) && interactive)
	  save (fp = interactive);
      }
    }
  }

  else if (!strcmp (s, "dump")) {
    char * file = NULL;
    parse_params ((Params[]){{"file", pstring, &file}, {NULL}});
    dump (file = file);
  }
  
  else if (!strcmp (s, "input_gfs")) {
    char * file = NULL;
    parse_params ((Params[]){{"file", pstring, &file}, {NULL}});
    if (file) {
      input_gfs (file = file, list = all);
      fields_stats();
      clear();
      // rebuild display list using history
      if (history->len && load (buf = history) && interactive)
	save (fp = interactive);
    }
  }

  else if (!strcmp (s, "save")) {
    char * file = NULL;
    parse_params ((Params[]){{"file", pstring, &file}, {NULL}});
    if (file)
      save (file = file, history = history);
  }

  else if (!strcmp (s, "load")) {
    char * file = NULL;
    parse_params ((Params[]){{"file", pstring, &file}, {NULL}});
    if (file && load (file = file, history = history) && interactive) {
      load (buf = history);
      save (fp = interactive);
    }
  }
        
  else if (!strcmp (s, "cells")) {
    struct _cells p = {{0}};
    _cells_get (&p);
    cells (p);
    draw_append (buf, history, interactive);
  }
    
  else if (!strcmp (s, "draw_vof")) {
    struct _draw_vof p = {0};
    _draw_vof_get (&p);
    if (draw_vof (p))
      draw_append (buf, history, interactive);
  }
    
  else if (!strcmp (s, "squares")) {
    struct _squares p = {0};
    _squares_get (&p);
    squares (p);
    draw_append (buf, history, interactive);
  }
    
  else if (!strcmp (s, "display")) {
    if (interactive && history->len && load (buf = history))
      save (fp = interactive);
  }

  else if (!strcmp (s, "clear")) {
    clear();
    if (interactive)
      save (fp = interactive);
    history->len = 0;
  }

  else if (!strcmp (s, "show")) {
    if (interactive && history->len)
      save (fp = stderr, format = "bv", history = history);
  }

  else if (!strcmp (s, "box")) {
    box();
    draw_append (buf, history, interactive);
  }

  else if (!strcmp (s, "view")) {
    struct _view_set p = {0};
    _view_set_get (&p);
    view (p);
    if (p.width || p.height || p.samples) {
      // rebuild display list using history
      if (history->len && load (buf = history) && p.samples && interactive)
	save (fp = interactive);
    }
  }

  else if (!strcmp (s, "quit")) {
    free (buf);
    return false; // quit
  }
  
  else if (s[0] != '\n')
    fprintf (stderr, "load(): syntax error: '%s'\n", s);

  free (buf);
  return true;
}

bool load (struct _load p) {
  if (p.file) {
    p.fp = fopen (p.file, "r");
    if (!p.fp) {
      perror (p.file);
      return false;
    }
  }

  Array * history = array_new();
  if (p.fp) { // read lines from file
    char line[256];
    while (fgets (line, 256, p.fp) && process_line (line, history, NULL));
  }
  else if (p.buf) { // read lines from buffer
    int i = 0;
    char * s = p.buf->p;
    while (i < p.buf->len) {
      char * start = s;
      while (i < p.buf->len && *s != '\n')
	s++, i++;
      if (*s == '\n' && ++s > start) {
	char line[s - start + 1];
	strncpy (line, start, s - start);
	line[s - start] = '\0';
	process_line (line, history, NULL);
      }
    }
  }
  if (p.history)
    array_append (p.history, history->p, history->len);
  array_free (history);

  return true;
}

Usage

Examples

Tests