diff --git a/examples/heat-equation/heatEquationSolver_impl.h b/examples/heat-equation/heatEquationSolver_impl.h
index 6b6f2ae774f0505e29ffd5aaa9698470b9b1cf90..2ad086e8e7bfd53501fb31875bff8099fcb5de33 100644
--- a/examples/heat-equation/heatEquationSolver_impl.h
+++ b/examples/heat-equation/heatEquationSolver_impl.h
@@ -144,7 +144,7 @@ bool heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide,TimeFunc
       return false;
    }
    
-   boundaryCondition.applyBoundaryConditions(mesh,numericalSolution,0.0,timeFunction,analyticSpaceFunction);
+   //boundaryCondition.applyBoundaryConditions(mesh,numericalSolution,0.0,timeFunction,analyticSpaceFunction);
    timeFunction.applyInitTimeValues( numericalSolution);
    
    return true;
@@ -184,7 +184,7 @@ makeSnapshot( const RealType& time,
    if(ifLaplaceCompare == 1)
    {
       analyticSolution.computeLaplace( mesh, time, analyticLaplace, timeFunction, analyticSpaceFunction );
-      diffusion.getExplicitRHS( mesh, numericalSolution, numericalLaplace );
+      //diffusion.getExplicitRHS( mesh, numericalSolution, numericalLaplace );
       
       tnlString fileName;
       FileNameBaseNumberEnding( "analyticLaplace", 0, 1, ".tnl", fileName );
diff --git a/examples/heat-equation/tnlDirichletBoundaryConditions.h b/examples/heat-equation/tnlDirichletBoundaryConditions.h
index 926a96cfd976be59f7fc0581379f8426f01c3fde..8cde1f213aad5651bc9d09c03758f77c5774100c 100644
--- a/examples/heat-equation/tnlDirichletBoundaryConditions.h
+++ b/examples/heat-equation/tnlDirichletBoundaryConditions.h
@@ -37,28 +37,6 @@ class tnlDirichletBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >,
                                const CoordinatesType& coordinates,
                                DofVectorType& u,
                                DofVectorType& fu );
-
-   //--------------------------------
-
-   template< typename AnalyticSpaceFunction,
-             typename TimeFunction,
-             typename Vector >
-   void applyBoundaryConditions( const MeshType& mesh,
-                                 Vector& u,
-                                 const RealType& time,
-                                 TimeFunction& timeFunction,
-                                 AnalyticSpaceFunction& analyticSpaceFunction );
-   
-   template< typename AnalyticSpaceFunction,
-             typename TimeFunction,
-             typename Vector >
-   void applyBoundaryTimeDerivation( const MeshType& mesh,
-                                     Vector& u,
-                                     const RealType& time,
-                                     TimeFunction& timeFunction,
-                                     AnalyticSpaceFunction& analyticSpaceFunction );
-   //-----------------------------------
-
 };
 
 template< typename MeshReal,
@@ -87,27 +65,6 @@ class tnlDirichletBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >,
                                const CoordinatesType& coordinates,
                                DofVectorType& u,
                                DofVectorType& fu );
-
-
-   //----------------------------------------
-   template< typename AnalyticSpaceFunction,
-             typename TimeFunction,
-             typename Vector >
-   void applyBoundaryConditions( const MeshType& mesh,
-                                 Vector& u,
-                                 const RealType& time,
-                                 TimeFunction& timeFunction,
-                                 AnalyticSpaceFunction& analyticSpaceFunction );
-
-   template< typename AnalyticSpaceFunction,
-             typename TimeFunction,
-             typename Vector >
-   void applyBoundaryTimeDerivation( const MeshType& mesh,
-                                     Vector& u,
-                                     const RealType& time,
-                                     TimeFunction& timeFunction,
-                                     AnalyticSpaceFunction& analyticSpaceFunction );
-   //----------------------------------
 };
 
 template< typename MeshReal,
@@ -137,25 +94,6 @@ class tnlDirichletBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >,
                                DofVectorType& u,
                                DofVectorType& fu );
 
