diff --git a/src/TNL/Functions/MeshFunctionGnuplotWriter.h b/src/TNL/Functions/MeshFunctionGnuplotWriter.h
index d747e84a75d6c75d4029b49761ad17efdaf72368..244146ff6d9eea06068f7ac1b61379236fac7e02 100644
--- a/src/TNL/Functions/MeshFunctionGnuplotWriter.h
+++ b/src/TNL/Functions/MeshFunctionGnuplotWriter.h
@@ -68,11 +68,10 @@ template< typename MeshFunction,
 class MeshFunctionGnuplotWriter
 : public MeshFunctionGnuplotWriterBase
 {
-   public:
-
-      using MeshType = typename MeshFunction::MeshType;
-      using EntityType = typename MeshType::template EntityType< MeshFunction::getEntitiesDimension() >;
-      using GlobalIndex = typename MeshType::GlobalIndexType;
+public:
+   using MeshType = typename MeshFunction::MeshType;
+   using EntityType = typename MeshType::template EntityType< MeshFunction::getEntitiesDimension() >;
+   using GlobalIndex = typename MeshType::GlobalIndexType;
 
    static bool write( const MeshFunction& function,
                       std::ostream& str,
@@ -99,11 +98,10 @@ template< typename MeshFunction,
 class MeshFunctionGnuplotWriter< MeshFunction, Meshes::Grid< 2, Real, Device, Index >, EntityDimension >
 : public MeshFunctionGnuplotWriterBase
 {
-   public:
-
-      using MeshType = typename MeshFunction::MeshType;
-      using EntityType = typename MeshType::template EntityType< MeshFunction::getEntitiesDimension() >;
-      using GlobalIndex = typename MeshType::GlobalIndexType;
+public:
+   using MeshType = typename MeshFunction::MeshType;
+   using EntityType = typename MeshType::template EntityType< MeshFunction::getEntitiesDimension() >;
+   using GlobalIndex = typename MeshType::GlobalIndexType;
 
    static bool write( const MeshFunction& function,
                       std::ostream& str,
@@ -137,11 +135,10 @@ template< typename MeshFunction,
 class MeshFunctionGnuplotWriter< MeshFunction, Meshes::Grid< 3, Real, Device, Index >, EntityDimension >
 : public MeshFunctionGnuplotWriterBase
 {
-   public:
-
-      using MeshType = typename MeshFunction::MeshType;
-      using EntityType = typename MeshType::template EntityType< MeshFunction::getEntitiesDimension() >;
-      using GlobalIndex = typename MeshType::GlobalIndexType;
+public:
+   using MeshType = typename MeshFunction::MeshType;
+   using EntityType = typename MeshType::template EntityType< MeshFunction::getEntitiesDimension() >;
+   using GlobalIndex = typename MeshType::GlobalIndexType;
 
    static bool write( const MeshFunction& function,
                       std::ostream& str,
@@ -167,6 +164,5 @@ class MeshFunctionGnuplotWriter< MeshFunction, Meshes::Grid< 3, Real, Device, In
    }
 };
 
-
 } // namespace Functions
 } // namespace TNL
diff --git a/src/TNL/Functions/MeshFunctionVTKWriter.h b/src/TNL/Functions/MeshFunctionVTKWriter.h
index 78608de7461dc5510d280d29b9f7c329836e3eb8..201178c61197da4941f3c06af7914ec2428245b6 100644
--- a/src/TNL/Functions/MeshFunctionVTKWriter.h
+++ b/src/TNL/Functions/MeshFunctionVTKWriter.h
@@ -13,7 +13,7 @@
 #include <TNL/Meshes/Writers/VTKWriter.h>
 
 namespace TNL {
-namespace Functions {   
+namespace Functions {
 
 template< typename MeshFunction >
 class MeshFunctionVTKWriter
diff --git a/src/TNL/Functions/VectorFieldGnuplotWriter.h b/src/TNL/Functions/VectorFieldGnuplotWriter.h
index 41b59d511d680568d62ed545a8f230efc43dd575..a1a63883e8387b1bbe94cb463ae2d39286105570 100644
--- a/src/TNL/Functions/VectorFieldGnuplotWriter.h
+++ b/src/TNL/Functions/VectorFieldGnuplotWriter.h
@@ -16,15 +16,15 @@ namespace TNL {
 namespace Functions {
 
 template< int, typename > class VectorField;
+template< typename, int, typename > class MeshFunction;
 
 template< typename VectorField >
 class VectorFieldGnuplotWriter
 {
-   public:
-
-      static bool write( const VectorField& function,
-                         std::ostream& str,
-                         const double& scale );
+public:
+   static bool write( const VectorField& function,
+                      std::ostream& str,
+                      const double& scale );
 };
 
 /***
@@ -37,14 +37,14 @@ template< typename MeshReal,
           int VectorFieldSize >
 class VectorFieldGnuplotWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 1, Real > > >
 {
-   public:
-      typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 1, RealType > > VectorFieldType;
-
-      static bool write( const VectorFieldType& function,
-                         std::ostream& str,
-                         const double& scale  );
+public:
+   using MeshType = Meshes::Grid< 1, MeshReal, Device, MeshIndex >;
+   using RealType = Real;
+   using VectorFieldType = Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 1, RealType > >;
+
+   static bool write( const VectorFieldType& function,
+                      std::ostream& str,
+                      const double& scale );
 };
 
 /***
@@ -57,14 +57,14 @@ template< typename MeshReal,
           int VectorFieldSize >
 class VectorFieldGnuplotWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 0, Real > > >
 {
-   public:
-      typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 0, RealType > > VectorFieldType;
-
-      static bool write( const VectorFieldType& function,
-                         std::ostream& str,
-                         const double& scale  );
+public:
+   using MeshType = Meshes::Grid< 1, MeshReal, Device, MeshIndex >;
+   using RealType = Real;
+   using VectorFieldType = Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 0, RealType > >;
+
+   static bool write( const VectorFieldType& function,
+                      std::ostream& str,
+                      const double& scale );
 };
 
 
@@ -78,14 +78,14 @@ template< typename MeshReal,
           int VectorFieldSize >
 class VectorFieldGnuplotWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 2, Real > > >
 {
-   public:
-      typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 2, RealType > > VectorFieldType;
-
-      static bool write( const VectorFieldType& function,
-                         std::ostream& str,
-                         const double& scale  );
+public:
+   using MeshType = Meshes::Grid< 2, MeshReal, Device, MeshIndex >;
+   using RealType = Real;
+   using VectorFieldType = Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 2, RealType > >;
+
+   static bool write( const VectorFieldType& function,
+                      std::ostream& str,
+                      const double& scale );
 };
 
 /***
@@ -98,14 +98,14 @@ template< typename MeshReal,
           int VectorFieldSize >
 class VectorFieldGnuplotWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 1, Real > > >
 {
-   public:
-      typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 1, RealType > > VectorFieldType;
-
-      static bool write( const VectorFieldType& function,
-                         std::ostream& str,
-                         const double& scale  );
+public:
+   using MeshType = Meshes::Grid< 2, MeshReal, Device, MeshIndex >;
+   using RealType = Real;
+   using VectorFieldType = Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 1, RealType > >;
+
+   static bool write( const VectorFieldType& function,
+                      std::ostream& str,
+                      const double& scale );
 };
 
 /***
@@ -118,14 +118,14 @@ template< typename MeshReal,
           int VectorFieldSize >
 class VectorFieldGnuplotWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 0, Real > > >
 {
-   public:
-      typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 0, RealType > > VectorFieldType;
-
-      static bool write( const VectorFieldType& function,
-                         std::ostream& str,
-                         const double& scale  );
+public:
+   using MeshType = Meshes::Grid< 2, MeshReal, Device, MeshIndex >;
+   using RealType = Real;
+   using VectorFieldType = Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 0, RealType > >;
+
+   static bool write( const VectorFieldType& function,
+                      std::ostream& str,
+                      const double& scale );
 };
 
 
@@ -139,14 +139,14 @@ template< typename MeshReal,
           int VectorFieldSize >
 class VectorFieldGnuplotWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 3, Real > > >
 {
-   public:
-      typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 3, RealType > > VectorFieldType;
-
-      static bool write( const VectorFieldType& function,
-                         std::ostream& str,
-                         const double& scale  );
+public:
+   using MeshType = Meshes::Grid< 3, MeshReal, Device, MeshIndex >;
+   using RealType = Real;
+   using VectorFieldType = Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 3, RealType > >;
+
+   static bool write( const VectorFieldType& function,
+                      std::ostream& str,
+                      const double& scale );
 };
 
 /***
@@ -159,14 +159,14 @@ template< typename MeshReal,
           int VectorFieldSize >
 class VectorFieldGnuplotWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 2, Real > > >
 {
-   public:
-      typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 2, RealType > > VectorFieldType;
-
-      static bool write( const VectorFieldType& function,
-                         std::ostream& str,
-                         const double& scale  );
+public:
+   using MeshType = Meshes::Grid< 3, MeshReal, Device, MeshIndex >;
+   using RealType = Real;
+   using VectorFieldType = Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 2, RealType > >;
+
+   static bool write( const VectorFieldType& function,
+                      std::ostream& str,
+                      const double& scale );
 };
 
 /***
@@ -179,14 +179,14 @@ template< typename MeshReal,
           int VectorFieldSize >
 class VectorFieldGnuplotWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 0, Real > > >
 {
-   public:
-      typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 0, RealType > > VectorFieldType;
-
-      static bool write( const VectorFieldType& function,
-                         std::ostream& str,
-                         const double& scale  );
+public:
+   using MeshType = Meshes::Grid< 3, MeshReal, Device, MeshIndex >;
+   using RealType = Real;
+   using VectorFieldType = Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 0, RealType > >;
+
+   static bool write( const VectorFieldType& function,
+                      std::ostream& str,
+                      const double& scale );
 };
 
 } // namespace Functions
diff --git a/src/TNL/Functions/VectorFieldGnuplotWriter_impl.h b/src/TNL/Functions/VectorFieldGnuplotWriter_impl.h
index 500bdc4d8fbb0d3ece8c94ec6938c4c81a51b1c8..ebda5972a6a0fbff85726d2ff5abc534afd112cc 100644
--- a/src/TNL/Functions/VectorFieldGnuplotWriter_impl.h
+++ b/src/TNL/Functions/VectorFieldGnuplotWriter_impl.h
@@ -43,9 +43,8 @@ write( const VectorFieldType& vectorField,
 {
    const MeshType& mesh = vectorField.getMesh();
    typename MeshType::Cell entity( mesh );
-   for( entity.getCoordinates().x() = 0;
-        entity.getCoordinates().x() < mesh.getDimensions().x();
-        entity.getCoordinates().x() ++ )
+   auto& c = entity.getCoordinates();
+   for( c.x() = 0; c.x() < mesh.getDimensions().x(); c.x()++ )
    {
       entity.refresh();
       typename MeshType::PointType v = entity.getCenter();
@@ -73,9 +72,8 @@ write( const VectorFieldType& vectorField,
 {
    const MeshType& mesh = vectorField.getMesh();
    typename MeshType::Vertex entity( mesh );
-   for( entity.getCoordinates().x() = 0;
-        entity.getCoordinates().x() <= mesh.getDimensions().x();
-        entity.getCoordinates().x() ++ )
+   auto& c = entity.getCoordinates();
+   for( c.x() = 0; c.x() <= mesh.getDimensions().x(); c.x()++ )
    {
       entity.refresh();
       typename MeshType::PointType v = entity.getCenter();
@@ -104,13 +102,10 @@ write( const VectorFieldType& vectorField,
 {
    const MeshType& mesh = vectorField.getMesh();
    typename MeshType::Cell entity( mesh );
-   for( entity.getCoordinates().y() = 0;
-        entity.getCoordinates().y() < mesh.getDimensions().y();
-        entity.getCoordinates().y() ++ )
+   auto& c = entity.getCoordinates();
+   for( c.y() = 0; c.y() < mesh.getDimensions().y(); c.y()++ )
    {
-      for( entity.getCoordinates().x() = 0;
-           entity.getCoordinates().x() < mesh.getDimensions().x();
-           entity.getCoordinates().x() ++ )
+      for( c.x() = 0; c.x() < mesh.getDimensions().x(); c.x()++ )
       {
          entity.refresh();
          typename MeshType::PointType v = entity.getCenter();
@@ -142,15 +137,12 @@ write( const VectorFieldType& vectorField,
    typedef typename MeshType::Face EntityType;
    typedef typename EntityType::EntityOrientationType EntityOrientation;
    EntityType entity( mesh );
+   auto& c = entity.getCoordinates();
 
    entity.setOrientation( EntityOrientation( 1.0, 0.0 ) );
-   for( entity.getCoordinates().y() = 0;
-        entity.getCoordinates().y() < mesh.getDimensions().y();
-        entity.getCoordinates().y() ++ )
+   for( c.y() = 0; c.y() < mesh.getDimensions().y(); c.y()++ )
    {
-      for( entity.getCoordinates().x() = 0;
-           entity.getCoordinates().x() <= mesh.getDimensions().x();
-           entity.getCoordinates().x() ++ )
+      for( c.x() = 0; c.x() <= mesh.getDimensions().x(); c.x()++ )
       {
          entity.refresh();
          typename MeshType::PointType v = entity.getCenter();
@@ -163,15 +155,9 @@ write( const VectorFieldType& vectorField,
    }
 
    entity.setOrientation( EntityOrientation( 0.0, 1.0 ) );
-         for( entity.getCoordinates().x() = 0;
-           entity.getCoordinates().x() < mesh.getDimensions().x();
-           entity.getCoordinates().x() ++ )
-
+   for( c.x() = 0; c.x() < mesh.getDimensions().x(); c.x()++ )
    {
-            for( entity.getCoordinates().y() = 0;
-        entity.getCoordinates().y() <= mesh.getDimensions().y();
-        entity.getCoordinates().y() ++ )
-
+      for( c.y() = 0; c.y() <= mesh.getDimensions().y(); c.y()++ )
       {
          entity.refresh();
          typename MeshType::PointType v = entity.getCenter();
@@ -202,13 +188,10 @@ write( const VectorFieldType& vectorField,
 {
    const MeshType& mesh = vectorField.getMesh();
    typename MeshType::Vertex entity( mesh );
-   for( entity.getCoordinates().y() = 0;
-        entity.getCoordinates().y() <= mesh.getDimensions().y();
-        entity.getCoordinates().y() ++ )
+   auto& c = entity.getCoordinates();
+   for( c.y() = 0; c.y() <= mesh.getDimensions().y(); c.y()++ )
    {
-      for( entity.getCoordinates().x() = 0;
-           entity.getCoordinates().x() <= mesh.getDimensions().x();
-           entity.getCoordinates().x() ++ )
+      for( c.x() = 0; c.x() <= mesh.getDimensions().x(); c.x()++ )
       {
          entity.refresh();
          typename MeshType::PointType v = entity.getCenter();
@@ -239,16 +222,11 @@ write( const VectorFieldType& vectorField,
 {
    const MeshType& mesh = vectorField.getMesh();
    typename MeshType::Cell entity( mesh );
-   for( entity.getCoordinates().z() = 0;
-        entity.getCoordinates().z() < mesh.getDimensions().z();
-        entity.getCoordinates().z() ++ )
-      for( entity.getCoordinates().y() = 0;
-           entity.getCoordinates().y() < mesh.getDimensions().y();
-           entity.getCoordinates().y() ++ )
+   auto& c = entity.getCoordinates();
+   for( c.z() = 0; c.z() < mesh.getDimensions().z(); c.z()++ )
+      for( c.y() = 0; c.y() < mesh.getDimensions().y(); c.y()++ )
       {
-         for( entity.getCoordinates().x() = 0;
-              entity.getCoordinates().x() < mesh.getDimensions().x();
-              entity.getCoordinates().x() ++ )
+         for( c.x() = 0; c.x() < mesh.getDimensions().x(); c.x()++ )
          {
             entity.refresh();
             typename MeshType::PointType v = entity.getCenter();
@@ -280,18 +258,13 @@ write( const VectorFieldType& vectorField,
    typedef typename MeshType::Face EntityType;
    typedef typename EntityType::EntityOrientationType EntityOrientation;
    EntityType entity( mesh );
+   auto& c = entity.getCoordinates();
 
    entity.setOrientation( EntityOrientation( 1.0, 0.0, 0.0 ) );
-   for( entity.getCoordinates().z() = 0;
-        entity.getCoordinates().z() < mesh.getDimensions().z();
-        entity.getCoordinates().z() ++ )
-      for( entity.getCoordinates().y() = 0;
-           entity.getCoordinates().y() < mesh.getDimensions().y();
-           entity.getCoordinates().y() ++ )
+   for( c.z() = 0; c.z() < mesh.getDimensions().z(); c.z()++ )
+      for( c.y() = 0; c.y() < mesh.getDimensions().y(); c.y()++ )
       {
-         for( entity.getCoordinates().x() = 0;
-              entity.getCoordinates().x() <= mesh.getDimensions().x();
-              entity.getCoordinates().x() ++ )
+         for( c.x() = 0; c.x() <= mesh.getDimensions().x(); c.x()++ )
          {
             entity.refresh();
             typename MeshType::PointType v = entity.getCenter();
@@ -304,16 +277,10 @@ write( const VectorFieldType& vectorField,
       }
 
    entity.setOrientation( EntityOrientation( 0.0, 1.0, 0.0 ) );
-   for( entity.getCoordinates().z() = 0;
-        entity.getCoordinates().z() < mesh.getDimensions().z();
-        entity.getCoordinates().z() ++ )
-      for( entity.getCoordinates().x() = 0;
-           entity.getCoordinates().x() < mesh.getDimensions().x();
-           entity.getCoordinates().x() ++ )
+   for( c.z() = 0; c.z() < mesh.getDimensions().z(); c.z()++ )
+      for( c.x() = 0; c.x() < mesh.getDimensions().x(); c.x()++ )
       {
-         for( entity.getCoordinates().y() = 0;
-              entity.getCoordinates().y() <= mesh.getDimensions().y();
-              entity.getCoordinates().y() ++ )
+         for( c.y() = 0; c.y() <= mesh.getDimensions().y(); c.y()++ )
          {
             entity.refresh();
             typename MeshType::PointType v = entity.getCenter();
@@ -326,16 +293,10 @@ write( const VectorFieldType& vectorField,
       }
 
    entity.setOrientation( EntityOrientation( 0.0, 0.0, 1.0 ) );
-   for( entity.getCoordinates().x() = 0;
-        entity.getCoordinates().x() < mesh.getDimensions().x();
-        entity.getCoordinates().x() ++ )
-      for( entity.getCoordinates().y() = 0;
-           entity.getCoordinates().y() <= mesh.getDimensions().y();
-           entity.getCoordinates().y() ++ )
+   for( c.x() = 0; c.x() < mesh.getDimensions().x(); c.x()++ )
+      for( c.y() = 0; c.y() <= mesh.getDimensions().y(); c.y()++ )
       {
-         for( entity.getCoordinates().z() = 0;
-              entity.getCoordinates().z() < mesh.getDimensions().z();
-              entity.getCoordinates().z() ++ )
+         for( c.z() = 0; c.z() < mesh.getDimensions().z(); c.z()++ )
          {
             entity.refresh();
             typename MeshType::PointType v = entity.getCenter();
@@ -366,16 +327,11 @@ write( const VectorFieldType& vectorField,
 {
    const MeshType& mesh = vectorField.getMesh();
    typename MeshType::Vertex entity( mesh );
-   for( entity.getCoordinates().z() = 0;
-        entity.getCoordinates().z() <= mesh.getDimensions().z();
-        entity.getCoordinates().z() ++ )
-      for( entity.getCoordinates().y() = 0;
-           entity.getCoordinates().y() <= mesh.getDimensions().y();
-           entity.getCoordinates().y() ++ )
+   auto& c = entity.getCoordinates();
+   for( c.z() = 0; c.z() <= mesh.getDimensions().z(); c.z()++ )
+      for( c.y() = 0; c.y() <= mesh.getDimensions().y(); c.y()++ )
       {
-         for( entity.getCoordinates().x() = 0;
-              entity.getCoordinates().x() <= mesh.getDimensions().x();
-              entity.getCoordinates().x() ++ )
+         for( c.x() = 0; c.x() <= mesh.getDimensions().x(); c.x()++ )
          {
             entity.refresh();
             typename MeshType::PointType v = entity.getCenter();
@@ -391,4 +347,3 @@ write( const VectorFieldType& vectorField,
 
 } // namespace Functions
 } // namespace TNL
-
diff --git a/src/TNL/Functions/VectorFieldVTKWriter.h b/src/TNL/Functions/VectorFieldVTKWriter.h
index 6d8b1a8535b25e076c83e688706f37490086c39d..5eceea57fe21d52055294ca3f22a437988e39701 100644
--- a/src/TNL/Functions/VectorFieldVTKWriter.h
+++ b/src/TNL/Functions/VectorFieldVTKWriter.h
@@ -2,7 +2,7 @@
                           VectorFieldVTKWriter.h  -  description
                              -------------------
     begin                : Jan 10, 2018
-    copyright            : (C) 2018 by oberhuber
+    copyright            : (C) 2018 by Tomas Oberhuber et al.
     email                : tomas.oberhuber@fjfi.cvut.cz
  ***************************************************************************/
 
@@ -10,255 +10,52 @@
 
 #pragma once
 
-#include <TNL/Meshes/Grid.h>
+#include <TNL/Meshes/Writers/VTKWriter.h>
 
 namespace TNL {
 namespace Functions {
 
-template< int, typename > class VectorField;
-
 template< typename VectorField >
 class VectorFieldVTKWriter
 {
-   public:
-
-      static bool write( const VectorField& vectorField,
-                         std::ostream& str,
-                         const double& scale );
-      
-      static void writeHeader( const VectorField& vectorField,
-                               std::ostream& str ){}
-      
-};
-
-/***
- * 1D grids cells
- */
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-class VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 1, Real > > >
-{
-   public:
-      typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 1, RealType > > VectorFieldType;
-      using VectorType = typename VectorFieldType::VectorType;
-
-      static bool write( const VectorFieldType& function,
-                         std::ostream& str,
-                         const double& scale  );
-      
-      static void writeHeader( const VectorFieldType& vectorField,
-                               std::ostream& str );
-      
-};
-
-/***
- * 1D grids vertices
- */
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-class VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 0, Real > > >
-{
-   public:
-      typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 0, RealType > > VectorFieldType;
-      using VectorType = typename VectorFieldType::VectorType;
-
-      static bool write( const VectorFieldType& function,
-                         std::ostream& str,
-                         const double& scale  );
-      
-      static void writeHeader( const VectorFieldType& vectorField,
-                               std::ostream& str );
-      
-};
-
-
-/***
- * 2D grids cells
- */
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-class VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 2, Real > > >
-{
-   public:
-      typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 2, RealType > > VectorFieldType;
-      using VectorType = typename VectorFieldType::VectorType;
-
-      static bool write( const VectorFieldType& function,
-                         std::ostream& str,
-                         const double& scale  );
-
-      static void writeHeader( const VectorFieldType& vectorField,
-                               std::ostream& str );
-      
-};
-
-/***
- * 2D grids faces
- */
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-class VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 1, Real > > >
-{
-   public:
-      typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 1, RealType > > VectorFieldType;
-      using VectorType = typename VectorFieldType::VectorType;
-
-      static bool write( const VectorFieldType& function,
-                         std::ostream& str,
-                         const double& scale  );
-      
-      static void writeHeader( const VectorFieldType& vectorField,
-                               std::ostream& str );
-      
-};
-
-/***
- * 2D grids vertices
- */
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-class VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 0, Real > > >
-{
-   public:
-      typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 0, RealType > > VectorFieldType;
-      using VectorType = typename VectorFieldType::VectorType;
-
-      static bool write( const VectorFieldType& function,
-                         std::ostream& str,
-                         const double& scale  );
-      
-      static void writeHeader( const VectorFieldType& vectorField,
-                               std::ostream& str );
-      
-};
-
-
-/***
- * 3D grids cells
- */
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-class VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 3, Real > > >
-{
-   public:
-      typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 3, RealType > > VectorFieldType;
-      using VectorType = typename VectorFieldType::VectorType;
-
-      static bool write( const VectorFieldType& function,
-                         std::ostream& str,
-                         const double& scale  );
-      
-      static void writeHeader( const VectorFieldType& vectorField,
-                               std::ostream& str );
-      
-};
-
-/***
- * 3D grids faces
- */
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-class VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 2, Real > > >
-{
-   public:
-      typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 2, RealType > > VectorFieldType;
-      using VectorType = typename VectorFieldType::VectorType;
-
-      static bool write( const VectorFieldType& function,
-                         std::ostream& str,
-                         const double& scale  );
-      
-      static void writeHeader( const VectorFieldType& vectorField,
-                               std::ostream& str );
-      
-};
-
-/***
- * 3D grids edges
- */
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-class VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 1, Real > > >
-{
-   public:
-      typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 1, RealType > > VectorFieldType;
-      using VectorType = typename VectorFieldType::VectorType;
-
-      static bool write( const VectorFieldType& function,
-                         std::ostream& str,
-                         const double& scale  );
-      
-      static void writeHeader( const VectorFieldType& vectorField,
-                               std::ostream& str );
-      
-};
-
-/***
- * 3D grids vertices
- */
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-class VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 0, Real > > >
-{
-   public:
-      typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::VectorField< VectorFieldSize, MeshFunction< MeshType, 0, RealType > > VectorFieldType;
-      using VectorType = typename VectorFieldType::VectorType;
-
-      static bool write( const VectorFieldType& function,
-                         std::ostream& str,
-                         const double& scale  );
-      
-      static void writeHeader( const VectorFieldType& vectorField,
-                               std::ostream& str );
-      
+   using MeshType = typename VectorField::MeshType;
+   using MeshWriter = Meshes::Writers::VTKWriter< MeshType >;
+   using EntityType = typename MeshType::template EntityType< VectorField::getEntitiesDimension() >;
+   using GlobalIndex = typename MeshType::GlobalIndexType;
+
+public:
+   static bool write( const VectorField& field,
+                      std::ostream& str,
+                      const double& scale = 1.0,
+                      const String& fieldName = "cellVectorFieldValues" )
+   {
+      const MeshType& mesh = field.getMesh();
+      MeshWriter::template writeEntities< VectorField::getEntitiesDimension() >( mesh, str );
+      appendField( field, str, fieldName, scale );
+      return true;
+   }
+
+   // VTK supports writing multiple fields into the same file.
+   // You can call this after 'write', which initializes the mesh entities,
+   // with different field name.
+   static void appendField( const VectorField& field,
+                            std::ostream& str,
+                            const String& fieldName,
+                            const double& scale = 1.0 )
+   {
+      const MeshType& mesh = field.getMesh();
+      const GlobalIndex entitiesCount = mesh.template getEntitiesCount< EntityType >();
+      str << std::endl << "CELL_DATA " << entitiesCount << std::endl;
+      str << "VECTORS " << fieldName << " " << getType< typename VectorField::RealType >() << " 1" << std::endl;
+      for( GlobalIndex i = 0; i < entitiesCount; i++ ) {
+         const typename VectorField::VectorType vector = field.getElement( i );
+         static_assert( VectorField::getVectorDimension() <= 3, "The VTK format supports only up to 3D vector fields." );
+         for( int i = 0; i < 3; i++ )
+            str << scale * ( i < vector.getSize() ? vector[ i ] : 0.0 ) << " ";
+         str << "\n";
+      }
+   }
 };
 
 } // namespace Functions
 } // namespace TNL
-
-#include <TNL/Functions/VectorFieldVTKWriter_impl.h>
diff --git a/src/TNL/Functions/VectorFieldVTKWriter_impl.h b/src/TNL/Functions/VectorFieldVTKWriter_impl.h
deleted file mode 100644
index 938227d22b57f61d2f6c4d5f4b7b13a9044d3aa8..0000000000000000000000000000000000000000
--- a/src/TNL/Functions/VectorFieldVTKWriter_impl.h
+++ /dev/null
@@ -1,881 +0,0 @@
-/***************************************************************************
-                          VectorFieldVTKWriter_impl.h  -  description
-                             -------------------
-    begin                : Jan 10, 2018
-    copyright            : (C) 2018 by oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/* See Copyright Notice in tnl/Copyright */
-
-#pragma once
-
-#include <TNL/Functions/VectorFieldVTKWriter.h>
-#include <TNL/Functions/VectorField.h>
-
-namespace TNL {
-namespace Functions {   
-
-template< typename VectorField >
-bool
-VectorFieldVTKWriter< VectorField >::
-write( const VectorField& vectorField,
-       std::ostream& str,
-       const double& scale )
-{
-   std::cerr << "VTK writer for vector field defined on mesh type " << VectorField::MeshType::getType() << " is not (yet) implemented." << std::endl;
-   return false;
-}
-
-/****
- * 1D grid, cells
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-void
-VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 1, Real > > >::
-writeHeader( const VectorFieldType& vectorField,
-             std::ostream& str )
-{
-    const MeshType& mesh = vectorField.getMesh();
-    const typename MeshType::PointType& origin = mesh.getOrigin();
-    const typename MeshType::PointType& proportions = mesh.getProportions();
-    str << "# vtk DataFile Version 2.0" << std::endl;
-    str << "TNL DATA" << std::endl;
-    str << "ASCII" << std::endl;
-    str << "DATASET UNSTRUCTURED_GRID" << std::endl;
-}
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-bool
-VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 1, Real > > >::
-write( const VectorFieldType& vectorField,
-       std::ostream& str,
-       const double& scale )
-{
-   writeHeader(vectorField, str);
- 
-   const MeshType& mesh = vectorField.getMesh();
-   const RealType origin = mesh.getOrigin().x();
-   const RealType spaceStep = mesh.getSpaceSteps().x();
- 
-   str << "POINTS " << mesh.getDimensions().x() + 1 << " float" << std::endl;
-   for (int i = 0; i <= mesh.getDimensions().x(); i++)
-   {
-       str << origin + i * spaceStep << " 0 0" << std::endl;
-   }
- 
-   str << std::endl << "CELLS " << mesh.getDimensions().x() << " " << mesh.getDimensions().x() * 3 << std::endl;
-   for (int i = 0; i < mesh.getDimensions().x(); i++)
-   {
-       str << "2 " << i << " " << i+1 << std::endl;
-   }
- 
-   str << std::endl << "CELL_TYPES " << mesh.getDimensions().x() << std::endl;
-   for (int i = 0; i < mesh.getDimensions().x(); i++)
-   {
-       str << "3 " << std::endl;
-   }
- 
-   str << std::endl << "CELL_DATA " << mesh.getDimensions().x() << std::endl;
-   str << "VECTORS cellVectorFieldValues " << getType< typename VectorFieldType::RealType >() << std::endl;
-
-   for( MeshIndex i = 0; i < mesh.template getEntitiesCount< typename MeshType::Cell >(); i++ )
-   {
-      typename MeshType::Cell entity = mesh.template getEntity< typename MeshType::Cell >( i );
-      entity.refresh();
-      const VectorType v = vectorField.getElement( entity.getIndex() );
-      for( int i = 0; i < 3; i++ )
-         str << scale * ( i < VectorFieldSize ? v[ i ] : 0.0 ) << " ";
-      str << std::endl;
-   }
- 
-   return true;
-}
-
-/****
- * 1D grid, vertices
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-void
-VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 0, Real > > >::
-writeHeader( const VectorFieldType& vectorField,
-             std::ostream& str )
-{
-    const MeshType& mesh = vectorField.getMesh();
-    const typename MeshType::PointType& origin = mesh.getOrigin();
-    const typename MeshType::PointType& proportions = mesh.getProportions();
-    str << "# vtk DataFile Version 2.0" << std::endl;
-    str << "TNL DATA" << std::endl;
-    str << "ASCII" << std::endl;
-    str << "DATASET UNSTRUCTURED_GRID" << std::endl;
-}
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-bool
-VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 0, Real > > >::
-write( const VectorFieldType& vectorField,
-       std::ostream& str,
-       const double& scale )
-{
-   writeHeader(vectorField, str);
- 
-   const MeshType& mesh = vectorField.getMesh();
-   const RealType origin = mesh.getOrigin().x();
-   const RealType spaceStep = mesh.getSpaceSteps().x();
- 
-   str << "POINTS " << mesh.getDimensions().x() + 1 << " float" << std::endl;
-   for (int i = 0; i < mesh.getDimensions().x() + 1; i++)
-   {
-       str << origin + i * spaceStep << " 0 0" << std::endl;
-   }
- 
-   str << std::endl << "CELLS " << mesh.getDimensions().x() + 1 << " " << ( mesh.getDimensions().x() + 1 ) * 2 << std::endl;
-   for (int i = 0; i < mesh.getDimensions().x() + 1; i++)
-   {
-       str << "1 " << i << std::endl;
-   }
- 
-   str << std::endl << "CELL_TYPES " << mesh.getDimensions().x() + 1 << std::endl;
-   for (int i = 0; i < mesh.getDimensions().x() + 1; i++)
-   {
-       str << "1 " << std::endl;
-   }
- 
-   str << std::endl << "CELL_DATA " << mesh.getDimensions().x() + 1 << std::endl;
-   str << "VECTORS VerticesVectorFieldValues " << getType< typename VectorFieldType::RealType >() << std::endl;
-
-   for( MeshIndex i = 0; i < mesh.template getEntitiesCount< typename MeshType::Vertex >(); i++ )
-   {
-      typename MeshType::Vertex entity = mesh.template getEntity< typename MeshType::Vertex >( i );
-      entity.refresh();
-      const VectorType v = vectorField.getElement( entity.getIndex() );
-      for( int i = 0; i < 3; i++ )
-         str << scale * ( i < VectorFieldSize ? v[ i ] : 0.0 ) << " ";
-      str << std::endl;
-   }
- 
-   return true;
-}
-
-/****
- * 2D grid, cells
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-void
-VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 2, Real > > >::
-writeHeader( const VectorFieldType& vectorField,
-             std::ostream& str )
-{
-    const MeshType& mesh = vectorField.getMesh();
-    const typename MeshType::PointType& origin = mesh.getOrigin();
-    const typename MeshType::PointType& proportions = mesh.getProportions();
-    str << "# vtk DataFile Version 2.0" << std::endl;
-    str << "TNL DATA" << std::endl;
-    str << "ASCII" << std::endl;
-    str << "DATASET UNSTRUCTURED_GRID" << std::endl;
-}
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-bool
-VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 2, Real > > >::
-write( const VectorFieldType& vectorField,
-       std::ostream& str,
-       const double& scale )
-{
-   writeHeader(vectorField, str);
- 
-   const MeshType& mesh = vectorField.getMesh();
-   const RealType originX = mesh.getOrigin().x();
-   const RealType spaceStepX = mesh.getSpaceSteps().x();
-   const RealType originY = mesh.getOrigin().y();
-   const RealType spaceStepY = mesh.getSpaceSteps().y();
-   const MeshIndex verticesCount = mesh.template getEntitiesCount< typename MeshType::Vertex >();
-   const MeshIndex entitiesCount = mesh.template getEntitiesCount< typename MeshType::Cell >();
-
-   str << "POINTS " << verticesCount << " " << getType< RealType >() << std::endl;
-   for (int j = 0; j < mesh.getDimensions().y() + 1; j++)
-   {
-        for (int i = 0; i < mesh.getDimensions().x() + 1; i++)
-        {
-             str << originX + i * spaceStepX << " " << originY + j * spaceStepY << " 0" << std::endl;
-        }
-   }
- 
-   str << std::endl << "CELLS " << entitiesCount << " " << entitiesCount * 5 << std::endl;
-   for (int j = 0; j < mesh.getDimensions().y(); j++)
-   {
-        for (int i = 0; i < mesh.getDimensions().x(); i++)
-        {
-            str << "4 " << j * ( mesh.getDimensions().x() + 1 ) + i << " " << j * ( mesh.getDimensions().x() + 1 )+ i + 1 <<
-                   " " << (j+1) * ( mesh.getDimensions().x() + 1 ) + i << " " << (j+1) * ( mesh.getDimensions().x() + 1 ) + i + 1 << std::endl;
-        }
-   }
- 
-   str << std::endl << "CELL_TYPES " << mesh.getDimensions().x() * mesh.getDimensions().y() << std::endl;
-   for (int i = 0; i < mesh.getDimensions().x()*mesh.getDimensions().y(); i++)
-   {
-       str << "8 " << std::endl;
-   }
- 
-   str << std::endl << "CELL_DATA " << entitiesCount << std::endl;
-   str << "VECTORS cellVectorFieldValues " << getType< RealType >() << std::endl;
-
-   for( MeshIndex i = 0; i < entitiesCount; i++ )
-   {
-      typename MeshType::Cell entity = mesh.template getEntity< typename MeshType::Cell >( i );
-      entity.refresh();
-      const VectorType v = vectorField.getElement( entity.getIndex() );
-      for( int i = 0; i < 3; i++ )
-      {
-         str << scale * ( i < VectorFieldSize ? v[ i ] : 0.0 );
-         if( i < 2 )
-            str << " ";
-      }
-      str << std::endl;
-   }
-
-   return true;
-}
-
-/****
- * 2D grid, faces
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-void
-VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 1, Real > > >::
-writeHeader( const VectorFieldType& vectorField,
-             std::ostream& str )
-{
-    const MeshType& mesh = vectorField.getMesh();
-    const typename MeshType::PointType& origin = mesh.getOrigin();
-    const typename MeshType::PointType& proportions = mesh.getProportions();
-    str << "# vtk DataFile Version 2.0" << std::endl;
-    str << "TNL DATA" << std::endl;
-    str << "ASCII" << std::endl;
-    str << "DATASET UNSTRUCTURED_GRID" << std::endl;
-}
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-bool
-VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 1, Real > > >::
-write( const VectorFieldType& vectorField,
-       std::ostream& str,
-       const double& scale )
-{
-   typedef typename MeshType::template EntityType< 0 > Vertex;
-   typedef typename MeshType::template EntityType< 1 > Face;
-   writeHeader(vectorField, str);
- 
-   const MeshType& mesh = vectorField.getMesh();
-   const RealType originX = mesh.getOrigin().x();
-   const RealType spaceStepX = mesh.getSpaceSteps().x();
-   const RealType originY = mesh.getOrigin().y();
-   const RealType spaceStepY = mesh.getSpaceSteps().y();
-   const MeshIndex verticesCount = mesh.template getEntitiesCount< typename MeshType::Vertex >();
-   const MeshIndex entitiesCount = mesh.template getEntitiesCount< typename MeshType::Face >();
- 
-   str << "POINTS " << verticesCount << " float" << std::endl;
-   for (int j = 0; j < ( mesh.getDimensions().y() + 1); j++)
-   {
-        for (int i = 0; i < ( mesh.getDimensions().x() + 1 ); i++)
-        {
-             str << originX + i * spaceStepX << " " << originY + j * spaceStepY << " 0" << std::endl;
-        }
-   }
- 
-   str << std::endl << "CELLS " << entitiesCount << " " << entitiesCount * 3 << std::endl;
-   for (int j = 0; j < mesh.getDimensions().y(); j++)
-   {
-        for (int i = 0; i < ( mesh.getDimensions().x() + 1 ); i++)
-        {
-            str << "2 " << j * ( mesh.getDimensions().x() + 1 ) + i << " " << (j+1) * ( mesh.getDimensions().x() + 1 ) + i << std::endl;
-        }
-   }
- 
-   for (int j = 0; j < (mesh.getDimensions().y()+1); j++)
-   {
-        for (int i = 0; i < mesh.getDimensions().x(); i++)
-        {
-            str << "2 " << j * ( mesh.getDimensions().x() + 1 ) + i << " " <<j * ( mesh.getDimensions().x() + 1 ) + i + 1<< std::endl;
-        }
-   }
- 
-   str << std::endl << "CELL_TYPES " << entitiesCount << std::endl;
-   for (int i = 0; i < entitiesCount; i++)
-   {
-       str << "3" << std::endl;
-   }
- 
-   str << std::endl << "CELL_DATA " << entitiesCount << std::endl;
-   str << "VECTORS FaceslVectorFieldValues " << getType< typename VectorFieldType::RealType >() << std::endl;
-
-   for( MeshIndex i = 0; i < entitiesCount; i++ )
-   {
-      typename MeshType::Face entity = mesh.template getEntity< typename MeshType::Face >( i );
-      entity.refresh();
-      const VectorType v = vectorField.getElement( entity.getIndex() );
-      for( int i = 0; i < 3; i++ )
-         str << scale * ( i < VectorFieldSize ? v[ i ] : 0.0 ) << " ";
-      str << std::endl;
-   }
-
-   return true;
-}
-
-/****
- * 2D grid, vertices
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-void
-VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 0, Real > > >::
-writeHeader( const VectorFieldType& vectorField,
-             std::ostream& str )
-{
-    const MeshType& mesh = vectorField.getMesh();
-    const typename MeshType::PointType& origin = mesh.getOrigin();
-    const typename MeshType::PointType& proportions = mesh.getProportions();
-    str << "# vtk DataFile Version 2.0" << std::endl;
-    str << "TNL DATA" << std::endl;
-    str << "ASCII" << std::endl;
-    str << "DATASET UNSTRUCTURED_GRID" << std::endl;
-}
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-bool
-VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 0, Real > > >::
-write( const VectorFieldType& vectorField,
-       std::ostream& str,
-       const double& scale )
-{
-   typedef typename MeshType::template EntityType< 0 > Vertex;
-   writeHeader(vectorField, str);
- 
-   const MeshType& mesh = vectorField.getMesh();
-   const RealType originX = mesh.getOrigin().x();
-   const RealType spaceStepX = mesh.getSpaceSteps().x();
-   const RealType originY = mesh.getOrigin().y();
-   const RealType spaceStepY = mesh.getSpaceSteps().y();
-   const MeshIndex verticesCount = mesh.template getEntitiesCount< typename MeshType::Vertex >();
- 
-   str << "POINTS " << verticesCount << " float" << std::endl;
-   for (int j = 0; j < ( mesh.getDimensions().y() + 1); j++)
-   {
-        for (int i = 0; i < ( mesh.getDimensions().x() + 1 ); i++)
-        {
-             str << originX + i * spaceStepX << " " << originY + j * spaceStepY << " 0" << std::endl;
-        }
-   }
- 
-   str << std::endl << "CELLS " << verticesCount << " " << verticesCount * 2 << std::endl;
-   for (int j = 0; j < ( mesh.getDimensions().y() + 1 ); j++)
-   {
-        for (int i = 0; i < ( mesh.getDimensions().x() + 1 ); i++)
-        {
-            str << "1 " << j * mesh.getDimensions().x() + i  << std::endl;
-        }
-   }
- 
-   str << std::endl << "CELL_TYPES " << verticesCount << std::endl;
-   for (int i = 0; i < verticesCount; i++)
-   {
-       str << "1" << std::endl;
-   }
- 
-   str << std::endl << "CELL_DATA " << verticesCount << std::endl;
-   str << "VECTORS VerticesVectorFieldValues " << getType< typename VectorFieldType::RealType >() << std::endl;
-
-   for( MeshIndex i = 0; i < verticesCount; i++ )
-   {
-      typename MeshType::Vertex entity = mesh.template getEntity< typename MeshType::Vertex >( i );
-      entity.refresh();
-      const VectorType v = vectorField.getElement( entity.getIndex() );
-      for( int i = 0; i < 3; i++ )
-         str << scale * ( i < VectorFieldSize ? v[ i ] : 0.0 ) << " ";
-      str << std::endl;
-   }
-
-   return true;
-}
-
-/****
- * 3D grid, cells
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-void
-VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 3, Real > > >::
-writeHeader( const VectorFieldType& vectorField,
-             std::ostream& str )
-{
-    const MeshType& mesh = vectorField.getMesh();
-    const typename MeshType::PointType& origin = mesh.getOrigin();
-    const typename MeshType::PointType& proportions = mesh.getProportions();
-    str << "# vtk DataFile Version 2.0" << std::endl;
-    str << "TNL DATA" << std::endl;
-    str << "ASCII" << std::endl;
-    str << "DATASET UNSTRUCTURED_GRID" << std::endl;
-}
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-bool
-VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 3, Real > > >::
-write( const VectorFieldType& vectorField,
-       std::ostream& str,
-       const double& scale )
-{
-   writeHeader(vectorField, str);
- 
-   const MeshType& mesh = vectorField.getMesh();
-   const RealType originX = mesh.getOrigin().x();
-   const RealType spaceStepX = mesh.getSpaceSteps().x();
-   const RealType originY = mesh.getOrigin().y();
-   const RealType spaceStepY = mesh.getSpaceSteps().y();
-   const RealType originZ = mesh.getOrigin().z();
-   const RealType spaceStepZ = mesh.getSpaceSteps().z();
-   const MeshIndex verticesCount = mesh.template getEntitiesCount< typename MeshType::Vertex >();
-   const MeshIndex entitiesCount = mesh.template getEntitiesCount< typename MeshType::Cell >();
- 
-   str << "POINTS " << verticesCount << " float" << std::endl;
-   for (int k = 0; k <= mesh.getDimensions().y(); k++)
-   {
-       for (int j = 0; j <= mesh.getDimensions().y(); j++)
-       {
-            for (int i = 0; i <= mesh.getDimensions().x(); i++)
-            {
-                 str << originX + i * spaceStepX << " " << originY + j * spaceStepY << " " <<
-                        originZ + k * spaceStepZ << std::endl;
-            }
-       }
-   }
- 
-   str << std::endl << "CELLS " << entitiesCount << " " <<
-          entitiesCount * 9 << std::endl;
-   for (int k = 0; k < mesh.getDimensions().z(); k++)
-   {
-        for (int j = 0; j < mesh.getDimensions().y(); j++)
-        {
-            for (int i = 0; i < mesh.getDimensions().x(); i++)
-            {
-                str << "8 " <<  k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << " "
-                    << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i + 1 << " "
-                    << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + (j+1) * ( mesh.getDimensions().x() + 1 ) + i << " "
-                    << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + (j+1) * ( mesh.getDimensions().x() + 1 ) + i + 1 << " "
-                    << (k+1) * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << " "
-                    << (k+1) * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i + 1 << " "
-                    << (k+1) * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + (j+1) * ( mesh.getDimensions().x() + 1 ) + i << " "
-                    << (k+1) * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + (j+1) * ( mesh.getDimensions().x() + 1 ) + i + 1 << std::endl;
-            }
-        }
-   }
- 
-   str << std::endl << "CELL_TYPES " << entitiesCount << std::endl;
-   for (int i = 0; i < entitiesCount; i++)
-   {
-       str << "11" << std::endl;
-   }
- 
-   str << std::endl << "CELL_DATA " << entitiesCount << std::endl;
-   str << "VECTORS cellVectorFieldValues " << getType< typename VectorFieldType::RealType >() << std::endl;
-
-   for( MeshIndex i = 0; i < entitiesCount; i++ )
-   {
-      typename MeshType::Cell entity = mesh.template getEntity< typename MeshType::Cell >( i );
-      entity.refresh();
-      const VectorType v = vectorField.getElement( entity.getIndex() );
-      for( int i = 0; i < 3; i++ )
-         str << scale * ( i < VectorFieldSize ? v[ i ] : 0.0 ) << " ";
-      str << std::endl;
-   }
-
-   return true;
-}
-
-/****
- * 3D grid, faces
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-void
-VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 2, Real > > >::
-writeHeader( const VectorFieldType& vectorField,
-             std::ostream& str )
-{
-    const MeshType& mesh = vectorField.getMesh();
-    const typename MeshType::PointType& origin = mesh.getOrigin();
-    const typename MeshType::PointType& proportions = mesh.getProportions();
-    str << "# vtk DataFile Version 2.0" << std::endl;
-    str << "TNL DATA" << std::endl;
-    str << "ASCII" << std::endl;
-    str << "DATASET UNSTRUCTURED_GRID" << std::endl;
-}
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-bool
-VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 2, Real > > >::
-write( const VectorFieldType& vectorField,
-       std::ostream& str,
-       const double& scale )
-{
-   writeHeader(vectorField, str);
- 
-   const MeshType& mesh = vectorField.getMesh();
-   const RealType originX = mesh.getOrigin().x();
-   const RealType spaceStepX = mesh.getSpaceSteps().x();
-   const RealType originY = mesh.getOrigin().y();
-   const RealType spaceStepY = mesh.getSpaceSteps().y();
-   const RealType originZ = mesh.getOrigin().z();
-   const RealType spaceStepZ = mesh.getSpaceSteps().z();
-   const MeshIndex verticesCount = mesh.template getEntitiesCount< typename MeshType::Vertex >();
-   const MeshIndex entitiesCount = mesh.template getEntitiesCount< typename MeshType::Face >();
- 
-   str << "POINTS " << verticesCount << " float" << std::endl;
-   for (int k = 0; k <= mesh.getDimensions().y(); k++)
-   {
-       for (int j = 0; j <= mesh.getDimensions().y(); j++)
-       {
-            for (int i = 0; i <= mesh.getDimensions().x(); i++)
-            {
-                 str << originX + i * spaceStepX << " " << originY + j * spaceStepY << " " <<
-                        originZ + k * spaceStepZ << std::endl;
-            }
-       }
-   }
- 
-   str << std::endl << "CELLS " << entitiesCount << " " << entitiesCount * 5 << std::endl;
-   for (int k = 0; k < mesh.getDimensions().z(); k++)
-   {
-        for (int j = 0; j < mesh.getDimensions().y(); j++)
-        {
-            for (int i = 0; i <= mesh.getDimensions().x(); i++)
-            {
-                str << "4 " << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << " "
-                    << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + (j+1) * ( mesh.getDimensions().x() + 1 ) + i << " "
-                    << (k+1) * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << " "
-                    << (k+1) * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + (j+1) * ( mesh.getDimensions().x() + 1 ) + i << std::endl;
-            }
-        }
-   }
- 
-   for (int k = 0; k < mesh.getDimensions().z(); k++)
-   {
-        for (int j = 0; j <= mesh.getDimensions().y(); j++)
-        {
-            for (int i = 0; i < mesh.getDimensions().x(); i++)
-            {
-                str << "4 " << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << " "
-                    << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i + 1 << " "
-                    << (k+1) * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << " "
-                    << (k+1) * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i + 1 << std::endl;
-            }
-        }
-   }
- 
-   for (int k = 0; k <= mesh.getDimensions().z(); k++)
-   {
-        for (int j = 0; j < mesh.getDimensions().y(); j++)
-        {
-            for (int i = 0; i < mesh.getDimensions().x(); i++)
-            {
-                str << "4 " << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << " "
-                    << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i + 1 << " "
-                    << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + (j+1) * ( mesh.getDimensions().x() + 1 ) + i << " "
-                    << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + (j+1) * ( mesh.getDimensions().x() + 1 ) + i + 1<< std::endl;
-            }
-        }
-   }
- 
-   str << std::endl << "CELL_TYPES " << entitiesCount << std::endl;
-   for (int i = 0; i < entitiesCount; i++)
-   {
-       str << "8" << std::endl;
-   }
- 
-   str << std::endl << "CELL_DATA " << entitiesCount << std::endl;
-   str << "VECTORS facesVectorFieldValues " << getType< typename VectorFieldType::RealType >() << std::endl;
-
-   for( MeshIndex i = 0; i < entitiesCount; i++ )
-   {
-      typename MeshType::Face entity = mesh.template getEntity< typename MeshType::Face >( i );
-      entity.refresh();
-      const VectorType v = vectorField.getElement( entity.getIndex() );
-      for( int i = 0; i < 3; i++ )
-         str << scale * ( i < VectorFieldSize ? v[ i ] : 0.0 ) << " ";
-      str << std::endl;
-   }
-
-   return true;
-}
-
-/****
- * 3D grid, edges
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-void
-VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 1, Real > > >::
-writeHeader( const VectorFieldType& vectorField,
-             std::ostream& str )
-{
-    const MeshType& mesh = vectorField.getMesh();
-    const typename MeshType::PointType& origin = mesh.getOrigin();
-    const typename MeshType::PointType& proportions = mesh.getProportions();
-    str << "# vtk DataFile Version 2.0" << std::endl;
-    str << "TNL DATA" << std::endl;
-    str << "ASCII" << std::endl;
-    str << "DATASET UNSTRUCTURED_GRID" << std::endl;
-}
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-bool
-VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 1, Real > > >::
-write( const VectorFieldType& vectorField,
-       std::ostream& str,
-       const double& scale )
-{
-   writeHeader(vectorField, str);
- 
-   const MeshType& mesh = vectorField.getMesh();
-   const RealType originX = mesh.getOrigin().x();
-   const RealType spaceStepX = mesh.getSpaceSteps().x();
-   const RealType originY = mesh.getOrigin().y();
-   const RealType spaceStepY = mesh.getSpaceSteps().y();
-   const RealType originZ = mesh.getOrigin().z();
-   const RealType spaceStepZ = mesh.getSpaceSteps().z();
-   const MeshIndex verticesCount = mesh.template getEntitiesCount< typename MeshType::Vertex >();
-   const MeshIndex entitiesCount = mesh.template getEntitiesCount< typename MeshType::Edge >();
- 
-   str << "POINTS " << verticesCount << " float" << std::endl;
-   for (int k = 0; k <= mesh.getDimensions().y(); k++)
-   {
-       for (int j = 0; j <= mesh.getDimensions().y(); j++)
-       {
-            for (int i = 0; i <= mesh.getDimensions().x(); i++)
-            {
-                 str << originX + i * spaceStepX << " " << originY + j * spaceStepY << " " <<
-                        originZ + k * spaceStepZ << std::endl;
-            }
-       }
-   }
- 
-   str << std::endl << "CELLS " << entitiesCount << " " << entitiesCount * 3 << std::endl;
-   for (int k = 0; k <= mesh.getDimensions().z(); k++)
-   {
-        for (int j = 0; j <= mesh.getDimensions().y(); j++)
-        {
-            for (int i = 0; i < mesh.getDimensions().x(); i++)
-            {
-                str << "3 " << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << " "
-                    << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i + 1 << std::endl;
-            }
-        }
-   }
- 
-   for (int k = 0; k <= mesh.getDimensions().z(); k++)
-   {
-        for (int j = 0; j < mesh.getDimensions().y(); j++)
-        {
-            for (int i = 0; i <= mesh.getDimensions().x(); i++)
-            {
-                str << "3 " << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << " "
-                    << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + (j+1) * ( mesh.getDimensions().x() + 1 ) + i << std::endl;
-            }
-        }
-   }
- 
-   for (int k = 0; k < mesh.getDimensions().z(); k++)
-   {
-        for (int j = 0; j <= mesh.getDimensions().y(); j++)
-        {
-            for (int i = 0; i <= mesh.getDimensions().x(); i++)
-            {
-                str << "3 " << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << " "
-                    << (k+1) * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i << std::endl;
-            }
-        }
-   }
- 
-   str << std::endl << "CELL_TYPES " << entitiesCount << std::endl;
-   for (int i = 0; i < entitiesCount; i++)
-   {
-       str << "3" << std::endl;
-   }
- 
-   str << std::endl << "CELL_DATA " << entitiesCount << std::endl;
-   str << "VECTORS edgesVectorFieldValues " << getType< typename VectorFieldType::RealType >() << std::endl;
-
-   for( MeshIndex i = 0; i < entitiesCount; i++ )
-   {
-      typename MeshType::Edge entity = mesh.template getEntity< typename MeshType::Edge >( i );
-      entity.refresh();
-      const VectorType v = vectorField.getElement( entity.getIndex() );
-      for( int i = 0; i < 3; i++ )
-         str << scale * ( i < VectorFieldSize ? v[ i ] : 0.0 ) << " ";
-      str << std::endl;
-   }
-
-   return true;
-}
-
-/****
- * 3D grid, vertices
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-void
-VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 0, Real > > >::
-writeHeader( const VectorFieldType& vectorField,
-             std::ostream& str )
-{
-    const MeshType& mesh = vectorField.getMesh();
-    const typename MeshType::PointType& origin = mesh.getOrigin();
-    const typename MeshType::PointType& proportions = mesh.getProportions();
-    str << "# vtk DataFile Version 2.0" << std::endl;
-    str << "TNL DATA" << std::endl;
-    str << "ASCII" << std::endl;
-    str << "DATASET UNSTRUCTURED_GRID" << std::endl;
-}
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          int VectorFieldSize >
-bool
-VectorFieldVTKWriter< VectorField< VectorFieldSize, MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 0, Real > > >::
-write( const VectorFieldType& vectorField,
-       std::ostream& str,
-       const double& scale )
-{
-   writeHeader(vectorField, str);
- 
-   const MeshType& mesh = vectorField.getMesh();
-   const RealType originX = mesh.getOrigin().x();
-   const RealType spaceStepX = mesh.getSpaceSteps().x();
-   const RealType originY = mesh.getOrigin().y();
-   const RealType spaceStepY = mesh.getSpaceSteps().y();
-   const RealType originZ = mesh.getOrigin().z();
-   const RealType spaceStepZ = mesh.getSpaceSteps().z();
-   const MeshIndex verticesCount = mesh.template getEntitiesCount< typename MeshType::Vertex >();
- 
-   str << "POINTS " << verticesCount << " float" << std::endl;
-   for (int k = 0; k <= mesh.getDimensions().y(); k++)
-   {
-       for (int j = 0; j <= mesh.getDimensions().y(); j++)
-       {
-            for (int i = 0; i <= mesh.getDimensions().x(); i++)
-            {
-                 str << originX + i * spaceStepX << " " << originY + j * spaceStepY << " " <<
-                        originZ + k * spaceStepZ << std::endl;
-            }
-       }
-   }
- 
-   str << std::endl << "CELLS " << verticesCount << " " << verticesCount * 2 << std::endl;
-   for (int k = 0; k < ( mesh.getDimensions().z() + 1 ); k++)
-   {
-        for (int j = 0; j < ( mesh.getDimensions().y() + 1 ); j++)
-        {
-            for (int i = 0; i < ( mesh.getDimensions().x() + 1 ); i++)
-            {
-                str << "1 " << k * ( mesh.getDimensions().y() + 1 ) * ( mesh.getDimensions().x() + 1 ) + j * ( mesh.getDimensions().x() + 1 ) + i  << std::endl;
-            }
-        }
-   }
- 
-   str << std::endl << "CELL_TYPES " << verticesCount << std::endl;
-   for (int i = 0; i < verticesCount; i++)
-   {
-       str << "1" << std::endl;
-   }
- 
-   str << std::endl << "CELL_DATA " << verticesCount << std::endl;
-   str << "VECTORS verticesVectorFieldValues " << getType< typename VectorFieldType::RealType >() << std::endl;
-
-   for( MeshIndex i = 0; i < verticesCount; i++ )
-   {
-      typename MeshType::Vertex entity = mesh.template getEntity< typename MeshType::Vertex >( i );
-      entity.refresh();
-      const VectorType v = vectorField.getElement( entity.getIndex() );
-      for( int i = 0; i < 3; i++ )
-         str << scale * ( i < VectorFieldSize ? v[ i ] : 0.0 ) << " ";
-      str << std::endl;
-   }
-
-   return true;
-}
-
-} // namespace Functions
-} // namespace TNL