"Fossies" - the Fresh Open Source Software Archive

Member "gfsview-snapshot-121130/gl/gfsgl3D.h" (30 Nov 2012, 32222 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 "gfsgl3D.h" 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 static gboolean plane_cuts_cell (FttCell * cell, gpointer data)
   22 {
   23   GfsGl2D * gl = data;
   24   FttVector p;
   25   ftt_cell_pos (cell, &p);
   26   gdouble r = ftt_cell_size (cell)*SLIGHTLY_LARGER*GFS_DIAGONAL;
   27   return (fabs ((p.x - gl->p[0].x)*gl->n.x + 
   28         (p.y - gl->p[0].y)*gl->n.y + 
   29         (p.z - gl->p[0].z)*gl->n.z)
   30       <= r);
   31 }
   32 
   33 /**
   34  * gfs_gl_cell_traverse_visible_plane:
   35  * @gl: a #GfsGl2D.
   36  * @f: a view frustum.
   37  * @func: a used-defined function.
   38  * @data: user data to pass to @func.
   39  *
   40  * Traverse the cells of @gl which are visible and whose bounding
   41  * sphere is intersected by the plane defined by @gl.
   42  */
   43 static void gfs_gl_cell_traverse_visible_plane (GfsGl * gl,
   44                         GfsFrustum * f,
   45                         FttCellTraverseFunc func,
   46                         gpointer data)
   47 {
   48   g_return_if_fail (gl != NULL);
   49   g_return_if_fail (f != NULL);
   50   g_return_if_fail (func != NULL);
   51 
   52   gfs_gl_cell_traverse_visible_condition (gl, f, plane_cuts_cell, gl, func, data);
   53 }
   54 
   55 /* GfsGlCells: Object */
   56 
   57 static void gl_cell (FttCell * cell, GfsGl * gl)
   58 {
   59   FttVector v[12];
   60   FttDirection d[12];
   61   guint nv = gfs_cut_cube_vertices (cell, gl->maxlevel,
   62                     &GFS_GL2D (gl)->p[0], &GFS_GL2D (gl)->n,
   63                     v, d,
   64                     NULL, NULL);
   65 
   66   if (nv > 2) {
   67     guint i;
   68     glBegin (GL_LINE_LOOP);
   69     for (i = 0; i < nv; i++)
   70       glVertex3d (v[i].x, v[i].y, v[i].z);
   71     glEnd ();
   72     gl->size++;
   73   }
   74 }
   75 
   76 /* GfsGlSquares: Object */
   77 
   78 static void gl_square (FttCell * cell, GfsGl * gl)
   79 {
   80   if (!GFS_HAS_DATA (cell, GFS_GL_SCALAR (gl)->v))
   81     return;
   82   FttVector v[12];
   83   FttDirection d[12];
   84   guint nv = gfs_cut_cube_vertices (cell, gl->maxlevel,
   85                     &GFS_GL2D (gl)->p[0], &GFS_GL2D (gl)->n,
   86                     v, d,
   87                     NULL, NULL);
   88   if (nv > 2) {
   89     GfsGlScalar * gls = GFS_GL_SCALAR (gl);
   90     GtsColor c;
   91     c = gfs_colormap_color (gls->cmap, gls->max > gls->min ?
   92                 (GFS_VALUE (cell, gls->v) - gls->min)/(gls->max - gls->min) :
   93                 0.5);
   94     glColor3f (c.r, c.g, c.b);
   95     guint i;
   96     glBegin (GL_POLYGON);
   97     for (i = 0; i < nv; i++)
   98       glVertex3d (v[i].x, v[i].y, v[i].z);
   99     glEnd ();
  100     gl->size++;
  101   }
  102 }
  103 
  104 static void gl_squares_draw (GfsGl * gl, GfsFrustum * f)
  105 {
  106   gl->size = 0;
  107   gfs_gl_normal (gl);
  108   gfs_gl_cell_traverse_visible_plane (gl, f, (FttCellTraverseFunc) gl_square, gl);
  109 
  110   (* GFS_GL_CLASS (GTS_OBJECT (gl)->klass->parent_class)->draw) (gl, f);
  111 }
  112 
  113 /* GfsGlLinear: Object */
  114 
  115 #define param(v) (gls->max > gls->min ? ((v) - gls->min)/(gls->max - gls->min) : 0.5)
  116 #define color(v) (gfs_colormap_color (gls->cmap, param (v)))
  117 
  118 static void gl_linear_color (FttCell * cell, GfsGl * gl)
  119 {
  120   GfsGlScalar * gls = GFS_GL_SCALAR (gl);
  121   if (!GFS_HAS_DATA (cell, gls->v))
  122     return;
  123   FttVector v[12];
  124   FttDirection d[12];
  125   gdouble val[12];
  126   guint nv = gfs_cut_cube_vertices (cell, gl->maxlevel,
  127                     &GFS_GL2D (gl)->p[0], &GFS_GL2D (gl)->n,
  128                     v, d,
  129                     gls->v, val);
  130 
  131   if (nv > 2) {
  132     guint i;
  133 
  134     glBegin (GL_POLYGON);
  135     for (i = 0; i < nv; i++) {
  136       GtsColor c = color (val[i]);
  137       glColor3f (c.r, c.g, c.b);
  138       glVertex3d (v[i].x, v[i].y, v[i].z);
  139     }
  140     glEnd ();
  141     gl->size++;
  142   }
  143 }
  144 
  145 static void gl_linear_texture (FttCell * cell, GfsGl * gl)
  146 {
  147   FttVector v[12];
  148   FttDirection d[12];
  149   gdouble val[12];
  150   GfsGlScalar * gls = GFS_GL_SCALAR (gl);
  151   guint nv = gfs_cut_cube_vertices (cell, gl->maxlevel,
  152                     &GFS_GL2D (gl)->p[0], &GFS_GL2D (gl)->n,
  153                     v, d,
  154                     gls->v, val);
  155 
  156   if (nv > 2) {
  157     guint i;
  158     FttVector c;
  159     gdouble vc = 0.;
  160 
  161     c.x = c.y = c.z = 0.;
  162     for (i = 0; i < nv; i++) {
  163       c.x += v[i].x; c.y += v[i].y; c.z += v[i].z;
  164       vc += val[i];
  165     }
  166 
  167     glBegin (GL_TRIANGLE_FAN);
  168     glTexCoord1d (param (vc/nv));
  169     glVertex3d (c.x/nv, c.y/nv, c.z/nv);
  170     for (i = 0; i < nv; i++) {
  171       glTexCoord1d (param (val[i]));
  172       glVertex3d (v[i].x, v[i].y, v[i].z);
  173     }
  174     glTexCoord1d (param (val[0]));
  175     glVertex3d (v[0].x, v[0].y, v[0].z);
  176     glEnd ();
  177     gl->size++;
  178   }
  179 }
  180 
  181 static void gl_linear_draw (GfsGl * gl, GfsFrustum * f)
  182 {
  183   gl->size = 0;
  184   glShadeModel (GL_SMOOTH);
  185   gfs_gl_normal (gl);
  186   if (gfs_gl_vector_format (gl))
  187     gfs_gl_cell_traverse_visible_plane (gl, f, (FttCellTraverseFunc) gl_linear_color, gl);
  188   else {
  189     glEnable (GL_TEXTURE_1D);
  190     gfs_colormap_texture (GFS_GL_SCALAR (gl)->cmap);
  191     glColor3f (1., 1., 1.);
  192     gfs_gl_cell_traverse_visible_plane (gl, f, (FttCellTraverseFunc) gl_linear_texture, gl);
  193     glDisable (GL_TEXTURE_1D);
  194   }
  195 
  196   (* GFS_GL_CLASS (GTS_OBJECT (gl)->klass->parent_class)->draw) (gl, f);
  197 }
  198 
  199 /* GfsGlIsoline: Object */
  200 
  201 static void gl_isoline (FttCell * cell, GfsGl * gl)
  202 {
  203   GfsGlScalar * gls = GFS_GL_SCALAR (gl);
  204   if (!GFS_HAS_DATA (cell, gls->v))
  205     return;
  206   FttVector v[12];
  207   FttDirection d[12];
  208   gdouble val[12];
  209   guint nv = gfs_cut_cube_vertices (cell, gl->maxlevel,
  210                     &GFS_GL2D (gl)->p[0], &GFS_GL2D (gl)->n,
  211                     v, d,
  212                     gls->v, val);
  213 
  214   if (nv > 2) {
  215     GArray * levels = GFS_GL_ISOLINE (gl)->levels;
  216     guint i;
  217 
  218     for (i = 0; i < levels->len; i++) {
  219       gdouble z = g_array_index (levels, gdouble, i);
  220       guint j;
  221       guint n = 0;
  222 
  223       val[nv] = val[0];
  224       v[nv] = v[0];
  225       for (j = 0; j < nv; j++)
  226     if ((val[j] > z && val[j + 1] <= z) || (val[j] <= z && val[j + 1] > z)) {
  227       gdouble a = (z - val[j])/(val[j + 1] - val[j]);
  228       glVertex3d (v[j].x + a*(v[j + 1].x - v[j].x),
  229               v[j].y + a*(v[j + 1].y - v[j].y),
  230               v[j].z + a*(v[j + 1].z - v[j].z));
  231       n++;
  232     }
  233       g_assert (n % 2 == 0);
  234     }
  235     gl->size++;
  236   }
  237 }
  238 
  239 static void gl_isoline_draw (GfsGl * gl, GfsFrustum * f)
  240 {
  241   gl_isoline_update_levels (gl);
  242   
  243   gl->size = 0;
  244   glMatrixMode (GL_PROJECTION);
  245   glPushMatrix ();
  246   glTranslatef (0., 0., gl->p->lc);
  247   glBegin (GL_LINES);
  248   gfs_gl_normal (gl);
  249   gfs_gl_cell_traverse_visible_plane (gl, f, (FttCellTraverseFunc) gl_isoline, gl);
  250   glEnd ();
  251   glPopMatrix ();
  252 }
  253 
  254 /* GfsGlVOF: Object */
  255 
  256 static void gl_vof (FttCell * cell, GfsGl * gl)
  257 {
  258   FttVector v[12], m;
  259   guint nv = gfs_vof_facet (cell, GFS_VARIABLE_TRACER_VOF (GFS_GL_VOF (gl)->vf->v), v, &m);
  260 
  261   if (nv > 2) {
  262     guint i;
  263     glBegin (GL_POLYGON);
  264     if (GFS_GL_VOF (gl)->use_scalar) {
  265       GfsGlScalar * gls = GFS_GL_SCALAR (gl);
  266       gdouble value;
  267       GtsColor c;
  268 
  269       if (GFS_GL_VOF (gl)->interpolate) {
  270     FttVector c;
  271     gfs_vof_center (cell, GFS_VARIABLE_TRACER_VOF (GFS_GL_VOF (gl)->vf->v), &c);
  272     value = gfs_interpolate (cell, c, gls->v);
  273       }
  274       else
  275     value = GFS_VALUE (cell, gls->v);
  276       c = gfs_colormap_color (gls->cmap, gls->max > gls->min ?
  277                   (value - gls->min)/(gls->max - gls->min) : 0.5);
  278       glColor3f (c.r, c.g, c.b);
  279     }
  280     if (GFS_GL_VOF (gl)->reversed)
  281       glNormal3d (-m.x, -m.y, -m.z);
  282     else
  283       glNormal3d (m.x, m.y, m.z);
  284     for (i = 0; i < nv; i++)
  285       glVertex3d (v[i].x, v[i].y, v[i].z);
  286     glEnd ();
  287     gl->size++;
  288   }
  289 }
  290 
  291 static void gl_vof_edges (FttCell * cell, GfsGl * gl)
  292 {
  293   FttVector v[12], m;
  294   guint nv = gfs_vof_facet (cell, GFS_VARIABLE_TRACER_VOF (GFS_GL_VOF (gl)->vf->v), v, &m);
  295 
  296   if (nv > 2) {
  297     guint i;
  298     glBegin (GL_LINE_LOOP);
  299     for (i = 0; i < nv; i++)
  300       glVertex3d (v[i].x, v[i].y, v[i].z);
  301     glEnd ();
  302   }
  303 }
  304 
  305 static gboolean is_vof (FttCell * cell, gpointer data)
  306 {
  307   gdouble f = GFS_VALUE (cell, GFS_GL_VOF (data)->vf->v);
  308 
  309   return !GFS_IS_FULL (f);
  310 }
  311 
  312 static void gl_vof_draw (GfsGl * gl, GfsFrustum * f)
  313 {
  314   if (GFS_IS_VARIABLE_TRACER_VOF (GFS_GL_VOF (gl)->vf->v)) {
  315     gl->size = 0;
  316     gfs_gl_cell_traverse_visible_condition (gl, f, is_vof, gl, (FttCellTraverseFunc) gl_vof, gl);
  317     if (GFS_GL_VOF (gl)->draw_edges) {
  318       glMatrixMode (GL_PROJECTION);
  319       glPushMatrix ();
  320       glTranslatef (0., 0., gl->p->lc);
  321       glColor3f (0., 0., 0.);
  322       gfs_gl_cell_traverse_visible_condition (gl, f, is_vof, gl,
  323                           (FttCellTraverseFunc) gl_vof_edges, gl);
  324       glPopMatrix ();
  325     }
  326   }
  327 }
  328 
  329 static gboolean plane_segment_intersection (GfsGl2D * plane, FttVector * A, FttVector * B,
  330                         FttVector * I)
  331 {
  332   GtsVector AB, AP0;
  333   gts_vector_init (AB, A, B);
  334   gts_vector_init (AP0, A, &plane->p[0]);
  335   gdouble lambda = gts_vector_scalar (AB, &plane->n.x);
  336   if (lambda != 0.) {
  337     lambda = gts_vector_scalar (AP0, &plane->n.x)/lambda;
  338     if (lambda >= 0. && lambda < 1.) {
  339       I->x = A->x + lambda*AB[0];
  340       I->y = A->y + lambda*AB[1];
  341       I->z = A->z + lambda*AB[2];
  342       return TRUE;
  343     }
  344   }
  345   return FALSE;
  346 }
  347 
  348 static void gl_vof_cut (GfsGl * gl, FttCell * cell, GfsGl2D * plane)
  349 {
  350   FttVector v[12], m;
  351   guint nv = gfs_vof_facet (cell, GFS_VARIABLE_TRACER_VOF (GFS_GL_VOF (gl)->vf->v), v, &m);
  352 
  353   if (nv > 2) {
  354     FttVector I[3];
  355     guint i, ni = 0;
  356     for (i = 0; i < nv - 1 && ni < 3; i++)
  357       if (plane_segment_intersection (plane, &v[i], &v[i+1], &I[ni]))
  358     ni++;
  359     if (plane_segment_intersection (plane, &v[nv - 1], &v[0], &I[ni]))
  360       ni++;
  361     if (ni == 2) {
  362       glVertex3d (I[0].x, I[0].y, I[0].z);
  363       glVertex3d (I[1].x, I[1].y, I[1].z);
  364     }
  365     GFS_GL (plane)->size++;
  366   }
  367 }
  368 
  369 /* GfsGlSolid: Object */
  370 
  371 typedef struct {
  372   guint m;
  373   FttVector * v;
  374   FttVector * n;
  375   gdouble * val;
  376 } polygon;
  377 
  378 typedef struct {
  379   polygon * p;
  380   guint n;
  381 } polygons;
  382 
  383 #define param(v) (gls->max > gls->min ? ((v) - gls->min)/(gls->max - gls->min) : 0.5)
  384 #define color(v) (gfs_colormap_color (gls->cmap, param (v)))
  385 
  386 static void polygon_draw (polygon * p, GfsGlSolid * gl)
  387 {
  388   guint i;
  389 
  390   glBegin (GL_POLYGON);
  391   for (i = 0; i < p->m; i++) {
  392     if (gl->use_scalar) {
  393       GfsGlScalar * gls = GFS_GL_SCALAR (gl);
  394 
  395       if (gfs_gl_vector_format (GFS_GL (gl))) {
  396     GtsColor c = color (p->val[i]);
  397     glColor3f (c.r, c.g, c.b);
  398       }
  399       else
  400     glTexCoord1d (param (p->val[i]));
  401     }
  402     if (gl->reversed)
  403       glNormal3d (-p->n[i].x, -p->n[i].y, -p->n[i].z);
  404     else
  405       glNormal3d (p->n[i].x, p->n[i].y, p->n[i].z);
  406     glVertex3d (p->v[i].x, p->v[i].y, p->v[i].z);
  407   }
  408   glEnd ();
  409 }
  410 
  411 static void polygons_draw (polygons * p, GfsGlSolid * gl)
  412 {
  413   guint i;
  414 
  415   for (i = 0; i < p->n; i++)
  416     polygon_draw (&p->p[i], gl);
  417 }
  418 
  419 static void polygons_destroy (polygons * p)
  420 {
  421   guint i;
  422   for (i = 0; i < p->n; i++) {
  423     g_free (p->p[i].v);
  424     g_free (p->p[i].n);
  425     g_free (p->p[i].val);
  426   }
  427   g_free (p->p);
  428   g_free (p);
  429 }
  430 
  431 static polygons * polygons_add (polygons * p, guint n)
  432 {
  433   if (p == NULL) {
  434     p = g_malloc (sizeof (polygons));
  435     p->p = g_malloc (sizeof (polygon));
  436     p->n = 0;
  437   }
  438   else
  439     p->p = g_realloc (p->p, sizeof (polygon)*(p->n + 1));
  440   p->p[p->n].v = g_malloc (n*sizeof (FttVector)); 
  441   p->p[p->n].n = g_malloc (n*sizeof (FttVector));
  442   p->p[p->n].val = g_malloc (n*sizeof (gdouble));
  443   p->p[p->n].m = n;
  444   p->n++;
  445   return p;
  446 }
  447 
  448 static void free_polygons (FttCell * cell, GfsGlSolid * gl)
  449 {
  450   gpointer po = GFS_DOUBLE_TO_POINTER (GFS_VALUE (cell, gl->p));
  451   gpointer s = GFS_DOUBLE_TO_POINTER (GFS_VALUE (cell, gl->s));
  452 
  453   if (po) {
  454     polygons_destroy (po);
  455     GFS_VALUE (cell, gl->p) = 0.;
  456   }
  457   if (s) {
  458     gts_object_destroy (s);
  459     GFS_VALUE (cell, gl->s) = 0.;
  460   }
  461 }
  462 
  463 void gfs_gl_solid_reset (GfsGlSolid * gl)
  464 {
  465   g_return_if_fail (gl != NULL);
  466   
  467   if (gl->p && gl->s && GFS_GL (gl)->sim)
  468     gfs_domain_cell_traverse (GFS_DOMAIN (GFS_GL (gl)->sim), 
  469                   FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
  470                   (FttCellTraverseFunc) free_polygons, gl);
  471   gl->needs_updating = TRUE;
  472 }
  473 
  474 static void gl_solid_destroy (GtsObject * o)
  475 {
  476   GfsGlSolid * gl = GFS_GL_SOLID (o);
  477 
  478   gfs_gl_solid_reset (gl);
  479   if (gl->p)
  480     gts_object_destroy (GTS_OBJECT (gl->p));
  481   if (gl->s)
  482     gts_object_destroy (GTS_OBJECT (gl->s));
  483   g_slist_free (gl->solids);
  484 
  485   (* GTS_OBJECT_CLASS (gfs_gl_solid_class ())->parent_class->destroy) (o);
  486 }
  487 
  488 static void gl_solid_read (GtsObject ** o, GtsFile * fp)
  489 {
  490   GfsGlSolid * gl = GFS_GL_SOLID (*o);
  491   GtsFileVariable var[] = {
  492     {GTS_INT, "reversed",   TRUE},
  493     {GTS_INT, "use_scalar", TRUE},
  494     {GTS_NONE}
  495   };
  496 
  497   (* GTS_OBJECT_CLASS (gfs_gl_solid_class ())->parent_class->read) (o, fp);
  498   if (fp->type == GTS_ERROR)
  499     return;
  500 
  501   var[0].data = &gl->reversed;
  502   var[1].data = &gl->use_scalar;
  503   gts_file_assign_variables (fp, var);
  504   if (fp->type == GTS_ERROR)
  505     return;
  506 }
  507 
  508 static void gl_solid_write (GtsObject * o, FILE * fp)
  509 {
  510   GfsGlSolid * gl = GFS_GL_SOLID (o);
  511 
  512   (* GTS_OBJECT_CLASS (gfs_gl_solid_class ())->parent_class->write) (o, fp);
  513 
  514   fprintf (fp, " {\n"
  515        "  reversed = %d\n"
  516        "  use_scalar = %d\n"
  517        "}",
  518        gl->reversed, (gl->use_scalar != NULL));
  519 }
  520 
  521 typedef struct {
  522   GtsPoint p[8];
  523   GfsSegment s[12];
  524   gdouble v[12];
  525   GtsVector n1[12];
  526 } CellCube;
  527 
  528 static void cell_size (FttCell * cell, FttVector * h)
  529 {
  530   h->x = h->y = ftt_cell_size (cell);
  531 #if !FTT_2D
  532   h->z = h->x;
  533 #endif /* 3D */
  534 }
  535 
  536 static gboolean cube_intersections (FttCell * cell,
  537                     GfsGenericSurface * s,
  538                     GfsVariable * v,
  539                     FttVector p[12],
  540                     FttVector n[12],
  541                     gint orient[12],
  542                     gdouble val[12],
  543                     gint max_level)
  544 {
  545   CellCube cube;
  546   FttVector o, h;
  547   guint i;
  548   gint inside[8] = {0,0,0,0,0,0,0,0};
  549   gboolean cut = FALSE;
  550 
  551   ftt_cell_pos (cell, &o);
  552   cell_size (cell, &h);
  553   for (i = 0; i < FTT_DIMENSION; i++)
  554     (&o.x)[i] -= (&h.x)[i]/2.;
  555   for (i = 0; i < 8; i++) { /* for each vertex of the cube */
  556     cube.p[i].x = o.x + h.x*vertex[i].x;
  557     cube.p[i].y = o.y + h.y*vertex[i].y;
  558     cube.p[i].z = o.z + h.z*vertex[i].z;
  559   }
  560 
  561   for (i = 0; i < 12; i++) {
  562     GfsSegment * e = &cube.s[i];
  563     e->E = &cube.p[edge1[i][0]];
  564     e->D = &cube.p[edge1[i][1]];
  565     gfs_surface_segment_intersection (s, cell, e);
  566     if (e->n % 2 != 0) {
  567       gfs_surface_segment_normal (s, cell, e, cube.n1[i]);
  568       cut = TRUE;
  569     }
  570   }
  571   
  572   if (!cut)
  573     return FALSE;
  574 
  575   if (v)
  576     for (i = 0; i < 8; i++) /* for each vertex of the cube */
  577       cube.v[i] = gfs_cell_corner_value (cell, corner[i], v, max_level);
  578   
  579   for (i = 0; i < 12; i++) { /* for each edge of the cube */
  580     GfsSegment * e = &cube.s[i];
  581     if (e->n % 2 != 0) { /* only for odd number of intersections */
  582       guint j = edge1[i][0], k = edge1[i][1];
  583 
  584       /* intersection vertex position is the average of all the n[i] intersections */
  585       e->x /= e->n;
  586 
  587       p[i].x = (1. - e->x)*cube.p[j].x + e->x*cube.p[k].x;
  588       p[i].y = (1. - e->x)*cube.p[j].y + e->x*cube.p[k].y;
  589       p[i].z = (1. - e->x)*cube.p[j].z + e->x*cube.p[k].z;
  590 
  591       n[i].x = cube.n1[i][0];
  592       n[i].y = cube.n1[i][1];
  593       n[i].z = cube.n1[i][2];
  594 
  595       if (v)
  596     val[i] = (1. - e->x)*cube.v[j] + e->x*cube.v[k];
  597 
  598       g_assert (inside[j] == 0 || inside[j] == e->inside);
  599       g_assert (inside[k] == 0 || inside[k] == - e->inside);
  600       inside[j] = e->inside;
  601       inside[k] = - e->inside;
  602       orient[i] = (inside[j] > 0);
  603     }
  604     else
  605       orient[i] = -1;
  606   }
  607 
  608   return TRUE;
  609 }
  610 
  611 #define solid_cube_vertices(cut,s,var,max_level,block) {\
  612   FttVector _p[12], _n[12], * v[12], * n[12];\
  613   gdouble _val[12], val[12];\
  614   gint _orient[12];\
  615   guint _i;\
  616   if (cube_intersections (cell, s, var, _p, _n, _orient, _val, max_level)) \
  617     for (_i = 0; _i < 12; _i++) {                   \
  618       guint nv = 0, _e = _i;                        \
  619       while (_orient[_e] >= 0) {                    \
  620     guint _m = 0, * _ne = connect[_e][_orient[_e]];         \
  621     n[nv] = &(_n[_e]);                      \
  622     val[nv] = _val[_e];                     \
  623     v[nv++] = &(_p[_e]);                        \
  624     _orient[_e] = -1;                       \
  625     while (_m < 3 && _orient[_e] < 0)               \
  626       _e = _ne[_m++];                       \
  627       }                                 \
  628       if (nv > 2) {                         \
  629     block                               \
  630       }                                     \
  631     }\
  632 }
  633 
  634 static void gl_solid (FttCell * cell, GfsGl * gl)
  635 {
  636   GfsGlSolid * gls = GFS_GL_SOLID (gl);
  637   polygons * p = GFS_DOUBLE_TO_POINTER (GFS_VALUE (cell, gls->p));
  638 
  639   if (p) {
  640     polygons_draw (p, gls);
  641     gl->size++;
  642   }
  643 }
  644 
  645 static void gl_solid_update (FttCell * cell, GfsGenericSurface * s, GfsGl * gl)
  646 {
  647   GfsGlSolid * gls = GFS_GL_SOLID (gl);
  648   polygons * p = GFS_DOUBLE_TO_POINTER (GFS_VALUE (cell, gls->p));
  649 
  650   g_assert (!p);
  651   solid_cube_vertices (cut, s, gls->use_scalar ? GFS_GL_SCALAR (gl)->v : NULL, gl->maxlevel, {
  652       guint i;
  653       polygon * q;
  654 
  655       p = polygons_add (p, nv);
  656       q = &p->p[p->n - 1];
  657 
  658       for (i = 0; i < nv; i++) {
  659     q->v[i] = *(v[i]);
  660     q->n[i] = *(n[i]);
  661     q->val[i] = val[i];
  662       }
  663     });
  664   if (p) {
  665     polygons_draw (p, gls);
  666     GFS_DOUBLE_TO_POINTER (GFS_VALUE (cell, gls->p)) = p;
  667     gl->size++;
  668   }
  669 }
  670 
  671 static gboolean is_solid (FttCell * cell, gpointer data)
  672 {
  673   return GFS_IS_MIXED (cell);
  674 }
  675 
  676 static void gl_solid_draw (GfsGl * gl, GfsFrustum * f)
  677 {
  678   GfsGlSolid * gls = GFS_GL_SOLID (gl);
  679 
  680   glShadeModel (GL_SMOOTH);
  681   if (gls->use_scalar && gls->use_scalar != GFS_GL_SCALAR (gl)->v) {
  682     gls->use_scalar = GFS_GL_SCALAR (gl)->v;
  683     gfs_gl_solid_reset (gls);
  684   }
  685   if (gls->use_scalar && !gfs_gl_vector_format (gl)) {
  686     glEnable (GL_TEXTURE_1D);
  687     gfs_colormap_texture (GFS_GL_SCALAR (gl)->cmap);
  688     glColor3f (1., 1., 1.);
  689   }
  690 
  691   gl->size = 0;
  692   if (gls->needs_updating) {
  693     GSList * i = gls->solids;
  694     while (i) {
  695       gfs_domain_traverse_cut (GFS_DOMAIN (gl->sim), GFS_SOLID (i->data)->s, 
  696                    FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS,
  697                    (FttCellTraverseCutFunc) gl_solid_update, gl);
  698       i = i->next;
  699     }
  700     gls->needs_updating = FALSE;
  701   }
  702   else {
  703     gdouble res = f->res;
  704     f->res = 0.;
  705     gfs_gl_cell_traverse_visible_condition (gl, f, is_solid, NULL,
  706                         (FttCellTraverseFunc) gl_solid, gl);
  707     f->res = res;
  708   }
  709 
  710   if (gls->use_scalar && !gfs_gl_vector_format (gl))
  711     glDisable (GL_TEXTURE_1D);
  712 }
  713 
  714 static void reset_p_s (FttCell * cell, GfsGlSolid * gl)
  715 {
  716   GFS_VALUE (cell, gl->p) = 0.;
  717   GFS_VALUE (cell, gl->s) = 0.;
  718 }
  719 
  720 static void gl_solid_set_simulation (GfsGl * object, GfsSimulation * sim)
  721 {
  722   GfsDomain * domain = GFS_DOMAIN (sim);
  723   GfsGlSolid * gls = GFS_GL_SOLID (object);
  724 
  725   gfs_gl_solid_reset (gls);
  726 
  727   (*GFS_GL_CLASS (GTS_OBJECT_CLASS (gfs_gl_solid_class ())->parent_class)->set_simulation)
  728     (object, sim);
  729 
  730   if (gls->p)
  731     gts_object_destroy (GTS_OBJECT (gls->p));
  732   gls->p = gfs_temporary_variable (domain);
  733   if (gls->s)
  734     gts_object_destroy (GTS_OBJECT (gls->s));
  735   gls->s = gfs_temporary_variable (domain);
  736   g_slist_free (gls->solids);
  737   gls->solids = gfs_simulation_get_solids (sim);
  738 
  739   gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
  740                 (FttCellTraverseFunc) reset_p_s, gls);
  741 }
  742 
  743 static gboolean gl_solid_relevant (GfsSimulation * sim)
  744 {
  745   return (sim->solids->items != NULL);
  746 }
  747 
  748 static void gl_solid_class_init (GfsGlClass * klass)
  749 {
  750   GTS_OBJECT_CLASS (klass)->destroy = gl_solid_destroy;
  751   GTS_OBJECT_CLASS (klass)->read = gl_solid_read;
  752   GTS_OBJECT_CLASS (klass)->write = gl_solid_write;
  753   klass->set_simulation = gl_solid_set_simulation;
  754   klass->draw = gl_solid_draw;
  755   klass->pick = NULL;
  756   klass->relevant = gl_solid_relevant;
  757 }
  758 
  759 static void gl_solid_init (GfsGl * gl)
  760 {
  761   GtsColor c = { 1., 1., 1. };
  762 
  763   gl->lc = c;
  764   GFS_GL_SOLID (gl)->needs_updating = TRUE;
  765 }
  766 
  767 GfsGlClass * gfs_gl_solid_class (void)
  768 {
  769   static GfsGlClass * klass = NULL;
  770 
  771   if (klass == NULL) {
  772     GtsObjectClassInfo gfs_gl_solid_info = {
  773       "GfsGlSolid",
  774       sizeof (GfsGlSolid),
  775       sizeof (GfsGlScalarClass),
  776       (GtsObjectClassInitFunc) gl_solid_class_init,
  777       (GtsObjectInitFunc) gl_solid_init,
  778       (GtsArgSetFunc) NULL,
  779       (GtsArgGetFunc) NULL
  780     };
  781     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_gl_scalar_class ()),
  782                   &gfs_gl_solid_info);
  783   }
  784 
  785   return klass;
  786 }
  787 
  788 /* GfsGlFractions: Object */
  789 
  790 static void gl_fractions (FttCell * cell, GfsGl * gl)
  791 {
  792   GfsSolidVector * s = GFS_STATE (cell)->solid;
  793   FttCellFace f;
  794   gdouble h = ftt_cell_size (cell)/2.;
  795 
  796   f.cell = cell;
  797   for (f.d = 0; f.d < FTT_NEIGHBORS; f.d++) 
  798     if (s->s[f.d] > 0. && s->s[f.d] < 1.) {
  799       FttComponent c = f.d/2;
  800       static FttVector o[FTT_DIMENSION][4] = {
  801     {{0.,1.,-1.}, {0.,1.,1.}, {0.,-1.,1.}, {0.,-1.,-1.}},
  802     {{1.,0.,-1.}, {1.,0.,1.}, {-1.,0.,1.}, {-1.,0.,-1.}},
  803     {{1.,-1.,0.}, {1.,1.,0.}, {-1.,1.,0.}, {-1.,-1.,0.}}
  804       };
  805       static FttVector n[FTT_NEIGHBORS] = {
  806     { 1., 0., 0. }, { -1., 0., 0.},
  807     { 0., 1., 0. }, { 0., -1., 0.},
  808     { 0., 0., 1. }, { 0., 0., -1.},
  809       };
  810       FttVector p;
  811       gdouble l = h*sqrt (s->s[f.d]);
  812 
  813       glNormal3d (n[f.d].x, n[f.d].y, n[f.d].z);
  814       ftt_face_pos (&f, &p);
  815       glVertex3d (p.x + l*o[c][0].x, p.y + l*o[c][0].y, p.z + l*o[c][0].z);
  816       glVertex3d (p.x + l*o[c][1].x, p.y + l*o[c][1].y, p.z + l*o[c][1].z);
  817       glVertex3d (p.x + l*o[c][2].x, p.y + l*o[c][2].y, p.z + l*o[c][2].z);
  818       glVertex3d (p.x + l*o[c][3].x, p.y + l*o[c][3].y, p.z + l*o[c][3].z);
  819     }
  820   gl->size++;
  821 }
  822 
  823 static void gl_fractions_draw (GfsGl * gl, GfsFrustum * f)
  824 {
  825   gl->size = 0;
  826   glShadeModel (GL_CONSTANT);
  827   glBegin (GL_QUADS);
  828   gfs_gl_cell_traverse_visible_mixed (gl, f, (FttCellTraverseFunc) gl_fractions, gl);
  829   glEnd ();
  830 }
  831 
  832 /* GfsGlBoundaries: Object */
  833 
  834 static void gl_boundaries (FttCell * cell, GfsGl * gl)
  835 {
  836   if (!GFS_IS_MIXED (cell)) {
  837     FttCellNeighbors n;
  838     gdouble h = ftt_cell_size (cell);
  839     FttVector p;
  840     guint i;
  841 
  842     ftt_cell_neighbors (cell, &n);
  843     ftt_cell_pos (cell, &p);
  844     p.x -= h/2.; p.y -= h/2.; p.z -= h/2.;
  845     for (i = 0; i < 12; i++) {
  846       static FttDirection edge2[12][2] = {
  847     {FTT_BOTTOM,FTT_BACK},{FTT_BOTTOM,FTT_FRONT},
  848     {FTT_TOP,FTT_FRONT},  {FTT_TOP,FTT_BACK},
  849     {FTT_LEFT,FTT_BACK},  {FTT_LEFT,FTT_FRONT},
  850     {FTT_RIGHT,FTT_FRONT},{FTT_RIGHT,FTT_BACK},
  851     {FTT_BOTTOM,FTT_LEFT},{FTT_BOTTOM,FTT_RIGHT},
  852     {FTT_TOP,FTT_RIGHT},   {FTT_TOP,FTT_LEFT}
  853       };
  854       guint j = edge2[i][0], k = edge2[i][1];
  855 
  856       if ((!n.c[j] || GFS_CELL_IS_BOUNDARY (n.c[j])) &&
  857       (!n.c[k] || GFS_CELL_IS_BOUNDARY (n.c[k]))) {
  858     gl->size++;
  859     glVertex3d (p.x + edge[i][0].x*h, 
  860             p.y + edge[i][0].y*h,
  861             p.z + edge[i][0].z*h);
  862     glVertex3d (p.x + edge[i][1].x*h, 
  863             p.y + edge[i][1].y*h,
  864             p.z + edge[i][1].z*h);
  865       }
  866     }
  867   }
  868 }
  869 
  870 /* GfsGlLevels: Object */
  871 
  872 static void gl_face (FttCell * cell, GfsGlLevels * gl)
  873 {
  874   FttVector v[12];
  875   FttDirection d[12];
  876   guint nv = gfs_cut_cube_vertices (cell, GFS_GL (gl)->maxlevel,
  877                     &GFS_GL2D (gl)->p[0], &GFS_GL2D (gl)->n,
  878                     v, d,
  879                     NULL, NULL);
  880   gboolean drawn = FALSE;
  881 
  882   if (nv > 2) {
  883     guint i;
  884     for (i = 0; i < nv - 1; i++) {
  885       FttCell * n = ftt_cell_neighbor (cell, d[i]);
  886       if (n && (floor (GFS_VALUE (cell, gl->v)) != 
  887         floor (GFS_VALUE (n, gl->v)))) {
  888     glVertex3d (v[i].x, v[i].y, v[i].z);
  889     glVertex3d (v[i+1].x, v[i+1].y, v[i+1].z);
  890     drawn = TRUE;
  891       }
  892     }
  893   }
  894   if (drawn)
  895     GFS_GL (gl)->size++;
  896 }
  897 
  898 /* GfsGlVectors: Object */
  899 
  900 static void gl_vector (FttCell * cell, GfsGl * gl)
  901 {
  902   GfsGlVectors * gls = GFS_GL_VECTORS (gl);
  903   if (gls->use_scalar && !GFS_HAS_DATA (cell, GFS_GL_SCALAR (gl)->v))
  904     return;
  905   if (gfs_plane_cuts_cell (GFS_GL2D (gl)->p, cell)) {
  906     GfsGl2D * gl2D = GFS_GL2D (gl);
  907     FttComponent c;
  908     FttVector pos, f;
  909     gdouble a;
  910     
  911     gl->size++;
  912     if (gls->use_scalar) {
  913       GfsGlScalar * gla = GFS_GL_SCALAR (gl);
  914       GtsColor c = gfs_colormap_color (gla->cmap, gla->max > gla->min ?
  915                        (GFS_VALUE (cell, gla->v) - gla->min)/
  916                        (gla->max - gla->min) :
  917                        0.5);
  918       glColor3f (c.r, c.g, c.b);
  919     }
  920     gfs_cell_cm (cell, &pos);
  921     a = ((gl2D->p[0].x - pos.x)*gl2D->n.x + 
  922      (gl2D->p[0].y - pos.y)*gl2D->n.y + 
  923      (gl2D->p[0].z - pos.z)*gl2D->n.z);
  924     pos.x += a*gl2D->n.x; pos.y += a*gl2D->n.y; pos.z += a*gl2D->n.z;
  925     for (c = 0; c < FTT_DIMENSION; c++)
  926       (&f.x)[c] = GFS_VALUE (cell, gls->v[c]);
  927     f.x *= gls->scale;
  928     f.y *= gls->scale;
  929     f.z *= gls->scale;
  930     glVertex3d (pos.x + f.x - (f.x - f.y/2.)/5., pos.y + f.y - (f.x/2. + f.y)/5., pos.z + f.z);
  931     glVertex3d (pos.x + f.x, pos.y + f.y, pos.z + f.z);
  932     glVertex3d (pos.x + f.x, pos.y + f.y, pos.z + f.z);
  933     glVertex3d (pos.x + f.x - (f.x + f.y/2.)/5., pos.y + f.y + (f.x/2. - f.y)/5., pos.z + f.z);
  934     glVertex3d (pos.x, pos.y, pos.z);
  935     glVertex3d (pos.x + f.x, pos.y + f.y, pos.z + f.z);
  936   }
  937 }
  938 
  939 /* GfsGlEllipses: Object */
  940 
  941 static void gl_ellipse (FttCell * cell, GfsGl * gl)
  942 {
  943   GfsGlEllipses * gls = GFS_GL_ELLIPSES (gl);
  944   if (gls->use_scalar && !GFS_HAS_DATA (cell, GFS_GL_SCALAR (gl)->v))
  945     return;
  946   if (gfs_plane_cuts_cell (GFS_GL2D (gl)->p, cell)) {
  947     GfsGl2D * gl2D = GFS_GL2D (gl);
  948     FttVector pos;
  949     gdouble a, t;
  950     
  951     gl->size++;
  952     if (gls->use_scalar) {
  953       GfsGlScalar * gla = GFS_GL_SCALAR (gl);
  954       GtsColor c = gfs_colormap_color (gla->cmap, gla->max > gla->min ?
  955                        (GFS_VALUE (cell, gla->v) - gla->min)/
  956                        (gla->max - gla->min) :
  957                        0.5);
  958       glColor3f (c.r, c.g, c.b);
  959     }
  960     gfs_cell_cm (cell, &pos);
  961     a = ((gl2D->p[0].x - pos.x)*gl2D->n.x + 
  962      (gl2D->p[0].y - pos.y)*gl2D->n.y + 
  963      (gl2D->p[0].z - pos.z)*gl2D->n.z);
  964     pos.x += a*gl2D->n.x; pos.y += a*gl2D->n.y; pos.z += a*gl2D->n.z;
  965     glBegin (GL_LINE_LOOP);
  966     for (t = 0.; t < 2.*M_PI; t += 2.*M_PI/20.) {
  967       gdouble cost = cos (t), sint = sin (t);
  968       
  969       glVertex3d (pos.x + gls->scale*(GFS_VALUE (cell, gls->v[0])*cost + 
  970                       GFS_VALUE (cell, gls->v[2])*sint),
  971           pos.y + gls->scale*(GFS_VALUE (cell, gls->v[1])*cost + 
  972                       GFS_VALUE (cell, gls->v[3])*sint),
  973           pos.z);
  974     }
  975     glEnd ();
  976   }
  977 }
  978 
  979 /* GfsGlLocate: Object */
  980 
  981 static void gl_locate (FttCell * cell, GfsGl * gl)
  982 {
  983   FttVector p;
  984   gdouble size = ftt_cell_size (cell)/2.;
  985   
  986   ftt_cell_pos (cell, &p);
  987   glBegin (GL_LINE_LOOP);
  988   glVertex3d (p.x - size, p.y - size, p.z - size);
  989   glVertex3d (p.x + size, p.y - size, p.z - size);
  990   glVertex3d (p.x + size, p.y + size, p.z - size);
  991   glVertex3d (p.x - size, p.y + size, p.z - size);
  992   glEnd ();
  993   glBegin (GL_LINE_LOOP);
  994   glVertex3d (p.x - size, p.y - size, p.z + size);
  995   glVertex3d (p.x + size, p.y - size, p.z + size);
  996   glVertex3d (p.x + size, p.y + size, p.z + size);
  997   glVertex3d (p.x - size, p.y + size, p.z + size);
  998   glEnd ();
  999   glBegin (GL_LINES);
 1000   glVertex3d (p.x - size, p.y - size, p.z - size);
 1001   glVertex3d (p.x - size, p.y - size, p.z + size);
 1002   glVertex3d (p.x + size, p.y - size, p.z - size);
 1003   glVertex3d (p.x + size, p.y - size, p.z + size);
 1004   glVertex3d (p.x + size, p.y + size, p.z - size);
 1005   glVertex3d (p.x + size, p.y + size, p.z + size);
 1006   glVertex3d (p.x - size, p.y + size, p.z - size);
 1007   glVertex3d (p.x - size, p.y + size, p.z + size);
 1008   glEnd ();  
 1009   gl->size++;
 1010 }
 1011 
 1012 /* GfsGlStreamline: Object */
 1013 
 1014 static GSList * circle_profile (GtsVertexClass * klass, 
 1015                 gdouble radius, guint np)
 1016 {
 1017   GSList * lp = NULL;
 1018   guint i;
 1019 
 1020   for (i = 0; i <= np; i++) {
 1021     gdouble a = 2.*M_PI*i/(gdouble) np;
 1022     gdouble cosa = cos (a), sina = sin (a);
 1023     GtsPoint * p = gts_point_new (GTS_POINT_CLASS (klass), radius*cosa, radius*sina, 0.);
 1024     gdouble * n = GTS_VERTEX_NORMAL (p)->n;
 1025 
 1026     n[0] = cosa; n[1] = sina; n[2] = 0.;
 1027     lp = g_slist_prepend (lp, p);
 1028   }
 1029   return lp;
 1030 }
 1031 
 1032 static void matrix_transpose (GtsMatrix * m)
 1033 {
 1034   guint i, j;
 1035 
 1036   for (i = 1; i < 3; i++)
 1037     for (j = 0; j < i; j++) {
 1038       gdouble t = m[i][j];
 1039 
 1040       m[i][j] = m[j][i];
 1041       m[j][i] = t;
 1042     }
 1043 }
 1044 
 1045 static void base (GtsMatrix * b, GtsPoint * p1, GtsPoint * p2)
 1046 {
 1047   GtsVector x, y;
 1048 
 1049   x[0] = b[0][0];
 1050   x[1] = b[1][0];
 1051   x[2] = b[2][0];
 1052   gts_vector_init (b[2], p2, p1);
 1053   gts_vector_normalize (b[2]);
 1054   gts_vector_cross (y, b[2], x);
 1055   if (gts_vector_norm (y) > 1e-2) {
 1056     b[1][0] = y[0];
 1057     b[1][1] = y[1];
 1058     b[1][2] = y[2];
 1059     gts_vector_normalize (b[1]);
 1060   }
 1061   gts_vector_cross (b[0], b[1], b[2]);
 1062   gts_vector_normalize (b[0]);
 1063   matrix_transpose (b);
 1064 }
 1065 
 1066 static void point_list (GtsMatrix * b, 
 1067             GtsMatrix * c,
 1068             GtsPoint * o,
 1069             GSList * profile,
 1070             GtsVertexNormal * pl, 
 1071             guint np)
 1072 {
 1073   guint i;
 1074 
 1075 #if 0
 1076   gboolean colored = FALSE;
 1077   if (GTS_IS_COLORED_VERTEX (o) && 
 1078       gts_object_class_is_from_class (GTS_OBJECT_CLASS (s->vertex_class),
 1079                       GTS_OBJECT_CLASS (gts_colored_vertex_class ())))
 1080     colored = TRUE;
 1081 #endif
 1082   if (FALSE/*GFS_IS_TWISTED_VERTEX (o)*/) {
 1083 #if 0
 1084     gdouble t = GFS_TWISTED_VERTEX (o)->theta;
 1085     gdouble sint = sin (t), cost = cos (t);
 1086     GtsMatrix * r = gts_matrix_new (cost, -sint, 0., 0.,
 1087                     sint,  cost, 0., 0.,
 1088                     0.,      0., 1., 0.,
 1089                     0.,      0., 0., 0.);
 1090     
 1091     c = gts_matrix_product (b, r);
 1092     gts_matrix_destroy (r);
 1093 #endif
 1094   }
 1095   else
 1096     gts_matrix_assign (c,
 1097                b[0][0], b[0][1], b[0][2], 0.,
 1098                b[1][0], b[1][1], b[1][2], 0.,
 1099                b[2][0], b[2][1], b[2][2], 0.,
 1100                0., 0., 0., 0.);
 1101 
 1102   for (i = 0; i < np; i++, profile = profile->next, pl++) {
 1103     GtsPoint * p = profile->data;
 1104     gdouble * n = GTS_VERTEX_NORMAL (p)->n;
 1105     GtsPoint n1;
 1106     
 1107     *GTS_POINT (pl) = *p;
 1108 #if 0
 1109     if (colored)
 1110       GTS_COLORED_VERTEX (v)->c = GTS_COLORED_VERTEX (o)->c;
 1111 #endif
 1112     gts_point_transform (GTS_POINT (pl), c);
 1113     GTS_POINT (pl)->x += o->x;
 1114     GTS_POINT (pl)->y += o->y;
 1115     GTS_POINT (pl)->z += o->z;
 1116 
 1117     n1.x = n[0]; n1.y = n[1]; n1.z = n[2];
 1118     gts_point_transform (&n1, c);
 1119     n = pl->n;
 1120     n[0] = n1.x; n[1] = n1.y; n[2] = n1.z;
 1121   }
 1122 }
 1123 
 1124 static void draw_point (GtsPoint * p)
 1125 {
 1126   glVertex3d (p->x, p->y, p->z);
 1127 }
 1128 
 1129 static void draw_normal (GtsVertexNormal * v)
 1130 {
 1131   glNormal3d (v->n[0], v->n[1], v->n[2]);
 1132 }
 1133 
 1134 static void draw_faces (GtsVertexNormal * v1, GtsVertexNormal * v2, guint np)
 1135 {
 1136   guint i;
 1137 
 1138   glBegin (GL_QUAD_STRIP);
 1139   for (i = 0; i < np; i++, v1++, v2++) {
 1140     draw_normal (v1);
 1141     draw_point (GTS_POINT (v1));
 1142     draw_normal (v2);
 1143     draw_point (GTS_POINT (v2));    
 1144   }
 1145   glEnd ();
 1146 }
 1147 
 1148 static GList * next_far_enough (GList * p, gdouble size)
 1149 {
 1150   GtsPoint * ps;
 1151   GList * pf = NULL;
 1152 
 1153   if (p == NULL)
 1154     return NULL;
 1155   ps = p->data;
 1156   p = p->next;
 1157   size *= size;
 1158   while (p && !pf) {
 1159     if (gts_point_distance2 (ps, p->data) > size)
 1160       pf = p;
 1161     p = p->next;
 1162   }
 1163   return pf;
 1164 }
 1165 
 1166 static void extrude_profile (GSList * profile, GList * path)
 1167 {
 1168   GtsMatrix * r, * c;
 1169   GtsPoint * p0, * p1, * p2;
 1170   GtsVertexNormal * pl1, * pl2, * tmp;
 1171   GtsBBox * bb;
 1172   gdouble size;
 1173   guint np;
 1174 
 1175   g_return_if_fail (profile != NULL);
 1176   g_return_if_fail (path != NULL);
 1177 
 1178   bb = gts_bbox_points (gts_bbox_class (), profile);
 1179   size = bb->x2 - bb->x1;
 1180   if (bb->y2 - bb->y1 > size)
 1181     size = bb->y2 - bb->y1;
 1182   gts_object_destroy (GTS_OBJECT (bb));
 1183 
 1184   size /= 4.;
 1185 
 1186   p0 = path->data;
 1187   path = next_far_enough (path, size);
 1188   if (path == NULL)
 1189     return;
 1190   p1 = path->data;
 1191   r = gts_matrix_identity (NULL);
 1192   c = gts_matrix_identity (NULL);
 1193   np = g_slist_length (profile);
 1194   pl1 = g_malloc (sizeof (GtsVertexNormal)*np);
 1195   pl2 = g_malloc (sizeof (GtsVertexNormal)*np);
 1196 
 1197   base (r, p0, p1);
 1198   point_list (r, c, p0, profile, pl1, np);
 1199   do {
 1200     path = next_far_enough (path, size);
 1201     p2 = path ? path->data : NULL;
 1202     if (p2)
 1203       base (r, p0, p2);
 1204     else
 1205       base (r, p0, p1);
 1206     point_list (r, c, p1, profile, pl2, np);
 1207     draw_faces (pl1, pl2, np);
 1208     tmp = pl1;
 1209     pl1 = pl2;
 1210     pl2 = tmp;
 1211     p0 = p1;
 1212     p1 = p2;
 1213   } while (p1);
 1214 
 1215   g_free (pl1);
 1216   g_free (pl2);
 1217   gts_matrix_destroy (c);
 1218   gts_matrix_destroy (r);
 1219 }
 1220 
 1221 static void gl_streamline_draw (GfsGlStreamline * s,
 1222                 GfsGlStreamlines * gls)
 1223 {
 1224   if (gls->radius > 0.) {
 1225     GSList * profile = circle_profile (gts_vertex_normal_class (), gls->radius, 10);
 1226 
 1227     glShadeModel (GL_SMOOTH);
 1228     extrude_profile (profile, s->l);
 1229 
 1230     g_slist_foreach (profile, (GFunc) gts_object_destroy, NULL);
 1231     g_slist_free (profile);
 1232   }
 1233   else {
 1234     glBegin (GL_LINE_STRIP);
 1235     g_list_foreach (s->l, (GFunc) draw_point, NULL);
 1236     glEnd ();
 1237   }
 1238 }
 1239