diff --git a/examples/navier-stokes/CMakeLists.txt b/examples/navier-stokes/CMakeLists.txt
index 06a2a8cdd2a9902a04edf101a8280f3fb20053e1..0cac83dd5eb8c39759e0ecad4ae5856c9370885a 100755
--- a/examples/navier-stokes/CMakeLists.txt
+++ b/examples/navier-stokes/CMakeLists.txt
@@ -1,10 +1 @@
 add_subdirectory( share )
-
-INSTALL( FILES Makefile
-               main.cpp
-               navierStokesSolver.h
-               navierStokesSolver_impl.h
-               navierStokesSetter.h
-               navierStokesSetter_impl.h
-               navier-stokes.cfg.desc
-         DESTINATION share/tnl-${tnlVersion}/examples/navier-stokes )
\ No newline at end of file
diff --git a/examples/navier-stokes/navierStokesSolver.h b/examples/navier-stokes/navierStokesSolver.h
index e112a5bc1bb8794d23e603b138f0052aa0ba6143..94aae01378bcca5b99dcedc1fc41c89e8dd53662 100644
--- a/examples/navier-stokes/navierStokesSolver.h
+++ b/examples/navier-stokes/navierStokesSolver.h
@@ -28,6 +28,7 @@
 #include "navierStokesSolverMonitor.h"
 #include <schemes/euler/fvm/tnlLaxFridrichs.h>
 #include <schemes/gradient/tnlCentralFDMGradient.h>
+#include <schemes/diffusion/tnlLinearDiffusion.h>
 
 template< typename Mesh,
           typename EulerScheme >
@@ -97,6 +98,8 @@ class navierStokesSolver
 
    EulerScheme eulerScheme;
 
+   tnlLinearDiffusion< MeshType > u1Viscosity, u2Viscosity;
+
    tnlCentralFDMGradient< MeshType > pressureGradient;
 
    navierStokesSolverMonitor< RealType, IndexType > solverMonitor;
diff --git a/examples/navier-stokes/navierStokesSolver_impl.h b/examples/navier-stokes/navierStokesSolver_impl.h
index 4f564f1cad7bf356345d3bc0556bcfbbf0b7d6b4..fed60964a307b211872866c028276ff2c2efe591 100644
--- a/examples/navier-stokes/navierStokesSolver_impl.h
+++ b/examples/navier-stokes/navierStokesSolver_impl.h
@@ -148,7 +148,10 @@ bool navierStokesSolver< Mesh, EulerScheme > :: init( const tnlParameterContaine
    pressureGradient. bindMesh( this -> mesh );
    this -> eulerScheme. bindMesh( this -> mesh );
    this -> eulerScheme. setPressureGradient( this -> pressureGradient );
-
+   this -> u1Viscosity. bindMesh( this -> mesh );
+   this -> u1Viscosity. setFunction( this -> u1 );
+   this -> u2Viscosity. bindMesh( this -> mesh );
+   this -> u2Viscosity. setFunction( this -> u2 );
    return true;
 }
 
@@ -427,33 +430,19 @@ void navierStokesSolver< Mesh, EulerScheme > :: GetExplicitRHS(  const RealType&
                continue;
             }
 
-            IndexType e = mesh. getElementIndex( j, i + 1 );
-            IndexType w = mesh. getElementIndex( j, i - 1 );
-            IndexType n = mesh. getElementIndex( j + 1, i );
-            IndexType s = mesh. getElementIndex( j - 1, i );
-
-            //const RealType& u = this -> u1[ c ];
-            //const RealType& v = this -> u2[ c ];
-            //const RealType u_sqr = u * u;
-            //const RealType v_sqr = v * v;
             eulerScheme. getExplicitRhs( c,
                                          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 );
-                             //- startUpCoefficient * this -> gravity * this -> rho[ c ];
+            //rho_u1_t[ c ] += ;
+            rho_u2_t[ c ] -= startUpCoefficient * this -> gravity * this -> rho[ c ];
 
             /***
              * Add the viscosity term
              */
