From 78393f314863312d3e962febfd96f51f38c415d5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jakub=20Klinkovsk=C3=BD?= <klinkjak@fjfi.cvut.cz>
Date: Sat, 4 Mar 2017 22:41:19 +0100
Subject: [PATCH] Refactoring mesh function writers

---
 src/TNL/Functions/CMakeLists.txt              |   1 -
 src/TNL/Functions/MeshFunctionVTKWriter.h     | 258 +----
 .../Functions/MeshFunctionVTKWriter_impl.h    | 976 ------------------
 src/TNL/Meshes/CMakeLists.txt                 |   1 +
 src/TNL/Meshes/GridDetails/GridEntity_impl.h  |  13 +
 src/TNL/Meshes/GridEntity.h                   |   4 +
 src/TNL/Meshes/Writers/CMakeLists.txt         |   5 +
 src/TNL/Meshes/Writers/VTKWriter.h            |  56 +
 src/TNL/Meshes/Writers/VTKWriter_impl.h       | 475 +++++++++
 9 files changed, 588 insertions(+), 1201 deletions(-)
 delete mode 100644 src/TNL/Functions/MeshFunctionVTKWriter_impl.h
 create mode 100755 src/TNL/Meshes/Writers/CMakeLists.txt
 create mode 100644 src/TNL/Meshes/Writers/VTKWriter.h
 create mode 100644 src/TNL/Meshes/Writers/VTKWriter_impl.h