-   //--------------------------------------------------
-   template< typename AnalyticSpaceFunction,
-             typename TimeFunction,
-             typename Vector >
-   void applyBoundaryConditions( const MeshType& mesh,
-                                 Vector& u,
-                                 const RealType& time,
-                                 TimeFunction& timeFunction,
-                                 AnalyticSpaceFunction& analyticSpaceFunction );
-   
-   template< typename AnalyticSpaceFunction,
-             typename TimeFunction,
-             typename Vector >
-   void applyBoundaryTimeDerivation( const MeshType& mesh,
-                                     Vector& u,
-                                     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 18ea0fd52c44d0d32f3429cd40a24a4e00c6dc30..662fbe1e933634af481042cba13ac42bb1d66200 100644
--- a/examples/heat-equation/tnlDirichletBoundaryConditions_impl.h
+++ b/examples/heat-equation/tnlDirichletBoundaryConditions_impl.h
@@ -22,70 +22,6 @@ setBoundaryConditions( const RealType& time,
    u[ index ] = 0;
 }
 
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          typename Index >
-   template< typename AnalyticSpaceFunction,
-             typename TimeFunction,
-             typename Vector >
-void tnlDirichletBoundaryConditions< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index >::
-applyBoundaryConditions( const MeshType& mesh,
-                         Vector& u,
-                         const RealType& time,
-                         TimeFunction& timeFunction,
-                         AnalyticSpaceFunction& analyticSpaceFunction )
-{
-   RealType timeFunctionValue = timeFunction.getTimeValue(time);
-   
-   CoordinatesType coordinates;
-   VertexType vertex;
-   
-   coordinates = mesh.getCellCoordinates( 0 );
-   vertex = mesh.getCellCenter( coordinates );
-   
-   u[0] = timeFunctionValue*analyticSpaceFunction.getF(vertex);
-   
-   coordinates = mesh.getCellCoordinates( u.getSize() - 1 );
-   vertex = mesh.getCellCenter( coordinates );
-
-   u[u.getSize()-1] = timeFunctionValue*analyticSpaceFunction.getF(vertex);
-}
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          typename Index >
-   template< typename AnalyticSpaceFunction,
-             typename TimeFunction,
-             typename Vector >
-void tnlDirichletBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
-applyBoundaryTimeDerivation( const MeshType& mesh,
-                             Vector& u,
-                             const RealType& time,
-                             TimeFunction& timeFunction,
-                             AnalyticSpaceFunction& analyticSpaceFunction )
-{
-   RealType timeFunctionDerivationValue = timeFunction.getDerivation(time);
-   
-   CoordinatesType coordinates;
-   VertexType vertex;
-   
-   coordinates = mesh.getCellCoordinates( 0 );
-   vertex = mesh.getCellCenter( coordinates );
-   
-   u[0] = timeFunctionDerivationValue*analyticSpaceFunction.getF(vertex);
-   
-   coordinates = mesh.getCellCoordinates( u.getSize()-1 );
-   vertex = mesh.getCellCenter( coordinates );
-
-   u[u.getSize()-1] = timeFunctionDerivationValue*analyticSpaceFunction.getF(vertex);
-}
-
-
 template< typename MeshReal,
           typename Device,
           typename MeshIndex,
@@ -105,118 +41,6 @@ setBoundaryConditions( const RealType& time,
    u[ index ] = 0;
 }
 
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          typename Index >
-   template< typename AnalyticSpaceFunction,
-             typename TimeFunction,
-             typename Vector >
-void tnlDirichletBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
-applyBoundaryConditions( const MeshType& mesh,
-                         Vector& u,
-                         const RealType& time,
-                         TimeFunction& timeFunction,
-                         AnalyticSpaceFunction& analyticSpaceFunction )
-{  
-   RealType timeFunctionValue = timeFunction.getTimeValue(time);   
-   
-   CoordinatesType dimensions = mesh.getDimensions();
-   CoordinatesType coordinates1;
-   CoordinatesType coordinates2;
-   VertexType vertex;
-   
-   coordinates1.y()=0;
-   coordinates2.y()=dimensions.y()-1;
-
-   #ifdef HAVE_OPENMP
-    #pragma omp parallel for firstprivate(coordinates1,coordinates2) private(vertex)
-   #endif
-   for(IndexType i=0; i<dimensions.x(); i++)
-   {
-      coordinates1.x()=coordinates2.x()=i;
-      vertex = mesh.getCellCenter( coordinates1 );
-      u[i] = timeFunctionValue*analyticSpaceFunction.getF(vertex);
-      
-      vertex = mesh.getCellCenter( coordinates2 );
-      u[(dimensions.y()-1)*mesh.getDimensions().x()+i] = timeFunctionValue*analyticSpaceFunction.getF(vertex);
-   }
-   
-   coordinates1.x()=0;
-   coordinates2.x()=dimensions.x()-1;
-
-   #ifdef HAVE_OPENMP
-    #pragma omp parallel for firstprivate(coordinates1,coordinates2) private(vertex)
-   #endif
-   for(IndexType i=0; i<dimensions.y(); i++)
-   {
-      coordinates1.y()=coordinates2.y()=i;
-           
-      vertex = mesh.getCellCenter( coordinates1 );
-      u[i*mesh.getDimensions().x()] = timeFunctionValue*analyticSpaceFunction.getF(vertex);
-      
-      vertex = mesh.getCellCenter( coordinates2 );
-      u[i*mesh.getDimensions().x()+dimensions.x()-1] = timeFunctionValue*analyticSpaceFunction.getF(vertex);
-   }
-}
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          typename Index >
-   template< typename AnalyticSpaceFunction,
-             typename TimeFunction,
-             typename Vector >
-void tnlDirichletBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
-applyBoundaryTimeDerivation( const MeshType& mesh,
-                             Vector& u,
-                             const RealType& time,
-                             TimeFunction& timeFunction,
-                             AnalyticSpaceFunction& analyticSpaceFunction )
-{
-   RealType timeFunctionDerivationValue = timeFunction.getDerivation(time);
-   
-   CoordinatesType dimensions = mesh.getDimensions();
-   CoordinatesType coordinates1;
-   CoordinatesType coordinates2;
-   VertexType vertex;
-
-   coordinates1.y()=0;
-   coordinates2.y()=dimensions.y()-1;  
-
-   #ifdef HAVE_OPENMP
-    #pragma omp parallel for firstprivate(coordinates1,coordinates2) private(vertex)
-   #endif
-   for(IndexType i=0; i<dimensions.x(); i++)
-   {
-      coordinates1.x()=coordinates2.x()=i;
-      vertex = mesh.getCellCenter( coordinates1 );
-      u[mesh.getCellIndex( CoordinatesType( i,0) )] = timeFunctionDerivationValue*analyticSpaceFunction.getF(vertex);
-      
-      vertex = mesh.getCellCenter( coordinates2 );
-      u[mesh.getCellIndex( CoordinatesType( i,dimensions.y()-1) )] = timeFunctionDerivationValue*analyticSpaceFunction.getF(vertex);
-   }
-   
-   coordinates1.x()=0;
-   coordinates2.x()=dimensions.x()-1; 
-
-   #ifdef HAVE_OPENMP
-    #pragma omp parallel for firstprivate(coordinates1,coordinates2) private(vertex)
-   #endif
-   for(IndexType i=0;i < dimensions.y(); i++)
-   {
-      coordinates1.y()=coordinates2.y()=i;
-           
-      vertex = mesh.getCellCenter( coordinates1 );
-      u[mesh.getCellIndex( CoordinatesType( 0,i) )] = timeFunctionDerivationValue*analyticSpaceFunction.getF(vertex);
-
-      vertex = mesh.getCellCenter( coordinates2 );
-      u[mesh.getCellIndex( CoordinatesType( dimensions.x()-1,i) )] = timeFunctionDerivationValue*analyticSpaceFunction.getF(vertex);
-   }
-}
-
 template< typename MeshReal,
           typename Device,
           typename MeshIndex,
@@ -236,170 +60,6 @@ setBoundaryConditions( const RealType& time,
    u[ index ] = 0;
 }
 
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          typename Index >
-   template< typename AnalyticSpaceFunction,
-             typename TimeFunction,
-             typename Vector >
-void tnlDirichletBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
-applyBoundaryConditions( const MeshType& mesh,
-                         Vector& u,
-                         const RealType& time,
-                         TimeFunction& timeFunction,
-                         AnalyticSpaceFunction& analyticSpaceFunction )
-{
-   RealType timeFunctionValue = timeFunction.getTimeValue(time);
-   
-   CoordinatesType dimensions = mesh.getDimensions();
-   CoordinatesType coordinates;
-   VertexType vertex;
-
-   #ifdef HAVE_OPENMP
-    #pragma omp parallel for firstprivate(coordinates) private(vertex)
-   #endif
-   for(IndexType i=0; i<dimensions.x(); i++)
-   {
-      for(IndexType j=0; j<dimensions.y(); j++)
-      {
-         coordinates.x()=i;
-         coordinates.y()=j;
-         
-         coordinates.z()=0;
-         vertex = mesh.getCellCenter( coordinates );
-         u[mesh.getCellIndex( CoordinatesType( i,j,0 ) )] = timeFunctionValue*analyticSpaceFunction.getF(vertex);
-         
-         coordinates.z()=dimensions.z()-1;
-         vertex = mesh.getCellCenter( coordinates );
-         u[mesh.getCellIndex( CoordinatesType( i,j,dimensions.z()-1) ) ] = timeFunctionValue*analyticSpaceFunction.getF(vertex);
-      }
-   }
-
-   #ifdef HAVE_OPENMP
-    #pragma omp parallel for firstprivate(coordinates) private(vertex)
-   #endif
-   for(IndexType i=0; i<dimensions.x(); i++)
-   {
-      for(IndexType j=0; j<dimensions.y(); j++)
-      {
-         coordinates.y()=i;
-         coordinates.z()=j;
-         
-         coordinates.x()=0;
-         vertex = mesh.getCellCenter( coordinates );
-         u[mesh.getCellIndex( CoordinatesType( 0,i,j ) ) ] = timeFunctionValue*analyticSpaceFunction.getF(vertex);
-         
-         coordinates.x()=dimensions.x()-1;
-         vertex = mesh.getCellCenter( coordinates );
-         u[mesh.getCellIndex( CoordinatesType( dimensions.x()-1,i,j ) ) ] = timeFunctionValue*analyticSpaceFunction.getF(vertex);
-      }
-   }
-
-   #ifdef HAVE_OPENMP
-    #pragma omp parallel for firstprivate(coordinates) private(vertex)
-   #endif
-   for(IndexType i=0; i<dimensions.x(); i++)
-   {
-      for(IndexType j=1; j<dimensions.y(); j++)
-      {
-         coordinates.x()=i;
-         coordinates.z()=j;
-         
-         coordinates.y()=0;
-         vertex = mesh.getCellCenter( coordinates );         
-         u[ mesh.getCellIndex( CoordinatesType( i,0,j ) ) ] = timeFunctionValue*analyticSpaceFunction.getF(vertex);
-         
-         coordinates.y()=dimensions.y()-1;
-         vertex = mesh.getCellCenter( coordinates );   
-         u[mesh.getCellIndex( CoordinatesType( i,dimensions.y()-1,j ) ) ] = timeFunctionValue*analyticSpaceFunction.getF(vertex);
-      }
-   } 
-}
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          typename Index >
-   template< typename AnalyticSpaceFunction,
-             typename TimeFunction,
-             typename Vector >
-void tnlDirichletBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
-applyBoundaryTimeDerivation( const MeshType& mesh,
-                             Vector& u,
-                             const RealType& time,
-                             TimeFunction& timeFunction,
-                             AnalyticSpaceFunction& analyticSpaceFunction )
-{
-   RealType timeFunctionDerivationValue = timeFunction.getDerivation(time);
-   
-   CoordinatesType dimensions = mesh.getDimensions();
-   CoordinatesType coordinates;
-   VertexType vertex;
-
-   #ifdef HAVE_OPENMP
-    #pragma omp parallel for firstprivate(coordinates) private(vertex)
-   #endif
-   for(IndexType i=0; i<dimensions.x(); i++)
-   {
-      for(IndexType j=0; j<dimensions.y(); j++)
-      {
-         coordinates.x()=i;
-         coordinates.y()=j;
-         
-         coordinates.z()=0;
-         vertex = mesh.getCellCenter( coordinates );
-         u[mesh.getCellIndex( CoordinatesType( i,j,0 ) ) ] = timeFunctionDerivationValue*analyticSpaceFunction.getF(vertex);
-         
-         coordinates.z()=dimensions.z()-1;
-         vertex = mesh.getCellCenter( coordinates );
-         u[mesh.getCellIndex( CoordinatesType( i,j,dimensions.z()-1) ) ] = timeFunctionDerivationValue*analyticSpaceFunction.getF(vertex);
-      }
-   }
-
-   #ifdef HAVE_OPENMP
-    #pragma omp parallel for firstprivate(coordinates) private(vertex)
-   #endif
-   for(IndexType i=0; i<dimensions.x(); i++)
-   {
-      for(IndexType j=0; j<dimensions.y(); j++)
-      {
-         coordinates.y()=i;
-         coordinates.z()=j;
-         
-         coordinates.x()=0;
-         vertex = mesh.getCellCenter( coordinates );
-         u[mesh.getCellIndex( CoordinatesType( 0,i,j) ) ] = timeFunctionDerivationValue*analyticSpaceFunction.getF(vertex);
-         
-         coordinates.x()=dimensions.x()-1;
-         vertex = mesh.getCellCenter( coordinates );
-         u[mesh.getCellIndex( CoordinatesType( dimensions.x()-1,i,j ) ) ] = timeFunctionDerivationValue*analyticSpaceFunction.getF(vertex);
-      }
-   }
-
-   #ifdef HAVE_OPENMP
-    #pragma omp parallel for firstprivate(coordinates) private(vertex)
-   #endif
-   for(IndexType i=0; i<dimensions.x(); i++)
-   {
-      for(IndexType j=0; j<dimensions.y(); j++)
-      {
-         coordinates.x()=i;
-         coordinates.z()=j;
-         
-         coordinates.y()=0;
-         vertex = mesh.getCellCenter( coordinates );         
-         u[mesh.getCellIndex( CoordinatesType( i,0,j) ) ] = timeFunctionDerivationValue*analyticSpaceFunction.getF(vertex);
-         
-         coordinates.y()=dimensions.y()-1;
-         vertex = mesh.getCellCenter( coordinates );   
-         u[mesh.getCellIndex( CoordinatesType( i,dimensions.y()-1,j) ) ] = timeFunctionDerivationValue*analyticSpaceFunction.getF(vertex);
-      }
-   } 
-}
-
 
 #endif	/* TNLDIRICHLETBOUNDARYCONDITIONS_IMPL_H */
 
diff --git a/examples/heat-equation/tnlLinearDiffusion.h b/examples/heat-equation/tnlLinearDiffusion.h
index 0263b3c2b6ec002fe3cf8ebd1be50a32462dbed9..2b5309c07e90956d1b7f0c2e3e07e5e3da9f2aff 100644
--- a/examples/heat-equation/tnlLinearDiffusion.h
+++ b/examples/heat-equation/tnlLinearDiffusion.h
@@ -26,27 +26,18 @@ class tnlLinearDiffusion< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index
    typedef Real RealType;
    typedef Device DeviceType;
    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,
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   void explicitUpdate( const RealType& time,
+                        const RealType& tau,
+                        const MeshType& mesh,
+                        const IndexType cellIndex,
                         const CoordinatesType& coordinates,
-                        Vector& _u,
-                        Vector& _fu );
-
-   template< typename Vector >
-   void getExplicitRHS( const MeshType& mesh,
-                        Vector& _u,
-                        Vector& _fu );
+                        Vector& u,
+                        Vector& fu );
    
 };
 
