diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt
index 04f786f263dc79b3af2a8f43954b1bdd40e22806..5dd8d96f69df4167aaa72ce28c11ae078f3c3430 100755
--- a/src/core/CMakeLists.txt
+++ b/src/core/CMakeLists.txt
@@ -6,6 +6,7 @@ set (headers tnlArray.h
       	    tnlCuda.h
   		       tnlDataElement.h
   		       tnlDevice.h
+  		       tnlFeature.h
   		       tnlFile.h 
   		       tnlFlopsCounter.h
    		    tnlHost.h 
diff --git a/src/core/tnlFeature.h b/src/core/tnlFeature.h
new file mode 100644
index 0000000000000000000000000000000000000000..1368bb87bdb579a36450d4642d3de2992aba59d4
--- /dev/null
+++ b/src/core/tnlFeature.h
@@ -0,0 +1,30 @@
+/***************************************************************************
+                          tnlFeature.h  -  description
+                             -------------------
+    begin                : May 18, 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 TNLFEATURE_H_
+#define TNLFEATURE_H_
+
+template< bool featureEnabled >
+class tnlFeature
+{
+   public:
+
+   enum{ enabled = featureEnabled };
+};
+
+
+#endif /* TNLFEATURE_H_ */
diff --git a/src/implementation/mesh/tnlGrid2D_impl.h b/src/implementation/mesh/tnlGrid2D_impl.h
index b475d06f44c2969dbb1a1545615ece0efd8bf9af..760c113d7bdb03e6b4d7efc409333650d6589419 100644
--- a/src/implementation/mesh/tnlGrid2D_impl.h
+++ b/src/implementation/mesh/tnlGrid2D_impl.h
@@ -72,6 +72,12 @@ void tnlGrid< 2, Real, Device, Index, Geometry > :: setDimensions( const Index x
    parametricStep. x() = proportions. x() / xSize;
    parametricStep. y() = proportions. y() / ySize;
    geometry. setParametricStep( parametricStep );
+   if( GeometryType :: ElementMeasureStorage :: enabled )
+   {
+      elementsMeasure. setSize( this -> getDofs() );
+      dualElementsMeasure. setSize( this -> getDofs() );
+      edgeNormals. setSize( this -> getNumberOfEdges() );
+   }
 }
 
 template< typename Real,
@@ -154,18 +160,6 @@ template< typename Real,
 const typename tnlGrid< 2, Real, Device, Index, Geometry > :: VertexType& 
    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 "
-                   << this -> getName() );
-   tnlAssert( dimensions. y() > 0,
-              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 > parametricStep;
-   parametricStep. x() =
-            this -> proportions. x() / ( Real ) ( this -> dimensions. x() - 1 );
-   parametricStep. y() =
-            this -> proportions. y() / ( Real ) ( this -> dimensions. y() - 1 );*/
    return geometry. getParametricStep();
 }
 
@@ -178,15 +172,80 @@ Index tnlGrid< 2, Real, Device, Index, Geometry > :: getElementIndex( const Inde
    tnlAssert( i < dimensions. x(),
               cerr << "Index i ( " << i
                    << " ) is out of range ( " << dimensions. x()
-                   << " ) in tnlGrid " << this -> getName(); )
+                   << " ) in tnlGrid " << this -> getName(); );
    tnlAssert( j < dimensions. y(),
               cerr << "Index j ( " << j
                    << " ) is out of range ( " << dimensions. y()
-                   << " ) in tnlGrid " << this -> getName(); )
+                   << " ) in tnlGrid " << this -> getName(); );
 
    return j * this -> dimensions. x() + i;
 }
 
