"Fossies" - the Fresh Open Source Software Archive

Member "gfsview-snapshot-121130/gl/gfsgl.c" (30 Nov 2012, 126980 Bytes) of package /linux/privat/gfsview-snapshot-121130.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "gfsgl.c" see the Fossies "Dox" file reference documentation.

    1 /* Gerris - The GNU Flow Solver
    2  * Copyright (C) 2004-2012 National Institute of Water and Atmospheric
    3  * Research
    4  *
    5  * This program is free software; you can redistribute it and/or
    6  * modify it under the terms of the GNU General Public License as
    7  * published by the Free Software Foundation; either version 2 of the
    8  * License, or (at your option) any later version.
    9  *
   10  * This program is distributed in the hope that it will be useful,
   11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
   12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
   13  * General Public License for more details.
   14  *
   15  * You should have received a copy of the GNU General Public License
   16  * along with this program; if not, write to the Free Software
   17  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
   18  * 02111-1307, USA.  
   19  */
   20 
   21 #include "gfsgl.h"
   22 #include <gerris/isocube.h>
   23 #include <gerris/river.h>
   24 #include "trackball.h"
   25 #include "config.h"
   26 #if HAVE_FTGL
   27 # include <FTGL/ftgl.h>
   28 #endif
   29 
   30 #define DEFAULT_FONT (PACKAGE_DATA_DIR "/fonts/Garuda.ttf")
   31 #define DEFAULT_FONT_SIZE 72.
   32 
   33 static gpointer default_font = NULL, default_raster_font = NULL;
   34 
   35 #ifndef GL_CONSTANT
   36 # define GL_CONSTANT GL_CONSTANT_EXT
   37 #endif
   38 
   39 static void gl_isoline_update_levels (GfsGl * gl);
   40 
   41 typedef struct {
   42   GfsVariable * min;
   43   GfsVariable * max;
   44   gdouble level;
   45   GfsGl * gl;
   46 } IsoParams;
   47 
   48 #if FTT_2D
   49 #  include "gfsgl2D.h"
   50 #else  /* 3D */
   51 #  include "gfsgl3D.h"
   52 #endif /* 3D */
   53 #if defined(__APPLE__)
   54 #  include <OpenGL/glu.h>
   55 #else
   56 #  include <GL/glu.h>
   57 #endif
   58 
   59 static void draw_face (GtsTriangle * t)
   60 {
   61   GtsVertex * v1, * v2, * v3;
   62   gts_triangle_vertices (t, &v1, &v2, &v3);
   63   glNormal3d (GTS_POINT (v1)->x, GTS_POINT (v1)->y, GTS_POINT (v1)->z);
   64   glVertex3d (GTS_POINT (v1)->x, GTS_POINT (v1)->y, GTS_POINT (v1)->z);
   65   glNormal3d (GTS_POINT (v2)->x, GTS_POINT (v2)->y, GTS_POINT (v2)->z);
   66   glVertex3d (GTS_POINT (v2)->x, GTS_POINT (v2)->y, GTS_POINT (v2)->z);
   67   glNormal3d (GTS_POINT (v3)->x, GTS_POINT (v3)->y, GTS_POINT (v3)->z);
   68   glVertex3d (GTS_POINT (v3)->x, GTS_POINT (v3)->y, GTS_POINT (v3)->z);
   69 }
   70 
   71 static void draw_sphere (void)
   72 {
   73   glShadeModel (GL_SMOOTH);
   74   glMatrixMode (GL_MODELVIEW);
   75   glPushMatrix ();
   76   glScalef (0.5, 0.5, 0.5);
   77   glBegin (GL_TRIANGLES);
   78   static GtsSurface * sphere = NULL;
   79   if (!sphere) {
   80     sphere = gts_surface_new (gts_surface_class (), gts_face_class (), 
   81                   gts_edge_class (), gts_vertex_class ());
   82     gts_surface_generate_sphere (sphere, 3);
   83   }
   84   gts_surface_foreach_face (sphere, (GtsFunc) draw_face, NULL);
   85   glEnd ();
   86   glPopMatrix ();
   87 }
   88 
   89 FILE * gfs_gl_popen (const gchar * fname)
   90 {
   91   g_return_val_if_fail (fname != NULL, NULL);
   92 
   93   FILE * fptr = fopen (fname, "r");
   94   if (fptr == NULL)
   95     return NULL;
   96   fclose (fptr);
   97 
   98   gchar * script = g_strconcat ("if gunzip -t \"", fname, "\" 2> /dev/null; then "
   99                 "  gunzip -c \"", fname, "\" 2> /dev/null; else ",
  100                 "  cat \"", fname, "\"; fi", NULL);
  101   fptr = popen (script, "r");
  102   g_free (script);
  103   return fptr;
  104 }
  105 
  106 static GtsColor * color_new (gdouble r, gdouble g, gdouble b)
  107 {
  108   GtsColor * c = g_malloc (sizeof (GtsColor));
  109   c->r = r; c->g = g; c->b = b;
  110   return c;
  111 }
  112 
  113 static void color_destroy (GtsColor * color)
  114 {
  115   g_return_if_fail (color != NULL);
  116 
  117   g_free (color);
  118 }
  119 
  120 static void colormap_set_texture (GfsColormap * cmap)
  121 {
  122   guint i;
  123 
  124   for (i = 0; i < GFS_COLORMAP_TEXTURE_SAMPLES; i++) {
  125     GtsColor c = gfs_colormap_color (cmap, i/(gdouble) (GFS_COLORMAP_TEXTURE_SAMPLES - 1));
  126     cmap->texture[3*i] = c.r;
  127     cmap->texture[3*i + 1] = c.g;
  128     cmap->texture[3*i + 2] = c.b;
  129   }
  130 }
  131 
  132 GfsColormap * gfs_colormap_jet (void)
  133 {
  134   GfsColormap * cmap = g_malloc (sizeof (GfsColormap));
  135   gint i;
  136 
  137   cmap->reversed = FALSE;
  138   cmap->colors = g_ptr_array_new ();
  139   for (i = 0; i < 127; i++) {
  140     gdouble r = 
  141       i <= 46 ? 0. : 
  142       i >= 111 ? -0.03125*(i - 111) + 1. :
  143       i >= 78 ? 1. : 
  144       0.03125*(i - 46);
  145     gdouble g = 
  146       i <= 14 || i >= 111 ? 0. : 
  147       i >= 79 ? -0.03125*(i - 111) : 
  148       i <= 46 ? 0.03125*(i - 14) : 
  149       1.;
  150     gdouble b =
  151       i >= 79 ? 0. :
  152       i >= 47 ? -0.03125*(i - 79) :
  153       i <= 14 ? 0.03125*(i - 14) + 1.:
  154       1.;
  155 
  156     g_ptr_array_add (cmap->colors, color_new (r, g, b));
  157   }
  158   colormap_set_texture (cmap);
  159   cmap->name = g_strdup ("Jet");
  160   return cmap;
  161 }
  162 
  163 GfsColormap * gfs_colormap_cool_warm (void)
  164 {
  165   GfsColormap * cmap = g_malloc (sizeof (GfsColormap));
  166   gint i;
  167 
  168   /* diverging cool-warm from:
  169    *  http://www.sandia.gov/~kmorel/documents/ColorMaps/CoolWarmFloat33.csv
  170    * see also:
  171    *  Diverging Color Maps for Scientific Visualization (Expanded)
  172    *  Kenneth Moreland
  173    */
  174   static double basemap[33][3] = {  
  175     {0.2298057,   0.298717966, 0.753683153},
  176     {0.26623388,  0.353094838, 0.801466763},
  177     {0.30386891,  0.406535296, 0.84495867},
  178     {0.342804478, 0.458757618, 0.883725899},
  179     {0.38301334,  0.50941904,  0.917387822},
  180     {0.424369608, 0.558148092, 0.945619588},
  181     {0.46666708,  0.604562568, 0.968154911},
  182     {0.509635204, 0.648280772, 0.98478814},
  183     {0.552953156, 0.688929332, 0.995375608},
  184     {0.596262162, 0.726149107, 0.999836203},
  185     {0.639176211, 0.759599947, 0.998151185},
  186     {0.681291281, 0.788964712, 0.990363227},
  187     {0.722193294, 0.813952739, 0.976574709},
  188     {0.761464949, 0.834302879, 0.956945269},
  189     {0.798691636, 0.849786142, 0.931688648},
  190     {0.833466556, 0.860207984, 0.901068838},
  191     {0.865395197, 0.86541021,  0.865395561},
  192     {0.897787179, 0.848937047, 0.820880546},
  193     {0.924127593, 0.827384882, 0.774508472},
  194     {0.944468518, 0.800927443, 0.726736146},
  195     {0.958852946, 0.769767752, 0.678007945},
  196     {0.96732803,  0.734132809, 0.628751763},
  197     {0.969954137, 0.694266682, 0.579375448},
  198     {0.966811177, 0.650421156, 0.530263762},
  199     {0.958003065, 0.602842431, 0.481775914},
  200     {0.943660866, 0.551750968, 0.434243684},
  201     {0.923944917, 0.49730856,  0.387970225},
  202     {0.89904617,  0.439559467, 0.343229596},
  203     {0.869186849, 0.378313092, 0.300267182},
  204     {0.834620542, 0.312874446, 0.259301199},
  205     {0.795631745, 0.24128379,  0.220525627},
  206     {0.752534934, 0.157246067, 0.184115123},
  207     {0.705673158, 0.01555616,  0.150232812} 
  208   };
  209   
  210   cmap->reversed = FALSE;
  211   cmap->colors = g_ptr_array_new ();
  212   for (i = 0; i < 33; i++)
  213     g_ptr_array_add (cmap->colors, color_new (basemap[i][0], basemap[i][1], basemap[i][2]));
  214   colormap_set_texture (cmap);
  215   cmap->name = g_strdup ("Cool");
  216   return cmap;
  217 }
  218 
  219 GfsColormap * gfs_colormap_gray (void)
  220 {
  221   GfsColormap * cmap = g_malloc (sizeof (GfsColormap));
  222   cmap->reversed = FALSE;
  223   cmap->colors = g_ptr_array_new ();
  224   g_ptr_array_add (cmap->colors, color_new (0,0,0));
  225   g_ptr_array_add (cmap->colors, color_new (1,1,1));
  226   colormap_set_texture (cmap);
  227   cmap->name = g_strdup ("Gray");
  228   return cmap;
  229 }
  230 
  231 void gfs_colormap_destroy (GfsColormap * colormap)
  232 {
  233   guint i;
  234 
  235   g_return_if_fail (colormap != NULL);
  236 
  237   for (i = 0; i < colormap->colors->len; i++)
  238     color_destroy (colormap->colors->pdata[i]);
  239   g_ptr_array_free (colormap->colors, TRUE);
  240   g_free (colormap->name);
  241   g_free (colormap);
  242 }
  243 
  244 GtsColor gfs_colormap_color (GfsColormap * cmap, gdouble val)
  245 {
  246   GtsColor c = {1., 1., 1.}, * c1, * c2;
  247   guint i, n;
  248   gdouble coef;
  249 
  250   g_return_val_if_fail (cmap != NULL, c);
  251 
  252   if (val > 1.0) val = 1.0;
  253   else if (val < 0.0) val = 0.0;
  254   if (cmap->reversed)
  255     val = 1.0 - val;
  256 
  257   n = cmap->colors->len;
  258   if (n == 0)
  259     return c;
  260   if (n == 1)
  261     return *((GtsColor *)cmap->colors->pdata[0]);
  262 
  263   i = floor ((gdouble)val*(gdouble)(n - 1));
  264   if (i == n - 1)
  265     return *((GtsColor *)cmap->colors->pdata[cmap->colors->len - 1]);
  266   coef = val*(gdouble)(n - 1) - (gdouble)i;
  267   c1 = cmap->colors->pdata[i];
  268   c2 = cmap->colors->pdata[i+1];
  269   c.r = c1->r + coef*(c2->r - c1->r);
  270   c.g = c1->g + coef*(c2->g - c1->g);
  271   c.b = c1->b + coef*(c2->b - c1->b);
  272 
  273   return c;
  274 }
  275 
  276 void gfs_colormap_texture (GfsColormap * cmap)
  277 {
  278   g_return_if_fail (cmap != NULL);
  279 
  280   glTexImage1D (GL_TEXTURE_1D, 0, GL_RGB, GFS_COLORMAP_TEXTURE_SAMPLES, 0, GL_RGB, GL_FLOAT,
  281         cmap->texture);
  282   glTexParameteri (GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
  283   glTexParameteri (GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
  284   glTexParameteri (GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
  285   glTexParameteri (GL_TEXTURE_1D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
  286 }
  287 
  288 #define RC(r,c) m[(r)+(c)*4]
  289 #define RCM(m,r,c) (m)[(r)+(c)*4]
  290 
  291 static void matrix_multiply (float * m, float * n)
  292 {
  293   float o[16];
  294   guint i;
  295   
  296   for (i = 0; i < 16; i++) o[i] = m[i];
  297   RC(0,0)=RCM(o,0,0)*RCM(n,0,0)+RCM(o,0,1)*RCM(n,1,0)+
  298           RCM(o,0,2)*RCM(n,2,0)+RCM(o,0,3)*RCM(n,3,0);
  299   RC(0,1)=RCM(o,0,0)*RCM(n,0,1)+RCM(o,0,1)*RCM(n,1,1)+
  300           RCM(o,0,2)*RCM(n,2,1)+RCM(o,0,3)*RCM(n,3,1);
  301   RC(0,2)=RCM(o,0,0)*RCM(n,0,2)+RCM(o,0,1)*RCM(n,1,2)+
  302           RCM(o,0,2)*RCM(n,2,2)+RCM(o,0,3)*RCM(n,3,2);
  303   RC(0,3)=RCM(o,0,0)*RCM(n,0,3)+RCM(o,0,1)*RCM(n,1,3)+
  304           RCM(o,0,2)*RCM(n,2,3)+RCM(o,0,3)*RCM(n,3,3);
  305   RC(1,0)=RCM(o,1,0)*RCM(n,0,0)+RCM(o,1,1)*RCM(n,1,0)+
  306           RCM(o,1,2)*RCM(n,2,0)+RCM(o,1,3)*RCM(n,3,0);
  307   RC(1,1)=RCM(o,1,0)*RCM(n,0,1)+RCM(o,1,1)*RCM(n,1,1)+
  308           RCM(o,1,2)*RCM(n,2,1)+RCM(o,1,3)*RCM(n,3,1);
  309   RC(1,2)=RCM(o,1,0)*RCM(n,0,2)+RCM(o,1,1)*RCM(n,1,2)+
  310           RCM(o,1,2)*RCM(n,2,2)+RCM(o,1,3)*RCM(n,3,2);
  311   RC(1,3)=RCM(o,1,0)*RCM(n,0,3)+RCM(o,1,1)*RCM(n,1,3)+
  312           RCM(o,1,2)*RCM(n,2,3)+RCM(o,1,3)*RCM(n,3,3);
  313   RC(2,0)=RCM(o,2,0)*RCM(n,0,0)+RCM(o,2,1)*RCM(n,1,0)+
  314           RCM(o,2,2)*RCM(n,2,0)+RCM(o,2,3)*RCM(n,3,0);
  315   RC(2,1)=RCM(o,2,0)*RCM(n,0,1)+RCM(o,2,1)*RCM(n,1,1)+
  316           RCM(o,2,2)*RCM(n,2,1)+RCM(o,2,3)*RCM(n,3,1);
  317   RC(2,2)=RCM(o,2,0)*RCM(n,0,2)+RCM(o,2,1)*RCM(n,1,2)+
  318           RCM(o,2,2)*RCM(n,2,2)+RCM(o,2,3)*RCM(n,3,2);
  319   RC(2,3)=RCM(o,2,0)*RCM(n,0,3)+RCM(o,2,1)*RCM(n,1,3)+
  320           RCM(o,2,2)*RCM(n,2,3)+RCM(o,2,3)*RCM(n,3,3);
  321   RC(3,0)=RCM(o,3,0)*RCM(n,0,0)+RCM(o,3,1)*RCM(n,1,0)+
  322           RCM(o,3,2)*RCM(n,2,0)+RCM(o,3,3)*RCM(n,3,0);
  323   RC(3,1)=RCM(o,3,0)*RCM(n,0,1)+RCM(o,3,1)*RCM(n,1,1)+
  324           RCM(o,3,2)*RCM(n,2,1)+RCM(o,3,3)*RCM(n,3,1);
  325   RC(3,2)=RCM(o,3,0)*RCM(n,0,2)+RCM(o,3,1)*RCM(n,1,2)+
  326           RCM(o,3,2)*RCM(n,2,2)+RCM(o,3,3)*RCM(n,3,2);
  327   RC(3,3)=RCM(o,3,0)*RCM(n,0,3)+RCM(o,3,1)*RCM(n,1,3)+
  328           RCM(o,3,2)*RCM(n,2,3)+RCM(o,3,3)*RCM(n,3,3);
  329 }
  330 
  331 static void vector_multiply (float * v, float * m)
  332 {
  333   float o[4];
  334   guint i;
  335   
  336   for (i = 0; i < 4; i++) o[i] = v[i];
  337   
  338   v[0]=RC(0,0)*o[0]+RC(0,1)*o[1]+RC(0,2)*o[2]+RC(0,3)*o[3];
  339   v[1]=RC(1,0)*o[0]+RC(1,1)*o[1]+RC(1,2)*o[2]+RC(1,3)*o[3];
  340   v[2]=RC(2,0)*o[0]+RC(2,1)*o[1]+RC(2,2)*o[2]+RC(2,3)*o[3];
  341   v[3]=RC(3,0)*o[0]+RC(3,1)*o[1]+RC(3,2)*o[2]+RC(3,3)*o[3];
  342 }
  343 
  344 void gfs_gl_get_frustum (GfsGlViewParams * vp, GList * symmetries, GfsFrustum * f)
  345 {
  346   GLint v[4];
  347   float p[16];
  348   int i;
  349 
  350   f->res = 2.*vp->res;
  351   f->symmetries = symmetries;
  352   guint n = 1;
  353   while (symmetries) {
  354     n *= 2;
  355     symmetries = symmetries->next;
  356   }
  357   f->s = g_malloc (n*sizeof (FttVector));
  358 
  359   glGetIntegerv (GL_VIEWPORT, v);
  360   f->width = v[2];
  361   glGetFloatv (GL_MODELVIEW_MATRIX, f->m);
  362   glGetFloatv (GL_PROJECTION_MATRIX, f->p);
  363   for (i = 0; i < 16; i++) p[i] = f->p[i];
  364   matrix_multiply (p, f->m);
  365 
  366   /* right */
  367   f->n[0][0] = p[3] - p[0];
  368   f->n[0][1] = p[7] - p[4];
  369   f->n[0][2] = p[11] - p[8];
  370   f->d[0]    = p[15] - p[12];
  371    
  372   /* left */
  373   f->n[1][0] = p[3] + p[0];
  374   f->n[1][1] = p[7] + p[4];
  375   f->n[1][2] = p[11] + p[8];
  376   f->d[1]    = p[15] + p[12];
  377   
  378   /* top */
  379   f->n[2][0] = p[3] - p[1];
  380   f->n[2][1] = p[7] - p[5];
  381   f->n[2][2] = p[11] - p[9];
  382   f->d[2]    = p[15] - p[13];
  383 
  384   /* bottom */
  385   f->n[3][0] = p[3] + p[1];
  386   f->n[3][1] = p[7] + p[5];
  387   f->n[3][2] = p[11] + p[9];
  388   f->d[3]    = p[15] + p[13];
  389   
  390   /* front */
  391   f->n[4][0] = p[3] + p[2];
  392   f->n[4][1] = p[7] + p[6];
  393   f->n[4][2] = p[11] + p[10];
  394   f->d[4]    = p[15] + p[14];
  395   
  396   /* back */
  397   f->n[5][0] = p[3] - p[2];
  398   f->n[5][1] = p[7] - p[6];
  399   f->n[5][2] = p[11] - p[10];
  400   f->d[5]    = p[15] - p[14];
  401   
  402   for (i = 0; i < 6; i++) {
  403     gdouble n = gts_vector_norm (f->n[i]);
  404     if (n > 0.) {
  405       f->n[i][0] /= n; f->n[i][1] /= n; f->n[i][2] /= n;
  406       f->d[i] /= n;
  407     }
  408   }
  409 }
  410 
  411 void gfs_gl_frustum_free (GfsFrustum * f)
  412 {
  413   g_return_if_fail (f != NULL);
  414   g_free (f->s);
  415 }
  416 
  417 static guint create_symmetries (FttVector * f, GList * s, FttVector * p)
  418 {
  419   guint j = 0, n;
  420 
  421   f[j++] = *p; n = j;
  422   while (s) {
  423     guint i;
  424     for (i = 0; i < n; i++)
  425       gfs_gl_symmetry_transform (s->data, &f[i], &f[j++]);
  426     n = j;
  427     s = s->next;
  428   }
  429   return n;
  430 }
  431 
  432 /**
  433  * gfs_sphere_in_frustum:
  434  * @p: the sphere center.
  435  * @r: the sphere radius.
  436  * @f: the view frustum.
  437  * 
  438  * Returns: GTS_OUT if the sphere is outside the view frustum, GTS_IN
  439  * if it is inside, GTS_ON if it is partly inside.
  440  */
  441 GtsIntersect gfs_sphere_in_frustum (FttVector * p, gdouble r, GfsFrustum * f)
  442 {
  443   GtsIntersect I1 = GTS_OUT;
  444 
  445   g_return_val_if_fail (p != NULL, GTS_OUT);
  446   g_return_val_if_fail (f != NULL, GTS_OUT);
  447 
  448   guint j, n = create_symmetries (f->s, f->symmetries, p);
  449   for (j = 0; j < n; j++) {
  450     p = &f->s[j];
  451     guint i;
  452     GtsIntersect I = GTS_IN;
  453     for (i = 0; i < 6; i++) {
  454       gdouble d = f->n[i][0]*p->x + f->n[i][1]*p->y + f->n[i][2]*p->z + f->d[i];
  455       if (d < -r) {
  456     I = GTS_OUT;
  457     break;
  458       }
  459       if (d < r)
  460     I = GTS_ON;
  461     }
  462     if (I == GTS_IN)
  463       return GTS_IN;
  464     if (I == GTS_ON)
  465       I1 = GTS_ON;
  466   }
  467   return I1;
  468 }
  469 
  470 /**
  471  * gfs_sphere_is_small:
  472  * @c: the sphere center.
  473  * @r: the sphere radius.
  474  * @f: the view frustum.
  475  * 
  476  * Returns: %TRUE if the screen size (in pixels) of the projected
  477  * sphere is smaller than the resolution of @f, %FALSE otherwise.
  478  */
  479 gboolean gfs_sphere_is_small (FttVector * c, gdouble r, GfsFrustum * f)
  480 {
  481   g_return_val_if_fail (c != NULL, FALSE);
  482   g_return_val_if_fail (f != NULL, FALSE);
  483 
  484   guint j, n = create_symmetries (f->s, f->symmetries, c);
  485   for (j = 0; j < n; j++) {
  486     c = &f->s[j];
  487     float v[4];
  488     v[0] = c->x; v[1] = c->y; v[2] = c->z; v[3] = 1.;
  489     vector_multiply (v, f->m);
  490     v[0] = r;
  491     vector_multiply (v, f->p);
  492     float rp = v[3] == 0. ? 0 : v[0]*f->width/v[3];
  493     if (rp >= f->res)
  494       return FALSE;
  495   }
  496   return TRUE;
  497 }
  498 
  499 static void cell_traverse_visible_no_check (FttCell * root,
  500                         GfsFrustum * f,
  501                         gint maxlevel,
  502                         FttCellTraverseFunc func,
  503                         gpointer data)
  504 {
  505   if (FTT_CELL_IS_LEAF (root) || ftt_cell_level (root) == maxlevel)
  506     (* func) (root, data);
  507   else {
  508     gdouble r = ftt_cell_size (root)*GFS_DIAGONAL;
  509     FttVector p;
  510     
  511     ftt_cell_pos (root, &p);
  512     if (gfs_sphere_is_small (&p, r, f))
  513       (* func) (root, data);
  514     else {
  515       struct _FttOct * children = root->children;
  516       guint n;
  517       
  518       for (n = 0; n < FTT_CELLS; n++) {
  519     FttCell * c = &(children->cell[n]);
  520     if (!FTT_CELL_IS_DESTROYED (c))
  521       cell_traverse_visible_no_check (c, f, maxlevel, func, data);
  522       }
  523     }
  524   }
  525 }
  526 
  527 static void cell_traverse_visible (FttCell * root,
  528                    GfsFrustum * f,
  529                    gint maxlevel,
  530                    FttCellTraverseFunc func,
  531                    gpointer data)
  532 {
  533   gdouble r = ftt_cell_size (root)*GFS_DIAGONAL;
  534   FttVector p;
  535   GtsIntersect i;
  536 
  537   ftt_cell_pos (root, &p);
  538   i = gfs_sphere_in_frustum (&p, r, f);
  539   if (i == GTS_OUT)
  540     return;
  541   if (FTT_CELL_IS_LEAF (root) ||
  542       ftt_cell_level (root) == maxlevel || 
  543       gfs_sphere_is_small (&p, r, f))
  544     (* func) (root, data);
  545   else if (i == GTS_IN)
  546     cell_traverse_visible_no_check (root, f, maxlevel, func, data);
  547   else {
  548     struct _FttOct * children = root->children;
  549     guint n;
  550 
  551     for (n = 0; n < FTT_CELLS; n++) {
  552       FttCell * c = &(children->cell[n]);
  553       if (!FTT_CELL_IS_DESTROYED (c))
  554     cell_traverse_visible (c, f, maxlevel, func, data);
  555     }
  556   }
  557 }
  558 
  559 static void box_traverse_visible (GfsBox * b, gpointer * data)
  560 {
  561   cell_traverse_visible (b->root, data[0], *((gint *)data[3]), data[1], data[2]);
  562 }
  563 
  564 /**
  565  * gfs_gl_cell_traverse_visible:
  566  * @gl: a #GfsGl.
  567  * @f: a view frustum.
  568  * @func: a used-defined function.
  569  * @data: user data to pass to @func.
  570  *
  571  * Traverse the cells of @gl which are visible.
  572  */
  573 void gfs_gl_cell_traverse_visible (GfsGl * gl,
  574                    GfsFrustum * f,
  575                    FttCellTraverseFunc func,
  576                    gpointer data)
  577 {
  578   gpointer dat[4];
  579 
  580   g_return_if_fail (gl != NULL);
  581   g_return_if_fail (f != NULL);
  582   g_return_if_fail (func != NULL);
  583 
  584   dat[0] = f;
  585   dat[1] = func;
  586   dat[2] = data;
  587   dat[3] = &gl->maxlevel;
  588   gts_container_foreach (GTS_CONTAINER (gl->sim), (GtsFunc) box_traverse_visible, dat);
  589 }
  590 
  591 typedef struct {
  592   GfsFrustum * f;
  593   gint maxlevel;
  594   gboolean (* condition) (FttCell *, gpointer);
  595   gpointer datacon;
  596   FttCellTraverseFunc func;
  597   gpointer data;
  598 } ConditionParams;
  599 
  600 static void cell_traverse_visible_condition_no_check (FttCell * root, ConditionParams * par)
  601 {
  602   if (!(* par->condition) (root, par->datacon))
  603     return;
  604   if (FTT_CELL_IS_LEAF (root) || ftt_cell_level (root) == par->maxlevel)
  605     (* par->func) (root, par->data);
  606   else {
  607     gdouble r = ftt_cell_size (root)*SLIGHTLY_LARGER*GFS_DIAGONAL;
  608     FttVector p;
  609     
  610     ftt_cell_pos (root, &p);
  611     if (gfs_sphere_is_small (&p, r, par->f))
  612       (* par->func) (root, par->data);
  613     else {
  614       struct _FttOct * children = root->children;
  615       guint n;
  616       
  617       for (n = 0; n < FTT_CELLS; n++) {
  618     FttCell * c = &(children->cell[n]);
  619     if (!FTT_CELL_IS_DESTROYED (c))
  620       cell_traverse_visible_condition_no_check (c, par);
  621       }
  622     }
  623   }
  624 }
  625 
  626 static void cell_traverse_visible_condition (FttCell * root, ConditionParams * par)
  627 {
  628   gdouble r = ftt_cell_size (root)*SLIGHTLY_LARGER*GFS_DIAGONAL;
  629   FttVector p;
  630   GtsIntersect i;
  631 
  632   if (!(* par->condition) (root, par->datacon))
  633     return;
  634   ftt_cell_pos (root, &p);
  635   i = gfs_sphere_in_frustum (&p, r, par->f);
  636   if (i == GTS_OUT)
  637     return;
  638   if (FTT_CELL_IS_LEAF (root) ||
  639       ftt_cell_level (root) == par->maxlevel || 
  640       gfs_sphere_is_small (&p, r, par->f))
  641     (* par->func) (root, par->data);
  642   else if (i == GTS_IN)
  643     cell_traverse_visible_condition_no_check (root, par);
  644   else {
  645     struct _FttOct * children = root->children;
  646     guint n;
  647 
  648     for (n = 0; n < FTT_CELLS; n++) {
  649       FttCell * c = &(children->cell[n]);
  650       if (!FTT_CELL_IS_DESTROYED (c))
  651     cell_traverse_visible_condition (c, par);
  652     }
  653   }
  654 }
  655 
  656 static void box_traverse_visible_condition (GfsBox * b, ConditionParams * par)
  657 {
  658   cell_traverse_visible_condition (b->root, par);
  659 }
  660 
  661 /**
  662  * gfs_gl_cell_traverse_visible_condition:
  663  * @gl: a #GfsGl.
  664  * @f: a view frustum.
  665  * @condition: the condition to verify.
  666  * @datacon: user data to pass to @condition.
  667  * @func: a used-defined function.
  668  * @data: user data to pass to @func.
  669  *
  670  * Traverse the cells of @gl which are visible and verify @condition.
  671  */
  672 void gfs_gl_cell_traverse_visible_condition (GfsGl * gl,
  673                          GfsFrustum * f,
  674                          gboolean (* condition) (FttCell *, gpointer),
  675                          gpointer datacon,
  676                          FttCellTraverseFunc func,
  677                          gpointer data)
  678 {
  679   ConditionParams par;
  680 
  681   g_return_if_fail (gl != NULL);
  682   g_return_if_fail (f != NULL);
  683   g_return_if_fail (condition != NULL);
  684   g_return_if_fail (func != NULL);
  685 
  686   par.f = f;
  687   par.condition = condition;
  688   par.datacon = datacon;
  689   par.func = func;
  690   par.data = data;
  691   par.maxlevel = gl->maxlevel;
  692   gts_container_foreach (GTS_CONTAINER (gl->sim), (GtsFunc) box_traverse_visible_condition, &par);
  693 }
  694 
  695 static gboolean is_mixed (FttCell * cell, gpointer data)
  696 {
  697   return GFS_IS_MIXED (cell);
  698 }
  699 
  700 /**
  701  * gfs_gl_cell_traverse_visible_mixed:
  702  * @gl: a #GfsGl.
  703  * @f: a view frustum.
  704  * @func: a used-defined function.
  705  * @data: user data to pass to @func.
  706  *
  707  * Traverse the cells of @gl which are visible and mixed.
  708  */
  709 void gfs_gl_cell_traverse_visible_mixed (GfsGl * gl,
  710                      GfsFrustum * f,
  711                      FttCellTraverseFunc func,
  712                      gpointer data)
  713 {
  714   g_return_if_fail (gl != NULL);
  715   g_return_if_fail (f != NULL);
  716   g_return_if_fail (func != NULL);
  717 
  718   gfs_gl_cell_traverse_visible_condition (gl, f, is_mixed, NULL, func, data);
  719 }
  720 
  721 static void cell_traverse_visible_boundary_no_check (FttCell * root,
  722                              GfsFrustum * f,
  723                              FttDirection d,
  724                              gint maxlevel,
  725                              FttCellTraverseFunc func,
  726                              gpointer data)
  727 {
  728   if (FTT_CELL_IS_LEAF (root) || ftt_cell_level (root) == maxlevel)
  729     (* func) (root, data);
  730   else {
  731     gdouble r = ftt_cell_size (root)*GFS_DIAGONAL;
  732     FttVector p;
  733     
  734     ftt_cell_pos (root, &p);
  735     if (gfs_sphere_is_small (&p, r, f))
  736       (* func) (root, data);
  737     else {
  738       FttCellChildren child;
  739       guint n;
  740       
  741       ftt_cell_children_direction (root, d, &child);
  742       for (n = 0; n < FTT_CELLS/2; n++)
  743     if (child.c[n])
  744       cell_traverse_visible_boundary_no_check (child.c[n], f, d, maxlevel, func, data);
  745     }
  746   }
  747 }
  748 
  749 static void cell_traverse_visible_boundary (FttCell * root,
  750                         GfsFrustum * f,
  751                         FttDirection d,
  752                         gint maxlevel,
  753                         FttCellTraverseFunc func,
  754                         gpointer data)
  755 {
  756   gdouble r = ftt_cell_size (root)*GFS_DIAGONAL;
  757   FttVector p;
  758   GtsIntersect i;
  759 
  760   ftt_cell_pos (root, &p);
  761   i = gfs_sphere_in_frustum (&p, r, f);
  762   if (i == GTS_OUT)
  763     return;
  764   if (FTT_CELL_IS_LEAF (root) ||
  765       ftt_cell_level (root) == maxlevel ||
  766       gfs_sphere_is_small (&p, r, f))
  767     (* func) (root, data);
  768   else if (i == GTS_IN)
  769     cell_traverse_visible_boundary_no_check (root, f, d, maxlevel, func, data);
  770   else {
  771     FttCellChildren child;
  772     guint n;
  773 
  774     ftt_cell_children_direction (root, d, &child);
  775     for (n = 0; n < FTT_CELLS/2; n++)
  776       if (child.c[n])
  777     cell_traverse_visible_boundary (child.c[n], f, d, maxlevel, func, data);
  778   }
  779 }
  780 
  781 static void box_traverse_visible_boundary (GfsBox * b, gpointer * data)
  782 {
  783   FttDirection d;
  784 
  785   for (d = 0; d < FTT_NEIGHBORS; d++)
  786     if (!b->neighbor[d] || GFS_IS_BOUNDARY (b->neighbor[d]))
  787       cell_traverse_visible_boundary (b->root, data[0], d, *((gint *)data[3]), data[1], data[2]);
  788 }
  789 
  790 /**
  791  * gfs_gl_cell_traverse_visible_boundary:
  792  * @gl: a #GfsGl.
  793  * @f: a view frustum.
  794  * @func: a used-defined function.
  795  * @data: user data to pass to @func.
  796  *
  797  * Traverse the boundary cells of @gl which are visible.
  798  */
  799 void gfs_gl_cell_traverse_visible_boundary (GfsGl * gl,
  800                         GfsFrustum * f,
  801                         FttCellTraverseFunc func,
  802                         gpointer data)
  803 {
  804   gpointer dat[4];
  805 
  806   g_return_if_fail (gl != NULL);
  807   g_return_if_fail (f != NULL);
  808   g_return_if_fail (func != NULL);
  809 
  810   dat[0] = f;
  811   dat[1] = func;
  812   dat[2] = data;
  813   dat[3] = &gl->maxlevel;
  814   gts_container_foreach (GTS_CONTAINER (gl->sim), (GtsFunc) box_traverse_visible_boundary, dat);
  815 }
  816 
  817 void gfs_gl_init (void)
  818 {
  819   gfs_gl_label_class ();
  820   gfs_gl_cells_class ();
  821   gfs_gl_fractions_class ();
  822   gfs_gl_boundaries_class ();
  823   gfs_gl_squares_class ();
  824   gfs_gl_linear_class ();
  825   gfs_gl_isoline_class ();
  826   gfs_gl_solid_class ();
  827   gfs_gl_vof_class ();
  828   gfs_gl_solid_class ();
  829   gfs_gl_levels_class ();
  830   gfs_gl_vectors_class ();
  831     gfs_gl_streamlines_class ();
  832   gfs_gl_ellipses_class ();
  833   gfs_gl_location_class ();
  834     gfs_gl_height_class ();
  835   gfs_gl_locate_class ();
  836   gfs_gl_pipes_class ();
  837   gfs_gl_clip_plane_class ();
  838   gfs_gl_clip_plane_class ();
  839   gfs_gl_cut_plane_class ();
  840 #if (!FTT_2D)
  841   gfs_gl_isosurface_class ();
  842 #endif /* 3D */
  843   gfs_gl_symmetry_class ();
  844     gfs_gl_periodic_class ();
  845 }
  846 
  847 void gfs_gl_init_gl (void)
  848 {
  849   GLfloat light0_pos[4]   = { 0.0, 0.0, 50.0, 0.0 };
  850   GLfloat light0_color[4] = { 1., 1., 1., 1.0 }; /* white light */
  851 
  852   glDisable (GL_CULL_FACE);
  853   glEnable (GL_DEPTH_TEST);
  854   glEnable (GL_NORMALIZE);
  855 
  856   /* speedups */
  857   glEnable (GL_DITHER);
  858   glShadeModel (GL_SMOOTH);
  859   glHint (GL_PERSPECTIVE_CORRECTION_HINT, GL_FASTEST);
  860   glHint (GL_POLYGON_SMOOTH_HINT, GL_FASTEST);
  861 
  862   /* light */
  863   glLightfv (GL_LIGHT0, GL_POSITION, light0_pos);
  864   glLightfv (GL_LIGHT0, GL_DIFFUSE,  light0_color);
  865   glEnable (GL_LIGHT0);
  866   glEnable (GL_LIGHTING);
  867 
  868   glColorMaterial (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE);
  869   glEnable (GL_COLOR_MATERIAL);
  870 
  871   /* fonts */
  872 #if HAVE_FTGL
  873   if (default_font)
  874     ftglDestroyFont (default_font);
  875   if (default_raster_font)
  876     ftglDestroyFont (default_raster_font);
  877   default_font = default_raster_font = NULL;
  878 #endif /* HAVE_FTGL */
  879 }
  880 
  881 /* GfsGl: Object */
  882 
  883 static void gl_read (GtsObject ** o, GtsFile * fp)
  884 {
  885   GfsGl * gl = GFS_GL (*o);
  886   gchar * shading = NULL;
  887   GtsFileVariable var[] = {
  888     {GTS_FLOAT,  "r",           TRUE, &gl->lc.r},
  889     {GTS_FLOAT,  "g",           TRUE, &gl->lc.g},
  890     {GTS_FLOAT,  "b",           TRUE, &gl->lc.b},
  891     {GTS_STRING, "shading",     TRUE, &shading},
  892     {GTS_INT,    "maxlevel",    TRUE, &gl->maxlevel},
  893     {GTS_FLOAT,  "font_size",   TRUE, &gl->font_size},
  894     {GTS_INT,    "raster_font", TRUE, &gl->use_raster_font},
  895     {GTS_FLOAT,  "line_width",  TRUE, &gl->line_width},
  896     {GTS_NONE}
  897   };
  898 
  899   if (fp->type != GTS_STRING) {
  900     gts_file_error (fp, "expecting a string (class)");
  901     return;
  902   }
  903   gts_file_next_token (fp);
  904 
  905   gts_file_assign_variables (fp, var);
  906   if (fp->type == GTS_ERROR) {
  907     g_free (shading);
  908     return;
  909   }
  910 
  911   if (var[3].set) {
  912     if (!strcmp (shading, "Constant"))
  913       gl->shading = GFS_GL_CONSTANT;
  914     else if (!strcmp (shading, "Flat"))
  915       gl->shading = GFS_GL_FLAT;
  916     else if (!strcmp (shading, "Smooth"))
  917       gl->shading = GFS_GL_SMOOTH;
  918     else if (!strcmp (shading, "CSmooth"))
  919       gl->shading = GFS_GL_CSMOOTH;
  920     else {
  921       gts_file_variable_error (fp, var, "shading", "unknown shading `%s'", shading);
  922       g_free (shading);
  923       return;
  924     }
  925     g_free (shading);
  926   }
  927 
  928   gfs_gl_set_raster_font (gl, gl->use_raster_font);
  929 }
  930 
  931 static void gl_write (GtsObject * o, FILE * fp)
  932 {
  933   GfsGl * gl = GFS_GL (o);
  934 
  935   g_assert (strlen (o->klass->info.name) > 5);
  936   fprintf (fp, "%s {\n"
  937        "  r = %g g = %g b = %g\n"
  938        "  shading = %s\n"
  939        "  maxlevel = %d\n"
  940        "  font_size = %g\n"
  941        "  raster_font = %d\n"
  942        "  line_width = %g\n"
  943        "}",
  944        &o->klass->info.name[5],
  945        gl->lc.r, gl->lc.g, gl->lc.b,
  946        gl->shading == GFS_GL_CONSTANT ? "Constant" :
  947        gl->shading == GFS_GL_FLAT ?     "Flat" :
  948        gl->shading == GFS_GL_SMOOTH ?   "Smooth" :
  949        gl->shading == GFS_GL_CSMOOTH ?  "CSmooth" : 
  950                                         "Unknown",
  951        gl->maxlevel,
  952        gl->font_size,
  953        gl->use_raster_font,
  954        gl->line_width);
  955 }
  956 
  957 static void gl_set_simulation (GfsGl * gl, GfsSimulation * sim)
  958 {
  959   gl->sim = sim;
  960 }
  961 
  962 static gboolean gl_relevant (GfsSimulation * sim)
  963 {
  964   return TRUE;
  965 }
  966 
  967 static void gl_class_init (GfsGlClass * klass)
  968 {
  969   klass->set_simulation = gl_set_simulation;
  970   klass->relevant = gl_relevant;
  971   GTS_OBJECT_CLASS (klass)->read = gl_read;
  972   GTS_OBJECT_CLASS (klass)->write = gl_write;
  973 }
  974 
  975 static void gl_init (GfsGl * gl)
  976 {
  977   GtsColor c = { 0., 0., 0. };
  978 
  979   gl->shading = GFS_GL_CONSTANT;
  980   gl->lc = c;
  981   gl->maxlevel = -1;
  982   gl->format = GFSGL_SCREEN;
  983   gl->font_size = 1.;
  984   gl->use_raster_font = 1;
  985   gl->line_width = 1.;
  986 }
  987 
  988 GfsGlClass * gfs_gl_class (void)
  989 {
  990   static GfsGlClass * klass = NULL;
  991 
  992   if (klass == NULL) {
  993     GtsObjectClassInfo gfs_gl_info = {
  994       "GfsGl",
  995       sizeof (GfsGl),
  996       sizeof (GfsGlClass),
  997       (GtsObjectClassInitFunc) gl_class_init,
  998       (GtsObjectInitFunc) gl_init,
  999       (GtsArgSetFunc) NULL,
 1000       (GtsArgGetFunc) NULL
 1001     };
 1002     klass = gts_object_class_new (GTS_OBJECT_CLASS (gts_object_class ()),
 1003                   &gfs_gl_info);
 1004   }
 1005 
 1006   return klass;
 1007 }
 1008 
 1009 GfsGl * gfs_gl_new (GfsGlClass * klass)
 1010 {
 1011   GfsGl * object;
 1012 
 1013   g_return_val_if_fail (klass != NULL, NULL);
 1014 
 1015   object = GFS_GL (gts_object_new (GTS_OBJECT_CLASS (klass)));
 1016 
 1017   return object;
 1018 }
 1019 
 1020 /**
 1021  * gfs_gl_new_from_file:
 1022  * @fp: a #GtsFile pointer.
 1023  *
 1024  * Returns: a new #GfsGl object read from @fp.
 1025  */
 1026 GfsGl * gfs_gl_new_from_file (GtsFile * fp)
 1027 {
 1028   GtsObjectClass * klass;
 1029   GfsGl * gl;
 1030   GtsObject * o;
 1031 
 1032   g_return_val_if_fail (fp != NULL, NULL);
 1033 
 1034   klass = gts_object_class_from_name (fp->token->str);
 1035   if (klass == NULL) {
 1036     gchar * ename = g_strconcat ("GfsGl", fp->token->str, NULL);
 1037     klass = gts_object_class_from_name (ename);
 1038     g_free (ename);
 1039   }
 1040   if (klass == NULL || !gts_object_class_is_from_class (klass, gfs_gl_class ()))
 1041     return NULL;
 1042   gl = gfs_gl_new (GFS_GL_CLASS (klass));
 1043   o = GTS_OBJECT (gl);
 1044   (* klass->read) (&o, fp);
 1045   if (fp->type == GTS_ERROR) {
 1046     gts_object_destroy (o);
 1047     return NULL;
 1048   }
 1049   else
 1050     return gl;
 1051 }
 1052 
 1053 void gfs_gl_set_simulation (GfsGl * gl, GfsSimulation * sim)
 1054 {
 1055   g_return_if_fail (gl != NULL);
 1056   g_return_if_fail (sim != NULL);
 1057 
 1058   (* GFS_GL_CLASS (GTS_OBJECT (gl)->klass)->set_simulation) (gl, sim);
 1059 }
 1060 
 1061 void gfs_gl_draw (GfsGl * gl, GfsFrustum * f)
 1062 {
 1063   g_return_if_fail (gl != NULL);
 1064 
 1065   if (gl->sim && GFS_GL_CLASS (GTS_OBJECT (gl)->klass)->draw) {
 1066     glColor3f (gl->lc.r, gl->lc.g, gl->lc.b);
 1067     gl2psLineWidth (gl->line_width*gl->p->lw);
 1068     glLineWidth (gl->line_width);
 1069     (* GFS_GL_CLASS (GTS_OBJECT (gl)->klass)->draw) (gl, f);
 1070   }
 1071 }
 1072 
 1073 void gfs_gl_set_raster_font (GfsGl * gl, gboolean raster)
 1074 {
 1075   g_return_if_fail (gl != NULL);
 1076 
 1077   gl->use_raster_font = raster;
 1078 }
 1079 
 1080 void gfs_gl_set_font_size (GfsGl * gl, gfloat size)
 1081 {
 1082   g_return_if_fail (gl != NULL);
 1083 
 1084   gl->font_size = size;
 1085 }
 1086 
 1087 /* GfsGlLabel: Object */
 1088 
 1089 static void gl_label_destroy (GtsObject * o)
 1090 {
 1091   GfsGlLabel * l = GFS_GL_LABEL (o);
 1092 
 1093   if (l->formatted_label != l->label)
 1094     g_free (l->formatted_label);
 1095   g_free (l->label);
 1096 
 1097   (* GTS_OBJECT_CLASS (gfs_gl_label_class ())->parent_class->destroy) (o);
 1098 }
 1099 
 1100 static void gl_label_read (GtsObject ** o, GtsFile * fp)
 1101 {
 1102   (* GTS_OBJECT_CLASS (gfs_gl_label_class ())->parent_class->read) (o, fp);
 1103   if (fp->type == GTS_ERROR)
 1104     return;
 1105 
 1106   GfsGlLabel * gl = GFS_GL_LABEL (*o);
 1107   GtsFileVariable var[] = {
 1108     {GTS_DOUBLE, "x", TRUE, &gl->p.x },
 1109     {GTS_DOUBLE, "y", TRUE, &gl->p.y },
 1110     {GTS_DOUBLE, "z", TRUE, &gl->p.z },
 1111     {GTS_STRING, "label", TRUE, &gl->label },
 1112     {GTS_INT,  "symbol", TRUE, &gl->symbol },
 1113     {GTS_NONE}
 1114   };
 1115   gts_file_assign_variables (fp, var);
 1116 }
 1117 
 1118 static void gl_label_write (GtsObject * o, FILE * fp)
 1119 {
 1120   (* GTS_OBJECT_CLASS (gfs_gl_label_class ())->parent_class->write) (o, fp);
 1121   GfsGlLabel * gl = GFS_GL_LABEL (o);
 1122   fprintf (fp, 
 1123        " {\n"
 1124        "  x = %g y = %g z = %g\n"
 1125        "  label = \"%s\"\n"
 1126        "  symbol = %d\n"
 1127        "}",
 1128        gl->p.x, gl->p.y, gl->p.z,
 1129        gl->label, gl->symbol);
 1130 }
 1131 
 1132 static void gl_text_setup (GfsGl * gl)
 1133 {
 1134 #if HAVE_FTGL
 1135   if (!default_font) {
 1136     default_font = ftglCreatePolygonFont (DEFAULT_FONT);
 1137     default_raster_font = ftglCreateTextureFont (DEFAULT_FONT);
 1138     if (!default_font || !default_raster_font)
 1139       g_warning ("cannot create default FTGL font: %s", DEFAULT_FONT);
 1140     else {
 1141       ftglSetFontFaceSize (default_font,        DEFAULT_FONT_SIZE, DEFAULT_FONT_SIZE);
 1142       ftglSetFontFaceSize (default_raster_font, DEFAULT_FONT_SIZE, DEFAULT_FONT_SIZE);
 1143       /* fixme: for some weird reason using display lists doesn't work for polygon fonts */
 1144       ftglSetFontDisplayList (default_font, 0);
 1145     }
 1146   }
 1147 #endif /* HAVE_FTGL */
 1148 }
 1149 
 1150 static void gl_draw_text (GfsGl * gl, const gchar * text, 
 1151               gdouble x, gdouble y, gdouble z, 
 1152               gdouble size)
 1153 {
 1154 #if HAVE_FTGL
 1155   if (text) {
 1156     gl_text_setup (gl);
 1157     glMatrixMode (GL_MODELVIEW);
 1158     glPushMatrix ();
 1159     glNormal3f (0., 0., 1.);
 1160     glTranslatef (x, y, z);
 1161     size /= DEFAULT_FONT_SIZE;
 1162     glScalef (size, size, size);
 1163     ftglRenderFont (gfs_gl_vector_format (gl) || !gl->use_raster_font ? 
 1164             default_font : default_raster_font,
 1165             text,
 1166             FTGL_RENDER_ALL);
 1167     glPopMatrix ();
 1168   }
 1169 #endif /* HAVE_FTGL */
 1170 }
 1171 
 1172 static void gl_text_bbox (GfsGl * gl, const gchar * text, 
 1173               gdouble size,
 1174               float bbox[6])
 1175 {
 1176 #if HAVE_FTGL
 1177   if (text) {
 1178     gl_text_setup (gl);
 1179     ftglGetFontBBox (gfs_gl_vector_format (gl) || !gl->use_raster_font ? 
 1180              default_font : default_raster_font,
 1181              text, strlen (text), bbox);
 1182     int i;
 1183     for (i = 0; i < 6; i++)
 1184       bbox[i] *= size/DEFAULT_FONT_SIZE;
 1185   }
 1186 #endif /* HAVE_FTGL */
 1187 }
 1188 
 1189 static void gl_label_draw (GfsGl * gl, GfsFrustum * f)
 1190 {
 1191   GfsGlLabel * l = GFS_GL_LABEL (gl);
 1192   FttVector p = l->p;
 1193   gfs_simulation_map (gl->sim, &p);
 1194   gdouble size;
 1195   if (gl->maxlevel >= 0)
 1196     size = ftt_level_size (gl->maxlevel);
 1197   else {
 1198     FttCell * cell = gfs_domain_locate (GFS_DOMAIN (gl->sim), p, -1, NULL);
 1199     size = cell ? ftt_cell_size (cell) : 1./64.;
 1200   }
 1201 
 1202   gl->size = 0;
 1203   glMatrixMode (GL_PROJECTION);
 1204   glPushMatrix ();
 1205   glTranslatef (0., 0., gl->p->lc);
 1206   if (l->symbol) {
 1207     glMatrixMode (GL_MODELVIEW);
 1208     glPushMatrix ();
 1209     glTranslated (p.x, p.y, p.z);
 1210     glScaled (size, size, size);
 1211     draw_sphere ();
 1212     glPopMatrix ();
 1213     p.x += size/2.; p.y += size/2.;
 1214   }
 1215   glNormal3f (0., 0., 1.);
 1216   gl_draw_text (gl, l->formatted_label, p.x, p.y, p.z, gl->font_size*size);
 1217   glMatrixMode (GL_PROJECTION);
 1218   glPopMatrix ();
 1219 }
 1220 
 1221 void gfs_gl_label_set_label (GfsGlLabel * gl, const gchar * label, GfsSimulation * sim)
 1222 {
 1223   g_return_if_fail (gl != NULL);
 1224   g_return_if_fail (label != NULL);
 1225   g_return_if_fail (sim != NULL);
 1226 
 1227   if (label != gl->label) {
 1228     if (gl->formatted_label != gl->label)
 1229       g_free (gl->formatted_label);
 1230     gl->formatted_label = NULL;
 1231     g_free (gl->label);
 1232     gl->label = g_strdup (label);
 1233   }
 1234 
 1235   gboolean dynamic = FALSE, parallel = FALSE;
 1236   GSList * formats = gfs_format_new (gl->label, NULL, &dynamic, &parallel);
 1237   if (dynamic || parallel) {
 1238     if (gl->formatted_label != gl->label)
 1239       g_free (gl->formatted_label);
 1240     gl->formatted_label = 
 1241       gfs_format_string (formats, GFS_DOMAIN (sim)->pid, sim->time.i, sim->time.t);
 1242   }
 1243   else
 1244     gl->formatted_label = gl->label;
 1245   gfs_format_destroy (formats);  
 1246 }
 1247 
 1248 static void gl_label_set_simulation (GfsGl * gl, GfsSimulation * sim)
 1249 {  
 1250   (*GFS_GL_CLASS (GTS_OBJECT_CLASS (gfs_gl_label_class ())->parent_class)->set_simulation) 
 1251     (gl, sim);
 1252 
 1253   gfs_gl_label_set_label (GFS_GL_LABEL (gl), GFS_GL_LABEL (gl)->label, sim);
 1254 }
 1255 
 1256 static void gl_label_class_init (GfsGlClass * klass)
 1257 {
 1258   GTS_OBJECT_CLASS (klass)->destroy = gl_label_destroy;
 1259   GTS_OBJECT_CLASS (klass)->read = gl_label_read;
 1260   GTS_OBJECT_CLASS (klass)->write = gl_label_write;
 1261   klass->set_simulation = gl_label_set_simulation;
 1262   klass->draw = gl_label_draw;
 1263 }
 1264 
 1265 static void label_init (GfsGl * gl)
 1266 {
 1267   gl->maxlevel = 4;
 1268   gl->font_size = 3;
 1269   GFS_GL_LABEL (gl)->label = g_strdup ("Label");
 1270 }
 1271 
 1272 GfsGlClass * gfs_gl_label_class (void)
 1273 {
 1274   static GfsGlClass * klass = NULL;
 1275 
 1276   if (klass == NULL) {
 1277     GtsObjectClassInfo gfs_gl_label_info = {
 1278       "GfsGlLabel",
 1279       sizeof (GfsGl2D),
 1280       sizeof (GfsGl2DClass),
 1281       (GtsObjectClassInitFunc) gl_label_class_init,
 1282       (GtsObjectInitFunc) label_init,
 1283       (GtsArgSetFunc) NULL,
 1284       (GtsArgGetFunc) NULL
 1285     };
 1286     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl_class ()), &gfs_gl_label_info);
 1287   }
 1288 
 1289   return klass;
 1290 }
 1291 
 1292 /* GfsGl2D: Object */
 1293 
 1294 static void gl2D_read (GtsObject ** o, GtsFile * fp)
 1295 {
 1296   GfsGl2D * gl = GFS_GL2D (*o);
 1297   GtsFileVariable var[] = {
 1298     {GTS_DOUBLE, "n.x", TRUE},
 1299     {GTS_DOUBLE, "n.y", TRUE},
 1300     {GTS_DOUBLE, "n.z", TRUE},
 1301     {GTS_DOUBLE, "pos", TRUE},
 1302     {GTS_NONE}
 1303   };
 1304 
 1305   (* GTS_OBJECT_CLASS (gfs_gl2D_class ())->parent_class->read) (o, fp);
 1306   if (fp->type == GTS_ERROR)
 1307     return;
 1308 
 1309   var[0].data = &gl->n.x;
 1310   var[1].data = &gl->n.y;
 1311   var[2].data = &gl->n.z;
 1312   var[3].data = &gl->pos;
 1313   gts_file_assign_variables (fp, var);
 1314   if (fp->type == GTS_ERROR)
 1315     return;
 1316 
 1317   gfs_gl2D_update_plane (gl);
 1318 }
 1319 
 1320 static void gl2D_write (GtsObject * o, FILE * fp)
 1321 {
 1322   GfsGl2D * gl = GFS_GL2D (o);
 1323 
 1324   (* GTS_OBJECT_CLASS (gfs_gl2D_class ())->parent_class->write) (o, fp);
 1325 
 1326   fprintf (fp, " {\n"
 1327        "  n.x = %g n.y = %g n.z = %g\n"
 1328        "  pos = %g\n"
 1329        "}",
 1330        gl->n.x, gl->n.y, gl->n.z,
 1331        gl->pos);
 1332 }
 1333 
 1334 static void gl2D_update_plane (GfsGl2D * gl)
 1335 {
 1336   gdouble n;
 1337   GtsVector Q0 = {0., 0., 0.};
 1338   GtsVector Q1 = {0., 0., 0.};
 1339   gdouble max;
 1340   guint d = 0;
 1341   
 1342   n = gts_vector_norm (&gl->n.x);
 1343   g_assert (n > 0.);
 1344   gts_vector_normalize (&gl->n.x);
 1345 
 1346   /* build a vector orthogonal to the constraint */
 1347   max = gl->n.x*gl->n.x;
 1348   if (gl->n.y*gl->n.y > max) { max = gl->n.y*gl->n.y; d = 1; }
 1349   if (gl->n.z*gl->n.z > max) { max = gl->n.z*gl->n.z; d = 2; }
 1350   switch (d) {
 1351   case 0: Q0[0] = - gl->n.z/gl->n.x; Q0[2] = 1.0; break;
 1352   case 1: Q0[1] = - gl->n.z/gl->n.y; Q0[2] = 1.0; break;
 1353   case 2: Q0[2] = - gl->n.x/gl->n.z; Q0[0] = 1.0; break;
 1354   }
 1355   
 1356   /* build a second vector orthogonal to the first and to the constraint */
 1357   gts_vector_cross (Q1, &gl->n.x, Q0);
 1358   
 1359   gl->p[0].x = gl->pos*gl->n.x;
 1360   gl->p[0].y = gl->pos*gl->n.y;
 1361   gl->p[0].z = gl->pos*gl->n.z;
 1362   
 1363   gl->p[1].x = gl->p[0].x + Q0[0];
 1364   gl->p[1].y = gl->p[0].y + Q0[1];
 1365   gl->p[1].z = gl->p[0].z + Q0[2];
 1366   
 1367   gl->p[2].x = gl->p[0].x + Q1[0];
 1368   gl->p[2].y = gl->p[0].y + Q1[1];
 1369   gl->p[2].z = gl->p[0].z + Q1[2];
 1370 }
 1371 
 1372 static gdouble gl2D_pick (GfsGl * gl, GfsGlRay * r)
 1373 {
 1374   GfsGl2D * gl2 = GFS_GL2D (gl);
 1375   GtsVector AB;
 1376   gdouble ABn;
 1377 
 1378   gts_vector_init (AB, &r->a, &r->b);
 1379   ABn = gts_vector_scalar (AB, &gl2->n.x);
 1380   if (fabs (ABn) < 1e-6) {
 1381     gl2->picked = NULL;
 1382     return GFS_NODATA;
 1383   }
 1384   else {
 1385     GtsVector AP;
 1386     gdouble a;
 1387 
 1388     AP[0] = gl2->n.x*gl2->pos - r->a.x;
 1389     AP[1] = gl2->n.y*gl2->pos - r->a.y;
 1390     AP[2] = gl2->n.z*gl2->pos - r->a.z;
 1391     a = gts_vector_scalar (AP, &gl2->n.x)/ABn;
 1392     gl2->pickedpos.x = r->a.x*(1. - a) + a*r->b.x;
 1393     gl2->pickedpos.y = r->a.y*(1. - a) + a*r->b.y;
 1394     gl2->pickedpos.z = r->a.z*(1. - a) + a*r->b.z;
 1395     gl2->picked = gfs_domain_locate (GFS_DOMAIN (gl->sim), gl2->pickedpos, gl->maxlevel, NULL);
 1396     return gl2->picked ? a : GFS_NODATA;
 1397   }
 1398 }
 1399 
 1400 static void gl2D_class_init (GfsGlClass * klass)
 1401 {
 1402   klass->pick = gl2D_pick;
 1403   GTS_OBJECT_CLASS (klass)->read = gl2D_read;
 1404   GTS_OBJECT_CLASS (klass)->write = gl2D_write;
 1405   GFS_GL2D_CLASS (klass)->update_plane = gl2D_update_plane;
 1406 }
 1407 
 1408 static void gl2D_init (GfsGl2D * object)
 1409 {
 1410   object->n.x = 0.; object->n.y = 0.; object->n.z = 1.;
 1411   object->pos = 0.;
 1412 
 1413   object->p[0].x = object->p[0].y = object->p[0].z = 0.;
 1414   object->p[1].x = 1.; object->p[1].y = object->p[1].z = 0.;
 1415   object->p[2].y = 1.; object->p[2].x = object->p[2].z = 0.;
 1416 
 1417   gfs_gl2D_update_plane (object);
 1418 }
 1419 
 1420 GfsGl2DClass * gfs_gl2D_class (void)
 1421 {
 1422   static GfsGl2DClass * klass = NULL;
 1423 
 1424   if (klass == NULL) {
 1425     GtsObjectClassInfo gfs_gl2D_info = {
 1426       "GfsGl2D",
 1427       sizeof (GfsGl2D),
 1428       sizeof (GfsGl2DClass),
 1429       (GtsObjectClassInitFunc) gl2D_class_init,
 1430       (GtsObjectInitFunc) gl2D_init,
 1431       (GtsArgSetFunc) NULL,
 1432       (GtsArgGetFunc) NULL
 1433     };
 1434     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl_class ()), &gfs_gl2D_info);
 1435   }
 1436 
 1437   return klass;
 1438 }
 1439 
 1440 void gfs_gl2D_update_plane (GfsGl2D * gl)
 1441 {
 1442   g_return_if_fail (gl != NULL);
 1443   
 1444   (* GFS_GL2D_CLASS (GTS_OBJECT (gl)->klass)->update_plane) (gl);
 1445 }
 1446 
 1447 /* GfsGlSymmetry: Object */
 1448 
 1449 void gfs_gl_symmetry_transform (GfsGl * gl, FttVector * p, FttVector * t)
 1450 {
 1451   g_return_if_fail (gl != NULL);
 1452   g_return_if_fail (p != NULL);
 1453   g_return_if_fail (t != NULL);
 1454 
 1455   GLfloat in[4];
 1456   in[0] = p->x; in[1] = p->y; in[2] = p->z; in[3] = 1.;
 1457   vector_multiply (in, GFS_GL_SYMMETRY (gl)->m);
 1458   t->x = in[0]; t->y = in[1]; t->z = in[2];
 1459 }
 1460 
 1461 void gfs_gl_symmetry_apply (GList * symmetry, GLuint display_list)
 1462 {
 1463   if (symmetry) {
 1464     guint length = g_list_length (symmetry);
 1465     GLuint symmetry_list = glGenLists (length);
 1466     GLuint list = symmetry_list;
 1467     GList * i = symmetry;
 1468     while (i) {
 1469       glNewList (list, GL_COMPILE);
 1470       glCallList (display_list);
 1471       glMatrixMode (GL_MODELVIEW);
 1472       glPushMatrix ();
 1473       glMultMatrixf (GFS_GL_SYMMETRY (i->data)->m);
 1474       glCallList (display_list);
 1475       glPopMatrix ();
 1476       glEndList ();
 1477       display_list = list++;
 1478       i = i->next;
 1479     }
 1480     glCallList (display_list);
 1481     glDeleteLists (symmetry_list, length);
 1482   }
 1483   else
 1484     glCallList (display_list);
 1485 }
 1486 
 1487 static void gl_symmetry_update (GfsGl2D * gl)
 1488 {
 1489   GLfloat * s = GFS_GL_SYMMETRY (gl)->m, t[16];
 1490 
 1491   gts_vector_normalize (&gl->n.x);
 1492   if (GFS_GL_SYMMETRY (gl)->periodic) {
 1493     s[0] = 1.;                    s[1] = 0.;                      s[2] = 0.;
 1494     s[3] = 0.;
 1495     s[4] = s[1];                  s[5] = 1.; s[6] = 0.;       
 1496     s[7] = 0.;
 1497     s[8] = s[2];                  s[9] = s[6];                    s[10] = 1.; 
 1498     s[11] = 0.;
 1499     s[12] = 0.;                   s[13] = 0.;                     s[14] = 0.;                    
 1500     s[15] = 1.;
 1501 
 1502     t[0] = 1.;  t[1] = 0.;   t[2] = 0.;  t[3] = 0.;
 1503     t[4] = 0.;  t[5] = 1.;   t[6] = 0.;  t[7] = 0.;
 1504     t[8] = 0.;  t[9] = 0.;   t[10] = 1.; t[11] = 0.;
 1505     t[12] = gl->n.x*gl->pos; 
 1506     t[13] = gl->n.y*gl->pos;  
 1507     t[14] = gl->n.z*gl->pos; 
 1508     t[15] = 1.;
 1509   }
 1510   else {
 1511     s[0] = 1. - 2.*gl->n.x*gl->n.x; s[1] = - 2.*gl->n.x*gl->n.y;  s[2] = - 2.*gl->n.x*gl->n.z;
 1512     s[3] = 0.;
 1513     s[4] = s[1];                  s[5] = 1. - 2.*gl->n.y*gl->n.y; s[6] = - 2.*gl->n.y*gl->n.z;
 1514     s[7] = 0.;
 1515     s[8] = s[2];                  s[9] = s[6];                    s[10] = 1. - 2.*gl->n.z*gl->n.z; 
 1516     s[11] = 0.;
 1517     s[12] = 0.;                   s[13] = 0.;                     s[14] = 0.;                    
 1518     s[15] = 1.;
 1519 
 1520     t[0] = 1.;  t[1] = 0.;   t[2] = 0.;  t[3] = 0.;
 1521     t[4] = 0.;  t[5] = 1.;   t[6] = 0.;  t[7] = 0.;
 1522     t[8] = 0.;  t[9] = 0.;   t[10] = 1.; t[11] = 0.;
 1523     t[12] = - 2.*gl->n.x*gl->pos; 
 1524     t[13] = - 2.*gl->n.y*gl->pos;  
 1525     t[14] = - 2.*gl->n.z*gl->pos; 
 1526     t[15] = 1.;
 1527   }
 1528   matrix_multiply (s, t);
 1529 }
 1530 
 1531 static void gl_symmetry_class_init (GfsGl2DClass * klass)
 1532 {
 1533   GFS_GL_CLASS (klass)->pick = NULL;
 1534   klass->update_plane = gl_symmetry_update;
 1535 }
 1536 
 1537 GfsGlClass * gfs_gl_symmetry_class (void)
 1538 {
 1539   static GfsGlClass * klass = NULL;
 1540 
 1541   if (klass == NULL) {
 1542     GtsObjectClassInfo gfs_gl_symmetry_info = {
 1543       "GfsGlSymmetry",
 1544       sizeof (GfsGlSymmetry),
 1545       sizeof (GfsGl2DClass),
 1546       (GtsObjectClassInitFunc) gl_symmetry_class_init,
 1547       (GtsObjectInitFunc) NULL,
 1548       (GtsArgSetFunc) NULL,
 1549       (GtsArgGetFunc) NULL
 1550     };
 1551     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl2D_class ()), &gfs_gl_symmetry_info);
 1552   }
 1553 
 1554   return klass;
 1555 }
 1556 
 1557 /* GfsGlPeriodic: Object */
 1558 
 1559 static void gl_periodic_init (GfsGlSymmetry * gl)
 1560 {
 1561   gl->periodic = TRUE;
 1562 }
 1563 
 1564 GfsGlClass * gfs_gl_periodic_class (void)
 1565 {
 1566   static GfsGlClass * klass = NULL;
 1567 
 1568   if (klass == NULL) {
 1569     GtsObjectClassInfo gfs_gl_periodic_info = {
 1570       "GfsGlPeriodic",
 1571       sizeof (GfsGlSymmetry),
 1572       sizeof (GfsGl2DClass),
 1573       (GtsObjectClassInitFunc) NULL,
 1574       (GtsObjectInitFunc) gl_periodic_init,
 1575       (GtsArgSetFunc) NULL,
 1576       (GtsArgGetFunc) NULL
 1577     };
 1578     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl_symmetry_class ()), 
 1579                   &gfs_gl_periodic_info);
 1580   }
 1581 
 1582   return klass;
 1583 }
 1584 
 1585 /* GfsGlCells: Object */
 1586 
 1587 static void gl_cells_draw (GfsGl * gl, GfsFrustum * f)
 1588 {
 1589   gl->size = 0;
 1590   glMatrixMode (GL_PROJECTION);
 1591   glPushMatrix ();
 1592   glTranslatef (0., 0., gl->p->lc);
 1593   gfs_gl_normal (gl);
 1594   gfs_gl_cell_traverse_visible_plane (gl, f, (FttCellTraverseFunc) gl_cell, gl);
 1595   glPopMatrix ();
 1596 }
 1597 
 1598 static void gl_cells_class_init (GfsGlClass * klass)
 1599 {
 1600   klass->draw = gl_cells_draw;
 1601 }
 1602 
 1603 GfsGlClass * gfs_gl_cells_class (void)
 1604 {
 1605   static GfsGlClass * klass = NULL;
 1606 
 1607   if (klass == NULL) {
 1608     GtsObjectClassInfo gfs_gl_cells_info = {
 1609       "GfsGlCells",
 1610       sizeof (GfsGl2D),
 1611       sizeof (GfsGl2DClass),
 1612       (GtsObjectClassInitFunc) gl_cells_class_init,
 1613       (GtsObjectInitFunc) NULL,
 1614       (GtsArgSetFunc) NULL,
 1615       (GtsArgGetFunc) NULL
 1616     };
 1617     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl2D_class ()), &gfs_gl_cells_info);
 1618   }
 1619 
 1620   return klass;
 1621 }
 1622 
 1623 /* GfsGlFractions: Object */
 1624 
 1625 static gboolean gl_fractions_relevant (GfsSimulation * sim)
 1626 {
 1627   return (sim->solids->items != NULL);
 1628 }
 1629 
 1630 static void gl_fractions_class_init (GfsGlClass * klass)
 1631 {
 1632   klass->draw = gl_fractions_draw;
 1633   klass->relevant = gl_fractions_relevant;
 1634 }
 1635 
 1636 static void gl_fractions_init (GfsGl * gl)
 1637 {
 1638   GtsColor c = { 0., 0., 1.};
 1639 
 1640   gl->lc = c;
 1641 }
 1642 
 1643 GfsGlClass * gfs_gl_fractions_class (void)
 1644 {
 1645   static GfsGlClass * klass = NULL;
 1646 
 1647   if (klass == NULL) {
 1648     GtsObjectClassInfo gfs_gl_fractions_info = {
 1649       "GfsGlFractions",
 1650       sizeof (GfsGl),
 1651       sizeof (GfsGlClass),
 1652       (GtsObjectClassInitFunc) gl_fractions_class_init,
 1653       (GtsObjectInitFunc) gl_fractions_init,
 1654       (GtsArgSetFunc) NULL,
 1655       (GtsArgGetFunc) NULL
 1656     };
 1657     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl_class ()), &gfs_gl_fractions_info);
 1658   }
 1659 
 1660   return klass;
 1661 }
 1662 
 1663 /* GfsGlBoundaries: Object */
 1664 
 1665 static void gl_boundaries_draw (GfsGl * gl, GfsFrustum * f)
 1666 {
 1667   gl->size = 0;
 1668   glMatrixMode (GL_PROJECTION);
 1669   glPushMatrix ();
 1670   glTranslatef (0., 0., gl->p->lc);
 1671   glBegin (GL_LINES);
 1672   gfs_gl_cell_traverse_visible_boundary (gl, f, (FttCellTraverseFunc) gl_boundaries, gl);
 1673   glEnd ();
 1674   glPopMatrix ();
 1675 }
 1676 
 1677 static void gl_boundaries_class_init (GfsGlClass * klass)
 1678 {
 1679   klass->draw = gl_boundaries_draw;
 1680 }
 1681 
 1682 static void gl_boundaries_init (GfsGl * gl)
 1683 {
 1684   GtsColor c = { 0., 0., 0.};
 1685 
 1686   gl->lc = c;
 1687 }
 1688 
 1689 GfsGlClass * gfs_gl_boundaries_class (void)
 1690 {
 1691   static GfsGlClass * klass = NULL;
 1692 
 1693   if (klass == NULL) {
 1694     GtsObjectClassInfo gfs_gl_boundaries_info = {
 1695       "GfsGlBoundaries",
 1696       sizeof (GfsGl),
 1697       sizeof (GfsGlClass),
 1698       (GtsObjectClassInitFunc) gl_boundaries_class_init,
 1699       (GtsObjectInitFunc) gl_boundaries_init,
 1700       (GtsArgSetFunc) NULL,
 1701       (GtsArgGetFunc) NULL
 1702     };
 1703     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl_class ()), &gfs_gl_boundaries_info);
 1704   }
 1705 
 1706   return klass;
 1707 }
 1708 
 1709 /* GfsGlLevels: Object */
 1710 
 1711 static void gl_levels_destroy (GtsObject * object)
 1712 {
 1713   GfsGlLevels * gl = GFS_GL_LEVELS (object);
 1714 
 1715   if (gl->v)
 1716     gts_object_destroy (GTS_OBJECT (gl->v));
 1717 
 1718   (* GTS_OBJECT_CLASS (gfs_gl_levels_class ())->parent_class->destroy) (object);
 1719 }
 1720 
 1721 static void set_level (FttCell * cell, GfsVariable * v)
 1722 {
 1723   GFS_VALUE (cell, v) = ftt_cell_level (cell);
 1724 }
 1725 
 1726 static void gl_levels_set_simulation (GfsGl * object, GfsSimulation * sim)
 1727 {
 1728   GfsDomain * domain = GFS_DOMAIN (sim);
 1729   GfsGlLevels * gl = GFS_GL_LEVELS (object);
 1730 
 1731   (*GFS_GL_CLASS (GTS_OBJECT_CLASS (gfs_gl_levels_class ())->parent_class)->set_simulation)
 1732     (object, sim);
 1733 
 1734   if (gl->v)
 1735     gts_object_destroy (GTS_OBJECT (gl->v));
 1736   gl->v = gfs_temporary_variable (domain);
 1737 
 1738   gfs_domain_cell_traverse (domain,
 1739                 FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
 1740                 (FttCellTraverseFunc) set_level, gl->v);
 1741   gfs_domain_cell_traverse (domain,
 1742                 FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
 1743                 (FttCellTraverseFunc) gfs_get_from_below_intensive, gl->v);
 1744   gfs_domain_bc (domain, FTT_TRAVERSE_ALL, -1, gl->v);
 1745 }
 1746 
 1747 static void gl_levels_draw (GfsGl * gl, GfsFrustum * f)
 1748 {
 1749   gl->size = 0;
 1750   glMatrixMode (GL_PROJECTION);
 1751   glPushMatrix ();
 1752   glTranslatef (0., 0., gl->p->lc);
 1753   glBegin (GL_LINES);
 1754   gfs_gl_normal (gl);
 1755   gfs_gl_cell_traverse_visible_plane (gl, f, (FttCellTraverseFunc) gl_face, gl);
 1756   glEnd ();
 1757   glPopMatrix ();
 1758 }
 1759 
 1760 static void min_max_level (FttCell * cell, gint * l)
 1761 {
 1762   guint level = ftt_cell_level (cell);
 1763   if (level < l[0]) l[0] = level;
 1764   if (level > l[1]) l[1] = level;
 1765 }
 1766 
 1767 static gboolean gl_levels_relevant (GfsSimulation * sim)
 1768 {
 1769   gint l[2] = { G_MAXINT, 0 };
 1770 
 1771   gfs_domain_cell_traverse (GFS_DOMAIN (sim), FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
 1772                 (FttCellTraverseFunc) min_max_level, l);
 1773   return (l[1] > l[0]);
 1774 }
 1775 
 1776 static void gl_levels_class_init (GfsGlClass * klass)
 1777 {
 1778   klass->set_simulation = gl_levels_set_simulation;
 1779   klass->draw = gl_levels_draw;
 1780   klass->pick = NULL;
 1781   klass->relevant = gl_levels_relevant;
 1782   GTS_OBJECT_CLASS (klass)->destroy = gl_levels_destroy;
 1783 }
 1784 
 1785 GfsGlClass * gfs_gl_levels_class (void)
 1786 {
 1787   static GfsGlClass * klass = NULL;
 1788 
 1789   if (klass == NULL) {
 1790     GtsObjectClassInfo gfs_gl_levels_info = {
 1791       "GfsGlLevels",
 1792       sizeof (GfsGlLevels),
 1793       sizeof (GfsGl2DClass),
 1794       (GtsObjectClassInitFunc) gl_levels_class_init,
 1795       (GtsObjectInitFunc) NULL,
 1796       (GtsArgSetFunc) NULL,
 1797       (GtsArgGetFunc) NULL
 1798     };
 1799     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl2D_class ()),
 1800                   &gfs_gl_levels_info);
 1801   }
 1802 
 1803   return klass;
 1804 }
 1805 
 1806 /* GfsGlVarFunc: object */
 1807 
 1808 GfsGlVarFunc * gfs_gl_var_func_new (void)
 1809 {
 1810   GfsGlVarFunc * vf = g_malloc0 (sizeof (GfsGlVarFunc));
 1811   return vf;
 1812 }
 1813 
 1814 void gfs_gl_var_func_destroy (GfsGlVarFunc * vf)
 1815 {
 1816   g_return_if_fail (vf != NULL);
 1817 
 1818   if (vf->v && vf->v != gfs_function_get_variable (vf->f))
 1819     gts_object_destroy (GTS_OBJECT (vf->v));
 1820   if (vf->f)
 1821     gts_object_destroy (GTS_OBJECT (vf->f));
 1822   g_free (vf);
 1823 }
 1824 
 1825 static void update_v (FttCell * cell, gpointer * data)
 1826 {
 1827   GFS_VALUE (cell, GFS_VARIABLE (data[1])) = gfs_function_value (GFS_FUNCTION (data[0]), cell);
 1828 }
 1829 
 1830 GtsFile * gfs_gl_var_func_set (GfsGlVarFunc * vf, 
 1831                    GfsSimulation * sim, 
 1832                    const gchar * func,
 1833                    GString * expr,
 1834                    GfsVariableClass * klass)
 1835 {
 1836   GfsDomain * domain;
 1837   GfsFunction * f;
 1838   GtsFile * fp;
 1839 
 1840   g_return_val_if_fail (vf != NULL, NULL);
 1841   g_return_val_if_fail (sim != NULL, NULL);
 1842   g_return_val_if_fail (func != NULL, NULL);
 1843   if (!klass)
 1844     klass = gfs_variable_class ();
 1845   domain = GFS_DOMAIN (sim);
 1846 
 1847   fp = gts_file_new_from_string (func);
 1848   f = gfs_function_new (gfs_function_class (), 0.);
 1849   gfs_function_read (f, domain, fp);
 1850   gfs_pending_functions_compilation (fp);
 1851   if (fp->type == GTS_ERROR) {
 1852     gts_object_destroy (GTS_OBJECT (f));
 1853     return fp;
 1854   }
 1855   gts_file_destroy (fp);
 1856   
 1857   GfsVariable * v;
 1858   if (!(v = gfs_function_get_variable (f)) ||
 1859       !gts_object_class_is_from_class (GTS_OBJECT (v)->klass, klass) ||
 1860       (gfs_variable_is_dimensional (v) && sim->physical_params.L != 1.)) {
 1861     v = gfs_variable_new (klass, domain, NULL, NULL);
 1862     gfs_catch_floating_point_exceptions ();
 1863     gpointer data[2] = { f, v };
 1864     gfs_domain_cell_traverse (domain,
 1865                   FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
 1866                   (FttCellTraverseFunc) update_v, data);
 1867     if (gfs_restore_floating_point_exceptions ()) {
 1868       fp = gts_file_new_from_string (func);
 1869       gts_file_error (fp, "Floating-point exception");
 1870       gts_object_destroy (GTS_OBJECT (v));
 1871       gts_object_destroy (GTS_OBJECT (f));
 1872       return fp;
 1873     }
 1874     gfs_event_init (GFS_EVENT (v), sim);
 1875     gfs_domain_cell_traverse (domain,
 1876                   FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
 1877                   (FttCellTraverseFunc) v->fine_coarse, v);
 1878     gfs_domain_bc (domain, FTT_TRAVERSE_ALL, -1, v);
 1879   }
 1880 
 1881   if (vf->v && vf->v != gfs_function_get_variable (vf->f))
 1882     gts_object_destroy (GTS_OBJECT (vf->v));
 1883   if (vf->f)
 1884     gts_object_destroy (GTS_OBJECT (vf->f));
 1885   vf->f = f;
 1886   vf->v = v;
 1887 
 1888   if (expr && expr->str != func) {
 1889     g_free (expr->str);
 1890     expr->str = g_strdup (func);
 1891     expr->len = strlen (expr->str);
 1892   }
 1893   return NULL;
 1894 }
 1895 
 1896 /* GfsGlScalar: Object */
 1897 
 1898 static void gl_scalar_read (GtsObject ** o, GtsFile * fp)
 1899 {
 1900   GfsGlScalar * gl = GFS_GL_SCALAR (*o);
 1901   gchar * cmap;
 1902   GtsFileVariable var[] = {
 1903     {GTS_INT,    "amin", TRUE, &gl->amin},
 1904     {GTS_DOUBLE, "min",  TRUE, &gl->min},
 1905     {GTS_INT,    "amax", TRUE, &gl->amax},
 1906     {GTS_DOUBLE, "max",  TRUE, &gl->max},
 1907     {GTS_STRING, "cmap", TRUE, &cmap},
 1908     {GTS_INT,    "show", TRUE, &gl->show},
 1909     {GTS_NONE}
 1910   };
 1911 
 1912   (* GTS_OBJECT_CLASS (gfs_gl_scalar_class ())->parent_class->read) (o, fp);
 1913   if (fp->type == GTS_ERROR)
 1914     return;
 1915 
 1916   g_string_free (gl->expr, TRUE);
 1917   if (!(gl->expr = gfs_function_expression (fp, NULL)))
 1918     return;
 1919   gts_file_next_token (fp);
 1920 
 1921   gts_file_assign_variables (fp, var);
 1922   if (fp->type == GTS_ERROR)
 1923     return;
 1924 
 1925   if (var[4].set) {
 1926     gfs_colormap_destroy (gl->cmap);
 1927     if (!strcmp (cmap, "Cool"))
 1928       gl->cmap = gfs_colormap_cool_warm ();
 1929     else if (!strcmp (cmap, "Gray"))
 1930       gl->cmap = gfs_colormap_gray ();
 1931     else if (!strcmp (cmap, "Jet"))
 1932       gl->cmap = gfs_colormap_jet ();
 1933     else
 1934       gts_file_error (fp, "unknown colormap '%s'", cmap);
 1935     g_free (cmap);
 1936   }
 1937 }
 1938 
 1939 static void gl_scalar_write (GtsObject * o, FILE * fp)
 1940 {
 1941   GfsGlScalar * gl = GFS_GL_SCALAR (o);
 1942   
 1943   (* GTS_OBJECT_CLASS (gfs_gl_scalar_class ())->parent_class->write) (o, fp);
 1944 
 1945   fprintf (fp, " %s {\n  amin = %d", gl->expr->str, gl->amin);
 1946   if (gl->amin)
 1947     fputc ('\n', fp);
 1948   else
 1949     fprintf (fp, " min = %g\n", gl->min);
 1950   fprintf (fp, "  amax = %d", gl->amax);
 1951   if (gl->amax)
 1952     fputc ('\n', fp);
 1953   else
 1954     fprintf (fp, " max = %g\n", gl->max);
 1955   if (gl->show)
 1956     fputs ("  show = 1\n", fp);
 1957   fprintf (fp, "  cmap = %s\n}", gl->cmap->name);
 1958 }
 1959 
 1960 static void gl_scalar_destroy (GtsObject * object)
 1961 {
 1962   GfsGlScalar * gl = GFS_GL_SCALAR (object);
 1963 
 1964   gfs_gl_var_func_destroy (gl->vf);
 1965   g_string_free (gl->expr, TRUE);
 1966   if (gl->cmap)
 1967     gfs_colormap_destroy (gl->cmap);
 1968 
 1969   (* GTS_OBJECT_CLASS (gfs_gl_scalar_class ())->parent_class->destroy) (object);
 1970 }
 1971 
 1972 static void min_max (FttCell * cell, GfsGlScalar * gl)
 1973 {
 1974   gdouble v = GFS_VALUE (cell, gl->v);
 1975   if (v < G_MAXDOUBLE && v > gl->amaxv)
 1976     gl->amaxv = v;
 1977   if (v < gl->aminv)
 1978     gl->aminv = v;
 1979 }
 1980 
 1981 static GtsFile * gl_scalar_set (GfsGlScalar * gl, const gchar * func)
 1982 {
 1983   GtsFile * fp;
 1984 
 1985   if ((fp = gfs_gl_var_func_set (gl->vf, GFS_GL (gl)->sim, func, gl->expr, NULL)))
 1986     return fp;
 1987 
 1988   gl->v = gl->vf->v;
 1989   gl->amaxv = -G_MAXDOUBLE;
 1990   gl->aminv =  G_MAXDOUBLE;
 1991   gfs_domain_cell_traverse (GFS_DOMAIN (GFS_GL (gl)->sim),
 1992                 FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
 1993                 (FttCellTraverseFunc) min_max, gl);
 1994   gfs_all_reduce (GFS_DOMAIN (GFS_GL (gl)->sim), gl->amaxv, MPI_DOUBLE, MPI_MAX);
 1995   gfs_all_reduce (GFS_DOMAIN (GFS_GL (gl)->sim), gl->aminv, MPI_DOUBLE, MPI_MIN);
 1996   if (gl->amax) gl->max = gl->amaxv;
 1997   if (gl->amin) gl->min = gl->aminv;
 1998 
 1999   return NULL;
 2000 }
 2001 
 2002 static void gl_scalar_set_simulation (GfsGl * object, GfsSimulation * sim)
 2003 {
 2004   GfsGlScalar * gls = GFS_GL_SCALAR (object);
 2005   GtsFile * fp = NULL;
 2006 
 2007   (*GFS_GL_CLASS (GTS_OBJECT_CLASS (gfs_gl_scalar_class ())->parent_class)->set_simulation)
 2008     (object, sim);
 2009 
 2010   if (gls->expr->str[0] == '\0' || (fp = gfs_gl_scalar_set (gls, gls->expr->str))) {
 2011     GfsDomain * domain = GFS_DOMAIN (sim);
 2012     
 2013     if (domain->variables) {
 2014       GfsVariable * v = domain->variables->data;
 2015       gfs_gl_scalar_set (gls, v->name);
 2016     }
 2017     else
 2018       gfs_gl_scalar_set (gls, "0");
 2019   }
 2020   if (fp)
 2021     gts_file_destroy (fp);
 2022 }
 2023 
 2024 static gchar * gl_scalar_properties (GfsGlScalar * gl, int len)
 2025 {
 2026   gchar * s = g_strndup (gl->expr->str, len + 3);
 2027   if (strlen (s) > len) {
 2028     guint i;
 2029     for (i = 0; i < 3; i++)
 2030       s[len + i] = '.';
 2031   }
 2032   return s;
 2033 }
 2034 
 2035 static void gl_scalar_draw (GfsGl * gl, GfsFrustum * f)
 2036 {
 2037   GfsGlScalar * gls = GFS_GL_SCALAR (gl);
 2038 
 2039   /* draw colorscale */
 2040   if (gls->show && gls->max > gls->min) {
 2041     glMatrixMode (GL_PROJECTION);
 2042     glPushMatrix ();
 2043     GLint viewport[4];
 2044     glGetIntegerv(GL_VIEWPORT, viewport);
 2045     int width = viewport[2], height = viewport[3];
 2046     glLoadIdentity ();
 2047     gluOrtho2D (0, width, 0, height);
 2048     glMatrixMode (GL_MODELVIEW);
 2049     glPushMatrix ();
 2050     glLoadIdentity ();
 2051 
 2052     int i, np = 10;
 2053     float border = gl->font_size*height/20.;
 2054     float cwidth = border;
 2055     float cheight = (gl->font_size*height - 3.*border)/np;
 2056 
 2057     gchar s[80];
 2058     snprintf (s, 80, "% 6.3e", M_PI);
 2059     float bbox[6];
 2060     gl_text_bbox (gl, s, cheight/2., bbox);
 2061     float twidth  = bbox[3] - bbox[0];
 2062     float theight = bbox[4] - bbox[1];
 2063 
 2064     glTranslatef (width - border/2. - cwidth - border/4. - twidth, border, 0.);
 2065 
 2066     glColor3f (gl->lc.r, gl->lc.g, gl->lc.b);
 2067     gchar * title = gl_scalar_properties (gls, strlen (s) - 3.);
 2068     gl_text_bbox (gl, title, cheight/2., bbox);
 2069     float Twidth  = bbox[3] - bbox[0];
 2070     gl_draw_text (gl, title, 
 2071           cwidth + border/4. + (twidth - Twidth)/2., 
 2072           np*cheight + cheight/10., 0., cheight/2.);
 2073     g_free (title);
 2074     
 2075     for (i = 0; i < np; i++) {
 2076       GtsColor c = gfs_colormap_color (gls->cmap, (0.5 + i)/(float)np);
 2077       glColor3f (c.r, c.g, c.b);
 2078       glRectf (0, i*cheight, cwidth, (i + 1)*cheight);
 2079       snprintf (s, 80, "% 6.3e", gls->min + (0.5 + i)*(gls->max - gls->min)/(float)np);
 2080       glColor3f (gl->lc.r, gl->lc.g, gl->lc.b);
 2081       gl_draw_text (gl, s, cwidth + border/4., i*cheight + (cheight - theight)/2., 0., cheight/2.);
 2082     }
 2083 
 2084     glPopMatrix ();
 2085     glMatrixMode (GL_PROJECTION);
 2086     glPopMatrix ();
 2087   }
 2088 }
 2089 
 2090 static void gl_scalar_class_init (GfsGlClass * klass)
 2091 {
 2092   klass->set_simulation = gl_scalar_set_simulation;
 2093   klass->draw = gl_scalar_draw;
 2094   GTS_OBJECT_CLASS (klass)->read = gl_scalar_read;
 2095   GTS_OBJECT_CLASS (klass)->write = gl_scalar_write;
 2096   GTS_OBJECT_CLASS (klass)->destroy = gl_scalar_destroy;
 2097   GFS_GL_SCALAR_CLASS (klass)->set_scalar = gl_scalar_set;
 2098 }
 2099 
 2100 static void gl_scalar_init (GfsGlScalar * object)
 2101 {
 2102   object->expr = g_string_new ("");
 2103   object->vf = gfs_gl_var_func_new ();
 2104   object->amax = object->amin = TRUE;
 2105   object->min = 0.;
 2106   object->max = 1.;
 2107   object->cmap = gfs_colormap_jet ();
 2108 }
 2109 
 2110 GfsGlScalarClass * gfs_gl_scalar_class (void)
 2111 {
 2112   static GfsGlScalarClass * klass = NULL;
 2113 
 2114   if (klass == NULL) {
 2115     GtsObjectClassInfo gfs_gl_scalar_info = {
 2116       "GfsGlScalar",
 2117       sizeof (GfsGlScalar),
 2118       sizeof (GfsGlScalarClass),
 2119       (GtsObjectClassInitFunc) gl_scalar_class_init,
 2120       (GtsObjectInitFunc) gl_scalar_init,
 2121       (GtsArgSetFunc) NULL,
 2122       (GtsArgGetFunc) NULL
 2123     };
 2124     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl2D_class ()),
 2125                   &gfs_gl_scalar_info);
 2126   }
 2127 
 2128   return klass;
 2129 }
 2130 
 2131 GtsFile * gfs_gl_scalar_set (GfsGlScalar * gl, const gchar * func)
 2132 {
 2133   g_return_val_if_fail (gl != NULL, NULL);
 2134   g_return_val_if_fail (func != NULL, NULL);
 2135 
 2136   return (* GFS_GL_SCALAR_CLASS (GTS_OBJECT (gl)->klass)->set_scalar) (gl, func);
 2137 }
 2138 
 2139 /* GfsGlSquares: Object */
 2140 
 2141 static void gl_squares_class_init (GfsGlClass * klass)
 2142 {
 2143   klass->draw = gl_squares_draw;
 2144 }
 2145 
 2146 GfsGlClass * gfs_gl_squares_class (void)
 2147 {
 2148   static GfsGlClass * klass = NULL;
 2149 
 2150   if (klass == NULL) {
 2151     GtsObjectClassInfo gfs_gl_squares_info = {
 2152       "GfsGlSquares",
 2153       sizeof (GfsGlSquares),
 2154       sizeof (GfsGlScalarClass),
 2155       (GtsObjectClassInitFunc) gl_squares_class_init,
 2156       (GtsObjectInitFunc) NULL,
 2157       (GtsArgSetFunc) NULL,
 2158       (GtsArgGetFunc) NULL
 2159     };
 2160     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl_scalar_class ()),
 2161                   &gfs_gl_squares_info);
 2162   }
 2163 
 2164   return klass;
 2165 }
 2166 
 2167 /* GfsGlLinear: Object */
 2168 
 2169 static void gl_linear_destroy (GtsObject * o)
 2170 {
 2171   GfsGlLinear * gl = GFS_GL_LINEAR (o);
 2172 
 2173   gfs_gl_var_func_destroy (gl->vf);
 2174   g_string_free (gl->expr, TRUE);
 2175   if (gl->nx) {
 2176     gts_object_destroy (GTS_OBJECT (gl->nx));
 2177     gts_object_destroy (GTS_OBJECT (gl->ny));
 2178   }
 2179 
 2180   (* GTS_OBJECT_CLASS (gfs_gl_linear_class ())->parent_class->destroy) (o);
 2181 }
 2182 
 2183 static void gl_linear_read (GtsObject ** o, GtsFile * fp)
 2184 {
 2185   GfsGlLinear * gl = GFS_GL_LINEAR (*o);
 2186 
 2187   (* GTS_OBJECT_CLASS (gfs_gl_linear_class ())->parent_class->read) (o, fp);
 2188   if (fp->type == GTS_ERROR)
 2189     return;
 2190 
 2191   if (fp->type != '{') { /* for backward compatibility */
 2192     GtsFileVariable var[] = {
 2193       {GTS_INT,    "reversed",   TRUE},
 2194       {GTS_INT,    "use_scalar", TRUE},
 2195       {GTS_NONE}
 2196     };
 2197 
 2198     g_string_free (gl->expr, TRUE);
 2199     if (!(gl->expr = gfs_function_expression (fp, NULL)))
 2200       return;
 2201     gts_file_next_token (fp);
 2202     
 2203     var[0].data = &gl->reversed;
 2204     var[1].data = &gl->use_scalar;
 2205     gts_file_assign_variables (fp, var);
 2206     if (fp->type == GTS_ERROR)
 2207       return;
 2208   }
 2209   else { /* this is the old format for linear */
 2210     g_warning ("obsolete GfsGlLinear/GfsGlIsoline syntax");
 2211     if (!GFS_IS_GL_ISOLINE (gl)) {
 2212       gdouble scale;
 2213       GtsFileVariable var[] = {
 2214     {GTS_DOUBLE, "scale", TRUE},
 2215     {GTS_NONE}
 2216       };
 2217       var[0].data = &scale;
 2218       gts_file_assign_variables (fp, var);
 2219       if (fp->type == GTS_ERROR)
 2220     return;
 2221     }
 2222     gl->use_scalar = (GfsVariable *) TRUE;
 2223   }
 2224 }
 2225 
 2226 static void gl_linear_write (GtsObject * o, FILE * fp)
 2227 {
 2228   GfsGlLinear * gl = GFS_GL_LINEAR (o);
 2229   
 2230   (* GTS_OBJECT_CLASS (gfs_gl_linear_class ())->parent_class->write) (o, fp);
 2231   fprintf (fp, " %s {\n"
 2232        "  reversed = %d\n"
 2233        "  use_scalar = %d\n"
 2234        "}",
 2235        gl->expr->str,
 2236        gl->reversed, 
 2237        (gl->use_scalar != NULL));
 2238 }
 2239 
 2240 static void gl_linear_set_simulation (GfsGl * object, GfsSimulation * sim)
 2241 {
 2242   GfsGlLinear * gls = GFS_GL_LINEAR (object);
 2243   GtsFile * fp = NULL;
 2244   
 2245   if (GFS_GL_CLASS (GTS_OBJECT_CLASS (gfs_gl_linear_class ())->parent_class)->set_simulation)
 2246     (*GFS_GL_CLASS (GTS_OBJECT_CLASS (gfs_gl_linear_class ())->parent_class)->set_simulation)
 2247       (object, sim);
 2248 
 2249   if (gls->nx) {
 2250     gts_object_destroy (GTS_OBJECT (gls->nx));
 2251     gts_object_destroy (GTS_OBJECT (gls->ny));
 2252     gls->nx = gls->ny = NULL;
 2253   }
 2254 
 2255   if (gls->expr->str[0] == '\0' || (fp = gfs_gl_linear_set (gls, gls->expr->str)))
 2256     gfs_gl_linear_set (gls, "0");
 2257   if (fp)
 2258     gts_file_destroy (fp);
 2259 }
 2260 
 2261 static void gl_linear_class_init (GfsGlClass * klass)
 2262 {
 2263   GTS_OBJECT_CLASS (klass)->destroy = gl_linear_destroy;
 2264   GTS_OBJECT_CLASS (klass)->read  = gl_linear_read;
 2265   GTS_OBJECT_CLASS (klass)->write = gl_linear_write;
 2266   klass->set_simulation = gl_linear_set_simulation;
 2267   klass->draw = gl_linear_draw;
 2268 }
 2269 
 2270 static void gl_linear_init (GfsGlLinear * gl)
 2271 {
 2272   GtsColor c = { 1., 1., 1. };
 2273 
 2274   gl->expr = g_string_new ("");
 2275   gl->vf = gfs_gl_var_func_new ();
 2276   GFS_GL (gl)->lc = c;
 2277   gl->use_scalar = (GfsVariable *) TRUE;
 2278 }
 2279 
 2280 GfsGlClass * gfs_gl_linear_class (void)
 2281 {
 2282   static GfsGlClass * klass = NULL;
 2283 
 2284   if (klass == NULL) {
 2285     GtsObjectClassInfo gfs_gl_linear_info = {
 2286       "GfsGlLinear",
 2287       sizeof (GfsGlLinear),
 2288       sizeof (GfsGlScalarClass),
 2289       (GtsObjectClassInitFunc) gl_linear_class_init,
 2290       (GtsObjectInitFunc) gl_linear_init,
 2291       (GtsArgSetFunc) NULL,
 2292       (GtsArgGetFunc) NULL
 2293     };
 2294     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl_scalar_class ()),
 2295                   &gfs_gl_linear_info);
 2296   }
 2297 
 2298   return klass;
 2299 }
 2300 
 2301 static void normals (FttCell * cell, GfsGlLinear * gl)
 2302 {
 2303   gdouble size = ftt_cell_size (cell);
 2304   GFS_VALUE (cell, gl->nx) = gfs_center_gradient (cell, FTT_X, gl->vf->v->i)/size;
 2305   GFS_VALUE (cell, gl->ny) = gfs_center_gradient (cell, FTT_Y, gl->vf->v->i)/size;
 2306 }
 2307 
 2308 GtsFile * gfs_gl_linear_set (GfsGlLinear * gl, const gchar * func)
 2309 {
 2310   GtsFile * fp;
 2311 
 2312   g_return_val_if_fail (gl != NULL, NULL);
 2313   g_return_val_if_fail (func != NULL, NULL);
 2314 
 2315   if ((fp = gfs_gl_var_func_set (gl->vf, GFS_GL (gl)->sim, func, gl->expr, NULL)))
 2316     return fp;
 2317 
 2318   if (gfs_function_get_constant_value (gl->vf->f) != 0.) {
 2319     GfsDomain * domain = GFS_DOMAIN (GFS_GL (gl)->sim);
 2320     if (!gl->nx) {
 2321       gl->nx = gfs_temporary_variable (domain);
 2322       gl->nx->component = FTT_X;
 2323       gl->ny = gfs_temporary_variable (domain);
 2324       gl->nx->component = FTT_Y;
 2325     }
 2326     gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
 2327                   (FttCellTraverseFunc) normals, gl);
 2328     gfs_domain_bc (domain, FTT_TRAVERSE_ALL, -1, gl->nx);
 2329     gfs_domain_bc (domain, FTT_TRAVERSE_ALL, -1, gl->ny);
 2330   }
 2331 
 2332   return NULL;
 2333 }
 2334 
 2335 /* GfsGlIsoline: Object */
 2336 
 2337 static gboolean iso_cuts_cell (FttCell * cell, gpointer data)
 2338 {
 2339   IsoParams * par = data;
 2340   return (par->level >= GFS_VALUE (cell, par->min) &&
 2341       par->level <= GFS_VALUE (cell, par->max));
 2342 }
 2343 
 2344 /**
 2345  * gfs_gl_cell_traverse_visible_iso:
 2346  * @gl: a #GfsGl2D.
 2347  * @f: a view frustum.
 2348  * @min: a variable storing the minimum level.
 2349  * @max: a variable storing the maximum level.
 2350  * @level: the level of the isoline/surface.
 2351  * @func: a used-defined function.
 2352  * @data: user data to pass to @func.
 2353  *
 2354  * Traverse the cells of @gl which are visible and whose bounding
 2355  * sphere is intersected by the isoline/surface defined by @min, @max
 2356  * and @level.
 2357  */
 2358 void gfs_gl_cell_traverse_visible_iso (GfsGl * gl,
 2359                        GfsFrustum * f,
 2360                        GfsVariable * min,
 2361                        GfsVariable * max,
 2362                        gdouble level,
 2363                        FttCellTraverseFunc func,
 2364                        gpointer data)
 2365 {
 2366   IsoParams p = { min, max, level };
 2367 
 2368   g_return_if_fail (gl != NULL);
 2369   g_return_if_fail (f != NULL);
 2370   g_return_if_fail (min != NULL);
 2371   g_return_if_fail (max != NULL);
 2372   g_return_if_fail (func != NULL);
 2373 
 2374   gfs_gl_cell_traverse_visible_condition (gl, f, iso_cuts_cell, &p, func, data);
 2375 }
 2376 
 2377 static void gl_isoline_destroy (GtsObject * o)
 2378 {
 2379   GfsGlIsoline * gl = GFS_GL_ISOLINE (o);
 2380 
 2381   if (gl->min)
 2382     gts_object_destroy (GTS_OBJECT (gl->min));
 2383   if (gl->max)
 2384     gts_object_destroy (GTS_OBJECT (gl->max));
 2385   g_array_free (gl->levels, TRUE);
 2386   g_free (gl->ls);
 2387 
 2388   (* GTS_OBJECT_CLASS (gfs_gl_isoline_class ())->parent_class->destroy) (o);
 2389 }
 2390 
 2391 static void gl_isoline_update_levels (GfsGl * gl)
 2392 {
 2393   GfsGlScalar * gls = GFS_GL_SCALAR (gl);
 2394   GfsGlIsoline * gli = GFS_GL_ISOLINE (gl);
 2395 
 2396   g_array_set_size (gli->levels, 0);
 2397   if (gli->n > 0.) {
 2398     guint i;
 2399     for (i = 0; i < gli->n; i++) {
 2400       gdouble l = gls->min + (i + 1)*(gls->max - gls->min)/(gli->n + 1.);
 2401       g_array_append_val (gli->levels, l);
 2402     }
 2403   }
 2404 
 2405   if (gli->ls) {
 2406     gchar ** list, ** s;
 2407 
 2408     s = list = g_strsplit (gli->ls, ",", -1);
 2409     while (*s) {
 2410       char * end;
 2411       gdouble l = strtod (*s, &end);
 2412       if (*end == '\0')
 2413     g_array_append_val (gli->levels, l);
 2414       s++;
 2415     }
 2416     g_strfreev (list);
 2417   }
 2418 }
 2419 
 2420 static void gl_isoline_read (GtsObject ** o, GtsFile * fp)
 2421 {
 2422   GtsFileVariable var[] = {
 2423     {GTS_DOUBLE, "n",      TRUE},
 2424     {GTS_STRING, "levels", TRUE},
 2425     {GTS_NONE}
 2426   };
 2427 
 2428   (* GTS_OBJECT_CLASS (gfs_gl_isoline_class ())->parent_class->read) (o, fp);
 2429   if (fp->type == GTS_ERROR)
 2430     return;
 2431 
 2432   var[0].data = &GFS_GL_ISOLINE (*o)->n;
 2433   var[1].data = &GFS_GL_ISOLINE (*o)->ls;
 2434   gts_file_assign_variables (fp, var);
 2435   if (fp->type == GTS_ERROR)
 2436     return;
 2437 
 2438   gl_isoline_update_levels (GFS_GL (*o));
 2439 }
 2440 
 2441 static void gl_isoline_write (GtsObject * o, FILE * fp)
 2442 {
 2443   const gchar * levels = GFS_GL_ISOLINE (o)->ls;
 2444 
 2445   (* GTS_OBJECT_CLASS (gfs_gl_isoline_class ())->parent_class->write) (o, fp);
 2446 
 2447   fprintf (fp, " {\n  n = %g", GFS_GL_ISOLINE (o)->n);
 2448   if (levels && levels[0] != '\0')
 2449     fprintf (fp, " levels = %s", levels);
 2450   fprintf (fp, "\n}");
 2451 }
 2452 
 2453 static void min_max_iso (FttCell * cell, GfsGlIsoline * gl)
 2454 {
 2455   gdouble min = G_MAXDOUBLE, max = -G_MAXDOUBLE;
 2456   guint i;
 2457 
 2458   if (FTT_CELL_IS_LEAF (cell)) {
 2459     for (i = 0; i < NCORNERS; i++) {
 2460       gdouble v = gfs_cell_corner_value (cell, corner[i], GFS_GL_SCALAR (gl)->v, -1);
 2461       if (v < min) min = v;
 2462       if (v > max) max = v;
 2463     }
 2464   }
 2465   else {
 2466     FttCellChildren child;
 2467 
 2468     ftt_cell_children (cell, &child);
 2469     for (i = 0; i < FTT_CELLS; i++)
 2470       if (child.c[i]) {
 2471     gdouble vmin = GFS_VALUE (child.c[i], gl->min);
 2472     gdouble vmax = GFS_VALUE (child.c[i], gl->max);
 2473     if (vmin < min) min = vmin;
 2474     if (vmax > max) max = vmax;
 2475       }
 2476   }
 2477   GFS_VALUE (cell, gl->min) = min;
 2478   GFS_VALUE (cell, gl->max) = max;
 2479 }
 2480 
 2481 static GtsFile * gl_isoline_set_scalar (GfsGlScalar * gl, const gchar * func)
 2482 {
 2483   GtsFile * fp = (*GFS_GL_SCALAR_CLASS 
 2484           (GTS_OBJECT_CLASS (gfs_gl_isoline_class ())->parent_class)->set_scalar)
 2485     (gl, func);
 2486   if (fp != NULL)
 2487     return fp;
 2488   
 2489   gfs_domain_cell_traverse (GFS_DOMAIN (GFS_GL (gl)->sim),
 2490                 FTT_POST_ORDER, FTT_TRAVERSE_ALL, -1,
 2491                 (FttCellTraverseFunc) min_max_iso, gl);  
 2492   return fp;
 2493 }
 2494 
 2495 static void gl_isoline_set_simulation (GfsGl * object, GfsSimulation * sim)
 2496 {
 2497   GfsDomain * domain = GFS_DOMAIN (sim);
 2498   GfsGlIsoline * gls = GFS_GL_ISOLINE (object);
 2499 
 2500   if (gls->min)
 2501     gts_object_destroy (GTS_OBJECT (gls->min));
 2502   gls->min = gfs_temporary_variable (domain);
 2503   if (gls->max)
 2504     gts_object_destroy (GTS_OBJECT (gls->max));
 2505   gls->max = gfs_temporary_variable (domain);
 2506 
 2507   (*GFS_GL_CLASS (GTS_OBJECT_CLASS (gfs_gl_isoline_class ())->parent_class)->set_simulation)
 2508     (object, sim);
 2509 }
 2510 
 2511 static void gl_isoline_class_init (GfsGlClass * klass)
 2512 {
 2513   GTS_OBJECT_CLASS (klass)->destroy = gl_isoline_destroy;
 2514   GTS_OBJECT_CLASS (klass)->read = gl_isoline_read;
 2515   GTS_OBJECT_CLASS (klass)->write = gl_isoline_write;
 2516   klass->set_simulation = gl_isoline_set_simulation;
 2517   GFS_GL_SCALAR_CLASS (klass)->set_scalar = gl_isoline_set_scalar;
 2518   klass->draw = gl_isoline_draw;
 2519   klass->pick = NULL;
 2520 }
 2521 
 2522 static void gl_isoline_init (GfsGlIsoline * gl)
 2523 {
 2524   GtsColor c = { 0., 0., 0. };
 2525   GFS_GL (gl)->lc = c;
 2526   gl->levels = g_array_new (FALSE, FALSE, sizeof (gdouble));
 2527   gl->n = 10.;
 2528 }
 2529 
 2530 GfsGlClass * gfs_gl_isoline_class (void)
 2531 {
 2532   static GfsGlClass * klass = NULL;
 2533 
 2534   if (klass == NULL) {
 2535     GtsObjectClassInfo gfs_gl_isoline_info = {
 2536       "GfsGlIsoline",
 2537       sizeof (GfsGlIsoline),
 2538       sizeof (GfsGlScalarClass),
 2539       (GtsObjectClassInitFunc) gl_isoline_class_init,
 2540       (GtsObjectInitFunc) gl_isoline_init,
 2541       (GtsArgSetFunc) NULL,
 2542       (GtsArgGetFunc) NULL
 2543     };
 2544     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl_linear_class ()),
 2545                   &gfs_gl_isoline_info);
 2546   }
 2547 
 2548   return klass;
 2549 }
 2550 
 2551 void gfs_gl_isoline_set_levels (GfsGlIsoline * gl, const gchar * levels)
 2552 {
 2553   g_return_if_fail (gl != NULL);
 2554   g_return_if_fail (levels != NULL);
 2555 
 2556   g_free (gl->ls);
 2557   gl->ls = g_strdup (levels);
 2558 }
 2559 
 2560 /* GfsGlVOF: Object */
 2561 
 2562 static void gl_vof_destroy (GtsObject * o)
 2563 {
 2564   GfsGlVOF * gl = GFS_GL_VOF (o);
 2565 
 2566   gfs_gl_var_func_destroy (gl->vf);
 2567   g_string_free (gl->expr, TRUE);
 2568 
 2569   (* GTS_OBJECT_CLASS (gfs_gl_vof_class ())->parent_class->destroy) (o);
 2570 }
 2571 
 2572 static void gl_vof_read (GtsObject ** o, GtsFile * fp)
 2573 {
 2574   GfsGlVOF * gl = GFS_GL_VOF (*o);
 2575   GtsFileVariable var[] = {
 2576     {GTS_INT, "reversed",    TRUE},
 2577     {GTS_INT, "use_scalar",  TRUE},
 2578     {GTS_INT, "draw_edges",  TRUE},
 2579     {GTS_INT, "interpolate", TRUE},
 2580     {GTS_NONE}
 2581   };
 2582 
 2583   (* GTS_OBJECT_CLASS (gfs_gl_vof_class ())->parent_class->read) (o, fp);
 2584   if (fp->type == GTS_ERROR)
 2585     return;
 2586 
 2587   g_string_free (gl->expr, TRUE);
 2588   if (!(gl->expr = gfs_function_expression (fp, NULL)))
 2589     return;
 2590   gts_file_next_token (fp);
 2591 
 2592   var[0].data = &gl->reversed;
 2593   var[1].data = &gl->use_scalar;
 2594   var[2].data = &gl->draw_edges;
 2595   var[3].data = &gl->interpolate;
 2596   gts_file_assign_variables (fp, var);
 2597 }
 2598 
 2599 static void gl_vof_write (GtsObject * o, FILE * fp)
 2600 {
 2601   GfsGlVOF * gl = GFS_GL_VOF (o);
 2602 
 2603   (* GTS_OBJECT_CLASS (gfs_gl_vof_class ())->parent_class->write) (o, fp);
 2604 
 2605   fprintf (fp, " %s {\n"
 2606        "  reversed = %d\n"
 2607        "  use_scalar = %d\n"
 2608        "  draw_edges = %d\n"
 2609        "  interpolate = %d\n"
 2610        "}",
 2611        gl->expr->str,
 2612        gl->reversed, 
 2613        (gl->use_scalar != NULL),
 2614        gl->draw_edges,
 2615        gl->interpolate);
 2616 }
 2617 
 2618 GtsFile * gfs_gl_vof_set (GfsGlVOF * gl, const gchar * func)
 2619 {
 2620   GtsFile * fp;
 2621 
 2622   g_return_val_if_fail (gl != NULL, NULL);
 2623   g_return_val_if_fail (func != NULL, NULL);
 2624 
 2625   if ((fp = gfs_gl_var_func_set (gl->vf, GFS_GL (gl)->sim, func, gl->expr, 
 2626                  GFS_VARIABLE_CLASS (gfs_variable_tracer_vof_class ()))))
 2627     return fp;
 2628   return NULL;
 2629 }
 2630 
 2631 static void gl_vof_set_simulation (GfsGl * object, GfsSimulation * sim)
 2632 {
 2633   GfsGlVOF * gls = GFS_GL_VOF (object);
 2634   GtsFile * fp = NULL;
 2635 
 2636   (*GFS_GL_CLASS (GTS_OBJECT_CLASS (gfs_gl_vof_class ())->parent_class)->set_simulation)
 2637     (object, sim);
 2638 
 2639   if (gls->expr->str[0] == '\0' || (fp = gfs_gl_vof_set (gls, gls->expr->str))) {
 2640     GfsDomain * domain = GFS_DOMAIN (sim);
 2641     
 2642     if (domain->variables) {
 2643       GfsVariable * v = NULL;
 2644       GSList * i = domain->variables;
 2645       
 2646       while (i && !v) {
 2647     if (GFS_IS_VARIABLE_TRACER_VOF (i->data))
 2648       v = i->data;
 2649     i = i->next;
 2650       }
 2651       if (!v)
 2652     v = domain->variables->data;
 2653       gfs_gl_vof_set (gls, v->name);
 2654     }
 2655     else
 2656       gfs_gl_vof_set (gls, "0");
 2657   }
 2658   if (fp)
 2659     gts_file_destroy (fp);
 2660 }
 2661 
 2662 static void gl_vof_class_init (GfsGlClass * klass)
 2663 {
 2664   GTS_OBJECT_CLASS (klass)->destroy = gl_vof_destroy;
 2665   GTS_OBJECT_CLASS (klass)->read = gl_vof_read;
 2666   GTS_OBJECT_CLASS (klass)->write = gl_vof_write;
 2667   klass->set_simulation = gl_vof_set_simulation;
 2668   klass->draw = gl_vof_draw;
 2669   klass->cut = gl_vof_cut;
 2670   klass->pick = NULL;
 2671 }
 2672 
 2673 static void gl_vof_init (GfsGlVOF * gl)
 2674 {
 2675   gl->expr = g_string_new ("");
 2676   gl->vf = gfs_gl_var_func_new ();
 2677 #if !FTT_2D
 2678   GFS_GL (gl)->lc.r = 1.; GFS_GL (gl)->lc.g = 1.; GFS_GL (gl)->lc.b = 1.;
 2679 #endif
 2680 }
 2681 
 2682 GfsGlClass * gfs_gl_vof_class (void)
 2683 {
 2684   static GfsGlClass * klass = NULL;
 2685 
 2686   if (klass == NULL) {
 2687     GtsObjectClassInfo gfs_gl_vof_info = {
 2688       "GfsGlVOF",
 2689       sizeof (GfsGlVOF),
 2690       sizeof (GfsGlScalarClass),
 2691       (GtsObjectClassInitFunc) gl_vof_class_init,
 2692       (GtsObjectInitFunc) gl_vof_init,
 2693       (GtsArgSetFunc) NULL,
 2694       (GtsArgGetFunc) NULL
 2695     };
 2696     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl_scalar_class ()),
 2697                   &gfs_gl_vof_info);
 2698   }
 2699 
 2700   return klass;
 2701 }
 2702 
 2703 /* GfsGlVectors: Object */
 2704 
 2705 static void gl_vectors_read (GtsObject ** o, GtsFile * fp)
 2706 {
 2707   GfsGlVectors * gl = GFS_GL_VECTORS (*o);
 2708   FttComponent c;
 2709   GtsFileVariable var[] = {
 2710     {GTS_DOUBLE, "scale",      TRUE},
 2711     {GTS_INT,    "use_scalar", TRUE},
 2712     {GTS_NONE}
 2713   };
 2714 
 2715   (* GTS_OBJECT_CLASS (gfs_gl_vectors_class ())->parent_class->read) (o, fp);
 2716   if (fp->type == GTS_ERROR)
 2717     return;
 2718 
 2719   for (c = 0; c < FTT_DIMENSION; c++) {
 2720     g_string_free (gl->expr[c], TRUE);
 2721     if (!(gl->expr[c] = gfs_function_expression (fp, NULL)))
 2722       return;
 2723     gts_file_next_token (fp);
 2724   }
 2725 
 2726   var[0].data = &gl->scale;
 2727   var[1].data = &gl->use_scalar;
 2728   gts_file_assign_variables (fp, var);
 2729   if (fp->type == GTS_ERROR)
 2730     return;
 2731   if (var[0].set)
 2732     gl->already_set = TRUE;
 2733 }
 2734 
 2735 static void gl_vectors_write (GtsObject * o, FILE * fp)
 2736 {
 2737   GfsGlVectors * gl = GFS_GL_VECTORS (o);
 2738   FttComponent c;
 2739 
 2740   (* GTS_OBJECT_CLASS (gfs_gl_vectors_class ())->parent_class->write) (o, fp);
 2741 
 2742   for (c = 0; c < FTT_DIMENSION; c++)
 2743     fprintf (fp, " %s", gl->expr[c]->str);
 2744   fprintf (fp, 
 2745        " {\n"
 2746        "  scale = %g\n"
 2747        "  use_scalar = %d\n"
 2748        "}",
 2749        gl->scale, gl->use_scalar);
 2750 }
 2751 
 2752 static void gl_vectors_destroy (GtsObject * object)
 2753 {
 2754   GfsGlVectors * gl = GFS_GL_VECTORS (object);
 2755   FttComponent c;
 2756 
 2757   for (c = 0; c < FTT_DIMENSION; c++) {
 2758     gfs_gl_var_func_destroy (gl->vf[c]);
 2759     g_string_free (gl->expr[c], TRUE);
 2760   }
 2761 
 2762   (* GTS_OBJECT_CLASS (gfs_gl_vectors_class ())->parent_class->destroy) (object);
 2763 }
 2764 
 2765 static void maxv (FttCell * cell, GfsGlVectors * gl)
 2766 {
 2767   FttComponent c;
 2768   gdouble n2 = 0.;
 2769   gdouble size = ftt_cell_size (cell);
 2770 
 2771   for (c = 0; c < FTT_DIMENSION; c++)
 2772     n2 += GFS_VALUE (cell, gl->v[c])*GFS_VALUE (cell, gl->v[c]);
 2773   if (n2 > gl->max)
 2774     gl->max = n2;
 2775   if (size < gl->h)
 2776     gl->h = size;
 2777 }
 2778 
 2779 static void update_norm (GfsGlVectors * gl)
 2780 {
 2781   gl->max = -G_MAXDOUBLE;
 2782   gl->h = G_MAXDOUBLE;
 2783   gfs_domain_cell_traverse (GFS_DOMAIN (GFS_GL (gl)->sim),
 2784                 FTT_POST_ORDER, FTT_TRAVERSE_LEAFS, -1,
 2785                 (FttCellTraverseFunc) maxv, gl);
 2786   gl->max = gl->max >= 0. ? sqrt (gl->max) : 1.;
 2787   if (!gl->already_set) {
 2788     gl->scale = gl->max > 0. ? gl->h/gl->max : 1.;
 2789     gl->already_set = TRUE;
 2790   }
 2791 }
 2792 
 2793 static void gl_vectors_set_simulation (GfsGl * object, GfsSimulation * sim)
 2794 {
 2795   GfsGlVectors * gls = GFS_GL_VECTORS (object);
 2796   GfsDomain * domain;
 2797   FttComponent c;
 2798 
 2799   (*GFS_GL_CLASS (GTS_OBJECT_CLASS (gfs_gl_vectors_class ())->parent_class)->set_simulation)
 2800     (object, sim);
 2801 
 2802   domain = GFS_DOMAIN (sim);
 2803   for (c = 0; c < FTT_DIMENSION; c++) {
 2804     GtsFile * fp = NULL;
 2805 
 2806     if (gls->expr[c]->str[0] == '\0' ||
 2807     (fp = gfs_gl_var_func_set (gls->vf[c], sim, gls->expr[c]->str, NULL, NULL))) {
 2808       GfsVariable ** u = gfs_domain_velocity (domain);
 2809 
 2810       if (u)
 2811     gfs_gl_var_func_set (gls->vf[c], sim, u[c]->name, gls->expr[c], NULL);
 2812       else if (domain->variables) {
 2813     GfsVariable * v = domain->variables->data;
 2814     gfs_gl_var_func_set (gls->vf[c], sim, v->name, gls->expr[c], NULL);
 2815       }
 2816       else
 2817     gfs_gl_var_func_set (gls->vf[c], sim, "0", gls->expr[c], NULL);
 2818     }
 2819     if (fp)
 2820       gts_file_destroy (fp);
 2821     gls->v[c] = gls->vf[c]->v;
 2822   }
 2823   update_norm (gls);
 2824 }
 2825 
 2826 static void gl_vectors_draw (GfsGl * gl, GfsFrustum * f)
 2827 {
 2828   GfsGlVectors * gls = GFS_GL_VECTORS (gl);
 2829   
 2830   gl->size = 0;
 2831   glMatrixMode (GL_PROJECTION);
 2832   glPushMatrix ();
 2833   glTranslatef (0., 0., gl->p->lc);
 2834   glBegin (GL_LINES);
 2835   gfs_gl_normal (gl);
 2836   gfs_gl_cell_traverse_visible_plane (gl, f, (FttCellTraverseFunc) gl_vector, gl);
 2837   glEnd ();
 2838   glPopMatrix ();
 2839 
 2840   if (gls->use_scalar)
 2841     (* GFS_GL_CLASS (GTS_OBJECT (gl)->klass->parent_class)->draw) (gl, f);  
 2842 }
 2843 
 2844 static void gl_vectors_class_init (GfsGlClass * klass)
 2845 {
 2846   GTS_OBJECT_CLASS (klass)->read = gl_vectors_read;
 2847   GTS_OBJECT_CLASS (klass)->write = gl_vectors_write;
 2848   GTS_OBJECT_CLASS (klass)->destroy = gl_vectors_destroy;
 2849   klass->set_simulation = gl_vectors_set_simulation;
 2850   klass->draw = gl_vectors_draw;
 2851   klass->pick = NULL;
 2852 }
 2853 
 2854 static void gl_vectors_init (GfsGlVectors * object)
 2855 {
 2856   FttComponent c;
 2857   
 2858   for (c = 0; c < FTT_DIMENSION; c++) {
 2859     object->vf[c] = gfs_gl_var_func_new ();
 2860     object->expr[c] = g_string_new ("");
 2861   }
 2862   object->scale = 1.;
 2863 }
 2864 
 2865 GfsGlClass * gfs_gl_vectors_class (void)
 2866 {
 2867   static GfsGlClass * klass = NULL;
 2868 
 2869   if (klass == NULL) {
 2870     GtsObjectClassInfo gfs_gl_vectors_info = {
 2871       "GfsGlVectors",
 2872       sizeof (GfsGlVectors),
 2873       sizeof (GfsGlScalarClass),
 2874       (GtsObjectClassInitFunc) gl_vectors_class_init,
 2875       (GtsObjectInitFunc) gl_vectors_init,
 2876       (GtsArgSetFunc) NULL,
 2877       (GtsArgGetFunc) NULL
 2878     };
 2879     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl_scalar_class ()),
 2880                   &gfs_gl_vectors_info);
 2881   }
 2882 
 2883   return klass;
 2884 }
 2885 
 2886 GtsFile * gfs_gl_vectors_set (GfsGlVectors * gl, FttComponent c, const gchar * func)
 2887 {
 2888   GtsFile * fp;
 2889 
 2890   g_return_val_if_fail (gl != NULL, NULL);
 2891   g_return_val_if_fail (c < FTT_DIMENSION, NULL);
 2892   g_return_val_if_fail (func != NULL, NULL);
 2893 
 2894   if ((fp = gfs_gl_var_func_set (gl->vf[c], GFS_GL (gl)->sim, func, gl->expr[c], NULL)))
 2895     return fp;
 2896 
 2897   gl->v[c] = gl->vf[c]->v;
 2898   update_norm (gl);
 2899 
 2900   return NULL;
 2901 }
 2902 
 2903 /* GfsGlStreamline: Object */
 2904 
 2905 static void gl_streamline_destroy (GtsObject * o)
 2906 {
 2907   gfs_streamline_destroy (GFS_GL_STREAMLINE (o)->l);
 2908   glDeleteLists (GFS_GL_STREAMLINE (o)->list, 1);
 2909 
 2910   (* GTS_OBJECT_CLASS (gfs_gl_streamline_class ())->parent_class->destroy) (o);
 2911 }
 2912 
 2913 static void gl_streamline_read (GtsObject ** o, GtsFile * fp)
 2914 {
 2915   GfsGlStreamline * gl = GFS_GL_STREAMLINE (*o);
 2916   GtsFileVariable var[] = {
 2917     {GTS_DOUBLE, "x", TRUE},
 2918     {GTS_DOUBLE, "y", TRUE},
 2919     {GTS_DOUBLE, "z", TRUE},
 2920     {GTS_NONE}
 2921   };
 2922 
 2923   var[0].data = &gl->c.x;
 2924   var[1].data = &gl->c.y;
 2925   var[2].data = &gl->c.z;
 2926   gts_file_assign_variables (fp, var);
 2927   if (fp->type == GTS_ERROR)
 2928     return;
 2929 }
 2930 
 2931 static void gl_streamline_write (GtsObject * o, FILE * fp)
 2932 {
 2933   GfsGlStreamline * gl = GFS_GL_STREAMLINE (o);
 2934 
 2935   fprintf (fp,
 2936        " {\n"
 2937        "    x = %.16f y = %.16f z = %.16f\n"
 2938        "  }",
 2939        gl->c.x, gl->c.y, gl->c.z);
 2940 }
 2941 
 2942 static void gl_streamline_class_init (GtsObjectClass * klass)
 2943 {
 2944   klass->destroy = gl_streamline_destroy;
 2945   klass->read = gl_streamline_read;
 2946   klass->write = gl_streamline_write;
 2947 }
 2948 
 2949 GtsObjectClass * gfs_gl_streamline_class (void)
 2950 {
 2951   static GtsObjectClass * klass = NULL;
 2952 
 2953   if (klass == NULL) {
 2954     GtsObjectClassInfo gfs_gl_streamline_info = {
 2955       "GfsGlStreamline",
 2956       sizeof (GfsGlStreamline),
 2957       sizeof (GtsObjectClass),
 2958       (GtsObjectClassInitFunc) gl_streamline_class_init,
 2959       (GtsObjectInitFunc) NULL,
 2960       (GtsArgSetFunc) NULL,
 2961       (GtsArgGetFunc) NULL
 2962     };
 2963     klass = gts_object_class_new (gts_object_class (), &gfs_gl_streamline_info);
 2964   }
 2965 
 2966   return klass;
 2967 }
 2968 
 2969 static void gfs_gl_streamline_write (GfsGlStreamline * s, FILE * fp)
 2970 {
 2971   (* GTS_OBJECT (s)->klass->write) (GTS_OBJECT (s), fp);
 2972 }
 2973 
 2974 /* GfsGlStreamlines: Object */
 2975 
 2976 static void gl_streamlines_destroy (GtsObject * o)
 2977 {
 2978   GfsGlStreamlines * gl = GFS_GL_STREAMLINES (o);
 2979 
 2980   gfs_gl_streamlines_reset (gl);
 2981   g_list_foreach (gl->stream, (GFunc) gts_object_destroy, NULL);
 2982   g_list_free (gl->stream);
 2983 
 2984   if (gl->s)
 2985     gts_object_destroy (GTS_OBJECT (gl->s));
 2986 
 2987   (* GTS_OBJECT_CLASS (gfs_gl_streamlines_class ())->parent_class->destroy) (o);
 2988 }
 2989 
 2990 static void gl_streamlines_read (GtsObject ** o, GtsFile * fp)
 2991 {
 2992   GfsGlStreamlines * gl = GFS_GL_STREAMLINES (*o);
 2993   GtsFileVariable var[] = {
 2994     {GTS_INT,    "show_cells", TRUE},
 2995     {GTS_DOUBLE, "dmin",       TRUE},
 2996     {GTS_DOUBLE, "radius",     TRUE},
 2997     {GTS_NONE}
 2998   };
 2999 
 3000   (* GTS_OBJECT_CLASS (gfs_gl_streamlines_class ())->parent_class->read) (o, fp);
 3001   if (fp->type == GTS_ERROR)
 3002     return;
 3003   
 3004   var[0].data = &gl->show_cells;
 3005   var[1].data = &gl->dmin;
 3006   var[2].data = &gl->radius;
 3007   gts_file_assign_variables (fp, var);
 3008   if (fp->type == GTS_ERROR)
 3009     return;
 3010 
 3011   if (fp->type != '{') {
 3012     gts_file_error (fp, "expecting an opening brace");
 3013     return;
 3014   }
 3015   fp->scope_max++;
 3016   gts_file_next_token (fp);
 3017   while (fp->type == '\n')
 3018     gts_file_next_token (fp);
 3019   while (fp->type == '{') {
 3020     GtsObject * o = gts_object_new (gfs_gl_streamline_class ());
 3021 
 3022     (* o->klass->read) (&o, fp);
 3023     if (fp->type == GTS_ERROR) {
 3024       gts_object_destroy (o);
 3025       return;
 3026     }
 3027     gl->stream = g_list_append (gl->stream, o);
 3028   }
 3029   while (fp->type == '\n')
 3030     gts_file_next_token (fp);
 3031   if (fp->type != '}') {
 3032     gts_file_error (fp, "expecting a closing brace");
 3033     return;
 3034   }
 3035   fp->scope_max--;
 3036   gts_file_next_token (fp);
 3037 }
 3038 
 3039 static void gl_streamlines_write (GtsObject * o, FILE * fp)
 3040 {
 3041   GfsGlStreamlines * gl = GFS_GL_STREAMLINES (o);
 3042 
 3043   (* GTS_OBJECT_CLASS (gfs_gl_streamlines_class ())->parent_class->write) (o, fp);
 3044 
 3045   fprintf (fp, " { show_cells = %d dmin = %g radius = %g } {\n ", 
 3046        gl->show_cells, gl->dmin, gl->radius);
 3047   g_list_foreach (gl->stream, (GFunc) gfs_gl_streamline_write, fp);
 3048   fputs ("\n}", fp);
 3049 }
 3050 
 3051 static void gl_streamlines_set_simulation (GfsGl * object, GfsSimulation * sim)
 3052 {
 3053   GfsGlStreamlines * gls = GFS_GL_STREAMLINES (object);
 3054   
 3055   gfs_gl_streamlines_reset (gls);
 3056 
 3057   (*GFS_GL_CLASS (GTS_OBJECT_CLASS (gfs_gl_streamlines_class ())->parent_class)->set_simulation)
 3058     (object, sim);
 3059 
 3060   if (gls->s)
 3061     gts_object_destroy (GTS_OBJECT (gls->s));
 3062   gls->s = gfs_temporary_variable (GFS_DOMAIN (sim));
 3063   gfs_domain_cell_traverse (GFS_DOMAIN (sim), FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
 3064                 (FttCellTraverseFunc) gfs_cell_reset, gls->s);
 3065 }
 3066 
 3067 static gboolean cell_overlaps_segment (FttCell * cell, gpointer data)
 3068 {
 3069   GList * l = data;
 3070   GtsBBox bb;
 3071   GtsSegment s;
 3072 
 3073   ftt_cell_bbox (cell, &bb);
 3074   s.v1 = l->data; s.v2 = l->next->data;
 3075   return gts_bbox_overlaps_segment (&bb, &s);
 3076 }
 3077 
 3078 static void add_segment (FttCell * cell, gpointer * data)
 3079 {
 3080   GfsVariable * v = data[0];
 3081   GList * l = data[1];
 3082 
 3083   if (FTT_CELL_IS_LEAF (cell))
 3084     GFS_DOUBLE_TO_POINTER (GFS_VALUE (cell, v)) = 
 3085       g_slist_prepend (GFS_DOUBLE_TO_POINTER (GFS_VALUE (cell, v)), l);
 3086   else
 3087     GFS_DOUBLE_TO_POINTER (GFS_VALUE (cell, v)) = l;
 3088 }
 3089 
 3090 static gdouble distance2 (GList * i, GtsPoint * p, GtsSegment * closest)
 3091 {
 3092   gdouble d1 = G_MAXDOUBLE, d2 = G_MAXDOUBLE;
 3093   GtsSegment s;
 3094 
 3095   s.v1 = i->data;
 3096   if (i->next) {
 3097     s.v2 = i->next->data;
 3098     d1 = gts_point_segment_distance2 (p, &s);
 3099   }
 3100   if (i->prev) {
 3101     s.v2 = i->prev->data;
 3102     d2 = gts_point_segment_distance2 (p, &s);
 3103   }
 3104   if (closest) {
 3105     closest->v1 = s.v1;
 3106     closest->v2 = d1 < d2 ? i->next->data : i->prev->data;
 3107   }
 3108   return MIN (d1, d2);
 3109 }
 3110 
 3111 static gboolean distance_is_larger (GList * a, GList * b, guint d)
 3112 {
 3113   while (a && a != b && --d)
 3114     a = a->next;
 3115   return !d;
 3116 }
 3117 
 3118 static gdouble cell_distance2 (FttCell * cell, GtsPoint * p, gpointer data)
 3119 {
 3120   GfsGlStreamlines * gls = ((gpointer *) data)[0];
 3121   GSList * i = GFS_DOUBLE_TO_POINTER (GFS_VALUE (cell, gls->s));
 3122 
 3123   if (i == NULL)
 3124     return G_MAXDOUBLE;
 3125   else if (!FTT_CELL_IS_LEAF (cell))
 3126     return ftt_cell_point_distance2_min (cell, p);
 3127   else {
 3128     GfsGlStreamline * s = ((gpointer *) data)[1];
 3129     GList * l = ((gpointer *) data)[2];
 3130     gdouble dmin = G_MAXDOUBLE, h = ftt_cell_size (cell);
 3131 
 3132     while (i) {
 3133       GList * j = i->data;
 3134       gboolean same_stream = (GTS_OBJECT (j->data)->reserved == s);
 3135       
 3136       if (!same_stream || distance_is_larger (l, j, 16)) {
 3137     gdouble d = distance2 (j, p, NULL);
 3138 
 3139     if (d < dmin)
 3140       dmin = d;
 3141     if (same_stream && d < h*h/100.)
 3142       return 0.;
 3143       }
 3144       i = i->next;
 3145     }
 3146     return dmin;
 3147   }
 3148 }
 3149 
 3150 static gboolean stop (FttCell * cell, GList * l, gpointer data)
 3151 {
 3152   GfsGlStreamlines * gls = ((gpointer *) data)[0];
 3153   GfsGlStreamline * s = ((gpointer *) data)[1];
 3154   gboolean stop = FALSE;
 3155 
 3156   ((gpointer *) data)[2] = l;
 3157   stop = (gfs_domain_cell_point_distance2 (GFS_DOMAIN (GFS_GL (gls)->sim), l->data,
 3158                        cell_distance2, data, NULL) 
 3159       <= gls->dmin*gls->dmin/4.);
 3160   if (l->next) {
 3161     gpointer datum[2];
 3162     
 3163     datum[0] = gls->s;
 3164     datum[1] = l;
 3165     gfs_domain_cell_traverse_condition (GFS_DOMAIN (GFS_GL (gls)->sim), 
 3166                     FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
 3167                     (FttCellTraverseFunc) add_segment, datum,
 3168                     cell_overlaps_segment, l);
 3169   }
 3170   GTS_OBJECT (l->data)->reserved = s;
 3171   return stop;
 3172 }
 3173 
 3174 static void gl_streamline_update_display_list (GfsGlStreamline * s,
 3175                            GfsGlStreamlines * gls)
 3176 {
 3177   if (s->list == 0) {
 3178     s->list = glGenLists (1);
 3179     if (s->list == 0)
 3180       g_warning ("No available OpenGL display list!");
 3181   }
 3182   glNewList (s->list, GL_COMPILE);
 3183   gl_streamline_draw (s, gls);
 3184   glEndList ();
 3185 }
 3186 
 3187 static void gl_streamline_update (GfsGlStreamline * s,
 3188                   GfsGlStreamlines * gls)
 3189 {
 3190   GtsPoint p;
 3191   gpointer data[3];
 3192   
 3193   if (s->l)
 3194     return;
 3195 
 3196   data[0] = gls;
 3197   data[1] = s;
 3198   p.x = s->c.x; p.y = s->c.y; p.z = s->c.z;
 3199   if (gfs_domain_cell_point_distance2 (GFS_DOMAIN (GFS_GL (gls)->sim), &p,
 3200                        cell_distance2, data, NULL) > gls->dmin*gls->dmin) {
 3201     s->l = gfs_streamline_new (GFS_DOMAIN (GFS_GL (gls)->sim),
 3202                    GFS_GL_VECTORS (gls)->v, s->c,
 3203                    GFS_GL_VECTORS (gls)->use_scalar ? GFS_GL_SCALAR (gls)->v : NULL,
 3204                    GFS_GL_SCALAR (gls)->min, GFS_GL_SCALAR (gls)->max, FALSE,
 3205                    stop, data);
 3206     if (s->l && !s->l->next) {
 3207       gfs_streamline_destroy (s->l);
 3208       s->l = NULL;
 3209     }
 3210     else
 3211       gl_streamline_update_display_list (s, gls);
 3212   }
 3213 }
 3214 
 3215 static void gl_streamlines_draw (GfsGl * gl, GfsFrustum * f)
 3216 {
 3217   GfsGlStreamlines * gls = GFS_GL_STREAMLINES (gl);
 3218   GList * i = gls->stream;
 3219 
 3220   gl->size = 0;
 3221   glMatrixMode (GL_PROJECTION);
 3222   glPushMatrix ();
 3223   glTranslatef (0., 0., gl->p->lc/2.);
 3224   gfs_gl_normal (gl);
 3225   if (gls->show_cells) {
 3226     glColor3f (0.3, 0.3, 0.3);
 3227     gfs_gl_cell_traverse_visible_plane (gl, f, (FttCellTraverseFunc) gl_cell, gl);
 3228   }
 3229   glTranslatef (0., 0., gl->p->lc/2.);
 3230   while (i) {
 3231     GfsGlStreamline * s = i->data;
 3232 
 3233     if (i == gls->selected)
 3234       glColor3f (1., 0., 0.);
 3235     else if (!GFS_GL_VECTORS (gl)->use_scalar)
 3236       glColor3f (gl->lc.r, gl->lc.g, gl->lc.b);
 3237     gl_streamline_update (s, gls);
 3238     if (s->l) {
 3239       if (gl->format == GFSGL_PPM_OFFSCREEN)
 3240     gl_streamline_draw (s, gls);
 3241       else
 3242     glCallList (s->list);
 3243     }
 3244     gl->size++;
 3245     i = i->next;
 3246   }
 3247   glPopMatrix ();
 3248 }
 3249 
 3250 static gdouble gl_streamlines_pick (GfsGl * gl, GfsGlRay * r)
 3251 {
 3252   gint maxlevel = gl->maxlevel;
 3253   gdouble a;
 3254 
 3255   gl->maxlevel = -1;
 3256   a = gl2D_pick (gl, r);
 3257   gl->maxlevel = maxlevel;
 3258   if (a != GFS_NODATA) {
 3259     GfsGl2D * gl2 = GFS_GL2D (gl);
 3260     GfsGlStreamlines * gls = GFS_GL_STREAMLINES (gl);
 3261     GSList * i = GFS_DOUBLE_TO_POINTER (GFS_VALUE (gl2->picked, gls->s));
 3262     gdouble dmin = G_MAXDOUBLE;
 3263     GfsGlStreamline * closest = NULL;
 3264     GtsPoint p;
 3265 
 3266     p.x = gl2->pickedpos.x; p.y = gl2->pickedpos.y; p.z = gl2->pickedpos.z;
 3267     while (i) {
 3268       GList * j = i->data;
 3269       gdouble d = distance2 (j, &p, NULL);
 3270 
 3271       if (d < dmin) {
 3272     dmin = d;
 3273     closest = GTS_OBJECT (j->data)->reserved;
 3274       }
 3275       i = i->next;
 3276     }
 3277     gls->candidate = closest ? g_list_find (gls->stream, closest) : NULL;
 3278   }
 3279   return a;
 3280 }
 3281 
 3282 static void gl_streamlines_class_init (GfsGlClass * klass)
 3283 {
 3284   GTS_OBJECT_CLASS (klass)->destroy = gl_streamlines_destroy;
 3285   GTS_OBJECT_CLASS (klass)->read = gl_streamlines_read;
 3286   GTS_OBJECT_CLASS (klass)->write = gl_streamlines_write;
 3287   klass->set_simulation = gl_streamlines_set_simulation;
 3288   klass->draw = gl_streamlines_draw;
 3289   klass->pick = gl_streamlines_pick;
 3290 }
 3291 
 3292 static void gl_streamlines_init (GfsGlStreamlines * gl)
 3293 {
 3294   GtsColor c = { 1., 1., 1. };
 3295 
 3296   GFS_GL (gl)->lc = c;
 3297   gl->show_cells = TRUE;
 3298 }
 3299 
 3300 GfsGlClass * gfs_gl_streamlines_class (void)
 3301 {
 3302   static GfsGlClass * klass = NULL;
 3303 
 3304   if (klass == NULL) {
 3305     GtsObjectClassInfo gfs_gl_streamlines_info = {
 3306       "GfsGlStreamlines",
 3307       sizeof (GfsGlStreamlines),
 3308       sizeof (GfsGlScalarClass),
 3309       (GtsObjectClassInitFunc) gl_streamlines_class_init,
 3310       (GtsObjectInitFunc) gl_streamlines_init,
 3311       (GtsArgSetFunc) NULL,
 3312       (GtsArgGetFunc) NULL
 3313     };
 3314     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl_vectors_class ()),
 3315                   &gfs_gl_streamlines_info);
 3316   }
 3317 
 3318   return klass;
 3319 }
 3320 
 3321 static void reset_segments (FttCell * cell, GfsGlStreamlines * gl)
 3322 {
 3323   if (FTT_CELL_IS_LEAF (cell)) {
 3324     GSList * l = GFS_DOUBLE_TO_POINTER (GFS_VALUE (cell, gl->s));
 3325     g_slist_free (l);
 3326   }
 3327   GFS_VALUE (cell, gl->s) = 0.;
 3328 }
 3329 
 3330 void gfs_gl_streamlines_reset (GfsGlStreamlines * gl)
 3331 {
 3332   GList * i;
 3333 
 3334   g_return_if_fail (gl != NULL);
 3335 
 3336   if (GFS_GL (gl)->sim)
 3337     gfs_domain_cell_traverse (GFS_DOMAIN (GFS_GL (gl)->sim), 
 3338                   FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
 3339                   (FttCellTraverseFunc) reset_segments, gl);
 3340   i = gl->stream;
 3341   while (i) {
 3342     GfsGlStreamline * s = i->data;
 3343     gfs_streamline_destroy (s->l);
 3344     s->l = NULL;    
 3345     i = i->next;
 3346   }
 3347 }
 3348 
 3349 static void remove_segment (FttCell * cell, gpointer * data)
 3350 {
 3351   GfsVariable * v = data[0];
 3352 
 3353   if (FTT_CELL_IS_LEAF (cell)) {
 3354     GList * i = data[1];
 3355     GSList * l = GFS_DOUBLE_TO_POINTER (GFS_VALUE (cell, v));
 3356     l = g_slist_remove (l, i);
 3357     l = g_slist_remove (l, i->next);
 3358     GFS_DOUBLE_TO_POINTER (GFS_VALUE (cell, v)) = l;
 3359   }
 3360   else {
 3361     FttCellChildren child;
 3362     gboolean empty = TRUE;
 3363     guint i;
 3364     
 3365     ftt_cell_children (cell, &child);
 3366     for (i = 0; i < FTT_CELLS && empty; i++)
 3367       if (child.c[i] && GFS_VALUE (child.c[i], v) != 0.)
 3368     empty = FALSE;
 3369     if (empty)
 3370       GFS_VALUE (cell, v) = 0.;
 3371   }
 3372 }
 3373 
 3374 void gfs_gl_streamlines_reset_selected (GfsGlStreamlines * gl)
 3375 {
 3376   GfsGlStreamline * s;
 3377   gpointer datum[2];
 3378   GList * i;
 3379 
 3380   g_return_if_fail (gl != NULL);
 3381   g_return_if_fail (gl->selected != NULL);
 3382 
 3383   s = gl->selected->data;
 3384   datum[0] = gl->s;
 3385   i = s->l;
 3386   while (i && i->next) {
 3387     datum[1] = i;
 3388     gfs_domain_cell_traverse_condition (GFS_DOMAIN (GFS_GL (gl)->sim), 
 3389                     FTT_POST_ORDER, FTT_TRAVERSE_ALL, -1,
 3390                     (FttCellTraverseFunc) remove_segment, datum,
 3391                     cell_overlaps_segment, i);
 3392     i = i->next;
 3393   }
 3394   gfs_streamline_destroy (s->l);
 3395   s->l = NULL;
 3396 }
 3397 
 3398 GfsGlStreamline * gfs_gl_streamlines_add (GfsGlStreamlines * gl, FttVector p)
 3399 {
 3400   GfsGlStreamline * s;
 3401 
 3402   g_return_val_if_fail (gl != NULL, NULL);
 3403 
 3404   s = GFS_GL_STREAMLINE (gts_object_new (gfs_gl_streamline_class ()));
 3405   s->c = p;
 3406   gl_streamline_update (s, gl);
 3407   if (s->l) {
 3408     gl->stream = g_list_append (gl->stream, s);
 3409     gl->selected = g_list_last (gl->stream);
 3410     return s;
 3411   }
 3412   else {
 3413     gts_object_destroy (GTS_OBJECT (s));
 3414     return NULL;
 3415   }
 3416 }
 3417 
 3418 gboolean gfs_gl_streamlines_remove_selected (GfsGlStreamlines * gl)
 3419 {
 3420   g_return_val_if_fail (gl != NULL, FALSE);
 3421 
 3422   if (gl->selected) {
 3423     gfs_gl_streamlines_reset_selected (gl);
 3424     gts_object_destroy (gl->selected->data);
 3425     gl->stream = g_list_remove_link (gl->stream, gl->selected);
 3426     g_list_free (gl->selected);
 3427     gl->selected = NULL;
 3428     return TRUE;
 3429   }
 3430   return FALSE;
 3431 }
 3432 
 3433 void gfs_gl_streamlines_update_display_lists (GfsGlStreamlines * gl)
 3434 {
 3435   g_return_if_fail (gl != NULL);
 3436 
 3437   g_list_foreach (gl->stream, (GFunc) gl_streamline_update_display_list, gl);
 3438 }
 3439 
 3440 gdouble gfs_gl_streamlines_closest (GfsGlStreamlines * gl, 
 3441                     FttVector * p,
 3442                     GtsPoint * closest)
 3443 {
 3444   FttCell * cell = NULL;
 3445   GtsPoint p1;
 3446   gdouble dmin;
 3447   gpointer data[3];
 3448 
 3449   g_return_val_if_fail (gl != NULL, G_MAXDOUBLE);
 3450   g_return_val_if_fail (p != NULL, G_MAXDOUBLE);
 3451   g_return_val_if_fail (closest != NULL, G_MAXDOUBLE);
 3452 
 3453   data[0] = gl; data[1] = data[2] = NULL;
 3454   p1.x = p->x; p1.y = p->y; p1.z = p->z;
 3455   dmin = gfs_domain_cell_point_distance2 (GFS_DOMAIN (GFS_GL (gl)->sim), &p1,
 3456                       cell_distance2, data, &cell);
 3457   if (cell) {
 3458     GSList * i = GFS_DOUBLE_TO_POINTER (GFS_VALUE (cell, gl->s));
 3459     GtsSegment segment;
 3460     gdouble dmin = G_MAXDOUBLE;
 3461     
 3462     while (i) {
 3463       GList * j = i->data;
 3464       GtsSegment s1;
 3465       gdouble d = distance2 (j, &p1, &s1);
 3466       
 3467       if (d < dmin) {
 3468     dmin = d;
 3469     segment = s1;
 3470       }
 3471       i = i->next;
 3472     }
 3473     gts_point_segment_closest (&p1, &segment, closest);
 3474   }
 3475   return dmin;
 3476 }
 3477 
 3478 static GList * streamline_step (GList * i, GtsPoint * p, gdouble ds)
 3479 {
 3480   if (!i || !i->next)
 3481     return NULL;
 3482   if (p->x == G_MAXDOUBLE) {
 3483     GtsPoint * p1 = i->data;
 3484     p->x = p1->x; p->y = p1->y; p->z = p1->z;
 3485     return i;
 3486   }
 3487   else {
 3488     gdouble l = gts_point_distance (i->data, i->next->data);
 3489     gdouble s;
 3490   
 3491     g_assert (l > 0.);
 3492     s = (gts_point_distance (i->data, p) + ds)/l;
 3493     while (s > 1. && i->next->next) {
 3494       s = l*(s - 1.);
 3495       i = i->next;
 3496       l = gts_point_distance (i->data, i->next->data);
 3497       g_assert (l > 0.);
 3498       s /= l;
 3499     }
 3500     if (s > 1.)
 3501       return NULL;
 3502     else {
 3503       GtsPoint * p1 = i->data;
 3504       GtsPoint * p2 = i->next->data;
 3505       p->x = p1->x + s*(p2->x - p1->x);
 3506       p->y = p1->y + s*(p2->y - p1->y);
 3507       p->z = p1->z + s*(p2->z - p1->z);
 3508       return i;
 3509     }
 3510   }
 3511 }
 3512 
 3513 static void push_stream (GfsGlStreamline * s, GtsFifo * fifo)
 3514 {
 3515   gts_fifo_push (fifo, s);
 3516 }
 3517 
 3518 void gfs_gl_streamlines_evenly_spaced (GfsGlStreamlines * gl,
 3519                        gboolean (* callback) (GfsGlStreamlines *, gpointer),
 3520                        gpointer data)
 3521 {
 3522   GtsFifo * fifo;
 3523   GfsGlStreamline * s;
 3524   gboolean stop = FALSE;
 3525 
 3526   g_return_if_fail (gl != NULL);
 3527   g_return_if_fail (gl->dmin > 0.);
 3528 
 3529   fifo = gts_fifo_new ();
 3530   g_list_foreach (gl->stream, (GFunc) push_stream, fifo);
 3531   while ((s = gts_fifo_pop (fifo)) && !stop) {
 3532     GtsPoint p;
 3533     GList * i = s->l;
 3534 
 3535     g_assert (i);
 3536     p.x = G_MAXDOUBLE;
 3537     while ((i = streamline_step (i, &p, gl->dmin/10.))) {
 3538       GtsVector n;
 3539       FttVector p1;
 3540       GfsGlStreamline * s1;
 3541 
 3542       gts_vector_init (n, GTS_POINT (i->data), GTS_POINT (i->next->data));
 3543       gts_vector_normalize (n);
 3544 
 3545       p1.x = p.x + n[1]*gl->dmin;
 3546       p1.y = p.y - n[0]*gl->dmin;
 3547       p1.z = 0.;
 3548       if ((s1 = gfs_gl_streamlines_add (gl, p1))) {
 3549     if (callback)
 3550       stop |= (* callback) (gl, data);
 3551     gts_fifo_push (fifo, s1);
 3552       }
 3553 
 3554       p1.x = p.x - n[1]*gl->dmin;
 3555       p1.y = p.y + n[0]*gl->dmin;
 3556       p1.z = 0.;
 3557       if ((s1 = gfs_gl_streamlines_add (gl, p1))) {
 3558     if (callback)
 3559       stop |= (* callback) (gl, data);  
 3560     gts_fifo_push (fifo, s1);
 3561       }
 3562     }
 3563   }
 3564   gts_fifo_destroy (fifo);
 3565   gl->selected = NULL;
 3566 }
 3567 
 3568 /* GfsGlEllipses: Object */
 3569 
 3570 static void gl_ellipses_read (GtsObject ** o, GtsFile * fp)
 3571 {
 3572   GfsGlEllipses * gl = GFS_GL_ELLIPSES (*o);
 3573   GtsFileVariable var[] = {
 3574     {GTS_DOUBLE, "scale",      TRUE},
 3575     {GTS_INT,    "use_scalar", TRUE},
 3576     {GTS_NONE}
 3577   };
 3578   guint i;
 3579 
 3580   (* GTS_OBJECT_CLASS (gfs_gl_ellipses_class ())->parent_class->read) (o, fp);
 3581   if (fp->type == GTS_ERROR)
 3582     return;
 3583 
 3584   for (i = 0; i < 4; i++) {
 3585     g_string_free (gl->expr[i], TRUE);
 3586     if (!(gl->expr[i] = gfs_function_expression (fp, NULL)))
 3587       return;
 3588     gts_file_next_token (fp);
 3589   }
 3590 
 3591   var[0].data = &gl->scale;
 3592   var[1].data = &gl->use_scalar;
 3593   gts_file_assign_variables (fp, var);
 3594   if (fp->type == GTS_ERROR)
 3595     return;
 3596   if (var[0].set)
 3597     gl->already_set = TRUE;
 3598 }
 3599 
 3600 static void gl_ellipses_write (GtsObject * o, FILE * fp)
 3601 {
 3602   GfsGlEllipses * gl = GFS_GL_ELLIPSES (o);
 3603   guint i;
 3604 
 3605   (* GTS_OBJECT_CLASS (gfs_gl_ellipses_class ())->parent_class->write) (o, fp);
 3606 
 3607   for (i = 0; i < 4; i++)
 3608     fprintf (fp, " %s", gl->expr[i]->str);
 3609   fprintf (fp, " {\n"
 3610        "  scale = %g\n"
 3611        "  use_scalar = %d\n"
 3612        "}",
 3613        gl->scale, gl->use_scalar);
 3614 }
 3615 
 3616 static void gl_ellipses_destroy (GtsObject * o)
 3617 {
 3618   GfsGlEllipses * gl = GFS_GL_ELLIPSES (o);
 3619   guint i;
 3620 
 3621   for (i = 0; i < 4; i++) {
 3622     gfs_gl_var_func_destroy (gl->vf[i]);
 3623     g_string_free (gl->expr[i], TRUE);
 3624   }
 3625 
 3626   (* GTS_OBJECT_CLASS (gfs_gl_ellipses_class ())->parent_class->destroy) (o);
 3627 }
 3628 
 3629 static void maxe (FttCell * cell, GfsGlEllipses * gl)
 3630 {
 3631   guint i;
 3632   gdouble n2 = 0.;
 3633   gdouble size = ftt_cell_size (cell);
 3634 
 3635   for (i = 0; i < 4; i++)
 3636     n2 += GFS_VALUE (cell, gl->v[i])*GFS_VALUE (cell, gl->v[i]);
 3637   if (n2 > gl->max)
 3638     gl->max = n2;
 3639   if (size < gl->h)
 3640     gl->h = size;
 3641 }
 3642 
 3643 static void update_ellipses_norm (GfsGlEllipses * gl)
 3644 {
 3645   gl->max = -G_MAXDOUBLE;
 3646   gl->h = G_MAXDOUBLE;
 3647   gfs_domain_cell_traverse (GFS_DOMAIN (GFS_GL (gl)->sim),
 3648                 FTT_POST_ORDER, FTT_TRAVERSE_LEAFS, -1,
 3649                 (FttCellTraverseFunc) maxe, gl);
 3650   gl->max = gl->max >= 0. ? sqrt (gl->max) : 1.;
 3651   if (!gl->already_set) {
 3652     gl->scale = gl->max > 0. ? gl->h/gl->max : 1.;
 3653     gl->already_set = TRUE;
 3654   }
 3655 }
 3656 
 3657 static void gl_ellipses_set_simulation (GfsGl * object, GfsSimulation * sim)
 3658 {
 3659   GfsGlEllipses * gle = GFS_GL_ELLIPSES (object);
 3660   GfsDomain * domain;
 3661   FttComponent c;
 3662 
 3663   (*GFS_GL_CLASS (GTS_OBJECT_CLASS (gfs_gl_ellipses_class ())->parent_class)->set_simulation)
 3664     (object, sim);
 3665 
 3666   domain = GFS_DOMAIN (sim);
 3667   for (c = 0; c < 4; c++) {
 3668     GtsFile * fp = NULL;
 3669 
 3670     if (gle->expr[c]->str[0] == '\0' ||
 3671     (fp = gfs_gl_var_func_set (gle->vf[c], sim, gle->expr[c]->str, NULL, NULL))) {
 3672       if (domain->variables) {
 3673     GfsVariable * v = domain->variables->data;
 3674     gfs_gl_var_func_set (gle->vf[c], sim, v->name, gle->expr[c], NULL);
 3675       }
 3676       else
 3677     gfs_gl_var_func_set (gle->vf[c], sim, "0", gle->expr[c], NULL);
 3678     }
 3679     if (fp)
 3680       gts_file_destroy (fp);
 3681     gle->v[c] = gle->vf[c]->v;
 3682   }
 3683   update_ellipses_norm (gle);
 3684 }
 3685 
 3686 static void gl_ellipses_draw (GfsGl * gl, GfsFrustum * f)
 3687 {
 3688   GfsGlEllipses * gls = GFS_GL_ELLIPSES (gl);
 3689   
 3690   gl->size = 0;
 3691   glMatrixMode (GL_PROJECTION);
 3692   glPushMatrix ();
 3693   glTranslatef (0., 0., gl->p->lc);
 3694   gfs_gl_normal (gl);
 3695   gfs_gl_cell_traverse_visible_plane (gl, f, (FttCellTraverseFunc) gl_ellipse, gl);
 3696   glPopMatrix ();
 3697 
 3698   if (gls->use_scalar)
 3699     (* GFS_GL_CLASS (GTS_OBJECT (gl)->klass->parent_class)->draw) (gl, f);  
 3700 }
 3701 
 3702 static void gl_ellipses_class_init (GfsGlClass * klass)
 3703 {
 3704   GTS_OBJECT_CLASS (klass)->read = gl_ellipses_read;
 3705   GTS_OBJECT_CLASS (klass)->write = gl_ellipses_write;
 3706   GTS_OBJECT_CLASS (klass)->destroy = gl_ellipses_destroy;
 3707   klass->set_simulation = gl_ellipses_set_simulation;
 3708   klass->draw = gl_ellipses_draw;
 3709   klass->pick = NULL;
 3710 }
 3711 
 3712 static void gl_ellipses_init (GfsGlEllipses * object)
 3713 {
 3714   FttComponent c;
 3715   
 3716   for (c = 0; c < 4; c++) {
 3717     object->vf[c] = gfs_gl_var_func_new ();
 3718     object->expr[c] = g_string_new ("");
 3719   }
 3720   object->scale = 1.;
 3721 }
 3722 
 3723 GfsGlClass * gfs_gl_ellipses_class (void)
 3724 {
 3725   static GfsGlClass * klass = NULL;
 3726 
 3727   if (klass == NULL) {
 3728     GtsObjectClassInfo gfs_gl_ellipses_info = {
 3729       "GfsGlEllipses",
 3730       sizeof (GfsGlEllipses),
 3731       sizeof (GfsGlScalarClass),
 3732       (GtsObjectClassInitFunc) gl_ellipses_class_init,
 3733       (GtsObjectInitFunc) gl_ellipses_init,
 3734       (GtsArgSetFunc) NULL,
 3735       (GtsArgGetFunc) NULL
 3736     };
 3737     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl_scalar_class ()),
 3738                   &gfs_gl_ellipses_info);
 3739   }
 3740 
 3741   return klass;
 3742 }
 3743 
 3744 GtsFile * gfs_gl_ellipses_set (GfsGlEllipses * gl, guint i, const gchar * func)
 3745 {
 3746   GtsFile * fp;
 3747   
 3748   g_return_val_if_fail (gl != NULL, NULL);
 3749   g_return_val_if_fail (i < 4, NULL);
 3750   g_return_val_if_fail (func != NULL, NULL);
 3751 
 3752   if ((fp = gfs_gl_var_func_set (gl->vf[i], GFS_GL (gl)->sim, func, gl->expr[i], NULL)))
 3753     return fp;
 3754 
 3755   gl->v[i] = gl->vf[i]->v;
 3756   update_ellipses_norm (gl);
 3757 
 3758   return NULL;
 3759 }
 3760 
 3761 /* GfsGlLocation: Object */
 3762 
 3763 static void gl_location_read (GtsObject ** o, GtsFile * fp)
 3764 {
 3765   GtsFileVariable var[] = {
 3766     {GTS_DOUBLE, "size",  TRUE},
 3767     {GTS_INT,    "label", TRUE},
 3768     {GTS_NONE}
 3769   };
 3770 
 3771   (* GTS_OBJECT_CLASS (gfs_gl_location_class ())->parent_class->read) (o, fp);
 3772   if (fp->type == GTS_ERROR)
 3773     return;
 3774   
 3775   var[0].data = &GFS_GL_LOCATION (*o)->size;
 3776   var[1].data = &GFS_GL_LOCATION (*o)->label;
 3777   gts_file_assign_variables (fp, var);
 3778   if (fp->type == GTS_ERROR)
 3779     return;
 3780 }
 3781 
 3782 static void gl_location_write (GtsObject * o, FILE * fp)
 3783 {
 3784   (* GTS_OBJECT_CLASS (gfs_gl_location_class ())->parent_class->write) (o, fp);
 3785 
 3786   fprintf (fp, " {\n  size = %g\n  label = %d\n}", 
 3787        GFS_GL_LOCATION (o)->size, GFS_GL_LOCATION (o)->label);
 3788 }
 3789 
 3790 static void gl_location_draw (GfsGl * gl, GfsFrustum * f)
 3791 {
 3792   GSList * i;
 3793   gdouble size;
 3794 
 3795   i = gl->sim->events->items;
 3796   while (i) {
 3797     if (GFS_IS_OUTPUT_LOCATION (i->data)) {
 3798       GfsOutputLocation * l = i->data;
 3799       guint j;
 3800 
 3801       for (j = 0; j < l->p->len; j++) {
 3802     FttVector p = g_array_index (l->p, FttVector, j);
 3803     gfs_simulation_map (gl->sim, &p);
 3804     FttCell * cell = gfs_domain_locate (GFS_DOMAIN (gl->sim), p, -1, NULL);
 3805 
 3806     if (gl->maxlevel >= 0)
 3807       size = GFS_GL_LOCATION (gl)->size*ftt_level_size (gl->maxlevel);
 3808     else if (cell)
 3809       size = GFS_GL_LOCATION (gl)->size*ftt_cell_size (cell);
 3810     else
 3811       size = GFS_GL_LOCATION (gl)->size/64.;
 3812 
 3813     if (cell)
 3814       glColor3f (gl->lc.r, gl->lc.g, gl->lc.b);
 3815     else
 3816       glColor3f (1., 0., 0.);
 3817 
 3818     glMatrixMode (GL_MODELVIEW);
 3819     glPushMatrix ();
 3820     glTranslated (p.x, p.y, p.z);
 3821     glScaled (size, size, size);
 3822     draw_sphere ();
 3823     glPopMatrix ();
 3824     
 3825     if (j == 0 && GFS_GL_LOCATION (gl)->label)
 3826       gl_draw_text (gl, l->label ? l->label : GFS_OUTPUT (l)->format,
 3827             p.x + size/2., p.y + size/2., p.z + size/2.,
 3828             size*gl->font_size);
 3829       }
 3830     }
 3831     i = i->next;
 3832   }
 3833 }
 3834 
 3835 static void gl_location_init (GfsGl * gl)
 3836 {
 3837   GtsColor c = { 0., 0., 1. };
 3838 
 3839   gl->lc = c;
 3840   gl->font_size = 2.;
 3841 
 3842   GFS_GL_LOCATION (gl)->size = 1.;
 3843 }
 3844 
 3845 static gboolean gl_location_relevant (GfsSimulation * sim)
 3846 {
 3847   GSList * i = sim->events->items;
 3848 
 3849   while (i) {
 3850     if (GFS_IS_OUTPUT_LOCATION (i->data))
 3851       return TRUE;
 3852     i = i->next;
 3853   }
 3854   return FALSE;
 3855 }
 3856 
 3857 static void gl_location_class_init (GfsGlClass * klass)
 3858 {
 3859   GTS_OBJECT_CLASS (klass)->read = gl_location_read;
 3860   GTS_OBJECT_CLASS (klass)->write = gl_location_write;
 3861   klass->draw = gl_location_draw;
 3862   klass->relevant = gl_location_relevant;
 3863 }
 3864 
 3865 GfsGlClass * gfs_gl_location_class (void)
 3866 {
 3867   static GfsGlClass * klass = NULL;
 3868 
 3869   if (klass == NULL) {
 3870     GtsObjectClassInfo gfs_gl_location_info = {
 3871       "GfsGlLocation",
 3872       sizeof (GfsGlLocation),
 3873       sizeof (GfsGlClass),
 3874       (GtsObjectClassInitFunc) gl_location_class_init,
 3875       (GtsObjectInitFunc) gl_location_init,
 3876       (GtsArgSetFunc) NULL,
 3877       (GtsArgGetFunc) NULL
 3878     };
 3879     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl_class ()), &gfs_gl_location_info);
 3880   }
 3881 
 3882   return klass;
 3883 }
 3884 
 3885 /* GfsGlHeight: Object */
 3886 
 3887 typedef struct {
 3888   GfsGl * gl;
 3889   GfsVariableTracerVOFHeight * t;
 3890 } HFData;
 3891 
 3892 static void draw_height (GfsGl * gl, FttComponent c, 
 3893              FttVector p, gdouble size,
 3894              gdouble h, const char * label)
 3895 {
 3896   glColor3f (gl->lc.r, gl->lc.g, gl->lc.b);
 3897   glMatrixMode (GL_MODELVIEW);
 3898   glPushMatrix ();
 3899   FttVector o = p;
 3900   (&o.x)[c] += h;
 3901   glTranslated (o.x, o.y, o.z);
 3902   glScaled (size, size, size);
 3903   draw_sphere ();
 3904   glPopMatrix ();
 3905 
 3906   if (GFS_GL_LOCATION (gl)->label)
 3907     gl_draw_text (gl, label,
 3908           o.x + size/2., o.y + size/2., o.z + size/2.,
 3909           size*gl->font_size);
 3910   
 3911   gl->size++;
 3912 }
 3913 
 3914 static void gl_height (FttCell * cell, HFData * hf)
 3915 {
 3916   GfsGl * gl = hf->gl;
 3917   FttVector p;
 3918   gdouble h = ftt_cell_size (cell);
 3919   ftt_cell_pos (cell, &p);
 3920   
 3921   gdouble size = GFS_GL_LOCATION (gl)->size;
 3922   if (gl->maxlevel >= 0)
 3923     size *= ftt_level_size (gl->maxlevel);
 3924   else
 3925     size *= h;
 3926 
 3927   FttComponent c;
 3928   for (c = 0; c < FTT_DIMENSION; c++) {
 3929     static gchar hb[3][4] = {"Hbx", "Hby", "Hbz"};
 3930     static gchar ht[3][4] = {"Htx", "Hty", "Htz"};
 3931     if (GFS_HAS_DATA (cell, hf->t->hb[c]))
 3932       draw_height (gl, c, p, size, GFS_VALUE (cell, hf->t->hb[c])*h, hb[c]);
 3933     if (GFS_HAS_DATA (cell, hf->t->ht[c]))
 3934       draw_height (gl, c, p, size, - GFS_VALUE (cell, hf->t->ht[c])*h, ht[c]);
 3935   }
 3936 }
 3937 
 3938 static void gl_height_draw (GfsGl * gl, GfsFrustum * f)
 3939 {
 3940   HFData hf = { gl, NULL };
 3941   GSList * i = GFS_DOMAIN (gl->sim)->variables;
 3942 
 3943   while (i && !hf.t) {
 3944     if (GFS_IS_VARIABLE_TRACER_VOF_HEIGHT (i->data))
 3945       hf.t = i->data;
 3946     i = i->next;
 3947   }
 3948   if (hf.t == NULL)
 3949     return;
 3950 
 3951   gl->size = 0;
 3952   gfs_gl_cell_traverse_visible (gl, f, (FttCellTraverseFunc) gl_height, &hf);
 3953 }
 3954 
 3955 static void gl_height_init (GfsGl * gl)
 3956 {
 3957   GtsColor c = { 0., 1., 0. };
 3958   gl->lc = c;
 3959   GFS_GL_LOCATION (gl)->size = 0.2;
 3960 }
 3961 
 3962 static gboolean gl_height_relevant (GfsSimulation * sim)
 3963 {
 3964   GSList * i = GFS_DOMAIN (sim)->variables;
 3965 
 3966   while (i) {
 3967     if (GFS_IS_VARIABLE_TRACER_VOF_HEIGHT (i->data))
 3968       return TRUE;
 3969     i = i->next;
 3970   }
 3971   return FALSE;
 3972 }
 3973 
 3974 static void gl_height_class_init (GfsGlClass * klass)
 3975 {
 3976   klass->draw = gl_height_draw;
 3977   klass->relevant = gl_height_relevant;
 3978 }
 3979 
 3980 GfsGlClass * gfs_gl_height_class (void)
 3981 {
 3982   static GfsGlClass * klass = NULL;
 3983 
 3984   if (klass == NULL) {
 3985     GtsObjectClassInfo info = {
 3986       "GfsGlHeight",
 3987       sizeof (GfsGlLocation),
 3988       sizeof (GfsGlClass),
 3989       (GtsObjectClassInitFunc) gl_height_class_init,
 3990       (GtsObjectInitFunc) gl_height_init,
 3991       (GtsArgSetFunc) NULL,
 3992       (GtsArgGetFunc) NULL
 3993     };
 3994     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl_location_class ()), &info);
 3995   }
 3996 
 3997   return klass;
 3998 }
 3999 
 4000 /* GfsGlLocate: Object */
 4001 
 4002 static void gl_locate_read (GtsObject ** o, GtsFile * fp)
 4003 {
 4004   GfsGlLocate * gl = GFS_GL_LOCATE (*o);
 4005   GtsFileVariable var[] = {
 4006     {GTS_DOUBLE, "x", TRUE},
 4007     {GTS_DOUBLE, "y", TRUE},
 4008     {GTS_DOUBLE, "z", TRUE},
 4009     {GTS_NONE}
 4010   };
 4011 
 4012   (* GTS_OBJECT_CLASS (gfs_gl_locate_class ())->parent_class->read) (o, fp);
 4013   if (fp->type == GTS_ERROR)
 4014     return;
 4015 
 4016   var[0].data = &gl->p.x;
 4017   var[1].data = &gl->p.y;
 4018   var[2].data = &gl->p.z;
 4019   gts_file_assign_variables (fp, var);
 4020   if (fp->type == GTS_ERROR)
 4021     return;
 4022 }
 4023 
 4024 static void gl_locate_write (GtsObject * o, FILE * fp)
 4025 {
 4026   GfsGlLocate * gl = GFS_GL_LOCATE (o);
 4027 
 4028   (* GTS_OBJECT_CLASS (gfs_gl_locate_class ())->parent_class->write) (o, fp);
 4029 
 4030   fprintf (fp, " { x = %g y = %g z = %g }", gl->p.x, gl->p.y, gl->p.z);
 4031 }
 4032 
 4033 static void gl_locate_draw (GfsGl * gl, GfsFrustum * f)
 4034 {
 4035   FttVector p = GFS_GL_LOCATE (gl)->p;
 4036   gfs_simulation_map (gl->sim, &p);
 4037   FttCell * cell = gfs_domain_locate (GFS_DOMAIN (gl->sim), p, gl->maxlevel, NULL);
 4038   if (cell) {
 4039     gl->size = 0;
 4040     glMatrixMode (GL_PROJECTION);
 4041     glPushMatrix ();
 4042     glTranslatef (0., 0., gl->p->lc);
 4043     glNormal3d (0., 0., 1.);
 4044     gl_locate (cell, gl);
 4045     glPopMatrix ();
 4046   }
 4047 }
 4048 
 4049 static void gl_locate_class_init (GfsGlClass * klass)
 4050 {
 4051   GTS_OBJECT_CLASS (klass)->read = gl_locate_read;
 4052   GTS_OBJECT_CLASS (klass)->write = gl_locate_write;  
 4053   klass->draw = gl_locate_draw;
 4054 }
 4055 
 4056 GfsGlClass * gfs_gl_locate_class (void)
 4057 {
 4058   static GfsGlClass * klass = NULL;
 4059 
 4060   if (klass == NULL) {
 4061     GtsObjectClassInfo gfs_gl_locate_info = {
 4062       "GfsGlLocate",
 4063       sizeof (GfsGlLocate),
 4064       sizeof (GfsGlClass),
 4065       (GtsObjectClassInitFunc) gl_locate_class_init,
 4066       (GtsObjectInitFunc) NULL,
 4067       (GtsArgSetFunc) NULL,
 4068       (GtsArgGetFunc) NULL
 4069     };
 4070     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl_class ()), &gfs_gl_locate_info);
 4071   }
 4072 
 4073   return klass;
 4074 }
 4075 
 4076 /* GfsGlPipes: Object */
 4077 
 4078 static void gl_pipes_draw (GfsGl * gl, GfsFrustum * f)
 4079 {
 4080   glMatrixMode (GL_PROJECTION);
 4081   glPushMatrix ();
 4082   glTranslatef (0., 0., gl->p->lc);
 4083   glNormal3d (0., 0., 1.);
 4084   glColor3f (gl->lc.r, gl->lc.g, gl->lc.b);
 4085   
 4086   GSList * i = gl->sim->events->items;
 4087   while (i) {
 4088     if (GFS_IS_SOURCE_PIPE (i->data)) {
 4089       GfsSourcePipe * p = i->data;
 4090       FttVector start = p->start, end = p->end;
 4091       gfs_simulation_map (gl->sim, &start);
 4092       gfs_simulation_map (gl->sim, &end);
 4093       FttVector t = { end.x - start.x, end.y - start.y };
 4094       gdouble dt = ftt_vector_norm (&t);
 4095       if (dt > 0.) {
 4096     t.x /= dt; t.y /= dt;
 4097     glMatrixMode (GL_MODELVIEW);
 4098     glPushMatrix ();
 4099     glTranslatef (start.x, start.y, 0.);
 4100     glRotatef (atan2 (t.y, t.x)*180./M_PI, 0., 0., 1.);
 4101     /* fixme: metric is not taken into account properly */
 4102     gdouble h = p->diameter/gl->sim->physical_params.L/2.;
 4103     glBegin (GL_QUADS);
 4104     glVertex2d (dt, h);
 4105     glVertex2d (0., h);
 4106     glVertex2d (0., -h);
 4107     glVertex2d (dt, -h);
 4108     glEnd ();
 4109     if (GFS_EVENT (p)->name) {
 4110       gchar * label = GFS_EVENT (p)->name;
 4111       float bbox[6];
 4112       float size = dt/strlen(label)*gl->font_size;
 4113       gl_text_bbox (gl, label, size, bbox);
 4114       float twidth  = bbox[3] - bbox[0];
 4115       float theight = bbox[4] - bbox[1];    
 4116       gl_draw_text (gl, label, dt/2. - twidth/2., theight/2., 0., size);
 4117     }
 4118     glPopMatrix ();
 4119     glMatrixMode (GL_PROJECTION);
 4120       }
 4121     }
 4122     i = i->next;
 4123   }
 4124 
 4125   glPopMatrix ();
 4126 }
 4127 
 4128 static void gl_pipes_init (GfsGl * gl)
 4129 {
 4130   GtsColor c = { 1., 0., 0. };
 4131   gl->lc = c;
 4132   gl->font_size = 1.;
 4133 }
 4134 
 4135 static gboolean gl_pipes_relevant (GfsSimulation * sim)
 4136 {
 4137   GSList * i = sim->events->items;
 4138 
 4139   while (i) {
 4140     if (GFS_IS_SOURCE_PIPE (i->data))
 4141       return TRUE;
 4142     i = i->next;
 4143   }
 4144   return FALSE;
 4145 }
 4146 
 4147 static void gl_pipes_class_init (GfsGlClass * klass)
 4148 {
 4149   klass->draw = gl_pipes_draw;
 4150   klass->relevant = gl_pipes_relevant;
 4151 }
 4152 
 4153 GfsGlClass * gfs_gl_pipes_class (void)
 4154 {
 4155   static GfsGlClass * klass = NULL;
 4156 
 4157   if (klass == NULL) {
 4158     GtsObjectClassInfo info = {
 4159       "GfsGlPipes",
 4160       sizeof (GfsGl),
 4161       sizeof (GfsGlClass),
 4162       (GtsObjectClassInitFunc) gl_pipes_class_init,
 4163       (GtsObjectInitFunc) gl_pipes_init,
 4164       (GtsArgSetFunc) NULL,
 4165       (GtsArgGetFunc) NULL
 4166     };
 4167     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl_class ()), &info);
 4168   }
 4169 
 4170   return klass;
 4171 }
 4172 
 4173 /* GfsGlClipPlane: Object */
 4174 
 4175 static void gl_clip_plane_destroy (GtsObject * o)
 4176 {
 4177   gint i = GFS_GL_CLIP_PLANE (o)->i;
 4178 
 4179   if (i >= 0) {
 4180     glDisable (GL_CLIP_PLANE0 + i);
 4181     GFS_GL (o)->p->cp[i] = FALSE;
 4182   }
 4183 
 4184   (* GTS_OBJECT_CLASS (gfs_gl_clip_plane_class ())->parent_class->destroy) (o);
 4185 }
 4186 
 4187 static void gl_clip_plane_draw (GfsGl * gl, GfsFrustum * f)
 4188 {
 4189   if (GFS_GL_CLIP_PLANE (gl)->i >= 0) {
 4190     guint cp = GL_CLIP_PLANE0 + GFS_GL_CLIP_PLANE (gl)->i;
 4191     
 4192     if (GFS_GL_CLIP_PLANE (gl)->disabled)
 4193       glDisable (cp);
 4194     else {
 4195       GfsGl2D * gl2 = GFS_GL2D (gl);
 4196       GLdouble eq[] = {-gl2->n.x, -gl2->n.y, -gl2->n.z, gl2->pos};
 4197       
 4198       glClipPlane (cp, eq);
 4199       glEnable (cp);
 4200     }
 4201   }
 4202 }
 4203 
 4204 static void gl_clip_plane_set_simulation (GfsGl * object, GfsSimulation * sim)
 4205 {
 4206   (*GFS_GL_CLASS (GTS_OBJECT_CLASS (gfs_gl_clip_plane_class ())->parent_class)->set_simulation)
 4207     (object, sim);
 4208 
 4209   if (GFS_GL_CLIP_PLANE (object)->i < 0) {
 4210     guint i;
 4211   
 4212     g_assert (object->p);
 4213     for (i = 0; i < 6; i++)
 4214       if (!object->p->cp[i]) {
 4215     object->p->cp[i] = TRUE;
 4216     GFS_GL_CLIP_PLANE (object)->i = i;
 4217     return;
 4218       }
 4219     g_warning ("too many clipping planes!");
 4220   }
 4221 }
 4222 
 4223 static void gl_clip_plane_class_init (GfsGlClass * klass)
 4224 {
 4225   GTS_OBJECT_CLASS (klass)->destroy = gl_clip_plane_destroy;
 4226   klass->set_simulation = gl_clip_plane_set_simulation;
 4227   klass->draw = gl_clip_plane_draw;
 4228   klass->pick = NULL;
 4229 }
 4230 
 4231 static void gl_clip_plane_init (GfsGlClipPlane * cp)
 4232 {
 4233   cp->i = -1;
 4234 }
 4235 
 4236 GfsGlClass * gfs_gl_clip_plane_class (void)
 4237 {
 4238   static GfsGlClass * klass = NULL;
 4239 
 4240   if (klass == NULL) {
 4241     GtsObjectClassInfo gfs_gl_clip_plane_info = {
 4242       "GfsGlClipPlane",
 4243       sizeof (GfsGlClipPlane),
 4244       sizeof (GfsGl2DClass),
 4245       (GtsObjectClassInitFunc) gl_clip_plane_class_init,
 4246       (GtsObjectInitFunc) gl_clip_plane_init,
 4247       (GtsArgSetFunc) NULL,
 4248       (GtsArgGetFunc) NULL
 4249     };
 4250     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl2D_class ()), &gfs_gl_clip_plane_info);
 4251   }
 4252 
 4253   return klass;
 4254 }
 4255 
 4256 void gfs_gl_clip_plane_disable (GfsGlClipPlane * gl)
 4257 {
 4258   gboolean disabled = gl->disabled;
 4259   gl->disabled = TRUE;
 4260   gfs_gl_draw (GFS_GL (gl), NULL);
 4261   gl->disabled = disabled;
 4262 }
 4263 
 4264 /* GfsGlCutPlane: Object */
 4265 
 4266 static void gl_cut_plane (FttCell * cell, GfsGl * gl)
 4267 {
 4268   GList * i = GFS_GL_CUT_PLANE (gl)->list;
 4269   while (i) {
 4270     if (GFS_GL_CLASS (GTS_OBJECT (i->data)->klass)->cut)
 4271       (* GFS_GL_CLASS (GTS_OBJECT (i->data)->klass)->cut) (i->data, cell, GFS_GL2D (gl));
 4272     i = i->next;
 4273   }
 4274 }
 4275 
 4276 static void gfs_gl_cut_plane_draw (GfsGl * gl, GfsFrustum * f)
 4277 {
 4278   gl->size = 0;
 4279   glMatrixMode (GL_PROJECTION);
 4280   glPushMatrix ();
 4281   glTranslatef (0., 0., gl->p->lc);
 4282   glBegin (GL_LINES);
 4283   gfs_gl_normal (gl);
 4284   gfs_gl_cell_traverse_visible_plane (gl, f, (FttCellTraverseFunc) gl_cut_plane, gl);
 4285   glEnd ();
 4286   glPopMatrix ();
 4287 }
 4288 
 4289 static void gl_cut_plane_class_init (GfsGlClass * klass)
 4290 {
 4291   klass->draw = gfs_gl_cut_plane_draw;
 4292   klass->pick = NULL;
 4293 }
 4294 
 4295 GfsGlClass * gfs_gl_cut_plane_class (void)
 4296 {
 4297   static GfsGlClass * klass = NULL;
 4298 
 4299   if (klass == NULL) {
 4300     GtsObjectClassInfo gfs_gl_cut_plane_info = {
 4301       "GfsGlCutPlane",
 4302       sizeof (GfsGlCutPlane),
 4303       sizeof (GfsGl2DClass),
 4304       (GtsObjectClassInitFunc) gl_cut_plane_class_init,
 4305       (GtsObjectInitFunc) NULL,
 4306       (GtsArgSetFunc) NULL,
 4307       (GtsArgGetFunc) NULL
 4308     };
 4309     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl2D_class ()), &gfs_gl_cut_plane_info);
 4310   }
 4311 
 4312   return klass;
 4313 }
 4314 
 4315 /* GfsGlViewParams */
 4316 
 4317 void gfs_gl_view_params_init (GfsGlViewParams * p)
 4318 {
 4319   g_return_if_fail (p != NULL);
 4320 
 4321   p->do_init = TRUE;
 4322   p->beginx = p->beginy = 0;
 4323   p->dx = p->dy = 0;
 4324   p->tx = p->ty = 0;
 4325   p->sx = p->sy = p->sz = 1.;
 4326   p->quat[0] = p->quat[1] = p->quat[2] = 0; p->quat[3] = 1;
 4327   p->dquat[0] = p->dquat[1] = p->dquat[2] = 0; p->dquat[3] = 1;
 4328   p->fov = 30.;
 4329   gfs_gl_trackball (p->quat , 0.0, 0.0, 0.0, 0.0);
 4330 
 4331   p->bg.r = 0.3; p->bg.g = 0.4; p->bg.b = 0.6;
 4332   p->res = p->base_res = 1.;
 4333   p->lc = 0.001;
 4334   p->reactivity = 0.1;
 4335   p->motion = FALSE;
 4336 
 4337   guint i;
 4338   for (i = 0; i < 6; i++)
 4339     p->cp[i] = FALSE;
 4340 }
 4341 
 4342 GfsGlViewParams * gfs_gl_view_params_new (void)
 4343 {
 4344   GfsGlViewParams * p = g_malloc (sizeof (GfsGlViewParams));
 4345   gfs_gl_view_params_init (p);
 4346   return p;
 4347 }
 4348 
 4349 void gfs_gl_view_params_write (GfsGlViewParams * p, FILE * fp)
 4350 {
 4351   g_return_if_fail (p != NULL);
 4352   g_return_if_fail (fp != NULL);
 4353   
 4354   fprintf (fp, 
 4355        "View {\n"
 4356        "  tx = %g ty = %g\n"
 4357        "  sx = %g sy = %g sz = %g\n"
 4358        "  q0 = %g q1 = %g q2 = %g q3 = %g\n"
 4359        "  fov = %g\n"
 4360        "  r = %g g = %g b = %g\n"
 4361        "  res = %g\n"
 4362        "  lc = %g\n"
 4363        "  reactivity = %g\n"
 4364        "}",
 4365        p->tx, p->ty,
 4366        p->sx, p->sy, p->sz,
 4367        p->quat[0], p->quat[1], p->quat[2], p->quat[3], 
 4368        p->fov,
 4369        p->bg.r, p->bg.g, p->bg.b,
 4370        p->base_res,
 4371        p->lc,
 4372        p->reactivity);
 4373 }
 4374 
 4375 void gfs_gl_view_params_read (GfsGlViewParams * p, GtsFile * fp)
 4376 {
 4377   GtsFileVariable var[] = {
 4378     {GTS_FLOAT, "tx",         TRUE},
 4379     {GTS_FLOAT, "ty",         TRUE},
 4380     {GTS_FLOAT, "q0",         TRUE},
 4381     {GTS_FLOAT, "q1",         TRUE},
 4382     {GTS_FLOAT, "q2",         TRUE},
 4383     {GTS_FLOAT, "q3",         TRUE},
 4384     {GTS_FLOAT, "fov",        TRUE},
 4385     {GTS_FLOAT, "r",          TRUE},
 4386     {GTS_FLOAT, "g",          TRUE},
 4387     {GTS_FLOAT, "b",          TRUE},
 4388     {GTS_FLOAT, "res",        TRUE},
 4389     {GTS_FLOAT, "lc",         TRUE},
 4390     {GTS_FLOAT, "reactivity", TRUE},
 4391     {GTS_FLOAT, "sx",         TRUE},
 4392     {GTS_FLOAT, "sy",         TRUE},
 4393     {GTS_FLOAT, "sz",         TRUE},
 4394     {GTS_NONE}
 4395   };
 4396 
 4397   g_return_if_fail (p != NULL);
 4398   g_return_if_fail (fp != NULL);
 4399   
 4400   if (fp->type != GTS_STRING) {
 4401     gts_file_error (fp, "expecting a string (\"View\")");
 4402     return;
 4403   }
 4404   if (strcmp (fp->token->str, "View")) {
 4405     gts_file_error (fp, "unknown keyword `%s'", fp->token->str);
 4406     return;
 4407   }
 4408   gts_file_next_token (fp);
 4409 
 4410   var[0].data = &p->tx;
 4411   var[1].data = &p->ty;
 4412   var[2].data = &p->quat[0];
 4413   var[3].data = &p->quat[1];
 4414   var[4].data = &p->quat[2];
 4415   var[5].data = &p->quat[3];
 4416   var[6].data = &p->fov;
 4417   var[7].data = &p->bg.r;
 4418   var[8].data = &p->bg.g;
 4419   var[9].data = &p->bg.b;
 4420   var[10].data = &p->base_res;
 4421   var[11].data = &p->lc;
 4422   var[12].data = &p->reactivity;
 4423   var[13].data = &p->sx;
 4424   var[14].data = &p->sy;
 4425   var[15].data = &p->sz;
 4426   gts_file_assign_variables (fp, var);
 4427   if (fp->type == GTS_ERROR)
 4428     return;
 4429 }
 4430 
 4431 void gfs_gl2ps_params_init (GfsGl2PSParams * p)
 4432 {
 4433   g_return_if_fail (p != NULL);
 4434 
 4435   p->format = GFSGL_PPM;
 4436   p->options = (GL2PS_SIMPLE_LINE_OFFSET |
 4437         GL2PS_SILENT |
 4438         GL2PS_BEST_ROOT |
 4439         GL2PS_OCCLUSION_CULL |
 4440         GL2PS_USE_CURRENT_VIEWPORT |
 4441         GL2PS_TIGHT_BOUNDING_BOX);
 4442   p->width = p->height = 0;
 4443   p->sort = GL2PS_SIMPLE_SORT;
 4444   p->lw = 1.;
 4445 }
 4446 
 4447 void gfs_gl2ps_params_read (GfsGl2PSParams * p, GtsFile * fp)
 4448 {
 4449   g_return_if_fail (p != NULL);
 4450   g_return_if_fail (fp != NULL);
 4451 
 4452   gchar * format, * orientation, * sort;
 4453   GtsFileVariable var[] = {
 4454     {GTS_STRING, "format",      TRUE, &format},
 4455     {GTS_STRING, "orientation", TRUE, &orientation},
 4456     {GTS_FLOAT,  "line_width",  TRUE, &p->lw},
 4457     {GTS_UINT,   "width",       TRUE, &p->width},
 4458     {GTS_UINT,   "height",      TRUE, &p->height},
 4459     {GTS_STRING, "sort",        TRUE, &sort},
 4460     {GTS_NONE}
 4461   };
 4462 
 4463   gfs_gl2ps_params_init (p);
 4464   gts_file_assign_variables (fp, var);
 4465   if (fp->type == GTS_ERROR)
 4466     return;
 4467 
 4468   if (var[0].set) {
 4469     if (!strcmp (format, "PPM"))
 4470       p->format = GFSGL_PPM;
 4471     else if (!strcmp (format, "PS"))
 4472       p->format = GL2PS_PS;
 4473     else if (!strcmp (format, "EPS"))
 4474       p->format = GL2PS_EPS;
 4475     else if (!strcmp (format, "PDF"))
 4476       p->format = GL2PS_PDF;
 4477     else if (!strcmp (format, "SVG"))
 4478       p->format = GL2PS_SVG;
 4479     else if (!strcmp (format, "Gnuplot"))
 4480       p->format = GFSGL_GNUPLOT;
 4481     else if (!strcmp (format, "OBJ"))
 4482       p->format = GFSGL_OBJ;
 4483     else if (!strcmp (format, "KML"))
 4484       p->format = GFSGL_KML;
 4485     else if (!strcmp (format, "Latex"))
 4486       p->format = GL2PS_TEX;
 4487     else {
 4488       gts_file_variable_error (fp, var, "format", "unknown format `%s'", format);
 4489       g_free (format);
 4490       return;
 4491     }
 4492     g_free (format);
 4493   }
 4494 
 4495   if (var[1].set) {
 4496     if (!strcmp (orientation, "Portrait"))
 4497       p->options &= ~GL2PS_LANDSCAPE;
 4498     else if (!strcmp (orientation, "Landscape"))
 4499       p->options |= GL2PS_LANDSCAPE;
 4500     else {
 4501       gts_file_variable_error (fp, var, "orientation", "unknown orientation `%s'", fp->token->str);
 4502       g_free (orientation);
 4503       return;
 4504     }
 4505     g_free (orientation);
 4506   }
 4507   else
 4508     p->options &= ~GL2PS_LANDSCAPE;  
 4509 
 4510   if (var[5].set) {
 4511     if (!strcmp (sort, "BSP"))
 4512       p->sort = GL2PS_BSP_SORT;
 4513     else if (!strcmp (sort, "Simple"))
 4514       p->sort = GL2PS_SIMPLE_SORT;
 4515     else if (!strcmp (sort, "None"))
 4516       p->sort = GL2PS_NO_SORT;
 4517     else {
 4518       gts_file_variable_error (fp, var, "sort", "unknown option `%s'", sort);
 4519       g_free (sort);
 4520       return;
 4521     }
 4522     g_free (sort);
 4523   }
 4524 }
 4525 
 4526 void gfs_gl2ps_params_write (GfsGl2PSParams * p, FILE * fp)
 4527 {
 4528   g_return_if_fail (p != NULL);
 4529   g_return_if_fail (fp != NULL);
 4530 
 4531   fprintf (fp, " { format = %s orientation = %s line_width = %g width = %d height = %d sort = %s }",
 4532        p->format == GFSGL_PPM ?     "PPM" :
 4533        p->format == GL2PS_PS ?      "PS" :
 4534        p->format == GL2PS_EPS ?     "EPS" :
 4535        p->format == GL2PS_PDF ?     "PDF" :
 4536        p->format == GL2PS_SVG ?     "SVG" :
 4537        p->format == GFSGL_GNUPLOT ? "Gnuplot" :
 4538        p->format == GFSGL_OBJ ?     "OBJ" :
 4539        p->format == GFSGL_KML ?     "KML" :
 4540        p->format == GL2PS_TEX ?     "Latex" : 
 4541        "?",
 4542        p->options & GL2PS_LANDSCAPE ? "Landscape" : "Portrait",
 4543        p->lw, p->width, p->height,
 4544        p->sort == GL2PS_SIMPLE_SORT ? "Simple" :
 4545        p->sort == GL2PS_BSP_SORT ?    "BSP" :
 4546        p->sort == GL2PS_NO_SORT ?     "None" :
 4547        "?");
 4548 }
 4549 
 4550 void gfs_gl_write_image (FILE * fp, const GLubyte * buffer, guint width, guint height)
 4551 {
 4552   gint i, x, y;
 4553   const GLubyte *ptr = buffer;
 4554 
 4555   g_return_if_fail (fp != NULL);
 4556   g_return_if_fail (buffer != NULL);
 4557 
 4558   fprintf (fp, "P6 %d %d 255\n", width, height);
 4559   for (y = height - 1; y >= 0; y--) {
 4560     for (x = 0; x < width; x++) {
 4561       i = (y*width + x)*4;
 4562       fputc (ptr[i], fp);   /* write red */
 4563       fputc (ptr[i+1], fp); /* write green */
 4564       fputc (ptr[i+2], fp); /* write blue */
 4565     }
 4566   }
 4567 }
 4568 
 4569 typedef struct {
 4570   GList * symmetries;
 4571   FttVector * s;
 4572   gdouble max;
 4573 } ExtentData;
 4574 
 4575 static void update_max (FttVector * p, ExtentData * d)
 4576 {
 4577   guint i, n = create_symmetries (d->s, d->symmetries, p);
 4578   for (i = 0; i < n; i++) {
 4579     FttComponent c;
 4580     gdouble m = 0.;
 4581     for (c = 0; c < FTT_DIMENSION; c++)
 4582       m += (&d->s[i].x)[c]*(&d->s[i].x)[c];
 4583     if (m > d->max) d->max = m;
 4584   }
 4585 }
 4586 
 4587 static void extent (GfsBox * b, ExtentData * d)
 4588 {
 4589   FttCell * cell = b->root;
 4590   gdouble h = ftt_cell_size (cell)/2.;
 4591   FttVector p, o;
 4592 
 4593   ftt_cell_pos (cell, &p);
 4594 #if FTT_2D
 4595   o.x = p.x + h; o.y = p.y + h; update_max (&o, d);
 4596   o.x = p.x + h; o.y = p.y - h; update_max (&o, d);
 4597   o.x = p.x - h; o.y = p.y + h; update_max (&o, d);
 4598   o.x = p.x - h; o.y = p.y - h; update_max (&o, d);
 4599 #else /* 3D */
 4600   o.x = p.x + h; o.y = p.y + h; o.z = p.z + h; update_max (&o, d);
 4601   o.x = p.x + h; o.y = p.y - h; o.z = p.z + h; update_max (&o, d);
 4602   o.x = p.x - h; o.y = p.y + h; o.z = p.z + h; update_max (&o, d);
 4603   o.x = p.x - h; o.y = p.y - h; o.z = p.z + h; update_max (&o, d);
 4604   o.x = p.x + h; o.y = p.y + h; o.z = p.z - h; update_max (&o, d);
 4605   o.x = p.x + h; o.y = p.y - h; o.z = p.z - h; update_max (&o, d);
 4606   o.x = p.x - h; o.y = p.y + h; o.z = p.z - h; update_max (&o, d);
 4607   o.x = p.x - h; o.y = p.y - h; o.z = p.z - h; update_max (&o, d);
 4608 #endif /* 3D */
 4609 }
 4610 
 4611 /**
 4612  * gfs_gl_domain_extent:
 4613  * @domain: a #GfsDomain.
 4614  * @symmetries: a list of #GfsGlSymmetry.
 4615  *
 4616  * Returns: the maximum extent along any dimension of @domain.
 4617  */
 4618 gdouble gfs_gl_domain_extent (GfsDomain * domain, GList * symmetries)
 4619 {
 4620   if (domain == NULL)
 4621     return 1.;
 4622 
 4623   ExtentData d;
 4624   d.max = 0.;
 4625   d.symmetries = symmetries;
 4626   guint n = 1;
 4627   while (symmetries) {
 4628     n *= 2;
 4629     symmetries = symmetries->next;
 4630   }
 4631   d.s = g_malloc (n*sizeof (FttVector));
 4632   gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) extent, &d);
 4633   g_free (d.s);
 4634   gfs_all_reduce (domain, d.max, MPI_DOUBLE, MPI_MAX);
 4635   return 2.*sqrt (d.max);
 4636 }
 4637 
 4638 /**
 4639  * gfs_gl_feedback_begin:
 4640  * @buffersize: the size of the feedback buffer.
 4641  *
 4642  * Returns: a newly setup openGL feedback buffer.
 4643  */
 4644 GfsGlFeedback * gfs_gl_feedback_begin (guint buffersize)
 4645 {
 4646   GfsGlFeedback * f;
 4647 
 4648   g_return_val_if_fail (buffersize > 0, NULL);
 4649 
 4650   f = g_malloc (sizeof (GfsGlFeedback));
 4651   f->feedback = g_malloc (sizeof (GLfloat)*buffersize);
 4652   glFeedbackBuffer (buffersize, GL_3D_COLOR, f->feedback);
 4653   glRenderMode (GL_FEEDBACK);
 4654   return f;
 4655 }
 4656 
 4657 /**
 4658  * gfs_gl_feedback_end:
 4659  * @f: the #GfsGlFeedback.
 4660  * @sim: a #GfsSimulation.
 4661  * @fp: a file pointer.
 4662  * @format: the file format.
 4663  *
 4664  * Writes to @fp the contents of @f.
 4665  *
 4666  * Returns: %FALSE if @f has overflowed, %TRUE otherwise.
 4667  */
 4668 gboolean gfs_gl_feedback_end (GfsGlFeedback * f, GfsSimulation * sim, FILE * fp, GfsGlFormat format)
 4669 {
 4670   GLint size;
 4671 
 4672   g_return_val_if_fail (f != NULL, FALSE);
 4673   g_return_val_if_fail (sim != NULL, FALSE);
 4674   g_return_val_if_fail (fp != NULL, FALSE);
 4675 
 4676   size = glRenderMode (GL_RENDER);
 4677   if (size >= 0) {
 4678     GLint used = size;
 4679     GLfloat * current = f->feedback;
 4680     GLdouble model[16], proj[16];
 4681     GLint viewport[4];
 4682     FttVector p, lastp = { G_MAXDOUBLE, G_MAXDOUBLE, G_MAXDOUBLE };
 4683     guint nps = 0, np = 0, linetoken = FALSE;
 4684 
 4685     /* Header */
 4686     switch (format) {
 4687     case GFSGL_KML:
 4688       fputs ("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n"
 4689          "<kml xmlns=\"http://www.opengis.net/kml/2.2\">\n"
 4690          "<Placemark>\n"
 4691          "<name>GfsView</name>\n"
 4692          "<MultiGeometry>\n", fp);
 4693       break;
 4694     default:
 4695       break;
 4696     }
 4697 
 4698     /* Body */
 4699     glGetDoublev (GL_MODELVIEW_MATRIX, model);
 4700     glGetDoublev (GL_PROJECTION_MATRIX, proj);
 4701     glGetIntegerv (GL_VIEWPORT, viewport);
 4702     while (used > 0)
 4703       switch ((GLint) *current) {
 4704 
 4705       case GL_POINT_TOKEN :
 4706     current++; used--;
 4707     gluUnProject (current[0], current[1], current[2], model, proj, viewport, &p.x, &p.y, &p.z);
 4708     gfs_simulation_map_inverse (sim, &p);
 4709     switch (format) {
 4710     case GFSGL_GNUPLOT:
 4711       fprintf (fp, "%g %g %g\n\n", p.x, p.y, p.z); break;
 4712     case GFSGL_OBJ:
 4713       fprintf (fp, "v %g %g %g\np -1\n", p.x, p.y, p.z); 
 4714       np++;
 4715       nps = 0;
 4716       break;
 4717     case GFSGL_KML:
 4718       if (linetoken) {
 4719         fputs ("</coordinates>\n</LineString>\n", fp);
 4720         linetoken = FALSE;
 4721       }
 4722       fprintf (fp, "<Point>\n<coordinates>%f,%f,%f</coordinates>\n</Point>\n", p.x, p.y, p.z);
 4723       break;
 4724     default:
 4725       g_assert_not_reached ();
 4726     }
 4727     current += 7; used -= 7;
 4728     break;
 4729 
 4730       case GL_LINE_TOKEN :
 4731       case GL_LINE_RESET_TOKEN :
 4732     current++; used--;
 4733     gluUnProject (current[0], current[1], current[2], model, proj, viewport, &p.x, &p.y, &p.z);
 4734     gfs_simulation_map_inverse (sim, &p);
 4735     if (p.x != lastp.x || p.y != lastp.y || p.z != lastp.z) {
 4736       switch (format) {
 4737       case GFSGL_GNUPLOT:
 4738         fprintf (fp, "\n%g %g %g\n", p.x, p.y, p.z);
 4739         break;
 4740       case GFSGL_OBJ:
 4741         if (nps > 0) {
 4742           fputc ('l', fp);
 4743           int i;
 4744           for (i = nps; i <= np; i++)
 4745         fprintf (fp, " %d", i);
 4746           fputc ('\n', fp);
 4747         }
 4748         nps = ++np;
 4749         fprintf (fp, "v %g %g %g\n", p.x, p.y, p.z);
 4750         break;
 4751       case GFSGL_KML:
 4752         if (linetoken)
 4753           fputs ("</coordinates>\n</LineString>\n", fp);
 4754         fprintf (fp, "<LineString>\n<coordinates>\n%f,%f,%f\n", p.x, p.y, p.z);
 4755         linetoken = TRUE;
 4756         break;
 4757       default:
 4758         g_assert_not_reached ();
 4759       }
 4760     }
 4761     current += 7; used -= 7;
 4762     gluUnProject (current[0], current[1], current[2], model, proj, viewport, &p.x, &p.y, &p.z);
 4763     gfs_simulation_map_inverse (sim, &p);
 4764     lastp = p;
 4765     switch (format) {
 4766     case GFSGL_GNUPLOT:
 4767       fprintf (fp, "%g %g %g\n", p.x, p.y, p.z); 
 4768       break;
 4769     case GFSGL_OBJ:
 4770       fprintf (fp, "v %g %g %g\n", p.x, p.y, p.z); 
 4771       np++;
 4772       break;
 4773     case GFSGL_KML:
 4774       fprintf (fp, "%f,%f,%f\n", p.x, p.y, p.z);
 4775       break;
 4776     default:
 4777       g_assert_not_reached ();
 4778     }
 4779     current += 7; used -= 7;
 4780     break;
 4781 
 4782       case GL_POLYGON_TOKEN : {
 4783     GLint count = (GLint) current[1], vcount = 0;
 4784     current += 2;
 4785     used -= 2;
 4786     if (format == GFSGL_KML) {
 4787