WriterVTK.cpp

Go to the documentation of this file.

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 Stéphane Del Pino

//  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, or (at your option)
//  any later version.

//  This program is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

//  $Id: WriterVTK.cpp,v 1.14 2008/05/25 11:56:36 delpinux Exp $

#include <WriterVTK.hpp>
#include <Mesh.hpp>

#include <FieldOfScalarFunction.hpp>

#include <ScalarFunctionBase.hpp>

#include <FEMFunctionBase.hpp>

#include <Mesh.hpp>

#include <Structured3DMesh.hpp>

#include <MeshOfTetrahedra.hpp>
#include <MeshOfHexahedra.hpp>

#include <SpectralMesh.hpp>

#include <OctreeMesh.hpp>

#include <SurfaceMeshOfTriangles.hpp>
#include <SurfaceMeshOfQuadrangles.hpp>

#include <XMLWriter.hpp>

#include <Vector.hpp>

#include <EndianConverter.hpp>
#include <fstream>

class BinarySerializerVTK
{
  std::vector<char> __data;

  template <typename T>
  struct DataTypeTraits {
    typedef T DataType;
  };

  template <typename T>
  void __add(T data)
  {
    littleEndianize(data);
    for (size_t i=0; i<sizeof(T); ++i) {
      __data.push_back(reinterpret_cast<char*>(&data)[i]);
    }
  }

public:
  size_t offset() const
  {
    return __data.size();
  }

  friend std::ostream& operator<<(std::ostream& os,
                                  const BinarySerializerVTK& serializer)
  {
    os.write(&serializer.__data[0], serializer.__data.size());
    return os;
  }

  template <typename T>
  void add(const Vector<T>& values)
  {
    typedef typename DataTypeTraits<T>::DataType DataType;

    int size = values.size()*sizeof(DataType);
    this->__add(size);
    for (size_t i=0; i<values.size(); ++i) {
      DataType v = values[i];
      this->__add(v);
    }
  }

  BinarySerializerVTK()
  {
    ;
  }
};

template <>
struct BinarySerializerVTK::DataTypeTraits<real_t>
{
  typedef double DataType;
};

BinarySerializerVTK serializer;



template <typename CellType>
struct WriterVTK::Traits {};

template <>
struct WriterVTK::Traits<Triangle>
{
  enum {
    VTKType = 5
  };
};

template <>
struct WriterVTK::Traits<Quadrangle>
{
  enum {
    VTKType = 9
  };
};

template <>
struct WriterVTK::Traits<Tetrahedron>
{
  enum {
    VTKType = 10
  };
};

template <>
struct WriterVTK::Traits<Hexahedron>
{
  enum {
    VTKType = 12
  };
};

template <>
struct WriterVTK::Traits<CartesianHexahedron>
{
  enum {
    VTKType = 12
  };
};