@@ -65,28 +56,18 @@ class tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index
    typedef Real RealType;
    typedef Device DeviceType;
    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,
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   void explicitUpdate( 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,
-                        Vector& _u,
-                        Vector& _fu);
+                        Vector& u,
+                        Vector& fu );
    
 };
 
@@ -105,28 +86,18 @@ class tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index
    typedef Real RealType;
    typedef Device DeviceType;
    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,
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   void explicitUpdate( 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,
-                        Vector& _u,
-                        Vector& _fu );
+                        Vector& u,
+                        Vector& fu );
    
 };
 
diff --git a/examples/heat-equation/tnlLinearDiffusion_impl.h b/examples/heat-equation/tnlLinearDiffusion_impl.h
index a9db826923f7f82bea455b22ce9dbfa0284dc60d..d588a992d8da47453ced99612a1bf6ab3050053b 100644
--- a/examples/heat-equation/tnlLinearDiffusion_impl.h
+++ b/examples/heat-equation/tnlLinearDiffusion_impl.h
@@ -13,66 +13,18 @@ template< typename MeshReal,
    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,
-          typename Real,
-          typename Index >
-   template< typename Vector >
-void tnlLinearDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
-getExplicitRHS( const MeshType& mesh,
+explicitUpdate( const RealType& time,
+                const RealType& tau,
+                const MeshType& mesh,
+                const IndexType cellIndex,
                 const CoordinatesType& coordinates,
-                Vector& _u,
-                Vector& _fu )
+                Vector& u,
+                Vector& fu )
 {
-   if((coordinates.x() <= 0) || (coordinates.x() >= (mesh.getDimensions().x()-1)))
-   {
-      cerr << "It is not possible to compute Laplace operator on the given dof. Coordinates are "
-           <<coordinates.x()<<"in x-axis. ";
-      return;
-   }
-   
-   _fu[mesh.getElementIndex(coordinates.x())]=((_u[mesh.getElementIndex(coordinates.x()-1)]-
-           2*_u[mesh.getElementIndex(coordinates.x())]+_u[mesh.getElementIndex(coordinates.x()+1)])/
-           (mesh.getParametricStep().x()*mesh.getParametricStep().x()));
+   fu[ cellIndex ] = ( u[ mesh.getCellXPredecessor( cellIndex ) ]
+                       - 2.0 * u[ cellIndex ]
+                       + u[ mesh.getCellXSuccessor( cellIndex ) ] ) * mesh.getHxSquareInverse();
 }
-   
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          typename Index >
-   template< typename Vector >
-void tnlLinearDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
-getExplicitRHS( const MeshType& mesh,
-                Vector& _u,
-                Vector& _fu )
-{
-   RealType stepXSquare = mesh.getCellProportions().x()*mesh.getCellProportions().x();
-   
-   CoordinatesType dimensions = mesh.getDimensions();
-
-   #ifdef HAVE_OPENMP
-    #pragma omp parallel for
-   #endif
-   for(IndexType i=1;i<(dimensions.x()-1);i++)
-   {
-      _fu[i]=(_u[i-1]-2.0*_u[i]+_u[i+1])/(stepXSquare);
-   }
-}
-
 
 template< typename MeshReal,
           typename Device,
@@ -82,75 +34,20 @@ template< typename MeshReal,
    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,
-          typename Real,
-          typename Index >
-   template< typename Vector >
-void tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
-getExplicitRHS( const MeshType& mesh,
+explicitUpdate( const RealType& time,
+                const RealType& tau,
+                const MeshType& mesh,
+                const IndexType cellIndex,
                 const CoordinatesType& coordinates,
-                Vector& _u,
-                Vector& _fu )
+                Vector& u,
+                Vector& fu )
 {
-   if((coordinates.x() <= 0) || (coordinates.x() >= (mesh.getDimensions().x()-1)) || 
-      (coordinates.y() <= 0 ) || (coordinates.y() >= (mesh.getDimensions().y()-1)))
-   {
-      cerr << "It is not possible to compute Laplace operator on the given dof. Coordinates are "
-           <<coordinates.x()<<"in x-axis, "<<coordinates.y()<<"in y-axis. ";
-      return;
-   }
-
-   _fu[mesh.getElementIndex(coordinates.x(),coordinates.y())]=((_u[mesh.getElementIndex(coordinates.x()-1,coordinates.y())]-
-           2*_u[mesh.getElementIndex(coordinates.x(),coordinates.y())]+_u[mesh.getElementIndex(coordinates.x()+1,coordinates.y())])/
-           (mesh.getParametricStep().x()*mesh.getParametricStep().x()));
-   _fu[mesh.getElementIndex(coordinates.x(),coordinates.y())]+=((_u[mesh.getElementIndex(coordinates.x(),coordinates.y()-1)]-
-           2*_u[mesh.getElementIndex(coordinates.x(),coordinates.y())]+_u[mesh.getElementIndex(coordinates.x(),coordinates.y()+1)])/
-           (mesh.getParametricStep().y()*mesh.getParametricStep().y()));
-}
-   
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          typename Index >
-   template< typename Vector >
-void tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index>::
-getExplicitRHS( const MeshType& mesh,
-                Vector& _u,
-                Vector& _fu )
-{
-   RealType stepXSquare = mesh.getCellProportions().x()*mesh.getCellProportions().x();
-   RealType stepYSquare = mesh.getCellProportions().y()*mesh.getCellProportions().y(); 
-   
-   CoordinatesType dimensions = mesh.getDimensions();
-
-   #ifdef HAVE_OPENMP
-    #pragma omp parallel for
-   #endif
-   for(IndexType i=1;i<(dimensions.x()-1);i++)
-   {
-      for(IndexType j=1;j<(dimensions.y()-1);j++)
-      {
-         _fu[j*mesh.getDimensions().x()+i]=(_u[j*mesh.getDimensions().x()+i-1]-2*_u[j*mesh.getDimensions().x()+i]+
-             _u[j*mesh.getDimensions().x()+i+1])/(stepXSquare)+(_u[(j-1)*mesh.getDimensions().x()+i]-2*_u[j*mesh.getDimensions().x()+i]+
-             _u[(j+1)*mesh.getDimensions().x()+i])/(stepYSquare);
-      }
-   }
+   fu[ cellIndex ] = ( u[ mesh.getCellXPredecessor( cellIndex ) ]
+                       - 2.0 * u[ cellIndex ]
+                       + u[ mesh.getCellXSuccessor( cellIndex ) ] ) * mesh.getHxSquareInverse() +
+                     ( u[ mesh.getCellYPredecessor( cellIndex ) ]
+                       - 2.0 * u[ cellIndex ]
+                       + u[ mesh.getCellYSuccessor( cellIndex ) ] ) * mesh.getHySquareInverse();
 }
 
 template< typename MeshReal,
@@ -161,94 +58,23 @@ template< typename MeshReal,
    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,