-            rho_u1_t[ c ] += this -> mu *
-                                   ( ( u1[ e ] - 2.0 * u1[ c ] + u1[ w ] ) / ( hx * hx ) +
-                                     ( u1[ n ] - 2.0 * u1[ c ] + u1[ s ] ) / ( hy * hy ) );
-            rho_u2_t[ c ] += this -> mu *
-                                   ( ( u2[ e ] - 2.0 * u2[ c ] + u2[ w ] ) / ( hx * hx ) +
-                                     ( u2[ n ] - 2.0 * u2[ c ] + u2[ s ] ) / ( hy * hy ) );
+            rho_u1_t[ c ] += this -> mu * u1Viscosity. getDiffusion( c );
+            rho_u2_t[ c ] += this -> mu * u2Viscosity. getDiffusion( c );
 
          }
    }
diff --git a/src/implementation/mesh/CMakeLists.txt b/src/implementation/mesh/CMakeLists.txt
index c9615e4dd0bcde87ccb5069838319ccfb8b89ea6..d91af8415a3596272f224115c11914db95adabcb 100755
--- a/src/implementation/mesh/CMakeLists.txt
+++ b/src/implementation/mesh/CMakeLists.txt
@@ -1,9 +1,10 @@
 SET( headers tnlGrid1D_impl.h
              tnlGrid2D_impl.h
-             tnlGrid3D_impl.h  )
+             tnlGrid3D_impl.h
+             tnlIdenticalGridGeometry_impl.h )
 
 SET( CURRENT_DIR ${CMAKE_SOURCE_DIR}/src/implementation/mesh )    
-set( tnl_implementation_mesh_SOURCES
+SET( tnl_implementation_mesh_SOURCES
      ${CURRENT_DIR}/tnlGrid_impl.cpp
      PARENT_SCOPE )
 
diff --git a/src/implementation/mesh/tnlGrid2D_impl.h b/src/implementation/mesh/tnlGrid2D_impl.h
index 99503c3515ac815d5ea0f4f426105a30fe69c3bc..bafa78139a195b849d88358a79f2b88fc4e463e6 100644
--- a/src/implementation/mesh/tnlGrid2D_impl.h
+++ b/src/implementation/mesh/tnlGrid2D_impl.h
@@ -68,6 +68,10 @@ void tnlGrid< 2, Real, Device, Index, Geometry > :: setDimensions( const Index y
    this -> dimensions. x() = xSize;
    this -> dimensions. y() = ySize;
    dofs = ySize * xSize;
+   tnlTuple< 2, Real > parametricStep;
+   parametricStep. x() = proportions. x() / ( xSize - 1 );
+   parametricStep. y() = proportions. y() / ( ySize - 1 );
+   geometry. setParametricStep( parametricStep );
 }
 
 template< typename Real,
@@ -76,13 +80,7 @@ template< typename Real,
           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." );
-   tnlAssert( dimensions. y() > 1,
-              cerr << "The number of Elements along y-axis must be larger than 1." );
-
-   this -> dimensions = dimensions;
-   dofs = this -> dimensions. x() * this -> dimensions. y();
+   this -> setDimensions( dimensions. y(), dimensions. x() );
 }
 
 template< typename Real,
@@ -119,6 +117,10 @@ template< typename Real,
 void tnlGrid< 2, Real, Device, Index, Geometry > :: setProportions( const tnlTuple< 2, Real >& proportions )
 {
    this -> proportions = proportions;
+   tnlTuple< 2, Real > parametricStep;
+   parametricStep. x() = proportions. x() / ( this -> dimensions. x() - 1 );
+   parametricStep. y() = proportions. y() / ( this -> dimensions. y() - 1 );
+   geometry. setParametricStep( parametricStep );
 }
 
 template< typename Real,