void WriterVTK::
__fillCrossedComponent(const FieldOfScalarFunction& field,
                       Vector<real_t>& values) const
{
  const size_t numberOfComponents = field.numberOfComponents();
  ASSERT(values.size() == __mesh->numberOfVertices()*numberOfComponents);

  switch (__mesh->type()) {
  case Mesh::cartesianHexahedraMesh: {
    const Structured3DMesh& mesh = dynamic_cast<const Structured3DMesh&>(*__mesh);
    for (size_t l = 0; l<numberOfComponents; ++l) {
      const ScalarFunctionBase& function = *field.function(l);

      switch (function.type()) {
      case ScalarFunctionBase::femfunction: {
        const FEMFunctionBase& fem
          = static_cast<const FEMFunctionBase&>(function);

        if ((fem.discretizationType() == ScalarDiscretizationTypeBase::lagrangianFEM1)
            and (fem.baseMesh() == __mesh)) {
          size_t position = 0;
          for (size_t k=0; k<mesh.shape().nz(); ++k)
            for (size_t j=0; j<mesh.shape().ny(); ++j)
              for (size_t i=0; i<mesh.shape().nx(); ++i) {
                const size_t pointNumber = mesh.shape()(i,j,k);
                values[position*numberOfComponents+l] = fem[pointNumber];
                position++;
              }
          break;
        }
      }
      default: {
        size_t position = 0;
        for (size_t k=0; k<mesh.shape().nz(); ++k)
          for (size_t j=0; j<mesh.shape().ny(); ++j)
            for (size_t i=0; i<mesh.shape().nx(); ++i) {
              const Vertex& x = mesh.vertex(i,j,k);
              values[position*numberOfComponents+l] = function(x);
              position++;
            }
      }
      }
    }
    break;
  }
  default: {
    for (size_t i = 0; i<numberOfComponents; ++i) {
      const ScalarFunctionBase& function = *field.function(i);

      switch (function.type()) {
      case ScalarFunctionBase::femfunction: {
        const FEMFunctionBase& fem
          = static_cast<const FEMFunctionBase&>(function);
        if ((fem.baseMesh() == __mesh) and
            ((fem.discretizationType() == ScalarDiscretizationTypeBase::lagrangianFEM1))) {
          for (size_t j=0; j<fem.values().size(); ++j) {
            values[numberOfComponents*j+i] = fem[j];
          }
          break;
        } // if not continues the standard method
      }
      default: {
        for (size_t j=0; j<__mesh->numberOfVertices(); ++j) {
          const TinyVector<3,real_t>& x = __mesh->vertex(j);
          values[numberOfComponents*j+i] = function(x);
        }
      }
      }
    }
  }
  }
}