+template< typename Real,
+          typename Device,
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+Index tnlGrid< 2, Real, Device, Index, Geometry > :: getEdgeIndex( const Index i,
+                                                                   const Index j,
+                                                                   const Index dx,
+                                                                   const Index dy ) const
+{
+   tnlAssert( i < dimensions. x(),
+              cerr << "Index i ( " << i
+                   << " ) is out of range ( " << dimensions. x()
+                   << " ) in tnlGrid " << this -> getName(); );
+   tnlAssert( j < dimensions. y(),
+              cerr << "Index j ( " << j
+                   << " ) is out of range ( " << dimensions. y()
+                   << " ) in tnlGrid " << this -> getName(); );
+   tnlAssert( dx == 0 && ( dy == 1 || dy == -1 ) ||
+              dy == 0 && ( dx == 1 || dx == -1 ),
+              cerr << "dx = " << dx << ", dy = " << dy << endl;);
+   return ( j + dy ) * this -> dimensions. x() + i + dx;
+}
+
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          template< int, typename, typename, typename > class Geometry >
+void tnlGrid< 2, Real, Device, Index, Geometry > :: refresh()
+{
+   if( GeometryType :: ElementMeasureStorage :: enabled )
+      for( IndexType j = 0; j < dimensions. y(); j ++ )
+         for( IndexType i = 0; i < dimensions. x(); i ++ )
+            elementsMeasure[ getElementIndex( i, j ) ] =
+                     geometry. getElementMeasure( CoordinatesType( i, j ) );
+
+   if( GeometryType :: DualElementMeasureStorage :: enabled )
+      for( IndexType j = 1; j < dimensions. y() - 1; j ++ )
+         for( IndexType i = 1; i < dimensions. x() - 1; i ++ )
+         {
+            dualElementsMeasure[ getElementIndex( i + 1, j ) ] =
+                     geometry. getDualElementMeasure<  1,  0 >( CoordinatesType( i, j ) );
+            dualElementsMeasure[ getElementIndex( i - 1, j ) ] =
+                     geometry. getDualElementMeasure< -1,  0 >( CoordinatesType( i, j ) );
+            dualElementsMeasure[ getElementIndex( i, j + 1 ) ] =
+                     geometry. getDualElementMeasure<  0,  1 >( CoordinatesType( i, j ) );
+            dualElementsMeasure[ getElementIndex( i, j - 1 ) ] =
+                     geometry. getDualElementMeasure<  0, -1 >( CoordinatesType( i, j ) );
+         }
+
+   if( GeometryType :: EdgeNormalStorage :: enabled )
+      for( IndexType j = 0; j < dimensions. y(); j ++ )
+         for( IndexType i = 0; i < dimensions. x(); i ++ )
+         {
+            geometry. getEdgeNormal<  1,  0 >( CoordinatesType( i, j ),
+                                               edgeNormals[ getEdgeIndex( i, j, 1, 0 ) ] );
+            geometry. getEdgeNormal< -1,  0 >( CoordinatesType( i, j ),
+                                                edgeNormals[ getEdgeIndex( i, j, -1, 0 ) ] );
+            geometry. getEdgeNormal<  0,  1 >( CoordinatesType( i, j ),
+                                               edgeNormals[ getEdgeIndex( i, j, 0, 1 ) ]  );
+            geometry. getEdgeNormal<  0, -1 >( CoordinatesType( i, j ),
+                                               edgeNormals[ getEdgeIndex( i, j, 0, -1 ) ] );
+         }
+}
+
 template< typename Real,
           typename Device,
           typename Index,
@@ -218,7 +277,6 @@ Index tnlGrid< 2, Real, Device, Index, Geometry > :: getElementNeighbour( const
    return Element + dy * this -> dimensions. x() + dx;
 }
 
-
 template< typename Real,
           typename Device,
           typename Index,
@@ -232,53 +290,44 @@ template< typename Real,
           typename Device,
           typename Index,
           template< int, typename, typename, typename > class Geometry >
-void tnlGrid< 2, Real, Device, Index, Geometry > :: getElementCenter( const CoordinatesType& coordinates,
-                                                                      VertexType& center ) const
+Index tnlGrid< 2, Real, Device, Index, Geometry > :: getNumberOfEdges() const
 {
-      geometry. getElementCenter( origin, coordinates, center );
+   return ( this -> dimensions. x() + 1 ) * ( this -> dimensions. y() + 1 );
 }
 
 template< typename Real,
           typename Device,
           typename Index,
           template< int, typename, typename, typename > class Geometry >
