diff --git a/src/TNL/Meshes/MeshOrdering.h b/src/TNL/Meshes/MeshOrdering.h
new file mode 100644
index 0000000000000000000000000000000000000000..e3646ea439a3050d45ad2f72513640eb8042f890
--- /dev/null
+++ b/src/TNL/Meshes/MeshOrdering.h
@@ -0,0 +1,50 @@
+// Copyright (c) 2004-2022 Tomáš Oberhuber et al.
+//
+// This file is part of TNL - Template Numerical Library (https://tnl-project.org/)
+//
+// SPDX-License-Identifier: MIT
+
+// Implemented by: Jakub KlinkovskĂ˝
+
+#pragma once
+
+#include <TNL/Meshes/Mesh.h>
+
+namespace TNL {
+namespace Meshes {
+
+// general implementation covering grids
+template< typename Mesh, typename Ordering >
+struct MeshOrdering
+{
+   void
+   reorder( Mesh& mesh )
+   {}
+};
+
+// reordering makes sense only for unstructured meshes
+template< typename MeshConfig, typename Device, typename Ordering >
+struct MeshOrdering< TNL::Meshes::Mesh< MeshConfig, Device >, Ordering >
+{
+   using Mesh = TNL::Meshes::Mesh< MeshConfig, Device >;
+
+   void
+   reorder( Mesh& mesh )
+   {
+      Algorithms::staticFor< int, 0, Mesh::getMeshDimension() >(
+         [ & ]( auto i )
+         {
+            // make sure to reorder cells first
+            constexpr int dim = Mesh::getMeshDimension() - i;
+            using PermutationArray = typename Mesh::GlobalIndexArray;
+            PermutationArray perm;
+            PermutationArray iperm;
+            using EntityType = typename Mesh::template EntityType< dim >;
+            Ordering::template getPermutations< EntityType >( mesh, perm, iperm );
+            mesh.template reorderEntities< dim >( perm, iperm );
+         } );
+   }
+};
+
+}  // namespace Meshes
+}  // namespace TNL
diff --git a/src/TNL/Meshes/MeshOrdering/CuthillMcKeeOrdering.h b/src/TNL/Meshes/MeshOrdering/CuthillMcKeeOrdering.h
new file mode 100644
index 0000000000000000000000000000000000000000..5bb757ba01cc33fbc83f28b5227b9a236f3a9ab7
--- /dev/null
+++ b/src/TNL/Meshes/MeshOrdering/CuthillMcKeeOrdering.h
@@ -0,0 +1,163 @@
+// Copyright (c) 2004-2022 Tomáš Oberhuber et al.
+//
+// This file is part of TNL - Template Numerical Library (https://tnl-project.org/)
+//
+// SPDX-License-Identifier: MIT
+
+// Implemented by: Jakub KlinkovskĂ˝
+
+#pragma once
+
+#include <type_traits>
+
+#include <TNL/Containers/Array.h>
+
+namespace TNL {
+namespace Meshes {
+
+template< bool reverse = true >
+struct CuthillMcKeeOrdering
+{
+   template< typename MeshEntity, typename Mesh, typename PermutationArray >
+   static void
+   getPermutations( const Mesh& mesh, PermutationArray& perm, PermutationArray& iperm )
+   {
+      static_assert( std::is_same< typename Mesh::DeviceType, TNL::Devices::Host >::value, "" );
+      static_assert( std::is_same< typename PermutationArray::DeviceType, TNL::Devices::Host >::value, "" );
+
+      // The reverse Cuthill-McKee ordering is implemented only for cells,
+      // other entities are ordered from the current order of cells exactly
+      // as if the mesh was re-initialized from the cell seeds.
+      Wrapper< MeshEntity, Mesh >::getPermutations( mesh, perm, iperm );
+   }
+
+private:
+   template< typename MeshEntity, typename Mesh, bool is_cell = MeshEntity::getEntityDimension() == Mesh::getMeshDimension() >
+   struct Wrapper
+   {
+      template< typename PermutationArray >
+      static void
+      getPermutations( const Mesh& mesh, PermutationArray& perm, PermutationArray& iperm )
+      {
+         reorderCells( mesh, perm, iperm );
+      }
+   };
+
+   template< typename MeshEntity, typename Mesh >
+   struct Wrapper< MeshEntity, Mesh, false >
+   {
+      template< typename PermutationArray >
+      static void
+      getPermutations( const Mesh& mesh, PermutationArray& perm, PermutationArray& iperm )
+      {
+         reorderEntities< MeshEntity >( mesh, perm, iperm );
+      }
+   };
+
+   template< typename Mesh, typename PermutationArray >
+   static void
+   reorderCells( const Mesh& mesh, PermutationArray& perm, PermutationArray& iperm )
+   {
+      using IndexType = typename Mesh::GlobalIndexType;
+      const IndexType numberOfCells = mesh.template getEntitiesCount< typename Mesh::Cell >();
+
+      // allocate permutation vectors
+      perm.setSize( numberOfCells );
+      iperm.setSize( numberOfCells );
+
+      // vector view for the neighbor counts
+      const auto neighborCounts = mesh.getNeighborCounts().getConstView();
+
+      // worker array - marker for inserted elements
+      TNL::Containers::Array< bool, TNL::Devices::Host, IndexType > marker( numberOfCells );
+      marker.setValue( false );
+      // worker vector for collecting neighbors
+      std::vector< IndexType > neighbors;
+
+      // comparator functor
+      auto comparator = [ & ]( IndexType a, IndexType b )
+      {
+         return neighborCounts[ a ] < neighborCounts[ b ];
+      };
+
+      // counter for assigning indices
+      IndexType permIndex = 0;
+
+      // modifier for the reversed variant
+      auto mod = [ numberOfCells ]( IndexType i )
+      {
+         if( reverse )
+            return numberOfCells - 1 - i;
+         else
+            return i;
+      };
+
+      // start with a peripheral node
+      const IndexType peripheral = argMin( neighborCounts ).first;
+      perm[ mod( permIndex ) ] = peripheral;
+      iperm[ peripheral ] = mod( permIndex );
+      permIndex++;
+      marker[ peripheral ] = true;
+
+      // Cuthill--McKee
+      IndexType i = 0;
+      while( permIndex < numberOfCells ) {
+         const IndexType k = perm[ mod( i ) ];
+         // collect all neighbors which were not marked yet
+         const IndexType count = neighborCounts[ k ];
+         for( IndexType n = 0; n < count; n++ ) {
+            const IndexType nk = mesh.getCellNeighborIndex( k, n );
+            if( marker[ nk ] == false ) {
+               neighbors.push_back( nk );
+               marker[ nk ] = true;
+            }
+         }
+         // sort collected neighbors with ascending neighbors count
+         std::sort( neighbors.begin(), neighbors.end(), comparator );
+         // assign an index to the neighbors in this order
+         for( auto nk : neighbors ) {
+            perm[ mod( permIndex ) ] = nk;
+            iperm[ nk ] = mod( permIndex );
+            permIndex++;
+         }
+         // next iteration
+         i++;
+         neighbors.clear();
+      }
+   }
+
+   template< typename MeshEntity, typename Mesh, typename PermutationArray >
+   static void
+   reorderEntities( const Mesh& mesh, PermutationArray& perm, PermutationArray& iperm )
+   {
+      using IndexType = typename Mesh::GlobalIndexType;
+      const IndexType numberOfEntities = mesh.template getEntitiesCount< MeshEntity >();
+      const IndexType numberOfCells = mesh.template getEntitiesCount< typename Mesh::Cell >();
+
+      // allocate permutation vectors
+      perm.setSize( numberOfEntities );
+      iperm.setSize( numberOfEntities );
+
+      // worker array - marker for numbered entities
+      TNL::Containers::Array< bool, TNL::Devices::Host, IndexType > marker( numberOfEntities );
+      marker.setValue( false );
+
+      IndexType permIndex = 0;
+      for( IndexType K = 0; K < numberOfCells; K++ ) {
+         const auto& cell = mesh.template getEntity< Mesh::getMeshDimension() >( K );
+         for( typename Mesh::LocalIndexType e = 0; e < cell.template getSubentitiesCount< MeshEntity::getEntityDimension() >();
+              e++ ) {
+            const auto E = cell.template getSubentityIndex< MeshEntity::getEntityDimension() >( e );
+            if( marker[ E ] == false ) {
+               marker[ E ] = true;
+               perm[ permIndex ] = E;
+               iperm[ E ] = permIndex;
+               permIndex++;
+            }
+         }
+      }
+   }
+};
+
+}  // namespace Meshes
+}  // namespace TNL
diff --git a/src/TNL/Meshes/MeshOrdering/HilbertOrdering.h b/src/TNL/Meshes/MeshOrdering/HilbertOrdering.h
new file mode 100644
index 0000000000000000000000000000000000000000..1d005bff427e3131cde7f85dec9b504ab5b22965
--- /dev/null
+++ b/src/TNL/Meshes/MeshOrdering/HilbertOrdering.h
@@ -0,0 +1,108 @@
+// Copyright (c) 2004-2022 Tomáš Oberhuber et al.
+//
+// This file is part of TNL - Template Numerical Library (https://tnl-project.org/)
+//
+// SPDX-License-Identifier: MIT
+
+// Implemented by: Jakub KlinkovskĂ˝
+
+#pragma once
+
+#ifdef HAVE_CGAL
+
+   #include <TNL/Meshes/Geometry/getEntityCenter.h>
+
+   // Documentation:
+   // - https://doc.cgal.org/latest/Spatial_sorting/index.html
+   #include <CGAL/hilbert_sort.h>
+
+namespace TNL {
+namespace Meshes {
+
+struct HilbertOrdering
+{
+   template< typename MeshEntity, typename Mesh, typename PermutationArray >
+   static void
+   getPermutations( const Mesh& mesh, PermutationArray& perm, PermutationArray& iperm )
+   {
+      static_assert( std::is_same< typename Mesh::DeviceType, TNL::Devices::Host >::value, "" );
+      static_assert( std::is_same< typename PermutationArray::DeviceType, TNL::Devices::Host >::value, "" );
+      using GlobalIndexType = typename Mesh::GlobalIndexType;
+      using PointType = typename Mesh::PointType;
+      using RealType = typename Mesh::RealType;
+
+      // wrappers for CGAL
+      using pair = std::pair< PointType, GlobalIndexType >;
+      struct Compute_d
+      {
+         RealType
+         operator()( const pair& p, int d ) const
+         {
+            return p.first[ d ];
+         }
+      };
+      struct Less_d
+      {
+         bool
+         operator()( const pair& p, const pair& q, int d ) const
+         {
+            return p.first[ d ] < q.first[ d ];
+         }
+      };
+      struct Point_dim_d
+      {
+         int
+         operator()( const pair& p ) const
+         {
+            return p.first.getSize();
+         }
+      };
+      struct SortingTraits
+      {
+         using Point_d = pair;
+         using Compute_coordinate_d = Compute_d;
+         using Less_coordinate_d = Less_d;
+         using Point_dimension_d = Point_dim_d;
+         Compute_coordinate_d
+         compute_coordinate_d_object() const
+         {
+            return {};
+         }
+         Less_coordinate_d
+         less_coordinate_d_object() const
+         {
+            return {};
+         }
+         Point_dimension_d
+         point_dimension_d_object() const
+         {
+            return {};
+         }
+      };
+
+      // create a vector with entity centers and initial indices
+      std::vector< pair > points;
+      const GlobalIndexType numberOfEntities = mesh.template getEntitiesCount< MeshEntity >();
+      for( GlobalIndexType i = 0; i < numberOfEntities; i++ ) {
+         const auto& entity = mesh.template getEntity< MeshEntity >( i );
+         const auto center = getEntityCenter( mesh, entity );
+         points.push_back( std::make_pair( center, i ) );
+      }
+
+      // sort the points
+      CGAL::hilbert_sort( points.begin(), points.end(), SortingTraits(), CGAL::Hilbert_sort_middle_policy() );
+
+      // build the permutations
+      perm.setSize( numberOfEntities );
+      iperm.setSize( numberOfEntities );
+      for( GlobalIndexType i = 0; i < numberOfEntities; i++ ) {
+         perm[ i ] = points[ i ].second;
+         iperm[ points[ i ].second ] = i;
+      }
+   }
+};
+
+}  // namespace Meshes
+}  // namespace TNL
+
+#endif  // HAVE_CGAL
diff --git a/src/TNL/Meshes/MeshOrdering/KdTreeOrdering.h b/src/TNL/Meshes/MeshOrdering/KdTreeOrdering.h
new file mode 100644
index 0000000000000000000000000000000000000000..a3912daf924e8aa11b6e8f2c93b8ab200b69af69
--- /dev/null
+++ b/src/TNL/Meshes/MeshOrdering/KdTreeOrdering.h
@@ -0,0 +1,95 @@
+// Copyright (c) 2004-2022 Tomáš Oberhuber et al.
+//
+// This file is part of TNL - Template Numerical Library (https://tnl-project.org/)
+//
+// SPDX-License-Identifier: MIT
+
+// Implemented by: Jakub KlinkovskĂ˝
+
+#pragma once
+
+#ifdef HAVE_CGAL
+
+   #include <TNL/Meshes/Geometry/getEntityCenter.h>
+
+   // Documentation:
+   // - https://doc.cgal.org/latest/Spatial_searching/index.html
+   // - https://doc.cgal.org/latest/Spatial_searching/group__PkgSpatialSearchingDRef.html
+   // - https://doc.cgal.org/latest/Spatial_searching/classCGAL_1_1Kd__tree.html
+   #include <CGAL/Search_traits.h>
+   #include <CGAL/Search_traits_adapter.h>
+   #include <CGAL/Kd_tree.h>
+
+namespace TNL {
+namespace Meshes {
+
+struct KdTreeOrdering
+{
+   template< typename MeshEntity, typename Mesh, typename PermutationArray >
+   static void
+   getPermutations( const Mesh& mesh, PermutationArray& perm, PermutationArray& iperm )
+   {
+      static_assert( std::is_same< typename Mesh::DeviceType, TNL::Devices::Host >::value, "" );
+      static_assert( std::is_same< typename PermutationArray::DeviceType, TNL::Devices::Host >::value, "" );
+      using GlobalIndexType = typename Mesh::GlobalIndexType;
+      using PointType = typename Mesh::PointType;
+      using RealType = typename Mesh::RealType;
+
+      const GlobalIndexType numberOfEntities = mesh.template getEntitiesCount< MeshEntity >();
+
+      // allocate permutation vectors
+      perm.setSize( numberOfEntities );
+      iperm.setSize( numberOfEntities );
+
+      // wrapper for CGAL
+      struct Construct_coord_iterator
+      {
+         using result_type = const RealType*;
+         const RealType*
+         operator()( const PointType& p ) const
+         {
+            return static_cast< const RealType* >( &p[ 0 ] );
+         }
+         const RealType*
+         operator()( const PointType& p, int ) const
+         {
+            return static_cast< const RealType* >( &p[ 0 ] + p.getSize() );
+         }
+      };
+
+      // CGAL types
+      using DimTag = CGAL::Dimension_tag< PointType::getSize() >;
+      using Point_and_index = std::tuple< PointType, GlobalIndexType >;
+      using Traits_base = CGAL::Search_traits< RealType, PointType, const RealType*, Construct_coord_iterator, DimTag >;
+      using TreeTraits =
+         CGAL::Search_traits_adapter< Point_and_index, CGAL::Nth_of_tuple_property_map< 0, Point_and_index >, Traits_base >;
+      // Note: splitters affect the quality of the ordering
+      // available splitters are: https://doc.cgal.org/latest/Spatial_searching/classSplitter.html
+      using Splitter = CGAL::Sliding_fair< TreeTraits >;
+      using Tree = CGAL::Kd_tree< TreeTraits, Splitter >;
+
+      // build a k-d tree
+      Tree tree;
+      for( GlobalIndexType i = 0; i < numberOfEntities; i++ ) {
+         const auto& entity = mesh.template getEntity< MeshEntity >( i );
+         const auto center = getEntityCenter( mesh, entity );
+         tree.insert( std::make_tuple( center, i ) );
+      }
+      tree.build();
+
+      GlobalIndexType permIndex = 0;
+
+      // in-order traversal of the k-d tree
+      for( auto iter : tree ) {
+         const GlobalIndexType i = std::get< 1 >( iter );
+         perm[ permIndex ] = i;
+         iperm[ i ] = permIndex;
+         permIndex++;
+      }
+   }
+};
+
+}  // namespace Meshes
+}  // namespace TNL
+
+#endif  // HAVE_CGAL
diff --git a/src/Tools/CMakeLists.txt b/src/Tools/CMakeLists.txt
index 8bfc42752c5e0daf73f57423c52b8748fb650dc2..bf7152e58d687b32b4fe4db0fb8c412b7083bf97 100644
--- a/src/Tools/CMakeLists.txt
+++ b/src/Tools/CMakeLists.txt
@@ -10,6 +10,8 @@ set( MESH_TOOLS
          tnl-planar-correct-mesh
          tnl-refine-mesh
          tnl-game-of-life
+         tnl-reorder-mesh
+         tnl-plot-mesh-ordering
 )
 set( CUDA_MESH_TOOLS )
 if( BUILD_CUDA )