-          typename MeshIndex,
-          typename Real,
-          typename Index >
-   template< typename Vector >
-void tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
-getExplicitRHS( const MeshType& mesh,
+explicitUpdate( const RealType& time,
+                const RealType& tau,
+                const MeshType& mesh,
+                const IndexType cellIndex,
                 const CoordinatesType& coordinates,
-                Vector& _u,
-                Vector& _fu )
-{
-   if((coordinates.x() <= 0) || (coordinates.x() >= (mesh.getDimensions().x()-1)) || 
-      (coordinates.y() <= 0 ) || (coordinates.y() >= (mesh.getDimensions().y()-1)) || 
-      (coordinates.y() <= 0 ) || (coordinates.z() >= (mesh.getDimensions().z()-1)))
-   {
-      cerr << "It is not possible to compute Laplace operator on given dof. Coordinates are "
-           <<coordinates.x()<<"in x-axis, "<<coordinates.y()<<"in y-axis, "<<coordinates.z()<<"in z-axis. ";
-      return;
-   }
-   
-   _fu[mesh.getElementIndex(coordinates.x(),coordinates.y(),coordinates.z())]=
-           ((_u[mesh.getElementIndex(coordinates.x()-1,coordinates.y(),coordinates.z())]-
-           2*_u[mesh.getElementIndex(coordinates.x(),coordinates.y(),coordinates.z())]+
-           _u[mesh.getElementIndex(coordinates.x()+1,coordinates.y(),coordinates.z())])/
-           (mesh.getParametricStep().x()*mesh.getParametricStep().x()));
-   _fu[mesh.getElementIndex(coordinates.x(),coordinates.y(),coordinates.z())]+=
-           ((_u[mesh.getElementIndex(coordinates.x(),coordinates.y()-1,coordinates.z())]-
-           2*_u[mesh.getElementIndex(coordinates.x(),coordinates.y(),coordinates.z())]+
-           _u[mesh.getElementIndex(coordinates.x(),coordinates.y()+1,coordinates.z())])/
-           (mesh.getParametricStep().y()*mesh.getParametricStep().y()));
-   _fu[mesh.getElementIndex(coordinates.x(),coordinates.y(),coordinates.z())]+=
-           ((_u[mesh.getElementIndex(coordinates.x(),coordinates.y(),coordinates.z()-1)]-
-           2*_u[mesh.getElementIndex(coordinates.x(),coordinates.y(),coordinates.z())]+
-           _u[mesh.getElementIndex(coordinates.x(),coordinates.y(),coordinates.z()+1)])/
-           (mesh.getParametricStep().z()*mesh.getParametricStep().z()));
-}
-   
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          typename Index >
-   template< typename Vector >
-void tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
-getExplicitRHS( const MeshType& mesh,
-                Vector& _u,
-                Vector& _fu )
+                Vector& u,
+                Vector& fu )
 {
-   RealType stepXSquare = mesh.getCellProportions().x()*mesh.getCellProportions().x();
-   RealType stepYSquare = mesh.getCellProportions().y()*mesh.getCellProportions().y();   
-   RealType stepZSquare = mesh.getCellProportions().z()*mesh.getCellProportions().z();  
-   
-   CoordinatesType dimensions = mesh.getDimensions();
-
-   #ifdef HAVE_OPENMP
-    #pragma omp parallel for
-   #endif
-   for(IndexType i=1;i<(dimensions.x()-1);i++)
-   {
-      for(IndexType j=1;j<(dimensions.y()-1);j++)
-      {
-         for(IndexType k=1;k<(dimensions.z()-1);k++)
-         {
-            _fu[ mesh.getCellIndex( CoordinatesType( i,j,k ) ) ] = (_u[mesh.getCellIndex( CoordinatesType( i-1,j,k) ) ]-
-                    2*_u[mesh.getCellIndex( CoordinatesType( i,j,k) ) ]+_u[mesh.getCellIndex( CoordinatesType( i+1,j,k) )])/
-                    (stepXSquare) + (_u[mesh.getCellIndex( CoordinatesType( i,j-1,k) ) ]-
-                    2*_u[mesh.getCellIndex( CoordinatesType( i,j,k) ) ]+_u[mesh.getCellIndex( CoordinatesType( i,j+1,k) )])/
-                    (stepYSquare) + (_u[mesh.getCellIndex( CoordinatesType( i,j,k-1) )]-
-                    2*_u[mesh.getCellIndex( CoordinatesType( i,j,k) )]+_u[mesh.getCellIndex( CoordinatesType( i,j,k+1) )])/
-                    (stepZSquare);
-         }
-      }
-   }   
+   fu[ cellIndex ] = ( u[ mesh.getCellXPredecessor( cellIndex ) ]
+                       - 2.0 * u[ cellIndex ]
+                       + u[ mesh.getCellXSuccessor( cellIndex ) ] ) * mesh.getHxSquareInverse() +
+                     ( u[ mesh.getCellYPredecessor( cellIndex ) ]
+                       - 2.0 * u[ cellIndex ]
+                       + u[ mesh.getCellYSuccessor( cellIndex ) ] ) * mesh.getHySquareInverse() +
+                     ( u[ mesh.getCellZPredecessor( cellIndex ) ]
+                       - 2.0 * u[ cellIndex ]
+                       + u[ mesh.getCellZSuccessor( cellIndex ) ] ) * mesh.getHzSquareInverse();
 }
 
 
diff --git a/src/config/tnlConfigDescription.cpp b/src/config/tnlConfigDescription.cpp
index 3a1b5a9e7bea0c32ad17253343017a81e6269641..6989ad5dd5096d2528ae151f7ad1176f7a0a2454 100644
--- a/src/config/tnlConfigDescription.cpp
+++ b/src/config/tnlConfigDescription.cpp
@@ -167,6 +167,7 @@ bool tnlConfigDescription :: checkMissingEntries( tnlParameterContainer& paramet
       cerr << "Add the following missing  parameters to the command line: " << endl << "   ";
       for( int i = 0; i < missingParameters.getSize(); i++ )
          cerr << "--" << missingParameters[ i ] << " ... ";
+      cerr << endl;
       return false;
    }
    return true;
diff --git a/src/implementation/mesh/CMakeLists.txt b/src/implementation/mesh/CMakeLists.txt
index 9d5f0578790f2d720a5bdc2bbb7727238bad9fb4..5a46422a511a93f5ce8268809ee135d12418bc1c 100755
--- a/src/implementation/mesh/CMakeLists.txt
+++ b/src/implementation/mesh/CMakeLists.txt
@@ -1,7 +1,6 @@
 SET( headers tnlGrid1D_impl.h
              tnlGrid2D_impl.h
              tnlGrid3D_impl.h
-             tnlGridCellNeighbours_impl.h
              tnlTraversal_Grid1D_impl.h
              tnlTraversal_Grid2D_impl.h
              tnlTraversal_Grid3D_impl.h )
diff --git a/src/implementation/mesh/tnlGrid1D_impl.h b/src/implementation/mesh/tnlGrid1D_impl.h
index c9acc9456655f72d23f37bb9fe9a357368e8c829..4b2e6752f3215f558ef1f2bd1fab7dba82097ddb 100644
--- a/src/implementation/mesh/tnlGrid1D_impl.h
+++ b/src/implementation/mesh/tnlGrid1D_impl.h
@@ -55,6 +55,21 @@ tnlString tnlGrid< 1, Real, Device, Index > :: getTypeVirtual() const
    return this -> getType();
 }
 
+template< typename Real,
+          typename Device,
+          typename Index >
+void tnlGrid< 1, Real, Device, Index > :: computeSpaceSteps()
+{
+   if( this->getDimensions().x() != 0 )
+   {
+      this->cellProportions = this->proportions.x() / ( Real ) this->getDimensions().x();
+      this->hx = this->proportions.x() / ( Real )  this->getDimensions().x();
+      this->hxSquare = this->hx * this->hx;
+      this->hxInverse = 1.0 / this->hx;
+      this->hxSquareInverse = this->hxInverse * this->hxInverse;
+   }
+}
+
 template< typename Real,
           typename Device,
           typename Index  >
@@ -64,7 +79,7 @@ void tnlGrid< 1, Real, Device, Index >::setDimensions( const Index xSize )
    this->dimensions.x() = xSize;
    this->numberOfCells = xSize;
    this->numberOfVertices = xSize + 1;
-   this->cellProportions = this->proportions.x() / ( Real ) xSize;
+   computeSpaceSteps();
 }
 
 template< typename Real,
