pot2surf.c - OpenGrok cross reference for /dports/science/v_sim/v_sim-3.8.0/src/extraFunctions/pot2surf.c

/*   EXTRAITS DE LA LICENCE
	Copyright CEA, contributeurs : Luc BILLARD, Damien
	CALISTE, Olivier D'Astier, laboratoire L_Sim, (2001-2005)

	Adresse mèl :
	BILLARD, non joignable par mèl ;
	CALISTE, damien P caliste AT cea P fr.
	D'ASTIER, dastier AT iie P cnam P fr.

	Ce logiciel est un programme informatique servant à visualiser des
	structures atomiques dans un rendu pseudo-3D.

	Ce logiciel est régi par la licence CeCILL soumise au droit français et
	respectant les principes de diffusion des logiciels libres. Vous pouvez
	utiliser, modifier et/ou redistribuer ce programme sous les conditions
	de la licence CeCILL telle que diffusée par le CEA, le CNRS et l'INRIA
	sur le site "http://www.cecill.info".

	Le fait que vous puissiez accéder à cet en-tête signifie que vous avez
	pris connaissance de la licence CeCILL, et que vous en avez accepté les
	termes (cf. le fichier Documentation/licence.fr.txt fourni avec ce logiciel).
*/

/*   LICENCE SUM UP
	Copyright CEA, contributors : Luc BILLARD and Damien
	CALISTE and Olivier D'Astier, laboratoire L_Sim, (2001-2005)

	E-mail addresses :
	BILLARD, not reachable any more ;
	CALISTE, damien P caliste AT cea P fr.
	D'ASTIER, dastier AT iie P cnam P fr.

	This software is a computer program whose purpose is to visualize atomic
	configurations in 3D.

	This software is governed by the CeCILL  license under French law and
	abiding by the rules of distribution of free software.  You can  use,
	modify and/ or redistribute the software under the terms of the CeCILL
	license as circulated by CEA, CNRS and INRIA at the following URL
	"http://www.cecill.info".

	The fact that you are presently reading this means that you have had
	knowledge of the CeCILL license and that you accept its terms. You can
	find a copy of this licence shipped with this software at Documentation/licence.en.txt.
*/
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <unistd.h>

#include <math.h>

#include <visu_tools.h>
#include "surfaces.h"
#include "pot2surf.h"

/**
 * SECTION:pot2surf
 * @short_description: Creates surfaces from scalar fields.
 *
 * <para>Originally written by Luc Billard for his program
 * VISUALISE. It has been since transformed to a single function and
 * integrated in V_Sim. .pot file are text files which specification
 * are the following :
 * <itemizedlist>
 * <listitem><para>1st line: full pathname of the potential file to
 * read</para></listitem>
 * <listitem><para>2nd line: full pathname of the surface file to
 * build</para></listitem>
 * <listitem><para>3rd line: an integer giving the nbr n of
 * isosurfaces to build </para></listitem>
 * <listitem><para>Each of the n following lines must match the
 * pattern [value name] where value is a real number for the isovalue
 * and name is the name given for the corresponding isosurface to
 * build. Each surface should be named surface_*.</para></listitem>
 * </itemizedlist>
 * The function will fail if it finds no isosurface corresponding to
 * some of the given isovalues. The <link
 * linkend="v-sim-panelSurfacesTools">panelSurfacesTools</link>
 * contains a frontend to build valid .instruc files. This panel is
 * originally integrated in V_Sim. You can access it through the
 * Convert tab in the Isosurfaces panel.</para>
 */

#define ISOSURFACES_FLAG_POTENTIAL "# potentialValue"


static int nx, ny, nz;
static double dxx, dyx, dyy, dzx, dzy, dzz;
static double exx, eyx, eyy, ezx, ezy, ezz;
static int nxm1, nym1, nzm1;
static int nz3, nyz3, nxyz3;
static double dxx1, dyx1, dyy1, dzx1, dzy1, dzz1;
static double ***f = NULL;
static double isovalue;
static FILE *in = NULL;

static int *itab = NULL;
static double *xs = NULL, *ys = NULL, *zs = NULL, *xns = NULL, *yns = NULL, *zns = NULL;
static gboolean create_fromSF_uniform(VisuSurface **surf, const VisuScalarField *field,
                                      double isoValue, const gchar *name);
static gboolean create_fromSF_nonuniform(VisuSurface **surf, const VisuScalarField *field,
                                         double isoValue, const gchar *name);
static gboolean Create_surf(int nelez3, int neleyz3, int nelexyz3,
                            guint sizem1[3], int *iTab, GArray *points,
                            double isoValue, const VisuScalarField *field, double scalarBox[6],
                            const gchar *name, VisuSurface **surf);

/* Tables from :
   http://www.swin.edu.au/astronomy/pbourke/modelling/polygonise/index.html
   (valid link on February 25, 2002)

   YOU WILL FIND MANY VERY INTERESTING PAGES at Paul Bourke's site:
   http://astronomy.swin.edu.au/pbourke/

   Other routines by Luc Billard, on March 1, 2002
*/

static int edgeTable[256]={
0x0  , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c,
0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac,
0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c,
0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc,
0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c,
0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc ,
0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460,
0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0,
0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230,
0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0   };

