diff --git a/examples/navier-stokes/navier-stokes.cfg.desc b/examples/navier-stokes/navier-stokes.cfg.desc
index 3b0d8920c55dc51ccbc3eef478f5e6eaa69db8a2..dc4775afc9ed1697bef233581f0bdd0696d69387 100644
--- a/examples/navier-stokes/navier-stokes.cfg.desc
+++ b/examples/navier-stokes/navier-stokes.cfg.desc
@@ -13,7 +13,7 @@ group Problem
    integer y-size(!)                    [Number of grid cells along y-axis.];   
    real max-inflow-velocity(1.0)        [];
    real max-outflow-velocity(0.8)       [];
-   real start-up(1.0)                   [Start-up time of the riser.];
+   real start-up(0.0)                   [Start-up time of the riser.];
    real mu(!)                           [Viscosity.];
    real p0(!)                           [Pressure for the initial condition.];
    real T(!)                            [Temperature.];
diff --git a/examples/navier-stokes/navierStokesSetter_impl.h b/examples/navier-stokes/navierStokesSetter_impl.h
index 8b1ebb132d7256604dcbc2d1094e1045446eee9c..b55d099496034c65941a02b4b5e45d5f5417bb49 100644
--- a/examples/navier-stokes/navierStokesSetter_impl.h
+++ b/examples/navier-stokes/navierStokesSetter_impl.h
@@ -20,6 +20,7 @@
 
 #include <mesh/tnlGrid.h>
 #include <schemes/euler/fvm/tnlLaxFridrichs.h>
+#include <schemes/gradient/tnlCentralFDMGradient.h>
 
 template< typename SolverStarter >
    template< typename RealType,
@@ -39,7 +40,13 @@ bool navierStokesSetter< SolverStarter > :: run( const tnlParameterContainer& pa
    {
       typedef tnlGrid< 2, RealType, DeviceType, IndexType > MeshType;
       if( schemeName == "lax-fridrichs" )
-         return solverStarter. run< navierStokesSolver< MeshType, tnlLaxFridrichs< MeshType > > >( parameters );
+         return solverStarter. run< navierStokesSolver< MeshType,
+                                                        tnlLaxFridrichs< MeshType,
+                                                                        tnlCentralFDMGradient< MeshType :: Dimensions,
+                                                                                               typename MeshType :: RealType,
+                                                                                               typename MeshType :: DeviceType,
+                                                                                               typename MeshType :: IndexType > > > >
+                                                        ( parameters );
    }
 }
 
diff --git a/examples/navier-stokes/navierStokesSolver.h b/examples/navier-stokes/navierStokesSolver.h
index 055dddb72fbf720d7f69bda67d28bed68e7026a5..82903ae36329f3ab7f72bb0f3d6864b9590a3d5e 100644
--- a/examples/navier-stokes/navierStokesSolver.h
+++ b/examples/navier-stokes/navierStokesSolver.h
@@ -27,8 +27,10 @@
 #include <solvers/tnlSolverMonitor.h>
 #include "navierStokesSolverMonitor.h"
 #include <schemes/euler/fvm/tnlLaxFridrichs.h>
+#include <schemes/gradient/tnlCentralFDMGradient.h>
 