-Real tnlGrid< 2, Real, Device, Index, Geometry > :: getElementMeasure( const CoordinatesType& coordinates ) const
+void tnlGrid< 2, Real, Device, Index, Geometry > :: getElementCenter( const CoordinatesType& coordinates,
+                                                                      VertexType& center ) const
 {
-   return geometry. getElementMeasure( coordinates );
+      geometry. getElementCenter( origin, coordinates, center );
 }
 
 template< typename Real,
           typename Device,
           typename Index,
           template< int, typename, typename, typename > class Geometry >
-   template< int dx, int dy >
-Real tnlGrid< 2, Real, Device, Index, Geometry > :: getElementCoVolumeMeasure( const CoordinatesType& coordinates ) const
+Real tnlGrid< 2, Real, Device, Index, Geometry > :: getElementMeasure( const CoordinatesType& coordinates ) const
 {
-   return geometry. getElementCoVolumeMeasure< dx, dy >( coordinates );
+   if( GeometryType :: ElementMeasureStorage :: enabled )
+      return elementsMeasure[ getElementIndex( coordinates. x(), coordinates. y() ) ];
+   return geometry. getElementMeasure( coordinates );
 }
 
 template< typename Real,
           typename Device,
           typename Index,
           template< int, typename, typename, typename > class Geometry >
-Real tnlGrid< 2, Real, Device, Index, Geometry > :: getElementsDistance( const CoordinatesType& c1,
-                                                                         const CoordinatesType& c2 ) const
+   template< int dx, int dy >
+Real tnlGrid< 2, Real, Device, Index, Geometry > :: getDualElementMeasure( const CoordinatesType& coordinates ) const
 {
-   return geometry. getElementsDistance( c1, c2 );
+   if( GeometryType :: DualElementMeasureStorage :: enabled )
+      return dualElementsMeasure[ getElementIndex( coordinates. x() + dx, coordinates. y() + dy ) ];
+   return geometry. getDualElementMeasure< dx, dy >( coordinates );
 }
 
