diff --git a/examples/incompressible-navier-stokes/tnlNSFastBuildConfig.h b/examples/incompressible-navier-stokes/tnlNSFastBuildConfig.h
new file mode 100644
index 0000000000000000000000000000000000000000..fc03ebd1ff9a5f274f89e92d99cc1b72c7c29df7
--- /dev/null
+++ b/examples/incompressible-navier-stokes/tnlNSFastBuildConfig.h
@@ -0,0 +1,71 @@
+/***************************************************************************
+                          tnlNSFastBuildConfig.h  -  description
+                             -------------------
+    begin                : Jul 7, 2014
+    copyright            : (C) 2014 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef TNLNSFASTBUILDCONFIG_H_
+#define TNLNSFASTBUILDCONFIG_H_
+
+#include <solvers/tnlConfigTags.h>
+
+class tnlNSFastBuildConfig
+{
+   public:
+
+      static void print() { cerr << "tnlNSFastBuildConfig" << endl; }
+};
+
+/****
+ * Turn off support for float and long double.
+ */
+template<> struct tnlConfigTagReal< tnlNSFastBuildConfig, float > { enum { enabled = false }; };
+template<> struct tnlConfigTagReal< tnlNSFastBuildConfig, long double > { enum { enabled = false }; };
+
+/****
+ * Turn off support for short int and long int indexing.
+ */
+template<> struct tnlConfigTagIndex< tnlNSFastBuildConfig, short int >{ enum { enabled = false }; };
+template<> struct tnlConfigTagIndex< tnlNSFastBuildConfig, long int >{ enum { enabled = false }; };
+
+/****
+ * 1, 2, and 3 dimensions are enabled by default
+ */
+template<> struct tnlConfigTagDimensions< tnlNSFastBuildConfig, 1 >{ enum { enabled = false }; };
+template<> struct tnlConfigTagDimensions< tnlNSFastBuildConfig, 2 >{ enum { enabled = true }; };
+template<> struct tnlConfigTagDimensions< tnlNSFastBuildConfig, 3 >{ enum { enabled = false }; };
+
+/****
+ * Use of tnlGrid is enabled for allowed dimensions and Real, Device and Index types.
+ */
+template< int Dimensions, typename Real, typename Device, typename Index >
+   struct tnlConfigTagMesh< tnlNSFastBuildConfig, tnlGrid< Dimensions, Real, Device, Index > >
+      { enum { enabled = tnlConfigTagDimensions< tnlNSFastBuildConfig, Dimensions >::enabled  &&
+                         tnlConfigTagReal< tnlNSFastBuildConfig, Real >::enabled &&
+                         tnlConfigTagDevice< tnlNSFastBuildConfig, Device >::enabled &&
+                         tnlConfigTagIndex< tnlNSFastBuildConfig, Index >::enabled }; };
+
+/****
+ * Please, chose your preferred time discretisation  here.
+ */
+template<> struct tnlConfigTagTimeDiscretisation< tnlNSFastBuildConfig, tnlExplicitTimeDiscretisationTag >{ enum { enabled = true }; };
+template<> struct tnlConfigTagTimeDiscretisation< tnlNSFastBuildConfig, tnlSemiImplicitTimeDiscretisationTag >{ enum { enabled = true }; };
+template<> struct tnlConfigTagTimeDiscretisation< tnlNSFastBuildConfig, tnlImplicitTimeDiscretisationTag >{ enum { enabled = false }; };
+
+/****
+ * Only the Runge-Kutta-Merson solver is enabled by default.
+ */
+template<> struct tnlConfigTagExplicitSolver< tnlNSFastBuildConfig, tnlExplicitEulerSolverTag >{ enum { enabled = false }; };
+
+#endif /* TNLNSFASTBUILDCONFIG_H_ */
diff --git a/examples/incompressible-navier-stokes/tnlPoissonProblem.h b/examples/incompressible-navier-stokes/tnlPoissonProblem.h
new file mode 100644
index 0000000000000000000000000000000000000000..74d2ff515eb46f0f07772c9509cfd4134835ee5d
--- /dev/null
+++ b/examples/incompressible-navier-stokes/tnlPoissonProblem.h
@@ -0,0 +1,110 @@
+#ifndef TNLPOISSONPROBLEM_H
+#define	TNLPOISSONPROBLEM_H
+
+#include <core/vectors/tnlVector.h>
+#include <mesh/tnlGrid.h>
+
+template< typename Mesh,
+          typename Real = typename Mesh::RealType,
+          typename Index = typename Mesh::IndexType >
+class tnlPoissonProblem
+{
+ 
+};
+
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+		  typename Index,
+		  typename MatrixType>
+class tnlPoissonProblem< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >
+{
+   public: 
+   
+   typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;
+   typedef typename MeshType::CoordinatesType CoordinatesType;
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+
+   static tnlString getType()
+   {
+	  return tnlString( "tnlPoisssonMatrix< " ) +
+			 MeshType::getType() + ", " +
+			 ::getType< Real >() + ", " +
+			 ::getType< Index >() + " >";
+   }
+   
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   Real getValue( const MeshType& mesh,
+                  const IndexType cellIndex,
+                  const CoordinatesType& coordinates,
+                  const Vector& u,
+                  const Real& time ) const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   Index getLinearSystemRowLength( const MeshType& mesh,
+                                   const IndexType& index,
+								   const CoordinatesType& coordinates ) const
+   {
+	   return 5;
+   }
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+	  void initLinearSystem(
+                               const MeshType& mesh,
+                               const IndexType& index,
+							   )
+   {
+	   _matrix.setElement(index,index,4);
+	   _matrix.setElement(index,mesh.getCellYPredecessor( index ),-1);
+	   _matrix.setElement(index,mesh.getCellXPredecessor( index ),-1);
+	   _matrix.setElement(index,mesh.getCellXSuccessor( index ),-1);
+	   _matrix.setElement(index,mesh.getCellYSuccessor( index ),-1);
+
+	  /*matrixRow.setElement( 0, mesh.getCellYPredecessor( index ), -1 );
+	  matrixRow.setElement( 1, mesh.getCellXPredecessor( index ), -1 );
+	  matrixRow.setElement( 2, index,                            4.0 );
+	  matrixRow.setElement( 3, mesh.getCellXSuccessor( index ),   -1 );
+	  matrixRow.setElement( 4, mesh.getCellYSuccessor( index ),   -1 );*/
+   }
+
+   void init(const MeshType& mesh)
+   {
+	   typename MatrixType::RowLengthsVector rowLenghts;
+	   IndexType N = mesh.getNumberOfCells();
+	   rowLenghts.setSize(N);
+	   rowLenghts.setValue(5);
+	   _matrix.setDimensions( N, N );
+	   _matrix.setRowLengths(rowLenghts);
+
+	   for (IndexType i = 0; i < N; i++)
+		   initLinearSystem(mesh, i);
+
+   }
+   template< typename Vector, typename Solver >
+   void solve(Solver & solver, const Vector & rhs, Vector & result)
+   {
+	   solver.setMatrix(_matrix);
+	   solver.solve(rhs,result);
+   }
+
+   MatrixType _matrix;
+};
+
+
+
+
+#include "tnlPoissonProblem_impl.h"
+
+
+#endif	/* TNLPOISSONPROBLEM_H */
diff --git a/examples/incompressible-navier-stokes/tnlPoissonProblem_impl.h b/examples/incompressible-navier-stokes/tnlPoissonProblem_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..43f8c29d951edc1a6aaf05c39191f965ff473838
--- /dev/null
+++ b/examples/incompressible-navier-stokes/tnlPoissonProblem_impl.h
@@ -0,0 +1,37 @@
+
+#ifndef TNLPOISSONPROBLEM_IMP_H
+#define	TNLPOISSONPROBLEM_IMP_H
+
+#include "tnlPoissonProblem.h"
+#include <mesh/tnlGrid.h>
+
+
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlPoissonProblem< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getValue( const MeshType& mesh,
+          const IndexType cellIndex,
+          const CoordinatesType& coordinates,
+          const Vector& u,
+          const Real& time ) const
+{
+   return ( u[ mesh.getCellXPredecessor( cellIndex ) ]
+            - 2.0 * u[ cellIndex ]
+            + u[ mesh.getCellXSuccessor( cellIndex ) ] ) * mesh.getHxSquareInverse() +
+           ( u[ mesh.getCellYPredecessor( cellIndex ) ]
+             - 2.0 * u[ cellIndex ]
+             + u[ mesh.getCellYSuccessor( cellIndex ) ] ) * mesh.getHySquareInverse();
+}
+
+
+
+#endif	/* TNLPOISSONPROBLEM_IMP_H */
diff --git a/examples/incompressible-navier-stokes/visit_writer.cpp b/examples/incompressible-navier-stokes/visit_writer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c27478eacd5849a2cd2c8425007c5c27cff3977a
--- /dev/null
+++ b/examples/incompressible-navier-stokes/visit_writer.cpp
@@ -0,0 +1,1066 @@
+/* ************************************************************************* //
+//                             visit_writer.c                                //
+// ************************************************************************* */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "visit_writer.h" //mozna uvozovky
+
+
+/*
+ * Globals.
+ */
+
+static FILE *fp = NULL;
+static int useBinary = 0;
+static int numInColumn = 0;
+
+
+/* ****************************************************************************
+ *  Function: end_line
+ *
+ *  Purpose:
+ *      If floats or ints have been written using the write_float or write_int
+ *      functions, this will issue a newline (if necessary) so that a new
+ *      heading can be placed.
+ *
+ *  Programmer: Hank Childs
+ *  Creation:   September 3, 2004
+ * 
+ * ************************************************************************* */
+
+static void end_line(void)
+{
+    if (!useBinary)
+    {
+        char str2[8] = "\n";
+        fprintf(fp, str2);
+        numInColumn = 0;
+    }
+}
+
+
+/* ****************************************************************************
+ *  Function: open_file
+ *
+ *  Purpose:
+ *      Opens a file for writing and assigns the handle to the global variable
+ *      "fp".
+ *
+ *  Programmer: Hank Childs
+ *  Creation:   September 3, 2004
+ * 
+ * ************************************************************************* */
+
+static void open_file(const char *filename)
+{
+    char full_filename[1024];
+    if (strstr(filename, ".vtk") != NULL)
+    {
+        strcpy(full_filename, filename);
+    }
+    else
+    {
+        sprintf(full_filename, "%s.vtk", filename);
+    }
+
+    fp = fopen(full_filename, "w+");
+}
+
+
+/* ****************************************************************************
+ *  Function: close_file
+ *
+ *  Purpose:
+ *      Closes the file with handle "fp" (a global variable).
+ *
+ *  Programmer: Hank Childs
+ *  Creation:   September 3, 2004
+ * 
+ * ************************************************************************* */
+
+static void close_file(void)
+{
+    end_line();
+    fclose(fp);
+    fp = NULL;
+}
+
+
+/* ****************************************************************************
+ *  Function: force_big_endian
+ *
+ *  Purpose:
+ *      Determines if the machine is little-endian.  If so, then, for binary
+ *      data, it will force the data to be big-endian.
+ *
+ *  Note:       This assumes that all inputs are 4 bytes long.
+ *
+ *  Programmer: Hank Childs
+ *  Creation:   September 3, 2004
+ * 
+ * ************************************************************************* */
+
+static void force_big_endian(unsigned char *bytes)
+{
+    static int doneTest = 0;
+    static int shouldSwap = 0;
+    if (!doneTest)
+    {
+        int tmp1 = 1;
+        unsigned char *tmp2 = (unsigned char *) &tmp1;
+        if (*tmp2 != 0)
+            shouldSwap = 1;
+        doneTest = 1;
+    }
+
+    if (shouldSwap & useBinary)
+    {
+        unsigned char tmp = bytes[0];
+        bytes[0] = bytes[3];
+        bytes[3] = tmp;
+        tmp = bytes[1];
+        bytes[1] = bytes[2];
+        bytes[2] = tmp;
+    }
+}
+
+/* ****************************************************************************
+ *  Function: force_double_big_endian
+ *
+ *  Purpose:
+ *      Determines if the machine is little-endian.  If so, then, for binary
+ *      data, it will force the data to be big-endian.
+ *
+ *  Note:       This assumes that all inputs are 8 bytes long.
+ *
+ *  Programmer: Hank Childs
+ *  Creation:   September 3, 2004
+ *
+ * ************************************************************************* */
+
+static void force_double_big_endian(unsigned char *bytes)
+{
+    static int doneTest = 0;
+    static int shouldSwap = 0;
+    if (!doneTest)
+    {
+        int tmp1 = 1;
+        unsigned char *tmp2 = (unsigned char *) &tmp1;
+        if (*tmp2 != 0)
+            shouldSwap = 1;
+        doneTest = 1;
+    }
+
+    if (shouldSwap & useBinary)
+    {
+        unsigned char tmp = bytes[0];
+        bytes[0] = bytes[7];
+        bytes[7] = tmp;
+        tmp = bytes[1];
+        bytes[1] = bytes[6];
+        bytes[6] = tmp;
+        tmp = bytes[2];
+        bytes[2] = bytes[5];
+        bytes[5] = tmp;
+        tmp = bytes[3];
+        bytes[3] = bytes[4];
+        bytes[4] = tmp;
+    }
+}
+
+
+/* ****************************************************************************
+ *  Function: write_string
+ *
+ *  Purpose:
+ *      Writes a character to the open file.
+ *
+ *  Programmer: Hank Childs
+ *  Creation:   September 3, 2004
+ * 
+ * ************************************************************************* */
+
+static void write_string(const char *str)
+{
+    fprintf(fp, str);
+}
+
+
+/* ****************************************************************************
+ *  Function: new_section
+ *
+ *  Purpose:
+ *      Adds a new line, provided we didn't already just do so and we are
+ *      writing an ASCII file.
+ *
+ *  Programmer: Hank Childs
+ *  Creation:   September 3, 2004
+ * 
+ * ************************************************************************* */
+
+static void new_section(void)
+{
+    if (numInColumn != 0)
+        end_line();
+    numInColumn = 0;
+}
+
+
+/* ****************************************************************************
+ *  Function: write_int
+ *
+ *  Purpose:
+ *      Writes an integer to the currently open file.  This routine takes
+ *      care of ASCII vs binary issues.
+ *
+ *  Programmer: Hank Childs
+ *  Creation:   September 3, 2004
+ * 
+ * ************************************************************************* */
+
+static void write_int(int val)
+{
+    if (useBinary)
+    {
+        force_big_endian((unsigned char *) &val);
+        fwrite(&val, sizeof(int), 1, fp);
+    }
+    else
+    {
+        char str[128];
+        sprintf(str, "%d ", val);
+        fprintf(fp, str);
+        if (((numInColumn++) % 9) == 8)
+        {
+            char str2[8] = "\n";
+            fprintf(fp, str2);
+            numInColumn = 0;
+        }
+    }
+}
+
+
+/* ****************************************************************************
+ *  Function: write_float
+ *
+ *  Purpose:
+ *      Writes an float to the currently open file.  This routine takes
+ *      care of ASCII vs binary issues.
+ *
+ *  Programmer: Hank Childs
+ *  Creation:   September 3, 2004
+ * 
+ *  Modifications:
+ *  
+ *    Hank Childs, Fri Apr 22 09:14:44 PDT 2005
+ *    Make precision changes suggested by Jeff McAninch
+ *
+ * ************************************************************************* */
+
+static void write_float(float val)
+{
+    if (useBinary)
+    {
+        force_big_endian((unsigned char *) &val);
+        fwrite(&val, sizeof(float), 1, fp);
+    }
+    else
+    {
+        char str[128];
+        sprintf(str, "%20.12e ", val);
+        fprintf(fp, str);
+        if (((numInColumn++) % 9) == 8)
+        {
+            end_line();
+        }
+    }
+}
+/* ****************************************************************************
+ *  Function: write_double
+ *
+ *  Purpose:
+ *      Writes a double to the currently open file.  This routine takes
+ *      care of ASCII vs binary issues.
+ *
+ *  Programmer: Hank Childs
+ *  Creation:   September 3, 2004
+ *
+ *  Modifications:
+ *
+ *    Hank Childs, Fri Apr 22 09:14:44 PDT 2005
+ *    Make precision changes suggested by Jeff McAninch
+ *
+ * ************************************************************************* */
+
+static void write_double(double val)
+{
+    if (useBinary)
+    {
+        force_double_big_endian((unsigned char *) &val);
+        fwrite(&val, sizeof(double), 1, fp);
+    }
+    else
+    {
+        char str[128];
+        sprintf(str, "%20.12e ", val);
+        fprintf(fp, str);
+        if (((numInColumn++) % 9) == 8)
+        {
+            end_line();
+        }
+    }
+}
+
+
+/* ****************************************************************************
+ *  Function: write_header
+ *
+ *  Purpose:
+ *      Writes the standard VTK header to the file.  This should be the first
+ *      thing written to the file.
+ *
+ *  Programmer: Hank Childs
+ *  Creation:   September 3, 2004
+ * 
+ * ************************************************************************* */
+
+static void write_header(void)
+{
+    fprintf(fp, "# vtk DataFile Version 2.0\n");
+    fprintf(fp, "Written using VisIt writer\n");
+    if (useBinary)
+        fprintf(fp, "BINARY\n");
+    else
+        fprintf(fp, "ASCII\n");
+}
+
+
+/* ****************************************************************************
+ *  Function: write_variables
+ *
+ *  Purpose:
+ *      Writes the variables to the file.  This can be a bit tricky.  The
+ *      cell data must be written first, followed by the point data.  When
+ *      writing the [point|cell] data, one variable must be declared the
+ *      primary scalar and another the primary vector (provided scalars
+ *      or vectors exist).  The rest of the arrays are added through the
+ *      "field data" mechanism.  Field data should support groups of arrays 
+ *      with different numbers of components (ie a scalar and a vector), but
+ *      there is a failure with the VTK reader.  So the scalars are all written
+ *      one group of field data and then the vectors as another.  If you don't
+ *      write it this way, the vectors do not show up.
+ *
+ *  Programmer: Hank Childs
+ *  Creation:   September 3, 2004
+ * 
+ * ************************************************************************* */
+
+void write_variables(int nvars, int *vardim, int *centering, 
+                     const char * const * varname, double **vars,
+                     int npts, int ncells)
+{
+    char str[1024];
+    int i, j, first_scalar, first_vector;
+    int num_scalars, num_vectors;
+    int num_field = 0;
+
+    new_section();
+    sprintf(str, "CELL_DATA %d\n", ncells);
+    write_string(str);
+
+    first_scalar = 0;
+    first_vector = 0;
+    num_scalars = 0;
+    num_vectors = 0;
+    /* The field data is where the non-primary scalars and vectors are 
+     * stored.  They must all be grouped together at the end of the point
+     * data.  So write out the primary scalars and vectors first.
+     */
+    for (i = 0 ; i < nvars ; i++)
+    {
+        if (centering[i] == 0)
+        {
+            int num_to_write = 0;
+            int should_write = 0;
+
+            if (vardim[i] == 1)
+            {
+                if (first_scalar == 0)
+                {
+                    should_write = 1;
+                    sprintf(str, "SCALARS %s double\n", varname[i]);
+                    write_string(str);
+                    write_string("LOOKUP_TABLE default\n");
+                    first_scalar = 1;
+                }
+                else 
+                    num_scalars++;
+            }
+            else if (vardim[i] == 3)
+            {
+                if (first_vector == 0)
+                {
+                    should_write = 1;
+                    sprintf(str, "VECTORS %s double\n", varname[i]);
+                    write_string(str);
+                    first_vector = 1;
+                }
+                else 
+                    num_vectors++;
+            }
+            else
+            {
+                printf("Only supported variable dimensions are 1 and 3.\n");
+                printf("Ignoring variable %s.\n", varname[i]);
+                continue;
+            }
+
+            if (should_write)
+            {
+                num_to_write = ncells*vardim[i];
+                for (j = 0 ; j < num_to_write ; j++)
+                {
+                    write_double(vars[i][j]);
+                }
+                end_line();
+            }
+        }
+    }
+
+    first_scalar = 0;
+    if (num_scalars > 0)
+    {
+        sprintf(str, "FIELD FieldData %d\n", num_scalars);
+        write_string(str);
+        for (i = 0 ; i < nvars ; i++)
+        {
+            int should_write = 0;
+            if (centering[i] == 0)
+            {
+                if (vardim[i] == 1)
+                {
+                    if (first_scalar == 0)
+                    {
+                        first_scalar = 1;
+                    }
+                    else
+                    {
+                        should_write = 1;
+                        sprintf(str, "%s 1 %d double\n", varname[i], ncells);
+                        write_string(str);
+                    }
+                }
+            }
+
+            if (should_write)
+            {
+                int num_to_write = ncells*vardim[i];
+                for (j = 0 ; j < num_to_write ; j++)
+                {
+                    write_double(vars[i][j]);
+                }
+                end_line();
+            }
+        }
+    }
+
+    first_vector = 0;
+    if (num_vectors > 0)
+    {
+        sprintf(str, "FIELD FieldData %d\n", num_vectors);
+        write_string(str);
+        for (i = 0 ; i < nvars ; i++)
+        {
+            int should_write = 0;
+            if (centering[i] == 0)
+            {
+                int num_to_write = 0;
+    
+                if (vardim[i] == 3)
+                {
+                    if (first_vector == 0)
+                    {
+                        first_vector = 1;
+                    }
+                    else
+                    {
+                        should_write = 1;
+                        sprintf(str, "%s 3 %d double\n", varname[i], ncells);
+                        write_string(str);
+                    }
+                }
+            }
+
+            if (should_write)
+            {
+                int num_to_write = ncells*vardim[i];
+                for (j = 0 ; j < num_to_write ; j++)
+                {
+                    write_double(vars[i][j]);
+                }
+                end_line();
+            }
+        }
+    }
+
+    new_section();
+    sprintf(str, "POINT_DATA %d\n", npts);
+    write_string(str);
+
+    first_scalar = 0;
+    first_vector = 0;
+    num_scalars = 0;
+    num_vectors = 0;
+    /* The field data is where the non-primary scalars and vectors are 
+     * stored.  They must all be grouped together at the end of the point
+     * data.  So write out the primary scalars and vectors first.
+     */
+    for (i = 0 ; i < nvars ; i++)
+    {
+        if (centering[i] != 0)
+        {
+            int num_to_write = 0;
+            int should_write = 0;
+
+            if (vardim[i] == 1)
+            {
+                if (first_scalar == 0)
+                {
+                    should_write = 1;
+                    sprintf(str, "SCALARS %s double\n", varname[i]);
+                    write_string(str);
+                    write_string("LOOKUP_TABLE default\n");
+                    first_scalar = 1;
+                }
+                else 
+                    num_scalars++;
+            }
+            else if (vardim[i] == 3)
+            {
+                if (first_vector == 0)
+                {
+                    should_write = 1;
+                    sprintf(str, "VECTORS %s double\n", varname[i]);
+                    write_string(str);
+                    first_vector = 1;
+                }
+                else 
+                    num_vectors++;
+            }
+            else
+            {
+                printf("Only supported variable dimensions are 1 and 3.\n");
+                printf("Ignoring variable %s.\n", varname[i]);
+                continue;
+            }
+
+            if (should_write)
+            {
+                num_to_write = npts*vardim[i];
+                for (j = 0 ; j < num_to_write ; j++)
+                {
+                    write_double(vars[i][j]);
+                }
+                end_line();
+            }
+        }
+    }
+
+    first_scalar = 0;
+    if (num_scalars > 0)
+    {
+        sprintf(str, "FIELD FieldData %d\n", num_scalars);
+        write_string(str);
+        for (i = 0 ; i < nvars ; i++)
+        {
+            int should_write = 0;
+            if (centering[i] != 0)
+            {
+                if (vardim[i] == 1)
+                {
+                    if (first_scalar == 0)
+                    {
+                        first_scalar = 1;
+                    }
+                    else
+                    {
+                        should_write = 1;
+                        sprintf(str, "%s 1 %d double\n", varname[i], npts);
+                        write_string(str);
+                    }
+                }
+            }
+
+            if (should_write)
+            {
+                int num_to_write = npts*vardim[i];
+                for (j = 0 ; j < num_to_write ; j++)
+                {
+                    write_double(vars[i][j]);
+                }
+                end_line();
+            }
+        }
+    }
+
+    first_vector = 0;
+    if (num_vectors > 0)
+    {
+        sprintf(str, "FIELD FieldData %d\n", num_vectors);
+        write_string(str);
+        for (i = 0 ; i < nvars ; i++)
+        {
+            int should_write = 0;
+            if (centering[i] != 0)
+            {
+                int num_to_write = 0;
+    
+                if (vardim[i] == 3)
+                {
+                    if (first_vector == 0)
+                    {
+                        first_vector = 1;
+                    }
+                    else
+                    {
+                        should_write = 1;
+                        sprintf(str, "%s 3 %d double\n", varname[i], npts);
+                        write_string(str);
+                    }
+                }
+            }
+
+            if (should_write)
+            {
+                int num_to_write = npts*vardim[i];
+                for (j = 0 ; j < num_to_write ; j++)
+                {
+                    write_double(vars[i][j]);
+                }
+                end_line();
+            }
+        }
+    }
+}
+
+
+/* ****************************************************************************
+//  Function: write_point_mesh
+//
+//  Purpose:
+//      Writes out a point mesh.
+//
+//  Arguments:
+//      filename   The name of the file to write.  If the extension ".vtk" is
+//                 not present, it will be added.
+//      useBinary  '0' to write ASCII, !0 to write binary
+//      npts       The number of points in the mesh.
+//      pts        The spatial locations of the points.  This array should
+//                 be size 3*npts.  The points should be encoded as:
+//                 <x1, y1, z1, x2, y2, z2, ..., xn, yn, zn>
+//      nvars      The number of variables.
+//      vardim     The dimension of each variable.  The size of vardim should
+//                 be nvars.  If var i is a scalar, then vardim[i] = 1.
+//                 If var i is a vector, then vardim[i] = 3.
+//      vars       An array of variables.  The size of vars should be nvars.
+//                 The size of vars[i] should be npts*vardim[i].
+//
+//  Programmer: Hank Childs
+//  Creation:   September 2, 2004
+//
+// ***************************************************************************/
+
+void write_point_mesh(const char *filename, int ub, int npts, float *pts,
+                      int nvars, int *vardim, const char * const *varnames,
+                      double **vars)
+{
+    int   i;
+    char  str[128];
+    int  *centering = NULL;
+
+    useBinary = ub;
+    open_file(filename);
+    write_header();
+
+    write_string("DATASET UNSTRUCTURED_GRID\n");
+    sprintf(str, "POINTS %d float\n", npts);
+    write_string(str);
+    for (i = 0 ; i < 3*npts ; i++)
+    {
+        write_float(pts[i]);
+    }
+
+    new_section();
+    sprintf(str, "CELLS %d %d\n", npts, 2*npts);
+    write_string(str);
+    for (i = 0 ; i < npts ; i++)
+    {
+        write_int(1);
+        write_int(i);
+        end_line();
+    }
+
+    new_section();
+    sprintf(str, "CELL_TYPES %d\n", npts);
+    write_string(str);
+    for (i = 0 ; i < npts ; i++)
+    {
+        write_int(VISIT_VERTEX);
+        end_line();
+    }
+
+    centering = (int *) malloc(nvars*sizeof(int));
+    for (i = 0 ; i < nvars ; i++)
+        centering[i] = 1;
+    write_variables(nvars, vardim, centering, varnames, vars, npts, npts);
+    free(centering);
+
+    close_file();
+}
+
+
+/* ****************************************************************************
+ *  Function: num_points_for_cell
+ *
+ *  Purpose:
+ *      Determines the number of points for the type of cell.
+ *
+ *  Programmer: Hank Childs
+ *  Creation:   September 3, 2004
+ * 
+ * ************************************************************************* */
+
+static int num_points_for_cell(int celltype)
+{
+    int npts = 0;
+    switch (celltype)
+    {
+       case VISIT_VERTEX:
+         npts = 1;
+         break;
+       case VISIT_LINE:
+         npts = 2;
+         break;
+       case VISIT_TRIANGLE:
+         npts = 3;
+         break;
+       case VISIT_QUAD:
+         npts = 4;
+         break;
+       case VISIT_TETRA:
+         npts = 4;
+         break;
+       case VISIT_HEXAHEDRON:
+         npts = 8;
+         break;
+       case VISIT_WEDGE:
+         npts = 6;
+         break;
+       case VISIT_PYRAMID:
+         npts = 5;
+         break;
+    }
+    return npts;
+}
+
+
+/* ****************************************************************************
+//  Function: write_unstructured_mesh
+//
+//  Purpose:
+//      Writes out a unstructured mesh.
+//
+//
+//  Arguments:
+//      filename   The name of the file to write.  If the extension ".vtk" is
+//                 not present, it will be added.
+//      useBinary  '0' to write ASCII, !0 to write binary
+//      npts       The number of points in the mesh.
+//      pts        The spatial locations of the points.  This array should
+//                 be size 3*npts.  The points should be encoded as:
+//                 <x1, y1, z1, x2, y2, z2, ..., xn, yn, zn>
+//      ncells     The number of cells.
+//      celltypes  The type of each cell.
+//      conn       The connectivity array.
+//      nvars      The number of variables.
+//      vardim     The dimension of each variable.  The size of vardim should
+//                 be nvars.  If var i is a scalar, then vardim[i] = 1.
+//                 If var i is a vector, then vardim[i] = 3.
+//      centering  The centering of each variable.  The size of centering
+//                 should be nvars.  If centering[i] == 0, then the variable
+//                 is cell-based.  If centering[i] != 0, then the variable
+//                 is point-based.
+//      vars       An array of variables.  The size of vars should be nvars.
+//                 The size of vars[i] should be npts*vardim[i].
+//
+//  Programmer: Hank Childs
+//  Creation:   September 2, 2004
+//
+// ***************************************************************************/
+
+void write_unstructured_mesh(const char *filename, int ub, int npts,
+                             float *pts, int ncells, int *celltypes, int *conn,
+                             int nvars, int *vardim, int *centering,
+                             const char * const *varnames, double **vars)
+{
+    int   i, j;
+    char  str[128];
+    int   conn_size = 0;
+    int  *curr_conn = conn;
+
+    useBinary = ub;
+    open_file(filename);
+    write_header();
+
+    write_string("DATASET UNSTRUCTURED_GRID\n");
+    sprintf(str, "POINTS %d float\n", npts);
+    write_string(str);
+    for (i = 0 ; i < 3*npts ; i++)
+    {
+        write_float(pts[i]);
+    }
+
+    new_section();
+    for (i = 0 ; i < ncells ; i++)
+    {
+        int npts = num_points_for_cell(celltypes[i]);
+         
+        conn_size += npts+1;
+    }
+    sprintf(str, "CELLS %d %d\n", ncells, conn_size);
+    write_string(str);
+    for (i = 0 ; i < ncells ; i++)
+    {
+        int npts = num_points_for_cell(celltypes[i]);
+        write_int(npts);
+        for (j = 0 ; j < npts ; j++)
+            write_int(*curr_conn++);
+        end_line();
+    }
+
+    new_section();
+    sprintf(str, "CELL_TYPES %d\n", ncells);
+    write_string(str);
+    for (i = 0 ; i < ncells ; i++)
+    {
+        write_int(celltypes[i]);
+        end_line();
+    }
+
+    write_variables(nvars, vardim, centering, varnames, vars, npts, ncells);
+
+    close_file();
+}
+
+
+/* ****************************************************************************
+//  Function: write_rectilinear_mesh
+//
+//  Purpose:
+//      Writes out a rectilinear mesh.
+//
+//
+//  Arguments:
+//      filename   The name of the file to write.  If the extension ".vtk" is
+//                 not present, it will be added.
+//      useBinary  '0' to write ASCII, !0 to write binary
+//      dims       An array of size 3 = { nX, nY, nZ }, where nX is the
+//                 number of points in the X-dimension, etc.
+//      x          An array of size dims[0] that contains the x-coordinates.
+//      y          An array of size dims[1] that contains the x-coordinates.
+//      z          An array of size dims[2] that contains the x-coordinates.
+//      nvars      The number of variables.
+//      vardim     The dimension of each variable.  The size of vardim should
+//                 be nvars.  If var i is a scalar, then vardim[i] = 1.
+//                 If var i is a vector, then vardim[i] = 3.
+//      centering  The centering of each variable.  The size of centering
+//                 should be nvars.  If centering[i] == 0, then the variable
+//                 is cell-based.  If centering[i] != 0, then the variable
+//                 is point-based.
+//      vars       An array of variables.  The size of vars should be nvars.
+//                 The size of vars[i] should be npts*vardim[i].
+//
+//
+//  Programmer: Hank Childs
+//  Creation:   September 2, 2004
+//
+//  Modifications:
+//
+//    Hank Childs, Wed Apr  6 16:22:57 PDT 2005
+//    Fix problem with 2D structured meshes and assessing cell count.
+//
+// ***************************************************************************/
+
+void write_rectilinear_mesh(const char *filename, int ub, int *dims,
+                            float *x, float *y, float *z,
+                            int nvars, int *vardim, int *centering,
+                            const char * const *varnames, double **vars)
+{
+    int   i, j;
+    char  str[128];
+    int npts = dims[0]*dims[1]*dims[2];
+    int ncX = (dims[0] - 1 < 1 ? 1 : dims[0] - 1);
+    int ncY = (dims[1] - 1 < 1 ? 1 : dims[1] - 1);
+    int ncZ = (dims[2] - 1 < 1 ? 1 : dims[2] - 1);
+    int ncells = ncX*ncY*ncZ;
+
+    useBinary = ub;
+    open_file(filename);
+    write_header();
+
+    write_string("DATASET RECTILINEAR_GRID\n");
+    sprintf(str, "DIMENSIONS %d %d %d\n", dims[0], dims[1], dims[2]);
+    write_string(str);
+    sprintf(str, "X_COORDINATES %d float\n", dims[0]);
+    write_string(str);
+    for (i = 0 ; i < dims[0] ; i++)
+        write_float(x[i]);
+    new_section();
+    sprintf(str, "Y_COORDINATES %d float\n", dims[1]);
+    write_string(str);
+    for (i = 0 ; i < dims[1] ; i++)
+        write_float(y[i]);
+    new_section();
+    sprintf(str, "Z_COORDINATES %d float\n", dims[2]);
+    write_string(str);
+    for (i = 0 ; i < dims[2] ; i++)
+        write_float(z[i]);
+
+    write_variables(nvars, vardim, centering, varnames, vars, npts, ncells);
+
+    close_file();
+}
+
+
+/* ****************************************************************************
+//  Function: write_regular_mesh
+//
+//  Purpose:
+//      Writes out a regular mesh.  A regular mesh is one where the data lies
+//      along regular intervals.  "Brick of bytes/doubles",
+//      "Block of bytes/double", and MRI data all are examples of data that
+//      lie on regular meshes.
+//
+//
+//  Arguments:
+//      filename   The name of the file to write.  If the extension ".vtk" is
+//                 not present, it will be added.
+//      useBinary  '0' to write ASCII, !0 to write binary
+//      dims       An array of size 3 = { nX, nY, nZ }, where nX is the
+//                 number of points in the X-dimension, etc.
+//      nvars      The number of variables.
+//      vardim     The dimension of each variable.  The size of vardim should
+//                 be nvars.  If var i is a scalar, then vardim[i] = 1.
+//                 If var i is a vector, then vardim[i] = 3.
+//      centering  The centering of each variable.  The size of centering
+//                 should be nvars.  If centering[i] == 0, then the variable
+//                 is cell-based.  If centering[i] != 0, then the variable
+//                 is point-based.
+//      vars       An array of variables.  The size of vars should be nvars.
+//                 The size of vars[i] should be npts*vardim[i].
+//
+//
+//  Programmer: Hank Childs
+//  Creation:   September 2, 2004
+//
+// ***************************************************************************/
+
+void write_regular_mesh(const char *filename, int ub, int *dims,
+                        int nvars, int *vardim, int *centering,
+                        const char * const *varnames, double **vars)
+{
+    int  i;
+
+    float *x = (float *) malloc(sizeof(float)*dims[0]);
+    float *y = (float *) malloc(sizeof(float)*dims[1]);
+    float *z = (float *) malloc(sizeof(float)*dims[2]);
+
+    for (i = 0 ; i < dims[0] ; i++)
+        x[i] = (float) i;
+    for (i = 0 ; i < dims[1] ; i++)
+        y[i] = (float) i;
+    for (i = 0 ; i < dims[2] ; i++)
+        z[i] = (float) i;
+
+    write_rectilinear_mesh(filename, ub, dims, x, y, z, nvars, vardim,
+                           centering, varnames, vars);
+
+    free(x);
+    free(y);
+    free(z);
+}
+
+
+/* ****************************************************************************
+//  Function: write_curvilinear_mesh
+//
+//  Purpose:
+//      Writes out a curvilinear mesh.
+//
+//
+//  Arguments:
+//      filename   The name of the file to write.  If the extension ".vtk" is
+//                 not present, it will be added.
+//      useBinary  '0' to write ASCII, !0 to write binary
+//      dims       An array of size 3 = { nI, nJ, nK }, where nI is the
+//                 number of points in the logical I dimension, etc.
+//      pts        An array of size nI*nJ*nK*3.  The array should be layed
+//                 out as (pt(i=0,j=0,k=0), pt(i=1,j=0,k=0), ...
+//                 pt(i=nI-1,j=0,k=0), pt(i=0,j=1,k=0), ...).
+//      nvars      The number of variables.
+//      vardim     The dimension of each variable.  The size of vardim should
+//                 be nvars.  If var i is a scalar, then vardim[i] = 1.
+//                 If var i is a vector, then vardim[i] = 3.
+//      centering  The centering of each variable.  The size of centering
+//                 should be nvars.  If centering[i] == 0, then the variable
+//                 is cell-based.  If centering[i] != 0, then the variable
+//                 is point-based.
+//      vars       An array of variables.  The size of vars should be nvars.
+//                 The size of vars[i] should be npts*vardim[i].
+//
+//
+//  Programmer: Hank Childs
+//  Creation:   September 2, 2004
+//
+//  Modifications:
+//
+//    Hank Childs, Wed Apr  6 16:22:57 PDT 2005
+//    Fix problem with 2D structured meshes and assessing cell count.
+//
+// ***************************************************************************/
+
+void write_curvilinear_mesh(const char *filename, int ub, int *dims,float *pts,
+                            int nvars, int *vardim, int *centering,
+                            const char * const *varnames, double **vars)
+{
+    int   i, j;
+    char  str[128];
+    int npts = dims[0]*dims[1]*dims[2];
+    int ncX = (dims[0] - 1 < 1 ? 1 : dims[0] - 1);
+    int ncY = (dims[1] - 1 < 1 ? 1 : dims[1] - 1);
+    int ncZ = (dims[2] - 1 < 1 ? 1 : dims[2] - 1);
+    int ncells = ncX*ncY*ncZ;
+
+    useBinary = ub;
+    open_file(filename);
+    write_header();
+
+    write_string("DATASET STRUCTURED_GRID\n");
+    sprintf(str, "DIMENSIONS %d %d %d\n", dims[0], dims[1], dims[2]);
+    write_string(str);
+    sprintf(str, "POINTS %d float\n", npts);
+    write_string(str);
+    for (i = 0 ; i < 3*npts ; i++)
+    {
+        write_float(pts[i]);
+    }
+
+    write_variables(nvars, vardim, centering, varnames, vars, npts, ncells);
+
+    close_file();
+}
+
+
diff --git a/examples/incompressible-navier-stokes/visit_writer.h b/examples/incompressible-navier-stokes/visit_writer.h
new file mode 100644
index 0000000000000000000000000000000000000000..ae67e3ccc8de624ccf7f027cddf223cc32cd476c
--- /dev/null
+++ b/examples/incompressible-navier-stokes/visit_writer.h
@@ -0,0 +1,271 @@
+/* ************************************************************************* //
+//                              visit_writer.h                               //
+// ************************************************************************* */
+
+/* 
+// This file contains function prototypes for writing out point meshes, 
+// unstructured meshes, rectilinear meshes, regular meshes, and 
+// structured/curvilinear meshes into files that can later be read by VisIt.
+//
+// Each routine assumes that the data being written is three-dimensional.
+// If the data is two-dimensional, you must still write out the data
+// as three-dimensional (ie pad arrays so that they are the correct size, etc).
+// However: the VisIt reader will determine that the data is truly two-
+// dimensional and visualize it as a two-dimensional dataset.
+//
+// All writers have an ASCII vs Binary decision.  The tradeoffs are the
+// standard ones: ASCII is human readable, but slow.  The
+// binary is much faster, but not human readable.  Note: the binary format
+// is portable, since it converts all data to be big-endian (this was a 
+// design decision for the format the visit_writer writes to -- the VTK
+// format).
+//
+// If you have multiple grids, you can write out one file for each grid.
+// There are potential pitfalls in doing this, where extra geometry and
+// interpolation problems appear along grid boundaries.  For additional
+// help with this issue, e-mail visit-help@llnl.gov
+*/
+
+
+/* ****************************************************************************
+//  Function: write_point_mesh
+//
+//  Purpose:
+//      Writes out a point mesh.
+//
+//  Arguments:
+//      filename   The name of the file to write.  If the extension ".vtk" is
+//                 not present, it will be added.
+//      useBinary  '0' to write ASCII, !0 to write binary
+//      npts       The number of points in the mesh.
+//      pts        The spatial locations of the points.  This array should
+//                 be size 3*npts.  The points should be encoded as:
+//                 <x1, y1, z1, x2, y2, z2, ..., xn, yn, zn>
+//      nvars      The number of variables.
+//      vardim     The dimension of each variable.  The size of vardim should
+//                 be nvars.  If var i is a scalar, then vardim[i] = 1.
+//                 If var i is a vector, then vardim[i] = 3.
+//      vars       An array of variables.  The size of vars should be nvars.
+//                 The size of vars[i] should be npts*vardim[i].
+//      
+//  Programmer: Hank Childs
+//  Creation:   September 2, 2004
+//
+// ***************************************************************************/
+
+void write_point_mesh(const char *filename, int useBinary, int npts, 
+                      float *pts, int nvars, int *vardim, 
+                      const char * const *varnames, double **vars);
+
+
+
+/* ****************************************************************************
+//  Function: write_unstructured_mesh
+//
+//  Purpose:
+//      Writes out a unstructured mesh.
+//     
+//
+//  Arguments:
+//      filename   The name of the file to write.  If the extension ".vtk" is
+//                 not present, it will be added.
+//      useBinary  '0' to write ASCII, !0 to write binary
+//      npts       The number of points in the mesh.
+//      pts        The spatial locations of the points.  This array should
+//                 be size 3*npts.  The points should be encoded as:
+//                 <x1, y1, z1, x2, y2, z2, ..., xn, yn, zn>
+//      ncells     The number of cells.
+//      celltypes  The type of each cell.
+//      conn       The connectivity array.
+//      nvars      The number of variables.
+//      vardim     The dimension of each variable.  The size of vardim should
+//                 be nvars.  If var i is a scalar, then vardim[i] = 1.
+//                 If var i is a vector, then vardim[i] = 3.
+//      centering  The centering of each variable.  The size of centering 
+//                 should be nvars.  If centering[i] == 0, then the variable
+//                 is cell-based.  If centering[i] != 0, then the variable
+//                 is point-based.
+//      vars       An array of variables.  The size of vars should be nvars.
+//                 The size of vars[i] should be npts*vardim[i].
+//      
+//  Example:
+//      You have two triangles.  The first has points (0,0,0), (0,1,0), and
+//      (1,1,0).  The second has points (0,0,0), (1,1,0), and (1,0,0).
+//
+//      There are four unique points. 
+//
+//      float pts[12] = { 0,0,0, 0,1,0, 1,1,0, 1,0,0 };
+//
+//      It is important the points list contain only unique points,
+//      because VisIt is not able to correctly determine the connectivity of a 
+//      dataset when points are duplicated.
+//
+//      There are two triangles.
+//      int ncells = 2;
+//
+//      The cells are both triangles.
+//      int celltypes[2] = { VISIT_TRIANGLE, VISIT_TRIANGLE };
+//
+//      The connectivity contains indices into the points list.  The indexing
+//      assumes that each point has size 3 (x,y,z).
+//
+//      int conn[6] = { 0, 1, 2, 0, 2, 3 };
+//
+//  Hint:  
+//      When writing an unstructured mesh, it is easy to get the orientation
+//      of a cell backwards.  VisIt typically does okay with this, but it
+//      can cause problems.  To test if this is happening, bring up VisIt on
+//      your newly outputted dataset and make a Pseudocolor plot of 
+//      "mesh_quality/volume" for 3D datasets or "mesh_quality/area" for 2D 
+//      datasets.  If the cells are inside-out, the volumes or areas will be
+//      negative.
+//      
+//
+//  Programmer: Hank Childs
+//  Creation:   September 2, 2004
+//
+// ***************************************************************************/
+
+#define VISIT_VERTEX         1
+#define VISIT_LINE           3
+#define VISIT_TRIANGLE       5
+#define VISIT_QUAD           9
+#define VISIT_TETRA         10
+#define VISIT_HEXAHEDRON    12
+#define VISIT_WEDGE         13
+#define VISIT_PYRAMID       14
+
+void write_unstructured_mesh(const char *filename, int useBinary, int npts,
+                             float *pts, int ncells, int *celltypes, int *conn, 
+                             int nvars, int *vardim, int *centering,
+                             const char * const *varnames, double **vars);
+
+
+
+/* ****************************************************************************
+//  Function: write_regular_mesh
+//
+//  Purpose:
+//      Writes out a regular mesh.  A regular mesh is one where the data lies
+//      along regular intervals.  "Brick of bytes/floats", 
+//      "Block of bytes/floats", and MRI data all are examples of data that
+//      lie on regular meshes.
+//     
+//
+//  Arguments:
+//      filename   The name of the file to write.  If the extension ".vtk" is
+//                 not present, it will be added.
+//      useBinary  '0' to write ASCII, !0 to write binary
+//      dims       An array of size 3 = { nX, nY, nZ }, where nX is the
+//                 number of points in the X-dimension, etc.
+//      nvars      The number of variables.
+//      vardim     The dimension of each variable.  The size of vardim should
+//                 be nvars.  If var i is a scalar, then vardim[i] = 1.
+//                 If var i is a vector, then vardim[i] = 3.
+//      centering  The centering of each variable.  The size of centering 
+//                 should be nvars.  If centering[i] == 0, then the variable
+//                 is cell-based.  If centering[i] != 0, then the variable
+//                 is point-based.
+//      vars       An array of variables.  The size of vars should be nvars.
+//                 The size of vars[i] should be npts*vardim[i].
+//      
+//
+//  Programmer: Hank Childs
+//  Creation:   September 2, 2004
+//
+// ***************************************************************************/
+
+void write_regular_mesh(const char *filename, int useBinary, int *dims, 
+                        int nvars, int *vardim, int *centering,
+                        const char * const *varnames, double **vars);
+
+
+
+
+/* ****************************************************************************
+//  Function: write_rectilinear_mesh
+//
+//  Purpose:
+//      Writes out a rectilinear mesh.
+//     
+//
+//  Arguments:
+//      filename   The name of the file to write.  If the extension ".vtk" is
+//                 not present, it will be added.
+//      useBinary  '0' to write ASCII, !0 to write binary
+//      dims       An array of size 3 = { nX, nY, nZ }, where nX is the
+//                 number of points in the X-dimension, etc.
+//      x          An array of size dims[0] that contains the x-coordinates.
+//      y          An array of size dims[1] that contains the x-coordinates.
+//      z          An array of size dims[2] that contains the x-coordinates.
+//      nvars      The number of variables.
+//      vardim     The dimension of each variable.  The size of vardim should
+//                 be nvars.  If var i is a scalar, then vardim[i] = 1.
+//                 If var i is a vector, then vardim[i] = 3.
+//      centering  The centering of each variable.  The size of centering 
+//                 should be nvars.  If centering[i] == 0, then the variable
+//                 is cell-based.  If centering[i] != 0, then the variable
+//                 is point-based.
+//      vars       An array of variables.  The size of vars should be nvars.
+//                 The size of vars[i] should be npts*vardim[i].
+//      
+//
+//  Example:
+//      You have a rectilinear mesh with x = { 0, 1, 2}, y = { 1, 1.5, 2, 3 },
+//      and z = { 2.5, 3.5 }.
+//
+//      Then dims = { 3, 4, 2 }.
+//      
+//  Programmer: Hank Childs
+//  Creation:   September 2, 2004
+//
+// ***************************************************************************/
+
+void write_rectilinear_mesh(const char *filename, int useBinary, 
+                            int *dims, float *x, float *y, float *z, 
+                            int nvars, int *vardim, int *centering, 
+                            const char * const *varnames, double **vars);
+
+
+
+
+/* ****************************************************************************
+//  Function: write_curvilinear_mesh
+//
+//  Purpose:
+//      Writes out a curvilinear mesh.
+//     
+//
+//  Arguments:
+//      filename   The name of the file to write.  If the extension ".vtk" is
+//                 not present, it will be added.
+//      useBinary  '0' to write ASCII, !0 to write binary
+//      dims       An array of size 3 = { nI, nJ, nK }, where nI is the
+//                 number of points in the logical I dimension, etc.
+//      pts        An array of size nI*nJ*nK*3.  The array should be layed
+//                 out as (pt(i=0,j=0,k=0), pt(i=1,j=0,k=0), ...
+//                 pt(i=nI-1,j=0,k=0), pt(i=0,j=1,k=0), ...).
+//      nvars      The number of variables.
+//      vardim     The dimension of each variable.  The size of vardim should
+//                 be nvars.  If var i is a scalar, then vardim[i] = 1.
+//                 If var i is a vector, then vardim[i] = 3.
+//      centering  The centering of each variable.  The size of centering 
+//                 should be nvars.  If centering[i] == 0, then the variable
+//                 is cell-based.  If centering[i] != 0, then the variable
+//                 is point-based.
+//      vars       An array of variables.  The size of vars should be nvars.
+//                 The size of vars[i] should be npts*vardim[i].
+//      
+//
+//  Programmer: Hank Childs
+//  Creation:   September 2, 2004
+//
+// ***************************************************************************/
+
+void write_curvilinear_mesh(const char *filename, int useBinary, 
+                            int *dims, float *pts,
+                            int nvars, int *vardim, int *centering, 
+                            const char * const *varnames, double **vars);
+
+
+