-template< typename Mesh, typename EulerScheme >
+template< typename Mesh,
+          typename EulerScheme >
 class navierStokesSolver
 {
    public:
diff --git a/examples/navier-stokes/navierStokesSolver_impl.h b/examples/navier-stokes/navierStokesSolver_impl.h
index 6fc4076a237b255db6bec7709c87dc933cdcde3c..420e28a6ab9c88f92f7c9869865c66e90e1b698a 100644
--- a/examples/navier-stokes/navierStokesSolver_impl.h
+++ b/examples/navier-stokes/navierStokesSolver_impl.h
@@ -222,6 +222,27 @@ bool navierStokesSolver< Mesh, EulerScheme > :: makeSnapshot( const RealType& t,
    FileNameBaseNumberEnding( "u-", step, 5, ".tnl", fileName );
    if( ! u. save( fileName ) )
       return false;
+   FileNameBaseNumberEnding( "rho-t-", step, 5, ".tnl", fileName );
+   if( ! rho_t. save( fileName ) )
+      return false;
+   FileNameBaseNumberEnding( "rho-u1-t-", step, 5, ".tnl", fileName );
+   if( ! rho_u1_t. save( fileName ) )
+      return false;
+
+   FileNameBaseNumberEnding( "rho-u2-t-", step, 5, ".tnl", fileName );
+   if( ! rho_u2_t. save( fileName ) )
+      return false;
+
+   FileNameBaseNumberEnding( "rho-", step, 5, ".tnl", fileName );
+   if( ! rho. save( fileName ) )
+      return false;
+   FileNameBaseNumberEnding( "rho-u1-", step, 5, ".tnl", fileName );
+   if( ! rho_u1. save( fileName ) )
+      return false;
+
+   FileNameBaseNumberEnding( "rho-u2-", step, 5, ".tnl", fileName );
+   if( ! rho_u2. save( fileName ) )
+      return false;
    return true;
 }
 
@@ -300,10 +321,16 @@ void navierStokesSolver< Mesh, EulerScheme > :: GetExplicitRHS(  const RealType&
    rho_u1. bind( & u. getData()[ dofs ], dofs );
    rho_u2. bind( & u. getData()[ 2 * dofs ], dofs );
 
+   eulerScheme. setRho( rho );
+   eulerScheme. setRhoU1( rho_u1 );
+   eulerScheme. setRhoU2( rho_u2 );
+
    rho_t. bind( & fu. getData()[ 0 ], dofs );
    rho_u1_t. bind( & fu. getData()[ dofs ], dofs );
    rho_u2_t. bind( & fu. getData()[ 2 * dofs ], dofs );
 
+
+
    updatePhysicalQuantities( rho, rho_u1, rho_u2 );
 
    /****
@@ -407,17 +434,13 @@ void navierStokesSolver< Mesh, EulerScheme > :: GetExplicitRHS(  const RealType&
             const RealType u_sqr = u * u;
             const RealType v_sqr = v * v;
             eulerScheme. getExplicitRhs( c,
-                                         rho,
-                                         rho_u1,
-                                         rho_u2,
-                                         rho_t,
-                                         rho_u1_t,
-                                         rho_u2_t );
+                                         rho_t[ c ],
+                                         rho_u1_t[ c ],
+                                         rho_u2_t[ c ] );
             
-            rho_u1_t[ c ] += -( p[ e ] - p[ w ] ) / ( 2.0 * hx );
-            rho_u2_t[ c ] += -( p[ n ] - p[ s ] ) / ( 2.0 * hy );
+            //rho_u1_t[ c ] += -( p[ e ] - p[ w ] ) / ( 2.0 * hx );
+            rho_u2_t[ c ] += //-( p[ n ] - p[ s ] ) / ( 2.0 * hy );
                              - startUpCoefficient * this -> gravity * this -> rho[ c ];
-                  break;
 
             /***
              * Add the viscosity term
diff --git a/src/implementation/mesh/tnlGrid1D_impl.h b/src/implementation/mesh/tnlGrid1D_impl.h
index 6b43365b79c5b75868639c9ba8bebc64ebfc9250..c9963b22fe21b73adab38b8d01b8fae39623fc77 100644
--- a/src/implementation/mesh/tnlGrid1D_impl.h
+++ b/src/implementation/mesh/tnlGrid1D_impl.h
@@ -26,16 +26,18 @@ using namespace std;
 
 template< typename Real,
           typename Device,
-          typename Index >
-tnlGrid< 1, Real, Device, Index> :: tnlGrid()
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+tnlGrid< 1, Real, Device, Index, Geometry > :: tnlGrid()
 : dofs( 0 )
 {
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-tnlString tnlGrid< 1, Real, Device, Index> :: getTypeStatic()
+          typename Index,
+          template< int, typename, typename, typename > class Geometry  >
+tnlString tnlGrid< 1, Real, Device, Index, Geometry > :: getTypeStatic()
 {
    return tnlString( "tnlGrid< " ) +
           tnlString( Dimensions ) + ", " +
@@ -46,16 +48,18 @@ tnlString tnlGrid< 1, Real, Device, Index> :: getTypeStatic()
 
 template< typename Real,
           typename Device,
-          typename Index >
-tnlString tnlGrid< 1, Real, Device, Index> :: getType() const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry  >
+tnlString tnlGrid< 1, Real, Device, Index, Geometry > :: getType() const
 {
    return this -> getTypeStatic();
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-void tnlGrid< 1, Real, Device, Index> :: setDimensions( const Index xSize )
+          typename Index,
+          template< int, typename, typename, typename > class Geometry  >
+void tnlGrid< 1, Real, Device, Index, Geometry > :: setDimensions( const Index xSize )
 {
    tnlAssert( xSize > 1,
               cerr << "The number of Elements along x-axis must be larger than 1." );
@@ -65,8 +69,9 @@ void tnlGrid< 1, Real, Device, Index> :: setDimensions( const Index xSize )
 
 template< typename Real,
           typename Device,
-          typename Index >
-void tnlGrid< 1, Real, Device, Index> :: setDimensions( const tnlTuple< 1, Index >& dimensions )
+          typename Index,
+          template< int, typename, typename, typename > class Geometry  >
+void tnlGrid< 1, Real, Device, Index, Geometry > :: setDimensions( const tnlTuple< 1, Index >& dimensions )
 {
    tnlAssert( dimensions. x() > 1,
               cerr << "The number of Elements along x-axis must be larger than 1." );
@@ -76,48 +81,54 @@ void tnlGrid< 1, Real, Device, Index> :: setDimensions( const tnlTuple< 1, Index
 
 template< typename Real,
           typename Device,
-          typename Index >
-const tnlTuple< 1, Index >& tnlGrid< 1, Real, Device, Index> :: getDimensions() const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry  >
+const tnlTuple< 1, Index >& tnlGrid< 1, Real, Device, Index, Geometry > :: getDimensions() const
 {
    return this -> dimensions;
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-void tnlGrid< 1, Real, Device, Index> :: setOrigin( const tnlTuple< 1, Real >& origin )
+          typename Index,
+          template< int, typename, typename, typename > class Geometry  >
+void tnlGrid< 1, Real, Device, Index, Geometry > :: setOrigin( const tnlTuple< 1, Real >& origin )
 {
    this -> origin = origin;
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-const tnlTuple< 1, Real >& tnlGrid< 1, Real, Device, Index> :: getOrigin() const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry  >
+const tnlTuple< 1, Real >& tnlGrid< 1, Real, Device, Index, Geometry > :: getOrigin() const
 {
    return this -> origin;
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-void tnlGrid< 1, Real, Device, Index> :: setProportions( const tnlTuple< 1, Real >& proportions )
+          typename Index,
+          template< int, typename, typename, typename > class Geometry  >
+void tnlGrid< 1, Real, Device, Index, Geometry > :: setProportions( const tnlTuple< 1, Real >& proportions )
 {
    this -> proportions = proportions;
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-const tnlTuple< 1, Real >& tnlGrid< 1, Real, Device, Index> :: getProportions() const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry  >
+const tnlTuple< 1, Real >& tnlGrid< 1, Real, Device, Index, Geometry > :: getProportions() const
 {
    return this -> proportions;
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-void tnlGrid< 1, Real, Device, Index> :: setSpaceStep( const tnlTuple< 1, Real >& spaceStep )
+          typename Index,
+          template< int, typename, typename, typename > class Geometry  >
+void tnlGrid< 1, Real, Device, Index, Geometry > :: setSpaceStep( const tnlTuple< 1, Real >& spaceStep )
 {
    this -> proportions. x() = this -> dimensions. x() *
                               spaceStep. x();
@@ -125,8 +136,9 @@ void tnlGrid< 1, Real, Device, Index> :: setSpaceStep( const tnlTuple< 1, Real >
 
 template< typename Real,
           typename Device,
-          typename Index >
-tnlTuple< 1, Real > tnlGrid< 1, Real, Device, Index> :: getSpaceStep() const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry  >
+tnlTuple< 1, Real > tnlGrid< 1, Real, Device, Index, Geometry > :: getSpaceStep() const
 {
    tnlAssert( dimensions. x() > 0,
               cerr << "Cannot get the space step hx since number of Elements along the x axis is not known in tnlGrid "
@@ -140,8 +152,9 @@ tnlTuple< 1, Real > tnlGrid< 1, Real, Device, Index> :: getSpaceStep() const
 
 template< typename Real,
           typename Device,
-          typename Index >
-Index tnlGrid< 1, Real, Device, Index> :: getElementIndex( const Index i ) const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+Index tnlGrid< 1, Real, Device, Index, Geometry > :: getElementIndex( const Index i ) const
 {
    tnlAssert( i < dimensions. x(),
               cerr << "Index i ( " << i
@@ -152,16 +165,18 @@ Index tnlGrid< 1, Real, Device, Index> :: getElementIndex( const Index i ) const
 
 template< typename Real,
           typename Device,
-          typename Index >
-Index tnlGrid< 1, Real, Device, Index> :: getDofs() const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+Index tnlGrid< 1, Real, Device, Index, Geometry > :: getDofs() const
 {
    return this -> dofs;
 };
 
 template< typename Real,
           typename Device,
-          typename Index >
-bool tnlGrid< 1, Real, Device, Index> :: save( tnlFile& file ) const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+bool tnlGrid< 1, Real, Device, Index, Geometry > :: save( tnlFile& file ) const
 {
    if( ! tnlObject :: save( file ) )
       return false;
@@ -178,8 +193,9 @@ bool tnlGrid< 1, Real, Device, Index> :: save( tnlFile& file ) const
 
 template< typename Real,
           typename Device,
-          typename Index >
-bool tnlGrid< 1, Real, Device, Index> :: load( tnlFile& file )
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+bool tnlGrid< 1, Real, Device, Index, Geometry > :: load( tnlFile& file )
 {
    if( ! tnlObject :: load( file ) )
       return false;
@@ -197,25 +213,28 @@ bool tnlGrid< 1, Real, Device, Index> :: load( tnlFile& file )
 
 template< typename Real,
           typename Device,
-          typename Index >
-bool tnlGrid< 1, Real, Device, Index> :: save( const tnlString& fileName ) const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+bool tnlGrid< 1, Real, Device, Index, Geometry > :: save( const tnlString& fileName ) const
 {
    return tnlObject :: save( fileName );
 };
 
 template< typename Real,
           typename Device,
-          typename Index >
-bool tnlGrid< 1, Real, Device, Index> :: load( const tnlString& fileName )
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+bool tnlGrid< 1, Real, Device, Index, Geometry > :: load( const tnlString& fileName )
 {
    return tnlObject :: load( fileName );
 };
 
 template< typename Real,
            typename Device,
-           typename Index >
+           typename Index,
+           template< int, typename, typename, typename > class Geometry >
    template< typename MeshFunction >
-bool tnlGrid< 1, Real, Device, Index> :: write( const MeshFunction& function,
+bool tnlGrid< 1, Real, Device, Index, Geometry > :: write( const MeshFunction& function,
                                                 const tnlString& fileName,
                                                 const tnlString& format ) const
 {
diff --git a/src/implementation/mesh/tnlGrid2D_impl.h b/src/implementation/mesh/tnlGrid2D_impl.h
index 381140d6678b74345b1fc6d8ba365e58e8949275..99503c3515ac815d5ea0f4f426105a30fe69c3bc 100644
--- a/src/implementation/mesh/tnlGrid2D_impl.h
+++ b/src/implementation/mesh/tnlGrid2D_impl.h
@@ -25,16 +25,18 @@ using namespace std;
 
 template< typename Real,
           typename Device,
-          typename Index >
-tnlGrid< 2, Real, Device, Index> :: tnlGrid()
+          typename Index,
+          template< int, typename, typename, typename > class Geometry  >
+tnlGrid< 2, Real, Device, Index, Geometry > :: tnlGrid()
 : dofs( 0 )
 {
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-tnlString tnlGrid< 2, Real, Device, Index> :: getTypeStatic()
+          typename Index,
+          template< int, typename, typename, typename > class Geometry  >
+tnlString tnlGrid< 2, Real, Device, Index, Geometry > :: getTypeStatic()
 {
    return tnlString( "tnlGrid< " ) +
           tnlString( Dimensions ) + ", " +
@@ -45,16 +47,18 @@ tnlString tnlGrid< 2, Real, Device, Index> :: getTypeStatic()
 
 template< typename Real,
            typename Device,
-           typename Index >
-tnlString tnlGrid< 2, Real, Device, Index> :: getType() const
+           typename Index,
+           template< int, typename, typename, typename > class Geometry >
+tnlString tnlGrid< 2, Real, Device, Index, Geometry > :: getType() const
 {
    return this -> getTypeStatic();
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-void tnlGrid< 2, Real, Device, Index> :: setDimensions( const Index ySize, const Index xSize )
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+void tnlGrid< 2, Real, Device, Index, Geometry > :: setDimensions( const Index ySize, const Index xSize )
 {
    tnlAssert( xSize > 1,
               cerr << "The number of Elements along x-axis must be larger than 1." );
@@ -68,8 +72,9 @@ void tnlGrid< 2, Real, Device, Index> :: setDimensions( const Index ySize, const
 
 template< typename Real,
           typename Device,
-          typename Index >
-void tnlGrid< 2, Real, Device, Index> :: setDimensions( const tnlTuple< 2, Index >& dimensions )
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+void tnlGrid< 2, Real, Device, Index, Geometry > :: setDimensions( const tnlTuple< 2, Index >& dimensions )
 {
    tnlAssert( dimensions. x() > 1,
               cerr << "The number of Elements along x-axis must be larger than 1." );
@@ -82,48 +87,54 @@ void tnlGrid< 2, Real, Device, Index> :: setDimensions( const tnlTuple< 2, Index
 
 template< typename Real,
           typename Device,
-          typename Index >
-const tnlTuple< 2, Index >& tnlGrid< 2, Real, Device, Index> :: getDimensions() const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+const tnlTuple< 2, Index >& tnlGrid< 2, Real, Device, Index, Geometry > :: getDimensions() const
 {
    return this -> dimensions;
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-void tnlGrid< 2, Real, Device, Index> :: setOrigin( const tnlTuple< 2, Real >& origin )
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+void tnlGrid< 2, Real, Device, Index, Geometry > :: setOrigin( const tnlTuple< 2, Real >& origin )
 {
    this -> origin = origin;
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-const tnlTuple< 2, Real >& tnlGrid< 2, Real, Device, Index> :: getOrigin() const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+const tnlTuple< 2, Real >& tnlGrid< 2, Real, Device, Index, Geometry > :: getOrigin() const
 {
    return this -> origin;
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-void tnlGrid< 2, Real, Device, Index> :: setProportions( const tnlTuple< 2, Real >& proportions )
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+void tnlGrid< 2, Real, Device, Index, Geometry > :: setProportions( const tnlTuple< 2, Real >& proportions )
 {
    this -> proportions = proportions;
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-const tnlTuple< 2, Real >& tnlGrid< 2, Real, Device, Index> :: getProportions() const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+const tnlTuple< 2, Real >& tnlGrid< 2, Real, Device, Index, Geometry > :: getProportions() const
 {
    return this -> proportions;
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-void tnlGrid< 2, Real, Device, Index> :: setSpaceStep( const tnlTuple< 2, Real >& spaceStep )
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+void tnlGrid< 2, Real, Device, Index, Geometry > :: setSpaceStep( const tnlTuple< 2, Real >& spaceStep )
 {
    this -> proportions. x() = this -> dimensions. x() *
                               spaceStep. x();
@@ -134,8 +145,9 @@ void tnlGrid< 2, Real, Device, Index> :: setSpaceStep( const tnlTuple< 2, Real >
 
 template< typename Real,
           typename Device,
-          typename Index >
-tnlTuple< 2, Real > tnlGrid< 2, Real, Device, Index> :: getSpaceStep() const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+tnlTuple< 2, Real > tnlGrid< 2, Real, Device, Index, Geometry > :: getSpaceStep() const
 {
    tnlAssert( dimensions. x() > 0,
               cerr << "Cannot get the space step hx since number of Elements along the x axis is not known in tnlGrid "
@@ -154,8 +166,9 @@ tnlTuple< 2, Real > tnlGrid< 2, Real, Device, Index> :: getSpaceStep() const
 
 template< typename Real,
           typename Device,
-          typename Index >
-Index tnlGrid< 2, Real, Device, Index> :: getElementIndex( const Index j, const Index i ) const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+Index tnlGrid< 2, Real, Device, Index, Geometry > :: getElementIndex( const Index j, const Index i ) const
 {
    tnlAssert( i < dimensions. x(),
               cerr << "Index i ( " << i
@@ -171,8 +184,9 @@ Index tnlGrid< 2, Real, Device, Index> :: getElementIndex( const Index j, const
 
 template< typename Real,
           typename Device,
-          typename Index >
-Index tnlGrid< 2, Real, Device, Index> :: getElementNeighbour( const Index Element,
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+Index tnlGrid< 2, Real, Device, Index, Geometry > :: getElementNeighbour( const Index Element,
                                                             const Index dy,
                                                             const Index dx ) const
 {
@@ -187,16 +201,18 @@ Index tnlGrid< 2, Real, Device, Index> :: getElementNeighbour( const Index Eleme
 
 template< typename Real,
            typename Device,
-           typename Index >
-Index tnlGrid< 2, Real, Device, Index> :: getDofs() const
+           typename Index,
+           template< int, typename, typename, typename > class Geometry >
+Index tnlGrid< 2, Real, Device, Index, Geometry > :: getDofs() const
 {
    return this -> dofs;
 };
 
 template< typename Real,
           typename Device,
-          typename Index >
-bool tnlGrid< 2, Real, Device, Index> :: save( tnlFile& file ) const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+bool tnlGrid< 2, Real, Device, Index, Geometry > :: save( tnlFile& file ) const
 {
    if( ! tnlObject :: save( file ) )
       return false;
@@ -213,8 +229,9 @@ bool tnlGrid< 2, Real, Device, Index> :: save( tnlFile& file ) const
 
 template< typename Real,
           typename Device,
-          typename Index >
-bool tnlGrid< 2, Real, Device, Index> :: load( tnlFile& file )
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+bool tnlGrid< 2, Real, Device, Index, Geometry > :: load( tnlFile& file )
 {
    if( ! tnlObject :: load( file ) )
       return false;
@@ -233,25 +250,28 @@ bool tnlGrid< 2, Real, Device, Index> :: load( tnlFile& file )
 
 template< typename Real,
           typename Device,
-          typename Index >
-bool tnlGrid< 2, Real, Device, Index> :: save( const tnlString& fileName ) const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+bool tnlGrid< 2, Real, Device, Index, Geometry > :: save( const tnlString& fileName ) const
 {
    return tnlObject :: save( fileName );
 };
 
 template< typename Real,
            typename Device,
-           typename Index >
-bool tnlGrid< 2, Real, Device, Index> :: load( const tnlString& fileName )
+           typename Index,
+           template< int, typename, typename, typename > class Geometry >
+bool tnlGrid< 2, Real, Device, Index, Geometry > :: load( const tnlString& fileName )
 {
    return tnlObject :: load( fileName );
 };
 
 template< typename Real,
            typename Device,
-           typename Index >
+           typename Index,
+           template< int, typename, typename, typename > class Geometry >
    template< typename MeshFunction >
-bool tnlGrid< 2, Real, Device, Index> :: write( const MeshFunction& function,
+bool tnlGrid< 2, Real, Device, Index, Geometry > :: write( const MeshFunction& function,
                                                    const tnlString& fileName,
                                                    const tnlString& format ) const
 {
diff --git a/src/implementation/mesh/tnlGrid3D_impl.h b/src/implementation/mesh/tnlGrid3D_impl.h
index 4f50005bb9e6571e1781f2a42be3542233bff909..c07d9839746b6424394dfbcb85ce6c893c7bf3a0 100644
--- a/src/implementation/mesh/tnlGrid3D_impl.h
+++ b/src/implementation/mesh/tnlGrid3D_impl.h
@@ -22,16 +22,18 @@
 
 template< typename Real,
           typename Device,
-          typename Index >
-tnlGrid< 3, Real, Device, Index> :: tnlGrid()
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+tnlGrid< 3, Real, Device, Index, Geometry > :: tnlGrid()
 : dofs( 0 )
 {
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-tnlString tnlGrid< 3, Real, Device, Index> :: getTypeStatic()
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+tnlString tnlGrid< 3, Real, Device, Index, Geometry > :: getTypeStatic()
 {
    return tnlString( "tnlGrid< " ) +
           tnlString( Dimensions ) + ", " +
@@ -42,16 +44,18 @@ tnlString tnlGrid< 3, Real, Device, Index> :: getTypeStatic()
 
 template< typename Real,
           typename Device,
-          typename Index >
-tnlString tnlGrid< 3, Real, Device, Index> :: getType() const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+tnlString tnlGrid< 3, Real, Device, Index, Geometry > :: getType() const
 {
    return this -> getTypeStatic();
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-void tnlGrid< 3, Real, Device, Index> :: setDimensions( const Index zSize, const Index ySize, const Index xSize )
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+void tnlGrid< 3, Real, Device, Index, Geometry > :: setDimensions( const Index zSize, const Index ySize, const Index xSize )
 {
    tnlAssert( xSize > 1,
               cerr << "The number of Elements along x-axis must be larger than 1." );
@@ -68,8 +72,9 @@ void tnlGrid< 3, Real, Device, Index> :: setDimensions( const Index zSize, const
 
 template< typename Real,
           typename Device,
-          typename Index >
-void tnlGrid< 3, Real, Device, Index> :: setDimensions( const tnlTuple< 3, Index >& dimensions )
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+void tnlGrid< 3, Real, Device, Index, Geometry > :: setDimensions( const tnlTuple< 3, Index >& dimensions )
 {
    tnlAssert( dimensions. x() > 1,
               cerr << "The number of Elements along x-axis must be larger than 1." );
@@ -84,48 +89,54 @@ void tnlGrid< 3, Real, Device, Index> :: setDimensions( const tnlTuple< 3, Index
 
 template< typename Real,
           typename Device,
-          typename Index >
-const tnlTuple< 3, Index >& tnlGrid< 3, Real, Device, Index> :: getDimensions() const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+const tnlTuple< 3, Index >& tnlGrid< 3, Real, Device, Index, Geometry > :: getDimensions() const
 {
    return this -> dimensions;
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-void tnlGrid< 3, Real, Device, Index> :: setOrigin( const tnlTuple< 3, Real >& origin )
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+void tnlGrid< 3, Real, Device, Index, Geometry > :: setOrigin( const tnlTuple< 3, Real >& origin )
 {
    this -> origin = origin;
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-const tnlTuple< 3, Real >& tnlGrid< 3, Real, Device, Index> :: getOrigin() const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+const tnlTuple< 3, Real >& tnlGrid< 3, Real, Device, Index, Geometry > :: getOrigin() const
 {
    return this -> origin;
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-void tnlGrid< 3, Real, Device, Index> :: setProportions( const tnlTuple< 3, Real >& proportions )
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+void tnlGrid< 3, Real, Device, Index, Geometry > :: setProportions( const tnlTuple< 3, Real >& proportions )
 {
    this -> proportions = proportions;
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-const tnlTuple< 3, Real >& tnlGrid< 3, Real, Device, Index> :: getProportions() const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+const tnlTuple< 3, Real >& tnlGrid< 3, Real, Device, Index, Geometry > :: getProportions() const
 {
    return this -> proportions;
 }
 
 template< typename Real,
           typename Device,
-          typename Index >
-void tnlGrid< 3, Real, Device, Index> :: setSpaceStep( const tnlTuple< 3, Real >& spaceStep )
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+void tnlGrid< 3, Real, Device, Index, Geometry > :: setSpaceStep( const tnlTuple< 3, Real >& spaceStep )
 {
    this -> proportions. x() = this -> dimensions. x() *
                               spaceStep. x();
@@ -137,8 +148,9 @@ void tnlGrid< 3, Real, Device, Index> :: setSpaceStep( const tnlTuple< 3, Real >
 
 template< typename Real,
           typename Device,
-          typename Index >
-tnlTuple< 3, Real > tnlGrid< 3, Real, Device, Index> :: getSpaceStep() const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+tnlTuple< 3, Real > tnlGrid< 3, Real, Device, Index, Geometry > :: getSpaceStep() const
 {
    tnlAssert( dimensions. x() > 0,
               cerr << "Cannot get the space step hx since number of Elements along the x axis is not known in tnlGrid "
@@ -163,8 +175,9 @@ tnlTuple< 3, Real > tnlGrid< 3, Real, Device, Index> :: getSpaceStep() const
 
 template< typename Real,
           typename Device,
-          typename Index >
-Index tnlGrid< 3, Real, Device, Index> :: getElementIndex( const Index k, const Index j, const Index i ) const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+Index tnlGrid< 3, Real, Device, Index, Geometry > :: getElementIndex( const Index k, const Index j, const Index i ) const
 {
    tnlAssert( i < dimensions. x(),
               cerr << "Index i ( " << i
@@ -184,16 +197,18 @@ Index tnlGrid< 3, Real, Device, Index> :: getElementIndex( const Index k, const
 
 template< typename Real,
           typename Device,
-          typename Index >
-Index tnlGrid< 3, Real, Device, Index> :: getDofs() const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+Index tnlGrid< 3, Real, Device, Index, Geometry > :: getDofs() const
 {
    return this -> dofs;
 };
 
 template< typename Real,
           typename Device,
-          typename Index >
-bool tnlGrid< 3, Real, Device, Index> :: save( tnlFile& file ) const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+bool tnlGrid< 3, Real, Device, Index, Geometry > :: save( tnlFile& file ) const
 {
    if( ! tnlObject :: save( file ) )
       return false;
@@ -210,8 +225,9 @@ bool tnlGrid< 3, Real, Device, Index> :: save( tnlFile& file ) const
 
 template< typename Real,
           typename Device,
-          typename Index >
-bool tnlGrid< 3, Real, Device, Index> :: load( tnlFile& file )
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+bool tnlGrid< 3, Real, Device, Index, Geometry > :: load( tnlFile& file )
 {
    if( ! tnlObject :: load( file ) )
       return false;
@@ -232,27 +248,30 @@ bool tnlGrid< 3, Real, Device, Index> :: load( tnlFile& file )
 
 template< typename Real,
           typename Device,
-          typename Index >
-bool tnlGrid< 3, Real, Device, Index> :: save( const tnlString& fileName ) const
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+bool tnlGrid< 3, Real, Device, Index, Geometry > :: save( const tnlString& fileName ) const
 {
    return tnlObject :: save( fileName );
 };
 
 template< typename Real,
           typename Device,
-          typename Index >
-bool tnlGrid< 3, Real, Device, Index> :: load( const tnlString& fileName )
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+bool tnlGrid< 3, Real, Device, Index, Geometry > :: load( const tnlString& fileName )
 {
    return tnlObject :: load( fileName );
 };
 
 template< typename Real,
-           typename Device,
-           typename Index >
+          typename Device,
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
    template< typename MeshFunction >
-bool tnlGrid< 3, Real, Device, Index> :: write( const MeshFunction& function,
-                                                   const tnlString& fileName,
-                                                   const tnlString& format ) const
+bool tnlGrid< 3, Real, Device, Index, Geometry > :: write( const MeshFunction& function,
+                                                                                     const tnlString& fileName,
+                                                                                     const tnlString& format ) const
 {
    if( this -> getDofs() != function. getSize() )
    {
diff --git a/src/implementation/schemes/CMakeLists.txt b/src/implementation/schemes/CMakeLists.txt
index 59780c95cbb86b67f4daa193a9b3e7e08423628d..528e0aabe0d082a9ed873793985902c2026fe669 100755
--- a/src/implementation/schemes/CMakeLists.txt
+++ b/src/implementation/schemes/CMakeLists.txt
@@ -1,16 +1,19 @@
 ADD_SUBDIRECTORY( euler )
+ADD_SUBDIRECTORY( gradient )
 
 SET( CURRENT_DIR ${CMAKE_SOURCE_DIR}/src/implementation/schemes )
 
 if( BUILD_CUDA)
       set( tnl_implementation_schemes_CUDA__SOURCES
         ${tnl_implementation_schemes_euler_CUDA__SOURCES}
+        ${tnl_implementation_schemes_gradient_CUDA__SOURCES}
         ${common_SOURCES}
         PARENT_SCOPE )
 endif()
 
 set( tnl_implementation_schemes_SOURCES
      ${tnl_implementation_schemes_euler_SOURCES}
+     ${tnl_implementation_schemes_gradient_SOURCES}
      ${common_SOURCES}
      PARENT_SCOPE )
    
diff --git a/src/implementation/schemes/euler/fvm/tnlLaxFridrichs_impl.h b/src/implementation/schemes/euler/fvm/tnlLaxFridrichs_impl.h
index 6b1ce8d6feea12df0399504f65e7abdb7a8ba2d0..472db7edeab9e507ee1707b0f39c7d71d36633fc 100644
--- a/src/implementation/schemes/euler/fvm/tnlLaxFridrichs_impl.h
+++ b/src/implementation/schemes/euler/fvm/tnlLaxFridrichs_impl.h
@@ -19,44 +19,98 @@
 #define TNLLAXFRIDRICHS_IMPL_H_
 
 
-template< typename MeshType >
-tnlLaxFridrichs< MeshType > :: tnlLaxFridrichs()
+template< typename Real,
+          typename Device,
+          typename Index >
+template< typename PressureGradient >
+tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >,
+                 PressureGradient > :: tnlLaxFridrichs()
 : regularizeEps( 0.0 ),
   viscosityCoefficient( 1.0 ),
-  mesh( 0 )
+  mesh( 0 ),
+  pressureGradient( 0 )
 {
 }
 
-template< typename MeshType >
-void tnlLaxFridrichs< MeshType > :: bindMesh( const MeshType& mesh )
+template< typename Real,
+          typename Device,
+          typename Index >
+template< typename PressureGradient >
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient  > :: bindMesh( const MeshType& mesh )
 {
    this -> mesh = &mesh;
 }
 
 
-template< typename MeshType >
-void tnlLaxFridrichs< MeshType > :: setRegularization( const RealType& epsilon )
+template< typename Real,
+          typename Device,
+          typename Index >
+template< typename PressureGradient >
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient  > :: setRegularization( const RealType& epsilon )
 {
    this -> regularizeEps = epsilon;
 }
 
-template< typename MeshType >
-void tnlLaxFridrichs< MeshType > :: setViscosityCoefficient( const RealType& v )
+template< typename Real,
+          typename Device,
+          typename Index >
+template< typename PressureGradient >
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient  > :: setViscosityCoefficient( const RealType& v )
 {
    this -> viscosityCoefficient = v;
 }
 
-template< typename MeshType >
+template< typename Real,
+          typename Device,
+          typename Index >
+template< typename PressureGradient >
    template< typename Vector >
-void tnlLaxFridrichs< MeshType > :: getExplicitRhs( const IndexType centralVolume,
-                                                    const Vector& rho,
-                                                    const Vector& rho_u1,
-                                                    const Vector& rho_u2,
-                                                    Vector& rho_t,
-                                                    Vector& rho_u1_t,
-                                                    Vector& rho_u2_t ) const
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient  > :: setRho( Vector& rho )
 {
-   tnlAssert( mesh, cerr << "No mesh has been binded with the Lax-Fridrichs scheme.");
+   this -> rho. bind( rho );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+template< typename PressureGradient >
+   template< typename Vector >
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient  > :: setRhoU1( Vector& rho_u1 )
+{
+   this -> rho_u1. bind( rho_u1 );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+template< typename PressureGradient >
+   template< typename Vector >
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient  > :: setRhoU2( Vector& rho_u2 )
+{
+   this -> rho_u2. bind( rho_u2 );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+template< typename PressureGradient >
+   template< typename Vector >
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient  > :: setPressureGradient( Vector& grad_p )
+{
+   this -> pressureGradient = &grad_p;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+template< typename PressureGradient >
+void tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient  > :: getExplicitRhs( const IndexType centralVolume,
+                                                                      RealType& rho_t,
+                                                                      RealType& rho_u1_t,
+                                                                      RealType& rho_u2_t ) const
+{
+   tnlAssert( mesh, cerr << "No mesh has been binded with the Lax-Fridrichs scheme." );
+   tnlAssert( pressureGradient, cerr << "No pressure gradient was set in the the Lax-Fridrichs scheme." )
 
    const IndexType& xSize = this -> mesh -> getDimensions(). x();
    const IndexType& ySize = this -> mesh -> getDimensions(). y();
@@ -76,23 +130,37 @@ void tnlLaxFridrichs< MeshType > :: getExplicitRhs( const IndexType centralVolum
    const RealType u1_w = rho_u1[ w ] / regularize( rho[ w ] );
    const RealType u2_n = rho_u2[ n ] / regularize( rho[ n ] );
    const RealType u2_s = rho_u2[ s ] / regularize( rho[ s ] );
-   rho_t[ c ]= this -> viscosityCoefficient * 0.25 * ( rho[ e ] + rho[ w ] + rho[ s ] + rho[ n ] - 4.0 * rho[ c ] )
-                - ( rho[ e ] * u1_e - rho[ w ] * u1_w ) / ( 2.0 * hx )
-                - ( rho[ n ] * u2_n - rho[ s ] * u2_s ) / ( 2.0 * hy );
-
-    /****
-     * ( rho * u1 )_t + ( rho * u1 * u1 )_x + ( rho * u1 * u2 )_y =  0
-     */
-    rho_u1_t[ c ] = this -> viscosityCoefficient * 0.25 * ( rho_u1[ e ] + rho_u1[ w ] + rho_u1[ s ] + rho_u1[ n ] - 4.0 * rho_u1[ c ] )
-                    - ( rho_u1[ e ] * u1_e - rho_u1[ w ] * u1_w ) / ( 2.0 * hx )
-                    - ( rho_u1[ n ] * u2_n - rho_u1[ s ] * u2_s ) / ( 2.0 * hy );
-    rho_u2_t[ c ] = this -> viscosityCoefficient * 0.25 * ( rho_u2[ e ] + rho_u2[ w ] + rho_u2[ s ] + rho_u2[ n ] - 4.0 * rho_u2[ c ] )
-                    - ( rho_u2[ e ] * u1_e - rho_u2[ w ] * u1_w ) / ( 2.0 * hx )
-                    - ( rho_u2[ n ] * u2_n - rho_u2[ s ] * u2_s ) / ( 2.0 * hy );
+   rho_t= this -> viscosityCoefficient * 0.25 * ( rho[ e ] + rho[ w ] + rho[ s ] + rho[ n ] - 4.0 * rho[ c ] )
+               - ( rho[ e ] * u1_e - rho[ w ] * u1_w ) / ( 2.0 * hx )
+               - ( rho[ n ] * u2_n - rho[ s ] * u2_s ) / ( 2.0 * hy );
+
+   /****
+    * Compute the pressure gradient
+    */
+   RealType p_x, p_y;
+   pressureGradient -> getGradient( c, p_x, p_y );
+
+   /****
+    * ( rho * u1 )_t + ( rho * u1 * u1 )_x + ( rho * u1 * u2 )_y - p_x =  0
+    */
+   rho_u1_t = this -> viscosityCoefficient * 0.25 * ( rho_u1[ e ] + rho_u1[ w ] + rho_u1[ s ] + rho_u1[ n ] - 4.0 * rho_u1[ c ] )
+                   - ( rho_u1[ e ] * u1_e - rho_u1[ w ] * u1_w ) / ( 2.0 * hx )
+                   - ( rho_u1[ n ] * u2_n - rho_u1[ s ] * u2_s ) / ( 2.0 * hy )
+                   - p_x;
+   /****
+    * ( rho * u1 )_t + ( rho * u1 * u1 )_x + ( rho * u1 * u2 )_y - p_y =  0
+    */
+   rho_u2_t = this -> viscosityCoefficient * 0.25 * ( rho_u2[ e ] + rho_u2[ w ] + rho_u2[ s ] + rho_u2[ n ] - 4.0 * rho_u2[ c ] )
+                   - ( rho_u2[ e ] * u1_e - rho_u2[ w ] * u1_w ) / ( 2.0 * hx )
+                   - ( rho_u2[ n ] * u2_n - rho_u2[ s ] * u2_s ) / ( 2.0 * hy )
+                   - p_y;
 }
 
-template< typename MeshType >
-typename MeshType :: RealType tnlLaxFridrichs< MeshType > :: regularize( const RealType& r ) const
+template< typename Real,
+          typename Device,
+          typename Index >
+template< typename PressureGradient >
+typename MeshType :: RealType tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient  > :: regularize( const RealType& r ) const
 {
    return r + ( ( r >= 0 ) - ( r < 0 ) ) * this -> regularizeEps;
 }
diff --git a/src/implementation/schemes/gradient/CMakeLists.txt b/src/implementation/schemes/gradient/CMakeLists.txt
new file mode 100755
index 0000000000000000000000000000000000000000..145da4a3371a7f6c7e451b47d166717de2540252
--- /dev/null
+++ b/src/implementation/schemes/gradient/CMakeLists.txt
@@ -0,0 +1,15 @@
+set( headers tnlCentralFDMGradient_impl.h )
+
+SET( CURRENT_DIR ${CMAKE_SOURCE_DIR}/src/implementation/schemes/gradient )
+
+if( BUILD_CUDA)
+      set( tnl_implementation_schemes_gradient_CUDA__SOURCES
+        ${common_SOURCES}
+        PARENT_SCOPE )
+endif()
+
+set( tnl_implementation_schemes_gradient_SOURCES
+     ${common_SOURCES}
+     PARENT_SCOPE )
+
+INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/implementation/schemes/gradient )
diff --git a/src/implementation/schemes/gradient/tnlCentralFDMGradient_impl.h b/src/implementation/schemes/gradient/tnlCentralFDMGradient_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..1f964ae44a77eb3f8f87e0684dd42039fad1f560
--- /dev/null
+++ b/src/implementation/schemes/gradient/tnlCentralFDMGradient_impl.h
@@ -0,0 +1,56 @@
+/***************************************************************************
+                          tnlCentralFDMGradient_impl.h  -  description
+                             -------------------
+    begin                : Apr 26, 2013
+    copyright            : (C) 2013 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 TNLCENTRALFDMGRADIENT_IMPL_H_
+#define TNLCENTRALFDMGRADIENT_IMPL_H_
+
+template< typename Real, typename Device, typename Index >
+tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index > > :: tnlCentralFDMGradient()
+: mesh( 0 )
+{
+}
+
+template< typename Real, typename Device, typename Index >
+void tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index > > :: bindMesh( const tnlGrid< 2, RealType, DeviceType, IndexType >& mesh )
+{
+   this -> mesh = mesh;
+}
+
+template< typename Real, typename Device, typename Index >
+   template< typename Vector >
+void tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index > > :: setFunction( Vector& f )
+{
+   this -> f. bind( f );
+}
+
+template< typename Real, typename Device, typename Index >
+void tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index > > :: getGradient( const Index& i,
+                                                                     RealType& f_x,
+                                                                     RealType& f_y ) const
+{
+   tnlAssert( this -> mesh, cerr << "No mesh was set in tnlCentralFDMGradient" );
+
+   const Index e = mesh -> getElementNeighbour( i,  1,  0 );
+   const Index w = mesh -> getElementNeighbour( i, -1,  0 );
+   const Index n = mesh -> getElementNeighbour( i,  0,  1 );
+   const Index s = mesh -> getElementNeighbour( i,  0, -1 );
+
+   f_x = ( f[ e ] - f[ w ] ) / ( 2.0 * mesh -> getSpaceStep(). x() );
+   f_y = ( f[ n ] - f[ s ] ) / ( 2.0 * mesh -> getSpaceStep(). y() );
+}
+
+#endif
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index f16bcfbc4b56ef59a5663680a843b60b3a4f790e..54e53cd107a26992972d34188255e20050d306ff 100755
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -8,6 +8,7 @@ ADD_SUBDIRECTORY( io )
 ADD_SUBDIRECTORY( topology )
 
 SET( headers tnlGrid.h
+             tnlIdenticalGridGeometry.h
              tnlDummyMesh.h
              config.h 
              information.h 
diff --git a/src/mesh/tnlGrid.h b/src/mesh/tnlGrid.h
index 956d2b3ef1881c5da0ac564aa7c73662e2a2260d..bea0d287d147198d8c0aff72f2799593f345f2f5 100644
--- a/src/mesh/tnlGrid.h
+++ b/src/mesh/tnlGrid.h
@@ -21,19 +21,22 @@
 #include <core/tnlObject.h>
 #include <core/tnlHost.h>
 #include <core/tnlTuple.h>
+#include <mesh/tnlIdenticalGridGeometry.h>
 
 template< int Dimensions,
           typename Real = double,
           typename Device = tnlHost,
-          typename Index = int >
+          typename Index = int,
+          template< int, typename, typename, typename > class Geometry = tnlIdenticalGridGeometry >
 class tnlGrid : public tnlObject
 {
 };
 
 template< typename Real,
           typename Device,
-          typename Index >
-class tnlGrid< 1, Real, Device, Index> : public tnlObject
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+class tnlGrid< 1, Real, Device, Index, Geometry > : public tnlObject
 {
    public:
 
@@ -97,8 +100,9 @@ class tnlGrid< 1, Real, Device, Index> : public tnlObject
 
 template< typename Real,
           typename Device,
-          typename Index >
-class tnlGrid< 2, Real, Device, Index> : public tnlObject
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+class tnlGrid< 2, Real, Device, Index, Geometry > : public tnlObject
 {
    public:
 
@@ -134,8 +138,8 @@ class tnlGrid< 2, Real, Device, Index> : public tnlObject
    Index getElementIndex( const Index j, const Index i ) const;
 
    Index getElementNeighbour( const Index Element,
-                           const Index dy,
-                           const Index dx ) const;
+                              const Index dy,
+                              const Index dx ) const;
 
    Index getDofs() const;
 
@@ -166,8 +170,9 @@ class tnlGrid< 2, Real, Device, Index> : public tnlObject
 
 template< typename Real,
           typename Device,
-          typename Index >
-class tnlGrid< 3, Real, Device, Index> : public tnlObject
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+class tnlGrid< 3, Real, Device, Index, Geometry > : public tnlObject
 {
    public:
 
diff --git a/src/mesh/tnlIdenticalGridGeometry.h b/src/mesh/tnlIdenticalGridGeometry.h
new file mode 100644
index 0000000000000000000000000000000000000000..73b68533563106841ac5e1bac2ce64fb32b37ea3
--- /dev/null
+++ b/src/mesh/tnlIdenticalGridGeometry.h
@@ -0,0 +1,38 @@
+/***************************************************************************
+                          tnlIdenticalGridGeometry.h  -  description
+                             -------------------
+    begin                : Apr 28, 2013
+    copyright            : (C) 2013 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 TNLIDENTICALGRIDGEOMETRY_H_
+#define TNLIDENTICALGRIDGEOMETRY_H_
+
+#include <core/tnlHost.h>
+
+template< int Dimensions,
+          typename Real = double,
+          typename Device = tnlHost,
+          typename Index = int >
+class tnlIdenticalGridGeometry
+{
+   public:
+
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+
+
+};
+
+#endif /* TNLIDENTICALGRIDGEOMETRY_H_ */
diff --git a/src/schemes/CMakeLists.txt b/src/schemes/CMakeLists.txt
index 90ee7741e321e9a27b8fde2fe6353ba43bd3b798..60a6021346d5d64ab6d7fa43cfca17980a97a396 100755
--- a/src/schemes/CMakeLists.txt
+++ b/src/schemes/CMakeLists.txt
@@ -1 +1,2 @@
 ADD_SUBDIRECTORY( euler )