-/*
-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,
@@ -287,7 +336,10 @@ template< int dx, int dy >
 void tnlGrid< 2, Real, Device, Index, Geometry > :: getEdgeNormal( const CoordinatesType& coordinates,
                                                                    VertexType& normal ) const
 {
-   return geometry. getEdgeNormal< dx, dy >( coordinates, normal );
+   //if( GeometryType :: EdgeNormalStorage :: enabled )
+   //   normal = edgeNormals[ getEdgeIndex( coordinates. x(), coordinates. y(), dx, dy ) ];
+   //else
+      return geometry. getEdgeNormal< dx, dy >( coordinates, normal );
 }
 
 template< typename Real,
diff --git a/src/implementation/mesh/tnlIdenticalGridGeometry_impl.h b/src/implementation/mesh/tnlIdenticalGridGeometry_impl.h
index fdc91b8605113d4c46e634c510e2f918e7c16041..aae0ff4b3e2f7f7e3ddf6af1fb82b757a84af9d5 100644
--- a/src/implementation/mesh/tnlIdenticalGridGeometry_impl.h
+++ b/src/implementation/mesh/tnlIdenticalGridGeometry_impl.h
@@ -167,55 +167,11 @@ template< typename Real,
           typename Device,
           typename Index >
    template< int dx, int dy >
-Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementCoVolumeMeasure( const CoordinatesType& coordinates ) const
+Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getDualElementMeasure( const CoordinatesType& coordinates ) const
 {
    return 0.5 * elementMeasure;
 }
 
-
-template< typename Real,
-          typename Device,
-          typename Index >
-Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementsDistance( const CoordinatesType& c1,
-                                                                                const CoordinatesType& c2 ) const
-{
-   CoordinatesType c( c2 );
-   c -= c1;
-   if( c. y() == 0 )
-      return parametricStep. x() * fabs( c. x() );
-   if( c. x() == 0 )
-      return parametricStep. y() * fabs( c. y() );
-   return sqrt( c. x() * c. x() + c. y() * c. y() );
-}
-
-/*
-template< typename Real,
-          typename Device,
-          typename Index >
-template< Index dx, Index dy >
-void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getEdgeCoordinates( const Index i,
-                                                                               const Index j,
-                                                                               const VertexType& origin,
-                                                                               VertexType& 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< Index dx, Index dy >
-Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getEdgeLength( const Index i,
-                                                                          const Index j ) 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 >
diff --git a/src/implementation/mesh/tnlLinearGridGeometry_impl.h b/src/implementation/mesh/tnlLinearGridGeometry_impl.h
index 498b3279bbef1030330187e1566d9863a7129a91..993991c6aa22519e4fcf0c43f1587b7c0de0c714 100644
--- a/src/implementation/mesh/tnlLinearGridGeometry_impl.h
+++ b/src/implementation/mesh/tnlLinearGridGeometry_impl.h
@@ -175,59 +175,33 @@ template< typename Real,
           typename Device,
           typename Index >
    template< int dx, int dy >
-Real tnlLinearGridGeometry< 2, Real, Device, Index > :: getElementCoVolumeMeasure( const CoordinatesType& coordinates ) const
+Real tnlLinearGridGeometry< 2, Real, Device, Index > :: getDualElementMeasure( const CoordinatesType& coordinates ) const
 {
    tnlAssert( ( dx == 0 && ( dy == 1 || dy == -1 ) ) ||
               ( dy == 0 && ( dx == 1 || dx == -1 ) ),
               cerr << " dx = " << dx << " dy = " << dy << endl );
-
-   return 0.5 * elementMeasure;
-}
-
-
-template< typename Real,
-          typename Device,
-          typename Index >
-Real tnlLinearGridGeometry< 2, Real, Device, Index > :: getElementsDistance( const CoordinatesType& c1,
-                                                                                const CoordinatesType& c2 ) const
-{
-   CoordinatesType c( c2 );
-   c -= c1;
-   if( c. y() == 0 )
-      return parametricStep. x() * fabs( c. x() );
-   if( c. x() == 0 )
-      return parametricStep. y() * fabs( c. y() );
-   return sqrt( c. x() * c. x() + c. y() * c. y() );
-}
-
-/*
-template< typename Real,
-          typename Device,
-          typename Index >
-template< Index dx, Index dy >
-void tnlLinearGridGeometry< 2, Real, Device, Index > :: getEdgeCoordinates( const Index i,
-                                                                               const Index j,
-                                                                               const VertexType& origin,
-                                                                               VertexType& 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();
+   VertexType v0, v1, v2, v3;
+   const VertexType origin( 0.0 );
+   this -> getElementCenter( origin, coordinates, v0 );
+   if( dy == 0 )
+   {
+      this -> template getVertex< dx, -1 >( coordinates, origin, v1 );
+      this -> template getVertex< dx, 1 >( coordinates, origin, v2 );
+      CoordinatesType c2( coordinates );
+      c2. x() += dx;
+      this -> getElementCenter( origin, c2, v3 );
+   }
+   if( dx == 0 )
+   {
+      this -> template getVertex< -1, dy >( coordinates, origin, v1 );
+      this -> template getVertex<  1, dy >( coordinates, origin, v2 );
+      CoordinatesType c2( coordinates );
+      c2. y() += dy;
+      this -> getElementCenter( origin, c2, v3 );
+   }
+   return tnlTriangleArea( v0, v1, v3 ) + tnlTriangleArea( v2, v3, v1 );
 }
 
-template< typename Real,
-          typename Device,
-          typename Index >
-template< Index dx, Index dy >
-Real tnlLinearGridGeometry< 2, Real, Device, Index > :: getEdgeLength( const Index i,
-                                                                          const Index j ) 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 >
diff --git a/src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h b/src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h
index 06bff8b96d4f26c92fc356e4c418e8d177485b30..0cb9beb46c859e5fa6cf0a4248a21d829df8323e 100644
--- a/src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h
+++ b/src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h
@@ -121,45 +121,45 @@ Real tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index, GridGeometry > > :: ge
    this -> mesh -> template getVertex< -1,  1 >( cCoordinates, wnVertex );
    this -> mesh -> template getVertex< -1, -1 >( cCoordinates, wsVertex );
 
-   const RealType f_x_e = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< 1, 0 >( cCoordinates ) *
+   const RealType f_x_e = 1.0 / this -> mesh -> template getDualElementMeasure< 1, 0 >( cCoordinates ) *
                           ( f_cen * ( cCenter. y() - enVertex. y() ) +
                             f_ces * ( esVertex. y() - cCenter. y() ) +
                             0.5 * ( f_es + f[ e ] ) * ( eCenter. y() - esVertex. y() ) +
                             0.5 * ( f_en + f[ e ] ) * ( enVertex. y() - eCenter. y() ) );
-   const RealType f_y_e = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< 1, 0 >( cCoordinates ) *
+   const RealType f_y_e = 1.0 / this -> mesh -> template getDualElementMeasure< 1, 0 >( cCoordinates ) *
                           ( f_cen * ( enVertex. x() - cCenter. x() ) +
                             f_ces * ( cCenter. x() - esVertex. x() ) +
                             0.5 * ( f_es + f[ e ] ) * ( esVertex. x() - eCenter. x() ) +
                             0.5 * ( f_en + f[ e ] ) * ( eCenter. x() ) - enVertex. x() );
 
-   const RealType f_x_w = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< -1, 0 >( cCoordinates ) *
+   const RealType f_x_w = 1.0 / this -> mesh -> template getDualElementMeasure< -1, 0 >( cCoordinates ) *
                           ( f_cwn * ( wnVertex. y() - cCenter. y() ) +
                             f_cws * ( cCenter. y() - wsVertex. y()  ) +
                             0.5 * ( f_ws + f[ w ] ) * ( wsVertex. y() - wCenter. y() ) +
                             0.5 * ( f_wn + f[ w ] ) * ( wCenter. y() - wnVertex. y() ) );
-   const RealType f_y_w = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< -1, 0 >( cCoordinates ) *
+   const RealType f_y_w = 1.0 / this -> mesh -> template getDualElementMeasure< -1, 0 >( cCoordinates ) *
                           ( f_cwn * ( cCenter. x() - wnVertex. x() ) +
                             f_cws * ( wsVertex. x() - cCenter. x() ) +
                             0.5 * ( f_ws + f[ w ] ) * ( wCenter. x() - wsVertex. x() ) +
                             0.5 * ( f_wn + f[ w ] ) * ( wnVertex. x() - wCenter. x() ) );
 
-   const RealType f_x_n = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< 0, 1 >( cCoordinates ) *
+   const RealType f_x_n = 1.0 / this -> mesh -> template getDualElementMeasure< 0, 1 >( cCoordinates ) *
                           ( f_cen * ( enVertex. y() - cCenter. y() ) +
                             f_cwn * ( cCenter. y() - wnVertex. y() ) +
                             0.5 * ( f_en + f[ n ] ) * ( nCenter. y() - enVertex. y() ) +
                             0.5 * ( f_wn + f[ n ] ) * ( wnVertex. y() - nCenter. y() ) );
-   const RealType f_y_n = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< 0, 1 >( cCoordinates ) *
+   const RealType f_y_n = 1.0 / this -> mesh -> template getDualElementMeasure< 0, 1 >( cCoordinates ) *
                           ( f_cen * ( cCenter. x() - enVertex. x() ) +
                             f_cwn * ( wnVertex. x() - cCenter. x() ) +
                             0.5 * ( f_en + f[ n ] ) * ( enVertex. x() - nCenter. x() ) +
                             0.5 * ( f_wn + f[ n ] ) * ( nCenter. x() - wnVertex. x() ) );
 
-   const RealType f_x_s = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< 0, -1 >( cCoordinates ) *
+   const RealType f_x_s = 1.0 / this -> mesh -> template getDualElementMeasure< 0, -1 >( cCoordinates ) *
                           ( f_ces * ( cCenter. y() - esVertex. y() ) +
                             f_cws * ( wsVertex. y() - cCenter. y() ) +
                             0.5 * ( f_es + f[ s ] ) * ( esVertex. y() - sCenter. y() ) +
                             0.5 * ( f_ws + f[ s ] ) * ( sCenter. y() - wsVertex. y() ) );
-   const RealType f_y_s = 1.0 / this -> mesh -> template getElementCoVolumeMeasure< 0, -1 >( cCoordinates ) *
+   const RealType f_y_s = 1.0 / this -> mesh -> template getDualElementMeasure< 0, -1 >( cCoordinates ) *
                           ( f_ces * ( esVertex. x() - cCenter. x() ) +
                             f_cws * ( cCenter. x() - wsVertex. x() ) +
                             0.5 * ( f_es + f[ s ] ) * ( sCenter. x() - esVertex. x() ) +
diff --git a/src/mesh/tnlGrid.h b/src/mesh/tnlGrid.h
index e0b7d2404cf13ba530ea4570fc76f37002b47c59..b2eabf16f08c1a2745eb97b2c6f5f9b15b1dc173 100644
--- a/src/mesh/tnlGrid.h
+++ b/src/mesh/tnlGrid.h
@@ -21,6 +21,7 @@
 #include <core/tnlObject.h>
 #include <core/tnlHost.h>
 #include <core/tnlTuple.h>
+#include <core/tnlVector.h>
 #include <mesh/tnlIdenticalGridGeometry.h>
 
 template< int Dimensions,
@@ -147,6 +148,13 @@ class tnlGrid< 2, Real, Device, Index, Geometry > : public tnlObject
    Index getElementIndex( const Index i,
                           const Index j ) const;
 
+   Index getEdgeIndex( const Index i,
+                       const Index j,
+                       const Index dx,
+                       const Index dy ) const;
+
+   void refresh();
+
    void getElementCoordinates( const Index i,
                                CoordinatesType& coordinates ) const;
 
@@ -156,13 +164,15 @@ class tnlGrid< 2, Real, Device, Index, Geometry > : public tnlObject
 
    Index getDofs() const;
 
+   Index getNumberOfEdges() const;
+
    void getElementCenter( const CoordinatesType& coordinates,
                           VertexType& center ) const;
 
    Real getElementMeasure( const CoordinatesType& coordinates ) const;
 
    template< int dx, int dy >
-   Real getElementCoVolumeMeasure( const CoordinatesType& coordinates ) const;
+   Real getDualElementMeasure( const CoordinatesType& coordinates ) const;
 
    template< int dx, int dy >
    void getEdgeNormal( const CoordinatesType& elementCoordinates,
@@ -172,13 +182,6 @@ class tnlGrid< 2, Real, Device, Index, Geometry > : public tnlObject
    void getVertex( const CoordinatesType& elementCoordinates,
                    VertexType& vertex ) const;
 
-   Real getElementsDistance( const CoordinatesType& c1,
-                             const CoordinatesType& c2 ) const;
-
-   /*template< int dy, int dx >
-   Real getEdgeLength( const Index j,
-                       const Index i ) const;*/
-
    //! Method for saving the object to a file as a binary data
    bool save( tnlFile& file ) const;
 
