"Fossies" - the Fresh Open Source Software Archive

Member "gfsview-snapshot-121130/gl/iso.c" (30 Nov 2012, 18222 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 "iso.c" see the Fossies "Dox" file reference documentation.

    1 /* Gerris - The GNU Flow Solver
    2  * Copyright (C) 2001-2004 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 
   24 /* GfsGlIsosurface: Object */
   25 
   26 typedef struct {
   27   guint m;
   28   FttVector * v;
   29   FttVector * n;
   30   gdouble * sv;
   31   gboolean coarse_fine;
   32 } polygon;
   33 
   34 typedef struct {
   35   polygon * p;
   36   guint n;
   37 } polygons;
   38 
   39 #define param(v) (gls->max > gls->min ? ((v) - gls->min)/(gls->max - gls->min) : 0.5)
   40 #define color(v) (gfs_colormap_color (gls->cmap, param (v)))
   41 
   42 static void polygon_draw (polygon * p, GfsGlIsosurface * gl)
   43 {
   44   guint i;
   45 
   46   glBegin (GL_POLYGON);
   47   for (i = 0; i < p->m; i++) {
   48     if (gl->use_scalar) {
   49       GfsGlScalar * gls = GFS_GL_SCALAR (gl);
   50 
   51       if (gfs_gl_vector_format (GFS_GL (gl))) {
   52     GtsColor c = color (p->sv[i]);
   53     glColor3f (c.r, c.g, c.b);
   54       }
   55       else
   56     glTexCoord1d (param (p->sv[i]));
   57     }
   58     if (gl->reversed)
   59       glNormal3d (-p->n[i].x, -p->n[i].y, -p->n[i].z);
   60     else
   61       glNormal3d (p->n[i].x, p->n[i].y, p->n[i].z);
   62     glVertex3d (p->v[i].x, p->v[i].y, p->v[i].z);
   63   }
   64   glEnd ();
   65 }
   66 
   67 static void polygons_draw (polygons * p, GfsGlIsosurface * gl)
   68 {
   69   guint i;
   70 
   71   for (i = 0; i < p->n; i++)
   72     polygon_draw (&p->p[i], gl);
   73 }
   74 
   75 static void polygons_destroy (polygons * p)
   76 {
   77   guint i;
   78   for (i = 0; i < p->n; i++) {
   79     g_free (p->p[i].v);
   80     g_free (p->p[i].n);
   81     g_free (p->p[i].sv);
   82   }
   83   g_free (p->p);
   84   g_free (p);
   85 }
   86 
   87 static polygons * polygons_add (polygons * p, guint n)
   88 {
   89   if (p == NULL) {
   90     p = g_malloc (sizeof (polygons));
   91     p->p = g_malloc (sizeof (polygon));
   92     p->n = 0;
   93   }
   94   else
   95     p->p = g_realloc (p->p, sizeof (polygon)*(p->n + 1));
   96   p->p[p->n].v = g_malloc (n*sizeof (FttVector)); 
   97   p->p[p->n].n = g_malloc (n*sizeof (FttVector));
   98   p->p[p->n].sv = g_malloc (n*sizeof (gdouble));
   99   p->p[p->n].m = n;
  100   p->n++;
  101   return p;
  102 }
  103 
  104 static void free_polygons (FttCell * cell, GfsVariable * p)
  105 {
  106   gpointer po = GFS_DOUBLE_TO_POINTER (GFS_VALUE (cell, p));
  107 
  108   if (po) {
  109     polygons_destroy (po);
  110     GFS_VALUE (cell, p) = 0.;
  111   }
  112 }
  113 
  114 void gfs_gl_isosurface_reset (GfsGlIsosurface * gl)
  115 {
  116   g_return_if_fail (gl != NULL);
  117   
  118   if (gl->p && GFS_GL (gl)->sim)
  119     gfs_domain_cell_traverse (GFS_DOMAIN (GFS_GL (gl)->sim), 
  120                   FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
  121                   (FttCellTraverseFunc) free_polygons, gl->p);
  122 }
  123 
  124 static void gl_isosurface_destroy (GtsObject * o)
  125 {
  126   GfsGlIsosurface * gl = GFS_GL_ISOSURFACE (o);
  127 
  128   gfs_gl_var_func_destroy (gl->vf);
  129   g_string_free (gl->expr, TRUE);
  130 
  131   if (gl->min)
  132     gts_object_destroy (GTS_OBJECT (gl->min));
  133   if (gl->max)
  134     gts_object_destroy (GTS_OBJECT (gl->max));
  135 
  136   if (gl->p) {
  137     gfs_gl_isosurface_reset (gl);
  138     gts_object_destroy (GTS_OBJECT (gl->p));
  139   }
  140 
  141   (* GTS_OBJECT_CLASS (gfs_gl_isosurface_class ())->parent_class->destroy) (o);
  142 }
  143 
  144 static void gl_isosurface_read (GtsObject ** o, GtsFile * fp)
  145 {
  146   GfsGlIsosurface * gl = GFS_GL_ISOSURFACE (*o);
  147   GtsFileVariable var[] = {
  148     {GTS_DOUBLE, "level",      TRUE},
  149     {GTS_INT,    "reversed",   TRUE},
  150     {GTS_INT,    "use_scalar", TRUE},
  151     {GTS_NONE}
  152   };
  153 
  154   (* GTS_OBJECT_CLASS (gfs_gl_isosurface_class ())->parent_class->read) (o, fp);
  155   if (fp->type == GTS_ERROR)
  156     return;
  157 
  158   g_string_free (gl->expr, TRUE);
  159   if (!(gl->expr = gfs_function_expression (fp, NULL)))
  160     return;
  161   gts_file_next_token (fp);
  162 
  163   var[0].data = &gl->level;
  164   var[1].data = &gl->reversed;
  165   var[2].data = &gl->use_scalar;
  166   gts_file_assign_variables (fp, var);
  167 }
  168 
  169 static void gl_isosurface_write (GtsObject * o, FILE * fp)
  170 {
  171   GfsGlIsosurface * gl = GFS_GL_ISOSURFACE (o);
  172 
  173   (* GTS_OBJECT_CLASS (gfs_gl_isosurface_class ())->parent_class->write) (o, fp);
  174 
  175   fprintf (fp, " %s {\n"
  176        "  level = %g\n"
  177        "  reversed = %d\n"
  178        "  use_scalar = %d\n"
  179        "}",
  180        gl->expr->str,
  181        gl->level,
  182        gl->reversed, 
  183        (gl->use_scalar != NULL));
  184 }
  185 
  186 static void interpolated_normal (FttVector * n, GfsInterpolator * inter, GfsVariable * v)
  187 {
  188   if (n->x == G_MAXDOUBLE) {
  189     guint i;
  190     
  191     n->x = n->y = n->z = 0.;
  192     for (i = 0; i < inter->n; i++) {
  193       gdouble w = inter->w[i];
  194       FttCell * c = inter->c[i];
  195       
  196       n->x += w*gfs_center_gradient (c, FTT_X, v->i);
  197       n->y += w*gfs_center_gradient (c, FTT_Y, v->i);
  198       n->z += w*gfs_center_gradient (c, FTT_Z, v->i);
  199     }
  200   }
  201 }
  202 
  203 static void cube_iso_intersection (FttCell * cell,
  204                    FttVector p[12],
  205                    FttVector n[12],
  206                    gint orient[12],
  207                    gdouble sv[12],
  208                    gdouble vc[8],
  209                    gdouble svc[8],
  210                    FttVector nc[8],
  211                    GfsInterpolator inter[8],
  212                    GfsVariable * var,
  213                    gdouble val,
  214                    GfsVariable * svar,
  215                    gint max_level)
  216 {
  217   FttVector o;
  218   gdouble h = ftt_cell_size (cell)*SLIGHTLY_LARGER;
  219   guint i;
  220 
  221   for (i = 0; i < 8; i++) {
  222     guint j;
  223 
  224     gfs_cell_corner_interpolator (cell, corner[i], var->centered, max_level, &inter[i]);
  225     vc[i] = svc[i] = 0.;
  226     for (j = 0; j < inter[i].n; j++) {
  227       vc[i] += inter[i].w[j]*GFS_VALUE (inter[i].c[j], var);
  228       if (svar)
  229     svc[i] += inter[i].w[j]*GFS_VALUE (inter[i].c[j], svar);
  230     }
  231     nc[i].x = G_MAXDOUBLE;
  232   }
  233 
  234   ftt_cell_pos (cell, &o);
  235   o.x -= h/2.; o.y -= h/2.; o.z -= h/2.;
  236   for (i = 0; i < 12; i++) {
  237     guint j = edge1[i][0], k = edge1[i][1];
  238     if ((vc[j] >= val && vc[k] < val) || (vc[j] < val && vc[k] >= val)) {
  239       gdouble a = (val - vc[j])/(vc[k] - vc[j]);
  240       FttVector e, d;
  241 
  242       d.x = o.x + h*edge[i][0].x; d.y = o.y + h*edge[i][0].y; d.z = o.z + h*edge[i][0].z;
  243       e.x = o.x + h*edge[i][1].x; e.y = o.y + h*edge[i][1].y; e.z = o.z + h*edge[i][1].z;
  244       p[i].x = d.x + a*(e.x - d.x); p[i].y = d.y + a*(e.y - d.y); p[i].z = d.z + a*(e.z - d.z);
  245       sv[i] = svc[j] + a*(svc[k] - svc[j]);
  246 
  247       interpolated_normal (&nc[j], &inter[j], var);
  248       interpolated_normal (&nc[k], &inter[k], var);
  249       n[i].x = nc[j].x + a*(nc[k].x - nc[j].x);
  250       n[i].y = nc[j].y + a*(nc[k].y - nc[j].y); 
  251       n[i].z = nc[j].z + a*(nc[k].z - nc[j].z);
  252       orient[i] = (vc[j] >= val);
  253     }
  254     else
  255       orient[i] = -1;
  256   }
  257 }
  258 
  259 static gboolean face_is_coarse_fine (FttCell * cell, guint e, gint orient, gint max_depth)
  260 {
  261   FttCell * n = ftt_cell_neighbor (cell, connect[e][orient][3]);
  262 
  263   if (n && !FTT_CELL_IS_LEAF (n)) {
  264     guint l = ftt_cell_level (n);
  265     return (l != max_depth && l == ftt_cell_level (cell));
  266   }
  267   return FALSE;
  268 }
  269 
  270 static gint next_vertex (guint s, 
  271              gdouble v0, gdouble v1, gdouble v2, gdouble v3, 
  272              FttVector p[4], FttVector n[4], gdouble sv[4],
  273              gdouble val,
  274              FttVector * p1, FttVector * n1, gdouble * sv1)
  275 {
  276   gdouble v[4];
  277   guint i;
  278   v[0] = v0; v[1] = v1; v[2] = v2; v[3] = v3;
  279   for (i = 1; i < 4; i++) {
  280     guint k1 = (s + i) % 4;
  281     guint k2 = (s + i + 1) % 4;
  282     if ((v[k1] >= val && v[k2] < val) || (v[k1] < val && v[k2] >= val)) {
  283       gdouble a = (val - v[k1])/(v[k2] - v[k1]);
  284       p1->x = p[k1].x + a*(p[k2].x - p[k1].x);
  285       p1->y = p[k1].y + a*(p[k2].y - p[k1].y);
  286       p1->z = p[k1].z + a*(p[k2].z - p[k1].z);
  287       n1->x = n[k1].x + a*(n[k2].x - n[k1].x);
  288       n1->y = n[k1].y + a*(n[k2].y - n[k1].y);
  289       n1->z = n[k1].z + a*(n[k2].z - n[k1].z);
  290       *sv1 = sv[k1] + a*(sv[k2] - sv[k1]);
  291       return k1;
  292     }
  293   }
  294   return -1;
  295 }
  296 
  297 /* see doc/coarse_fine.fig */
  298 static gint next_square (guint square, guint edge, 
  299              gdouble v[9], gdouble sv[9], gdouble val,
  300              FttVector p[9], FttVector n[9],
  301              FttVector v3[24], FttVector n3[24], gdouble sv3[24],
  302              guint * nv)
  303 {
  304   static guint sqv[4][4] = {
  305     {3, 4, 8, 7}, {4, 0, 5, 8}, {7, 8, 6, 2}, {8, 5, 1, 6}
  306   };
  307   static gint sqn[4][4] = {
  308     {-1, 1, 2, -4}, {-1, -2, 3, 0}, {0, 3, -3, -4}, {1, -2, -3, 2}
  309   };
  310   static gint edn[4][4] = {
  311     {0, 3, 0, 0}, {0, 0, 0, 1}, {2, 3, 0, 0}, {2, 0, 0, 1}
  312   };
  313   FttVector p2[4], n2[4], p1, n1;  
  314   gdouble sv2[4], sv1;
  315   gint u, nsquare, nedge;
  316 
  317   p2[0] = p[sqv[square][0]];
  318   p2[1] = p[sqv[square][1]];
  319   p2[2] = p[sqv[square][2]];
  320   p2[3] = p[sqv[square][3]];
  321 
  322   n2[0] = n[sqv[square][0]];
  323   n2[1] = n[sqv[square][1]];
  324   n2[2] = n[sqv[square][2]];
  325   n2[3] = n[sqv[square][3]];
  326 
  327   sv2[0] = sv[sqv[square][0]];
  328   sv2[1] = sv[sqv[square][1]];
  329   sv2[2] = sv[sqv[square][2]];
  330   sv2[3] = sv[sqv[square][3]];
  331 
  332   u = next_vertex (edge,
  333            v[sqv[square][0]], v[sqv[square][1]],
  334            v[sqv[square][2]], v[sqv[square][3]], p2, n2, sv2, val, &p1, &n1, &sv1);
  335   if (u < 0) {
  336     *nv = 0;
  337     return -1;
  338   }
  339   nsquare = sqn[square][u];
  340   if (nsquare < 0)
  341     return - nsquare - 1;
  342   n3[*nv] = n1;
  343   sv3[*nv] = sv1;
  344   v3[(*nv)++] = p1;
  345   g_assert (*nv <= 24); 
  346   nedge = edn[square][u];
  347   return next_square (nsquare, nedge, v, sv, val, p, n, v3, n3, sv3, nv);
  348 }
  349 
  350 static void gl_isosurface (FttCell * cell, GfsGl * gl)
  351 {
  352   GfsGlIsosurface * gls = GFS_GL_ISOSURFACE (gl);
  353   polygons * p = GFS_DOUBLE_TO_POINTER (GFS_VALUE (cell, gls->p));
  354 
  355   if (p) {
  356     polygons_draw (p, gls);
  357     gl->size++;
  358   }
  359   else {
  360     gboolean coarse_fine = FALSE;
  361     
  362     {
  363       FttVector _p[12], _n[12], v[24], n[24];
  364       gint _orient[12];
  365       gdouble sv[24], _sv[12], _svc[8], _vc[8], _h = ftt_cell_size (cell)*SLIGHTLY_LARGER;
  366       FttVector _nc[8];
  367       GfsInterpolator _inter[8];
  368       guint _i;
  369       cube_iso_intersection (cell, _p, _n, _orient, _sv, _vc, _svc, _nc, _inter, 
  370                  gls->v, gls->level,
  371                  gls->use_scalar ? GFS_GL_SCALAR (gls)->v : NULL,
  372                  GFS_GL (gl)->maxlevel);
  373       for (_i = 0; _i < 12; _i++) {
  374     guint nv = 0, _e = _i;
  375     while (_orient[_e] >= 0) {
  376       guint _m = 0, * _ne = connect[_e][_orient[_e]];
  377       n[nv] = _n[_e];
  378       sv[nv] = _sv[_e];
  379       v[nv++] = _p[_e];
  380       if (face_is_coarse_fine (cell, _e, _orient[_e], gl->maxlevel)) {
  381         guint _v0 = connectv[_e][_orient[_e]][0];
  382         guint _v1 = connectv[_e][_orient[_e]][1];
  383         guint _v2 = connectv[_e][_orient[_e]][2];
  384         guint _v3 = connectv[_e][_orient[_e]][3];
  385         FttVector _o, _p2[9], _n2[9];
  386         gdouble _v[9], _vs[9];
  387         guint _q, j;
  388 
  389         _v[0] = _vc[_v0]; _v[1] = _vc[_v1]; _v[2] = _vc[_v2]; _v[3] = _vc[_v3];
  390         _v[4] = (_v[3] + _v[0])/2.;
  391         _v[5] = (_v[1] + _v[0])/2.;
  392         _v[6] = (_v[1] + _v[2])/2.;
  393         _v[7] = (_v[3] + _v[2])/2.;
  394         _v[8] = (_v[0] + _v[1] + _v[2] + _v[3])/4.;
  395 
  396         _vs[0] = _svc[_v0]; _vs[1] = _svc[_v1]; _vs[2] = _svc[_v2]; _vs[3] = _svc[_v3];
  397         _vs[4] = (_vs[3] + _vs[0])/2.;
  398         _vs[5] = (_vs[1] + _vs[0])/2.;
  399         _vs[6] = (_vs[1] + _vs[2])/2.;
  400         _vs[7] = (_vs[3] + _vs[2])/2.;
  401         _vs[8] = (_vs[0] + _vs[1] + _vs[2] + _vs[3])/4.;
  402 
  403         interpolated_normal (&_nc[_v0], &_inter[_v0], gls->v);
  404         _n2[0] = _nc[_v0];
  405         interpolated_normal (&_nc[_v1], &_inter[_v1], gls->v);
  406         _n2[1] = _nc[_v1];
  407         interpolated_normal (&_nc[_v2], &_inter[_v2], gls->v);
  408         _n2[2] = _nc[_v2];
  409         interpolated_normal (&_nc[_v3], &_inter[_v3], gls->v);
  410         _n2[3] = _nc[_v3];
  411         ftt_cell_pos (cell, &_o);
  412         _o.x -= _h/2.; _o.y -= _h/2.; _o.z -= _h/2;
  413         for (j = 0; j < 3; j++) {
  414           (&_p2[0].x)[j] = (&_o.x)[j] + _h*(&cvertex[_v0].x)[j];
  415           (&_p2[1].x)[j] = (&_o.x)[j] + _h*(&cvertex[_v1].x)[j];
  416           (&_p2[2].x)[j] = (&_o.x)[j] + _h*(&cvertex[_v2].x)[j];
  417           (&_p2[3].x)[j] = (&_o.x)[j] + _h*(&cvertex[_v3].x)[j];
  418           (&_p2[4].x)[j] = ((&_p2[3].x)[j] + (&_p2[0].x)[j])/2.;
  419           (&_p2[5].x)[j] = ((&_p2[1].x)[j] + (&_p2[0].x)[j])/2.;
  420           (&_p2[6].x)[j] = ((&_p2[1].x)[j] + (&_p2[2].x)[j])/2.;
  421           (&_p2[7].x)[j] = ((&_p2[3].x)[j] + (&_p2[2].x)[j])/2.;
  422           (&_p2[8].x)[j] = ((&_p2[0].x)[j] + (&_p2[1].x)[j] + 
  423                 (&_p2[2].x)[j] + (&_p2[3].x)[j])/4.;
  424           (&_n2[4].x)[j] = ((&_n2[3].x)[j] + (&_n2[0].x)[j])/2.;
  425           (&_n2[5].x)[j] = ((&_n2[1].x)[j] + (&_n2[0].x)[j])/2.;
  426           (&_n2[6].x)[j] = ((&_n2[1].x)[j] + (&_n2[2].x)[j])/2.;
  427           (&_n2[7].x)[j] = ((&_n2[3].x)[j] + (&_n2[2].x)[j])/2.;
  428           (&_n2[8].x)[j] = ((&_n2[0].x)[j] + (&_n2[1].x)[j] + 
  429                 (&_n2[2].x)[j] + (&_n2[3].x)[j])/4.;
  430         }
  431         _q = next_square ((gls->level - _v[3])/(_v[0] - _v[3]) > 0.5, 0, _v, _vs, gls->level,
  432                   _p2, _n2, v, n, sv, &nv);
  433         if (_q <= 0)
  434           break;
  435         _orient[_e] = -1;
  436         _e = _ne[_q - 1];
  437         coarse_fine = FALSE;
  438       }
  439       else {
  440         _orient[_e] = -1;
  441         while (_m < 3 && _orient[_e] < 0)
  442           _e = _ne[_m++];
  443       }
  444     }
  445     if (nv > 2) {
  446       guint i;
  447       polygon * q;
  448 
  449       p = polygons_add (p, nv);
  450       q = &p->p[p->n - 1];
  451 
  452       for (i = 0; i < nv; i++) {
  453         q->v[i] = v[i];
  454         q->n[i] = n[i];
  455         q->sv[i] = sv[i];
  456       }
  457       q->coarse_fine = coarse_fine;
  458     }
  459     coarse_fine = FALSE;
  460       }
  461     }
  462     if (p) {
  463       polygons_draw (p, gls);
  464       GFS_DOUBLE_TO_POINTER (GFS_VALUE (cell, gls->p)) = p;
  465       gl->size++;
  466     }
  467   }
  468 }
  469 
  470 static void gl_isosurface_draw (GfsGl * gl, GfsFrustum * f)
  471 {
  472   GfsGlIsosurface * gls = GFS_GL_ISOSURFACE (gl);
  473   
  474   if (gls->use_scalar && gls->use_scalar != GFS_GL_SCALAR (gl)->v) {
  475     gls->use_scalar = GFS_GL_SCALAR (gl)->v;
  476     gfs_gl_isosurface_reset (gls);
  477   }
  478 
  479   gl->size = 0;
  480   glShadeModel (GL_SMOOTH);
  481   if (gls->use_scalar && !gfs_gl_vector_format (gl)) {
  482     glEnable (GL_TEXTURE_1D);
  483     gfs_colormap_texture (GFS_GL_SCALAR (gl)->cmap);
  484     glColor3f (1., 1., 1.);
  485     gfs_gl_cell_traverse_visible_iso (gl, f, gls->min, gls->max, gls->level,
  486                       (FttCellTraverseFunc) gl_isosurface, gl);
  487     glDisable (GL_TEXTURE_1D);
  488   }
  489   else
  490     gfs_gl_cell_traverse_visible_iso (gl, f, gls->min, gls->max, gls->level,
  491                       (FttCellTraverseFunc) gl_isosurface, gl);
  492 
  493   if (gls->use_scalar)
  494     (* GFS_GL_CLASS (GTS_OBJECT (gl)->klass->parent_class)->draw) (gl, f);  
  495 }
  496 
  497 static void min_max_iso (FttCell * cell, GfsGlIsosurface * gl)
  498 {
  499   gdouble min = G_MAXDOUBLE, max = -G_MAXDOUBLE;
  500   guint i;
  501   gpointer p;
  502 
  503   if (FTT_CELL_IS_LEAF (cell))
  504     for (i = 0; i < 8; i++) {
  505       gdouble v = gfs_cell_corner_value (cell, corner[i], gl->v, -1);
  506       if (v < min) min = v;
  507       if (v > max) max = v;
  508     }
  509   else {
  510     FttCellChildren child;
  511 
  512     ftt_cell_children (cell, &child);
  513     for (i = 0; i < FTT_CELLS; i++)
  514       if (child.c[i]) {
  515     gdouble vmin = GFS_VALUE (child.c[i], gl->min);
  516     gdouble vmax = GFS_VALUE (child.c[i], gl->max);
  517     if (vmin < min) min = vmin;
  518     if (vmax > max) max = vmax;
  519       }
  520   }
  521   GFS_VALUE (cell, gl->min) = min;
  522   GFS_VALUE (cell, gl->max) = max;
  523   if ((p = GFS_DOUBLE_TO_POINTER (GFS_VALUE (cell, gl->p)))) {
  524     polygons_destroy (p);
  525     GFS_VALUE (cell, gl->p) = 0.;
  526   }
  527 }
  528 
  529 static void min_max_iso1 (GfsBox * b, GfsGlIsosurface * gl)
  530 {
  531   gdouble min = GFS_VALUE (b->root, gl->min);
  532   gdouble max = GFS_VALUE (b->root, gl->max);
  533   if (min < gl->minv) gl->minv = min;
  534   if (max > gl->maxv) gl->maxv = max;
  535 }
  536 
  537 GtsFile * gfs_gl_isosurface_set (GfsGlIsosurface * gl, const gchar * func)
  538 {
  539   GtsFile * fp;
  540 
  541   g_return_val_if_fail (gl != NULL, NULL);
  542   g_return_val_if_fail (func != NULL, NULL);
  543 
  544   if ((fp = gfs_gl_var_func_set (gl->vf, GFS_GL (gl)->sim, func, gl->expr, NULL)))
  545     return fp;
  546 
  547   gl->v = gl->vf->v;
  548   gfs_domain_cell_traverse (GFS_DOMAIN (GFS_GL (gl)->sim),
  549                 FTT_POST_ORDER, FTT_TRAVERSE_ALL, -1,
  550                 (FttCellTraverseFunc) min_max_iso, gl);
  551   gl->minv = G_MAXDOUBLE;
  552   gl->maxv = -G_MAXDOUBLE;
  553   gts_container_foreach (GTS_CONTAINER (GFS_GL (gl)->sim), (GtsFunc) min_max_iso1, gl);
  554 
  555   return NULL;
  556 }
  557 
  558 static void gl_isosurface_set_simulation (GfsGl * object, GfsSimulation * sim)
  559 {
  560   GfsDomain * domain = GFS_DOMAIN (sim);
  561   GfsGlIsosurface * gls = GFS_GL_ISOSURFACE (object);
  562   GtsFile * fp = NULL;
  563   
  564   gfs_gl_isosurface_reset (gls);
  565 
  566   if (GFS_GL_CLASS (GTS_OBJECT_CLASS (gfs_gl_isosurface_class ())->parent_class)->set_simulation)
  567     (*GFS_GL_CLASS (GTS_OBJECT_CLASS (gfs_gl_isosurface_class ())->parent_class)->set_simulation)
  568       (object, sim);
  569 
  570   if (gls->min)
  571     gts_object_destroy (GTS_OBJECT (gls->min));
  572   gls->min = gfs_temporary_variable (domain);
  573   if (gls->max)
  574     gts_object_destroy (GTS_OBJECT (gls->max));
  575   gls->max = gfs_temporary_variable (domain);
  576   if (gls->p)
  577     gts_object_destroy (GTS_OBJECT (gls->p));
  578   gls->p = gfs_temporary_variable (domain);
  579 
  580   gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
  581                 (FttCellTraverseFunc) gfs_cell_reset, gls->p);
  582 
  583   if (gls->expr->str[0] == '\0' || (fp = gfs_gl_isosurface_set (gls, gls->expr->str))) {
  584     GfsDomain * domain = GFS_DOMAIN (sim);
  585     
  586     if (domain->variables) {
  587       GfsVariable * v = domain->variables->data;
  588       gfs_gl_isosurface_set (gls, v->name);
  589     }
  590     else
  591       gfs_gl_isosurface_set (gls, "0");
  592   }
  593   if (fp)
  594     gts_file_destroy (fp);
  595 }
  596 
  597 static void gl_isosurface_class_init (GfsGlClass * klass)
  598 {
  599   GTS_OBJECT_CLASS (klass)->destroy = gl_isosurface_destroy;
  600   GTS_OBJECT_CLASS (klass)->read = gl_isosurface_read;
  601   GTS_OBJECT_CLASS (klass)->write = gl_isosurface_write;
  602   klass->set_simulation = gl_isosurface_set_simulation;
  603   klass->draw = gl_isosurface_draw;
  604   klass->pick = NULL;
  605 }
  606 
  607 static void gl_isosurface_init (GfsGlIsosurface * gl)
  608 {
  609   GtsColor c = { 1., 1., 1. };
  610 
  611   gl->expr = g_string_new ("");
  612   gl->vf = gfs_gl_var_func_new ();
  613   GFS_GL (gl)->lc = c;
  614 }
  615 
  616 GfsGlClass * gfs_gl_isosurface_class (void)
  617 {
  618   static GfsGlClass * klass = NULL;
  619 
  620   if (klass == NULL) {
  621     GtsObjectClassInfo gfs_gl_isosurface_info = {
  622       "GfsGlIsosurface",
  623       sizeof (GfsGlIsosurface),
  624       sizeof (GfsGlScalarClass),
  625       (GtsObjectClassInitFunc) gl_isosurface_class_init,
  626       (GtsObjectInitFunc) gl_isosurface_init,
  627       (GtsArgSetFunc) NULL,
  628       (GtsArgGetFunc) NULL
  629     };
  630     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl_scalar_class ()),
  631                   &gfs_gl_isosurface_info);
  632   }
  633 
  634   return klass;
  635 }
  636