@@ -42,21 +44,33 @@ endforeach()
 
 find_package( METIS QUIET )
 if( METIS_FOUND )
-   add_executable(tnl-decompose-mesh tnl-decompose-mesh.cpp)
-   target_include_directories(tnl-decompose-mesh PUBLIC ${METIS_INCLUDE_DIRS})
-   target_link_libraries(tnl-decompose-mesh ${METIS_LIBRARIES})
+   add_executable( tnl-decompose-mesh tnl-decompose-mesh.cpp )
+   target_include_directories( tnl-decompose-mesh PUBLIC ${METIS_INCLUDE_DIRS} )
+   target_link_libraries( tnl-decompose-mesh ${METIS_LIBRARIES} )
    if( ZLIB_FOUND )
-      target_compile_definitions(tnl-decompose-mesh PUBLIC "-DHAVE_ZLIB")
-      target_include_directories(tnl-decompose-mesh PUBLIC ${ZLIB_INCLUDE_DIRS})
-      target_link_libraries(tnl-decompose-mesh ${ZLIB_LIBRARIES})
+      target_compile_definitions( tnl-decompose-mesh PUBLIC "-DHAVE_ZLIB" )
+      target_include_directories( tnl-decompose-mesh PUBLIC ${ZLIB_INCLUDE_DIRS} )
+      target_link_libraries( tnl-decompose-mesh ${ZLIB_LIBRARIES} )
    endif()
    if( tinyxml2_FOUND )