template <typename MeshType>
void WriterVTK::
__proceed() const
{
  typedef typename MeshType::CellType CellType;

  std::string filename = __filename;
  filename += ".vtu";

  const MeshType& mesh
    = static_cast<const MeshType&>(*__mesh);

//   size_t surfaceMeshNumberOfCells = 0;
//   if (mesh.hasaceMesh()) {
//     surfaceMeshNumberOfCells = mesh.surfaceMesh()->numberOfCells();
//   }

  std::ofstream file(filename.c_str());
  file.precision(15);

  XMLWriter xmlWriter(file);

  xmlWriter.writeHeader();
  {
    XMLTag vtkFile("VTKFile");
    ReferenceCounting<XMLAttribute> attribute
      = new XMLAttribute("type","UnstructuredGrid");
    vtkFile.add(attribute);

    if (this->__binary) {
      ReferenceCounting<XMLAttribute> endiannes
        = new XMLAttribute("byte_order","LittleEndian");
      vtkFile.add(endiannes);
    }

    xmlWriter.add(vtkFile);
  }

  {
    XMLTag unstructuredGrid("UnstructuredGrid");
    xmlWriter.add(unstructuredGrid);
  }

  {
    XMLTag piece("Piece");
    ReferenceCounting<XMLAttribute>
      nbPoints = new XMLAttribute("NumberOfPoints",
                                  stringify(mesh.numberOfVertices()));
    piece.add(nbPoints);
    ReferenceCounting<XMLAttribute>
      nbCells = new XMLAttribute("NumberOfCells",
                                 stringify(mesh.numberOfCells()));
    piece.add(nbCells);
    xmlWriter.add(piece);

    {
      XMLTag points("Points");
      xmlWriter.add(points);
   
      {  
        XMLTag dataarray("DataArray");
        ReferenceCounting<XMLAttribute>
          type = new XMLAttribute("type","Float64");
        dataarray.add(type);
        ReferenceCounting<XMLAttribute>
          name = new XMLAttribute("Name","Position");
        dataarray.add(name);
        ReferenceCounting<XMLAttribute>
          nbComponent = new XMLAttribute("NumberOfComponents","3");
        dataarray.add(nbComponent);

        if (not this->__binary) {
          ReferenceCounting<XMLAttribute>
            format = new XMLAttribute("format","ascii");
          dataarray.add(format);

          xmlWriter.add(dataarray);

          for (size_t i=0; i<mesh.numberOfVertices(); ++i) {
            const Vertex& x = mesh.vertex(i); 
            xmlWriter.insert(x[0]);
            xmlWriter.insert(x[1]);
            xmlWriter.insert(x[2]);
          }
          xmlWriter.insertNewLine();
          xmlWriter.closeTag();
        } else {
          ReferenceCounting<XMLAttribute>
            format = new XMLAttribute("format","appended");
          dataarray.add(format);

          ReferenceCounting<XMLAttribute>
            offset = new XMLAttribute("offset",stringify(serializer.offset()));
          dataarray.add(offset);

          xmlWriter.add(dataarray);

          Vector<real_t> coords(3*mesh.numberOfVertices());
          for (size_t i=0; i<mesh.numberOfVertices(); ++i) {
            const Vertex& x = mesh.vertex(i); 
            coords[3*i+0] = x[0];
            coords[3*i+1] = x[1];
            coords[3*i+2] = x[2];
          }

          serializer.add(coords);
          xmlWriter.closeTag();
        }
      }
      xmlWriter.closeTag();
    }

    {
      XMLTag cells("Cells");
      xmlWriter.add(cells);

      {
        XMLTag dataarray("DataArray");
        ReferenceCounting<XMLAttribute>
          type = new XMLAttribute("type","Int32");
        dataarray.add(type);
        ReferenceCounting<XMLAttribute>
          name = new XMLAttribute("Name","connectivity");
        dataarray.add(name);
        ReferenceCounting<XMLAttribute>
          nbComponent = new XMLAttribute("NumberOfComponents","1");
        dataarray.add(nbComponent);

        if (not this->__binary) {
          ReferenceCounting<XMLAttribute>
            format = new XMLAttribute("format","ascii");
          dataarray.add(format);
          xmlWriter.add(dataarray);

          for (size_t i=0;i<mesh.numberOfCells(); ++i) {
            const CellType& cell = mesh.cell(i);
            for (size_t j=0; j<CellType::NumberOfVertices; ++j) {
              xmlWriter.insert(mesh.vertexNumber(cell(j)));
            }
          }
          xmlWriter.insertNewLine();
          xmlWriter.closeTag();
        } else {
          ReferenceCounting<XMLAttribute>
            format = new XMLAttribute("format","appended");
          dataarray.add(format);

          ReferenceCounting<XMLAttribute>
            offset = new XMLAttribute("offset",stringify(serializer.offset()));
          dataarray.add(offset);

          Vector<int> cellDescription(CellType::NumberOfVertices*mesh.numberOfCells());

          for (size_t i=0;i<mesh.numberOfCells(); ++i) {
            const CellType& cell = mesh.cell(i);
            for (size_t j=0; j<CellType::NumberOfVertices; ++j) {
              cellDescription[i*CellType::NumberOfVertices+j]
                = mesh.vertexNumber(cell(j));
            }
          }

          serializer.add(cellDescription);

          xmlWriter.add(dataarray);
          xmlWriter.closeTag();
        }
      }

      {
        XMLTag dataarray("DataArray");
        ReferenceCounting<XMLAttribute>
          type = new XMLAttribute("type","Int32");
        dataarray.add(type);
        ReferenceCounting<XMLAttribute>
          name = new XMLAttribute("Name","offsets");
        dataarray.add(name);
        ReferenceCounting<XMLAttribute>
          nbComponent = new XMLAttribute("NumberOfComponents","1");
        dataarray.add(nbComponent);
        if (not this->__binary) {
          ReferenceCounting<XMLAttribute>
            format = new XMLAttribute("format","ascii");
          dataarray.add(format);
          xmlWriter.add(dataarray);
     
          for (size_t i=1; i<=mesh.numberOfCells(); ++i) {
            xmlWriter.insert(i*CellType::NumberOfVertices);
          }

          xmlWriter.insertNewLine();
          xmlWriter.closeTag();      
        } else {
          ReferenceCounting<XMLAttribute>
            format = new XMLAttribute("format","appended");
          dataarray.add(format);

          ReferenceCounting<XMLAttribute>
            offset = new XMLAttribute("offset",stringify(serializer.offset()));
          dataarray.add(offset);

          Vector<int> offsets(mesh.numberOfCells());
          for (size_t i=0; i<mesh.numberOfCells(); ++i) {
            offsets[i] = ((i+1)*CellType::NumberOfVertices);
          }

          serializer.add(offsets);

          xmlWriter.add(dataarray);
          xmlWriter.closeTag();
        }
      }

      {
        XMLTag dataarray("DataArray");
        ReferenceCounting<XMLAttribute>
          type = new XMLAttribute("type","UInt8");
        dataarray.add(type);
        ReferenceCounting<XMLAttribute>
          name = new XMLAttribute("Name","types");
        dataarray.add(name);
        ReferenceCounting<XMLAttribute>
          nbComponent = new XMLAttribute("NumberOfComponents","1");
        dataarray.add(nbComponent);

        if (not this->__binary) {
          ReferenceCounting<XMLAttribute>
            format = new XMLAttribute("format","ascii");
          dataarray.add(format);
          xmlWriter.add(dataarray);
     
          for (size_t i=0; i<mesh.numberOfCells(); ++i) {
            xmlWriter.insert(Traits<CellType>::VTKType);
          }
          xmlWriter.insertNewLine();
          xmlWriter.closeTag();
        } else {
          ReferenceCounting<XMLAttribute>
            format = new XMLAttribute("format","appended");
          dataarray.add(format);

          ReferenceCounting<XMLAttribute>
            offset = new XMLAttribute("offset",stringify(serializer.offset()));
          dataarray.add(offset);

          Vector<unsigned char> types(mesh.numberOfCells());
          for (size_t i=0; i<mesh.numberOfCells(); ++i) {
            types[i] = Traits<CellType>::VTKType;
          }

          serializer.add(types);

          xmlWriter.add(dataarray);
          xmlWriter.closeTag();
        }
      }
    }
    xmlWriter.closeTag();
  }
  {
    XMLTag pointdata("PointData");
    xmlWriter.add(pointdata);

    for (size_t i=0; i<__scalarFunctionList.size(); ++i) {

      const ScalarFunctionBase& function = *__scalarFunctionList[i];

      if (not this->__binary) {
        XMLTag dataarray("DataArray");
        ReferenceCounting<XMLAttribute>
          type = new XMLAttribute("type","Float64");
        dataarray.add(type);
        ReferenceCounting<XMLAttribute>
          name = new XMLAttribute("Name",stringify(function));
        dataarray.add(name);
        ReferenceCounting<XMLAttribute>
          nbComponent = new XMLAttribute("NumberOfComponents","1");
        dataarray.add(nbComponent);
        ReferenceCounting<XMLAttribute>
          format = new XMLAttribute("format","ascii");
        dataarray.add(format);
    
        xmlWriter.add(dataarray);

        switch (function.type()) {
        case ScalarFunctionBase::femfunction: {
          const FEMFunctionBase& fem
            = static_cast<const FEMFunctionBase&>(function);

          if ((fem.discretizationType() == ScalarDiscretizationTypeBase::lagrangianFEM1)
              and (fem.baseMesh() == __mesh)) {
            for (size_t i=0; i<mesh.numberOfVertices(); ++i) {
              xmlWriter.insert(fem[i]);
            }
            break;
          }
        }
        default: {
          for (size_t i=0; i<mesh.numberOfVertices(); ++i) {
            const Vertex& x = mesh.vertex(i);
            xmlWriter.insert(function(x));
          }
        }
        }

        xmlWriter.insertNewLine();
        xmlWriter.closeTag();
      } else {
        XMLTag dataArray("DataArray");
        ReferenceCounting<XMLAttribute> type
          = new XMLAttribute("type","Float64");
        dataArray.add(type);
        ReferenceCounting<XMLAttribute> name
          = new XMLAttribute("Name",stringify(function));
        dataArray.add(name);
        ReferenceCounting<XMLAttribute>
          nbComponent = new XMLAttribute("NumberOfComponents","1");
        dataArray.add(nbComponent);

        ReferenceCounting<XMLAttribute>
          format = new XMLAttribute("format","appended");
        dataArray.add(format);

        ReferenceCounting<XMLAttribute>
          offset = new XMLAttribute("offset",stringify(serializer.offset()));
        dataArray.add(offset);

        xmlWriter.add(dataArray);
        xmlWriter.closeTag();

        Vector<real_t> values(mesh.numberOfVertices());

        switch (function.type()) {
        case ScalarFunctionBase::femfunction: {
          const FEMFunctionBase& fem
            = static_cast<const FEMFunctionBase&>(function);

          if ((fem.discretizationType() == ScalarDiscretizationTypeBase::lagrangianFEM1)
              and (fem.baseMesh() == __mesh)) {
            for (size_t i=0; i<mesh.numberOfVertices(); ++i) {
              values[i] = fem[i];
            }
            break;
          }
        }
        default: {
          for (size_t i=0; i<mesh.numberOfVertices(); ++i) {
            const Vertex& x = mesh.vertex(i);
            values[i] = function(x);
          }
        }
        }
        serializer.add(values);
      }
    }

    for (size_t i=0; i<__fieldList.size(); ++i) {

      const FieldOfScalarFunction& field = *__fieldList[i];

      XMLTag dataArray("DataArray");
      ReferenceCounting<XMLAttribute> type
        = new XMLAttribute("type","Float64");
      dataArray.add(type);
      ReferenceCounting<XMLAttribute> name
        = new XMLAttribute("Name",stringify(field));
      dataArray.add(name);
      ReferenceCounting<XMLAttribute>
        nbComponent = new XMLAttribute("NumberOfComponents",
                                       stringify(field.numberOfComponents()));
      dataArray.add(nbComponent);

      Vector<real_t> values(field.numberOfComponents()*__mesh->numberOfVertices());

      this->__fillCrossedComponent(field,values);
      if (not this->__binary) {
        ReferenceCounting<XMLAttribute> format
          = new XMLAttribute("format","ascii");
        dataArray.add(format);
        xmlWriter.add(dataArray);


        for (size_t i=0; i<values.size(); ++i) {
          xmlWriter.insert(values[i]);
        }

        xmlWriter.insertNewLine();
        xmlWriter.closeTag();
      } else {
        ReferenceCounting<XMLAttribute>
          format = new XMLAttribute("format","appended");
        dataArray.add(format);

        ReferenceCounting<XMLAttribute>
          offset = new XMLAttribute("offset",stringify(serializer.offset()));
        dataArray.add(offset);

        xmlWriter.add(dataArray);
        xmlWriter.closeTag();

        serializer.add(values);
      }
    }

    xmlWriter.closeTag();
  }

  if (this->__binary) {
    xmlWriter.closeTag();
    xmlWriter.closeTag();
    {
      XMLTag appended("AppendedData");

      ReferenceCounting<XMLAttribute>
        encoding = new XMLAttribute("encoding","raw");

      appended.add(encoding);
      xmlWriter.add(appended);
 
      xmlWriter.insertUnderscore();

      file << serializer;

      xmlWriter.insertNewLine();
    }
  }
}

