diff --git a/examples/heat-equation/heatEquationSolver.h b/examples/heat-equation/heatEquationSolver.h
index 3cc97549e7012ae69923f07f4f2a6a1221f42076..9f483faa06bd81b3d1d9d7515bd03d2ebf8227d6 100644
--- a/examples/heat-equation/heatEquationSolver.h
+++ b/examples/heat-equation/heatEquationSolver.h
@@ -24,6 +24,7 @@
 #include <core/tnlLogger.h>
 #include <core/vectors/tnlVector.h>
 #include <core/vectors/tnlSharedVector.h>
+#include <solvers/pde/tnlExplicitUpdater.h>
 #include "heatEquationSolver.h"
 #include "tnlAnalyticSolution.h"
 
@@ -86,7 +87,9 @@ class heatEquationSolver
                                                       exactSolution,
                                                       analyticLaplace,
                                                       numericalLaplace;
-   //MeshType mesh;
+
+   tnlExplicitUpdater< Mesh, DofVectorType, BoundaryCondition, Diffusion > explicitUpdater;
+
    AnalyticSpaceFunction analyticSpaceFunction;
    TimeFunction timeFunction;
    AnalyticSolution< MeshType, RealType, IndexType > analyticSolution;
diff --git a/examples/heat-equation/heatEquationSolver_impl.h b/examples/heat-equation/heatEquationSolver_impl.h
index 68d6f2621125fa3fddd2497b1fec9c364c5cf4cd..6b6f2ae774f0505e29ffd5aaa9698470b9b1cf90 100644
--- a/examples/heat-equation/heatEquationSolver_impl.h
+++ b/examples/heat-equation/heatEquationSolver_impl.h
@@ -21,6 +21,7 @@
 #include <core/mfilename.h>
 #include "heatEquationSolver.h"
 
+
 template< typename Mesh, typename Diffusion, typename BoundaryCondition, typename RightHandSide, 
           typename TimeFunction, typename AnalyticSpaceFunction>
 tnlString heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide,TimeFunction,AnalyticSpaceFunction> 
@@ -219,25 +220,13 @@ void heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide,TimeFunc
     */
 
    this->bindDofs( mesh, _u );
-   if( DeviceType :: getDevice() == tnlHostDevice )
-   {
-      boundaryCondition.applyBoundaryTimeDerivation(mesh, _fu, time, timeFunction, analyticSpaceFunction);
-      
-      boundaryCondition.applyBoundaryConditions(mesh, _u, time, timeFunction, analyticSpaceFunction);
-      
-      diffusion.getExplicitRHS(mesh,_u,_fu);  
-      
-      //RHS.applyRHSValues(mesh, time, _fu, timeFunction, analyticSpaceFunction);
-      
-   }
-#ifdef HAVE_CUDA
-   if( DeviceType :: getDevice() == tnlCudaDevice )
-   {
-      /****
-       * Write the CUDA solver here.
-       */
-   }
-#endif
+   explicitUpdater.template update< Mesh::Dimensions >( time,
+                                                        tau,
+                                                        mesh,
+                                                        this->boundaryCondition,
+                                                        this->diffusion,
+                                                        _u,
+                                                        _fu );
 }
 
 template< typename Mesh, typename Diffusion, typename BoundaryCondition, typename RightHandSide, 
diff --git a/examples/heat-equation/tnlDirichletBoundaryConditions.h b/examples/heat-equation/tnlDirichletBoundaryConditions.h
index b2a5c430703a2cf7e5066afdb83c9e6f5b289866..926a96cfd976be59f7fc0581379f8426f01c3fde 100644
--- a/examples/heat-equation/tnlDirichletBoundaryConditions.h
+++ b/examples/heat-equation/tnlDirichletBoundaryConditions.h
@@ -30,6 +30,16 @@ class tnlDirichletBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >,
    typedef tnlStaticVector< 1, RealType > VertexType;
    typedef typename MeshType::CoordinatesType CoordinatesType;
             
+   void setBoundaryConditions( const RealType& time,
+                               const RealType& tau,
+                               const MeshType& mesh,
+                               const IndexType index,
+                               const CoordinatesType& coordinates,
+                               DofVectorType& u,
+                               DofVectorType& fu );
+
+   //--------------------------------
+
    template< typename AnalyticSpaceFunction,
              typename TimeFunction,
              typename Vector >
@@ -47,6 +57,8 @@ class tnlDirichletBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >,
                                      const RealType& time,
                                      TimeFunction& timeFunction,
                                      AnalyticSpaceFunction& analyticSpaceFunction );
+   //-----------------------------------
+
 };
 
 template< typename MeshReal,
@@ -67,7 +79,17 @@ class tnlDirichletBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >,
    typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
    typedef tnlStaticVector< 2, RealType > VertexType;
    typedef typename MeshType::CoordinatesType CoordinatesType;
-            
+
+   void setBoundaryConditions( const RealType& time,
+                               const RealType& tau,
+                               const MeshType& mesh,
+                               const IndexType index,
+                               const CoordinatesType& coordinates,
+                               DofVectorType& u,
+                               DofVectorType& fu );
+
+
+   //----------------------------------------
    template< typename AnalyticSpaceFunction,
              typename TimeFunction,
              typename Vector >
@@ -85,6 +107,7 @@ class tnlDirichletBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >,
                                      const RealType& time,
                                      TimeFunction& timeFunction,
                                      AnalyticSpaceFunction& analyticSpaceFunction );
+   //----------------------------------
 };
 
 template< typename MeshReal,
@@ -105,7 +128,16 @@ class tnlDirichletBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >,
    typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
    typedef tnlStaticVector< 3, RealType > VertexType;
    typedef typename MeshType::CoordinatesType CoordinatesType;
-            
+
+   void setBoundaryConditions( const RealType& time,
+                               const RealType& tau,
+                               const MeshType& mesh,
+                               const IndexType index,
+                               const CoordinatesType& coordinates,
+                               DofVectorType& u,
+                               DofVectorType& fu );
+
+   //--------------------------------------------------
    template< typename AnalyticSpaceFunction,
              typename TimeFunction,
              typename Vector >
@@ -123,6 +155,7 @@ class tnlDirichletBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >,
                                      const RealType& time,
                                      TimeFunction& timeFunction,
                                      AnalyticSpaceFunction& analyticSpaceFunction );
+   //--------------------------
 };
 
 #include "tnlDirichletBoundaryConditions_impl.h"
diff --git a/examples/heat-equation/tnlDirichletBoundaryConditions_impl.h b/examples/heat-equation/tnlDirichletBoundaryConditions_impl.h
index d013d7201533196ba7023e8a532945aa6ec9ab5e..18ea0fd52c44d0d32f3429cd40a24a4e00c6dc30 100644
--- a/examples/heat-equation/tnlDirichletBoundaryConditions_impl.h
+++ b/examples/heat-equation/tnlDirichletBoundaryConditions_impl.h
@@ -3,6 +3,26 @@
 
 #include "tnlDirichletBoundaryConditions.h"
 
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+void
+tnlDirichletBoundaryConditions< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index >::
+setBoundaryConditions( const RealType& time,
+                       const RealType& tau,
+                       const MeshType& mesh,
+                       const IndexType index,
+                       const CoordinatesType& coordinates,
+                       DofVectorType& u,
+                       DofVectorType& fu )
+{
+   fu[ index ] = 0;
+   u[ index ] = 0;
+}
+
+
 template< typename MeshReal,
           typename Device,
           typename MeshIndex,
@@ -65,6 +85,26 @@ applyBoundaryTimeDerivation( const MeshType& mesh,
    u[u.getSize()-1] = timeFunctionDerivationValue*analyticSpaceFunction.getF(vertex);
 }
 
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+void
+tnlDirichletBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+setBoundaryConditions( const RealType& time,
+                       const RealType& tau,
+                       const MeshType& mesh,
+                       const IndexType index,
+                       const CoordinatesType& coordinates,
+                       DofVectorType& u,
+                       DofVectorType& fu )
+{
+   fu[ index ] = 0;
+   u[ index ] = 0;
+}
+
 template< typename MeshReal,
           typename Device,
           typename MeshIndex,
@@ -177,6 +217,25 @@ applyBoundaryTimeDerivation( const MeshType& mesh,
    }
 }
 
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+void
+tnlDirichletBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
+setBoundaryConditions( const RealType& time,
+                       const RealType& tau,
+                       const MeshType& mesh,
+                       const IndexType index,
+                       const CoordinatesType& coordinates,
+                       DofVectorType& u,
+                       DofVectorType& fu )
+{
+   fu[ index ] = 0;
+   u[ index ] = 0;
+}
+
 template< typename MeshReal,
           typename Device,
           typename MeshIndex,
diff --git a/examples/heat-equation/tnlLinearDiffusion.h b/examples/heat-equation/tnlLinearDiffusion.h
index e35bdedf831e5d44fdae4d15b7a6013178eae65b..0263b3c2b6ec002fe3cf8ebd1be50a32462dbed9 100644
--- a/examples/heat-equation/tnlLinearDiffusion.h
+++ b/examples/heat-equation/tnlLinearDiffusion.h
@@ -28,6 +28,15 @@ class tnlLinearDiffusion< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index
    typedef Index IndexType;
    //typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
 