@@ -134,20 +136,19 @@ template< typename Real,
           typename Device,
           typename Index,
           template< int, typename, typename, typename > class Geometry >
-void tnlGrid< 2, Real, Device, Index, Geometry > :: setSpaceStep( const tnlTuple< 2, Real >& spaceStep )
+void tnlGrid< 2, Real, Device, Index, Geometry > :: setParametricStep( const tnlTuple< 2, Real >& spaceStep )
 {
    this -> proportions. x() = this -> dimensions. x() *
-                              spaceStep. x();
+                              parametricStep. x();
    this -> proportions. y() = this -> dimensions. y() *
-                              spaceStep. y();
-
+                              parametricStep. y();
 }
 
 template< typename Real,
           typename Device,
           typename Index,
           template< int, typename, typename, typename > class Geometry >
-tnlTuple< 2, Real > tnlGrid< 2, Real, Device, Index, Geometry > :: getSpaceStep() const
+tnlTuple< 2, Real > tnlGrid< 2, Real, Device, Index, Geometry > :: getParametricStep() 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 "
@@ -156,12 +157,12 @@ tnlTuple< 2, Real > tnlGrid< 2, Real, Device, Index, Geometry > :: getSpaceStep(
               cerr << "Cannot get the space step hy since number of Elements along the y axis is not known in tnlGrid "
                    << this -> getName() );
 
-   tnlTuple< 2, RealType > spaceStep;
-   spaceStep. x() =
+   tnlTuple< 2, RealType > parametricStep;
+   parametricStep. x() =
             this -> proportions. x() / ( Real ) ( this -> dimensions. x() - 1 );
-   spaceStep. y() =
+   parametricStep. y() =
             this -> proportions. y() / ( Real ) ( this -> dimensions. y() - 1 );
-   return spaceStep;
+   return parametricStep;
 }
 
 template< typename Real,
@@ -208,6 +209,50 @@ Index tnlGrid< 2, Real, Device, Index, Geometry > :: getDofs() const
    return this -> dofs;
 };
 
+template< typename Real,
+          typename Device,
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+Real tnlGrid< 2, Real, Device, Index, Geometry > :: getElementMeasure( const Index j,
+                                                                       const Index i ) const
+{
+   return geometry. getElementMeasure( j, i );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+Real tnlGrid< 2, Real, Device, Index, Geometry > :: getElementsDistance( const Index j,
+                                                                         const Index i,
+                                                                         const Index dy,
+                                                                         const Index dx ) const
+{
+   return geometry. getElementsDistance( j, i, dy, dx );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+template< int dy, int dx >
+Real tnlGrid< 2, Real, Device, Index, Geometry > :: getEdgeLength( const Index j,
+                                                                   const Index i ) const
+{
+   return geometry. getEdgeLength< dy, dx >( j, i );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+template< int dy, int dx >
+tnlTuple< 2, Real > tnlGrid< 2, Real, Device, Index, Geometry > :: getEdgeNormal( const Index j,
+                                                                                  const Index i ) const
+{
+   return geometry. getEdgeNormal< dy, dx >( j, i );
+}
+
 template< typename Real,
           typename Device,
           typename Index,
@@ -224,6 +269,8 @@ bool tnlGrid< 2, Real, Device, Index, Geometry > :: save( tnlFile& file ) const
            << this -> getName() << endl;
       return false;
    }
+   if( ! geometry. save( file ) )
+      return false;
    return true;
 };
 
@@ -243,8 +290,14 @@ bool tnlGrid< 2, Real, Device, Index, Geometry > :: load( tnlFile& file )
            << this -> getName() << endl;
       return false;
    }
+   if( ! geometry. load( file ) )
+      return false;
    this -> dofs = this -> getDimensions(). x() *