+ADD_SUBDIRECTORY( gradient )
diff --git a/src/schemes/euler/fvm/tnlLaxFridrichs.h b/src/schemes/euler/fvm/tnlLaxFridrichs.h
index d89a535a27c4864ff759d8c6d4b9ccb102a37b98..fc1c7f011c68bcd975a1a68a240ba1d82b1fd86c 100644
--- a/src/schemes/euler/fvm/tnlLaxFridrichs.h
+++ b/src/schemes/euler/fvm/tnlLaxFridrichs.h
@@ -18,25 +18,35 @@
 #ifndef TNLLAXFRIDRICHS_H_
 #define TNLLAXFRIDRICHS_H_
 
-template< typename MeshType >
+#include <core/tnlSharedVector.h>
+#include <mesh/tnlGrid.h>
+#include <schemes/gradient/tnlCentralFDMGradient.h>
+
+template< typename MeshType,
+          typename PressureGradient >
 class tnlLaxFridrichs
+{
+};
+
+template< typename Real,
+          typename Device,
+          typename Index >
+template< typename PressureGradient >
+class tnlLaxFridrichs< tnlGrid< 2, Real, Device, Index >, PressureGradient >
 {
    public:
 
-   typedef typename MeshType :: RealType RealType;
-   typedef typename MeshType :: DeviceType DeviceType;
-   typedef typename MeshType :: IndexType IndexType;
+   typedef tnlGrid< 2, Real, Device, Index > MeshType;
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
 
    tnlLaxFridrichs();
 
-   template< typename Vector >
    void getExplicitRhs( const IndexType centralVolume,
-                        const Vector& rho,
-                        const Vector& rho_u1,
-                        const Vector& rho_u2,
-                        Vector& rho_t,
-                        Vector& rho_u1_t,
-                        Vector& rho_u2_t ) const;
+                        RealType& rho_t,
+                        RealType& rho_u1_t,
+                        RealType& rho_u2_t ) const;
 
    void setRegularization( const RealType& epsilon );
 