template <>
void WriterVTK::
__proceed<Structured3DMesh>() const
{
  std::string filename = __filename;
  filename += ".vti";

  const Structured3DMesh& mesh
    = static_cast<const Structured3DMesh&>(*__mesh);

  std::ofstream file(filename.c_str());
  file.precision(15);

  XMLWriter xmlWriter(file);

  xmlWriter.writeHeader();
  {
    XMLTag vtkFile("VTKFile");
    ReferenceCounting<XMLAttribute> attribute
      = new XMLAttribute("type","ImageData");
    vtkFile.add(attribute);

    if (this->__binary) {
      ReferenceCounting<XMLAttribute> endiannes
        = new XMLAttribute("byte_order","LittleEndian");
      vtkFile.add(endiannes);
    }

    xmlWriter.add(vtkFile);
  }

  {
    XMLTag rectilinearGrid("ImageData");
    ReferenceCounting<XMLAttribute> wholeExtent
      = new XMLAttribute("WholeExtent",
                         "0 "+stringify(mesh.shape().nx()-1)+
                         " 0 "+stringify(mesh.shape().ny()-1)+
                         " 0 "+stringify(mesh.shape().nz()-1));
    rectilinearGrid.add(wholeExtent);
    ReferenceCounting<XMLAttribute> origin
      = new XMLAttribute("Origin",
                         stringify(mesh.shape().a(0))+" "+
                         stringify(mesh.shape().a(1))+" "+
                         stringify(mesh.shape().a(2)));
    rectilinearGrid.add(origin);
    ReferenceCounting<XMLAttribute> spacing
      = new XMLAttribute("Spacing",
                         stringify(mesh.shape().hx())+" "+
                         stringify(mesh.shape().hy())+" "+
                         stringify(mesh.shape().hz()));
    rectilinearGrid.add(spacing);
    xmlWriter.add(rectilinearGrid);
  }

  {
    XMLTag piece("Piece");
    ReferenceCounting<XMLAttribute> extent
      = new XMLAttribute("Extent",
                         "0 "+stringify(mesh.shape().nx()-1)
                         +" 0 "+stringify(mesh.shape().ny()-1)
                         +" 0 "+stringify(mesh.shape().nz()-1));
    piece.add(extent);
    xmlWriter.add(piece);
  }

  {
    XMLTag pointData("PointData");
    XMLAttribute scalars("Scalars",__filename);
    xmlWriter.add(pointData);

    for (size_t i=0; i<__scalarFunctionList.size(); ++i) {

      const ScalarFunctionBase& function = *__scalarFunctionList[i];

      if (not this->__binary) {
        XMLTag dataArray("DataArray");
        ReferenceCounting<XMLAttribute> type
          = new XMLAttribute("type","Float64");
        dataArray.add(type);
        ReferenceCounting<XMLAttribute> name
          = new XMLAttribute("Name",stringify(function));
        dataArray.add(name);
        ReferenceCounting<XMLAttribute> format
          = new XMLAttribute("format","ascii");
        dataArray.add(format);
        xmlWriter.add(dataArray);

        switch (function.type()) {
        case ScalarFunctionBase::femfunction: {
          const FEMFunctionBase& fem
            = static_cast<const FEMFunctionBase&>(function);

          if ((fem.discretizationType() == ScalarDiscretizationTypeBase::lagrangianFEM1)
              and (fem.baseMesh() == __mesh)) {
            for (size_t k=0; k<mesh.shape().nz(); ++k)
              for (size_t j=0; j<mesh.shape().ny(); ++j)
                for (size_t i=0; i<mesh.shape().nx(); ++i) {
                  const size_t pointNumber = mesh.shape()(i,j,k);
                  xmlWriter.insert(fem[pointNumber]);
                }
            break;
          }
        }
        default: {
          for (size_t k=0; k<mesh.shape().nz(); ++k)
            for (size_t j=0; j<mesh.shape().ny(); ++j)
              for (size_t i=0; i<mesh.shape().nx(); ++i) {
                const Vertex& x = mesh.vertex(i,j,k);
                xmlWriter.insert(function(x));
              }
        }
        }
        xmlWriter.insertNewLine();
        xmlWriter.closeTag();
      } else {
        XMLTag dataArray("DataArray");
        ReferenceCounting<XMLAttribute> type
          = new XMLAttribute("type","Float64");
        dataArray.add(type);
        ReferenceCounting<XMLAttribute> name
          = new XMLAttribute("Name",stringify(function));
        dataArray.add(name);
        ReferenceCounting<XMLAttribute>
          nbComponent = new XMLAttribute("NumberOfComponents","1");
        dataArray.add(nbComponent);

        ReferenceCounting<XMLAttribute>
          format = new XMLAttribute("format","appended");
        dataArray.add(format);

        ReferenceCounting<XMLAttribute>
          offset = new XMLAttribute("offset",stringify(serializer.offset()));
        dataArray.add(offset);

        xmlWriter.add(dataArray);
        xmlWriter.closeTag();

        Vector<real_t> values(mesh.numberOfVertices());

        switch (function.type()) {
        case ScalarFunctionBase::femfunction: {
          const FEMFunctionBase& fem
            = static_cast<const FEMFunctionBase&>(function);

          if ((fem.discretizationType() == ScalarDiscretizationTypeBase::lagrangianFEM1)
              and (fem.baseMesh() == __mesh)) {
            size_t l=0;
            for (size_t k=0; k<mesh.shape().nz(); ++k)
              for (size_t j=0; j<mesh.shape().ny(); ++j)
                for (size_t i=0; i<mesh.shape().nx(); ++i) {
                  const size_t pointNumber = mesh.shape()(i,j,k);
                  values[l] = fem[pointNumber];
                  l++;
                }
            break;
          }
        }
        default: {
          size_t l=0;
          for (size_t k=0; k<mesh.shape().nz(); ++k)
            for (size_t j=0; j<mesh.shape().ny(); ++j)
              for (size_t i=0; i<mesh.shape().nx(); ++i) {
                const Vertex& x = mesh.vertex(i,j,k);
                values[l] = function(x);
                l++;
              }
        }
        }

        serializer.add(values);
      }
    }

    for (size_t i=0; i<__fieldList.size(); ++i) {

      const FieldOfScalarFunction& field = *__fieldList[i];

      XMLTag dataArray("DataArray");
      ReferenceCounting<XMLAttribute> type
        = new XMLAttribute("type","Float64");
      dataArray.add(type);
      ReferenceCounting<XMLAttribute> name
        = new XMLAttribute("Name",stringify(field));
      dataArray.add(name);
      ReferenceCounting<XMLAttribute>
        nbComponent = new XMLAttribute("NumberOfComponents",
                                       stringify(field.numberOfComponents()));
      dataArray.add(nbComponent);

      Vector<real_t> values(field.numberOfComponents()*__mesh->numberOfVertices());
      this->__fillCrossedComponent(field,values);

      if (not this->__binary) {
        ReferenceCounting<XMLAttribute> format
          = new XMLAttribute("format","ascii");
        dataArray.add(format);
        xmlWriter.add(dataArray);


        for (size_t i=0; i<values.size(); ++i) {
          xmlWriter.insert(values[i]);
        }

        xmlWriter.insertNewLine();
        xmlWriter.closeTag();
      } else {
        ReferenceCounting<XMLAttribute>
          format = new XMLAttribute("format","appended");
        dataArray.add(format);

        ReferenceCounting<XMLAttribute>
          offset = new XMLAttribute("offset",stringify(serializer.offset()));
        dataArray.add(offset);

        xmlWriter.add(dataArray);
        xmlWriter.closeTag();

        serializer.add(values);
      }
    }

    xmlWriter.closeTag(); // Point Data
  }

  xmlWriter.closeTag(); // Piece
  xmlWriter.closeTag(); // ImageData

  if (this->__binary) {
    XMLTag appended("AppendedData");

    ReferenceCounting<XMLAttribute>
      encoding = new XMLAttribute("encoding","raw");

    appended.add(encoding);
    xmlWriter.add(appended);
 
    xmlWriter.insertUnderscore();
    file << serializer;
    xmlWriter.insertNewLine();
  }
}