+   template< typename Vector >
+   void update( const RealType& time,
+                const RealType& tau,
+                const MeshType& mesh,
+                const IndexType index,
+                const CoordinatesType& coordinates,
+                Vector& u,
+                Vector& fu );
+
    template< typename Vector >
    void getExplicitRHS( const MeshType& mesh,
                         const CoordinatesType& coordinates,
@@ -58,6 +67,16 @@ class tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index
    typedef Index IndexType;
    //typedef tnlVector< RealType, DeviceType, IndexType > DofVectorType;
 
+   template< typename Vector >
+   void update( const RealType& time,
+                const RealType& tau,
+                const MeshType& mesh,
+                const IndexType index,
+                const CoordinatesType& coordinates,
+                Vector& u,
+                Vector& fu );
+
+
    template< typename Vector >
    void getExplicitRHS( const MeshType& mesh,
                         const CoordinatesType& coordinates,
@@ -88,6 +107,16 @@ class tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index
    typedef Index IndexType;
    //typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
 
+   template< typename Vector >
+   void update( const RealType& time,
+                const RealType& tau,
+                const MeshType& mesh,
+                const IndexType index,
+                const CoordinatesType& coordinates,
+                Vector& u,
+                Vector& fu );
+
+
    template< typename Vector >
    void getExplicitRHS( const MeshType& mesh,
                         const CoordinatesType& coordinates,
diff --git a/examples/heat-equation/tnlLinearDiffusion_impl.h b/examples/heat-equation/tnlLinearDiffusion_impl.h
index 669d2e2ef0538583c2d4c54ee098296c9e57e072..a9db826923f7f82bea455b22ce9dbfa0284dc60d 100644
--- a/examples/heat-equation/tnlLinearDiffusion_impl.h
+++ b/examples/heat-equation/tnlLinearDiffusion_impl.h
@@ -5,6 +5,26 @@
 #include "tnlLinearDiffusion.h"
 #include <mesh/tnlGrid.h>
 
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+   template< typename Vector >
+void
+tnlLinearDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
+update( const RealType& time,
+        const RealType& tau,
+        const MeshType& mesh,
+        const IndexType index,
+        const CoordinatesType& coordinates,
+        Vector& u,
+        Vector& fu )
+{
+
+}
+
+
 template< typename MeshReal,
           typename Device,
           typename MeshIndex,
@@ -53,6 +73,27 @@ getExplicitRHS( const MeshType& mesh,
    }
 }
 
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+   template< typename Vector >
+void
+tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+update( const RealType& time,
+        const RealType& tau,
+        const MeshType& mesh,
+        const IndexType index,
+        const CoordinatesType& coordinates,
+        Vector& u,
+        Vector& fu )
+{
+
+}
+
+
 template< typename MeshReal,
           typename Device,
           typename MeshIndex,
@@ -112,6 +153,26 @@ getExplicitRHS( const MeshType& mesh,
    }
 }
 
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+   template< typename Vector >
+void
+tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
+update( const RealType& time,
+        const RealType& tau,
+        const MeshType& mesh,
+        const IndexType index,
+        const CoordinatesType& coordinates,
+        Vector& u,
+        Vector& fu )
+{
+
+}
+
+
 
 template< typename MeshReal,
           typename Device,
diff --git a/src/implementation/mesh/CMakeLists.txt b/src/implementation/mesh/CMakeLists.txt
index e26aab5112a17e54d1b31915d94c009de3e8fdf0..5a46422a511a93f5ce8268809ee135d12418bc1c 100755
--- a/src/implementation/mesh/CMakeLists.txt
+++ b/src/implementation/mesh/CMakeLists.txt
@@ -1,6 +1,9 @@
 SET( headers tnlGrid1D_impl.h
              tnlGrid2D_impl.h
-             tnlGrid3D_impl.h )
+             tnlGrid3D_impl.h
+             tnlTraversal_Grid1D_impl.h
+             tnlTraversal_Grid2D_impl.h
+             tnlTraversal_Grid3D_impl.h )
 
 SET( CURRENT_DIR ${CMAKE_SOURCE_DIR}/src/implementation/mesh )    
 SET( tnl_implementation_mesh_SOURCES
diff --git a/src/implementation/mesh/tnlGrid1D_impl.h b/src/implementation/mesh/tnlGrid1D_impl.h
index 19dc6aba7b5d226a0418f324aae4bccc84b71c6e..c9acc9456655f72d23f37bb9fe9a357368e8c829 100644
--- a/src/implementation/mesh/tnlGrid1D_impl.h
+++ b/src/implementation/mesh/tnlGrid1D_impl.h
@@ -146,7 +146,7 @@ Index tnlGrid< 1, Real, Device, Index > :: getCellIndex( const CoordinatesType&
    tnlAssert( cellCoordinates.x() >= 0 && cellCoordinates.x() < this->getDimensions().x(),
               cerr << "cellCoordinates.x() = " << cellCoordinates.x()
                    << " this->getDimensions().x() = " << this->getDimensions().x()
-                   << " this->getName() = " << this->getName(); )
+                   << " this->getName() = " << this->getName(); );
    return cellCoordinates.x();
 }
 