@@ -95,8 +110,7 @@ void tnlGrid< 1, Real, Device, Index > :: setDomain( const VertexType& origin,
 {
    this->origin = origin;
    this->proportions = proportions;
-   if( this->dimensions.x() != 0 )
-      this->cellProportions.x() = this->proportions.x() / ( Real ) this->dimensions.x();
+   computeSpaceSteps();
 }
 
 template< typename Real,
@@ -197,6 +211,80 @@ tnlGrid< 1, Real, Device, Index > :: getVertexCoordinates( const Index vertexInd
    return CoordinatesType( vertexIndex );
 }
 
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+Index tnlGrid< 1, Real, Device, Index > :: getCellXPredecessor( const IndexType& cellIndex ) const
+{
+   tnlAssert( cellIndex > 0 && cellIndex < this->getNumberOfCells(),
+              cerr << " cellIndex = " << cellIndex
+                   << " this->getNumberOfCells() = " << this->getNumberOfCells()
+                   << " this->getName() " << this->getName(); );
+   return cellIndex - 1;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+Index tnlGrid< 1, Real, Device, Index > :: getCellXSuccessor( const IndexType& cellIndex ) const
+{
+   tnlAssert( cellIndex >= 0 && cellIndex < this->getNumberOfCells() - 1,
+              cerr << " cellIndex = " << cellIndex
+                   << " this->getNumberOfCells() - 1 = " << this->getNumberOfCells() - 1
+                   << " this->getName() " << this->getName(); );
+   return cellIndex + 1;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 1, Real, Device, Index > :: getHx() const
+{
+   return this->hx;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 1, Real, Device, Index > :: getHxSquare() const
+{
+   return this->hxSquare;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 1, Real, Device, Index > :: getHxInverse() const
+{
+   return this->hxInverse;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 1, Real, Device, Index > :: getHxSquareInverse() const
+{
+   return this->hxSquareInverse;
+}
+
 template< typename Real,
           typename Device,
           typename Index >
diff --git a/src/implementation/mesh/tnlGrid2D_impl.h b/src/implementation/mesh/tnlGrid2D_impl.h
index afc827381f57ee9d02ce6f02ac75752850455eb7..b7bfdde4cff46a49f7241fcd0de800e3e359ad46 100644
--- a/src/implementation/mesh/tnlGrid2D_impl.h
+++ b/src/implementation/mesh/tnlGrid2D_impl.h
@@ -57,6 +57,28 @@ tnlString tnlGrid< 2, Real, Device, Index > :: getTypeVirtual() const
    return this -> getType();
 }
 
+template< typename Real,
+          typename Device,
+          typename Index >
+void tnlGrid< 2, Real, Device, Index > :: computeSpaceSteps()
+{
+   if( this->getDimensions().x() > 0 && this->getDimensions().y() > 0 )
+   {
+      this->cellProportions.x() = this->proportions.x() / ( Real ) this->getDimensions().x();
+      this->cellProportions.y() = this->proportions.y() / ( Real ) this->getDimensions().y();
+      this->hx = this->proportions.x() / ( Real ) this->getDimensions().x();
+      this->hxSquare = this->hx * this->hx;
+      this->hxInverse = 1.0 / this->hx;
+      this->hxSquareInverse = this->hxInverse * this->hxInverse;
+      this->hy = this->proportions.y() / ( Real ) this->getDimensions().y();
+      this->hySquare = this->hy * this->hy;
+      this->hyInverse = 1.0 / this->hy;
+      this->hySquareInverse = this->hyInverse * this->hyInverse;
+      this->hxhy = this->hx * this->hy;
+      this->hxhyInverse = 1.0 / this->hxhy;
+   }
+}
+
 template< typename Real,
           typename Device,
           typename Index >
@@ -72,8 +94,7 @@ void tnlGrid< 2, Real, Device, Index > :: setDimensions( const Index xSize, cons
    this->numberOfNyFaces = xSize * ( ySize + 1 );
    this->numberOfFaces = this->numberOfNxFaces + this->numberOfNyFaces;
    this->numberOfVertices = ( xSize + 1 ) * ( ySize + 1 );
-   this->cellProportions.x() = this->proportions.x() / ( Real ) xSize;
-   this->cellProportions.y() = this->proportions.y() / ( Real ) ySize;
+   computeSpaceSteps();
 }
 
 template< typename Real,
@@ -104,11 +125,7 @@ void tnlGrid< 2, Real, Device, Index > :: setDomain( const VertexType& origin,
 {
    this->origin = origin;
    this->proportions = proportions;
-   if( this->getDimensions().x() > 0 && this->getDimensions().y() > 0 )
-   {
-      this->cellProportions.x() = proportions.x() / ( Real ) this->getDimensions().x();
-      this->cellProportions.y() = proportions.y() / ( Real ) this->getDimensions().y();
-   }
+   computeSpaceSteps();
 }
 
 template< typename Real,
@@ -279,6 +296,176 @@ tnlGrid< 2, Real, Device, Index > :: getVertexCoordinates( const Index vertexInd
    return CoordinatesType( vertexIndex % aux, vertexIndex / aux );
 }
 
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+Index tnlGrid< 2, Real, Device, Index > :: getCellXPredecessor( const IndexType& cellIndex ) const
+{
+   tnlAssert( cellIndex > 0 && cellIndex < this->getNumberOfCells(),
+              cerr << " cellIndex = " << cellIndex
+                   << " this->getNumberOfCells() = " << this->getNumberOfCells()
+                   << " this->getName() " << this->getName(); );
+   return cellIndex - 1;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+Index tnlGrid< 2, Real, Device, Index > :: getCellXSuccessor( const IndexType& cellIndex ) const
+{
+   tnlAssert( cellIndex >= 0 && cellIndex < this->getNumberOfCells() - 1,
+              cerr << " cellIndex = " << cellIndex
+                   << " this->getNumberOfCells() - 1 = " << this->getNumberOfCells() - 1
+                   << " this->getName() " << this->getName(); );
+   return cellIndex + 1;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+Index tnlGrid< 2, Real, Device, Index > :: getCellYPredecessor( const IndexType& cellIndex ) const
+{
+   tnlAssert( cellIndex >= this->getDimensions().x() && cellIndex < this->getNumberOfCells(),
+              cerr << " cellIndex = " << cellIndex
+                   << " this->getNumberOfCells() = " << this->getNumberOfCells()
+                   << " this->getName() " << this->getName(); );
+   return cellIndex - this->getDimensions().x();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+Index tnlGrid< 2, Real, Device, Index > :: getCellYSuccessor( const IndexType& cellIndex ) const
+{
+   tnlAssert( cellIndex >= 0 && cellIndex < this->getNumberOfCells() - this->getDimensions().x(),
+              cerr << " cellIndex = " << cellIndex
+                   << " this->getNumberOfCells() - this->getDimensions().x() = " << this->getNumberOfCells() - this->getDimensions().x()
+                   << " this->getName() " << this->getName(); );
+   return cellIndex + this->getDimensions().x();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 2, Real, Device, Index > :: getHx() const
+{
+   return this->hx;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 2, Real, Device, Index > :: getHxSquare() const
+{
+   return this->hxSquare;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 2, Real, Device, Index > :: getHxInverse() const
+{
+   return this->hxInverse;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 2, Real, Device, Index > :: getHxSquareInverse() const
+{
+   return this->hxSquareInverse;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 2, Real, Device, Index > :: getHy() const
+{
+   return this->hy;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 2, Real, Device, Index > :: getHySquare() const
+{
+   return this->hySquare;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 2, Real, Device, Index > :: getHyInverse() const
+{
+   return this->hyInverse;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 2, Real, Device, Index > :: getHySquareInverse() const
+{
+   return this->hySquareInverse;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 2, Real, Device, Index > :: getHxHy() const
+{
+   return this->hxhy;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 2, Real, Device, Index > :: getHxHyInverse() const
+{
+   return this->hxhyInverse;
+}
+
 template< typename Real,
           typename Device,
           typename Index >
diff --git a/src/implementation/mesh/tnlGrid3D_impl.h b/src/implementation/mesh/tnlGrid3D_impl.h
index 3828e4b571d1b25325b4db923ea871290627eb5e..9be4b24f04239b7c7c85e6458100f1948309e276 100644
--- a/src/implementation/mesh/tnlGrid3D_impl.h
+++ b/src/implementation/mesh/tnlGrid3D_impl.h
@@ -60,6 +60,39 @@ tnlString tnlGrid< 3, Real, Device, Index > :: getTypeVirtual() const
    return this -> getType();
 }
 
+template< typename Real,
+          typename Device,
+          typename Index >
+void tnlGrid< 3, Real, Device, Index > :: computeSpaceSteps()
+{
+   if( this->getDimensions().x() > 0 &&
+       this->getDimensions().y() > 0 &&
+       this->getDimensions().z() > 0 )
+   {
+      this->cellProportions.x() = this->proportions.x() / ( Real ) this->getDimensions().x();
+      this->cellProportions.y() = this->proportions.y() / ( Real ) this->getDimensions().y();
+      this->cellProportions.z() = this->proportions.z() / ( Real ) this->getDimensions().z();
+      this->hx = this->proportions.x() / ( Real ) this->getDimensions().x();
+      this->hxSquare = this->hx * this->hx;
+      this->hxInverse = 1.0 / this->hx;
+      this->hxSquareInverse = this->hxInverse * this->hxInverse;
+      this->hy = this->proportions.y() / ( Real ) this->getDimensions().y();
+      this->hySquare = this->hy * this->hy;
+      this->hyInverse = 1.0 / this->hy;
+      this->hySquareInverse = this->hyInverse * this->hyInverse;
+      this->hz = this->proportions.z() / ( Real ) this->getDimensions().z();
+      this->hzSquare = this->hz * this->hz;
+      this->hzInverse = 1.0 / this->hz;
+      this->hzSquareInverse = this->hzInverse * this->hzInverse;
+      this->hxhy = this->hx * this->hy;
+      this->hxhz = this->hx * this->hz;
+      this->hyhz = this->hy * this->hz;
+      this->hxhyInverse = 1.0 / this->hxhy;
+      this->hxhzInverse = 1.0 / this->hxhz;
+      this->hyhzInverse = 1.0 / this->hyhz;
+   }
+}
+
 template< typename Real,
           typename Device,
           typename Index >
@@ -89,9 +122,9 @@ void tnlGrid< 3, Real, Device, Index > :: setDimensions( const Index xSize, cons
                          this->numberOfDzEdges;
    this->numberOfVertices = ( xSize + 1 ) * ( ySize + 1 ) * ( zSize + 1 );
    
-   this->cellProportions.x() = this->proportions.x() / ( Real ) xSize;
-   this->cellProportions.y() = this->proportions.y() / ( Real ) ySize;
-   this->cellProportions.z() = this->proportions.z() / ( Real ) zSize;
+   this->cellZNeighboursStep = xSize * ySize;
+
+   computeSpaceSteps();
 }
 
 template< typename Real,
@@ -122,14 +155,7 @@ void tnlGrid< 3, Real, Device, Index > :: setDomain( const VertexType& origin,
 {
    this->origin = origin;
    this->proportions = proportions;
-   if( this->getDimensions().x() > 0 &&
-       this->getDimensions().y() > 0 &&
-       this->getDimensions().z() > 0 )
-   {
-      this->cellProportions.x() = proportions.x() / ( Real ) this->getDimensions().x();
-      this->cellProportions.y() = proportions.y() / ( Real ) this->getDimensions().y();
-      this->cellProportions.z() = proportions.z() / ( Real ) this->getDimensions().z();
-   }
+   computeSpaceSteps();
 }
 
 template< typename Real,
@@ -463,6 +489,294 @@ tnlGrid< 3, Real, Device, Index > :: getVertexCoordinates( const Index vertexInd
                            vertexIndex / ( auxX * auxY ) );
 }
 
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+Index tnlGrid< 3, Real, Device, Index > :: getCellXPredecessor( const IndexType& cellIndex ) const
+{
+   tnlAssert( cellIndex > 0 && cellIndex < this->getNumberOfCells(),
+              cerr << " cellIndex = " << cellIndex
+                   << " this->getNumberOfCells() = " << this->getNumberOfCells()
+                   << " this->getName() " << this->getName(); );
+   return cellIndex - 1;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+Index tnlGrid< 3, Real, Device, Index > :: getCellXSuccessor( const IndexType& cellIndex ) const
+{
+   tnlAssert( cellIndex >= 0 && cellIndex < this->getNumberOfCells() - 1,
+              cerr << " cellIndex = " << cellIndex
+                   << " this->getNumberOfCells() - 1 = " << this->getNumberOfCells() - 1
+                   << " this->getName() " << this->getName(); );
+   return cellIndex + 1;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+Index tnlGrid< 3, Real, Device, Index > :: getCellYPredecessor( const IndexType& cellIndex ) const
+{
+   tnlAssert( cellIndex >= this->getDimensions().x() && cellIndex < this->getNumberOfCells(),
+              cerr << " cellIndex = " << cellIndex
+                   << " this->getNumberOfCells() = " << this->getNumberOfCells()
+                   << " this->getName() " << this->getName(); );
+   return cellIndex - this->getDimensions().x();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+Index tnlGrid< 3, Real, Device, Index > :: getCellYSuccessor( const IndexType& cellIndex ) const
+{
+   tnlAssert( cellIndex >= 0 && cellIndex < this->getNumberOfCells() - this->getDimensions().x(),
+              cerr << " cellIndex = " << cellIndex
+                   << " this->getNumberOfCells() - this->getDimensions().x() = " << this->getNumberOfCells() - this->getDimensions().x()
+                   << " this->getName() " << this->getName(); );
+   return cellIndex + this->getDimensions().x();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+Index tnlGrid< 3, Real, Device, Index > :: getCellZPredecessor( const IndexType& cellIndex ) const
+{
+   tnlAssert( cellIndex >= this->cellZNeighboursStep && cellIndex < this->getNumberOfCells(),
+              cerr << " cellIndex = " << cellIndex
+                   << " this->getNumberOfCells() = " << this->getNumberOfCells()
+                   << " this->getName() " << this->getName(); );
+   return cellIndex - this->cellZNeighboursStep;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+Index tnlGrid< 3, Real, Device, Index > :: getCellZSuccessor( const IndexType& cellIndex ) const
+{
+   tnlAssert( cellIndex >= 0 && cellIndex < this->getNumberOfCells() - this->cellZNeighboursStep,
+              cerr << " cellIndex = " << cellIndex
+                   << " this->getNumberOfCells() - this->cellZNeighboursStep = " << this->getNumberOfCells() - this->cellZNeighboursStep
+                   << " this->getName() " << this->getName(); );
+   return cellIndex + this->cellZNeighboursStep;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 3, Real, Device, Index > :: getHx() const
+{
+   return this->hx;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 3, Real, Device, Index > :: getHxSquare() const
+{
+   return this->hxSquare;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 3, Real, Device, Index > :: getHxInverse() const
+{
+   return this->hxInverse;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 3, Real, Device, Index > :: getHxSquareInverse() const
+{
+   return this->hxSquareInverse;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 3, Real, Device, Index > :: getHy() const
+{
+   return this->hy;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 3, Real, Device, Index > :: getHySquare() const
+{
+   return this->hySquare;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 3, Real, Device, Index > :: getHyInverse() const
+{
+   return this->hyInverse;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 3, Real, Device, Index > :: getHySquareInverse() const
+{
+   return this->hySquareInverse;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 3, Real, Device, Index > :: getHz() const
+{
+   return this->hz;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 3, Real, Device, Index > :: getHzSquare() const
+{
+   return this->hzSquare;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 3, Real, Device, Index > :: getHzInverse() const
+{
+   return this->hzInverse;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 3, Real, Device, Index > :: getHzSquareInverse() const
+{
+   return this->hzSquareInverse;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 3, Real, Device, Index > :: getHxHy() const
+{
+   return this->hxhy;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 3, Real, Device, Index > :: getHxHz() const
+{
+   return this->hxhz;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 3, Real, Device, Index > :: getHyHz() const
+{
+   return this->hyhz;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 3, Real, Device, Index > :: getHxHyInverse() const
+{
+   return this->hxhyInverse;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 3, Real, Device, Index > :: getHxHzInverse() const
+{
+   return this->hxhzInverse;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+const Real& tnlGrid< 3, Real, Device, Index > :: getHyHzInverse() const
+{
+   return this->hyhzInverse;
+}
+
 template< typename Real,
           typename Device,
           typename Index >
diff --git a/src/implementation/mesh/tnlGridCellNeighbours_impl.h b/src/implementation/mesh/tnlGridCellNeighbours_impl.h
deleted file mode 100644
index b91f082976799a8ffa9e9fb002f281df1d26e84e..0000000000000000000000000000000000000000
--- a/src/implementation/mesh/tnlGridCellNeighbours_impl.h
+++ /dev/null
@@ -1,110 +0,0 @@
-/***************************************************************************
-                          tnlGridCellNeighbours_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 TNLGRIDCELLNEIGHBOURS_IMPL_H_
-#define TNLGRIDCELLNEIGHBOURS_IMPL_H_
-
-template< typename Real,
-          typename Device,
-          typename Index >
-tnlGridCellNeighbours< tnlGrid< 1, Real, Device, Index > >::
-tnlGridCellNeighbours( const GridType& grid,
-                       const CoordinatesType& cellCoordinates )
-{
-   const IndexType cellIndex = grid.getCellIndex( cellCoordinates );
-   this->xPredecessor = cellIndex - 1;
-   this->xSuccessor = cellIndex + 1;
-}
-
-template< typename Real,
-          typename Device,
-          typename Index >
-Index
-tnlGridCellNeighbours< tnlGrid< 1, Real, Device, Index > >::
-getXPredecessor() const
-{
-   return this->xPredecessor;
-}
-
-template< typename Real,
-          typename Device,
-          typename Index >
-Index
-tnlGridCellNeighbours< tnlGrid< 1, Real, Device, Index > >::
-getXSuccessor() const
-{
-   return this->xSuccessor;
-}
-
-template< typename Real,
-          typename Device,
-          typename Index >
-tnlGridCellNeighbours< tnlGrid< 2, Real, Device, Index > >::
-tnlGridCellNeighbours( const GridType& grid,
-                       const CoordinatesType& cellCoordinates )
-{
-   const IndexType cellIndex = grid.getCellIndex( cellCoordinates );
-   this->xPredecessor = cellIndex - 1;
-   this->xSuccessor = cellIndex + 1;
-   this->yPredecessor = cellIndex - grid.getDimensions().x();
-   this->ySuccessor = cellIndex + grid.getDimensions().x();
-}
-
-template< typename Real,
-          typename Device,
-          typename Index >
-Index
-tnlGridCellNeighbours< tnlGrid< 2, Real, Device, Index > >::
-getXPredecessor() const
-{
-   return this->xPredecessor;
-}
-
-template< typename Real,
-          typename Device,
-          typename Index >
-Index
-tnlGridCellNeighbours< tnlGrid< 2, Real, Device, Index > >::
-getXSuccessor() const
-{
-   return this->xSuccessor;
-}
-
-template< typename Real,
-          typename Device,
-          typename Index >
-Index
-tnlGridCellNeighbours< tnlGrid< 2, Real, Device, Index > >::
-getYPredecessor() const
-{
-   return this->yPredecessor;
-}
-
-template< typename Real,
-          typename Device,
-          typename Index >
-Index
-tnlGridCellNeighbours< tnlGrid< 2, Real, Device, Index > >::
-getYSuccessor() const
-{
-   return this->ySuccessor;
-}
-
-
-
-
-#endif /* TNLGRIDCELLNEIGHBOURS_IMPL_H_ */
diff --git a/src/implementation/mesh/tnlTraversal_Grid2D_impl.h b/src/implementation/mesh/tnlTraversal_Grid2D_impl.h
index 2130a6beb93adabade4fa274bc95aadffb25e573..d5c6cf749c74c151a145e22eac20a0e71162e544 100644
--- a/src/implementation/mesh/tnlTraversal_Grid2D_impl.h
+++ b/src/implementation/mesh/tnlTraversal_Grid2D_impl.h
@@ -29,12 +29,23 @@ tnlTraversal< tnlGrid< 2, Real, tnlHost, Index >, 2 >::
 processEntities( const GridType& grid,
                  UserData& userData,
                  BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
-                 InteriorEntitiesProcessor& interiorEntitesProcessor ) const
+                 InteriorEntitiesProcessor& interiorEntitiesProcessor ) const
 {
    /****
     * Traversing cells
     */
-
+   CoordinatesType coordinates;
+   const IndexType& xSize = grid.getDimensions().x();
+   const IndexType& ySize = grid.getDimensions().y();
+   for( coordinates.y() = 0; coordinates.y() < ySize; coordinates.y() ++ )
+      for( coordinates.x() = 0; coordinates.x() < xSize; coordinates.x() ++ )
+      {
+         const IndexType index = grid.getCellIndex( coordinates );
+         if( grid.isBoundaryCell( coordinates ) )
+            boundaryEntitiesProcessor.template processEntity< 2 >( grid, userData, index, coordinates );
+         else
+            interiorEntitiesProcessor.template processEntity< 2 >( grid, userData, index, coordinates );
+      }
 }
 
 template< typename Real,
diff --git a/src/implementation/mesh/tnlTraversal_Grid3D_impl.h b/src/implementation/mesh/tnlTraversal_Grid3D_impl.h
index a97d04a40392fd1e8f983231569148187cb066f1..8cc2a0db8bcc7d6d2d75e2a0f8221cde8363ec94 100644
--- a/src/implementation/mesh/tnlTraversal_Grid3D_impl.h
+++ b/src/implementation/mesh/tnlTraversal_Grid3D_impl.h
@@ -28,11 +28,25 @@ tnlTraversal< tnlGrid< 3, Real, tnlHost, Index >, 3 >::
 processEntities( const GridType& grid,
                  UserData& userData,
                  BoundaryEntitiesProcessor& boundaryEntitiesProcessor,
-                 InteriorEntitiesProcessor& interiorEntitesProcessor ) const
+                 InteriorEntitiesProcessor& interiorEntitiesProcessor ) const
 {
    /****
     * Traversing cells
     */
+   CoordinatesType coordinates;
+   const IndexType& xSize = grid.getDimensions().x();
+   const IndexType& ySize = grid.getDimensions().y();
+   const IndexType& zSize = grid.getDimensions().z();
+   for( coordinates.z() = 0; coordinates.z() < zSize; coordinates.z() ++ )
+      for( coordinates.y() = 0; coordinates.y() < ySize; coordinates.y() ++ )
+         for( coordinates.x() = 0; coordinates.x() < xSize; coordinates.x() ++ )
+         {
+            const IndexType index = grid.getCellIndex( coordinates );
+            if( grid.isBoundaryCell( coordinates ) )
+               boundaryEntitiesProcessor.template processEntity< 3 >( grid, userData, index, coordinates );
+            else
+               interiorEntitiesProcessor.template processEntity< 3 >( grid, userData, index, coordinates );
+         }
 }
 
 template< typename Real,
diff --git a/src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h b/src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h
index c86af510de6b55c294245cdbd9b765e2e8b92b6c..3bf8028bb3813c7ef0e601765cbb58a32bebef20 100644
--- a/src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h
+++ b/src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h
@@ -18,6 +18,7 @@
 #ifndef TNLLINEARDIFFUSION_IMPL_H_
 #define TNLLINEARDIFFUSION_IMPL_H_
 
+#ifdef UNDEF
 template< typename Real,
           typename Device,
           typename Index,
@@ -242,6 +243,8 @@ Real tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index, GridGeometry > > :: ge
             );
 }
 
+#endif /* UNDEF */
+
 /****
  * Specialization for the grids with no deformations (Identical grid geometry)
  */
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index afc30628d3c6bc1941358274981d6f4bff42ba52..1a55833a80446a379cf3f66d5e98b01e5821ac5d 100755
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -4,7 +4,6 @@ ADD_SUBDIRECTORY( traits )
 ADD_SUBDIRECTORY( topologies )
 
 SET( headers tnlGrid.h
-             tnlGridCellNeighbours.h
              tnlDummyMesh.h
              tnlGnuplotWriter.h
              tnlMesh.h
diff --git a/src/mesh/tnlGrid.h b/src/mesh/tnlGrid.h
index fc5e934b0b5ef897264a738f3d1aa368a680700c..b6de81e6473d33fa8945ef3541afd9e606f3f467 100644
--- a/src/mesh/tnlGrid.h
+++ b/src/mesh/tnlGrid.h
@@ -99,6 +99,36 @@ class tnlGrid< 1, Real, Device, Index > : public tnlObject
 #endif
    CoordinatesType getVertexCoordinates( const Index vertexCoordinates ) const;
 
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   IndexType getCellXPredecessor( const IndexType& cellIndex ) const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   IndexType getCellXSuccessor( const IndexType& cellIndex ) const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHx() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHxSquare() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHxInverse() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHxSquareInverse() const;
+
    /****
     * The type Vertex can have different Real type.
     */
@@ -167,12 +197,16 @@ class tnlGrid< 1, Real, Device, Index > : public tnlObject
 
    protected:
 
+   void computeSpaceSteps();
+
    CoordinatesType dimensions;
 
    VertexType origin, proportions, cellProportions;
 
    IndexType numberOfCells, numberOfVertices;
 
+   RealType hx, hxSquare, hxInverse, hxSquareInverse;
+
 };
 
 template< typename Real,
@@ -250,6 +284,76 @@ class tnlGrid< 2, Real, Device, Index > : public tnlObject
 #endif
    CoordinatesType getVertexCoordinates( const Index vertexIndex ) const;
 
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   IndexType getCellXPredecessor( const IndexType& cellIndex ) const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   IndexType getCellXSuccessor( const IndexType& cellIndex ) const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   IndexType getCellYPredecessor( const IndexType& cellIndex ) const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   IndexType getCellYSuccessor( const IndexType& cellIndex ) const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHx() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHxSquare() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHxInverse() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHxSquareInverse() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHy() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHySquare() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHyInverse() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHySquareInverse() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHxHy() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHxHyInverse() const;
+
    /****
     * The type Vertex can have different Real type.
     */
@@ -338,12 +442,18 @@ template< int nx, int ny, typename Vertex = VertexType >
 
    protected:
 
+   void computeSpaceSteps();
+
    CoordinatesType dimensions;
 
    VertexType origin, proportions, cellProportions;
 
    IndexType numberOfCells, numberOfNxFaces, numberOfNyFaces, numberOfFaces, numberOfVertices;
 
+   RealType hx, hxSquare, hxInverse, hxSquareInverse,
+            hy, hySquare, hyInverse, hySquareInverse,
+            hxhy, hxhyInverse;
+
 
 };
 
@@ -435,6 +545,126 @@ class tnlGrid< 3, Real, Device, Index > : public tnlObject
 #endif
    CoordinatesType getVertexCoordinates( const Index vertexIndex ) const;
 
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   IndexType getCellXPredecessor( const IndexType& cellIndex ) const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   IndexType getCellXSuccessor( const IndexType& cellIndex ) const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   IndexType getCellYPredecessor( const IndexType& cellIndex ) const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   IndexType getCellYSuccessor( const IndexType& cellIndex ) const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   IndexType getCellZPredecessor( const IndexType& cellIndex ) const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   IndexType getCellZSuccessor( const IndexType& cellIndex ) const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHx() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHxSquare() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHxInverse() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHxSquareInverse() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHy() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHySquare() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHyInverse() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHySquareInverse() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHz() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHzSquare() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHzInverse() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHzSquareInverse() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHxHy() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHxHz() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHyHz() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHxHyInverse() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHxHzInverse() const;
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   const RealType& getHyHzInverse() const;
+
    /****
     * The type Vertex can have different Real type.
     */
@@ -540,6 +770,8 @@ template< int dx, int dy, int dz, typename Vertex = VertexType >
 
    protected:
 
+   void computeSpaceSteps();
+
    CoordinatesType dimensions;
 
    VertexType origin, proportions, cellProportions;
@@ -548,6 +780,15 @@ template< int dx, int dy, int dz, typename Vertex = VertexType >
              numberOfNxFaces, numberOfNyFaces, numberOfNzFaces, numberOfNxAndNyFaces, numberOfFaces,
              numberOfDxEdges, numberOfDyEdges, numberOfDzEdges, numberOfDxAndDyEdges, numberOfEdges,
              numberOfVertices;
+   IndexType cellZNeighboursStep;
+
+   RealType hx, hxSquare, hxInverse, hxSquareInverse,
+            hy, hySquare, hyInverse, hySquareInverse,
+            hz, hzSquare, hzInverse, hzSquareInverse,
+            hxhy, hxhz, hyhz,
+            hxhyInverse, hxhzInverse, hyhzInverse;
+
+
 
 };
 
diff --git a/src/mesh/tnlGridCellNeighbours.h b/src/mesh/tnlGridCellNeighbours.h
deleted file mode 100644
index fd180486ec62c90e84f7210236dd372d855da1b7..0000000000000000000000000000000000000000
--- a/src/mesh/tnlGridCellNeighbours.h
+++ /dev/null
@@ -1,116 +0,0 @@
-/***************************************************************************
-                          tnlGridCellNeighbours.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 TNLGRIDCELLNEIGHBOURS_H_
-#define TNLGRIDCELLNEIGHBOURS_H_
-
-#include <mesh/tnlGrid.h>
-
-template< typename Grid >
-class tnlGridCellNeighbours{};
-
-template< typename Real,
-          typename Device,
-          typename Index >
-class tnlGridCellNeighbours< tnlGrid< 1, Real, Device, Index > >
-{
-   public:
-      typedef tnlGrid< 1, Real, Device, Index > GridType;
-      typedef Real RealType;
-      typedef Device DeviceType;
-      typedef Index IndexType;
-      typedef typename GridType::CoordinatesType CoordinatesType;
-      enum { Dimensions = GridType::Dimensions };
-
-      tnlGridCellNeighbours( const GridType& grid,
-                             const CoordinatesType& cellCoordinates );
-
-      const IndexType& getXPredecessor() const;
-
-      const IndexType& getXSuccessor() const;
-
-   protected:
-      IndexType xPredecessor, xSuccessor;
-};
-
-template< typename Real,
-          typename Device,
-          typename Index >
-class tnlGridCellNeighbours< tnlGrid< 2, Real, Device, Index > >
-{
-   public:
-      typedef tnlGrid< 2, Real, Device, Index > GridType;
-      typedef Real RealType;
-      typedef Device DeviceType;
-      typedef Index IndexType;
-      typedef typename GridType::CoordinatesType CoordinatesType;
-      enum { Dimensions = GridType::Dimensions };
-
-      tnlGridCellNeighbours( const GridType& grid,
-                             const CoordinatesType& cellCoordinates );
-
-      const IndexType& getXPredecessor() const;
-
-      const IndexType& getXSuccessor() const;
-
-      const IndexType& getYPredecessor() const;
-
-      const IndexType& getYSuccessor() const;
-
-
-   protected:
-      IndexType xPredecessor, xSuccessor,
-                yPredecessor, ySuccessor;
-};
-
-template< typename Real,
-          typename Device,
-          typename Index >
-class tnlGridCellNeighbours< tnlGrid< 2, Real, Device, Index > >
-{
-   public:
-      typedef tnlGrid< 2, Real, Device, Index > GridType;
-      typedef Real RealType;
-      typedef Device DeviceType;
-      typedef Index IndexType;
-      typedef typename GridType::CoordinatesType CoordinatesType;
-      enum { Dimensions = GridType::Dimensions };
-
-      tnlGridCellNeighbours( const GridType& grid,
-                             const CoordinatesType& cellCoordinates );
-
-      const IndexType& getXPredecessor() const;
-
-      const IndexType& getXSuccessor() const;
-
-      const IndexType& getYPredecessor() const;
-
-      const IndexType& getYSuccessor() const;
-
-      const IndexType& getZPredecessor() const;
-
-      const IndexType& getZSuccessor() const;
-
-   protected:
-      IndexType xPredecessor, xSuccessor,
-                yPredecessor, ySuccessor,
-                zPredecessor, zSuccessor;
-};
-
-#include <implementation/mesh/tnlGridCellNeighbours_impl.h>
-
-#endif /* TNLGRIDCELLNEIGHBOURS_H_ */
diff --git a/src/solvers/pde/tnlExplicitUpdater.h b/src/solvers/pde/tnlExplicitUpdater.h
index 0a04ef81252cd608b689049645927655c5606fb1..fde11e326a1965184cc14bc1919647d4a0a85559 100644
--- a/src/solvers/pde/tnlExplicitUpdater.h
+++ b/src/solvers/pde/tnlExplicitUpdater.h
@@ -177,13 +177,13 @@ class tnlExplicitUpdater< tnlGrid< Dimensions, Real, Device, Index >,
                                 const IndexType index,
                                 const CoordinatesType& coordinates )
             {
-               userData.interiorUpdater.update( userData.time,
-                                                userData.tau,
-                                                mesh,
-                                                index,
-                                                coordinates,
-                                                userData.u,
-                                                userData.fu );
+               userData.interiorUpdater.explicitUpdate( userData.time,
+                                                        userData.tau,
+                                                        mesh,
+                                                        index,
+                                                        coordinates,
+                                                        userData.u,
+                                                        userData.fu );
             }
 
       };
diff --git a/tools/src/CMakeLists.txt b/tools/src/CMakeLists.txt
index 4534b3d022d93fc3a8ae99d4d4ded77247528eda..1001e38f54f52ffaceb7e77a8a0dfc8cf48e9b1d 100755
--- a/tools/src/CMakeLists.txt
+++ b/tools/src/CMakeLists.txt
@@ -13,8 +13,8 @@ set( tnl_mesh_convert_SOURCES tnl-mesh-convert.cpp
                             tnl-mesh-convert.h )
                             
                             
-set( tnl_discrete_SOURCES tnl-discrete.cpp
-                          tnl-discrete.h )
+set( tnl_init_SOURCES tnl-init.cpp
+                      tnl-init.h )
 
 set( tnl_view_SOURCES tnl-view.cpp
                       tnl-view.h )
@@ -28,8 +28,8 @@ target_link_libraries (tnl-grid-setup${debugExt} tnl${debugExt}-${tnlVersion} )
 ADD_EXECUTABLE(tnl-mesh-convert${debugExt} ${tnl_mesh_convert_SOURCES})
 target_link_libraries (tnl-mesh-convert${debugExt} tnl${debugExt}-${tnlVersion} )
 
-ADD_EXECUTABLE(tnl-discrete${debugExt} ${tnl_discrete_SOURCES})
-target_link_libraries (tnl-discrete${debugExt} tnl${debugExt}-${tnlVersion} )
+ADD_EXECUTABLE(tnl-init${debugExt} ${tnl_init_SOURCES})
+target_link_libraries (tnl-init${debugExt} tnl${debugExt}-${tnlVersion} )
 
 ADD_EXECUTABLE(tnl-view${debugExt} ${tnl_view_SOURCES})
 target_link_libraries (tnl-view${debugExt} tnl${debugExt}-${tnlVersion} )
@@ -46,7 +46,7 @@ target_link_libraries (tnl-curve2gnuplot${debugExt} tnl${debugExt}-${tnlVersion}
 #ADD_EXECUTABLE( tnl-matrix-convert${debugExt} ${tnlmatrixconvertsources} )
 #target_link_libraries( tnl-matrix-convert${debugExt} tnl${debugExt}-${tnlVersion} )
 
-INSTALL( TARGETS tnl-discrete${debugExt}
+INSTALL( TARGETS tnl-init${debugExt}
                  tnl-view${debugExt}
                  tnl-diff${debugExt}
                  tnl-grid-setup${debugExt}
diff --git a/tools/src/tnl-discrete.cpp b/tools/src/tnl-init.cpp
similarity index 96%
rename from tools/src/tnl-discrete.cpp
rename to tools/src/tnl-init.cpp
index 035e7287ebe32645752ff6255272b0272945804f..38db422c07dc789d96bf69391b6d7b85ccfce339 100644
--- a/tools/src/tnl-discrete.cpp
+++ b/tools/src/tnl-init.cpp
@@ -1,5 +1,5 @@
 /***************************************************************************
-                          tnl-discrete.cpp  -  description
+                          tnl-init.cpp  -  description
                              -------------------
     begin                : Nov 23, 2013
     copyright            : (C) 2013 by Tomas Oberhuber
@@ -15,7 +15,7 @@
  *                                                                         *
  ***************************************************************************/
 
-#include "tnl-discrete.h"
+#include "tnl-init.h"
 
 #include <core/tnlFile.h>
 #include <debug/tnlDebug.h>
diff --git a/tools/src/tnl-discrete.h b/tools/src/tnl-init.h
similarity index 96%
rename from tools/src/tnl-discrete.h
rename to tools/src/tnl-init.h
index 70a5de8aa47d0df050bec1048554e7469cb13c37..e4273ea632e539758a1b33fd57e61084517aa402 100644
--- a/tools/src/tnl-discrete.h
+++ b/tools/src/tnl-init.h
@@ -1,5 +1,5 @@
 /***************************************************************************
-                          tnl-discrete.h  -  description
+                          tnl-init.h  -  description
                              -------------------
     begin                : Nov 23, 2013
     copyright            : (C) 2013 by Tomas Oberhuber
@@ -15,8 +15,8 @@
  *                                                                         *
  ***************************************************************************/
 
-#ifndef TNL_DISCRETE_H_
-#define TNL_DISCRETE_H_
+#ifndef TNL_INIT_H_
+#define TNL_INIT_H_
 
 #include <config/tnlParameterContainer.h>
 #include <core/vectors/tnlVector.h>
@@ -190,12 +190,6 @@ bool resolveMesh( const tnlList< tnlString >& parsedMeshType,
    cout << "+ -> Setting mesh type to " << parsedMeshType[ 0 ] << " ... " << endl;
    if( parsedMeshType[ 0 ] == "tnlGrid" )
    {
-      tnlList< tnlString > parsedGeometryType;
-      if( ! parseObjectType( parsedMeshType[ 5 ], parsedGeometryType ) )
-      {
-         cerr << "Unable to parse the geometry type " << parsedMeshType[ 5 ] << "." << endl;
-         return false;
-      }
       typedef tnlGrid< Dimensions, RealType, tnlHost, IndexType > MeshType;
       return resolveFunction< MeshType >( parameters );
    }
@@ -247,4 +241,4 @@ bool resolveMeshType( const tnlList< tnlString >& parsedMeshType,
       return resolveRealType< 3 >( parsedMeshType, parameters );
 
 }
-#endif /* TNL_DISCRETE_H_ */
+#endif /* TNL_INIT_H_ */
diff --git a/tools/src/tnl-view.cpp b/tools/src/tnl-view.cpp
index 88b6807156412d6ab1b29d397cf753de8ada4d04..d5875a2de91dd23a9f4ee1848ca122c6170d969d 100644
--- a/tools/src/tnl-view.cpp
+++ b/tools/src/tnl-view.cpp
@@ -36,7 +36,7 @@ void setupConfig( tnlConfigDescription& config )
    config.addEntry        < tnlString >           ( "mesh", "Mesh file.", "mesh.tnl" );
    config.addRequiredList < tnlString >           ( "input-files", "Input files." );
    config.addList         < tnlString >           ( "output-files", "Output files." );
-   config.addEntry        < bool >                ( "check-output-file", "If the output file already exists, do not recreate it.", "false" );
+   config.addEntry        < bool >                ( "check-output-file", "If the output file already exists, do not recreate it.", false );
 
    config.addDelimiter( "Grid settings:");
    config.addList         < double >              ( "level-lines", "List of level sets which will be drawn." );
@@ -93,13 +93,7 @@ int main( int argc, char* argv[] )
    }
    if( parsedMeshType[ 0 ] == "tnlGrid" )
    {
-      tnlList< tnlString > parsedGeometryType;
-      if( ! parseObjectType( parsedMeshType[ 5 ], parsedGeometryType ) )
-      {
-         cerr << "Unable to parse the geometry type " << parsedMeshType[ 5 ] << "." << endl;
-         return false;
-      }
-      int dimensions = atoi( parsedGeometryType[ 1 ].getString() );
+      int dimensions = atoi( parsedMeshType[ 1 ].getString() );
       if( dimensions == 1 )
       {
          typedef tnlGrid< 1, double, tnlHost, int > MeshType;