static int triTable[256][16] =
{{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1},
{3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1},
{3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
{3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1},
{9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1},
{9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
{2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1},
{8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1},
{9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
{4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1},
{3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1},
{1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1},
{4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1},
{4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
{5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1},
{2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1},
{9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
{0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
{2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1},
{10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1},
{5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1},
{5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1},
{9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1},
{0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1},
{10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1},
{8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1},
{2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1},
{7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1},
{2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1},
{11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1},
{5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1},
{11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1},
{11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
{1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1},
{9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1},
{5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1},
{2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
{5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1},
{6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1},
{3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1},
{6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1},
{5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1},
{1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
{10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1},
{6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1},
{8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1},
{7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1},
{3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
{5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1},
{0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1},
{9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1},
{8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1},
{5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1},
{0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1},
{6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1},
{10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1},
{10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1},
{8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1},
{1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1},
{0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1},
{10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1},
{3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1},
{6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1},
{9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1},
{8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1},
{3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1},
{6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1},
{0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1},
{10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1},
{10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1},
{2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1},
{7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1},
{7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1},
{2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1},
{1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1},
{11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1},
{8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1},
{0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1},
{7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
{10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
{2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
{6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1},
{7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1},
{2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1},
{1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1},
{10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1},
{10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1},
{0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1},
{7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1},
{6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1},
{8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1},
{9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1},
{6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1},
{4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1},
{10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1},
{8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1},
{0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1},
{1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1},
{8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1},
{10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1},
{4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1},
{10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
{5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
{11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1},
{9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
{6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1},
{7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1},
{3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1},
{7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1},
{3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1},
{6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1},
{9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1},
{1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1},
{4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1},
{7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1},
{6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1},
{3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1},
{0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1},
{6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1},
{0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1},
{11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1},
{6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1},
{5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1},
{9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1},
{1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1},
{10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1},
{0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1},
{5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1},
{10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1},
{11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1},
{9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1},
{7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1},
{2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1},
{8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1},
{9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1},
{9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1},
{1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1},
{9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1},
{9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1},
{5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1},
{0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1},
{10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1},
{2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1},
{0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1},
{0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1},
{9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1},
{5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1},
{3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1},
{5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1},
{8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1},
{0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1},
{9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1},
{1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1},
{3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1},
{4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1},
{9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1},
{11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1},
{11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1},
{2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1},
{9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1},
{3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1},
{1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1},
{4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
{3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1},
{0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1},
{9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1},
{1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}};

/******************************************************************************/

static double gx(int i, int j, int k) {
   if(i > 0 && i < nxm1) {
      return (f[i+1][j][k]-f[i-1][j][k]) / 2.0;
   }
   else if(i == 0) {
      return (f[i+1][j][k]-f[i][j][k]);
   }
   else {
      return (f[i][j][k]-f[i-1][j][k]);
   }
}

/******************************************************************************/

static double gy(int i, int j, int k) {
   if(j > 0 && j < nym1) {
      return (f[i][j+1][k]-f[i][j-1][k]) / 2.0;
   }
   else if(j == 0) {
      return (f[i][j+1][k]-f[i][j][k]);
   }
   else {
      return (f[i][j][k]-f[i][j-1][k]);
   }
}

/******************************************************************************/

static double gz(int i, int j, int k) {
   if(k > 0 && k < nzm1) {
      return (f[i][j][k+1]-f[i][j][k-1]) / 2.0;
   }
   else if(k == 0) {
      return (f[i][j][k+1]-f[i][j][k]);
   }
   else {
      return (f[i][j][k]-f[i][j][k-1]);
   }
}

/******************************************************************************/

static GString* calc(guint *n1, guint *n2) {

/*  I label local cell vertices as in ref. above
    and assigns x y z names like this :

                 y
                 |

                 4----------5
                /.         /|
               7----------6 |
               | .        | |
               | .        | |
               | 0........|.1  --x
               |/         |/
               3----------2

             /
            z

*/
   register int i, j, k;
   guint n;
   double fac;
   double xu, yu, zu, xui, xuij, yuj;
   int cubeindex;
   int e;
   int vertlist[12];
   int m1, m2, m3;
   double ux, uy, uz, vx, vy, vz, wx, wy, wz, s;
   double xnu, ynu, znu;
   GString *out;

   /* edges are labelled from each corner:
      x+ then z+ then y+
      for k varying fastest
      then j varying fastest
      then i varying fastest
   */
   n = 0;
   for(i=0; i<nx; i++) {
      xui = i*dxx1;
      for(j=0; j<ny; j++) {
         xuij = xui+j*dyx1;
         yuj = j*dyy1;
         for(k=0; k<nz; k++) {
            xu = xuij+k*dzx1;
            yu = yuj+k*dzy1;
            zu = k*dzz1;
            xnu = gx(i, j, k);
            ynu = gy(i, j, k);
            znu = gz(i, j, k);
            if(i<nxm1) {
               if((isovalue-f[i][j][k] < 0.0 && isovalue-f[i+1][j][k] >= 0) ||
                  (isovalue-f[i][j][k] >= 0.0 && isovalue-f[i+1][j][k] < 0)) {
                  fac = (isovalue-f[i][j][k])/(f[i+1][j][k]-f[i][j][k]);
                  itab[n] = *n2;
                  xs[*n2] = xu + fac*dxx1;
                  ys[*n2] = yu;
                  zs[*n2] = zu;
                  ux = xnu + fac*(gx(i+1, j, k) - xnu);
                  uy = ynu + fac*(gy(i+1, j, k) - ynu);
                  uz = znu + fac*(gz(i+1, j, k) - znu);
                  ux = ux*exx+uy*eyx+uz*ezx;
                  uy =        uy*eyy+uz*ezy;
                  uz =               uz*ezz;
                  s = sqrt(ux*ux+uy*uy+uz*uz);
                  if(s <= 0.0) {
                    g_warning("normal direction assumed");
                     xns[*n2] = 1.0;
                     yns[*n2] = 0.0;
                     zns[*n2] = 0.0;
                  }
                  else {
                     xns[*n2] = ux/s;
                     yns[*n2] = uy/s;
                     zns[*n2] = uz/s;
                  }
                  *n2 += 1;
               }
               else {
                  itab[n] = -1;
               }
            }
            else {
               itab[n] = -1;
            }
            n++;
            if(k<nzm1) {
               if((isovalue-f[i][j][k] < 0.0 && isovalue-f[i][j][k+1] >= 0) ||
                  (isovalue-f[i][j][k] >= 0.0 && isovalue-f[i][j][k+1] < 0)) {
                  fac = (isovalue-f[i][j][k])/(f[i][j][k+1]-f[i][j][k]);
                  itab[n] = *n2;
                  xs[*n2] = xu + fac*dzx1;
                  ys[*n2] = yu + fac*dzy1;
                  zs[*n2] = zu + fac*dzz1;
                  ux = xnu + fac*(gx(i, j, k+1) - xnu);
                  uy = ynu + fac*(gy(i, j, k+1) - ynu);
                  uz = znu + fac*(gz(i, j, k+1) - znu);
                  ux = ux*exx+uy*eyx+uz*ezx;
                  uy =        uy*eyy+uz*ezy;
                  uz =               uz*ezz;
                  s = sqrt(ux*ux+uy*uy+uz*uz);
                  if(s <= 0.0) {
                    g_warning("normal direction assumed");
                     xns[*n2] = 1.0;
                     yns[*n2] = 0.0;
                     zns[*n2] = 0.0;
                  }
                  else {
                     xns[*n2] = ux/s;
                     yns[*n2] = uy/s;
                     zns[*n2] = uz/s;
                  }
                  *n2 += 1;
               }
               else {
                  itab[n] = -1;
               }
            }
            else {
               itab[n] = -1;
            }
            n++;
            if(j<nym1) {
               if((isovalue-f[i][j][k] < 0.0 && isovalue-f[i][j+1][k] >= 0) ||
                  (isovalue-f[i][j][k] >= 0.0 && isovalue-f[i][j+1][k] < 0)) {
                  fac = (isovalue-f[i][j][k])/(f[i][j+1][k]-f[i][j][k]);
                  itab[n] = *n2;
                  xs[*n2] = xu + fac*dyx1;
                  ys[*n2] = yu + fac*dyy1;
                  zs[*n2] = zu;
                  ux = xnu + fac*(gx(i, j+1, k) - xnu);
                  uy = ynu + fac*(gy(i, j+1, k) - ynu);
                  uz = znu + fac*(gz(i, j+1, k) - znu);
                  ux = ux*exx+uy*eyx+uz*ezx;
                  uy =        uy*eyy+uz*ezy;
                  uz =               uz*ezz;
                  s = sqrt(ux*ux+uy*uy+uz*uz);
                  if(s <= 0.0) {
                    g_warning("normal direction assumed");
                     xns[*n2] = 1.0;
                     yns[*n2] = 0.0;
                     zns[*n2] = 0.0;
                  }
                  else {
                     xns[*n2] = ux/s;
                     yns[*n2] = uy/s;
                     zns[*n2] = uz/s;
                  }
                  *n2 += 1;
               }
               else {
                  itab[n] = -1;
               }
            }
            else {
               itab[n] = -1;
            }
            n++;
         }
      }
   }

   if(*n2 == 0) {
     g_warning("no isosurfaces found.");
     return (GString*)0;
   }

   out = g_string_new("");
   for(i=0; i<nxm1; i++) {
      for(j=0; j<nym1; j++) {
         for(k=0; k<nzm1; k++) {
            cubeindex = 0;
            if (f[i  ][j  ][k  ] < isovalue) cubeindex |=   1;
            if (f[i+1][j  ][k  ] < isovalue) cubeindex |=   2;
            if (f[i+1][j  ][k+1] < isovalue) cubeindex |=   4;
            if (f[i  ][j  ][k+1] < isovalue) cubeindex |=   8;
            if (f[i  ][j+1][k  ] < isovalue) cubeindex |=  16;
            if (f[i+1][j+1][k  ] < isovalue) cubeindex |=  32;
            if (f[i+1][j+1][k+1] < isovalue) cubeindex |=  64;
            if (f[i  ][j+1][k+1] < isovalue) cubeindex |= 128;

            e = edgeTable[cubeindex];

            if (e == 0) continue;

            n = i*nyz3+j*nz3+k*3;

            if (e & 1   ) vertlist[0 ] = itab[n];
            if (e & 2   ) vertlist[1 ] = itab[n + nyz3 + 1];
            if (e & 4   ) vertlist[2 ] = itab[n + 3];
            if (e & 8   ) vertlist[3 ] = itab[n + 1];
            if (e & 16  ) vertlist[4 ] = itab[n + nz3];
            if (e & 32  ) vertlist[5 ] = itab[n + nyz3 + 1 + nz3];
            if (e & 64  ) vertlist[6 ] = itab[n + 3 + nz3];
            if (e & 128 ) vertlist[7 ] = itab[n + 1 + nz3];
            if (e & 256 ) vertlist[8 ] = itab[n + 2];
            if (e & 512 ) vertlist[9 ] = itab[n + nyz3 + 2];
            if (e & 1024) vertlist[10] = itab[n + nyz3 + 3 + 2];
            if (e & 2048) vertlist[11] = itab[n + 3 + 2];

            for (n = 0; triTable[cubeindex][n]!= -1; n += 3) {
               m1 = vertlist[triTable[cubeindex][n  ]];
               if(m1 == -1)
                 {
                   g_warning("m1 %d %d %d", i, j, k);
                   g_string_free(out, TRUE);
                   return (GString*)0;
                 }
               m2 = vertlist[triTable[cubeindex][n+1]];
               if(m2 == -1)
                 {
                   g_warning("m2 %d %d %d", i, j, k);
                   g_string_free(out, TRUE);
                   return (GString*)0;
                 }
               m3 = vertlist[triTable[cubeindex][n+2]];
               if(m3 == -1)
                 {
                   g_warning("m3 %d %d %d", i, j, k);
                   g_string_free(out, TRUE);
                   return (GString*)0;
                 }
               ux = xs[m2] - xs[m1];
               uy = ys[m2] - ys[m1];
               uz = zs[m2] - zs[m1];
               vx = xs[m3] - xs[m1];
               vy = ys[m3] - ys[m1];
               vz = zs[m3] - zs[m1];
               wx = uy*vz-uz*vy;
               wy = uz*vx-ux*vz;
               wz = ux*vy-uy*vx;
               s = sqrt(wx*wx+wy*wy+wz*wz);
               if(s <= 0.0) continue;
               g_string_append_printf(out, "3 %d %d %d\n", m1+1, m2+1, m3+1);
               *n1 += 1;
            }
         }
      }
   }

   if(*n1 == 0) {
     g_warning("no isosurfaces found.");
     return (GString*)0;
   }


   for(n=0; n<*n2; n++)
     g_string_append_printf(out, "%g %g %g  %g %g %g\n",
                            xs[n], ys[n], zs[n], xns[n], yns[n], zns[n]);

   return out;
}

/**
 * visu_surface_createFromPotentialFile:
 * @surf_file_to_write: target surf file
 * @pot_file_to_read: source pot file
 * @nsurfs_to_build: number of surfaces to build
 * @surf_value: an array with iso-values
 * @surf_name: an array with iso-surface names.
 *
 * Read the given pot file and produce an associated surf
 * file. WARNING, it may be removed later.
 *
 * Deprecated: 3.7
 *
 * Returns: 0 in case of success, n!=0 otherwise.
 */
#ifndef V_SIM_DISABLE_DEPRECATED
int visu_surface_createFromPotentialFile(const gchar *surf_file_to_write, const gchar *pot_file_to_read, int nsurfs_to_build,
		    const float *surf_value, const gchar **surf_name)
{
   register int i, j, k;
   char rep[TOOL_MAX_LINE_LENGTH], info[TOOL_MAX_LINE_LENGTH];
   double v;
   GString *out, *local;
   guint n1, n2;
   int n1_tot, n2_tot;
   gchar *comment;
   GError *error;

   n1_tot = 0;
   n2_tot = 0;

   in = fopen(pot_file_to_read,"r");
   if (!in) {
      (void)fprintf(stderr, "\aERROR: cannot read potential file %s\n", pot_file_to_read);
      return 1;
   }
   DBG_fprintf(stderr, "pot2surf : Reading potential file %s\n", pot_file_to_read);

   out = g_string_new("");
   DBG_fprintf(stderr, "pot2surf : Writing surface file %s ...\n", surf_file_to_write);

   /* 1st line (comment) */
   g_return_val_if_fail(fgets(rep, TOOL_MAX_LINE_LENGTH, in), 1);
   rep[strlen(rep)-1] = 0; /* skipping \n */
   DBG_fprintf(stderr, "Pot2surf : read comment line '%s'\n", rep);
   comment = g_strdup(rep);

   g_return_val_if_fail(fgets(rep, TOOL_MAX_LINE_LENGTH, in), 1);
   g_return_val_if_fail(sscanf(rep, "%d %d %d", &nx, &ny, &nz) != 3, 1);

   f = g_malloc(nx*sizeof(double **));
   for(i=0; i<nx; i++) {
      f[i] = g_malloc(ny*sizeof(double *));
      for(j=0; j<ny; j++) {
         f[i][j] = g_malloc(nz*sizeof(double));
      }
   }

   nxm1 = nx - 1;
   nym1 = ny - 1;
   nzm1 = nz - 1;

   nz3 = 3*nz;
   nyz3 = ny*nz3;
   nxyz3 = nx*nyz3;

   itab = g_malloc(nxyz3*sizeof(int));

   xs = g_malloc(nxyz3*sizeof(double));
   ys = g_malloc(nxyz3*sizeof(double));
   zs = g_malloc(nxyz3*sizeof(double));
   xns = g_malloc(nxyz3*sizeof(double));
   yns = g_malloc(nxyz3*sizeof(double));
   zns = g_malloc(nxyz3*sizeof(double));

   g_return_val_if_fail(fgets(rep, TOOL_MAX_LINE_LENGTH, in), 1);
   g_return_val_if_fail(sscanf(rep, "%lf %lf %lf", &dxx, &dyx, &dyy) != 3, 1);
   g_return_val_if_fail(fgets(rep, TOOL_MAX_LINE_LENGTH, in), 1);
   g_return_val_if_fail(sscanf(rep, "%lf %lf %lf", &dzx, &dzy, &dzz) != 3, 1);;

   v = dxx*dyy*dzz;
   exx = dyy*dzz/v;
   eyx = -dyx*dzz/v;
   ezx = (dyx*dzy-dyy*dzx)/v;
   eyy = dxx*dzz/v;
   ezy = -dxx*dzy/v;
   ezz = dxx*dyy/v;

   dxx1 = dxx/nx;
   dyx1 = dyx/ny;
   dyy1 = dyy/ny;
   dzx1 = dzx/nz;
   dzy1 = dzy/nz;
   dzz1 = dzz/nz;

   g_return_val_if_fail(fgets(rep, TOOL_MAX_LINE_LENGTH, in), 1);
   if(sscanf(rep, "%s", info) != 1) {
     (void)fprintf(stderr, "\aERROR: missing chain xyz OR zyx\n");
     return 1;
   }
   if(!strcmp(info, "xyz")) {
     for ( k = 0; k < nz; k++ )
       for ( j = 0; j < ny; j++ )
         for ( i = 0; i < nx; i++ )
           g_return_val_if_fail(fscanf(in, "%lf", &f[i][j][k]) != 1, 1);
   }
   else if(!strcmp(info, "zyx")) {
     for ( i = 0; i < nx; i++ )
       for ( j = 0; j < ny; j++ )
         for ( k = 0; k < nz; k++ )
           g_return_val_if_fail(fscanf(in, "%lf", &f[i][j][k]) != 1, 1);
   }
   else {
      (void)fprintf(stderr,
         "\aERROR: 5th line should contain chain xyz OR zyx\n");
      return 1;
   }

   for(k=0; k<nsurfs_to_build; k++) {

     isovalue = surf_value[k];

     DBG_fprintf(stderr, "Building isosurface for isovalue = %f (%s)\n",
                 isovalue, surf_name[k]);
     g_string_append_printf(out, "%s %f\n", ISOSURFACES_FLAG_POTENTIAL, isovalue);
     if (surf_name[k])
       g_string_append(out, surf_name[k]);
     g_string_append(out, "\n");

     n1 = 0;
     n2 = 0;
     local = calc(&n1, &n2);
     if(!local)
       return 1;
     DBG_fprintf(stderr, "  Found %d facets and %d points\n", n1, n2);
     g_string_append_printf(out, "%12d%12d", n1, n2);
     g_string_append(out, local->str);
     g_string_free(local, TRUE);

     n1_tot += n1;
     n2_tot += n2;
   }

   local = g_string_new(comment);
   g_free(comment);
   g_string_append_printf(local, "%f %f %f\n", dxx, dyx, dyy);
   g_string_append_printf(local, "%f %f %f\n", dzx, dzy, dzz);
   g_string_append_printf(local, "%12d%12d%12d", nsurfs_to_build, n1_tot, n2_tot);
   g_string_prepend(out, local->str);
   g_string_free(local, TRUE);

   DBG_fprintf(stderr,
      "Found a total of %d facets and %d points\n", n1_tot, n2_tot);
   DBG_fprintf(stderr, "Done.\n");

   for(i=0; i<nx; i++)
     {
       for(j=0; j<ny; j++)
	 g_free(f[i][j]);
       g_free(f[i]);
     }
   g_free(f);

   g_free(itab);

   g_free(xs);
   g_free(ys);
   g_free(zs);
   g_free(xns);
   g_free(yns);
   g_free(zns);

   error = (GError*)0;
   if (!g_file_set_contents(surf_file_to_write, out->str, -1, &error))
     {
       g_warning("%s", error->message);
       g_error_free(error);
       return 1;
     }
   g_string_free(out, TRUE);

   return 0;
}
#endif

/******************************************************************************/
/******************************************************************************/


/* DEVELOPMENT AREA
   Methods are copied to be callable without all static stuff everywhere. */

static double getDeltaMesh(int i, double *mesh, int size, gboolean periodic)
{
  if (i >= size - 1)
    {
      if (periodic)
	return 1. + mesh[(i + 1)%size] - mesh[i%size];
      else
	return (mesh[size] - mesh[size - 1]);
    }
  return (mesh[i + 1] - mesh[i]);
}

/**
 * visu_surface_new_fromScalarField:
 * @field: the scalar field to create the surface from ;
 * @isoValue: the value of the isosurface ;
 * @name: (allow-none): the name of the surface to use (can be NULL).
 *
 * Create on the fly a surface from the scalar field @field. If
 * @name is given, the surface is created with it, if not, "Isosurface
 * @id + 1" is used. @surf can already contains several surfaces, in
 * that case, the new surface is added. If @surf is NULL, then a new
 * #VisuSurface object is created and returned.
 *
 * Returns: (transfer full): the newly created surface, if any.
 */
VisuSurface* visu_surface_new_fromScalarField(const VisuScalarField *field,
                                                double isoValue,
                                                const gchar *name)
{
  VisuSurface *surf;

  surf = (VisuSurface*)0;
  switch (visu_scalar_field_getMeshtype(field))
    {
    case VISU_SCALAR_FIELD_MESH_UNIFORM:
      create_fromSF_uniform   (&surf, field, isoValue, name);
      break;
    case VISU_SCALAR_FIELD_MESH_NON_UNIFORM:
      create_fromSF_nonuniform(&surf, field, isoValue, name);
      break;
    default:
      g_warning("Wrong value for 'meshtype'.");
      return FALSE;
    }
  if (surf)
    g_object_set_data(G_OBJECT(surf), "origin", (gpointer)field);
  return surf;
}

/**
 * visu_surface_new_defaultFromScalarField:
 * @field: a #VisuScalarField object.
 * @neg: (transfer full) (out): a #VisuSurface location.
 * @pos: (transfer full) (out): a #VisuSurface location.
 *
 * Creates one or two surface from @field. If @field is mostly
 * positive, one surface is generated at the mean value. If @field has
 * equivalent negative and positive values, two surfaces are
 * generated, one for th mean negative value and one for the mean
 * positive value.
 *
 * Since: 3.8
 **/
void visu_surface_new_defaultFromScalarField(const VisuScalarField *field,
                                             VisuSurface **neg,
                                             VisuSurface **pos)
{
  double mimax[2];
  VisuSurface *surf;
  VisuSurfaceResource *res;
  gchar *label;
  static guint id = 0;
  int blue[4] = {0, 24, 185, 196}, red[4] = {185, 24, 0, 196};
  gboolean new;

  if (neg)
    *neg = (VisuSurface*)0;
  if (pos)
    *pos = (VisuSurface*)0;
  id += 1;
  visu_scalar_field_getMinMax(field, mimax);
  /* For densities, this create only one surface. For
     wavefunctions, it create two and set their style
     accordingly. */
  if (mimax[0] <= mimax[1] && mimax[0] * mimax[1] < 0.f &&
      MIN(ABS(mimax[0]), ABS(mimax[1])) / MAX(ABS(mimax[0]), ABS(mimax[1])) > 0.2)
    {
      label = g_strdup_printf(_("Negative (%d)"), id);
      res = visu_surface_resource_new_fromName(label, &new);
      if (new)
        g_object_set(res, "color", tool_color_addIntRGBA(blue),
                     "rendered", TRUE, NULL);
      g_object_unref(res);
      surf = visu_surface_new_fromScalarField(field, mimax[0] / 2., label);
      g_free(label);
      if (surf && neg)
        *neg = surf;
      else if (surf)
        g_object_unref(surf);

      label = g_strdup_printf(_("Positive (%d)"), id);
      res = visu_surface_resource_new_fromName(label, &new);
      if (new)
        g_object_set(res, "color", tool_color_addIntRGBA(red),
                     "rendered", TRUE, NULL);
      g_object_unref(res);
      surf = visu_surface_new_fromScalarField(field, mimax[1] / 2., label);
      g_free(label);
      if (surf && pos)
        *pos = surf;
      else if (surf)
        g_object_unref(surf);
    }
  else if (mimax[0] <= mimax[1])
    {
      surf = visu_surface_new_fromScalarField(field, (mimax[0] + mimax[1]) / 2.,
                                              (const gchar*)0);
      if (surf && pos)
        *pos = surf;
      else if (surf)
        g_object_unref(surf);
    }
  else
    {
      surf = visu_surface_new_fromScalarField(field, 0., (gchar*)0);
      if (surf && pos)
        *pos = surf;
      else if (surf)
        g_object_unref(surf);
    }
}

static void _setNormal(double nrm[3], double s, double vux, double vuy, double vuz)
{
  if (s <= 0.0)
    {
      g_warning("normal direction assumed");
      nrm[0] = 1.;
      nrm[1] = 0.;
      nrm[2] = 0.;
    }
  else
    {
      nrm[0] = vux/s;
      nrm[1] = vuy/s;
      nrm[2] = vuz/s;
    }
}

/**
 * create_fromSF_uniform:
 * @surf: a location on a #VisuSurface pointer ;
 * @field: the scalar field to create the surface from ;
 * @isoValue: the value of the isosurface ;
 * @id: an integer identifying the surface ;
 * @name: the name of the surface to use (can be NULL).
 *
 * Create on the fly a surface from the scalar field @field. If @name is given, the surface
 * is created with it, if not, "Isosurface @id + 1" is used. @surf can already
 * contains several surfaces, in that case, the new surface is added. If @surf is
 * NULL, then a new #VisuSurface object is created and returned.
 *
 * Returns: TRUE if the surface is created.
 */

static gboolean create_fromSF_uniform(VisuSurface **surf, const VisuScalarField *field,
                                      double isoValue, const gchar *name)
{

/*  I label local cell vertices as in ref. above
    and assigns x y z names like this :

                 y
                 |

                 4----------5
                /.         /|
               7----------6 |
               | .        | |
               | .        | |
               | 0........|.1  --x
               |/         |/
               3----------2

             /
            z

*/

/* Local variables */

  register guint i, j, k;
  int n;
  int nelez3, neleyz3, nelexyz3;
  guint sizem1[3]; /*nx-1, ny-1 ,nz-1 : number of mesh intervals in the x,y and z direction*/
  guint scalarGrid[3]; /*nx, ny ,nz : number of mesh nodes in the x,y and z direction*/
  int *iTab;
  double fac;
  double xu, yu, zu, xui, xuij, yuj;
  double vux, vuy, vuz, s;
  double vxnu, vynu, vznu, v;
  double orthoBox[6];  /* inverse length of the box*/
  double dBox[6];   /*step of the mesh*/
  GArray *points;
  VisuSurfacePoint at;
  double val1, val2;
  VisuBox *box;
  double scalarBox[VISU_BOX_N_VECTORS]; /*define the box edges and topology*/
  gboolean per[3], status;
  float scalarShift[3];

/* Routine code*/

/* Ensure that the pointer is not zero */
  g_return_val_if_fail(surf, FALSE);

  /* edges are labelled from each corner:
     x+ then z+ then y+
     for k varying fastest
     then j varying fastest
     then i varying fastest
  */

  /* Get data from the field. */
  box = visu_boxed_getBox(VISU_BOXED(field));
  for (i = 0; i < VISU_BOX_N_VECTORS; i++)
    scalarBox[i] = visu_box_getGeometry(box, i);
  visu_box_getPeriodicity(box, per);
  visu_scalar_field_getGridSize(field, scalarGrid);
  visu_pointset_getTranslation(VISU_POINTSET(field), scalarShift);

  DBG_fprintf(stderr, "Pot2surf : compute isosurface at value %g from"
	      " field %p.\n", isoValue, (gpointer)field);
  DBG_fprintf(stderr, " | periodic %d, %d, %d.\n", per[0], per[1], per[2]);
  DBG_fprintf(stderr, " | value %g.\n", isoValue);
  DBG_fprintf(stderr, " | size %d %d %d.\n", scalarGrid[0],
	      scalarGrid[1], scalarGrid[2]);
  DBG_fprintf(stderr, " | box %g %g %g\n", scalarBox[0],
	      scalarBox[1], scalarBox[2]);
  DBG_fprintf(stderr, " |     %g %g %g.\n", scalarBox[3],
	      scalarBox[4], scalarBox[5]);
  DBG_fprintf(stderr, " | shift %g %g %g\n", scalarShift[0],
	      scalarShift[1], scalarShift[2]);

  /* Compute some local values. */
  if (per[0])
    scalarGrid[0] += 1;
  if (per[1])
    scalarGrid[1] += 1;
  if (per[2])
    scalarGrid[2] += 1;
  for (i = 0; i < 3; i++)
    sizem1[i] = scalarGrid[i] - 1;
  v   = scalarBox[0] * scalarBox[2] * scalarBox[5];  /* volume of the box = Lx*Ly*Lz */
  orthoBox[0] = scalarBox[2] * scalarBox[5] / v;     /*  1/Lx */
  orthoBox[1] = -scalarBox[1] * scalarBox[5] / v;
  orthoBox[2] = scalarBox[0] * scalarBox[5] / v;     /*  1/Ly */
  orthoBox[3] = (scalarBox[1] * scalarBox[4] - scalarBox[2] * scalarBox[3]) / v;
  orthoBox[4] = -scalarBox[0] * scalarBox[4] / v;
  orthoBox[5] = scalarBox[0] * scalarBox[2] / v;     /*  1/Lz */
  dBox[0] = scalarBox[0] / sizem1[0];
  dBox[1] = scalarBox[1] / sizem1[1];
  dBox[2] = scalarBox[2] / sizem1[1];
  dBox[3] = scalarBox[3] / sizem1[2];
  dBox[4] = scalarBox[4] / sizem1[2];
  dBox[5] = scalarBox[5] / sizem1[2];
  nelez3 = 3 * scalarGrid[2];
  neleyz3 = scalarGrid[1] * nelez3;
  nelexyz3 = scalarGrid[0] * neleyz3; /* 3*nx*ny*nz */

  /* Allocate buffers. */
  points = g_array_sized_new(FALSE, FALSE, sizeof(VisuSurfacePoint), nelexyz3 / 10);
  /* Indexes. */
  iTab = g_malloc(nelexyz3 * sizeof(int));

  n = 0;
  for(i=0; i<scalarGrid[0]; i++) {
    xui = i*dBox[0];
    for(j=0; j<scalarGrid[1]; j++) {
      xuij = xui+j*dBox[1];
      yuj = j*dBox[2];
      for(k=0; k<scalarGrid[2]; k++) {
	xu = xuij+k*dBox[3] - scalarShift[0];
	yu = yuj+k*dBox[4] - scalarShift[1];
	zu = k*dBox[5] - scalarShift[2];               /* xu,yu,zu coordinate in the mesh corresponding to index i,j,k */
	vxnu = visu_scalar_field_getGradAt(field, i, j, k, TOOL_XYZ_X);
	vynu = visu_scalar_field_getGradAt(field, i, j, k, TOOL_XYZ_Y);
	vznu = visu_scalar_field_getGradAt(field, i, j, k, TOOL_XYZ_Z);
	val1 = visu_scalar_field_getAt(field, i, j, k); /* v(i,j,k) */
	if(i<sizem1[0]) {
	  val2 = visu_scalar_field_getAt(field, i + 1, j, k); /* v(i+1,j,k)*/
	  if((isoValue - val1 < 0.0 && isoValue - val2 >= 0.0) ||
	     (isoValue - val1 >= 0.0 && isoValue - val2 < 0.0)) {
	    fac = (isoValue - val1) / (val2 - val1);
	    iTab[n] = points->len;
            at.at[0] = xu + fac*dBox[0];  /* surface intersection
                                          coordinates*/
            at.at[1] = yu;
            at.at[2] = zu;
	    vux = vxnu + fac*(visu_scalar_field_getGradAt(field, i+1, j, k, TOOL_XYZ_X) - vxnu);
	    vuy = vynu + fac*(visu_scalar_field_getGradAt(field, i+1, j, k, TOOL_XYZ_Y) - vynu);
	    vuz = vznu + fac*(visu_scalar_field_getGradAt(field, i+1, j, k, TOOL_XYZ_Z) - vznu);
	    vux = vux*orthoBox[0]+vuy*orthoBox[1]+vuz*orthoBox[3];
	    vuy =                 vuy*orthoBox[2]+vuz*orthoBox[4];
	    vuz =                                 vuz*orthoBox[5];
	    s = sqrt(vux*vux+vuy*vuy+vuz*vuz);
            _setNormal(at.normal, s, vux, vuy, vuz);
            g_array_append_val(points, at);
	  }
	  else {
	    iTab[n] = -1;
	  }
	}
	else {
	  iTab[n] = -1;
	}
	n++;
	if(k<sizem1[2]) {
	  val2 = visu_scalar_field_getAt(field, i, j, k + 1);
	  if((isoValue - val1 < 0.0 && isoValue - val2 >= 0.0) ||
	     (isoValue - val1 >= 0.0 && isoValue - val2 < 0.0)) {
	    fac = (isoValue - val1) / (val2 - val1);
	    iTab[n] = points->len;
            at.at[0] = xu + fac*dBox[3];
            at.at[1] = yu + fac*dBox[4];
            at.at[2] = zu + fac*dBox[5];
	    vux = vxnu + fac*(visu_scalar_field_getGradAt(field, i, j, k+1, TOOL_XYZ_X) - vxnu);
	    vuy = vynu + fac*(visu_scalar_field_getGradAt(field, i, j, k+1, TOOL_XYZ_Y) - vynu);
	    vuz = vznu + fac*(visu_scalar_field_getGradAt(field, i, j, k+1, TOOL_XYZ_Z) - vznu);
	    vux = vux*orthoBox[0]+vuy*orthoBox[1]+vuz*orthoBox[3];
	    vuy =                 vuy*orthoBox[2]+vuz*orthoBox[4];
	    vuz =                                 vuz*orthoBox[5];
	    s = sqrt(vux*vux+vuy*vuy+vuz*vuz);
            _setNormal(at.normal, s, vux, vuy, vuz);
            g_array_append_val(points, at);
	  }
	  else {
	    iTab[n] = -1;
	  }
	}
	else {
	  iTab[n] = -1;
	}
	n++;
	if(j<sizem1[1]) {
	  val2 = visu_scalar_field_getAt(field, i, j + 1, k);
	  if((isoValue - val1 < 0.0 && isoValue - val2 >= 0.0) ||
	     (isoValue - val1 >= 0.0 && isoValue - val2 < 0.0)) {
	    fac = (isoValue - val1) / (val2 - val1);
	    iTab[n] = points->len;
            at.at[0] = xu + fac*dBox[1];
            at.at[1] = yu + fac*dBox[2];
            at.at[2] = zu;
	    vux = vxnu + fac*(visu_scalar_field_getGradAt(field, i, j+1, k, TOOL_XYZ_X) - vxnu);
	    vuy = vynu + fac*(visu_scalar_field_getGradAt(field, i, j+1, k, TOOL_XYZ_Y) - vynu);
	    vuz = vznu + fac*(visu_scalar_field_getGradAt(field, i, j+1, k, TOOL_XYZ_Z) - vznu);
	    vux = vux*orthoBox[0]+vuy*orthoBox[1]+vuz*orthoBox[3];
	    vuy =                 vuy*orthoBox[2]+vuz*orthoBox[4];
	    vuz =                                 vuz*orthoBox[5];
	    s = sqrt(vux*vux+vuy*vuy+vuz*vuz);
            _setNormal(at.normal, s, vux, vuy, vuz);
            g_array_append_val(points, at);
	  }
	  else {
	    iTab[n] = -1;
	  }
	}
	else {
	  iTab[n] = -1;
	}
	n++;
      }
    }
  }
  DBG_fprintf(stderr, " | found %d points.\n", points->len);

  if (points->len == 0)
    {
      g_warning("no isosurface found.");
      g_array_unref(points);
      g_free(iTab);
      return FALSE;
    }
  status = Create_surf(nelez3, neleyz3, nelexyz3, sizem1, iTab,
                       points, isoValue, field, scalarBox, name, surf);

  DBG_fprintf(stderr, " | Surface successfully created.\n");

  return status;
}


/**
 * create_fromSF_nonuniform:
 * @surf: a location on a #VisuSurface pointer ;
 * @field: the scalar field to create the surface from ;
 * @isoValue: the value of the isosurface ;
 * @id: an integer identifying the surface ;
 * @name: the name of the surface to use (can be NULL).
 *
 * Create on the fly a surface from the scalar field @field. If @name is given, the surface
 * is created with it, if not, "Isosurface @id + 1" is used. @surf can already
 * contains several surfaces, in that case, the new surface is added. If @surf is
 * NULL, then a new #VisuSurface object is created and returned.
 *
 * Returns: TRUE if the surface is created.
 */

static gboolean create_fromSF_nonuniform(VisuSurface **surf, const VisuScalarField *field,
                                       double isoValue, const gchar *name)
{

/*  I label local cell vertices as in ref. above
    and assigns x y z names like this :

                 y
                 |

                 4---------------5
                /.              /|
               7---------------6 |
               | .             | |
               | .             | |
               | 0.............|.1  --x
               |/              |/
               3---------------2

             /
            z

*/

/* Local variables */

  register guint i, j, k;
  int n;
  int nelez3, neleyz3, nelexyz3;
  guint sizem1[3]; /*nx-1, ny-1 ,nz-1 : number of mesh intervals in the x,y and z direction*/
  guint scalarGrid[3]; /*nx, ny ,nz : number of mesh nodes in the x,y and z direction*/
  int *iTab;
  double fac;
  double xu, yu, zu, xui, xuij, yuj, dmesh;
  double vux, vuy, vuz, s;
  double vxnu, vynu, vznu, v;
  double orthoBox[6];  /* inverse length of the box*/
  GArray *points;
  VisuSurfacePoint at;
  double *meshx, *meshy, *meshz;
  double val1, val2;
  VisuBox *box;
  double scalarBox[6]; /*define the box edges and topology*/
  gboolean per[3], status;

/* Routine code*/

/* Ensure that the pointer is not zero */
  g_return_val_if_fail(surf, FALSE);

  /* edges are labelled from each corner:
     x+ then z+ then y+
     for k varying fastest
     then j varying fastest
     then i varying fastest
  */

  /* Get data from the field. */
  box = visu_boxed_getBox(VISU_BOXED(field));
  for (i = 0; i < VISU_BOX_N_VECTORS; i++)
    scalarBox[i] = visu_box_getGeometry(box, i);
  visu_box_getPeriodicity(box, per);
  meshx = visu_scalar_field_getMesh(field, TOOL_XYZ_X);
  meshy = visu_scalar_field_getMesh(field, TOOL_XYZ_Y);
  meshz = visu_scalar_field_getMesh(field, TOOL_XYZ_Z);
  visu_scalar_field_getGridSize(field, scalarGrid);

  DBG_fprintf(stderr, "Pot2surf : compute isosurface at value %g from"
	      " field %p.\n", isoValue, (gpointer)field);
  DBG_fprintf(stderr, " | periodic %d, %d, %d.\n", per[0], per[1], per[2]);
  DBG_fprintf(stderr, " | value %g.\n", isoValue);
  DBG_fprintf(stderr, " | size %d %d %d.\n", scalarGrid[0],
	      scalarGrid[1], scalarGrid[2]);
  DBG_fprintf(stderr, " | box %g %g %g\n", scalarBox[0],
	      scalarBox[1], scalarBox[2]);
  DBG_fprintf(stderr, " |     %g %g %g.\n", scalarBox[3],
	      scalarBox[4], scalarBox[5]);

  /* Compute some local values. */
  if (per[0])
    scalarGrid[0] += 1;
  if (per[1])
    scalarGrid[1] += 1;
  if (per[2])
    scalarGrid[2] += 1;
  for (i = 0; i < 3; i++)
    sizem1[i] = scalarGrid[i] - 1;
  v   = scalarBox[0] * scalarBox[2] * scalarBox[5];  /* volume of the box = Lx*Ly*Lz */
  orthoBox[0] = scalarBox[2] * scalarBox[5] / v;     /*  1/Lx */
  orthoBox[1] = -scalarBox[1] * scalarBox[5] / v;
  orthoBox[2] = scalarBox[0] * scalarBox[5] / v;     /*  1/Ly */
  orthoBox[3] = (scalarBox[1] * scalarBox[4] - scalarBox[2] * scalarBox[3]) / v;
  orthoBox[4] = -scalarBox[0] * scalarBox[4] / v;
  orthoBox[5] = scalarBox[0] * scalarBox[2] / v;     /*  1/Lz */
  nelez3 = 3 * scalarGrid[2];
  neleyz3 = scalarGrid[1] * nelez3;
  nelexyz3 = scalarGrid[0] * neleyz3; /* 3*nx*ny*nz */

  /* Allocate buffers. */
  points = g_array_sized_new(FALSE, FALSE, sizeof(VisuSurfacePoint), nelexyz3 / 10);
  /* Indexes. */
  iTab = g_malloc(nelexyz3 * sizeof(int));

  n = 0;
  for(i=0; i<scalarGrid[0]; i++) {
    xui   = ((per[0] && i == sizem1[0])?1.:meshx[i]) * scalarBox[0];
    for(j=0; j<scalarGrid[1]; j++) {
      xuij   = xui  + ((per[1] && j == sizem1[1])?1.:meshy[j]) * scalarBox[1];
      yuj    = ((per[1] && j == sizem1[1])?1.:meshy[j]) * scalarBox[2];
      for(k=0; k<scalarGrid[2]; k++) {
        xu   = xuij  + ((per[2] && k == sizem1[2])?1.:meshz[k]) * scalarBox[3];
	yu   = yuj   + ((per[2] && k == sizem1[2])?1.:meshz[k]) * scalarBox[4];
	zu   = ((per[2] && k == sizem1[2])?1.:meshz[k]) * scalarBox[5];
	/* xu  ,yu  ,zu coordinate in the mesh corresponding to index i ,j ,k */
        vxnu = visu_scalar_field_getGradAt(field, i, j, k, TOOL_XYZ_X);
        vynu = visu_scalar_field_getGradAt(field, i, j, k, TOOL_XYZ_Y);
        vznu = visu_scalar_field_getGradAt(field, i, j, k, TOOL_XYZ_Z);
        val1 = visu_scalar_field_getAt(field, i, j, k);       /*  v(i,j,k) */
        if(i<sizem1[0]) {
          val2 = visu_scalar_field_getAt(field, i + 1, j, k); /* v(i+1,j,k)*/
          if((isoValue - val1 < 0.0 && isoValue - val2 >= 0.0) ||
             (isoValue - val1 >= 0.0 && isoValue - val2 < 0.0)) {
            fac = (isoValue - val1) / (val2 - val1);
            iTab[n] = points->len;
	    dmesh = getDeltaMesh(i, meshx, sizem1[0], per[0]);
            at.at[0] = xu + fac * dmesh * scalarBox[0];
            at.at[1] = yu;
            at.at[2] = zu;
            vux = vxnu + fac*(visu_scalar_field_getGradAt(field, i+1, j, k, TOOL_XYZ_X) - vxnu);
            vuy = vynu + fac*(visu_scalar_field_getGradAt(field, i+1, j, k, TOOL_XYZ_Y) - vynu);
            vuz = vznu + fac*(visu_scalar_field_getGradAt(field, i+1, j, k, TOOL_XYZ_Z) - vznu);
            vux = vux*orthoBox[0]+vuy*orthoBox[1]+vuz*orthoBox[3];
            vuy =                 vuy*orthoBox[2]+vuz*orthoBox[4];
            vuz =                                 vuz*orthoBox[5];
            s = sqrt(vux*vux+vuy*vuy+vuz*vuz);
            _setNormal(at.normal, s, vux, vuy, vuz);
            g_array_append_val(points, at);
	  }
	  else {
	    iTab[n] = -1;
	  }
	}
	else {
	  iTab[n] = -1;
	}
	n++;
	if(k<sizem1[2]) {
	  val2 = visu_scalar_field_getAt(field, i, j, k + 1);
	  if((isoValue - val1 < 0.0 && isoValue - val2 >= 0.0) ||
	     (isoValue - val1 >= 0.0 && isoValue - val2 < 0.0)) {
	    fac = (isoValue - val1) / (val2 - val1);
	    iTab[n] = points->len;
	    dmesh = getDeltaMesh(k, meshz, sizem1[2], per[2]);
            at.at[0] = xu + fac * dmesh * scalarBox[3];
            at.at[1] = yu + fac * dmesh * scalarBox[4];
            at.at[2] = zu + fac * dmesh * scalarBox[5];
	    vux = vxnu + fac*(visu_scalar_field_getGradAt(field, i, j, k+1, TOOL_XYZ_X) - vxnu);
	    vuy = vynu + fac*(visu_scalar_field_getGradAt(field, i, j, k+1, TOOL_XYZ_Y) - vynu);
	    vuz = vznu + fac*(visu_scalar_field_getGradAt(field, i, j, k+1, TOOL_XYZ_Z) - vznu);
	    vux = vux*orthoBox[0]+vuy*orthoBox[1]+vuz*orthoBox[3];
	    vuy =                 vuy*orthoBox[2]+vuz*orthoBox[4];
	    vuz =                                 vuz*orthoBox[5];
	    s = sqrt(vux*vux+vuy*vuy+vuz*vuz);
            _setNormal(at.normal, s, vux, vuy, vuz);
            g_array_append_val(points, at);
	  }
	  else {
	    iTab[n] = -1;
	  }
	}
	else {
	  iTab[n] = -1;
	}
	n++;
	if(j<sizem1[1]) {
	  val2 = visu_scalar_field_getAt(field, i, j + 1, k);
	  if((isoValue - val1 < 0.0 && isoValue - val2 >= 0.0) ||
	     (isoValue - val1 >= 0.0 && isoValue - val2 < 0.0)) {
	    fac = (isoValue - val1) / (val2 - val1);
	    iTab[n] = points->len;
	    dmesh = getDeltaMesh(j, meshy, sizem1[1], per[1]);
            at.at[0] = xu + fac * dmesh * scalarBox[1];
            at.at[1] = yu + fac * dmesh * scalarBox[2];
            at.at[2] = zu;
	    vux = vxnu + fac*(visu_scalar_field_getGradAt(field, i, j+1, k, TOOL_XYZ_X) - vxnu);
	    vuy = vynu + fac*(visu_scalar_field_getGradAt(field, i, j+1, k, TOOL_XYZ_Y) - vynu);
	    vuz = vznu + fac*(visu_scalar_field_getGradAt(field, i, j+1, k, TOOL_XYZ_Z) - vznu);
	    vux = vux*orthoBox[0]+vuy*orthoBox[1]+vuz*orthoBox[3];
	    vuy =                 vuy*orthoBox[2]+vuz*orthoBox[4];
	    vuz =                                 vuz*orthoBox[5];
	    s = sqrt(vux*vux+vuy*vuy+vuz*vuz);
            _setNormal(at.normal, s, vux, vuy, vuz);
            g_array_append_val(points, at);
	  }
	  else {
	    iTab[n] = -1;
	  }
	}
	else {
	  iTab[n] = -1;
	}
	n++;
      }
    }
  }
  DBG_fprintf(stderr, " | found %d points.\n", points->len);

  if (points->len == 0)
    {
      g_warning("no isosurface found.");
      g_array_unref(points);
      g_free(iTab);
      return FALSE;
    }
  status = Create_surf(nelez3, neleyz3, nelexyz3, sizem1, iTab,
                       points, isoValue, field, scalarBox, name, surf);

  DBG_fprintf(stderr, " | Surface successfully created.\n");

  return status;
}

/**
 * Create_surf:
 * @id: an integer identifying the surface ;
 * @nelez3: an integer giving the number of element times 3 along z;
 * @neleyz3: an integer giving the number of element times 3 in plane y by z;
 * @nelexyz3: an integer giving the number of element times 3 in volume x by y by z;
 * @nPoints: an integer giving the number points of the surface;
 * @sizem1: the number of intervals in each direction;
 * @iTab: table containing the number of surface points if surface exists at this place otherwise -1;
 * @xsurf: x coordinate of the surface points;
 * @ysurf: y coordinate of the surface points;
 * @zsurf: z coordinate of the surface points;
 * @xnsurf: normal value of the surface points;
 * @ynsurf: normal of the surface points;
 * @znsurf: normal of the surface points;
 * @isoValue: the value of the isosurface;
 * @scalarData: the potential values;
 * @scalarBox: the box parameters;
 * @name: the name of the surface to use (can be NULL).
 * @per: boolean for periodicity;
 * @surf: a location on a #VisuSurface pointer ;
 *
 * Create on the fly a surface from the scalar field @field. If @name is given, the surface
 * is created with it, if not, "Isosurface @id + 1" is used. @surf can already
 * contains several surfaces, in that case, the new surface is added. If @surf is
 * NULL, then a new #VisuSurface object is created and returned.
 *
 * Returns: TRUE if the surface is created.
 */

static gboolean Create_surf(int nelez3, int neleyz3, int nelexyz3 _U_,
                            guint sizem1[3], int *iTab, GArray *points,
                            double isoValue, const VisuScalarField *field, double scalarBox[6],
                            const gchar *name, VisuSurface **surf)
{

/* Local variables */

  register guint i, j, k;
  int cubeindex;
  int e, n;
  int m1, m2, m3;
  int vertlist[12];
  GArray *vTab;
  VisuSurfacePoly poly;
  VisuSurfacePoint *at1, *at2, *at3;
  double vx, vy, vz, wx, wy, wz;
  double ux, uy, uz, s;
  float* densityData;
  VisuBox *box;

  /* Routine code */

  /* Ensure that the pointer is not zero */
  g_return_val_if_fail(surf, FALSE);

  /* Allocate buffers. The 16 coefficient is for the maximum number of
     polygons per grid points. */
  vTab = g_array_sized_new(FALSE, FALSE, sizeof(VisuSurfacePoly), points->len / 3);
  poly.nvertices = 3;

  /* Main loop */
  for(i=0; i<sizem1[0]; i++) {
    for(j=0; j<sizem1[1]; j++) {
      for(k=0; k<sizem1[2]; k++) {
	cubeindex = 0;
	if (visu_scalar_field_getAt(field, i  , j  , k  ) <= isoValue) cubeindex |=   1;
	if (visu_scalar_field_getAt(field, i+1, j  , k  ) <= isoValue) cubeindex |=   2;
	if (visu_scalar_field_getAt(field, i+1, j  , k+1) <= isoValue) cubeindex |=   4;
	if (visu_scalar_field_getAt(field, i  , j  , k+1) <= isoValue) cubeindex |=   8;
	if (visu_scalar_field_getAt(field, i  , j+1, k  ) <= isoValue) cubeindex |=  16;
	if (visu_scalar_field_getAt(field, i+1, j+1, k  ) <= isoValue) cubeindex |=  32;
	if (visu_scalar_field_getAt(field, i+1, j+1, k+1) <= isoValue) cubeindex |=  64;
	if (visu_scalar_field_getAt(field, i  , j+1, k+1) <= isoValue) cubeindex |= 128;

	e = edgeTable[cubeindex];

	if (e == 0) continue;

	n = i*neleyz3+j*nelez3+k*3;

	if (e & 1   ) vertlist[0 ] = iTab[n];
	if (e & 2   ) vertlist[1 ] = iTab[n + neleyz3 + 1];
	if (e & 4   ) vertlist[2 ] = iTab[n + 3];
	if (e & 8   ) vertlist[3 ] = iTab[n + 1];
	if (e & 16  ) vertlist[4 ] = iTab[n + nelez3];
	if (e & 32  ) vertlist[5 ] = iTab[n + neleyz3 + 1 + nelez3];
	if (e & 64  ) vertlist[6 ] = iTab[n + 3 + nelez3];
	if (e & 128 ) vertlist[7 ] = iTab[n + 1 + nelez3];
	if (e & 256 ) vertlist[8 ] = iTab[n + 2];
	if (e & 512 ) vertlist[9 ] = iTab[n + neleyz3 + 2];
	if (e & 1024) vertlist[10] = iTab[n + neleyz3 + 3 + 2];
	if (e & 2048) vertlist[11] = iTab[n + 3 + 2];

	for (n = 0; triTable[cubeindex][n]!= -1; n += 3) {
	  m1 = vertlist[triTable[cubeindex][n  ]];
	  if(m1 == -1)
	    {
	      g_warning("m1 %d %d %d.", i, j, k);
              g_array_unref(points);
              g_free(iTab);
              g_array_unref(vTab);
	      return FALSE;
	    }
	  m2 = vertlist[triTable[cubeindex][n+1]];
	  if(m2 == -1)
	    {
	      g_warning("m2 %d %d %d.", i, j, k);
              g_array_unref(points);
              g_free(iTab);
              g_array_unref(vTab);
	      return FALSE;
	    }
	  m3 = vertlist[triTable[cubeindex][n+2]];
	  if(m3 == -1)
	    {
	      g_warning("m3 %d %d %d.", i, j, k);
              g_array_unref(points);
              g_free(iTab);
              g_array_unref(vTab);
	      return FALSE;
	    }
          at1 = &g_array_index(points, VisuSurfacePoint, m1);
          at2 = &g_array_index(points, VisuSurfacePoint, m2);
          at3 = &g_array_index(points, VisuSurfacePoint, m3);
	  ux = at2->at[0] - at1->at[0];
	  uy = at2->at[1] - at1->at[1];
	  uz = at2->at[2] - at1->at[2];
          vx = at3->at[0] - at1->at[0];
	  vy = at3->at[1] - at1->at[1];
	  vz = at3->at[2] - at1->at[2];
	  wx = uy*vz-uz*vy;
	  wy = uz*vx-ux*vz;
	  wz = ux*vy-uy*vx;
	  s = sqrt(wx*wx+wy*wy+wz*wz);
	  if(s <= 0.0) continue;
          poly.indices[0] = m1;
          poly.indices[1] = m2;
          poly.indices[2] = m3;
          g_array_append_val(vTab, poly);
	}
      }
    }
  }
  DBG_fprintf(stderr, " | found %d polygons.\n", vTab->len);

  if (vTab->len == 0)
    {
      g_warning("no isosurface found.");
      g_array_unref(points);
      g_free(iTab);
      g_array_unref(vTab);
      return FALSE;
    }

  /* Ok, try to create a VisuSurface object from here. */
  g_return_val_if_fail(*surf == (VisuSurface*)0, FALSE);

  DBG_fprintf(stderr, " | create new VisuSurface.\n");
  *surf = visu_surface_new(name, points, vTab);
  box = visu_box_new(scalarBox, VISU_BOX_PERIODIC);
  visu_boxed_setBox(VISU_BOXED(*surf), VISU_BOXED(box));
  g_object_unref(box);
  densityData = visu_surface_addPropertyFloat(*surf, VISU_SURFACE_PROPERTY_POTENTIAL);
  densityData[0] = isoValue;

#if DEBUG == 1
  visu_surface_checkConsistency(*surf);
#endif

  g_array_unref(points);
  g_array_unref(vTab);
  g_free(iTab);
  return TRUE;
}


/*****************************/
/* XML files for iso-values. */
/*****************************/

/*
   <surfaces>
     <surface rendered="yes" iso-value="0.45" name="foo">
       <hidden-by-planes status="yes" />
       <color rgba="0.5 0.5 0.5 1." material="0.2 1. 0.5 0.5 0." />
     </surface>
   </surfaces>
*/

/* Known elements. */
#define SURF_PARSER_ELEMENT_SURFACES  "surfaces"
#define SURF_PARSER_ELEMENT_SURFACE   "surface"
#define SURF_PARSER_ELEMENT_HIDE      "hidden-by-planes"
#define SURF_PARSER_ELEMENT_COLOR     "color"
/* Known attributes. */
#define SURF_PARSER_ATTRIBUTES_RENDERED "rendered"
#define SURF_PARSER_ATTRIBUTES_VALUE    "iso-value"
#define SURF_PARSER_ATTRIBUTES_NAME     "name"
#define SURF_PARSER_ATTRIBUTES_STATUS   "status"
#define SURF_PARSER_ATTRIBUTES_RGBA     "rgba"
#define SURF_PARSER_ATTRIBUTES_MATERIAL "material"

struct _surfaces_xml
{
  gchar *name;
  float iso;
  gboolean rendered, masked;
  gboolean colorSet, materialSet;
  float color[4], material[5];
};

static gboolean startVisuSurface;

/* This method is called for every element that is parsed.
   The user_data must be a GList of _surfaces_xml. When a 'surface'
   element, a new struct instance is created and prepend in the list.
   When 'hidden-by-planes' or other qualificative elements are
   found, the first surface of the list is modified accordingly. */
static void surfacesXML_element(GMarkupParseContext *context _U_,
                          const gchar         *element_name,
                          const gchar        **attribute_names,
                          const gchar        **attribute_values,
                          gpointer             user_data,
                          GError             **error)
{
  GList **surfacesList;
  struct _surfaces_xml *surf;
  int i;

  g_return_if_fail(user_data);
  surfacesList = (GList **)user_data;

  DBG_fprintf(stderr, "Surf parser: found '%s' element.\n", element_name);
  if (!strcmp(element_name, SURF_PARSER_ELEMENT_SURFACES))
    {
      /* Should have no attributes. */
      if (attribute_names[0])
	{
	  g_set_error(error, G_MARKUP_ERROR, G_MARKUP_ERROR_UNKNOWN_ATTRIBUTE,
		      _("Unexpected attribute '%s' for element '%s'."),
		      attribute_names[0], SURF_PARSER_ELEMENT_SURFACES);
	  return;
	}
      /* Initialise the surfacesList. */
      if (*surfacesList)
	{
	  g_set_error(error, G_MARKUP_ERROR, G_MARKUP_ERROR_PARSE,
		      _("DTD error: element '%s' should appear only once."),
		      SURF_PARSER_ELEMENT_SURFACES);
	  return;
	}
      *surfacesList = (GList*)0;
      startVisuSurface = TRUE;
    }
  else if (!strcmp(element_name, SURF_PARSER_ELEMENT_SURFACE))
    {
      surf = g_malloc(sizeof(struct _surfaces_xml));
      surf->name        = (gchar*)0;
      surf->rendered    = TRUE;
      surf->masked      = TRUE;
      surf->iso         = 12345.6789f;
      surf->colorSet    = FALSE;
      surf->materialSet = FALSE;

      /* We parse the attributes. */
      for (i = 0; attribute_names[i]; i++)
	{
	  if (!strcmp(attribute_names[i], SURF_PARSER_ATTRIBUTES_NAME))
	    surf->name = g_strdup(attribute_values[i]);
	  else if (!strcmp(attribute_names[i], SURF_PARSER_ATTRIBUTES_RENDERED))
	    {
	      if (!strcmp(attribute_values[i], "yes"))
		surf->rendered = TRUE;
	      else if (!strcmp(attribute_values[i], "no"))
		surf->rendered = FALSE;
	      else
		g_set_error(error, G_MARKUP_ERROR, G_MARKUP_ERROR_INVALID_CONTENT,
			    _("Invalid value '%s' for attribute '%s'."),
			    attribute_values[i], SURF_PARSER_ATTRIBUTES_RENDERED);
	    }
	  else if (!strcmp(attribute_names[i], SURF_PARSER_ATTRIBUTES_VALUE))
	    {
	      if (!(sscanf(attribute_values[i], "%f", &surf->iso) == 1))
		g_set_error(error, G_MARKUP_ERROR, G_MARKUP_ERROR_INVALID_CONTENT,
			    _("Invalid value '%s' for attribute '%s'."),
			    attribute_values[i], SURF_PARSER_ATTRIBUTES_VALUE);
	    }
	  else
	    g_set_error(error, G_MARKUP_ERROR, G_MARKUP_ERROR_UNKNOWN_ATTRIBUTE,
			_("Unexpected attribute '%s' for element '%s'."),
			attribute_names[i], SURF_PARSER_ELEMENT_SURFACE);
	  if (*error)
	    {
	      g_free(surf->name);
	      g_free(surf);
	      return;
	    }
	}
      /* We check the mandatory attribute. */
      if (surf->iso == 12345.6789f)
	{
	  g_set_error(error, G_MARKUP_ERROR, G_MARKUP_ERROR_PARSE,
		      _("Missing attribute '%s' for element '%s'."),
		      SURF_PARSER_ATTRIBUTES_VALUE, SURF_PARSER_ELEMENT_SURFACE);
	  g_free(surf->name);
	  g_free(surf);
	  return;
	}
      DBG_fprintf(stderr, "Surf parser: add a new surface '%s' %f %d.\n",
		  (surf->name)?surf->name:"None", surf->iso, surf->rendered);
      *surfacesList = g_list_prepend(*surfacesList, (gpointer)surf);
    }
  else if (startVisuSurface && !strcmp(element_name, SURF_PARSER_ELEMENT_HIDE))
    {
      if (!*surfacesList)
	{
	  g_set_error(error, G_MARKUP_ERROR, G_MARKUP_ERROR_INVALID_CONTENT,
		      _("DTD error: parent element '%s' of element '%s' is missing."),
		      SURF_PARSER_ELEMENT_SURFACE, SURF_PARSER_ELEMENT_HIDE);
	  return;
	}
      surf = (struct _surfaces_xml*)((*surfacesList)->data);

      /* We read the mandatory attribute. */
      if (attribute_names[0])
	{
	  if (!strcmp(attribute_names[0], SURF_PARSER_ATTRIBUTES_STATUS))
	    {
	      if (!strcmp(attribute_values[0], "yes"))
		surf->masked = TRUE;
	      else if (!strcmp(attribute_values[0], "no"))
		surf->masked = FALSE;
	      else
		g_set_error(error, G_MARKUP_ERROR, G_MARKUP_ERROR_INVALID_CONTENT,
			    _("Invalid value '%s' for attribute '%s'."),
			    attribute_values[0], SURF_PARSER_ATTRIBUTES_STATUS);
	    }
	  else
	    g_set_error(error, G_MARKUP_ERROR, G_MARKUP_ERROR_UNKNOWN_ATTRIBUTE,
			_("Unexpected attribute '%s' for element '%s'."),
			attribute_names[0], SURF_PARSER_ELEMENT_HIDE);
	}
      else
	g_set_error(error, G_MARKUP_ERROR, G_MARKUP_ERROR_PARSE,
		    _("Missing attribute '%s' for element '%s'."),
		    SURF_PARSER_ATTRIBUTES_VALUE, SURF_PARSER_ELEMENT_SURFACE);
      if (*error)
	return;

    }
  else if (startVisuSurface && !strcmp(element_name, SURF_PARSER_ELEMENT_COLOR))
    {
      if (!*surfacesList)
	{
	  g_set_error(error, G_MARKUP_ERROR, G_MARKUP_ERROR_INVALID_CONTENT,
		      _("DTD error: parent element '%s' of element '%s' is missing."),
		      SURF_PARSER_ELEMENT_SURFACE, SURF_PARSER_ELEMENT_COLOR);
	  return;
	}
      surf = (struct _surfaces_xml*)((*surfacesList)->data);

      for(i = 0; attribute_names[i]; i++)
	{
	  if (!strcmp(attribute_names[i], SURF_PARSER_ATTRIBUTES_RGBA))
	    {
	      if (sscanf(attribute_values[i], "%g %g %g %g",
			 surf->color, surf->color + 1,
			 surf->color + 2, surf->color + 3) != 4)
		g_set_error(error, G_MARKUP_ERROR, G_MARKUP_ERROR_INVALID_CONTENT,
			    _("Invalid value '%s' for attribute '%s'."),
			    attribute_values[i], SURF_PARSER_ATTRIBUTES_RGBA);
	      else
		surf->colorSet = TRUE;
	    }
	  else if (!strcmp(attribute_names[i], SURF_PARSER_ATTRIBUTES_MATERIAL))
	    {
	      if (sscanf(attribute_values[i], "%g %g %g %g %g",
			 surf->material, surf->material + 1, surf->material + 2,
			 surf->material + 3, surf->material + 4) != 5)
		g_set_error(error, G_MARKUP_ERROR, G_MARKUP_ERROR_INVALID_CONTENT,
			    _("Invalid value '%s' for attribute '%s'."),
			    attribute_values[i], SURF_PARSER_ATTRIBUTES_MATERIAL);
	      else
		surf->materialSet = TRUE;
	    }
	  else
	    g_set_error(error, G_MARKUP_ERROR, G_MARKUP_ERROR_UNKNOWN_ATTRIBUTE,
			_("Unexpected attribute '%s' for element '%s'."),
			attribute_names[i], SURF_PARSER_ELEMENT_COLOR);
	  if (*error)
	    return;
	}
    }
  else if (startVisuSurface)
    {
      /* We silently ignore the element if surfacesList is unset, but
	 raise an error if surfacesList has been set. */
      g_set_error(error, G_MARKUP_ERROR, G_MARKUP_ERROR_UNKNOWN_ELEMENT,
		  _("Unexpected element '%s'."), element_name);
    }
}

/* Check when a element is closed that everything required has been set. */
void surfacesXML_end(GMarkupParseContext *context _U_,
		      const gchar         *element_name,
		      gpointer             user_data _U_,
		      GError             **error _U_)
{
  if (!strcmp(element_name, SURF_PARSER_ELEMENT_SURFACES))
    startVisuSurface = FALSE;
}

/* What to do when an error is raised. */
static void surfacesXML_error(GMarkupParseContext *context _U_,
			GError              *error,
			gpointer             user_data)
{
  GList *tmpLst;

  DBG_fprintf(stderr, "Surf parser: error raised '%s'.\n", error->message);
  g_return_if_fail(user_data);

  /* We free the current list of surfaces. */
  tmpLst = *(GList**)user_data;
  while (tmpLst)
    {
      g_free(((struct _surfaces_xml*)tmpLst->data)->name);
      g_free(tmpLst->data);
      tmpLst = g_list_next(tmpLst);
    }
  g_list_free(*(GList**)user_data);
}

/**
 * visu_surface_parseXMLFile:
 * @filename: a path to a file.
 * @surfaces: (inout) (element-type VisuSurface*): a location on a
 * #VisuSurface pointer.
 * @field: the scalar field to create the surface from.
 * @error: a location to store a possible error.
 *
 * Parse the given XML file, looking for the &lt;surfaces&gt; tag and create
 * the given surfaces.
 *
 * Returns: FALSE if a error occured.
 */
gboolean visu_surface_parseXMLFile(const gchar* filename, GList **surfaces,
			       VisuScalarField *field, GError **error)
{
  GMarkupParseContext* xmlContext;
  GMarkupParser parser;
  gboolean status;
  gsize size;
  gchar *buffer;
  GList *list, *tmpLst;
  VisuSurface *surface;
  struct _surfaces_xml *surf;
  VisuSurfaceResource *res;

  g_return_val_if_fail(filename && surfaces && field, FALSE);

  buffer = (gchar*)0;
  if (!g_file_get_contents(filename, &buffer, &size, error))
    return FALSE;

  /* Create context. */
  list = (GList*)0;
  parser.start_element = surfacesXML_element;
  parser.end_element   = surfacesXML_end;
  parser.text          = NULL;
  parser.passthrough   = NULL;
  parser.error         = surfacesXML_error;
  xmlContext = g_markup_parse_context_new(&parser, 0, &list, NULL);

  /* Parse data. */
  startVisuSurface = FALSE;
  status = g_markup_parse_context_parse(xmlContext, buffer, size, error);

  /* Free buffers. */
  g_markup_parse_context_free(xmlContext);
  g_free(buffer);

  if (!status)
    return FALSE;

  if (!list)
    {
      *error = g_error_new(G_MARKUP_ERROR, G_MARKUP_ERROR_EMPTY,
			  _("No iso-value found."));
      return FALSE;
    }

  /* Need to reverse the list since elements have been prepended. */
  list = g_list_reverse(list);

  /* Convert the list to new isosurfaces. */
  DBG_fprintf(stderr, "Surf parser: create %d new surfaces for %p.\n",
	      g_list_length(list), (gpointer)(*surfaces));
  tmpLst = list;
  while(tmpLst)
    {
      surf = (struct _surfaces_xml*)tmpLst->data;
      surface = visu_surface_new_fromScalarField(field, surf->iso, surf->name);
      if (surface)
	{
	  res = visu_surface_getResource(surface);
          g_object_set(G_OBJECT(res), "rendered", surf->rendered,
                       "maskable", surf->masked, NULL);
	  if (surf->colorSet)
            g_object_set(G_OBJECT(res), "color", tool_color_addFloatRGBA(surf->color, NULL), NULL);
	  if (surf->materialSet)
            g_object_set(G_OBJECT(res), "material", surf->material, NULL);
          *surfaces = g_list_append(*surfaces, surface);
	}
      g_free(surf->name);
      g_free(surf);
      tmpLst = g_list_next(tmpLst);
    }
  g_list_free(list);

  return TRUE;
}

/**
 * visu_surface_exportXMLFile:
 * @filename: a path to a file.
 * @values: an array of @n values.
 * @res: an array of @n #VisuSurfaceResource.
 * @n: number of surface resources to export.
 * @error: a location to store a possible error.
 *
 * Export the surface resources into an XML file.
 *
 * Returns: FALSE if a error occured.
 */
gboolean visu_surface_exportXMLFile(const gchar* filename, float *values,
                                     VisuSurfaceResource **res, int n, GError **error)
{
  gboolean valid;
  GString *output;
  int i;
  gchar *lbl;
  gboolean rendered, maskable;
  ToolColor *color;
  float *mat;

  /* We actually output the values. */
  output = g_string_new("  <surfaces>\n");
  for (i = 0; i < n; i++)
    {
      g_object_get(G_OBJECT(res[i]), "label", &lbl, "color", &color,
                   "material", &mat, "rendered", &rendered,
                   "maskable", &maskable, NULL);
      g_string_append_printf(output, "    <surface rendered=\"%s\" iso-value=\"%g\"",
			     (rendered)?"yes":"no", values[i]);
      if (lbl && lbl[0])
	g_string_append_printf(output, " name=\"%s\"", lbl);
      g_string_append(output, ">\n");
      g_string_append_printf(output, "      <hidden-by-planes status=\"%s\" />\n",
			     (maskable)?"yes":"no");
      g_string_append_printf(output, "      <color rgba=\"%g %g %g %g\""
			     " material=\"%g %g %g %g %g\" />\n",
			     color->rgba[0], color->rgba[1],
			     color->rgba[2], color->rgba[3],
			     mat[0], mat[1], mat[2], mat[3], mat[4]);
      g_string_append(output, "    </surface>\n");
    }
  g_string_append(output, "  </surfaces>");

  valid = tool_XML_substitute(output, filename, "surfaces", error);
  if (!valid)
    {
      g_string_free(output, TRUE);
      return FALSE;
    }

  valid = g_file_set_contents(filename, output->str, -1, error);
  g_string_free(output, TRUE);
  return valid;
}