@@ -196,14 +199,17 @@ class tnlGrid< 2, Real, Device, Index, Geometry > : public tnlObject
 
    protected:
 
-   tnlTuple< 2, IndexType > dimensions;
+   CoordinatesType dimensions;
 
-   tnlTuple< 2, RealType > origin, proportions;
+   VertexType origin, proportions;
 
    GeometryType geometry;
 
    IndexType dofs;
 
+   tnlVector< Real, Device, Index > elementsMeasure, dualElementsMeasure;
+   tnlVector< VertexType, Device, Index > edgeNormals;
+
 };
 
 template< typename Real,
diff --git a/src/mesh/tnlIdenticalGridGeometry.h b/src/mesh/tnlIdenticalGridGeometry.h
index 13012cfdaec71caa4c2a4f09983730b8546d89ad..41464693cf8a8e03345a478c485c0d874a74659c 100644
--- a/src/mesh/tnlIdenticalGridGeometry.h
+++ b/src/mesh/tnlIdenticalGridGeometry.h
@@ -19,6 +19,7 @@
 #define TNLIDENTICALGRIDGEOMETRY_H_
 
 #include <core/tnlHost.h>
+#include <core/tnlFeature.h>
 
 template< int Dimensions,
           typename Real = double,
@@ -93,6 +94,9 @@ class tnlIdenticalGridGeometry< 2, Real, Device, Index >
    typedef Index IndexType;
    typedef tnlTuple< 2, Index > CoordinatesType;
    typedef tnlTuple< 2, Real > VertexType;