@@ -232,6 +232,9 @@ Vertex tnlGrid< 1, Real, Device, Index >::getVertex( const CoordinatesType& vert
 template< typename Real,
           typename Device,
           typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
 Index tnlGrid< 1, Real, Device, Index > :: getNumberOfCells() const
 {
    return this->numberOfCells;
@@ -240,11 +243,48 @@ Index tnlGrid< 1, Real, Device, Index > :: getNumberOfCells() const
 template< typename Real,
           typename Device,
           typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
 Index tnlGrid< 1, Real, Device, Index > :: getNumberOfVertices() const
 {
    return this->numberOfVertices;
 };
 
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+bool tnlGrid< 1, Real, Device, Index > :: isBoundaryCell( const CoordinatesType& cellCoordinates ) const
+{
+   tnlAssert( cellCoordinates.x() >= 0 && cellCoordinates.x() < this->getDimensions().x(),
+              cerr << "cellCoordinates.x() = " << cellCoordinates.x()
+                   << " this->getDimensions().x() = " << this->getDimensions().x()
+                   << " this->getName() = " << this->getName(); );
+   if( cellCoordinates.x() == 0 || cellCoordinates.x() == this->getDimensions().x() - 1 )
+      return true;
+   return false;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+bool tnlGrid< 1, Real, Device, Index > :: isBoundaryVertex( const CoordinatesType& vertexCoordinates ) const
+{
+   tnlAssert( vertexCoordinates.x() >= 0 && vertexCoordinates.x() < this->getDimensions().x() + 1,
+              cerr << "vertexCoordinates.x() = " << vertexCoordinates.x()
+                   << " this->getDimensions().x() + 1 = " << this->getDimensions().x() + 1
+                   << " this->getName() = " << this->getName(); );
+   if( vertexCoordinates.x() == 0 || vertexCoordinates.x() == this->getDimensions().x() )
+      return true;
+   return false;
+}
+
 template< typename Real,
           typename Device,
           typename Index >
diff --git a/src/implementation/mesh/tnlGrid2D_impl.h b/src/implementation/mesh/tnlGrid2D_impl.h
index 28fd0d766891b5c7083e7b8182ce79c78f9f083d..afc827381f57ee9d02ce6f02ac75752850455eb7 100644
--- a/src/implementation/mesh/tnlGrid2D_impl.h
+++ b/src/implementation/mesh/tnlGrid2D_impl.h
@@ -365,6 +365,9 @@ Vertex tnlGrid< 2, Real, Device, Index >::getVertex( const CoordinatesType& vert
 template< typename Real,
           typename Device,
           typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
 Index tnlGrid< 2, Real, Device, Index > :: getNumberOfCells() const
 {
    return this->numberOfCells;
@@ -373,6 +376,9 @@ Index tnlGrid< 2, Real, Device, Index > :: getNumberOfCells() const
 template< typename Real,
           typename Device,
           typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
 Index tnlGrid< 2, Real, Device, Index > :: getNumberOfFaces() const
 {
    return this->numberOfFaces;
@@ -381,9 +387,96 @@ Index tnlGrid< 2, Real, Device, Index > :: getNumberOfFaces() const
 template< typename Real,
           typename Device,
           typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
 Index tnlGrid< 2, Real, Device, Index > :: getNumberOfVertices() const
 {
-   return numberOfVertices;
+   return this->numberOfVertices;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+bool tnlGrid< 2, Real, Device, Index > :: isBoundaryCell( const CoordinatesType& cellCoordinates ) const
+{
+   tnlAssert( cellCoordinates.x() >= 0 && cellCoordinates.x() < this->getDimensions().x(),
+              cerr << "cellCoordinates.x() = " << cellCoordinates.x()
+                   << " this->getDimensions().x() = " << this->getDimensions().x()
+                   << " this->getName() = " << this->getName(); );
+   tnlAssert( cellCoordinates.y() >= 0 && cellCoordinates.y() < this->getDimensions().y(),
+              cerr << "cellCoordinates.y() = " << cellCoordinates.y()
+                   << " this->getDimensions().y() = " << this->getDimensions().y()
+                   << " this->getName() = " << this->getName(); );
+
+   if( cellCoordinates.x() == 0 || cellCoordinates.x() == this->getDimensions().x() - 1 ||
+       cellCoordinates.y() == 0 || cellCoordinates.y() == this->getDimensions().y() - 1 )
+      return true;
+   return false;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+   template< int nx, int ny >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+bool tnlGrid< 2, Real, Device, Index > :: isBoundaryFace( const CoordinatesType& faceCoordinates ) const
+{
+   tnlStaticAssert( nx >= 0 && ny >= 0 && nx + ny == 1, "Wrong template parameters nx or ny." );
+   if( nx )
+   {
+      tnlAssert( faceCoordinates.x() >= 0 && faceCoordinates.x() < this->getDimensions().x() + 1,
+                 cerr << "faceCoordinates.x() = " << faceCoordinates.x()
+                      << " this->getDimensions().x() + 1 = " << this->getDimensions().x() + 1
+                      << " this->getName() = " << this->getName(); );
+      tnlAssert( faceCoordinates.y() >= 0 && faceCoordinates.y() < this->getDimensions().y(),
+                 cerr << "faceCoordinates.y() = " << faceCoordinates.y()
+                      << " this->getDimensions().y() = " << this->getDimensions().y()
+                      << " this->getName() = " << this->getName(); );
+      if( faceCoordinates.x() == 0 || faceCoordinates.x() == this->getDimensions().x() )
+         return true;
+      return false;
+   }
+   tnlAssert( faceCoordinates.x() >= 0 && faceCoordinates.x() < this->getDimensions().x(),
+              cerr << "faceCoordinates.x() = " << faceCoordinates.x()
+                   << " this->getDimensions().x() = " << this->getDimensions().x()
+                   << " this->getName() = " << this->getName(); );
+   tnlAssert( faceCoordinates.y() >= 0 && faceCoordinates.y() < this->getDimensions().y() + 1,
+              cerr << "faceCoordinates.y() = " << faceCoordinates.y()
+                   << " this->getDimensions().y() + 1 = " << this->getDimensions().y() + 1
+                   << " this->getName() = " << this->getName(); );
+   if( faceCoordinates.y() == 0 || faceCoordinates.y() == this->getDimensions().y() )
+      return true;
+   return false;
+}
+
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+bool tnlGrid< 2, Real, Device, Index > :: isBoundaryVertex( const CoordinatesType& vertexCoordinates ) const
+{
+   tnlAssert( vertexCoordinates.x() >= 0 && vertexCoordinates.x() < this->getDimensions().x() + 1,
+              cerr << "vertexCoordinates.x() = " << vertexCoordinates.x()
+                   << " this->getDimensions().x() + 1 = " << this->getDimensions().x() + 1
+                   << " this->getName() = " << this->getName(); );
+   tnlAssert( vertexCoordinates.y() >= 0 && vertexCoordinates.y() < this->getDimensions().y() + 1,
+              cerr << "vertexCoordinates.y() = " << vertexCoordinates.y()
+                   << " this->getDimensions().y() + 1 = " << this->getDimensions().y() + 1
+                   << " this->getName() = " << this->getName(); );
+
+   if( vertexCoordinates.x() == 0 || vertexCoordinates.x() == this->getDimensions().x() ||
+       vertexCoordinates.y() == 0 || vertexCoordinates.y() == this->getDimensions().y() )
+      return true;
+   return false;
 }
 
 template< typename Real,
diff --git a/src/implementation/mesh/tnlGrid3D_impl.h b/src/implementation/mesh/tnlGrid3D_impl.h
index e773aad74afc05f708d572aa42808c84cd9af7e2..3828e4b571d1b25325b4db923ea871290627eb5e 100644
--- a/src/implementation/mesh/tnlGrid3D_impl.h
+++ b/src/implementation/mesh/tnlGrid3D_impl.h
@@ -654,6 +654,9 @@ Vertex tnlGrid< 3, Real, Device, Index >::getVertex( const CoordinatesType& vert
 template< typename Real,
           typename Device,
           typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
 Index tnlGrid< 3, Real, Device, Index > :: getNumberOfCells() const
 {
    return this->numberOfCells;
@@ -662,6 +665,9 @@ Index tnlGrid< 3, Real, Device, Index > :: getNumberOfCells() const
 template< typename Real,
           typename Device,
           typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
 Index tnlGrid< 3, Real, Device, Index > :: getNumberOfFaces() const
 {
    return this->numberOfFaces;
@@ -670,6 +676,9 @@ Index tnlGrid< 3, Real, Device, Index > :: getNumberOfFaces() const
 template< typename Real,
           typename Device,
           typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
 Index tnlGrid< 3, Real, Device, Index > :: getNumberOfEdges() const
 {
    return this->numberOfEdges;
@@ -678,11 +687,199 @@ Index tnlGrid< 3, Real, Device, Index > :: getNumberOfEdges() const
 template< typename Real,
           typename Device,
           typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
 Index tnlGrid< 3, Real, Device, Index > :: getNumberOfVertices() const
 {
    return numberOfVertices;
 }
 
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+bool tnlGrid< 3, Real, Device, Index > :: isBoundaryCell( const CoordinatesType& cellCoordinates ) const
+{
+   tnlAssert( cellCoordinates.x() >= 0 && cellCoordinates.x() < this->getDimensions().x(),
+              cerr << "cellCoordinates.x() = " << cellCoordinates.x()
+                   << " this->getDimensions().x() = " << this->getDimensions().x()
+                   << " this->getName() = " << this->getName(); );
+   tnlAssert( cellCoordinates.y() >= 0 && cellCoordinates.y() < this->getDimensions().y(),
+              cerr << "cellCoordinates.y() = " << cellCoordinates.y()
+                   << " this->getDimensions().y() = " << this->getDimensions().y()
+                   << " this->getName() = " << this->getName(); );
+   tnlAssert( cellCoordinates.z() >= 0 && cellCoordinates.z() < this->getDimensions().z(),
+              cerr << "cellCoordinates.z() = " << cellCoordinates.z()
+                   << " this->getDimensions().z() = " << this->getDimensions().z()
+                   << " this->getName() = " << this->getName(); );
+
+
+   if( cellCoordinates.x() == 0 || cellCoordinates.x() == this->getDimensions().x() - 1 ||
+       cellCoordinates.y() == 0 || cellCoordinates.y() == this->getDimensions().y() - 1 ||
+       cellCoordinates.z() == 0 || cellCoordinates.z() == this->getDimensions().z() - 1 )
+      return true;
+   return false;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+   template< int nx, int ny, int nz >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+bool tnlGrid< 3, Real, Device, Index > :: isBoundaryFace( const CoordinatesType& faceCoordinates ) const
+{
+   tnlStaticAssert( nx >= 0 && ny >= 0 && nz >=0 && nx + ny == 1, "Wrong template parameters nx or ny." );
+   if( nx )
+   {
+      tnlAssert( faceCoordinates.x() >= 0 && faceCoordinates.x() < this->getDimensions().x() + 1,
+                 cerr << "faceCoordinates.x() = " << faceCoordinates.x()
+                      << " this->getDimensions().x() + 1 = " << this->getDimensions().x() + 1
+                      << " this->getName() = " << this->getName(); );
+      tnlAssert( faceCoordinates.y() >= 0 && faceCoordinates.y() < this->getDimensions().y(),
+                 cerr << "faceCoordinates.y() = " << faceCoordinates.y()
+                      << " this->getDimensions().y() = " << this->getDimensions().y()
+                      << " this->getName() = " << this->getName(); );
+      if( faceCoordinates.x() == 0 || faceCoordinates.x() == this->getDimensions().y() )
+         return true;
+      return false;
+   }
+   if( ny )
+   {
+      tnlAssert( faceCoordinates.x() >= 0 && faceCoordinates.x() < this->getDimensions().x(),
+                 cerr << "faceCoordinates.x() = " << faceCoordinates.x()
+                      << " this->getDimensions().x() = " << this->getDimensions().x()
+                      << " this->getName() = " << this->getName(); );
+      tnlAssert( faceCoordinates.y() >= 0 && faceCoordinates.y() < this->getDimensions().y() + 1,
+                 cerr << "faceCoordinates.y() = " << faceCoordinates.y()
+                      << " this->getDimensions().y() + 1 = " << this->getDimensions().y() + 1
+                      << " this->getName() = " << this->getName(); );
+      if( faceCoordinates.y() == 0 || faceCoordinates.y() == this->getDimensions().y() )
+         return true;
+      return false;
+   }
+   if( nz )
+   {
+      tnlAssert( faceCoordinates.x() >= 0 && faceCoordinates.x() < this->getDimensions.x(),
+                 cerr << "faceCoordinates.x() = " << faceCoordinates.x()
+                      << " this->getDimensions().x() = " << this->getDimensions().x()
+                      << " this->getName() = " << this->getName(); );
+      tnlAssert( faceCoordinates.y() >= 0 && faceCoordinates.y() < this->getDimensions.y(),
+                 cerr << "faceCoordinates.y() = " << faceCoordinates.y()
+                      << " this->getDimensions().y()= " << this->getDimensions().y()
+                      << " this->getName() = " << this->getName(); );
+      tnlAssert( faceCoordinates.z() >= 0 && faceCoordinates.z() < this->getDimensions.z() + 1,
+                 cerr << "faceCoordinates.z() = " << faceCoordinates.z()
+                      << " this->getDimensions().z() + 1 = " << this->getDimensions().z() + 1
+                      << " this->getName() = " << this->getName(); );
+      if( faceCoordinates.z() == 0 || faceCoordinates.z() == this->getDimensions().z() )
+         return true;
+      return false;
+   }
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+   template< int dx, int dy, int dz >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+bool tnlGrid< 3, Real, Device, Index >::isBoundaryEdge( const CoordinatesType& edgeCoordinates ) const
+{
+   tnlStaticAssert( dx >= 0 && dy >= 0 && dz >= 0 && dx + dy + dz = 1, "Wrong template parameters nx or ny or nz." );
+   if( dx )
+   {
+      tnlAssert( edgeCoordinates.x() >= 0 && edgeCoordinates.x() < this->getDimensions.x(),
+                 cerr << "edgeCoordinates.x() = " << edgeCoordinates.x()
+                      << " this->getDimensions().x() = " << this->getDimensions().x()
+                      << " this->getName() = " << this->getName(); );
+      tnlAssert( edgeCoordinates.y() >= 0 && edgeCoordinates.y() < this->getDimensions.y() + 1,
+                 cerr << "edgeCoordinates.y() = " << edgeCoordinates.y()
+                      << " this->getDimensions().y() + 1 = " << this->getDimensions().y() + 1
+                      << " this->getName() = " << this->getName(); );
+      tnlAssert( edgeCoordinates.z() >= 0 && edgeCoordinates.z() < this->getDimensions.z() + 1,
+                 cerr << "edgeCoordinates.z() = " << edgeCoordinates.z()
+                      << " this->getDimensions().z() + 1 = " << this->getDimensions().z() + 1
+                      << " this->getName() = " << this->getName(); );
+      if( edgeCoordinates.y() == 0 || edgeCoordinates.y() == this->getDimensions().y() ||
+          edgeCoordinates.z() == 0 || edgeCoordinates.z() == this->getDimensions().z() )
+         return true;
+      return false;
+   }
+   if( dy )
+   {
+      tnlAssert( edgeCoordinates.x() >= 0 && edgeCoordinates.x() < this->getDimensions.x() + 1,
+                 cerr << "edgeCoordinates.x() = " << edgeCoordinates.x()
+                      << " this->getDimensions().x() + 1 = " << this->getDimensions().x() + 1
+                      << " this->getName() = " << this->getName(); );
+      tnlAssert( edgeCoordinates.y() >= 0 && edgeCoordinates.y() < this->getDimensions.y(),
+                 cerr << "edgeCoordinates.y() = " << edgeCoordinates.y()
+                      << " this->getDimensions().y() = " << this->getDimensions().y()
+                      << " this->getName() = " << this->getName(); );
+      tnlAssert( edgeCoordinates.z() >= 0 && edgeCoordinates.z() < this->getDimensions.z() + 1,
+                 cerr << "edgeCoordinates.z() = " << edgeCoordinates.z()
+                      << " this->getDimensions().z() + 1 = " << this->getDimensions().z() + 1
+                      << " this->getName() = " << this->getName(); );
+      if( edgeCoordinates.x() == 0 || edgeCoordinates.x() == this->getDimensions().x() ||
+          edgeCoordinates.z() == 0 || edgeCoordinates.z() == this->getDimensions().z() )
+         return true;
+      return false;
+
+   }
+   if( dz )
+   {
+      tnlAssert( edgeCoordinates.x() >= 0 && edgeCoordinates.x() < this->getDimensions.x() + 1,
+                 cerr << "edgeCoordinates.x() = " << edgeCoordinates.x()
+                      << " this->getDimensions().x() = " << this->getDimensions().x()
+                      << " this->getName() = " << this->getName(); );
+      tnlAssert( edgeCoordinates.y() >= 0 && edgeCoordinates.y() < this->getDimensions.y() + 1,
+                 cerr << "edgeCoordinates.y() = " << edgeCoordinates.y()
+                      << " this->getDimensions().y() + 1 = " << this->getDimensions().y() + 1
+                      << " this->getName() = " << this->getName(); );
+      tnlAssert( edgeCoordinates.z() >= 0 && edgeCoordinates.z() < this->getDimensions.z(),
+                 cerr << "edgeCoordinates.z() = " << edgeCoordinates.z()
+                      << " this->getDimensions().z() = " << this->getDimensions().z()
+                      << " this->getName() = " << this->getName(); );
+      if( edgeCoordinates.x() == 0 || edgeCoordinates.x() == this->getDimensions().x() ||
+          edgeCoordinates.y() == 0 || edgeCoordinates.y() == this->getDimensions().y() )
+         return true;
+      return false;
+   }
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+bool tnlGrid< 3, Real, Device, Index > :: isBoundaryVertex( const CoordinatesType& vertexCoordinates ) const
+{
+   tnlAssert( vertexCoordinates.x() >= 0 && vertexCoordinates.x() < this->getDimensions().x() + 1,
+              cerr << "vertexCoordinates.x() = " << vertexCoordinates.x()
+                   << " this->getDimensions().x() + 1 = " << this->getDimensions().x() + 1
+                   << " this->getName() = " << this->getName(); );
+   tnlAssert( vertexCoordinates.y() >= 0 && vertexCoordinates.y() < this->getDimensions().y() + 1,
+              cerr << "vertexCoordinates.y() = " << vertexCoordinates.y()
+                   << " this->getDimensions().y() + 1 = " << this->getDimensions().y() + 1
+                   << " this->getName() = " << this->getName(); );
+   tnlAssert( vertexCoordinates.z() >= 0 && vertexCoordinates.z() < this->getDimensions().z() + 1,
+              cerr << "vertexCoordinates.z() = " << vertexCoordinates.z()
+                   << " this->getDimensions().z() + 1 = " << this->getDimensions().z() + 1
+                   << " this->getName() = " << this->getName(); );
+
+   if( vertexCoordinates.x() == 0 || vertexCoordinates.x() == this->getDimensions().x() ||
+       vertexCoordinates.y() == 0 || vertexCoordinates.y() == this->getDimensions().y() ||
+       vertexCoordinates.z() == 0 || vertexCoordinates.z() == this->getDimensions().z() )
+      return true;
+   return false;
+}
+
 template< typename Real,
           typename Device,
           typename Index >
diff --git a/src/implementation/mesh/tnlTraversal_Grid1D_impl.h b/src/implementation/mesh/tnlTraversal_Grid1D_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..4978fb15f9780b9022a3ba6d66b0c94284b7ad30
--- /dev/null
+++ b/src/implementation/mesh/tnlTraversal_Grid1D_impl.h
@@ -0,0 +1,114 @@
+/***************************************************************************
+                          tnlTraversal_Grid1D_impl.h  -  description
+                             -------------------
+    begin                : Jul 29, 2014
+    copyright            : (C) 2014 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef TNLTRAVERSAL_GRID1D_IMPL_H_
+#define TNLTRAVERSAL_GRID1D_IMPL_H_
+
+template< typename Real,
+          typename Index >
+   template< typename UserData,
+             typename BoundaryEntitiesProcessor,
+             typename InteriorEntitiesProcessor >
+void
+tnlTraversal< tnlGrid< 1, Real, tnlHost, Index >, 1 >::
+processEntities( const GridType& grid,
+                 UserData& userData,
+                 BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                 InteriorEntitiesProcessor& interiorEntitiesProcessor ) const
+{
+   /****
+    * Traversing cells
+    */
+   CoordinatesType coordinates;
+   const IndexType& xSize = grid.getDimensions().x();
+   for( coordinates.x() = 0; coordinates.x() < xSize; coordinates.x() ++ )
+   {
+      const IndexType index = grid.getCellIndex( coordinates );
+      if( grid.isBoundaryCell( coordinates ) )
+         boundaryEntitiesProcessor.template processEntity< 1 >( grid, userData, index, coordinates );
+      else
+         interiorEntitiesProcessor.template processEntity< 1 >( grid, userData, index, coordinates );
+   }
+}
+
+
+template< typename Real,
+          typename Index >
+   template< typename UserData,
+             typename BoundaryEntitiesProcessor,
+             typename InteriorEntitiesProcessor >
+void
+tnlTraversal< tnlGrid< 1, Real, tnlHost, Index >, 0 >::
+processEntities( const GridType& grid,
+                 UserData& userData,
+                 BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                 InteriorEntitiesProcessor& interiorEntitiesProcessor ) const
+{
+   /****
+    * Traversing vertices
+    */
+   CoordinatesType coordinates;
+   const IndexType& xSize = grid.getDimensions().x();
+   for( coordinates.x() = 0; coordinates.x() <= xSize; coordinates.x() ++ )
+   {
+      const IndexType index = grid.getVertexIndex( coordinates );
+      if( grid.isBoundaryVertex( coordinates ) )
+         boundaryEntitiesProcessor.template processEntity< 0 >( grid, userData, index, coordinates );
+      else
+         interiorEntitiesProcessor.template processEntity< 0 >( grid, userData, index, coordinates );
+   }
+}
+
+template< typename Real,
+          typename Index >
+   template< typename UserData,
+             typename BoundaryEntitiesProcessor,
+             typename InteriorEntitiesProcessor >
+void
+tnlTraversal< tnlGrid< 1, Real, tnlCuda, Index >, 1 >::
+processEntities( const GridType& grid,
+                 UserData& userData,
+                 BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                 InteriorEntitiesProcessor& interiorEntitesProcessor ) const
+{
+   /****
+    * Traversing cells
+    */
+
+}
+
+template< typename Real,
+          typename Index >
+   template< typename UserData,
+             typename BoundaryEntitiesProcessor,
+             typename InteriorEntitiesProcessor >
+void
+tnlTraversal< tnlGrid< 1, Real, tnlCuda, Index >, 0 >::
+processEntities( const GridType& grid,
+                 UserData& userData,
+                 BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                 InteriorEntitiesProcessor& interiorEntitesProcessor ) const
+{
+   /****
+    * Traversing vertices
+    */
+
+}
+
+
+
+#endif /* TNLTRAVERSAL_GRID1D_IMPL_H_ */
diff --git a/src/implementation/mesh/tnlTraversal_Grid2D_impl.h b/src/implementation/mesh/tnlTraversal_Grid2D_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..2130a6beb93adabade4fa274bc95aadffb25e573
--- /dev/null
+++ b/src/implementation/mesh/tnlTraversal_Grid2D_impl.h
@@ -0,0 +1,127 @@
+/***************************************************************************
+                          tnlTraversal_Grid2D_impl.h  -  description
+                             -------------------
+    begin                : Jul 29, 2014
+    copyright            : (C) 2014 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+
+#ifndef TNLTRAVERSAL_GRID2D_IMPL_H_
+#define TNLTRAVERSAL_GRID2D_IMPL_H_
+
+template< typename Real,
+          typename Index >
+   template< typename UserData,
+             typename BoundaryEntitiesProcessor,
+             typename InteriorEntitiesProcessor >
+void
+tnlTraversal< tnlGrid< 2, Real, tnlHost, Index >, 2 >::
+processEntities( const GridType& grid,
+                 UserData& userData,
+                 BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                 InteriorEntitiesProcessor& interiorEntitesProcessor ) const
+{
+   /****
+    * Traversing cells
+    */
+
+}
+
+template< typename Real,
+          typename Index >
+   template< typename UserData,
+             typename BoundaryEntitiesProcessor,
+             typename InteriorEntitiesProcessor >
+void
+tnlTraversal< tnlGrid< 2, Real, tnlHost, Index >, 1 >::
+processEntities( const GridType& grid,
+                 UserData& userData,
+                 BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                 InteriorEntitiesProcessor& interiorEntitesProcessor ) const
+{
+   /****
+    * Traversing faces
+    */
+}
+
+template< typename Real,
+          typename Index >
+   template< typename UserData,
+             typename BoundaryEntitiesProcessor,
+             typename InteriorEntitiesProcessor >
+void
+tnlTraversal< tnlGrid< 2, Real, tnlHost, Index >, 0 >::
+processEntities( const GridType& grid,
+                 UserData& userData,
+                 BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                 InteriorEntitiesProcessor& interiorEntitesProcessor ) const
+{
+   /****
+    * Traversing vertices
+    */
+}
+
+
+template< typename Real,
+          typename Index >
+   template< typename UserData,
+             typename BoundaryEntitiesProcessor,
+             typename InteriorEntitiesProcessor >
+void
+tnlTraversal< tnlGrid< 2, Real, tnlCuda, Index >, 2 >::
+processEntities( const GridType& grid,
+                 UserData& userData,
+                 BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                 InteriorEntitiesProcessor& interiorEntitesProcessor ) const
+{
+   /****
+    * Traversing cells
+    */
+}
+
+template< typename Real,
+          typename Index >
+   template< typename UserData,
+             typename BoundaryEntitiesProcessor,
+             typename InteriorEntitiesProcessor >
+void
+tnlTraversal< tnlGrid< 2, Real, tnlCuda, Index >, 1 >::
+processEntities( const GridType& grid,
+                 UserData& userData,
+                 BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                 InteriorEntitiesProcessor& interiorEntitesProcessor ) const
+{
+   /****
+    * Traversing faces
+    */
+}
+
+template< typename Real,
+          typename Index >
+   template< typename UserData,
+             typename BoundaryEntitiesProcessor,
+             typename InteriorEntitiesProcessor >
+void
+tnlTraversal< tnlGrid< 2, Real, tnlCuda, Index >, 0 >::
+processEntities( const GridType& grid,
+                 UserData& userData,
+                 BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                 InteriorEntitiesProcessor& interiorEntitesProcessor ) const
+{
+   /****
+    * Traversing edges
+    */
+}
+
+
+#endif /* TNLTRAVERSAL_GRID2D_IMPL_H_ */
diff --git a/src/implementation/mesh/tnlTraversal_Grid3D_impl.h b/src/implementation/mesh/tnlTraversal_Grid3D_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..a97d04a40392fd1e8f983231569148187cb066f1
--- /dev/null
+++ b/src/implementation/mesh/tnlTraversal_Grid3D_impl.h
@@ -0,0 +1,163 @@
+/***************************************************************************
+                          tnlTraversal_Grid3D_impl.h  -  description
+                             -------------------
+    begin                : Jul 29, 2014
+    copyright            : (C) 2014 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef TNLTRAVERSAL_GRID3D_IMPL_H_
+#define TNLTRAVERSAL_GRID3D_IMPL_H_
+
+template< typename Real,
+          typename Index >
+   template< typename UserData,
+             typename BoundaryEntitiesProcessor,
+             typename InteriorEntitiesProcessor >
+void
+tnlTraversal< tnlGrid< 3, Real, tnlHost, Index >, 3 >::
+processEntities( const GridType& grid,
+                 UserData& userData,
+                 BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                 InteriorEntitiesProcessor& interiorEntitesProcessor ) const
+{
+   /****
+    * Traversing cells
+    */
+}
+
+template< typename Real,
+          typename Index >
+   template< typename UserData,
+             typename BoundaryEntitiesProcessor,
+             typename InteriorEntitiesProcessor >
+void
+tnlTraversal< tnlGrid< 3, Real, tnlHost, Index >, 2 >::
+processEntities( const GridType& grid,
+                 UserData& userData,
+                 BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                 InteriorEntitiesProcessor& interiorEntitesProcessor ) const
+{
+   /****
+    * Traversing faces
+    */
+}
+
+template< typename Real,
+          typename Index >
+   template< typename UserData,
+             typename BoundaryEntitiesProcessor,
+             typename InteriorEntitiesProcessor >
+void
+tnlTraversal< tnlGrid< 3, Real, tnlHost, Index >, 1 >::
+processEntities( const GridType& grid,
+                 UserData& userData,
+                 BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                 InteriorEntitiesProcessor& interiorEntitesProcessor ) const
+{
+   /****
+    * Traversing edges
+    */
+}
+
+template< typename Real,
+          typename Index >
+   template< typename UserData,
+             typename BoundaryEntitiesProcessor,
+             typename InteriorEntitiesProcessor >
+void
+tnlTraversal< tnlGrid< 3, Real, tnlHost, Index >, 0 >::
+processEntities( const GridType& grid,
+                 UserData& userData,
+                 BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                 InteriorEntitiesProcessor& interiorEntitesProcessor ) const
+{
+   /****
+    * Traversing vertices
+    */
+}
+
+
+/****
+ * CUDA travelsal
+ */
+template< typename Real,
+          typename Index >
+   template< typename UserData,
+             typename BoundaryEntitiesProcessor,
+             typename InteriorEntitiesProcessor >
+void
+tnlTraversal< tnlGrid< 3, Real, tnlCuda, Index >, 3 >::
+processEntities( const GridType& grid,
+                 UserData& userData,
+                 BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                 InteriorEntitiesProcessor& interiorEntitesProcessor ) const
+{
+   /****
+    * Traversing cells
+    */
+}
+
+template< typename Real,
+          typename Index >
+   template< typename UserData,
+             typename BoundaryEntitiesProcessor,
+             typename InteriorEntitiesProcessor >
+void
+tnlTraversal< tnlGrid< 3, Real, tnlCuda, Index >, 2 >::
+processEntities( const GridType& grid,
+                 UserData& userData,
+                 BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                 InteriorEntitiesProcessor& interiorEntitesProcessor ) const
+{
+   /****
+    * Traversing faces
+    */
+}
+
+template< typename Real,
+          typename Index >
+   template< typename UserData,
+             typename BoundaryEntitiesProcessor,
+             typename InteriorEntitiesProcessor >
+void
+tnlTraversal< tnlGrid< 3, Real, tnlCuda, Index >, 1 >::
+processEntities( const GridType& grid,
+                 UserData& userData,
+                 BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                 InteriorEntitiesProcessor& interiorEntitesProcessor ) const
+{
+   /****
+    * Traversing edges
+    */
+}
+
+template< typename Real,
+          typename Index >
+   template< typename UserData,
+             typename BoundaryEntitiesProcessor,
+             typename InteriorEntitiesProcessor >
+void
+tnlTraversal< tnlGrid< 3, Real, tnlCuda, Index >, 0 >::
+processEntities( const GridType& grid,
+                 UserData& userData,
+                 BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                 InteriorEntitiesProcessor& interiorEntitesProcessor ) const
+{
+   /****
+    * Traversing vertices
+    */
+}
+
+
+
+#endif /* TNLTRAVERSAL_GRID3D_IMPL_H_ */
diff --git a/src/implementation/solvers/pde/CMakeLists.txt b/src/implementation/solvers/pde/CMakeLists.txt
index e4d2f3aa6cefb798c23a739d8cdce7b7339c329e..6dd0d2c0d1c569dcc05546ea6601b401a4ea2743 100755
--- a/src/implementation/solvers/pde/CMakeLists.txt
+++ b/src/implementation/solvers/pde/CMakeLists.txt
@@ -1,5 +1,6 @@
 SET( headers tnlExplicitTimeStepper_impl.h
-             tnlPDESolver_impl.h  )
+             tnlPDESolver_impl.h
+             tnlExplicitUpdater_impl.h  )
 
 INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/implementation/solvers/pde )
 
diff --git a/src/implementation/solvers/pde/tnlExplicitUpdater_impl.h b/src/implementation/solvers/pde/tnlExplicitUpdater_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..30ffc12c8d7e31747f6e61c9962296870687cae2
--- /dev/null
+++ b/src/implementation/solvers/pde/tnlExplicitUpdater_impl.h
@@ -0,0 +1,85 @@
+/***************************************************************************
+                          tnlExplicitUpdater_impl.h  -  description
+                             -------------------
+    begin                : Jul 29, 2014
+    copyright            : (C) 2014 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef TNLEXPLICITUPDATER_IMPL_H_
+#define TNLEXPLICITUPDATER_IMPL_H_
+
+#include <mesh/tnlTraversal_Grid1D.h>
+#include <mesh/tnlTraversal_Grid2D.h>
+#include <mesh/tnlTraversal_Grid3D.h>
+
+template< typename Mesh,
+          typename DofVector,
+          typename BoundaryConditions,
+          typename InteriorUpdater >
+   template< int EntityDimensions >
+void
+tnlExplicitUpdater< Mesh, DofVector, BoundaryConditions, InteriorUpdater >::
+update( const RealType& time,
+        const RealType& tau,
+        const Mesh& mesh,
+        BoundaryConditions& boundaryConditions,
+        InteriorUpdater& interiorUpdater,
+        DofVector& u,
+        DofVector& fu ) const
+{
+   TraversalUserData userData( time, tau, boundaryConditions, interiorUpdater, u, fu );
+   TraversalBoundaryEntitiesProcessor boundaryEntitiesProcessor;
+   TraversalInteriorEntitiesProcessor interiorEntitiesProcessor;
+   tnlTraversal< MeshType, EntityDimensions > meshTraversal;
+   meshTraversal.template processEntities< TraversalUserData,
+                                           TraversalBoundaryEntitiesProcessor,
+                                           TraversalInteriorEntitiesProcessor >
+                                          ( mesh,
+                                            userData,
+                                            boundaryEntitiesProcessor,
+                                            interiorEntitiesProcessor );
+}
+
+template< int Dimensions,
+          typename Real,
+          typename Device,
+          typename Index,
+          typename DofVector,
+          typename BoundaryConditions,
+          typename InteriorUpdater >
+   template< int EntityDimensions >
+void
+tnlExplicitUpdater< tnlGrid< Dimensions, Real, Device, Index >, DofVector, BoundaryConditions, InteriorUpdater >::
+update( const RealType& time,
+        const RealType& tau,
+        const tnlGrid< Dimensions, Real, Device, Index >& mesh,
+        BoundaryConditions& boundaryConditions,
+        InteriorUpdater& interiorUpdater,
+        DofVector& u,
+        DofVector& fu ) const
+{
+   TraversalUserData userData( time, tau, boundaryConditions, interiorUpdater, u, fu );
+   TraversalBoundaryEntitiesProcessor boundaryEntitiesProcessor;
+   TraversalInteriorEntitiesProcessor interiorEntitiesProcessor;
+   tnlTraversal< MeshType, EntityDimensions > meshTraversal;
+   meshTraversal.template processEntities< TraversalUserData,
+                                           TraversalBoundaryEntitiesProcessor,
+                                           TraversalInteriorEntitiesProcessor >
+                                          ( mesh,
+                                            userData,
+                                            boundaryEntitiesProcessor,
+                                            interiorEntitiesProcessor );
+}
+
+
+#endif /* TNLEXPLICITUPDATER_IMPL_H_ */
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 8574d3d0142ff29e162bf104461d618453ffcc4b..1a55833a80446a379cf3f66d5e98b01e5821ac5d 100755
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -15,6 +15,9 @@ SET( headers tnlGrid.h
              tnlMeshInitializer.h
              tnlMeshEntityInitializer.h
              tnlMeshSuperentityInitializerLayer.h
-             tnlTraversal.h )
+             tnlTraversal.h
+             tnlTraversal_Grid1D.h 
+             tnlTraversal_Grid2D.h 
+             tnlTraversal_Grid3D.h )
 
 INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/mesh )
\ No newline at end of file
diff --git a/src/mesh/tnlGrid.h b/src/mesh/tnlGrid.h
index d0a6f43630323571ca200a1917863a2e821ef758..fc5e934b0b5ef897264a738f3d1aa368a680700c 100644
--- a/src/mesh/tnlGrid.h
+++ b/src/mesh/tnlGrid.h
@@ -124,6 +124,16 @@ class tnlGrid< 1, Real, Device, Index > : public tnlObject
 #endif
    Index getNumberOfVertices() const;
 
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   bool isBoundaryCell( const CoordinatesType& cellCoordinates ) const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   bool isBoundaryVertex( const CoordinatesType& vertexCoordinates ) const;
+
    template< typename GridFunction >
    typename GridFunction::RealType getDifferenceAbsMax( const GridFunction& f1,
                                                         const GridFunction& f2 ) const;
@@ -276,6 +286,22 @@ template< int nx, int ny, typename Vertex = VertexType >
 #endif
    Index getNumberOfVertices() const;
 
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   bool isBoundaryCell( const CoordinatesType& cellCoordinates ) const;
+
+   template< int nx, int ny >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   bool isBoundaryFace( const CoordinatesType& faceCoordinates ) const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   bool isBoundaryVertex( const CoordinatesType& vertexCoordinates ) const;
+
    template< typename GridFunction >
    typename GridFunction::RealType getAbsMax( const GridFunction& f ) const;
 
@@ -456,6 +482,28 @@ template< int dx, int dy, int dz, typename Vertex = VertexType >
 #endif
    Index getNumberOfVertices() const;
 
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   bool isBoundaryCell( const CoordinatesType& cellCoordinates ) const;
+
+   template< int nx, int ny, int nz >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   bool isBoundaryFace( const CoordinatesType& faceCoordinates ) const;
+
+   template< int dx, int dy, int dz >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   bool isBoundaryEdge( const CoordinatesType& edgeCoordinates ) const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   bool isBoundaryVertex( const CoordinatesType& vertexCoordinates ) const;
+
    template< typename GridFunction >
    typename GridFunction::RealType getAbsMax( const GridFunction& f ) const;
 
diff --git a/src/mesh/tnlGridTraversal.h b/src/mesh/tnlGridTraversal.h
deleted file mode 100644
index e694ddb662b2f49420cba5272b55d5bdbc8f9918..0000000000000000000000000000000000000000
--- a/src/mesh/tnlGridTraversal.h
+++ /dev/null
@@ -1,91 +0,0 @@
-/***************************************************************************
-                          tnlGridTraversal.h  -  description
-                             -------------------
-    begin                : Jul 28, 2014
-    copyright            : (C) 2014 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/***************************************************************************
- *                                                                         *
- *   This program is free software; you can redistribute it and/or modify  *
- *   it under the terms of the GNU General Public License as published by  *
- *   the Free Software Foundation; either version 2 of the License, or     *
- *   (at your option) any later version.                                   *
- *                                                                         *
- ***************************************************************************/
-
-#ifndef TNLGRIDTRAVERSAL_H_
-#define TNLGRIDTRAVERSAL_H_
-
-#include <mesh/tnlTraversal.h>
-
-template< typename Real,
-          typename Index,
-          typename UserData >
-class tnlTraversal< tnlGrid< 1, Real, tnlHost, Index >, UserData >
-{
-   public:
-      typedef tnlGrid< 1, Real, tnlHost, Index > GridType;
-      typedef Real RealType;
-      typedef tnlHost DeviceType;
-      typedef Index IndexType;
-      typedef UserData UserDataType;
-
-      static void processBoundaryCells( const GridType& grid,
-                                        const UserDataType& userData );
-
-      static void processInteriorCells( const GridType& grid,
-                                        const UserDataType& userData );
-
-      static void processAllCells( const GridType& grid,
-                                   const UserDataType& userData );
-
-      static void processBoundaryVertices( const GridType& grid,
-                                           const UserDataType& userData );
-
-      static void processInteriorVertices( const GridType& grid,
-                                           const UserDataType& userData );
-
-      static void processAllVertices( const GridType& grid,
-                                      const UserDataType& userData );
-
-};
-
-template< typename Real,
-          typename Index,
-          typename UserData >
-class tnlTraversal< tnlGrid< 1, Real, tnlCuda, Index >, UserData >
-{
-   public:
-      typedef tnlGrid< 1, Real, tnlCuda, Index > GridType;
-      typedef Real RealType;
-      typedef tnlCuda DeviceType;
-      typedef Index IndexType;
-      typedef UserData UserDataType;
-
-      static void processBoundaryCells( const GridType& grid,
-                                        const UserDataType& userData );
-
-      static void processInteriorCells( const GridType& grid,
-                                        const UserDataType& userData );
-
-      static void processAllCells( const GridType& grid,
-                                   const UserDataType& userData );
-
-      static void processBoundaryVertices( const GridType& grid,
-                                           const UserDataType& userData );
-
-      static void processInteriorVertices( const GridType& grid,
-                                           const UserDataType& userData );
-
-      static void processAllVertices( const GridType& grid,
-                                      const UserDataType& userData );
-
-};
-
-
-
-
-
-#endif /* TNLGRIDTRAVERSAL_H_ */
diff --git a/src/mesh/tnlTraversal.h b/src/mesh/tnlTraversal.h
index a43addc3b7875b4ce599ee71f93f30420a4580b8..64a433dc4182d268d11787c8fcba904bc84bea95 100644
--- a/src/mesh/tnlTraversal.h
+++ b/src/mesh/tnlTraversal.h
@@ -19,7 +19,7 @@
 #define TNLTRAVERSAL_H_
 
 template< typename Mesh,
-          typename UserData >
+          int EntitiesDimensions = Mesh::Dimensions >
 class tnlTraversal{};
 
 
diff --git a/src/mesh/tnlTraversal_Grid1D.h b/src/mesh/tnlTraversal_Grid1D.h
new file mode 100644
index 0000000000000000000000000000000000000000..870c93460af795808eb40401fc4126ee51089ffc
--- /dev/null
+++ b/src/mesh/tnlTraversal_Grid1D.h
@@ -0,0 +1,111 @@
+/***************************************************************************
+                          tnlTraversal_Grid1D.h  -  description
+                             -------------------
+    begin                : Jul 28, 2014
+    copyright            : (C) 2014 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef TNLTRAVERSAL_GRID1D_H_
+#define TNLTRAVERSAL_GRID1D_H_
+
+#include <mesh/tnlTraversal.h>
+
+template< typename Real,
+          typename Index >
+class tnlTraversal< tnlGrid< 1, Real, tnlHost, Index >, 1 >
+{
+   public:
+      typedef tnlGrid< 1, Real, tnlHost, Index > GridType;
+      typedef Real RealType;
+      typedef tnlHost DeviceType;
+      typedef Index IndexType;
+      typedef typename GridType::CoordinatesType CoordinatesType;
+
+      template< typename UserData,
+                typename BoundaryEntitiesProcessor,
+                typename InteriorEntitiesProcessor >
+      void processEntities( const GridType& grid,
+                            UserData& userData,
+                            BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                            InteriorEntitiesProcessor& interiorEntitesProcessor ) const;
+};
+
+
+template< typename Real,
+          typename Index >
+class tnlTraversal< tnlGrid< 1, Real, tnlHost, Index >, 0 >
+{
+   public:
+      typedef tnlGrid< 1, Real, tnlHost, Index > GridType;
+      typedef Real RealType;
+      typedef tnlHost DeviceType;
+      typedef Index IndexType;
+      typedef typename GridType::CoordinatesType CoordinatesType;
+
+      template< typename UserData,
+                typename BoundaryEntitiesProcessor,
+                typename InteriorEntitiesProcessor >
+      void processEntities( const GridType& grid,
+                            UserData& userData,
+                            BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                            InteriorEntitiesProcessor& interiorEntitesProcessor ) const;
+};
+
+/****
+ * CUDA traversals
+ */
+template< typename Real,
+          typename Index >
+class tnlTraversal< tnlGrid< 1, Real, tnlCuda, Index >, 1 >
+{
+   public:
+      typedef tnlGrid< 1, Real, tnlCuda, Index > GridType;
+      typedef Real RealType;
+      typedef tnlCuda DeviceType;
+      typedef Index IndexType;
+      typedef typename GridType::CoordinatesType CoordinatesType;
+
+      template< typename UserData,
+                typename BoundaryEntitiesProcessor,
+                typename InteriorEntitiesProcessor >
+      void processEntities( const GridType& grid,
+                            UserData& userData,
+                            BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                            InteriorEntitiesProcessor& interiorEntitesProcessor ) const;
+
+};
+
+template< typename Real,
+          typename Index >
+class tnlTraversal< tnlGrid< 1, Real, tnlCuda, Index >, 0 >
+{
+   public:
+      typedef tnlGrid< 1, Real, tnlCuda, Index > GridType;
+      typedef Real RealType;
+      typedef tnlCuda DeviceType;
+      typedef Index IndexType;
+      typedef typename GridType::CoordinatesType CoordinatesType;
+
+      template< typename UserData,
+                typename BoundaryEntitiesProcessor,
+                typename InteriorEntitiesProcessor >
+      void processEntities( const GridType& grid,
+                            UserData& userData,
+                            BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                            InteriorEntitiesProcessor& interiorEntitesProcessor ) const;
+
+};
+
+#include <implementation/mesh/tnlTraversal_Grid1D_impl.h>
+
+#endif /* TNLTRAVERSAL_GRID1D_H_ */
diff --git a/src/mesh/tnlTraversal_Grid2D.h b/src/mesh/tnlTraversal_Grid2D.h
new file mode 100644
index 0000000000000000000000000000000000000000..33d4a7468220a77c241fbcf90ee30508abb148fe
--- /dev/null
+++ b/src/mesh/tnlTraversal_Grid2D.h
@@ -0,0 +1,150 @@
+/***************************************************************************
+                          tnlTraversal_Grid2D.h  -  description
+                             -------------------
+    begin                : Jul 29, 2014
+    copyright            : (C) 2014 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+
+#ifndef TNLTRAVERSAL_GRID2D_H_
+#define TNLTRAVERSAL_GRID2D_H_
+
+
+template< typename Real,
+          typename Index >
+class tnlTraversal< tnlGrid< 2, Real, tnlHost, Index >, 2 >
+{
+   public:
+      typedef tnlGrid< 2, Real, tnlHost, Index > GridType;
+      typedef Real RealType;
+      typedef tnlHost DeviceType;
+      typedef Index IndexType;
+      typedef typename GridType::CoordinatesType CoordinatesType;
+
+      template< typename UserData,
+                typename BoundaryEntitiesProcessor,
+                typename InteriorEntitiesProcessor >
+      void processEntities( const GridType& grid,
+                            UserData& userData,
+                            BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                            InteriorEntitiesProcessor& interiorEntitesProcessor ) const;
+};
+
+template< typename Real,
+          typename Index >
+class tnlTraversal< tnlGrid< 2, Real, tnlHost, Index >, 1 >
+{
+   public:
+      typedef tnlGrid< 2, Real, tnlHost, Index > GridType;
+      typedef Real RealType;
+      typedef tnlHost DeviceType;
+      typedef Index IndexType;
+      typedef typename GridType::CoordinatesType CoordinatesType;
+
+      template< typename UserData,
+                typename BoundaryEntitiesProcessor,
+                typename InteriorEntitiesProcessor >
+      void processEntities( const GridType& grid,
+                            UserData& userData,
+                            BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                            InteriorEntitiesProcessor& interiorEntitesProcessor ) const;
+};
+
+template< typename Real,
+          typename Index >
+class tnlTraversal< tnlGrid< 2, Real, tnlHost, Index >, 0 >
+{
+   public:
+      typedef tnlGrid< 2, Real, tnlHost, Index > GridType;
+      typedef Real RealType;
+      typedef tnlHost DeviceType;
+      typedef Index IndexType;
+      typedef typename GridType::CoordinatesType CoordinatesType;
+
+      template< typename UserData,
+                typename BoundaryEntitiesProcessor,
+                typename InteriorEntitiesProcessor >
+      void processEntities( const GridType& grid,
+                            UserData& userData,
+                            BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                            InteriorEntitiesProcessor& interiorEntitesProcessor ) const;
+};
+
+
+template< typename Real,
+          typename Index >
+class tnlTraversal< tnlGrid< 2, Real, tnlCuda, Index >, 2 >
+{
+   public:
+      typedef tnlGrid< 2, Real, tnlCuda, Index > GridType;
+      typedef Real RealType;
+      typedef tnlCuda DeviceType;
+      typedef Index IndexType;
+      typedef typename GridType::CoordinatesType CoordinatesType;
+
+      template< typename UserData,
+                typename BoundaryEntitiesProcessor,
+                typename InteriorEntitiesProcessor >
+      void processEntities( const GridType& grid,
+                            UserData& userData,
+                            BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                            InteriorEntitiesProcessor& interiorEntitesProcessor ) const;
+
+};
+
+template< typename Real,
+          typename Index >
+class tnlTraversal< tnlGrid< 2, Real, tnlCuda, Index >, 1 >
+{
+   public:
+      typedef tnlGrid< 2, Real, tnlCuda, Index > GridType;
+      typedef Real RealType;
+      typedef tnlCuda DeviceType;
+      typedef Index IndexType;
+      typedef typename GridType::CoordinatesType CoordinatesType;
+
+      template< typename UserData,
+                typename BoundaryEntitiesProcessor,
+                typename InteriorEntitiesProcessor >
+      void processEntities( const GridType& grid,
+                            UserData& userData,
+                            BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                            InteriorEntitiesProcessor& interiorEntitesProcessor ) const;
+
+};
+
+template< typename Real,
+          typename Index >
+class tnlTraversal< tnlGrid< 2, Real, tnlCuda, Index >, 0 >
+{
+   public:
+      typedef tnlGrid< 2, Real, tnlCuda, Index > GridType;
+      typedef Real RealType;
+      typedef tnlCuda DeviceType;
+      typedef Index IndexType;
+      typedef typename GridType::CoordinatesType CoordinatesType;
+
+      template< typename UserData,
+                typename BoundaryEntitiesProcessor,
+                typename InteriorEntitiesProcessor >
+      void processEntities( const GridType& grid,
+                            UserData& userData,
+                            BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                            InteriorEntitiesProcessor& interiorEntitesProcessor ) const;
+
+};
+
+
+#include <implementation/mesh/tnlTraversal_Grid2D_impl.h>
+
+#endif /* TNLTRAVERSAL_GRID2D_H_ */
diff --git a/src/mesh/tnlTraversal_Grid3D.h b/src/mesh/tnlTraversal_Grid3D.h
new file mode 100644
index 0000000000000000000000000000000000000000..b2bade1a3bbab3df8803203bec32ef64cbbd7c2b
--- /dev/null
+++ b/src/mesh/tnlTraversal_Grid3D.h
@@ -0,0 +1,193 @@
+/***************************************************************************
+                          tnlTraversal_Grid3D.h  -  description
+                             -------------------
+    begin                : Jul 29, 2014
+    copyright            : (C) 2014 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef TNLTRAVERSAL_GRID3D_H_
+#define TNLTRAVERSAL_GRID3D_H_
+
+
+template< typename Real,
+          typename Index >
+class tnlTraversal< tnlGrid< 3, Real, tnlHost, Index >, 3 >
+{
+   public:
+      typedef tnlGrid< 3, Real, tnlHost, Index > GridType;
+      typedef Real RealType;
+      typedef tnlHost DeviceType;
+      typedef Index IndexType;
+      typedef typename GridType::CoordinatesType CoordinatesType;
+
+      template< typename UserData,
+                typename BoundaryEntitiesProcessor,
+                typename InteriorEntitiesProcessor >
+      void processEntities( const GridType& grid,
+                            UserData& userData,
+                            BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                            InteriorEntitiesProcessor& interiorEntitesProcessor ) const;
+};
+
+template< typename Real,
+          typename Index >
+class tnlTraversal< tnlGrid< 3, Real, tnlHost, Index >, 2 >
+{
+   public:
+      typedef tnlGrid< 3, Real, tnlHost, Index > GridType;
+      typedef Real RealType;
+      typedef tnlHost DeviceType;
+      typedef Index IndexType;
+      typedef typename GridType::CoordinatesType CoordinatesType;
+
+      template< typename UserData,
+                typename BoundaryEntitiesProcessor,
+                typename InteriorEntitiesProcessor >
+      void processEntities( const GridType& grid,
+                            UserData& userData,
+                            BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                            InteriorEntitiesProcessor& interiorEntitesProcessor ) const;
+};
+
+template< typename Real,
+          typename Index >
+class tnlTraversal< tnlGrid< 3, Real, tnlHost, Index >, 1 >
+{
+   public:
+      typedef tnlGrid< 3, Real, tnlHost, Index > GridType;
+      typedef Real RealType;
+      typedef tnlHost DeviceType;
+      typedef Index IndexType;
+      typedef typename GridType::CoordinatesType CoordinatesType;
+
+      template< typename UserData,
+                typename BoundaryEntitiesProcessor,
+                typename InteriorEntitiesProcessor >
+      void processEntities( const GridType& grid,
+                            UserData& userData,
+                            BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                            InteriorEntitiesProcessor& interiorEntitesProcessor ) const;
+};
+
+template< typename Real,
+          typename Index >
+class tnlTraversal< tnlGrid< 3, Real, tnlHost, Index >, 0 >
+{
+   public:
+      typedef tnlGrid< 3, Real, tnlHost, Index > GridType;
+      typedef Real RealType;
+      typedef tnlHost DeviceType;
+      typedef Index IndexType;
+      typedef typename GridType::CoordinatesType CoordinatesType;
+
+      template< typename UserData,
+                typename BoundaryEntitiesProcessor,
+                typename InteriorEntitiesProcessor >
+      void processEntities( const GridType& grid,
+                            UserData& userData,
+                            BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                            InteriorEntitiesProcessor& interiorEntitesProcessor ) const;
+};
+
+/****
+ * CUDA traversal
+ */
+
+template< typename Real,
+          typename Index >
+class tnlTraversal< tnlGrid< 3, Real, tnlCuda, Index >, 3 >
+{
+   public:
+      typedef tnlGrid< 3, Real, tnlCuda, Index > GridType;
+      typedef Real RealType;
+      typedef tnlCuda DeviceType;
+      typedef Index IndexType;
+      typedef typename GridType::CoordinatesType CoordinatesType;
+
+      template< typename UserData,
+                typename BoundaryEntitiesProcessor,
+                typename InteriorEntitiesProcessor >
+      void processEntities( const GridType& grid,
+                            UserData& userData,
+                            BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                            InteriorEntitiesProcessor& interiorEntitesProcessor ) const;
+
+};
+
+template< typename Real,
+          typename Index >
+class tnlTraversal< tnlGrid< 3, Real, tnlCuda, Index >, 2 >
+{
+   public:
+      typedef tnlGrid< 3, Real, tnlCuda, Index > GridType;
+      typedef Real RealType;
+      typedef tnlCuda DeviceType;
+      typedef Index IndexType;
+      typedef typename GridType::CoordinatesType CoordinatesType;
+
+      template< typename UserData,
+                typename BoundaryEntitiesProcessor,
+                typename InteriorEntitiesProcessor >
+      void processEntities( const GridType& grid,
+                            UserData& userData,
+                            BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                            InteriorEntitiesProcessor& interiorEntitesProcessor ) const;
+
+};
+
+template< typename Real,
+          typename Index >
+class tnlTraversal< tnlGrid< 3, Real, tnlCuda, Index >, 1 >
+{
+   public:
+      typedef tnlGrid< 3, Real, tnlCuda, Index > GridType;
+      typedef Real RealType;
+      typedef tnlCuda DeviceType;
+      typedef Index IndexType;
+      typedef typename GridType::CoordinatesType CoordinatesType;
+
+      template< typename UserData,
+                typename BoundaryEntitiesProcessor,
+                typename InteriorEntitiesProcessor >
+      void processEntities( const GridType& grid,
+                            UserData& userData,
+                            BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                            InteriorEntitiesProcessor& interiorEntitesProcessor ) const;
+
+};
+
+template< typename Real,
+          typename Index >
+class tnlTraversal< tnlGrid< 3, Real, tnlCuda, Index >, 0 >
+{
+   public:
+      typedef tnlGrid< 3, Real, tnlCuda, Index > GridType;
+      typedef Real RealType;
+      typedef tnlCuda DeviceType;
+      typedef Index IndexType;
+      typedef typename GridType::CoordinatesType CoordinatesType;
+
+      template< typename UserData,
+                typename BoundaryEntitiesProcessor,
+                typename InteriorEntitiesProcessor >
+      void processEntities( const GridType& grid,
+                            UserData& userData,
+                            BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
+                            InteriorEntitiesProcessor& interiorEntitesProcessor ) const;
+
+};
+
+
+#include <implementation/mesh/tnlTraversal_Grid3D_impl.h>
+
+#endif /* TNLTRAVERSAL_GRID3D_H_ */
diff --git a/src/solvers/pde/CMakeLists.txt b/src/solvers/pde/CMakeLists.txt
index 2304b7a455881495b102f14738c3aff1120cc5c0..3c3da45fef7d0f4aa2efe17df89e1333c26f945f 100755
--- a/src/solvers/pde/CMakeLists.txt
+++ b/src/solvers/pde/CMakeLists.txt
@@ -1,5 +1,6 @@
 SET( headers tnlPDESolver.h
              tnlExplicitTimeStepper.h
+             tnlExplicitUpdater.h
    )
    
 INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/solvers/pde )
\ No newline at end of file
diff --git a/src/solvers/pde/tnlExplicitUpdater.h b/src/solvers/pde/tnlExplicitUpdater.h
new file mode 100644
index 0000000000000000000000000000000000000000..0a04ef81252cd608b689049645927655c5606fb1
--- /dev/null
+++ b/src/solvers/pde/tnlExplicitUpdater.h
@@ -0,0 +1,194 @@
+/***************************************************************************
+                          tnlExplicitUpdater.h  -  description
+                             -------------------
+    begin                : Jul 29, 2014
+    copyright            : (C) 2014 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef TNLEXPLICITUPDATER_H_
+#define TNLEXPLICITUPDATER_H_
+
+template< typename Real,
+          typename DofVector,
+          typename BoundaryConditions,
+          typename InteriorUpdater >
+class tnlExplicitUpdaterTraversalUserData
+{
+   public:
+      const Real &time, &tau;
+
+      BoundaryConditions& boundaryConditions;
+
+      InteriorUpdater& interiorUpdater;
+
+      DofVector &u, &fu;
+
+      tnlExplicitUpdaterTraversalUserData( const Real& time,
+                                           const Real& tau,
+                                           BoundaryConditions& boundaryConditions,
+                                           InteriorUpdater& interiorUpdater,
+                                           DofVector& u,
+                                           DofVector& fu )
+      : time( time ),
+        tau( tau ),
+        boundaryConditions( boundaryConditions ),
+        interiorUpdater( interiorUpdater ),
+        u( u ),
+        fu( fu )
+      {};
+};
+
+template< typename Mesh,
+          typename DofVector,
+          typename BoundaryConditions,
+          typename InteriorUpdater >
+class tnlExplicitUpdater
+{
+   public:
+      typedef Mesh MeshType;
+      typedef typename DofVector::RealType RealType;
+      typedef typename DofVector::DeviceType DeviceType;
+      typedef typename DofVector::IndexType IndexType;
+      typedef tnlExplicitUpdaterTraversalUserData< RealType,
+                                                   DofVector,
+                                                   BoundaryConditions,
+                                                   InteriorUpdater > TraversalUserData;
+
+      template< int EntityDimensions >
+      void update( const RealType& time,
+                   const RealType& tau,
+                   const MeshType& mesh,
+                   BoundaryConditions& boundaryConditions,
+                   InteriorUpdater& interiorUpdater,
+                   DofVector& u,
+                   DofVector& fu ) const;
+
+      class TraversalBoundaryEntitiesProcessor
+      {
+         public:
+
+            template< int EntityDimension >
+            void processEntity( const MeshType& mesh,
+                                TraversalUserData& userData,
+                                const IndexType index )
+            {
+               userData.boundaryConditions.setBoundaryConditions( userData.time,
+                                                                  userData.tau,
+                                                                  mesh,
+                                                                  index,
+                                                                  userData.u,
+                                                                  userData.fu );
+            }
+
+      };
+
+      class TraversalInteriorEntitiesProcessor
+      {
+         public:
+
+            template< int EntityDimension >
+            void processEntity( const MeshType& mesh,
+                                TraversalUserData& userData,
+                                const IndexType index )
+            {
+               userData.boundaryConditions.update( userData.time,
+                                                   userData.tau,
+                                                   mesh,
+                                                   index,
+                                                   userData.u,
+                                                   userData.fu );
+            }
+
+      };
+
+};
+
+template< int Dimensions,
+          typename Real,
+          typename Device,
+          typename Index,
+          typename DofVector,
+          typename BoundaryConditions,
+          typename InteriorUpdater >
+class tnlExplicitUpdater< tnlGrid< Dimensions, Real, Device, Index >,
+                          DofVector,
+                          BoundaryConditions,
+                          InteriorUpdater >
+{
+   public:
+
+      typedef tnlGrid< Dimensions, Real, Device, Index > MeshType;
+      typedef typename MeshType::RealType RealType;
+      typedef typename MeshType::DeviceType DeviceType;
+      typedef typename MeshType::IndexType IndexType;
+      typedef typename MeshType::CoordinatesType CoordinatesType;
+      typedef tnlExplicitUpdaterTraversalUserData< RealType,
+                                                   DofVector,
+                                                   BoundaryConditions,
+                                                   InteriorUpdater > TraversalUserData;
+      
+      template< int EntityDimensions >
+      void update( const RealType& time,
+                   const RealType& tau,
+                   const MeshType& mesh,
+                   BoundaryConditions& boundaryConditions,
+                   InteriorUpdater& interiorUpdater,
+                   DofVector& u,
+                   DofVector& fu ) const;
+
+      class TraversalBoundaryEntitiesProcessor
+      {
+         public:
+
+            template< int EntityDimension >
+            void processEntity( const MeshType& mesh,
+                                TraversalUserData& userData,
+                                const IndexType index,
+                                const CoordinatesType& coordinates )
+            {
+               userData.boundaryConditions.setBoundaryConditions( userData.time,
+                                                                  userData.tau,
+                                                                  mesh,
+                                                                  index,
+                                                                  coordinates,
+                                                                  userData.u,
+                                                                  userData.fu );
+            }
+
+      };
+
+      class TraversalInteriorEntitiesProcessor
+      {
+         public:
+
+            template< int EntityDimension >
+            void processEntity( const MeshType& mesh,
+                                TraversalUserData& userData,
+                                const IndexType index,
+                                const CoordinatesType& coordinates )
+            {
+               userData.interiorUpdater.update( userData.time,
+                                                userData.tau,
+                                                mesh,
+                                                index,
+                                                coordinates,
+                                                userData.u,
+                                                userData.fu );
+            }
+
+      };
+
+};
+
+#include <implementation/solvers/pde/tnlExplicitUpdater_impl.h>
+#endif /* TNLEXPLICITUPDATER_H_ */