void WriterVTK::
proceed() const
{
  switch (__mesh->type()) {
  case Mesh::cartesianHexahedraMesh: {
    this->__proceed<Structured3DMesh>();
    break;
  }
  case Mesh::hexahedraMesh: {
    this->__proceed<MeshOfHexahedra>();
    break;
  }
  case Mesh::octreeMesh: {
    this->__proceed<OctreeMesh>();
    break;
  }
  case Mesh::tetrahedraMesh: {
    this->__proceed<MeshOfTetrahedra>();
    break;
  }
  case Mesh::spectralMesh: {
    this->__proceed<SpectralMesh>();
    break;
  }
  case Mesh::surfaceMeshTriangles: {
    this->__proceed<SurfaceMeshOfTriangles>();
    break;
  }
  case Mesh::surfaceMeshQuadrangles: {
    this->__proceed<SurfaceMeshOfQuadrangles>();
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
                       "cannot save this mesh type to VTK format",
                       ErrorHandler::normal);
  }
  }
}

WriterVTK::
WriterVTK(ConstReferenceCounting<Mesh> mesh,
          const std::string& filename,
          const FileDescriptor& fileDescriptor)
  : WriterBase(mesh,
               filename,
               fileDescriptor),
    __binary((fileDescriptor.type() == FileDescriptor::binary) or 
             (fileDescriptor.type() == FileDescriptor::formatDefault))
{
  ;
}

WriterVTK::
~WriterVTK()
{
  ;
}

---