-                   this -> getDimensions(). y();
+                  this -> getDimensions(). y();
+   tnlTuple< 2, Real > parametricStep;
+   parametricStep. x() = proportions. x() / ( this -> dimensions. x() - 1 );
+   parametricStep. y() = proportions. y() / ( this -> dimensions. y() - 1 );
+   geometry. setParametricStep( parametricStep );
    return true;
 };
 
diff --git a/src/implementation/mesh/tnlIdenticalGridGeometry_impl.h b/src/implementation/mesh/tnlIdenticalGridGeometry_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..0003d8036cfe3b4059b6018b595e09ecac4b4ccb
--- /dev/null
+++ b/src/implementation/mesh/tnlIdenticalGridGeometry_impl.h
@@ -0,0 +1,126 @@
+/***************************************************************************
+                          tnlIdenticalGridGeometry_impl.h  -  description
+                             -------------------
+    begin                : May 1, 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_IMPL_H_
+#define TNLIDENTICALGRIDGEOMETRY_IMPL_H_
+
+#include <core/tnlFile.h>
+#include <core/tnlAssert.h>
+
+template< typename Real,
+          typename Device,
+          typename Index >
+void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: setParametricStep( const tnlTuple< 2, Real >& parametricStep )
+{
+   this -> parametricStep = parametricStep;
+   this -> elementMeasure - this -> parametricStep. x() * this -> parametricStep. y();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementCoordinates( const Index j,
+                                                                                  const Index i,
+                                                                                  tnlTuple< 2, Real >& coordinates ) const
+{
+   coordinates. x() = i * parametricStep. x();
+   coordinates. y() = j * parametricStep. y();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementMeasure( const Index j,
+                                                                              const Index i ) const
+{
+   return elementMeasure;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+   template< Index dy, Index dx >
+Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementsDistance( const Index j,
+                                                                                const Index i ) const
+{
+   if( dy == 0 && dx == 1 )
+      return parametricStep. x();
+   if( dy == 1 && dx == 0 )
+      return parametricStep. y();
+   const Real x = dx * parametricStep. x();
+   const Real y = dy * parametricStep. y();
+   return sqrt( dx * dx + dy * dy );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+template< Index dy, Index dx >
+void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getEdgeCoordinates( const Index j,
+                                                                               const Index i,
+                                                                               tnlTuple< 2, Real >& coordinates ) const
+{
+   coordinates. x() = origin. x() + ( i + 0.5 * ( 1.0 + dx ) ) * parametricStep. x();
+   coordinates. y() = origin. y() + ( j + 0.5 * ( 1.0 + dy ) ) * parametricStep. y();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+template< int dy, int dx >
+Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getEdgeLength( const Index j,
+                                                                          const Index i ) const
+{
+   if( dy == 0 && dx == 1 )
+      return parametricStep. y();
+   if( dy == 1 && dx == 0 )
+      return parametricStep. x();
+   tnlAssert( false, cerr << "Bad values of dx and dy - dx = " << dx << " dy = " << dy );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+template< int dy, int dx >
+void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getEdgeNormal( const Index j,
+                                                                          const Index i,
+                                                                          tnlTuple< 2, Real >& normal ) const
+{
+   tnlAssert( ( dx == 0 || dx == 1 || dx == -1 ||
+                dy == 0 || dy == 1 || dy == -1 ) &&
+               dx * dy == 0, cerr << " dx = " << dx << " dy = " << dy << endl );
+   normal. x() = dx;
+   normal. y() = dy;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+bool tnlIdenticalGridGeometry< 2, Real, Device, Index > :: save( tnlFile& file ) const
+{
+   return true;
+};
+
+template< typename Real,
+          typename Device,
+          typename Index >
+bool tnlIdenticalGridGeometry< 2, Real, Device, Index > :: load( tnlFile& file )
+{
+   return true;
+};
+
+#endif /* TNLIDENTICALGRIDGEOMETRY_IMPL_H_ */
diff --git a/src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h b/src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h
index c708676805fa9c4f9734c1e9d670d3b678edf051..0df20164f55e44906cd7bd689a7a1c11820639e0 100644
--- a/src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h
+++ b/src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h
@@ -39,21 +39,20 @@ void tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index > > :: setFunction( Vec
 }
 
 template< typename Real, typename Device, typename Index >