diff --git a/src/TNL/Functions/CMakeLists.txt b/src/TNL/Functions/CMakeLists.txt
index 75a34d3181..806f716635 100644
--- a/src/TNL/Functions/CMakeLists.txt
+++ b/src/TNL/Functions/CMakeLists.txt
@@ -11,7 +11,6 @@ SET( headers Domain.h
              MeshFunctionGnuplotWriter_impl.h
              MeshFunctionNormGetter.h
              MeshFunctionVTKWriter.h
-             MeshFunctionVTKWriter_impl.h             
              OperatorFunction.h
              TestFunction.h             
              TestFunction_impl.h
diff --git a/src/TNL/Functions/MeshFunctionVTKWriter.h b/src/TNL/Functions/MeshFunctionVTKWriter.h
index a910cc6d60..674202ea1f 100644
--- a/src/TNL/Functions/MeshFunctionVTKWriter.h
+++ b/src/TNL/Functions/MeshFunctionVTKWriter.h
@@ -2,7 +2,7 @@
                           MeshFunctionVTKWriter.h  -  description
                              -------------------
     begin                : Jan 28, 2016
-    copyright            : (C) 2016 by oberhuber
+    copyright            : (C) 2016 by Tomas Oberhuber et al.
     email                : tomas.oberhuber@fjfi.cvut.cz
  ***************************************************************************/
 
@@ -10,237 +10,47 @@
 
 #pragma once
 
-#include <TNL/Meshes/Grid.h>
+#include <TNL/Meshes/Writers/VTKWriter.h>
 
 namespace TNL {
 namespace Functions {   
 
-template< typename, int, typename > class MeshFunction;
-
 template< typename MeshFunction >
 class MeshFunctionVTKWriter
 {
-   public:
- 
-      static bool write( const MeshFunction& function,
-                         std::ostream& str );
-      static void writeHeader( const MeshFunction& function,
-                               std::ostream& str );
-};
-
-template< typename Mesh,
-          typename Real >
-class MeshFunctionVTKWriter< MeshFunction< Mesh, 0, Real > >
-{
-   using MeshFunctionType = MeshFunction< Mesh, 0, Real >;
-
-   public:
- 
-      static bool write( const MeshFunctionType& function,
-                         std::ostream& str );
-      static void writeHeader( const MeshFunctionType& function,
-                               std::ostream& str );
-};
-
-/***
- * 1D grid, cells
- */
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real >
-class MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 1, Real > >
-{
-   public:
-      typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::MeshFunction< MeshType, 1, RealType > MeshFunctionType;
-
-      static bool write( const MeshFunctionType& function,
-                         std::ostream& str,
-                         const double& scale );
-      
-      static void writeHeader(const MeshFunctionType& function,
-                         std::ostream& str );
-};
- 
-/***
- * 1D grid, vertices
- */
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real >
-class MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 0, Real > >
-{
-   public:
-      typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::MeshFunction< MeshType, 0, RealType > MeshFunctionType;
-
-      static bool write( const MeshFunctionType& function,
-                         std::ostream& str,
-                         const double& scale );
-      
-      static void writeHeader(const MeshFunctionType& function,
-                         std::ostream& str );
-};
-
-/***
- * 2D grid, cells
- */
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real >
-class MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 2, Real > >
-{
-   public:
-      typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::MeshFunction< MeshType, 2, RealType > MeshFunctionType;
-
-      static bool write( const MeshFunctionType& function,
-                         std::ostream& str,
-                         const double& scale );
-      
-      static void writeHeader(const MeshFunctionType& function,
-                         std::ostream& str );
-};
-
-/***
- * 2D grid, faces
- */
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real >
-class MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 1, Real > >
-{
-   public:
-      typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::MeshFunction< MeshType, 1, RealType > MeshFunctionType;
-
-      static bool write( const MeshFunctionType& function,
-                         std::ostream& str,
-                         const double& scale );
-      
-      static void writeHeader(const MeshFunctionType& function,
-                         std::ostream& str );
-};
-
-/***
- * 2D grid, vertices
- */
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real >
-class MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 0, Real > >
-{
-   public:
-      typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::MeshFunction< MeshType, 0, RealType > MeshFunctionType;
-
-      static bool write( const MeshFunctionType& function,
-                         std::ostream& str,
-                         const double& scale );
-      
-      static void writeHeader(const MeshFunctionType& function,
-                         std::ostream& str );
-};
-
-/***
- * 3D grid, cells
- */
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real >
-class MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 3, Real > >
-{
-   public:
-      typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::MeshFunction< MeshType, 3, RealType > MeshFunctionType;
-
-      static bool write( const MeshFunctionType& function,
-                         std::ostream& str,
-                         const double& scale );
-      
-      static void writeHeader(const MeshFunctionType& function,
-                         std::ostream& str );
-};
-
-/***
- * 3D grid, faces
- */
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real >
-class MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 2, Real > >
-{
-   public:
-      typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::MeshFunction< MeshType, 2, RealType > MeshFunctionType;
-
-      static bool write( const MeshFunctionType& function,
-                         std::ostream& str,
-                         const double& scale );
-      
-      static void writeHeader(const MeshFunctionType& function,
-                         std::ostream& str );
-};
-
-/***
- * 3D grid, edges
- */
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real >
-class MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 1, Real > >
-{
-   public:
-      typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::MeshFunction< MeshType, 1, RealType > MeshFunctionType;
-
-      static bool write( const MeshFunctionType& function,
-                         std::ostream& str,
-                         const double& scale );
-      
-      static void writeHeader(const MeshFunctionType& function,
-                         std::ostream& str );
-};
-
-/***
- * 3D grid, vertices
- */
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real >
-class MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 0, Real > >
-{
-   public:
-      typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Functions::MeshFunction< MeshType, 0, RealType > MeshFunctionType;
-
-      static bool write( const MeshFunctionType& function,
-                         std::ostream& str,
-                         const double& scale );
-      
-      static void writeHeader(const MeshFunctionType& function,
-                         std::ostream& str );
+   using MeshType = typename MeshFunction::MeshType;
+   using MeshWriter = Meshes::Writers::VTKWriter< MeshType >;
+   using EntityType = typename MeshType::template EntityType< MeshFunction::getEntitiesDimension() >;
+   using GlobalIndex = typename MeshType::GlobalIndexType;
+
+public:
+   static bool write( const MeshFunction& function,
+                      std::ostream& str,
+                      const String& functionName = "cellFunctionValues" )
+   {
+      const MeshType& mesh = function.getMesh();
+      MeshWriter::template writeEntities< MeshFunction::getEntitiesDimension() >( mesh, str );
+      appendFunction( function, str, functionName );
+      return true;
+   }
+
+   // VTK supports writing multiple functions into the same file.
+   // You can call this after 'write', which initializes the mesh entities,
+   // with different function name.
+   static void appendFunction( const MeshFunction& function,
+                               std::ostream& str,
+                               const String& functionName )
+   {
+      const MeshType& mesh = function.getMesh();
+      const GlobalIndex entitiesCount = mesh.template getEntitiesCount< EntityType >();
+      str << std::endl << "CELL_DATA " << entitiesCount << std::endl;
+      str << "SCALARS " << functionName << " " << getType< typename MeshFunction::RealType >() << " 1" << std::endl;
+      str << "LOOKUP_TABLE default" << std::endl;
+      for( GlobalIndex i = 0; i < entitiesCount; i++ ) {
+         str << function.getData().getElement( i ) << std::endl;
+      }
+   }
 };
 
 } // namespace Functions
 } // namespace TNL
-
-#include <TNL/Functions/MeshFunctionVTKWriter_impl.h>
diff --git a/src/TNL/Functions/MeshFunctionVTKWriter_impl.h b/src/TNL/Functions/MeshFunctionVTKWriter_impl.h
deleted file mode 100644
index 2c8777c76c..0000000000
--- a/src/TNL/Functions/MeshFunctionVTKWriter_impl.h
+++ /dev/null
@@ -1,976 +0,0 @@
-/***************************************************************************
-                          MeshFunctionVTKWriter_impl.h  -  description
-                             -------------------
-    begin                : Jan 28, 2016
-    copyright            : (C) 2016 by oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/* See Copyright Notice in tnl/Copyright */
-
-#pragma once
-
-#include <limits>
-
-#include <TNL/Functions/MeshFunctionVTKWriter.h>
-#include <TNL/Meshes/Readers/VTKEntityType.h>
-
-namespace TNL {
-namespace Functions {   
-
-template< typename MeshFunction >
-void
-MeshFunctionVTKWriter< MeshFunction >::
-writeHeader( const MeshFunction& function,
-             std::ostream& str )
-{
-    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 MeshFunction >
-bool
-MeshFunctionVTKWriter< MeshFunction >::
-write( const MeshFunction& function,
-       std::ostream& str )
-{
-   writeHeader( function, str );
-
-   using MeshType = typename MeshFunction::MeshType;
-   using EntityType = typename MeshType::template EntityType< MeshFunction::getEntitiesDimensions() >;
-   using GlobalIndex = typename MeshType::GlobalIndexType;
-   using LocalIndex = typename MeshType::LocalIndexType;
-
-   const MeshType& mesh = function.getMesh();
-   const GlobalIndex verticesCount = mesh.template getEntitiesCount< 0 >();
-   const GlobalIndex entitiesCount = mesh.template getEntitiesCount< MeshFunction::getEntitiesDimensions() >();
-   const LocalIndex verticesPerEntity = EntityType::getVerticesCount();
-
-   str << "POINTS " << verticesCount << " " << getType< typename MeshType::RealType >() << std::endl;
-   str.precision( std::numeric_limits< typename MeshType::RealType >::digits10 );
-   for( GlobalIndex i = 0; i < verticesCount; i++ ) {
-      const auto& vertex = mesh.template getEntity< 0 >( i );
-      const auto& point = vertex.getPoint();
-      for( LocalIndex j = 0; j < point.size; j++ ) {
-         str << point[ j ];
-         if( j < point.size - 1 )
-            str << " ";
-      }
-      // VTK needs zeros for unused dimensions
-      for( LocalIndex j = 0; j < 3 - point.size; j++ )
-         str << " 0";
-      str << std::endl;
-   }
- 
-   const GlobalIndex cellsListSize = entitiesCount * ( verticesPerEntity + 1 );
-   str << std::endl << "CELLS " << entitiesCount << " " << cellsListSize << std::endl;
-   for( GlobalIndex i = 0; i < entitiesCount; i++ ) {
-      const auto& entity = mesh.template getEntity< MeshFunction::getEntitiesDimensions() >( i );
-      str << verticesPerEntity;
-      for( LocalIndex j = 0; j < verticesPerEntity; j++ )
-         str << " " << entity.template getSubentityIndex< 0 >( j );
-      str << std::endl;
-   }
- 
-   str << std::endl << "CELL_TYPES " << entitiesCount << std::endl;
-   for( GlobalIndex i = 0; i < entitiesCount; i++ ) {
-      const int type = (int) Meshes::Readers::TopologyToVTKMap< typename EntityType::EntityTopology >::type;
-      str << type << std::endl;
-   }
- 
-   str << std::endl << "CELL_DATA " << entitiesCount << std::endl;
-   str << "SCALARS cellFunctionValues float 1" << std::endl;
-   str << "LOOKUP_TABLE default" << std::endl;
-   for( GlobalIndex i = 0; i < entitiesCount; i++ ) {
-      str << function.getData().getElement( i ) << std::endl;
-   }
- 
-   return true;
-}
-
-template< typename Mesh,
-          typename Real >
-void
-MeshFunctionVTKWriter< MeshFunction< Mesh, 0, Real > >::
-writeHeader( const MeshFunctionType& function,
-             std::ostream& str )
-{
-    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 Mesh,
-          typename Real >
-bool
-MeshFunctionVTKWriter< MeshFunction< Mesh, 0, Real > >::
-write( const MeshFunctionType& function,
-       std::ostream& str )
-{
-   writeHeader( function, str );
-
-   using MeshType = typename MeshFunctionType::MeshType;
-   using EntityType = typename MeshType::template EntityType< MeshFunctionType::getEntitiesDimensions() >;
-   using GlobalIndex = typename MeshType::GlobalIndexType;
-   using LocalIndex = typename MeshType::LocalIndexType;
-
-   const MeshType& mesh = function.getMesh();
-   const GlobalIndex verticesCount = mesh.template getEntitiesCount< 0 >();
-   const GlobalIndex entitiesCount = mesh.template getEntitiesCount< MeshFunctionType::getEntitiesDimensions() >();
-   const LocalIndex verticesPerEntity = 1;
-
-   str << "POINTS " << verticesCount << " float" << std::endl;
-   for( GlobalIndex i = 0; i < verticesCount; i++ ) {
-      const auto& vertex = mesh.template getEntity< 0 >( i );
-      const auto& point = vertex.getPoint();
-      for( LocalIndex j = 0; j < point.size; j++ ) {
-         str << point[ j ];
-         if( j < point.size - 1 )
-            str << " ";
-      }
-      // VTK needs zeros for unused dimensions
-      for( LocalIndex j = 0; j < 3 - point.size; j++ )
-         str << " 0";
-      str << std::endl;
-   }
- 
-   const GlobalIndex cellsListSize = entitiesCount * ( verticesPerEntity + 1 );
-   str << std::endl << "CELLS " << entitiesCount << " " << cellsListSize << std::endl;
-   for( GlobalIndex i = 0; i < entitiesCount; i++ ) {
-      const auto& entity = mesh.template getEntity< MeshFunctionType::getEntitiesDimensions() >( i );
-      str << verticesPerEntity << i << std::endl;
-   }
- 
-   str << std::endl << "CELL_TYPES " << entitiesCount << std::endl;
-   for( GlobalIndex i = 0; i < entitiesCount; i++ ) {
-      const int type = (int) Meshes::Readers::TopologyToVTKMap< typename EntityType::EntityTopology >::type;
-      str << type << std::endl;
-   }
- 
-   str << std::endl << "CELL_DATA " << entitiesCount << std::endl;
-   str << "SCALARS cellFunctionValues float 1" << std::endl;
-   str << "LOOKUP_TABLE default" << std::endl;
-   for( GlobalIndex i = 0; i < entitiesCount; i++ ) {
-      str << function.getData().getElement( i ) << std::endl;
-   }
- 
-   return true;
-}
-
-
-/****
- * 1D grid, cells
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real >
-void
-MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 1, Real > >::
-writeHeader( const MeshFunctionType& function,
-             std::ostream& str )
-{
-    const MeshType& mesh = function.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 >
-bool
-MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 1, Real > >::
-write( const MeshFunctionType& function,
-       std::ostream& str,
-       const double& scale )
-{
-   writeHeader(function, str);
- 
-   const MeshType& mesh = function.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 << "SCALARS cellFunctionValues float 1" << std::endl;
-   str << "LOOKUP_TABLE default" << 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();
-      str << scale * function.getData().getElement( entity.getIndex() ) << std::endl;
-   }
- 
-   return true;
-}
-
-/****
- * 1D grid, vertices
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real >
-void
-MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 0, Real > >::
-writeHeader( const MeshFunctionType& function,
-             std::ostream& str )
-{
-    const MeshType& mesh = function.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 >
-bool
-MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 0, Real > >::
-write( const MeshFunctionType& function,
-       std::ostream& str,
-       const double& scale )
-{
-   writeHeader(function, str);
- 
-   const MeshType& mesh = function.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 << "SCALARS VerticesFunctionValues float 1" << std::endl;
-   str << "LOOKUP_TABLE default" << 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();
-      str << scale * function.getData().getElement( entity.getIndex() ) << std::endl;
-   }
- 
-   return true;
-}
-
-/****
- * 2D grid, cells
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real >
-void
-MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 2, Real > >::
-writeHeader( const MeshFunctionType& function,
-             std::ostream& str )
-{
-    const MeshType& mesh = function.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 >
-bool
-MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 2, Real > >::
-write( const MeshFunctionType& function,
-       std::ostream& str,
-       const double& scale )
-{
-   writeHeader(function, str);
- 
-   const MeshType& mesh = function.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 << " 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 * 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 << "SCALARS cellFunctionValues float 1" << std::endl;
-   str << "LOOKUP_TABLE default" << std::endl;
-
-   for( MeshIndex i = 0; i < entitiesCount; i++ )
-   {
-      typename MeshType::Cell entity = mesh.template getEntity< typename MeshType::Cell >( i );
-      entity.refresh();
-      str << scale * function.getData().getElement( entity.getIndex() ) << std::endl;
-   }
-
-   return true;
-}
-
-/****
- * 2D grid, faces
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real >
-void
-MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 1, Real > >::
-writeHeader( const MeshFunctionType& function,
-             std::ostream& str )
-{
-    const MeshType& mesh = function.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 >
-bool
-MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 1, Real > >::
-write( const MeshFunctionType& function,
-       std::ostream& str,
-       const double& scale )
-{
-   typedef typename MeshType::template EntityType< 0 > Vertex;
-   typedef typename MeshType::template EntityType< 1 > Face;
-   writeHeader(function, str);
- 
-   const MeshType& mesh = function.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 << "SCALARS FaceslFunctionValues float 1" << std::endl;
-   str << "LOOKUP_TABLE default" << std::endl;
-
-   for( MeshIndex i = 0; i < entitiesCount; i++ )
-   {
-      typename MeshType::Face entity = mesh.template getEntity< typename MeshType::Face >( i );
-      entity.refresh();
-      str << scale * function.getData().getElement( entity.getIndex() ) << std::endl;
-   }
-
-   return true;
-}
-
-/****
- * 2D grid, vertices
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real >
-void
-MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 0, Real > >::
-writeHeader( const MeshFunctionType& function,
-             std::ostream& str )
-{
-    const MeshType& mesh = function.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 >
-bool
-MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 0, Real > >::
-write( const MeshFunctionType& function,
-       std::ostream& str,
-       const double& scale )
-{
-   typedef typename MeshType::template EntityType< 0 > Vertex;
-   writeHeader(function, str);
- 
-   const MeshType& mesh = function.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 << "SCALARS VerticesFunctionValues float 1" << std::endl;
-   str << "LOOKUP_TABLE default" << std::endl;
-
-   for( MeshIndex i = 0; i < verticesCount; i++ )
-   {
-      typename MeshType::Vertex entity = mesh.template getEntity< typename MeshType::Vertex >( i );
-      entity.refresh();
-      str << scale * function.getData().getElement( entity.getIndex() ) << std::endl;
-   }
-
-   return true;
-}
-
-/****
- * 3D grid, cells
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real >
-void
-MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 3, Real > >::
-writeHeader( const MeshFunctionType& function,
-             std::ostream& str )
-{
-    const MeshType& mesh = function.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 >
-bool
-MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 3, Real > >::
-write( const MeshFunctionType& function,
-       std::ostream& str,
-       const double& scale )
-{
-   writeHeader(function, str);
- 
-   const MeshType& mesh = function.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 << "SCALARS cellFunctionValues float 1" << std::endl;
-   str << "LOOKUP_TABLE default" << std::endl;
-
-   for( MeshIndex i = 0; i < entitiesCount; i++ )
-   {
-      typename MeshType::Cell entity = mesh.template getEntity< typename MeshType::Cell >( i );
-      entity.refresh();
-      str << scale * function.getData().getElement( entity.getIndex() ) << std::endl;
-   }
-
-   return true;
-}
-
-/****
- * 3D grid, faces
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real >
-void
-MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 2, Real > >::
-writeHeader( const MeshFunctionType& function,
-             std::ostream& str )
-{
-    const MeshType& mesh = function.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 >
-bool
-MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 2, Real > >::
-write( const MeshFunctionType& function,
-       std::ostream& str,
-       const double& scale )
-{
-   writeHeader(function, str);
- 
-   const MeshType& mesh = function.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 << "SCALARS facesFunctionValues float 1" << std::endl;
-   str << "LOOKUP_TABLE default" << std::endl;
-
-   for( MeshIndex i = 0; i < entitiesCount; i++ )
-   {
-      typename MeshType::Face entity = mesh.template getEntity< typename MeshType::Face >( i );
-      entity.refresh();
-      str << scale * function.getData().getElement( entity.getIndex() ) << std::endl;
-   }
-
-   return true;
-}
-
-/****
- * 3D grid, edges
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real >
-void
-MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 1, Real > >::
-writeHeader( const MeshFunctionType& function,
-             std::ostream& str )
-{
-    const MeshType& mesh = function.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 >
-bool
-MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 1, Real > >::
-write( const MeshFunctionType& function,
-       std::ostream& str,
-       const double& scale )
-{
-   writeHeader(function, str);
- 
-   const MeshType& mesh = function.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 << "SCALARS edgesFunctionValues float 1" << std::endl;
-   str << "LOOKUP_TABLE default" << std::endl;
-
-   for( MeshIndex i = 0; i < entitiesCount; i++ )
-   {
-      typename MeshType::Edge entity = mesh.template getEntity< typename MeshType::Edge >( i );
-      entity.refresh();
-      str << scale * function.getData().getElement( entity.getIndex() ) << std::endl;
-   }
-
-   return true;
-}
-
-/****
- * 3D grid, vertices
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real >
-void
-MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 0, Real > >::
-writeHeader( const MeshFunctionType& function,
-             std::ostream& str )
-{
-    const MeshType& mesh = function.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 >
-bool
-MeshFunctionVTKWriter< MeshFunction< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 0, Real > >::
-write( const MeshFunctionType& function,
-       std::ostream& str,
-       const double& scale )
-{
-   writeHeader(function, str);
- 
-   const MeshType& mesh = function.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 << "SCALARS verticesFunctionValues float 1" << std::endl;
-   str << "LOOKUP_TABLE default" << std::endl;
-
-   for( MeshIndex i = 0; i < verticesCount; i++ )
-   {
-      typename MeshType::Vertex entity = mesh.template getEntity< typename MeshType::Vertex >( i );
-      entity.refresh();
-      str << scale * function.getData().getElement( entity.getIndex() ) << std::endl;
-   }
-
-   return true;
-}
-
-} // namespace Functions
-} // namespace TNL
-
diff --git a/src/TNL/Meshes/CMakeLists.txt b/src/TNL/Meshes/CMakeLists.txt
index 82746c61d5..eddb9e7f20 100644
--- a/src/TNL/Meshes/CMakeLists.txt
+++ b/src/TNL/Meshes/CMakeLists.txt
@@ -3,6 +3,7 @@ ADD_SUBDIRECTORY( MeshDetails )
 ADD_SUBDIRECTORY( Readers )
 ADD_SUBDIRECTORY( Topologies )
 ADD_SUBDIRECTORY( TypeResolver )
+ADD_SUBDIRECTORY( Writers )
 
 SET( headers BuildConfigTags.h
              DimensionTag.h
diff --git a/src/TNL/Meshes/GridDetails/GridEntity_impl.h b/src/TNL/Meshes/GridDetails/GridEntity_impl.h
index 8703c064fd..7c3820f6b6 100644
--- a/src/TNL/Meshes/GridDetails/GridEntity_impl.h
+++ b/src/TNL/Meshes/GridDetails/GridEntity_impl.h
@@ -675,6 +675,19 @@ getCenter() const
    return GridEntityCenterGetter< ThisType >::getEntityCenter( *this );
 }
 
+template< int Dimension,
+          typename Real,
+          typename Device,
+          typename Index,
+          typename Config >
+__cuda_callable__ inline
+typename Meshes::Grid< Dimension, Real, Device, Index >::PointType
+GridEntity< Meshes::Grid< Dimension, Real, Device, Index >, 0, Config >::
+getPoint() const
+{
+   return getCenter();
+}
+
 template< int Dimension,
           typename Real,
           typename Device,
diff --git a/src/TNL/Meshes/GridEntity.h b/src/TNL/Meshes/GridEntity.h
index 4db5017328..abad0a57d5 100644
--- a/src/TNL/Meshes/GridEntity.h
+++ b/src/TNL/Meshes/GridEntity.h
@@ -362,6 +362,10 @@ class GridEntity< Meshes::Grid< Dimension, Real, Device, Index >, 0, Config >
       __cuda_callable__ inline
       PointType getCenter() const;
 
+      // compatibility with meshes, equivalent to getCenter
+      __cuda_callable__ inline
+      PointType getPoint() const;
+
       __cuda_callable__ inline
       const RealType getMeasure() const;
  
diff --git a/src/TNL/Meshes/Writers/CMakeLists.txt b/src/TNL/Meshes/Writers/CMakeLists.txt
new file mode 100755
index 0000000000..73354ff1a5
--- /dev/null
+++ b/src/TNL/Meshes/Writers/CMakeLists.txt
@@ -0,0 +1,5 @@
+SET( headers VTKWriter.h
+             VTKWriter_impl.h
+)
+
+INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/TNL/Meshes/Writers )
diff --git a/src/TNL/Meshes/Writers/VTKWriter.h b/src/TNL/Meshes/Writers/VTKWriter.h
new file mode 100644
index 0000000000..9b47f4efd7
--- /dev/null
+++ b/src/TNL/Meshes/Writers/VTKWriter.h
@@ -0,0 +1,56 @@
+/***************************************************************************
+                          VTKWriter.h  -  description
+                             -------------------
+    begin                : Mar 04, 2017
+    copyright            : (C) 2017 by Tomas Oberhuber et al.
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#pragma once
+
+#include <limits>
+
+#include <TNL/Meshes/Grid.h>
+#include <TNL/Meshes/Mesh.h>
+
+namespace TNL {
+namespace Meshes {
+namespace Writers {
+
+namespace __impl {
+
+template< typename Mesh, int EntityDimension > struct MeshEntitiesVTKWriter;
+template< typename Mesh, int EntityDimension > struct MeshEntityTypesVTKWriter;
+
+} // namespace __impl
+
+template< typename Mesh >
+class VTKWriter
+{
+   template< int EntityDimension >
+   using EntitiesWriter = __impl::MeshEntitiesVTKWriter< Mesh, EntityDimension >;
+
+   template< int EntityDimension >
+   using EntityTypesWriter = __impl::MeshEntityTypesVTKWriter< Mesh, EntityDimension >;
+
+public:
+   using Index = typename Mesh::GlobalIndexType;
+
+   static void writeAllEntities( const Mesh& mesh, std::ostream& str );
+
+   template< int EntityDimension = Mesh::getMeshDimension() >
+   static void writeEntities( const Mesh& mesh, std::ostream& str );
+
+protected:
+   static void writeHeader( const Mesh& mesh, std::ostream& str );
+
+   static void writePoints( const Mesh& mesh, std::ostream& str );
+};
+
+} // namespace Writers
+} // namespace Meshes
+} // namespace TNL
+
+#include <TNL/Meshes/Writers/VTKWriter_impl.h>
diff --git a/src/TNL/Meshes/Writers/VTKWriter_impl.h b/src/TNL/Meshes/Writers/VTKWriter_impl.h
new file mode 100644
index 0000000000..119d305715
--- /dev/null
+++ b/src/TNL/Meshes/Writers/VTKWriter_impl.h
@@ -0,0 +1,475 @@
+/***************************************************************************
+                          VTKWriter.h  -  description
+                             -------------------
+    begin                : Mar 04, 2017
+    copyright            : (C) 2017 by Tomas Oberhuber et al.
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#pragma once
+
+#include <TNL/Meshes/Writers/VTKWriter.h>
+#include <TNL/Meshes/Readers/VTKEntityType.h>
+
+namespace TNL {
+namespace Meshes {
+namespace Writers {
+
+namespace __impl {
+
+template< typename Entity >
+struct VerticesPerEntity
+{
+   static constexpr int count = Entity::getVerticesCount();
+};
+
+template< typename MeshConfig, typename Device >
+struct VerticesPerEntity< MeshEntity< MeshConfig, Device, MeshVertexTopology > >
+{
+   static constexpr int count = 1;
+};
+
+template< typename Grid, typename Config >
+struct VerticesPerEntity< GridEntity< Grid, 0, Config > >
+{
+   static constexpr int count = 1;
+};
+
+template< typename Grid, typename Config >
+struct VerticesPerEntity< GridEntity< Grid, 1, Config > >
+{
+   static constexpr int count = 2;
+};
+
+template< typename Grid, typename Config >
+struct VerticesPerEntity< GridEntity< Grid, 2, Config > >
+{
+   static constexpr int count = 4;
+};
+
+template< typename Grid, typename Config >
+struct VerticesPerEntity< GridEntity< Grid, 3, Config > >
+{
+   static constexpr int count = 8;
+};
+
+
+template< typename GridEntity >
+struct GridEntityToVTKType {};
+
+template< typename Grid, typename Config >
+struct GridEntityToVTKType< GridEntity< Grid, 0, Config > >
+{
+   static constexpr Readers::VTKEntityType type = Readers::VTKEntityType::Vertex;
+};
+
+template< typename Grid, typename Config >
+struct GridEntityToVTKType< GridEntity< Grid, 1, Config > >
+{
+   static constexpr Readers::VTKEntityType type = Readers::VTKEntityType::Line;
+};
+
+template< typename Grid, typename Config >
+struct GridEntityToVTKType< GridEntity< Grid, 2, Config > >
+{
+   static constexpr Readers::VTKEntityType type = Readers::VTKEntityType::Pixel;
+};
+
+template< typename Grid, typename Config >
+struct GridEntityToVTKType< GridEntity< Grid, 3, Config > >
+{
+   static constexpr Readers::VTKEntityType type = Readers::VTKEntityType::Voxel;
+};
+
+
+template< typename Mesh >
+typename Mesh::GlobalIndexType
+getAllMeshEntitiesCount( const Mesh& mesh, DimensionTag< 0 > )
+{
+   using EntityType = typename Mesh::template EntityType< 0 >;
+   return mesh.template getEntitiesCount< EntityType >();
+}
+
+// TODO: specialization for disabled entities
+template< typename Mesh,
+          typename DimensionTag = Meshes::DimensionTag< Mesh::getMeshDimension() > >
+typename Mesh::GlobalIndexType
+getAllMeshEntitiesCount( const Mesh& mesh, DimensionTag = DimensionTag() )
+{
+   using EntityType = typename Mesh::template EntityType< DimensionTag::value >;
+   return mesh.template getEntitiesCount< EntityType >() +
+          getAllMeshEntitiesCount( mesh, typename DimensionTag::Decrement() );
+}
+
+
+template< typename Mesh >
+typename Mesh::GlobalIndexType
+getCellsListSize( const Mesh& mesh, DimensionTag< 0 > )
+{
+   using EntityType = typename Mesh::template EntityType< 0 >;
+   return mesh.template getEntitiesCount< EntityType >() * 2;
+}
+
+// TODO: specialization for disabled entities
+template< typename Mesh,
+          typename DimensionTag = Meshes::DimensionTag< Mesh::getMeshDimension() > >
+typename Mesh::GlobalIndexType
+getCellsListSize( const Mesh& mesh, DimensionTag = DimensionTag() )
+{
+   using EntityType = typename Mesh::template EntityType< DimensionTag::value >;
+   const auto verticesPerEntity = VerticesPerEntity< EntityType >::count;
+   return ( mesh.template getEntitiesCount< EntityType >() * ( verticesPerEntity + 1 ) ) +
+          getCellsListSize( mesh, typename DimensionTag::Decrement() );
+}
+
+
+// TODO: specialization for disabled entities
+// Unstructured meshes, entities
+template< typename Mesh, int EntityDimension >
+struct MeshEntitiesVTKWriter
+{
+   static void exec( const Mesh& mesh, std::ostream& str )
+   {
+      using EntityType = typename Mesh::template EntityType< EntityDimension >;
+      using Index = typename Mesh::GlobalIndexType;
+
+      const Index entitiesCount = mesh.template getEntitiesCount< EntityType >();
+      const Index verticesPerEntity = VerticesPerEntity< EntityType >::count;;
+      for( Index i = 0; i < entitiesCount; i++ ) {
+         const auto& entity = mesh.template getEntity< EntityType >( i );
+         str << verticesPerEntity;
+         for( Index j = 0; j < verticesPerEntity; j++ )
+            str << " " << entity.template getSubentityIndex< 0 >( j );
+         str << "\n";
+      }
+   }
+};
+
+// Unstructured meshes, vertices
+template< typename Mesh >
+struct MeshEntitiesVTKWriter< Mesh, 0 >
+{
+   static void exec( const Mesh& mesh, std::ostream& str )
+   {
+      using EntityType = typename Mesh::template EntityType< 0 >;
+      using Index = typename Mesh::GlobalIndexType;
+
+      const Index entitiesCount = mesh.template getEntitiesCount< EntityType >();
+      const Index verticesPerEntity = 1;
+      for( Index i = 0; i < entitiesCount; i++ ) {
+         str << verticesPerEntity << " " << i << "\n";
+      }
+   }
+};
+
+// 1D grids, cells
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex >
+struct MeshEntitiesVTKWriter< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 1 >
+{
+   using MeshType = Meshes::Grid< 1, MeshReal, Device, MeshIndex >;
+
+   static void exec( const MeshType& mesh, std::ostream& str )
+   {
+      for( MeshIndex i = 0; i < mesh.getDimensions().x(); i++ )
+         str << "2 " << i << " " << i+1 << "\n";
+   }
+};
+
+// 1D grids, vertices
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex >
+struct MeshEntitiesVTKWriter< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 0 >
+{
+   using MeshType = Meshes::Grid< 1, MeshReal, Device, MeshIndex >;
+
+   static void exec( const MeshType& mesh, std::ostream& str )
+   {
+      for( MeshIndex i = 0; i < mesh.getDimensions().x() + 1; i++ )
+         str << "1 " << i << "\n";
+   }
+};
+
+// 2D grids, cells
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex >
+struct MeshEntitiesVTKWriter< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 2 >
+{
+   using MeshType = Meshes::Grid< 2, MeshReal, Device, MeshIndex >;
+
+   static void exec( const MeshType& mesh, std::ostream& str )
+   {
+      for( MeshIndex j = 0; j < mesh.getDimensions().y(); j++ )
+         for( MeshIndex 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 << "\n";
+   }
+};
+
+// 2D grids, faces
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex >
+struct MeshEntitiesVTKWriter< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 1 >
+{
+   using MeshType = Meshes::Grid< 2, MeshReal, Device, MeshIndex >;
+
+   static void exec( const MeshType& mesh, std::ostream& str )
+   {
+      for( MeshIndex j = 0; j < mesh.getDimensions().y(); j++ )
+         for( MeshIndex i = 0; i < ( mesh.getDimensions().x() + 1 ); i++ )
+            str << "2 " << j * ( mesh.getDimensions().x() + 1 ) + i << " "
+                        << (j+1) * ( mesh.getDimensions().x() + 1 ) + i << "\n";
+
+      for( MeshIndex j = 0; j < (mesh.getDimensions().y()+1); j++ )
+         for( MeshIndex i = 0; i < mesh.getDimensions().x(); i++ )
+            str << "2 " << j * ( mesh.getDimensions().x() + 1 ) + i << " "
+                        << j * ( mesh.getDimensions().x() + 1 ) + i + 1 << "\n";
+   }
+};
+
+// 2D grids, vertices
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex >
+struct MeshEntitiesVTKWriter< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 0 >
+{
+   using MeshType = Meshes::Grid< 2, MeshReal, Device, MeshIndex >;
+
+   static void exec( const MeshType& mesh, std::ostream& str )
+   {
+      for( MeshIndex j = 0; j < ( mesh.getDimensions().y() + 1 ); j++ )
+         for( MeshIndex i = 0; i < ( mesh.getDimensions().x() + 1 ); i++ )
+            str << "1 " << j * mesh.getDimensions().x() + i << "\n";
+   }
+};
+
+// 3D grids, cells
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex >
+struct MeshEntitiesVTKWriter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 3 >
+{
+   using MeshType = Meshes::Grid< 3, MeshReal, Device, MeshIndex >;
+
+   static void exec( const MeshType& mesh, std::ostream& str )
+   {
+      for( MeshIndex k = 0; k < mesh.getDimensions().z(); k++ )
+         for( MeshIndex j = 0; j < mesh.getDimensions().y(); j++ )
+            for( MeshIndex 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 << "\n";
+   }
+};
+
+// 3D grids, faces
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex >
+struct MeshEntitiesVTKWriter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 2 >
+{
+   using MeshType = Meshes::Grid< 3, MeshReal, Device, MeshIndex >;
+
+   static void exec( const MeshType& mesh, std::ostream& str )
+   {
+      for( MeshIndex k = 0; k < mesh.getDimensions().z(); k++ )
+         for( MeshIndex j = 0; j < mesh.getDimensions().y(); j++ )
+            for( MeshIndex 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 << "\n";
+
+      for( MeshIndex k = 0; k < mesh.getDimensions().z(); k++ )
+         for( MeshIndex j = 0; j <= mesh.getDimensions().y(); j++ )
+            for( MeshIndex 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 << "\n";
+
+      for( MeshIndex k = 0; k <= mesh.getDimensions().z(); k++ )
+         for( MeshIndex j = 0; j < mesh.getDimensions().y(); j++ )
+            for( MeshIndex 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 << "\n";
+   }
+};
+
+// 3D grids, edges
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex >
+struct MeshEntitiesVTKWriter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 1 >
+{
+   using MeshType = Meshes::Grid< 3, MeshReal, Device, MeshIndex >;
+
+   static void exec( const MeshType& mesh, std::ostream& str )
+   {
+      for( MeshIndex k = 0; k <= mesh.getDimensions().z(); k++ )
+         for( MeshIndex j = 0; j <= mesh.getDimensions().y(); j++ )
+            for( MeshIndex i = 0; i < mesh.getDimensions().x(); i++ )
+               str << "2 " << 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 << "\n";
+
+      for( MeshIndex k = 0; k <= mesh.getDimensions().z(); k++ )
+         for( MeshIndex j = 0; j < mesh.getDimensions().y(); j++ )
+            for( MeshIndex i = 0; i <= mesh.getDimensions().x(); i++ )
+               str << "2 " << 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 << "\n";
+
+      for( MeshIndex k = 0; k < mesh.getDimensions().z(); k++ )
+         for( MeshIndex j = 0; j <= mesh.getDimensions().y(); j++ )
+            for( MeshIndex i = 0; i <= mesh.getDimensions().x(); i++ )
+               str << "2 " << 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 << "\n";
+   }
+};
+
+// 3D grids, vertices
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex >
+struct MeshEntitiesVTKWriter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 0 >
+{
+   using MeshType = Meshes::Grid< 3, MeshReal, Device, MeshIndex >;
+
+   static void exec( const MeshType& mesh, std::ostream& str )
+   {
+      for( MeshIndex k = 0; k < ( mesh.getDimensions().z() + 1 ); k++ )
+         for( MeshIndex j = 0; j < ( mesh.getDimensions().y() + 1 ); j++ )
+            for( MeshIndex 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  << "\n";
+   }
+};
+
+
+// TODO: specialization for disabled entities
+template< typename Mesh, int EntityDimension >
+struct MeshEntityTypesVTKWriter
+{
+   static void exec( const Mesh& mesh, std::ostream& str )
+   {
+      using EntityType = typename Mesh::template EntityType< EntityDimension >;
+      using Index = typename Mesh::GlobalIndexType;
+
+      const Index entitiesCount = mesh.template getEntitiesCount< EntityType >();
+      for( Index i = 0; i < entitiesCount; i++ ) {
+         const int type = (int) Meshes::Readers::TopologyToVTKMap< typename EntityType::EntityTopology >::type;
+         str << type << "\n";
+      }
+   }
+};
+
+template< int Dimension,
+          typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          int EntityDimension >
+struct MeshEntityTypesVTKWriter< Grid< Dimension, MeshReal, Device, MeshIndex >, EntityDimension >
+{
+   using MeshType = Grid< Dimension, MeshReal, Device, MeshIndex >;
+
+   static void exec( const MeshType& mesh, std::ostream& str )
+   {
+      using EntityType = typename MeshType::template EntityType< EntityDimension >;
+
+      const MeshIndex entitiesCount = mesh.template getEntitiesCount< EntityType >();
+      for( MeshIndex i = 0; i < entitiesCount; i++ ) {
+         const int type = (int) __impl::GridEntityToVTKType< EntityType >::type;
+         str << type << "\n";
+      }
+   }
+};
+
+} // namespace __impl
+
+template< typename Mesh >
+void
+VTKWriter< Mesh >::writeAllEntities( const Mesh& mesh, std::ostream& str )
+{
+   writeHeader( mesh, str );
+   writePoints( mesh, str );
+
+   const Index allEntitiesCount = __impl::getAllMeshEntitiesCount( mesh );
+   const Index cellsListSize = __impl::getCellsListSize( mesh );
+
+   str << std::endl << "CELLS " << allEntitiesCount << " " << cellsListSize << std::endl;
+   StaticFor< int, 0, Mesh::getMeshDimension() + 1, EntitiesWriter >::exec( mesh, str );
+
+   str << std::endl << "CELL_TYPES " << allEntitiesCount << std::endl;
+   StaticFor< int, 0, Mesh::getMeshDimension() + 1, EntityTypesWriter >::exec( mesh, str );
+}
+
+template< typename Mesh >
+   template< int EntityDimension >
+void
+VTKWriter< Mesh >::writeEntities( const Mesh& mesh, std::ostream& str )
+{
+   writeHeader( mesh, str );
+   writePoints( mesh, str );
+
+   using EntityType = typename Mesh::template EntityType< EntityDimension >;
+   const Index entitiesCount = mesh.template getEntitiesCount< EntityType >();
+   const Index verticesPerEntity = __impl::VerticesPerEntity< EntityType >::count;
+   const Index cellsListSize = entitiesCount * ( verticesPerEntity + 1 );
+
+   str << std::endl << "CELLS " << entitiesCount << " " << cellsListSize << std::endl;
+   EntitiesWriter< EntityDimension >::exec( mesh, str );
+
+   str << std::endl << "CELL_TYPES " << entitiesCount << std::endl;
+   EntityTypesWriter< EntityDimension >::exec( mesh, str );
+}
+
+template< typename Mesh >
+void
+VTKWriter< Mesh >::writeHeader( const Mesh& mesh, std::ostream& str )
+{
+    str << "# vtk DataFile Version 2.0\n"
+        << "TNL DATA\n"
+        << "ASCII\n"
+        << "DATASET UNSTRUCTURED_GRID\n";
+}
+
+template< typename Mesh >
+void
+VTKWriter< Mesh >::writePoints( const Mesh& mesh, std::ostream& str )
+{
+   const Index verticesCount = mesh.template getEntitiesCount< typename Mesh::Vertex >();
+
+   str << "POINTS " << verticesCount << " " << getType< typename Mesh::RealType >() << std::endl;
+   str.precision( std::numeric_limits< typename Mesh::RealType >::digits10 );
+
+   for( Index i = 0; i < verticesCount; i++ ) {
+      const auto& vertex = mesh.template getEntity< typename Mesh::Vertex >( i );
+      const auto& point = vertex.getPoint();
+      for( Index j = 0; j < point.size; j++ ) {
+         str << point[ j ];
+         if( j < point.size - 1 )
+            str << " ";
+      }
+      // VTK needs zeros for unused dimensions
+      for( Index j = 0; j < 3 - point.size; j++ )
+         str << " 0";
+      str << "\n";
+   }
+}
+
+} // namespace Writers
+} // namespace Meshes
+} // namespace TNL
-- 
GitLab