@@ -44,6 +54,18 @@ class tnlLaxFridrichs
 
    void bindMesh( const MeshType& mesh );
 
+   template< typename Vector >
+   void setRho( Vector& rho ); // TODO: add const
+
+   template< typename Vector >
+   void setRhoU1( Vector& rho_u1 ); // TODO: add const
+
+   template< typename Vector >
+   void setRhoU2( Vector& rho_u2 ); // TODO: add const
+
+   template< typename Vector >
+   void setPressureGradient( Vector& grad_p ); // TODO: add const
+
    protected:
 
    RealType regularize( const RealType& r ) const;
@@ -51,6 +73,10 @@ class tnlLaxFridrichs
    RealType regularizeEps, viscosityCoefficient;
 
    const MeshType* mesh;
+
+   const PressureGradient* pressureGradient;
+
+   tnlSharedVector< RealType, DeviceType, IndexType > rho, rho_u1, rho_u2;
 };
 
 #include <implementation/schemes/euler/fvm/tnlLaxFridrichs_impl.h>
diff --git a/src/schemes/gradient/CMakeLists.txt b/src/schemes/gradient/CMakeLists.txt
new file mode 100755
index 0000000000000000000000000000000000000000..bac1e377a48b3d9fda916691bd0c8e4190eab612
--- /dev/null
+++ b/src/schemes/gradient/CMakeLists.txt
@@ -0,0 +1,4 @@
+SET( headers tnlCentralFDMGradient.h
+   )
+   
+INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/schemes/gradient )
diff --git a/src/schemes/gradient/tnlCentralFDMGradient.h b/src/schemes/gradient/tnlCentralFDMGradient.h
new file mode 100644
index 0000000000000000000000000000000000000000..e77e5b7f6352692aeda58f7407cc39ed9a09d2d6
--- /dev/null
+++ b/src/schemes/gradient/tnlCentralFDMGradient.h
@@ -0,0 +1,58 @@
+/***************************************************************************
+                          tnlCentralFDMGradient.h  -  description
+                             -------------------
+    begin                : Apr 26, 2013
+    copyright            : (C) 2013 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 TNLCENTRALFDMGRADIENT_H_
+#define TNLCENTRALFDMGRADIENT_H_
+
+#include <mesh/tnlGrid.h>
+#include <core/tnlHost.h>
+
+template< typename Mesh >
+class tnlCentralFDMGradient
+{
+};
+
+template< typename Real, typename Device, typename Index >
+class tnlCentralFDMGradient< tnlGrid< 2, Real, Device, Index > >
+{
+   public:
+
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+
+   tnlCentralFDMGradient();
+
+   void bindMesh( const tnlGrid< 2, RealType, DeviceType, IndexType >& mesh );
+
+   template< typename Vector >
+   void setFunction( Vector& f ); // TODO: add const
+
+   void getGradient( const Index& i,
+                     RealType& f_x,
+                     RealType& f_y ) const;
+   protected:
+
+   // TODO: change to ConstSharedVector
+   tnlSharedVector< RealType, DeviceType, IndexType > f;
+
+   const tnlGrid< 2, RealType, DeviceType, IndexType >* mesh;
+};
+
+#include <implementation/schemes/gradient/tnlCentralFDMGradient_impl.h>
+
+#endif