+   typedef tnlFeature< false > ElementMeasureStorage;
+   typedef tnlFeature< false > DualElementMeasureStorage;
+   typedef tnlFeature< false > EdgeNormalStorage;
 
    void setParametricStep( const VertexType& parametricStep );
 
@@ -105,20 +109,7 @@ class tnlIdenticalGridGeometry< 2, Real, Device, Index >
    Real getElementMeasure( const tnlTuple< 2, Index >& coordinates ) const;
 
    template< int dx, int dy >
-   Real getElementCoVolumeMeasure( const CoordinatesType& coordinates ) const;
-
-   Real getElementsDistance( const CoordinatesType& c1,
-                             const CoordinatesType& c2 ) const;
-
-   /*template< Index dx, Index dy >
-   void getEdgeCoordinates( const Index i,
-                            const Index j,
-                            const VertexType& origin,
-                            VertexType& coordinates ) const;
-
-   template< Index dx, Index dy >
-   Real getEdgeLength( const Index i,
-                       const Index j ) const;*/
+   Real getDualElementMeasure( const CoordinatesType& coordinates ) const;
 
    template< Index dx, Index dy >
    void getEdgeNormal( const CoordinatesType& coordinates,
diff --git a/src/mesh/tnlLinearGridGeometry.h b/src/mesh/tnlLinearGridGeometry.h
index 7c8263b2dc34da1449f2f3e9fbe0bb493656ef6d..c105aeb4e49de9eef665e61fa789576a49ee1d23 100644
--- a/src/mesh/tnlLinearGridGeometry.h
+++ b/src/mesh/tnlLinearGridGeometry.h
@@ -94,6 +94,10 @@ class tnlLinearGridGeometry< 2, Real, Device, Index >
    typedef Index IndexType;
    typedef tnlTuple< 2, Index > CoordinatesType;
    typedef tnlTuple< 2, Real > VertexType;
+   typedef tnlFeature< true > ElementMeasureStorage;
+   typedef tnlFeature< true > DualElementMeasureStorage;
+   typedef tnlFeature< true > EdgeNormalStorage;
+
 
    void setParametricStep( const VertexType& parametricStep );
 
@@ -106,20 +110,7 @@ class tnlLinearGridGeometry< 2, Real, Device, Index >
    Real getElementMeasure( const tnlTuple< 2, Index >& coordinates ) const;
 
    template< int dx, int dy >
-   Real getElementCoVolumeMeasure( const CoordinatesType& coordinates ) const;
-
-   Real getElementsDistance( const CoordinatesType& c1,
-                             const CoordinatesType& c2 ) const;
-
-   /*template< Index dx, Index dy >
-   void getEdgeCoordinates( const Index i,
-                            const Index j,
-                            const VertexType& origin,
-                            VertexType& coordinates ) const;
-
-   template< Index dx, Index dy >
-   Real getEdgeLength( const Index i,
-                       const Index j ) const;*/
+   Real getDualElementMeasure( const CoordinatesType& coordinates ) const;
 
    template< Index dx, Index dy >
    void getEdgeNormal( const CoordinatesType& coordinates,