-      target_compile_definitions(tnl-decompose-mesh PUBLIC "-DHAVE_TINYXML2")
-      target_link_libraries(tnl-decompose-mesh tinyxml2::tinyxml2)
+      target_compile_definitions( tnl-decompose-mesh PUBLIC "-DHAVE_TINYXML2" )
+      target_link_libraries( tnl-decompose-mesh tinyxml2::tinyxml2 )
    endif()
    install( TARGETS tnl-decompose-mesh DESTINATION bin )
 endif()
 
+set( CGAL_DATA_DIR "/" )  # to disable stupid warning: https://github.com/CGAL/cgal/pull/6649
+find_package( CGAL )
+if( CGAL_FOUND )
+   option( CGAL_DO_NOT_WARN_ABOUT_CMAKE_BUILD_TYPE "disable warning" ON )
+
+   target_compile_definitions( tnl-reorder-mesh PUBLIC "-DHAVE_CGAL" )
+   target_link_libraries( tnl-reorder-mesh CGAL::CGAL )
+
+   target_compile_definitions( tnl-plot-mesh-ordering PUBLIC "-DHAVE_CGAL" )
+   target_link_libraries( tnl-plot-mesh-ordering CGAL::CGAL )
+endif()
+
 
 add_executable( tnl-grid-setup tnl-grid-setup.cpp )
 