-void tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index > > :: getGradient( const Index& i,
-                                                                             RealType& diffusion ) const
+Real tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index > > :: getDiffusion( const Index& c ) const
 {
    tnlAssert( this -> mesh, cerr << "No mesh was set in tnlLinearDiffusion. Use the bindMesh method." );
 
-   const Real hx = mesh -> getSpaceStep(). x() );
-   const Real hy = mesh -> getSpaceStep(). y() );
+   const Real hx = mesh -> getSpaceStep(). x();
+   const Real hy = mesh -> getSpaceStep(). y();
 
-   const Index e = mesh -> getElementNeighbour( i,  0,  1 );
-   const Index w = mesh -> getElementNeighbour( i,  0, -1 );
-   const Index n = mesh -> getElementNeighbour( i,  1,  0 );
-   const Index s = mesh -> getElementNeighbour( i, -1,  0 );
+   const Index e = mesh -> getElementNeighbour( c,  0,  1 );
+   const Index w = mesh -> getElementNeighbour( c,  0, -1 );
+   const Index n = mesh -> getElementNeighbour( c,  1,  0 );
+   const Index s = mesh -> getElementNeighbour( c, -1,  0 );
 
-   diffusion = ( f[ e ] - 2.0 * f[ c ] + f[ w ] ) / ( hx * hx ) +
-               ( f[ n ] - 2.0 * f[ c ] + f[ s ] ) / ( hy * hy );
+   return ( f[ e ] - 2.0 * f[ c ] + f[ w ] ) / ( hx * hx ) +
+          ( f[ n ] - 2.0 * f[ c ] + f[ s ] ) / ( hy * hy );
 }
 
 #endif
diff --git a/src/mesh/tnlGrid.h b/src/mesh/tnlGrid.h
index bea0d287d147198d8c0aff72f2799593f345f2f5..194fc6b2d14bf974fced1c2c1a3f1c93ea71e03a 100644
--- a/src/mesh/tnlGrid.h
+++ b/src/mesh/tnlGrid.h
@@ -109,6 +109,7 @@ class tnlGrid< 2, Real, Device, Index, Geometry > : public tnlObject
    typedef Real RealType;
    typedef Device DeviceType;
    typedef Index IndexType;
+   typedef Geometry GeometryType;
    enum { Dimensions = 2};
 
    tnlGrid();
@@ -131,11 +132,12 @@ class tnlGrid< 2, Real, Device, Index, Geometry > : public tnlObject
 
    const tnlTuple< 2, Real >& getProportions() const;
 
-   void setSpaceStep( const tnlTuple< 2, Real >& spaceStep );
+   void setParametricStep( const tnlTuple< 2, Real >& spaceStep );
 
-   tnlTuple< 2, Real > getSpaceStep() const;
+   tnlTuple< 2, Real > getParametricStep() const;
 
-   Index getElementIndex( const Index j, const Index i ) const;
+   Index getElementIndex( const Index j,
+                          const Index i ) const;
 
    Index getElementNeighbour( const Index Element,
                               const Index dy,
@@ -143,6 +145,23 @@ class tnlGrid< 2, Real, Device, Index, Geometry > : public tnlObject
 
    Index getDofs() const;
 
+   Real getElementMeasure( const Index j,
+                           const Index i ) const;
+
+   Real getElementsDistance( const Index j,
+                             const Index i,
+                             const Index dy,
+                             const Index dx ) const;
+
+   template< int dy, int dx >
+   Real getEdgeLength( const Index j,
+                       const Index i ) const;
+
+   template< int dy, int dx >
+   tnlTuple< 2, Real > getEdgeNormal( const Index j,
+                                      const Index i ) const;
+
+
    //! Method for saving the object to a file as a binary data
    bool save( tnlFile& file ) const;
 
@@ -155,8 +174,8 @@ class tnlGrid< 2, Real, Device, Index, Geometry > : public tnlObject
 
    template< typename MeshFunction >
    bool write( const MeshFunction& function,
-                const tnlString& fileName,
-                const tnlString& format ) const;
+               const tnlString& fileName,
+               const tnlString& format ) const;
 
    protected:
 