diff --git a/src/Tools/tnl-plot-mesh-ordering.cpp b/src/Tools/tnl-plot-mesh-ordering.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1c063f7b857fc673f7f4e5367764f8564f8e6848
--- /dev/null
+++ b/src/Tools/tnl-plot-mesh-ordering.cpp
@@ -0,0 +1,151 @@
+#include <TNL/Meshes/Mesh.h>
+#include <TNL/Meshes/MeshOrdering.h>
+#include <TNL/Meshes/MeshOrdering/CuthillMcKeeOrdering.h>
+#include <TNL/Meshes/MeshOrdering/KdTreeOrdering.h>
+#include <TNL/Meshes/MeshOrdering/HilbertOrdering.h>
+#include <TNL/Meshes/Geometry/getEntityCenter.h>
+#include <TNL/Meshes/TypeResolver/resolveMeshType.h>
+#include <TNL/Meshes/Writers/VTUWriter.h>
+
+using namespace TNL;
+using namespace TNL::Meshes;
+
+struct MyConfigTag
+{};
+
+namespace TNL {
+namespace Meshes {
+namespace BuildConfigTags {
+
+/****
+ * Turn off support for float and long double.
+ */
+template<> struct GridRealTag< MyConfigTag, float >{ static constexpr bool enabled = false; };
+template<> struct GridRealTag< MyConfigTag, long double >{ static constexpr bool enabled = false; };
+
+/****
+ * Turn off support for short int and long int indexing.
+ */
+template<> struct GridIndexTag< MyConfigTag, short int >{ static constexpr bool enabled = false; };
+template<> struct GridIndexTag< MyConfigTag, long int >{ static constexpr bool enabled = false; };
+
+/****
+ * Unstructured meshes.
+ */
+template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Edge >{ static constexpr bool enabled = true; };
+template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Triangle >{ static constexpr bool enabled = true; };
+template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Quadrangle >{ static constexpr bool enabled = true; };
+template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Tetrahedron >{ static constexpr bool enabled = true; };
+template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Hexahedron >{ static constexpr bool enabled = true; };
+
+// Meshes are enabled only for the world dimension equal to the cell dimension.
+template< typename CellTopology, int WorldDimension >
+struct MeshSpaceDimensionTag< MyConfigTag, CellTopology, WorldDimension >
+{ static constexpr bool enabled = WorldDimension == CellTopology::dimension; };
+
+// Meshes are enabled only for types explicitly listed below.
+template<> struct MeshRealTag< MyConfigTag, float >{ static constexpr bool enabled = true; };
+template<> struct MeshRealTag< MyConfigTag, double >{ static constexpr bool enabled = true; };
+template<> struct MeshGlobalIndexTag< MyConfigTag, int >{ static constexpr bool enabled = true; };
+template<> struct MeshGlobalIndexTag< MyConfigTag, long int >{ static constexpr bool enabled = true; };
+template<> struct MeshLocalIndexTag< MyConfigTag, short int >{ static constexpr bool enabled = true; };
+
+}  // namespace BuildConfigTags
+}  // namespace Meshes
+}  // namespace TNL
+
+template< typename MeshEntity, typename MeshType >
+void
+writeMeshOrdering( const MeshType& mesh, const String& outputFileName )
+{
+   std::ofstream outputFile( outputFileName.getString() );
+
+   for( typename MeshType::GlobalIndexType i = 0; i < mesh.template getEntitiesCount< MeshEntity >(); i++ ) {
+      const auto cell = mesh.template getEntity< MeshEntity >( i );
+      const auto center = getEntityCenter( mesh, cell );
+      for( int j = 0; j < center.getSize(); j++ )
+         outputFile << center[ j ] << " ";
+      outputFile << std::endl;
+   }
+}
+
+template< typename Mesh >
+bool
+plotOrdering( const Mesh& mesh, const String& fileName )
+{
+   std::cerr << "plot-ordering is not implemented for meshes of type " << TNL::getType< Mesh >() << std::endl;
+   return false;
+}
+
+template< typename Mesh >
+void
+writeMeshVTU( const Mesh& mesh, const std::string& fileName )
+{
+   std::ofstream file( fileName );
+   using MeshWriter = TNL::Meshes::Writers::VTUWriter< Mesh >;
+   MeshWriter writer( file );
+   writer.template writeEntities< Mesh::getMeshDimension() >( mesh );
+}
+
+template< typename MeshConfig >
+bool
+plotOrdering( Mesh< MeshConfig, Devices::Host >& mesh, const String& fileName )
+{
+   using MeshType = Mesh< MeshConfig, Devices::Host >;
+   using MeshEntity = typename MeshType::Cell;
+
+   // strip the ".tnl" or ".vtk" or ".vtu" suffix
+   const String baseName( fileName.getString(), 0, fileName.getSize() - 4 );
+
+   writeMeshOrdering< MeshEntity >( mesh, baseName + ".original.gplt" );
+
+   // reverse Cuthill-McKee ordering
+   MeshOrdering< MeshType, CuthillMcKeeOrdering<> > rcm_ordering;
+   rcm_ordering.reorder( mesh );
+   writeMeshOrdering< MeshEntity >( mesh, baseName + ".rcm.gplt" );
+   writeMeshVTU( mesh, baseName + ".rcm.vtu" );
+
+#ifdef HAVE_CGAL
+   // k-d tree ordering
+   MeshOrdering< MeshType, KdTreeOrdering > kdtree_ordering;
+   kdtree_ordering.reorder( mesh );
+   writeMeshOrdering< MeshEntity >( mesh, baseName + ".kdtree.gplt" );
+   writeMeshVTU( mesh, baseName + ".kdtree.vtu" );
+#else
+   std::cerr << "CGAL support is missing. Skipping KdTreeOrdering." << std::endl;
+#endif
+
+#ifdef HAVE_CGAL
+   // Hilbert curve ordering
+   MeshOrdering< MeshType, HilbertOrdering > hilbert_ordering;
+   hilbert_ordering.reorder( mesh );
+   writeMeshOrdering< MeshEntity >( mesh, baseName + ".hilbert.gplt" );
+   writeMeshVTU( mesh, baseName + ".hilbert.vtu" );
+#else
+   std::cerr << "CGAL support is missing. Skipping HilbertOrdering." << std::endl;
+#endif
+
+   return true;
+}
+
+int
+main( int argc, char* argv[] )
+{
+   if( argc < 2 ) {
+      std::cerr << "Usage: " << argv[ 0 ] << " filename.[tnl|ng|vtk|vtu] ..." << std::endl;
+      return EXIT_FAILURE;
+   }
+
+   bool result = true;
+
+   for( int i = 1; i < argc; i++ ) {
+      String fileName = argv[ i ];
+      auto wrapper = [ & ]( auto& reader, auto&& mesh ) -> bool
+      {
+         return plotOrdering( mesh, fileName );
+      };
+      result &= resolveAndLoadMesh< MyConfigTag, Devices::Host >( wrapper, fileName );
+   }
+
+   return static_cast< int >( ! result );
+}
diff --git a/src/Tools/tnl-reorder-mesh.cpp b/src/Tools/tnl-reorder-mesh.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..31bfc95e43a7196053648002b9ce8b9a7e20b835
--- /dev/null
+++ b/src/Tools/tnl-reorder-mesh.cpp
@@ -0,0 +1,160 @@
+#include <experimental/filesystem>
+
+#include <TNL/Config/parseCommandLine.h>
+#include <TNL/Meshes/Mesh.h>
+#include <TNL/Meshes/MeshOrdering.h>
+#include <TNL/Meshes/MeshOrdering/CuthillMcKeeOrdering.h>
+#include <TNL/Meshes/MeshOrdering/KdTreeOrdering.h>
+#include <TNL/Meshes/MeshOrdering/HilbertOrdering.h>
+#include <TNL/Meshes/TypeResolver/resolveMeshType.h>
+#include <TNL/Meshes/Writers/VTKWriter.h>
+#include <TNL/Meshes/Writers/VTUWriter.h>
+
+using namespace TNL;
+using namespace TNL::Meshes;
+
+struct MyConfigTag
+{};
+
+namespace TNL {
+namespace Meshes {
+namespace BuildConfigTags {
+
+/****
+ * Turn off support for float and long double.
+ */
+template<> struct GridRealTag< MyConfigTag, float >{ static constexpr bool enabled = false; };
+template<> struct GridRealTag< MyConfigTag, long double >{ static constexpr bool enabled = false; };
+
+/****
+ * Turn off support for short int and long int indexing.
+ */
+template<> struct GridIndexTag< MyConfigTag, short int >{ static constexpr bool enabled = false; };
+template<> struct GridIndexTag< MyConfigTag, long int >{ static constexpr bool enabled = false; };
+
+/****
+ * Unstructured meshes.
+ */
+template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Edge >{ static constexpr bool enabled = true; };
+template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Triangle >{ static constexpr bool enabled = true; };
+template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Quadrangle >{ static constexpr bool enabled = true; };
+template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Tetrahedron >{ static constexpr bool enabled = true; };
+template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Hexahedron >{ static constexpr bool enabled = true; };
+
+// Meshes are enabled only for the world dimension equal to the cell dimension.
+template< typename CellTopology, int WorldDimension >
+struct MeshSpaceDimensionTag< MyConfigTag, CellTopology, WorldDimension >
+{ static constexpr bool enabled = WorldDimension == CellTopology::dimension; };
+
+// Meshes are enabled only for types explicitly listed below.
+template<> struct MeshRealTag< MyConfigTag, float >{ static constexpr bool enabled = true; };
+template<> struct MeshRealTag< MyConfigTag, double >{ static constexpr bool enabled = true; };
+template<> struct MeshGlobalIndexTag< MyConfigTag, int >{ static constexpr bool enabled = true; };
+template<> struct MeshGlobalIndexTag< MyConfigTag, long int >{ static constexpr bool enabled = true; };
+template<> struct MeshLocalIndexTag< MyConfigTag, short int >{ static constexpr bool enabled = true; };
+
+}  // namespace BuildConfigTags
+}  // namespace Meshes
+}  // namespace TNL
+
+template< typename Mesh >
+bool
+reorder( Mesh&& mesh, const std::string& ordering, const std::string& outputFileName, std::string outputFileFormat )
+{
+   if( ordering == "rcm" ) {
+      // reverse Cuthill-McKee ordering
+      MeshOrdering< Mesh, CuthillMcKeeOrdering<> > rcm_ordering;
+      rcm_ordering.reorder( mesh );
+   }
+   else if( ordering == "kdtree" ) {
+#ifdef HAVE_CGAL
+      // k-d tree ordering
+      MeshOrdering< Mesh, KdTreeOrdering > kdtree_ordering;
+      kdtree_ordering.reorder( mesh );
+#else
+      std::cerr << "CGAL support is missing. Recompile with -DHAVE_CGAL and try again." << std::endl;
+      return false;
+#endif
+   }
+   else if( ordering == "hilbert" ) {
+#ifdef HAVE_CGAL
+      // Hilbert curve ordering
+      MeshOrdering< Mesh, HilbertOrdering > hilbert_ordering;
+      hilbert_ordering.reorder( mesh );
+#else
+      std::cerr << "CGAL support is missing. Recompile with -DHAVE_CGAL and try again." << std::endl;
+      return false;
+#endif
+   }
+   else {
+      std::cerr << "unknown ordering algorithm: '" << ordering << "'. Available options are: rcm, kdtree, hilbert."
+                << std::endl;
+      return false;
+   }
+
+   namespace fs = std::experimental::filesystem;
+   if( outputFileFormat == "auto" ) {
+      outputFileFormat = fs::path( outputFileName ).extension();
+      if( outputFileFormat.length() > 0 )
+         // remove dot from the extension
+         outputFileFormat = outputFileFormat.substr( 1 );
+   }
+
+   if( outputFileFormat == "vtk" ) {
+      std::ofstream file( outputFileName );
+      using MeshWriter = TNL::Meshes::Writers::VTKWriter< Mesh >;
+      MeshWriter writer( file );
+      writer.template writeEntities< Mesh::getMeshDimension() >( mesh );
+   }
+   else if( outputFileFormat == "vtu" ) {
+      std::ofstream file( outputFileName );
+      using MeshWriter = TNL::Meshes::Writers::VTUWriter< Mesh >;
+      MeshWriter writer( file );
+      writer.template writeEntities< Mesh::getMeshDimension() >( mesh );
+   }
+   else {
+      std::cerr << "unknown output file format: '" << outputFileFormat << "'. Available options are: vtk, vtu." << std::endl;
+      return false;
+   }
+
+   return true;
+}
+
+void
+configSetup( Config::ConfigDescription& config )
+{
+   config.addDelimiter( "General settings:" );
+   config.addRequiredEntry< String >( "input-mesh", "Input file with the distributed mesh." );
+   config.addEntry< String >( "input-mesh-format", "Input mesh file format.", "auto" );
+   config.addRequiredEntry< String >( "ordering", "Algorithm to reorder the mesh." );
+   config.addEntryEnum( "rcm" );
+   config.addEntryEnum( "kdtree" );
+   config.addRequiredEntry< String >( "output-mesh", "Output file with the distributed mesh." );
+   config.addEntry< String >( "output-mesh-format", "Output mesh file format.", "auto" );
+}
+
+int
+main( int argc, char* argv[] )
+{
+   Config::ParameterContainer parameters;
+   Config::ConfigDescription conf_desc;
+
+   configSetup( conf_desc );
+
+   if( ! parseCommandLine( argc, argv, conf_desc, parameters ) )
+      return EXIT_FAILURE;
+
+   const String inputFileName = parameters.getParameter< String >( "input-mesh" );
+   const String inputFileFormat = parameters.getParameter< String >( "input-mesh-format" );
+   const String ordering = parameters.getParameter< String >( "ordering" );
+   const String outputFileName = parameters.getParameter< String >( "output-mesh" );
+   const String outputFileFormat = parameters.getParameter< String >( "output-mesh-format" );
+
+   auto wrapper = [ & ]( auto& reader, auto&& mesh ) -> bool
+   {
+      using MeshType = std::decay_t< decltype( mesh ) >;
+      return reorder( std::forward< MeshType >( mesh ), ordering, outputFileName, outputFileFormat );
+   };
+   const bool status = resolveAndLoadMesh< MyConfigTag, Devices::Host >( wrapper, inputFileName, inputFileFormat );
+   return static_cast< int >( ! status );
+}