@@ -164,6 +183,8 @@ class tnlGrid< 2, Real, Device, Index, Geometry > : public tnlObject
 
    tnlTuple< 2, RealType > origin, proportions;
 
+   GeometryType geometry;
+
    IndexType dofs;
 
 };
diff --git a/src/mesh/tnlIdenticalGridGeometry.h b/src/mesh/tnlIdenticalGridGeometry.h
index 73b68533563106841ac5e1bac2ce64fb32b37ea3..98cba4fec3ffe1f2e3723c1282d8c1b42e72379e 100644
--- a/src/mesh/tnlIdenticalGridGeometry.h
+++ b/src/mesh/tnlIdenticalGridGeometry.h
@@ -32,7 +32,48 @@ class tnlIdenticalGridGeometry
    typedef Device DeviceType;
    typedef Index IndexType;
 
+   void setParametricStep( const tnlTuple< 2, Real >& parametricStep );
+
+   void getElementCoordinates( const Index j,
+                               const Index i,
+                               tnlTuple< 2, Real >& coordinates ) const;
+
+   Real getElementMeasure( const Index j,
+                           const Index i ) const;
+
+   template< Index dy, Index dx >
+   Real getElementsDistance( const Index j,
+                             const Index i ) const;
+
+   template< Index dy, Index dx >
+   void getEdgeCoordinates( const Index j,
+                            const Index i,
+                            tnlTuple< 2, Real >& coordinates ) const
+
+   template< Index dy, Index dx >
+   Real getEdgeLength( const Index j,
+                       const Index i ) const;
+
+   template< Index dy, Index dx >
+   tnlTuple< 2, Real > getEdgeNormal( const Index j,
+                                      const Index i ) const;
+
+   void getVertexCoordinates( const Index j,
+                              const Index i,
+                              tnlTuple< 2, Real >& coordinates ) const;
+
+   bool save( tnlFile& file ) const;
+
+   bool load( tnlFile& file );
+
+   protected:
+
+   tnlTuple< Dimensions, Real > parametricStep;
+
+   Real elementMeasure;
 
 };
 
+#include <implementation/mesh/tnlIdenticalGridGeometry_impl.h>
+
 #endif /* TNLIDENTICALGRIDGEOMETRY_H_ */
diff --git a/src/schemes/diffusion/tnlLinearDiffusion.h b/src/schemes/diffusion/tnlLinearDiffusion.h
index 92174fa08cf1e82c111f6a8d3697e3988fd3e005..7c767e06bbf18fb7638e4331ca96807a7edd99ee 100644
--- a/src/schemes/diffusion/tnlLinearDiffusion.h
+++ b/src/schemes/diffusion/tnlLinearDiffusion.h
@@ -42,8 +42,7 @@ class tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index > >
    template< typename Vector >
    void setFunction( Vector& f ); // TODO: add const
 
-   void getGradient( const Index& i,
-                     RealType& diffusion ) const;
+   RealType getDiffusion( const Index& i ) const;
    protected:
 
    // TODO: change to ConstSharedVector
@@ -52,6 +51,6 @@ class tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index > >
    const tnlGrid< 2, RealType, DeviceType, IndexType >* mesh;
 };
 
-#include <implementation/schemes/gradient/tnlLinearDiffusion_impl.h>
+#include <implementation/schemes/diffusion/tnlLinearDiffusion_impl.h>
 
 #endif