diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index ab3efc8d7485cd0c8c7ada67626fafd0ffd48b04..6063a109864aaa2e0aec4e9e4b017def62c202a9 100755
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -2,3 +2,4 @@ add_subdirectory( make-project )
 add_subdirectory( simple-solver )
 add_subdirectory( heat-equation )
 add_subdirectory( navier-stokes )
+add_subdirectory( mean-curvature-flow )
diff --git a/examples/mean-curvature-flow/CMakeLists.txt b/examples/mean-curvature-flow/CMakeLists.txt
new file mode 100755
index 0000000000000000000000000000000000000000..9497b2f7a059be5008ffc4600a8d88ab7af6e489
--- /dev/null
+++ b/examples/mean-curvature-flow/CMakeLists.txt
@@ -0,0 +1,23 @@
+INSTALL( FILES tnlBackwardFiniteDifference.h
+					tnlBackwardFiniteDifference_impl.h
+					tnlExactFlowDiffusion.h
+					tnlExactFlowDiffusion_impl.h
+					tnlForwardFiniteDifference.h
+					tnlForwardFiniteDifference_impl.h
+					tnl-mean-curvature-flow.cpp
+					tnl-mean-curvature-flow.h
+					tnlMeanCurvatureFlowDiffusion.h
+					tnlMeanCurvatureFlowDiffusion_impl.h
+					tnl-mean-curvature-flow-eoc.cpp
+					tnl-mean-curvature-flow-eoc.h
+					tnlMeanCurvatureFlowEocProblem.h
+					tnlMeanCurvatureFlowEocProblem_impl.h
+					tnlMeanCurvatureFlowEocRhs.h
+					tnlMeanCurvatureFlowEocRhs_impl.h
+					tnlMeanCurvatureFlowProblem.h
+					tnlMeanCurvatureFlowProblem_impl.h
+					tnlOperatorQ.h
+					tnlOperatorQ_impl.h
+					tnl-run-mean-curvature-flow
+					tnl-run-mean-curvature-flow-eoc-test
+         DESTINATION share/tnl-${tnlVersion}/examples/mean-curvature-flow )
diff --git a/examples/mean-curvature-flow/Makefile b/examples/mean-curvature-flow/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..9ca9bc675cf55868456640d8162d6f572cb7230a
--- /dev/null
+++ b/examples/mean-curvature-flow/Makefile
@@ -0,0 +1,36 @@
+TNL_VERSION=0.1
+TNL_INSTALL_DIR=${HOME}/local/lib
+TNL_INCLUDE_DIR=${HOME}/local/include/tnl-${TNL_VERSION}
+
+TARGET = mean-curvature-flow-eoc
+INSTALL_DIR = ${HOME}/local
+CXX = g++
+CUDA_CXX = nvcc
+CXX_FLAGS = -std=gnu++0x -I$(TNL_INCLUDE_DIR)  -g -O0
+LD_FLAGS = -L$(TNL_INSTALL_DIR) -ltnl-dbg-0.1
+
+SOURCES = tnl-mean-curvature-flow-eoc.cpp
+HEADERS = 
+OBJECTS = tnl-mean-curvature-flow-eoc.o
+DIST = $(SOURCES) Makefile
+
+all: $(TARGET)
+clean: 
+	rm -f $(OBJECTS)	
+
+dist: $(DIST)
+	tar zcvf $(TARGET).tgz $(DIST) 
+
+install: $(TARGET)
+	cp $(TARGET) $(INSTALL_DIR)/bin
+	cp $(INSTALL_DIR)/share/tnl-0.1/examples/mean-curvature-flow-eoc
+
+uninstall: $(TARGET)
+	rm -f $(INSTALL_DIR)/bin/$(TARGET) 
+	rm -f $(INSTALL_DIR)/share/tnl-0.1/examples/mean-curvature-flow-eoc
+
+$(TARGET): $(OBJECTS)
+	$(CXX) -o $(TARGET) $(OBJECTS) $(LD_FLAGS)
+
+%.o: %.cpp $(HEADERS)
+	$(CXX) -c -o $@ $(CXX_FLAGS) $<
diff --git a/examples/mean-curvature-flow/tnl-mean-curvature-flow-eoc.cpp b/examples/mean-curvature-flow/tnl-mean-curvature-flow-eoc.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..256e235dc7fb854dcacbbba29ee8e99e8dadb721
--- /dev/null
+++ b/examples/mean-curvature-flow/tnl-mean-curvature-flow-eoc.cpp
@@ -0,0 +1,20 @@
+/***************************************************************************
+                          tnl-heat-equation-eoc.cpp  -  description
+                             -------------------
+    begin                : Sep 7, 2014
+    copyright            : (C) 2014 by 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.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#include "tnl-mean-curvature-flow-eoc.h"
+
+
diff --git a/examples/mean-curvature-flow/tnl-mean-curvature-flow-eoc.h b/examples/mean-curvature-flow/tnl-mean-curvature-flow-eoc.h
new file mode 100644
index 0000000000000000000000000000000000000000..0b09c6059b15ad9592dfa332093dbee17107593f
--- /dev/null
+++ b/examples/mean-curvature-flow/tnl-mean-curvature-flow-eoc.h
@@ -0,0 +1,86 @@
+/***************************************************************************
+                          tnl-heat-equation-eoc.h  -  description
+                             -------------------
+    begin                : Nov 29, 2014
+    copyright            : (C) 2014 by 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 TNL_MEAN_CURVATURE_FLOW_EOC_H_
+#define TNL_MEAN_CURVATURE_FLOW_EOC_H_
+
+#include <solvers/tnlSolver.h>
+#include <solvers/tnlFastBuildConfig.h>
+#include <solvers/tnlConfigTags.h>
+#include <functions/tnlTestFunction.h>
+#include <operators/tnlAnalyticDirichletBoundaryConditions.h>
+#include <operators/tnlAnalyticNeumannBoundaryConditions.h>
+#include "tnlMeanCurvatureFlowEocRhs.h"
+#include "tnlMeanCurvatureFlowEocProblem.h"
+#include "tnlExactFlowDiffusion.h"
+#include "tnlMeanCurvatureFlowDiffusion.h"
+
+//typedef tnlDefaultConfigTag BuildConfig;
+typedef tnlFastBuildConfig BuildConfig;
+
+template< typename ConfigTag >
+class meanCurvatureFlowEocConfig
+{
+   public:
+      static void configSetup( tnlConfigDescription& config )
+      {
+         config.addDelimiter( "Mean Curvature Flow EOC settings:" );
+         config.addDelimiter( "Tests setting::" );
+         tnlTestFunction< 2, double >::configSetup( config );
+      }
+};
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MeshType,
+          typename ConfigTag,
+          typename SolverStarter >
+class meanCurvatureFlowEocSetter
+{
+   public:
+
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+
+   typedef tnlStaticVector< MeshType::Dimensions, Real > Vertex;
+
+   static bool run( const tnlParameterContainer& parameters )
+   {
+      enum { Dimensions = MeshType::Dimensions };
+      typedef tnlMeanCurvatureFlowDiffusion< MeshType, Real, Index > ApproximateOperator;
+      typedef tnlExactFlowDiffusion< Dimensions > ExactOperator;
+      typedef tnlTestFunction< MeshType::Dimensions, Real, Device > TestFunction;
+      typedef tnlMeanCurvatureFlowEocRhs< ExactOperator, TestFunction > RightHandSide;
+      typedef tnlStaticVector < MeshType::Dimensions, Real > Vertex;
+      typedef tnlAnalyticNeumannBoundaryConditions< MeshType, TestFunction, Real, Index > BoundaryConditions;
+      typedef tnlMeanCurvatureFlowEocProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Solver;
+      SolverStarter solverStarter;
+      return solverStarter.template run< Solver >( parameters );
+   };
+};
+
+int main( int argc, char* argv[] )
+{
+   tnlSolver< meanCurvatureFlowEocSetter, meanCurvatureFlowEocConfig, BuildConfig > solver;
+   if( ! solver. run( argc, argv ) )
+      return EXIT_FAILURE;
+   return EXIT_SUCCESS;
+}
+
+#endif /* TNL_MEAN_CURVATURE_FLOW_EOC_H_ */
diff --git a/examples/mean-curvature-flow/tnl-mean-curvature-flow.cpp b/examples/mean-curvature-flow/tnl-mean-curvature-flow.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..78dc54d4cdad74eb75b7881fce6f274d9b7c6874
--- /dev/null
+++ b/examples/mean-curvature-flow/tnl-mean-curvature-flow.cpp
@@ -0,0 +1,18 @@
+/***************************************************************************
+                          tnl-heat-equation.cpp  -  description
+                             -------------------
+    begin                : Jan 12, 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.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#include "tnl-mean-curvature-flow.h"
diff --git a/examples/mean-curvature-flow/tnl-mean-curvature-flow.h b/examples/mean-curvature-flow/tnl-mean-curvature-flow.h
new file mode 100644
index 0000000000000000000000000000000000000000..a363b26c469e0bb87014c6f56ba2827acd628958
--- /dev/null
+++ b/examples/mean-curvature-flow/tnl-mean-curvature-flow.h
@@ -0,0 +1,117 @@
+/***************************************************************************
+                          tnl-heat-equation.h  -  description
+                             -------------------
+    begin                : Nov 29, 2014
+    copyright            : (C) 2014 by 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 TNL_MEAN_CURVATIVE_FLOW_H_
+#define TNL_MEAN_CURVATIVE_FLOW_H_
+
+#include <solvers/tnlSolver.h>
+#include <solvers/tnlFastBuildConfig.h>
+#include <solvers/tnlConfigTags.h>
+#include <operators/diffusion/tnlLinearDiffusion.h>
+#include <operators/tnlAnalyticDirichletBoundaryConditions.h>
+#include <operators/tnlDirichletBoundaryConditions.h>
+#include <operators/tnlAnalyticNeumannBoundaryConditions.h>
+#include <operators/tnlNeumannBoundaryConditions.h>
+#include <functions/tnlConstantFunction.h>
+#include <problems/tnlHeatEquationProblem.h>
+
+#include "tnlMeanCurvatureFlowProblem.h"
+#include "tnlMeanCurvatureFlowDiffusion.h"
+
+//typedef tnlDefaultConfigTag BuildConfig;
+typedef tnlFastBuildConfig BuildConfig;
+
+template< typename ConfigTag >
+class meanCurvatureFlowConfig
+{
+   public:
+      static void configSetup( tnlConfigDescription& config )
+      {
+         config.addDelimiter( "Mean Curvature Flow settings:" );
+         config.addEntry< tnlString >( "boundary-conditions-type", "Choose the boundary conditions type.", "dirichlet");
+            config.addEntryEnum< tnlString >( "dirichlet" );
+            config.addEntryEnum< tnlString >( "neumann" );
+
+         config.addEntry< tnlString >( "boundary-conditions-file", "File with the values of the boundary conditions.", "boundary.tnl" );
+         config.addEntry< double >( "boundary-conditions-constant", "This sets a value in case of the constant boundary conditions." );
+         config.addEntry< tnlString >( "initial-condition", "File with the initial condition.", "initial.tnl");
+      };
+};
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MeshType,
+          typename ConfigTag,
+          typename SolverStarter >
+class meanCurvatureFlowSetter
+{
+   public:
+
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+
+   typedef tnlStaticVector< MeshType::Dimensions, Real > Vertex;
+
+   static bool run( const tnlParameterContainer& parameters )
+   {
+      enum { Dimensions = MeshType::Dimensions };
+      typedef tnlMeanCurvatureFlowDiffusion< MeshType, Real, Index > ApproximateOperator;
+      typedef tnlConstantFunction< Dimensions, Real > RightHandSide;
+      typedef tnlStaticVector < MeshType::Dimensions, Real > Vertex;
+
+      tnlString boundaryConditionsType = parameters.getParameter< tnlString >( "boundary-conditions-type" );
+      if( parameters.checkParameter( "boundary-conditions-constant" ) )
+      {
+         typedef tnlConstantFunction< Dimensions, Real > ConstantFunction;
+         if( boundaryConditionsType == "dirichlet" )
+         {
+            typedef tnlAnalyticDirichletBoundaryConditions< MeshType, ConstantFunction, Real, Index > BoundaryConditions;
+            typedef tnlMeanCurvatureFlowProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Solver;
+            SolverStarter solverStarter;
+            return solverStarter.template run< Solver >( parameters );
+         }
+         typedef tnlAnalyticNeumannBoundaryConditions< MeshType, ConstantFunction, Real, Index > BoundaryConditions;
+         typedef tnlMeanCurvatureFlowProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Solver;
+         SolverStarter solverStarter;
+         return solverStarter.template run< Solver >( parameters );
+      }
+      typedef tnlVector< Real, Device, Index > VectorType;
+      if( boundaryConditionsType == "dirichlet" )
+      {
+         typedef tnlDirichletBoundaryConditions< MeshType, VectorType, Real, Index > BoundaryConditions;
+         typedef tnlMeanCurvatureFlowProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Solver;
+         SolverStarter solverStarter;
+         return solverStarter.template run< Solver >( parameters );
+      }
+      typedef tnlNeumannBoundaryConditions< MeshType, VectorType, Real, Index > BoundaryConditions;
+      typedef tnlMeanCurvatureFlowProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Solver;
+      SolverStarter solverStarter;
+      return solverStarter.template run< Solver >( parameters );
+   };
+};
+
+int main( int argc, char* argv[] )
+{
+   tnlSolver< meanCurvatureFlowSetter, meanCurvatureFlowConfig, BuildConfig > solver;
+   if( ! solver. run( argc, argv ) )
+      return EXIT_FAILURE;
+   return EXIT_SUCCESS;
+}
+
+#endif /* TNL_MEAN_CURVATIVE_FLOW_H_ */
diff --git a/examples/mean-curvature-flow/tnl-run-mean-curvature-flow b/examples/mean-curvature-flow/tnl-run-mean-curvature-flow
new file mode 100644
index 0000000000000000000000000000000000000000..c44ad55786ea7790f1e9bacc670b917cb1faacf7
--- /dev/null
+++ b/examples/mean-curvature-flow/tnl-run-mean-curvature-flow
@@ -0,0 +1,90 @@
+#!/bin/bash
+
+dofSize=64
+dimension=2;
+proportions=1
+
+analyticFunction="exp-bump"
+timeFunction="cosinus"
+
+amplitude=1.0
+waveLength=1.0
+waveLengthX=1.0
+waveLengthY=1.0
+waveLengthZ=1.0
+wavesNumber=0.0
+wavesNumberX=0.0
+wavesNumberY=0.0
+wavesNumberZ=0.0
+phase=0.0
+phaseX=0.0
+phaseY=0.0
+phaseZ=0.0
+sigma=1.0
+
+tnl-grid-setup --dimensions ${dimension} \
+               --proportions-x ${proportions} \
+               --proportions-y ${proportions} \
+               --proportions-z ${proportions} \
+               --origin-x 0 \
+               --origin-y 0 \
+               --origin-z 0 \
+               --size-x ${dofSize} \
+               --size-y ${dofSize} \
+               --size-z ${dofSize} \
+               
+tnl-init --mesh mesh.tnl \
+         --test-function ${analyticFunction} \
+         --output-file initial.tnl \
+         --amplitude ${amplitude} \
+         --wave-length ${waveLength} \
+         --wave-length-x ${waveLengthX} \
+         --wave-length-y ${waveLengthY} \
+         --wave-length-z ${waveLengthZ} \
+             --waves-number ${wavesNumber} \
+             --waves-number-x ${wavesNumberX} \
+             --waves-number-y ${wavesNumberY} \
+             --waves-number-z ${wavesNumberZ} \
+             --phase ${phase} \
+             --phase-x ${phaseX} \
+             --phase-y ${phaseY} \
+             --phase-z ${phaseZ} \
+             --sigma ${sigma} \
+
+tnl-mean-curvature-flow --time-discretisation explicit \
+                  --boundary-conditions-type dirichlet \
+                  --boundary-conditions-constant 0.5 \
+                  --discrete-solver merson \
+                  --snapshot-period 0.0005 \
+                  --final-time 0.1 \
+              
+tnl-view --mesh mesh.tnl \
+         --input-files *.tnl \ 
+
+seznam=`ls u-*.gplt`
+
+for fname in $seznam ; do
+   echo "Drawing $fname"
+gnuplot << EOF
+    set terminal unknown
+    #set view 33,33 #3D
+    #unset xtics 
+    #unset ytics
+    #unset ztics
+    unset border
+    set output '$fname.png'
+    set yrange [-1.2:1.2]
+    set zrange [0.4:1.1]    
+    set terminal png
+    set title "Numerical solution" 
+    splot '$fname' with line 
+EOF
+done
+
+
+mencoder "mf://u-*.png" -mf fps=22 -o diffusion.avi -ovc lavc -lavcopts vcodec=mpeg4
+
+#rm *.png
+#rm *.tnl         
+#rm *.gplt      
+              
diff --git a/examples/mean-curvature-flow/tnl-run-mean-curvature-flow-eoc-test b/examples/mean-curvature-flow/tnl-run-mean-curvature-flow-eoc-test
new file mode 100755
index 0000000000000000000000000000000000000000..abf3827cd5110cd42073eaba0b680a1343d45d3f
--- /dev/null
+++ b/examples/mean-curvature-flow/tnl-run-mean-curvature-flow-eoc-test
@@ -0,0 +1,200 @@
+#!/bin/bash
+
+device="host"
+dimensions="2D"
+sizes2D="16 32 64"
+testFunctions="sin-bumps"
+snapshotPeriod=0.1
+finalTime=1.0
+timeDependence="cosine"
+solverName="mean-curvature-flow-eoc"
+#solverName="gdb --args tnl-heat-equation-eoc-test-dbg"
+#
+
+setupTestFunction()
+{
+   testFunction=$1
+#   if test x${testFunction} = "xexp-bump";
+#   then
+      origin=0
+      proportions=1
+      amplitude=1.0
+      waveLength=1.0
+      waveLengthX=1.0
+      waveLengthY=1.0
+      waveLengthZ=1.0
+      wavesNumber=0.0
+      wavesNumberX=0.0
+      wavesNumberY=0.0
+      wavesNumberZ=0.0
+      phase=0.0
+      phaseX=0.0
+      phaseY=0.0
+      phaseZ=0.0
+      sigma=0.5
+#   fi
+}
+
+setupGrid()
+{
+   dimensions=$1
+   gridSize=$2
+   tnl-grid-setup --dimensions ${dimensions} \
+                  --origin-x ${origin} \
+                  --origin-y ${origin} \
+                  --origin-z ${origin} \
+                  --proportions-x ${proportions} \
+                  --proportions-y ${proportions} \
+                  --proportions-z ${proportions} \
+                  --size-x ${gridSize} \
+                  --size-y ${gridSize} \
+                  --size-z ${gridSize} 
+}
+
+setInitialCondition()
+{
+   testFunction=$1
+   tnl-init --test-function ${testFunction} \
+            --output-file exact-u.tnl \
+            --amplitude ${amplitude} \
+            --wave-length ${waveLength} \
+            --wave-length-x ${waveLengthX} \
+            --wave-length-y ${waveLengthY} \
+            --wave-length-z ${waveLengthZ} \
+            --waves-number ${wavesNumber} \
+            --waves-number-x ${wavesNumberX} \
+            --waves-number-y ${wavesNumberY} \
+            --waves-number-z ${wavesNumberZ} \
+            --phase ${phase} \
+            --phase-x ${phaseX} \
+            --phase-y ${phaseY} \
+            --phase-z ${phaseZ} \
+            --sigma ${sigma} \
+            --time-dependence ${timeDependence} \
+            --snapshot-period ${snapshotPeriod} \
+            --final-time ${finalTime}
+}
+
+solve()
+{
+   timeDiscretisation=$1
+   discreteSolver=$2
+   ${solverName} --device ${device} \
+                 --mesh mesh.tnl \
+                 --initial-condition exact-u-00000.tnl \
+                 --time-discretisation ${timeDiscretisation} \
+                 --time-step 10 \
+                 --time-step-order 2 \
+                 --discrete-solver ${discreteSolver} \
+                 --merson-adaptivity 1.0e-7 \
+                 --sor-omega 1.95 \
+                 --gmres-restarting 10 \
+                 --min-iterations 20 \
+                 --convergence-residue 1.0e-12 \
+                 --test-function ${testFunction}\
+                 --amplitude ${amplitude} \
+                 --wave-length ${waveLength} \
+                 --wave-length-x ${waveLengthX} \
+                 --wave-length-y ${waveLengthY} \
+                 --wave-length-z ${waveLengthZ} \
+                 --waves-number ${wavesNumber} \
+                 --waves-number-x ${wavesNumberX} \
+                 --waves-number-y ${wavesNumberY} \
+                 --waves-number-z ${wavesNumberZ} \
+                 --phase ${phase} \
+                 --phase-x ${phaseX} \
+                 --phase-y ${phaseY} \
+                 --phase-z ${phaseZ} \
+                 --sigma ${sigma} \
+                 --time-dependence ${timeDependence} \
+                 --snapshot-period ${snapshotPeriod} \
+                 --final-time ${finalTime}
+}
+               
+computeError()
+{
+   tnl-diff --mesh mesh.tnl \
+            --input-files exact-u-*.tnl u-*.tnl \
+            --mode halves \
+            --snapshot-period ${snapshotPeriod} \
+            --output-file errors.txt \
+            --write-difference yes
+}
+
+runTest()
+{
+   for testFunction in ${testFunctions};
+   do
+      mkdir -p ${testFunction}
+      cd ${testFunction}
+      setupTestFunction ${testFunction}
+      
+      for dim in ${dimensions};
+      do
+         mkdir -p $dim
+         cd ${dim}
+         if test $dim = 1D;
+         then 
+            sizes=$sizes1D
+         fi
+         if test $dim = 2D;
+         then 
+            sizes=$sizes2D
+         fi
+         if test $dim = 3D;
+         then 
+            sizes=$sizes3D
+         fi
+         
+         lastSize=""
+         for size in $sizes;
+         do
+            mkdir -p $size
+            cd $size
+            echo ""
+            echo ""
+            echo ""
+            if test ! -f computation-done;
+            then
+               touch computation-in-progress
+               echo "========================================================================="
+               echo "===                   SETTING UP THE GRID                             ==="
+               echo "========================================================================="
+               setupGrid $dim $size 
+               echo "========================================================================="
+               echo "===                WRITING THE EXACT SOLUTION                         ==="
+               echo "========================================================================="
+               setInitialCondition $testFunction
+               echo "========================================================================="
+               echo "===                   STARTING THE SOLVER                             ==="
+               echo "========================================================================="
+               solve explicit merson
+               #solve semi-implicit gmres
+               mv computation-in-progress computation-done
+            fi            
+            echo "========================================================================="
+            echo "===                   COMPUTING THE ERROR                             ==="
+            echo "========================================================================="
+            computeError
+            echo "========================================================================="
+            echo "===                     COMPUTING THE EOC                             ==="            
+            echo "========================================================================="
+            if test ! x$lastSize = x;
+            then
+               tnl-err2eoc ../$lastSize/errors.txt errors.txt
+            fi
+            echo "========================================================================="
+            echo "===                     COMPUTATION DONE                              ==="            
+            echo "========================================================================="
+            cd ..
+            lastSize=$size
+         done
+
+         cd ..
+      done
+      cd ..
+   done
+}
+
+runTest
+ 
diff --git a/examples/mean-curvature-flow/tnlBackwardFiniteDifference.h b/examples/mean-curvature-flow/tnlBackwardFiniteDifference.h
new file mode 100644
index 0000000000000000000000000000000000000000..b685d9df8ff092ecc9470c8139e774125a10fc88
--- /dev/null
+++ b/examples/mean-curvature-flow/tnlBackwardFiniteDifference.h
@@ -0,0 +1,128 @@
+#ifndef TNLBACKWARDFINITEDIFFERENCE_H
+#define	TNLBACKWARDFINITEDIFFERENCE_H
+
+#include <core/vectors/tnlVector.h>
+#include <mesh/tnlGrid.h>
+
+template< typename Mesh,
+          typename Real = typename Mesh::RealType,
+          typename Index = typename Mesh::IndexType >
+class tnlBackwardFiniteDifference
+{
+ 
+};
+
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+class tnlBackwardFiniteDifference< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index >
+{
+   public: 
+   
+   typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;
+   typedef typename MeshType::CoordinatesType CoordinatesType;
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+
+   static tnlString getType();
+
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif   
+   Real getValueX( const MeshType& mesh,
+                   const IndexType cellIndex,
+                   const Vector& u) const 
+   {}
+   
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif   
+   Real getValueY( const MeshType& mesh,
+                   const IndexType cellIndex,
+                   const Vector& u) const
+   {}
+};
+
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+class tnlBackwardFiniteDifference< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >
+{
+   public: 
+   
+   typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;
+   typedef typename MeshType::CoordinatesType CoordinatesType;
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+
+   static tnlString getType(); 
+
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif   
+   Real getValueX( const MeshType& mesh,
+                   const IndexType cellIndex,
+                   const Vector& u) const;
+   
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif   
+   Real getValueY( const MeshType& mesh,
+                   const IndexType cellIndex,
+                   const Vector& u) const;
+};
+
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+class tnlBackwardFiniteDifference< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >
+{
+   public: 
+   
+   typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;
+   typedef typename MeshType::CoordinatesType CoordinatesType;
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+
+   static tnlString getType();
+
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif   
+   Real getValueX( const MeshType& mesh,
+                   const IndexType cellIndex,
+                   const Vector& u) const 
+   {}
+   
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif   
+   Real getValueY( const MeshType& mesh,
+                   const IndexType cellIndex,
+                   const Vector& u) const
+   {}
+};
+
+
+#include "tnlBackwardFiniteDifference_impl.h"
+
+
+#endif	/* TNLBACKWARDFINITEDIFFERENCE_H */
diff --git a/examples/mean-curvature-flow/tnlBackwardFiniteDifference_impl.h b/examples/mean-curvature-flow/tnlBackwardFiniteDifference_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..fe3c3e0774fb74cdbfab26d501d0585ba5e4bc00
--- /dev/null
+++ b/examples/mean-curvature-flow/tnlBackwardFiniteDifference_impl.h
@@ -0,0 +1,91 @@
+
+#ifndef TNLBACKWARDFINITEDIFFERENCE_IMPL_H
+#define	TNLBACKWARDFINITEDIFFERENCE_IMPL_H
+
+#include "tnlBackwardFiniteDifference.h"
+#include <mesh/tnlGrid.h>
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+tnlString
+tnlBackwardFiniteDifference< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+   return tnlString( "tnlBackwardFiniteDifference< " ) +
+          MeshType::getType() + ", " +
+          ::getType< Real >() + ", " +
+          ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+tnlString
+tnlBackwardFiniteDifference< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+   return tnlString( "tnlBackwardFiniteDifference< " ) +
+          MeshType::getType() + ", " +
+          ::getType< Real >() + ", " +
+          ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlBackwardFiniteDifference< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getValueX( const MeshType& mesh,
+          const IndexType cellIndex,
+          const Vector& u) const
+{
+   return ( - u[ mesh.template getCellNextToCell< -1,0 >( cellIndex ) ]
+            + u[ cellIndex ]) * mesh.getHxInverse();   
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlBackwardFiniteDifference< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getValueY( const MeshType& mesh,
+          const IndexType cellIndex,
+          const Vector& u) const
+{
+   return ( - u[ mesh.template getCellNextToCell< 0,-1 >( cellIndex ) ]
+            + u[ cellIndex ]) * mesh.getHyInverse();   
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+tnlString
+tnlBackwardFiniteDifference< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+   return tnlString( "tnlBackwardFiniteDifference< " ) +
+          MeshType::getType() + ", " +
+          ::getType< Real >() + ", " +
+          ::getType< Index >() + " >";
+}
+
+#endif	/* TNLBACKWARDFINITEDIFFERENCE_IMPL_H */
diff --git a/examples/mean-curvature-flow/tnlExactFlowDiffusion.h b/examples/mean-curvature-flow/tnlExactFlowDiffusion.h
new file mode 100644
index 0000000000000000000000000000000000000000..43d206f255ee49766be207596e61ea0a780e6b0b
--- /dev/null
+++ b/examples/mean-curvature-flow/tnlExactFlowDiffusion.h
@@ -0,0 +1,108 @@
+/***************************************************************************
+                          tnlExactLinearDiffusion.h  -  description
+                             -------------------
+    begin                : Aug 8, 2014
+    copyright            : (C) 2014 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef TNLEXACTFLOWDIFFUSION_H_
+#define TNLEXACTFLOWDIFFUSION_H_
+
+#include <functions/tnlFunctionType.h>
+
+template< int Dimensions >
+class tnlExactFlowDiffusion
+{};
+
+template<>
+class tnlExactFlowDiffusion< 1 >
+{
+   public:
+
+      enum { Dimensions = 1 };
+
+      static tnlString getType();
+   
+#ifdef HAVE_NOT_CXX11      
+      template< typename Function, typename Vertex, typename Real >
+#else   
+      template< typename Function, typename Vertex, typename Real = typename Vertex::RealType >
+#endif      
+#ifdef HAVE_CUDA
+      __device__ __host__
+#endif
+      static Real getValue( const Function& function,
+                            const Vertex& v,
+                            const Real& time = 0.0 )
+      {
+         return 0;
+      }
+};
+
+template<>
+class tnlExactFlowDiffusion< 2 >
+{
+   public:
+
+      enum { Dimensions = 2 };
+
+      static tnlString getType();
+
+#ifdef HAVE_NOT_CXX11      
+      template< typename Function, typename Vertex, typename Real >
+#else   
+      template< typename Function, typename Vertex, typename Real = typename Vertex::RealType >
+#endif
+#ifdef HAVE_CUDA
+      __device__ __host__
+#endif      
+      static Real getValue( const Function& function,
+                            const Vertex& v,
+                            const Real& time = 0.0 );
+};
+
+template<>
+class tnlExactFlowDiffusion< 3 >
+{
+   public:
+
+      enum { Dimensions = 3 };
+
+      static tnlString getType();
+
+#ifdef HAVE_NOT_CXX11      
+      template< typename Function, typename Vertex, typename Real >
+#else   
+      template< typename Function, typename Vertex, typename Real = typename Vertex::RealType >
+#endif
+#ifdef HAVE_CUDA
+      __device__ __host__
+#endif
+      static Real getValue( const Function& function,
+                            const Vertex& v,
+                            const Real& time = 0.0 )
+      {
+         return 0;
+      }
+};
+
+template< int Dimensions >
+class tnlFunctionType< tnlExactFlowDiffusion< Dimensions > >
+{
+   public:
+      enum { Type = tnlAnalyticFunction };
+};
+
+#include "tnlExactFlowDiffusion_impl.h"
+
+#endif /* TNLEXACTFLOWDIFFUSION_H_ */
diff --git a/examples/mean-curvature-flow/tnlExactFlowDiffusion_impl.h b/examples/mean-curvature-flow/tnlExactFlowDiffusion_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..05d623831eab89af3364a532c027efe3b3e3dbca
--- /dev/null
+++ b/examples/mean-curvature-flow/tnlExactFlowDiffusion_impl.h
@@ -0,0 +1,59 @@
+/***************************************************************************
+                          tnlExactLinearDiffusion_impl.h  -  description
+                             -------------------
+    begin                : Aug 8, 2014
+    copyright            : (C) 2014 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef TNLEXACTFLOWDIFFUSION_IMPL_H_
+#define TNLEXACTFLOWDIFFUSION_IMPL_H_
+
+tnlString
+tnlExactFlowDiffusion< 1 >::
+getType()
+{
+   return "tnlExactFlowDiffusion< 1 >";
+}
+
+tnlString
+tnlExactFlowDiffusion< 2 >::
+getType()
+{
+   return "tnlExactFlowDiffusion< 2 >";
+}
+
+template< typename Function, typename Vertex, typename Real >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+Real
+tnlExactFlowDiffusion< 2 >::
+getValue( const Function& function,
+          const Vertex& v,
+          const Real& time )
+{
+   return (function.template getValue< 2, 0, 0, Vertex >( v, time )*(1.0+ function.template getValue< 0, 1, 0, Vertex >( v, time )*function.template getValue< 0, 1, 0, Vertex >( v, time ))
+          - 2.0*function.template getValue< 0, 1, 0, Vertex >( v, time )*function.template getValue< 1, 0, 0, Vertex >( v, time )*function.template getValue< 1, 1, 0, Vertex >( v, time )
+          + function.template getValue< 0, 2, 0, Vertex >( v, time )*(1.0+ function.template getValue< 1, 0, 0, Vertex >( v, time )*function.template getValue< 1, 0, 0, Vertex >( v, time )))
+          / (1.0+function.template getValue< 1, 0, 0, Vertex >( v, time )*function.template getValue< 1, 0, 0, Vertex >( v, time )
+          +function.template getValue< 0, 1, 0, Vertex >( v, time )*function.template getValue< 0, 1, 0, Vertex >( v, time ));
+}
+
+tnlString
+tnlExactFlowDiffusion< 3 >::
+getType()
+{
+   return "tnlExactFlowDiffusion< 3 >";
+}
+
+#endif /* TNLEXACTFLOWDIFFUSION_IMPL_H_ */
diff --git a/examples/mean-curvature-flow/tnlForwardFiniteDifference.h b/examples/mean-curvature-flow/tnlForwardFiniteDifference.h
new file mode 100644
index 0000000000000000000000000000000000000000..f207879ffffdf3face39fcb0315a760671d6b682
--- /dev/null
+++ b/examples/mean-curvature-flow/tnlForwardFiniteDifference.h
@@ -0,0 +1,128 @@
+#ifndef TNLFORWARDFINITEDIFFERENCE_H
+#define	TNLFORWARDFINITEDIFFERENCE_H
+
+#include <core/vectors/tnlVector.h>
+#include <mesh/tnlGrid.h>
+
+template< typename Mesh,
+          typename Real = typename Mesh::RealType,
+          typename Index = typename Mesh::IndexType >
+class tnlForwardFiniteDifference
+{
+ 
+};
+
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+class tnlForwardFiniteDifference< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index >
+{
+   public: 
+   
+   typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;
+   typedef typename MeshType::CoordinatesType CoordinatesType;
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+
+   static tnlString getType();
+
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif   
+   Real getValueX( const MeshType& mesh,
+                   const IndexType cellIndex,
+                   const Vector& u) const 
+   {}
+   
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif   
+   Real getValueY( const MeshType& mesh,
+                   const IndexType cellIndex,
+                   const Vector& u) const
+   {}
+};
+
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+class tnlForwardFiniteDifference< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >
+{
+   public: 
+   
+   typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;
+   typedef typename MeshType::CoordinatesType CoordinatesType;
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+
+   static tnlString getType(); 
+
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif   
+   Real getValueX( const MeshType& mesh,
+                   const IndexType cellIndex,
+                   const Vector& u) const;
+   
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif   
+   Real getValueY( const MeshType& mesh,
+                   const IndexType cellIndex,
+                   const Vector& u) const;
+};
+
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+class tnlForwardFiniteDifference< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >
+{
+   public: 
+   
+   typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;
+   typedef typename MeshType::CoordinatesType CoordinatesType;
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+
+   static tnlString getType();
+
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif   
+   Real getValueX( const MeshType& mesh,
+                   const IndexType cellIndex,
+                   const Vector& u) const 
+   {}
+   
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif   
+   Real getValueY( const MeshType& mesh,
+                   const IndexType cellIndex,
+                   const Vector& u) const
+   {}
+};
+
+
+#include "tnlForwardFiniteDifference_impl.h"
+
+
+#endif	/* TNLFORWARDFINITEDIFFERENCE_H */
diff --git a/examples/mean-curvature-flow/tnlForwardFiniteDifference_impl.h b/examples/mean-curvature-flow/tnlForwardFiniteDifference_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..419adb4786b90eba4c4689aaa5bfd3323d5c0aef
--- /dev/null
+++ b/examples/mean-curvature-flow/tnlForwardFiniteDifference_impl.h
@@ -0,0 +1,91 @@
+
+#ifndef TNLFORWARDFINITEDIFFERENCE_IMPL_H
+#define	TNLFORWARDFINITEDIFFERENCE_IMPL_H
+
+#include "tnlForwardFiniteDifference.h"
+#include <mesh/tnlGrid.h>
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+tnlString
+tnlForwardFiniteDifference< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+   return tnlString( "tnlForwardFiniteDifference< " ) +
+          MeshType::getType() + ", " +
+          ::getType< Real >() + ", " +
+          ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+tnlString
+tnlForwardFiniteDifference< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+   return tnlString( "tnlForwardFiniteDifference< " ) +
+          MeshType::getType() + ", " +
+          ::getType< Real >() + ", " +
+          ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlForwardFiniteDifference< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getValueX( const MeshType& mesh,
+          const IndexType cellIndex,
+          const Vector& u) const
+{
+   return ( u[ mesh.template getCellNextToCell< 1,0 >( cellIndex ) ]
+            - u[ cellIndex ]) * mesh.getHxInverse();
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlForwardFiniteDifference< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getValueY( const MeshType& mesh,
+          const IndexType cellIndex,
+          const Vector& u) const
+{
+   return ( u[ mesh.template getCellNextToCell< 0,1 >( cellIndex ) ]
+            - u[ cellIndex ]) * mesh.getHyInverse();   
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+tnlString
+tnlForwardFiniteDifference< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+   return tnlString( "tnlForwardFiniteDifference< " ) +
+          MeshType::getType() + ", " +
+          ::getType< Real >() + ", " +
+          ::getType< Index >() + " >";
+}
+
+#endif	/* TNLFORWARDFINITEDIFFERENCE_IMPL_H */
diff --git a/examples/mean-curvature-flow/tnlMeanCurvatureFlowDiffusion.h b/examples/mean-curvature-flow/tnlMeanCurvatureFlowDiffusion.h
new file mode 100644
index 0000000000000000000000000000000000000000..2800f91511ac79ffa27ba82207c4016acc8549bd
--- /dev/null
+++ b/examples/mean-curvature-flow/tnlMeanCurvatureFlowDiffusion.h
@@ -0,0 +1,211 @@
+#ifndef TNLMEANCURVATIVEFLOWDIFFUSION_H
+#define	TNLMEANCURVATIVEFLOWDIFFUSION_H
+
+#include <core/vectors/tnlVector.h>
+#include "tnlForwardFiniteDifference.h"
+#include "tnlBackwardFiniteDifference.h"
+#include "tnlOperatorQ.h"
+#include <mesh/tnlGrid.h>
+
+template< typename Mesh,
+          typename Real = typename Mesh::RealType,
+          typename Index = typename Mesh::IndexType >
+class tnlMeanCurvatureFlowDiffusion
+{
+ 
+};
+
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+class tnlMeanCurvatureFlowDiffusion< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index >
+{
+   public: 
+   
+   typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;
+   typedef typename MeshType::CoordinatesType CoordinatesType;
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+
+   static tnlString getType();
+   
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   Real getValue( const MeshType& mesh,
+                  const IndexType cellIndex,
+                  const CoordinatesType& coordinates,
+                  const Vector& u,
+                  const RealType& time) const 
+   {}
+   
+//   void setupDofs ( const MeshType& mesh )
+//   {}
+//   
+//   template< typename Vector >
+//   void computeFirstGradient( const MeshType& mesh,
+//                              const RealType& time,
+//                              const Vector u)
+//   {}
+   
+
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   Index getLinearSystemRowLength( const MeshType& mesh,
+                                   const IndexType& index,
+                                   const CoordinatesType& coordinates ) const
+   {}
+
+   template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+      void updateLinearSystem( const RealType& time,
+                               const RealType& tau,
+                               const MeshType& mesh,
+                               const IndexType& index,
+                               const CoordinatesType& coordinates,
+                               Vector& u,
+                               Vector& b,
+                               MatrixRow& matrixRow ) const
+   {}
+
+};
+
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+class tnlMeanCurvatureFlowDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >
+{
+   public: 
+   
+   typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;
+   typedef typename MeshType::CoordinatesType CoordinatesType;
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+//   typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
+   
+   static tnlString getType();
+   
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   Real getValue( const MeshType& mesh,
+                  const IndexType cellIndex,
+                  const CoordinatesType& coordinates,
+                  const Vector& u,
+                  const RealType& time) const;
+   
+//   void setupDofs ( const MeshType& mesh );
+//   
+//   template< typename Vector >
+//   void computeFirstGradient( const MeshType& mesh,
+//                              const RealType& time,
+//                              const Vector u);
+   
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   Index getLinearSystemRowLength( const MeshType& mesh,
+                                   const IndexType& index,
+                                   const CoordinatesType& coordinates ) const
+   {}
+
+   template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+      void updateLinearSystem( const RealType& time,
+                               const RealType& tau,
+                               const MeshType& mesh,
+                               const IndexType& index,
+                               const CoordinatesType& coordinates,
+                               Vector& u,
+                               Vector& b,
+                               MatrixRow& matrixRow ) const
+   {}
+   
+   protected:
+      
+   tnlForwardFiniteDifference<MeshType> fDifference;
+   tnlBackwardFiniteDifference<MeshType> bDifference;
+   tnlOperatorQ<MeshType> operatorQ;
+};
+
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+class tnlMeanCurvatureFlowDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >
+{
+   public: 
+   
+   typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;
+   typedef typename MeshType::CoordinatesType CoordinatesType;
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+
+   static tnlString getType();
+   
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   Real getValue( const MeshType& mesh,
+                  const IndexType cellIndex,
+                  const CoordinatesType& coordinates,
+                  const Vector& u,
+                  const RealType& time) const
+   {}
+   
+//   void setupDofs ( const MeshType& mesh )
+//   {}
+//   
+//   template< typename Vector >
+//   void computeFirstGradient( const MeshType& mesh,
+//                              const RealType& time,
+//                              const Vector u )
+//   {}
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+   Index getLinearSystemRowLength( const MeshType& mesh,
+                                   const IndexType& index,
+                                   const CoordinatesType& coordinates ) const 
+   {}
+
+   template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif
+      void updateLinearSystem( const RealType& time,
+                               const RealType& tau,
+                               const MeshType& mesh,
+                               const IndexType& index,
+                               const CoordinatesType& coordinates,
+                               Vector& u,
+                               Vector& b,
+                               MatrixRow& matrixRow ) const 
+   {}
+
+};
+
+
+#include "tnlMeanCurvatureFlowDiffusion_impl.h"
+
+
+#endif	/* TNLMEANCURVATIVEFLOWDIFFUSION_H */
diff --git a/examples/mean-curvature-flow/tnlMeanCurvatureFlowDiffusion_impl.h b/examples/mean-curvature-flow/tnlMeanCurvatureFlowDiffusion_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..061af10f4000e73c79101460c41fdc4089b9bd68
--- /dev/null
+++ b/examples/mean-curvature-flow/tnlMeanCurvatureFlowDiffusion_impl.h
@@ -0,0 +1,104 @@
+
+#ifndef TNLMEANCURVATIVEFLOWDIFFUSION_IMPL_H
+#define	TNLMEANCURVATIVEFLOWDIFFUSION_IMPL_H
+
+#include "tnlMeanCurvatureFlowDiffusion.h"
+#include "tnlForwardFiniteDifference.h"
+#include "tnlBackwardFiniteDifference.h"
+
+#include <mesh/tnlGrid.h>
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+tnlString
+tnlMeanCurvatureFlowDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+   return tnlString( "tnlForwardFiniteDifference< " ) +
+          MeshType::getType() + ", " +
+          ::getType< Real >() + ", " +
+          ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+tnlString
+tnlMeanCurvatureFlowDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+   return tnlString( "tnlForwardFiniteDifference< " ) +
+          MeshType::getType() + ", " +
+          ::getType< Real >() + ", " +
+          ::getType< Index >() + " >";
+}
+
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlMeanCurvatureFlowDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getValue( const MeshType& mesh,
+          const IndexType cellIndex,
+          const CoordinatesType& coordinates,
+          const Vector& u,
+          const Real& time ) const
+{
+   return operatorQ.getValueStriped(mesh,cellIndex,u)*((fDifference.getValueX(mesh,cellIndex,u)/operatorQ.getValue(mesh,cellIndex,u)
+          -bDifference.getValueX(mesh,cellIndex,u)/operatorQ.getValue(mesh,mesh.template getCellNextToCell<-1,0>(cellIndex),u))
+          *mesh.getHxInverse()+(fDifference.getValueY(mesh,cellIndex,u)/operatorQ.getValue(mesh,cellIndex,u)
+          -bDifference.getValueY(mesh,cellIndex,u)/operatorQ.getValue(mesh,mesh.template getCellNextToCell<0,-1>(cellIndex),u))*mesh.getHyInverse());
+}
+//template< typename MeshReal,
+//          typename Device,
+//          typename MeshIndex,
+//          typename Real,
+//          typename Index >
+//void tnlMeanCurvatureFlowDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::setupDofs ( const MeshType& mesh )
+//{
+//   dofs.setSize(mesh.getNumberOfCells());
+//}
+//
+//template< typename MeshReal,
+//          typename Device,
+//          typename MeshIndex,
+//          typename Real,
+//          typename Index >
+//template< typename Vector>
+//void 
+//tnlMeanCurvatureFlowDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+//computeFirstGradient( const MeshType& mesh,
+//                      const RealType& time,
+//                      const Vector u)
+//{
+//   
+//}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+tnlString
+tnlMeanCurvatureFlowDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+   return tnlString( "tnlForwardFiniteDifference< " ) +
+          MeshType::getType() + ", " +
+          ::getType< Real >() + ", " +
+          ::getType< Index >() + " >";
+}
+
+#endif	/* TNLMEANCURVATIVEFLOWDIFFUSION_IMPL_H */
diff --git a/examples/mean-curvature-flow/tnlMeanCurvatureFlowEocProblem.h b/examples/mean-curvature-flow/tnlMeanCurvatureFlowEocProblem.h
new file mode 100644
index 0000000000000000000000000000000000000000..5c94ec71000b1e0c075e280bbed7de0c122ae936
--- /dev/null
+++ b/examples/mean-curvature-flow/tnlMeanCurvatureFlowEocProblem.h
@@ -0,0 +1,40 @@
+/***************************************************************************
+                          tnlHeatEquationEocProblem.h  -  description
+                             -------------------
+    begin                : Nov 22, 2014
+    copyright            : (C) 2014 by 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 TNLMEANCURVATUREFLOWEOCPROBLEM_H_
+#define TNLMEANCURVATUREFLOWEOCPROBLEM_H_
+
+#include "tnlMeanCurvatureFlowProblem.h"
+
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator = tnlMeanCurvatureFlowDiffusion< Mesh,
+                                                              typename BoundaryCondition::RealType > >
+class tnlMeanCurvatureFlowEocProblem : public tnlMeanCurvatureFlowProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >
+{
+   public:
+
+      static tnlString getTypeStatic();
+
+      bool setup( const tnlParameterContainer& parameters );
+};
+
+#include "tnlMeanCurvatureFlowEocProblem_impl.h"
+
+#endif /* TNLMEANCURVATUREFLOWEOCPROBLEM_H_ */
diff --git a/examples/mean-curvature-flow/tnlMeanCurvatureFlowEocProblem_impl.h b/examples/mean-curvature-flow/tnlMeanCurvatureFlowEocProblem_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..dd7e1fb23036a5f3cde940c888df476d6a36e731
--- /dev/null
+++ b/examples/mean-curvature-flow/tnlMeanCurvatureFlowEocProblem_impl.h
@@ -0,0 +1,46 @@
+/***************************************************************************
+                          tnlHeatEquationEocProblem_impl.h  -  description
+                             -------------------
+    begin                : Nov 22, 2014
+    copyright            : (C) 2014 by 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 TNLMEANCURVATUREFLOWEOCPROBLEM_IMPL_H_
+#define TNLMEANCURVATUREFLOWEOCPROBLEM_IMPL_H_
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+tnlString
+tnlMeanCurvatureFlowEocProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+getTypeStatic()
+{
+   return tnlString( "tnlHeatEquationEocProblem< " ) + Mesh :: getTypeStatic() + " >";
+}
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+bool
+tnlMeanCurvatureFlowEocProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator  >::
+setup( const tnlParameterContainer& parameters )
+{
+   if( ! this->boundaryCondition.setup( parameters ) ||
+       ! this->rightHandSide.setup( parameters ) )
+      return false;
+   return true;
+}
+
+#endif /* TNLMEANCURVATUREFLOWEOCPROBLEM_IMPL_H_ */
diff --git a/examples/mean-curvature-flow/tnlMeanCurvatureFlowEocRhs.h b/examples/mean-curvature-flow/tnlMeanCurvatureFlowEocRhs.h
new file mode 100644
index 0000000000000000000000000000000000000000..60910448ee9da96d8fe413b8f90b943b39c5a761
--- /dev/null
+++ b/examples/mean-curvature-flow/tnlMeanCurvatureFlowEocRhs.h
@@ -0,0 +1,67 @@
+/***************************************************************************
+                          tnlHeatEquationEocRhs.h  -  description
+                             -------------------
+    begin                : Sep 8, 2014
+    copyright            : (C) 2014 by 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 TNLMEANCURVATUREFLOWEOCRHS_H_
+#define TNLMEANCURVATUREFLOWEOCRHS_H_
+
+#include <functions/tnlFunctionType.h>
+
+template< typename ExactOperator,
+          typename TestFunction >
+class tnlMeanCurvatureFlowEocRhs
+{
+   public:
+
+      typedef ExactOperator ExactOperatorType;
+      typedef TestFunction TestFunctionType;
+
+      bool setup( const tnlParameterContainer& parameters,
+                  const tnlString& prefix = "" )
+      {
+         if( ! testFunction.setup( parameters, prefix ) )
+            return false;
+         return true;
+      };
+
+      template< typename Vertex,
+                typename Real >
+#ifdef HAVE_CUDA
+      __device__ __host__
+#endif
+      Real getValue( const Vertex& vertex,
+                     const Real& time ) const
+      {
+         return testFunction.getTimeDerivative( vertex, time )
+                - exactOperator.getValue( testFunction, vertex, time );
+      };
+
+   protected:
+      ExactOperator exactOperator;
+
+      TestFunction testFunction;
+};
+
+template< typename ExactOperator,
+          typename TestFunction >
+class tnlFunctionType< tnlMeanCurvatureFlowEocRhs< ExactOperator, TestFunction > >
+{
+   public:
+
+      enum { Type = tnlAnalyticFunction };
+};
+
+#endif /* TNLMEANCURVATUREFLOWEOCRHS_H_ */
diff --git a/examples/mean-curvature-flow/tnlMeanCurvatureFlowProblem.h b/examples/mean-curvature-flow/tnlMeanCurvatureFlowProblem.h
new file mode 100644
index 0000000000000000000000000000000000000000..f2686b2dc2a927b7b84194e34192485e4cb85471
--- /dev/null
+++ b/examples/mean-curvature-flow/tnlMeanCurvatureFlowProblem.h
@@ -0,0 +1,102 @@
+/***************************************************************************
+                          tnlHeatEquationProblem.h  -  description
+                             -------------------
+    begin                : Feb 23, 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 TNLMEANCURVATUREFLOWPROBLEM_H_
+#define TNLMEANCURVATUREFLOWPROBLEM_H_
+
+#include "tnlMeanCurvatureFlowDiffusion.h"
+#include <problems/tnlPDEProblem.h>
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator = tnlMeanCurvatureFlowDiffusion< Mesh,
+                                                              typename BoundaryCondition::RealType > >
+class tnlMeanCurvatureFlowProblem : public tnlPDEProblem< Mesh,
+                                                     typename DifferentialOperator::RealType,
+                                                     typename Mesh::DeviceType,
+                                                     typename DifferentialOperator::IndexType  >
+{
+   public:
+
+      typedef typename DifferentialOperator::RealType RealType;
+      typedef typename Mesh::DeviceType DeviceType;
+      typedef typename DifferentialOperator::IndexType IndexType;
+      typedef tnlPDEProblem< Mesh, RealType, DeviceType, IndexType > BaseType;
+
+      using typename BaseType::MeshType;
+      using typename BaseType::DofVectorType;
+
+      static tnlString getTypeStatic();
+
+      tnlString getPrologHeader() const;
+
+      void writeProlog( tnlLogger& logger,
+                        const tnlParameterContainer& parameters ) const;
+
+      bool setup( const tnlParameterContainer& parameters );
+
+      bool setInitialCondition( const tnlParameterContainer& parameters,
+                                const MeshType& mesh,
+                                DofVectorType& dofs,
+                                DofVectorType& auxDofs );
+
+      template< typename MatrixType >
+      bool setupLinearSystem( const MeshType& mesh,
+                              MatrixType& matrix );
+
+      bool makeSnapshot( const RealType& time,
+                         const IndexType& step,
+                         const MeshType& mesh,
+                         DofVectorType& dofs,
+                         DofVectorType& auxDofs );
+
+      IndexType getDofs( const MeshType& mesh ) const;
+
+      void bindDofs( const MeshType& mesh,
+                     DofVectorType& dofs );
+
+      void getExplicitRHS( const RealType& time,
+                           const RealType& tau,
+                           const MeshType& mesh,
+                           DofVectorType& _u,
+                           DofVectorType& _fu );
+
+      template< typename MatrixType >
+      void assemblyLinearSystem( const RealType& time,
+                                 const RealType& tau,
+                                 const MeshType& mesh,
+                                 DofVectorType& dofs,
+                                 DofVectorType& auxDofs,
+                                 MatrixType& matrix,
+                                 DofVectorType& rightHandSide );
+
+
+      protected:
+
+      tnlSharedVector< RealType, DeviceType, IndexType > solution;
+
+      DifferentialOperator differentialOperator;
+
+      BoundaryCondition boundaryCondition;
+   
+      RightHandSide rightHandSide;
+};
+
+#include "tnlMeanCurvatureFlowProblem_impl.h"
+
+#endif /* TNLMEANCURVATUREFLOWPROBLEM_H_ */
diff --git a/examples/mean-curvature-flow/tnlMeanCurvatureFlowProblem_impl.h b/examples/mean-curvature-flow/tnlMeanCurvatureFlowProblem_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..d03ff415046d9cf37ed317891ffa9f0528ecd88f
--- /dev/null
+++ b/examples/mean-curvature-flow/tnlMeanCurvatureFlowProblem_impl.h
@@ -0,0 +1,247 @@
+/***************************************************************************
+                          tnlHeatEquationProblem_impl.h  -  description
+                             -------------------
+    begin                : Mar 10, 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 TNLMEANCURVATUREFLOWPROBLEM_IMPL_H_
+#define TNLMEANCURVATUREFLOWPROBLEM_IMPL_H_
+
+#include <core/mfilename.h>
+#include <matrices/tnlMatrixSetter.h>
+#include <matrices/tnlMultidiagonalMatrixSetter.h>
+#include <core/tnlLogger.h>
+#include <solvers/pde/tnlExplicitUpdater.h>
+#include <solvers/pde/tnlLinearSystemAssembler.h>
+
+#include "tnlMeanCurvatureFlowProblem.h"
+
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+tnlString
+tnlMeanCurvatureFlowProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+getTypeStatic()
+{
+   return tnlString( "tnlMeanCurvativeFlowProblem< " ) + Mesh :: getTypeStatic() + " >";
+}
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+tnlString
+tnlMeanCurvatureFlowProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+getPrologHeader() const
+{
+   return tnlString( "Mean Curvative Flow" );
+}
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+void
+tnlMeanCurvatureFlowProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+writeProlog( tnlLogger& logger, const tnlParameterContainer& parameters ) const
+{
+}
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+bool
+tnlMeanCurvatureFlowProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+setup( const tnlParameterContainer& parameters )
+{
+   if( ! this->boundaryCondition.setup( parameters, "boundary-conditions-" ) ||
+       ! this->rightHandSide.setup( parameters, "right-hand-side-" ) )
+      return false;
+   return true;
+}
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+typename tnlMeanCurvatureFlowProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::IndexType
+tnlMeanCurvatureFlowProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+getDofs( const MeshType& mesh ) const
+{
+   /****
+    * Set-up DOFs and supporting grid functions
+    */
+   return mesh.getNumberOfCells();
+}
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+void
+tnlMeanCurvatureFlowProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+bindDofs( const MeshType& mesh,
+          DofVectorType& dofVector )
+{
+   const IndexType dofs = mesh.getNumberOfCells();
+   this->solution.bind( dofVector.getData(), dofs );
+//   this->differentialOperator.setupDofs(mesh);
+}
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+bool
+tnlMeanCurvatureFlowProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+setInitialCondition( const tnlParameterContainer& parameters,
+                     const MeshType& mesh,
+                     DofVectorType& dofs,
+                     DofVectorType& auxiliaryDofs )
+{
+   this->bindDofs( mesh, dofs );
+   const tnlString& initialConditionFile = parameters.getParameter< tnlString >( "initial-condition" );
+   if( ! this->solution.load( initialConditionFile ) )
+   {
+      cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." << endl;
+      return false;
+   }
+   return true;
+}
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+  template< typename MatrixType >          
+bool
+tnlMeanCurvatureFlowProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+setupLinearSystem( const MeshType& mesh,
+                   MatrixType& matrix )
+{
+//   const IndexType dofs = this->getDofs( mesh );
+//   typedef typename MatrixType::RowLengthsVector RowLengthsVectorType;
+//   RowLengthsVectorType rowLengths;
+//   if( ! rowLengths.setSize( dofs ) )
+//      return false;
+//   tnlMatrixSetter< MeshType, DifferentialOperator, BoundaryCondition, RowLengthsVectorType > matrixSetter;
+//   matrixSetter.template getRowLengths< Mesh::Dimensions >( mesh,
+//                                                            differentialOperator,
+//                                                            boundaryCondition,
+//                                                            rowLengths );
+//   matrix.setDimensions( dofs, dofs );
+//   if( ! matrix.setRowLengths( rowLengths ) )
+//      return false;
+//   return true;
+//   //return tnlMultidiagonalMatrixSetter< Mesh >::setupMatrix( mesh, matrix );
+}
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+bool
+tnlMeanCurvatureFlowProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+makeSnapshot( const RealType& time,
+              const IndexType& step,
+              const MeshType& mesh,
+              DofVectorType& dofs,
+              DofVectorType& auxiliaryDofs )
+{
+   cout << endl << "Writing output at time " << time << " step " << step << "." << endl;
+
+   this->bindDofs( mesh, dofs );
+   //cout << "dofs = " << dofs << endl;
+   tnlString fileName;
+   FileNameBaseNumberEnding( "u-", step, 5, ".tnl", fileName );
+   if( ! this->solution.save( fileName ) )
+      return false;
+   return true;
+}
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+void
+tnlMeanCurvatureFlowProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+getExplicitRHS( const RealType& time,
+                const RealType& tau,
+                const MeshType& mesh,
+                DofVectorType& u,
+                DofVectorType& fu )
+{
+   /****
+    * If you use an explicit solver like tnlEulerSolver or tnlMersonSolver, you
+    * need to implement this method. Compute the right-hand side of
+    *
+    *   d/dt u(x) = fu( x, u )
+    *
+    * You may use supporting vectors again if you need.
+    */
+
+//   this->differentialOperator.computeFirstGradient(mesh,time,u);
+   
+   //cout << "u = " << u << endl;
+   this->bindDofs( mesh, u );
+   tnlExplicitUpdater< Mesh, DofVectorType, DifferentialOperator, BoundaryCondition, RightHandSide > explicitUpdater;
+   explicitUpdater.template update< Mesh::Dimensions >( time,
+                                                        mesh,
+                                                        this->differentialOperator,
+                                                        this->boundaryCondition,
+                                                        this->rightHandSide,
+                                                        u,
+                                                        fu );
+   //cout << "u = " << u << endl;
+   //cout << "fu = " << fu << endl;
+   //_u.save( "u.tnl" );
+   //_fu.save( "fu.tnl" );
+   //getchar();
+}
+
+template< typename Mesh,
+          typename BoundaryCondition,
+          typename RightHandSide,
+          typename DifferentialOperator >
+    template< typename MatrixType >          
+void
+tnlMeanCurvatureFlowProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+assemblyLinearSystem( const RealType& time,
+                      const RealType& tau,
+                      const MeshType& mesh,
+                      DofVectorType& u,
+                      DofVectorType& auxDofs,
+                      MatrixType& matrix,
+                      DofVectorType& b )
+{
+//   tnlLinearSystemAssembler< Mesh, DofVectorType, DifferentialOperator, BoundaryCondition, RightHandSide, MatrixType > systemAssembler;
+//   systemAssembler.template assembly< Mesh::Dimensions >( time,
+//                                                          tau,
+//                                                          mesh,
+//                                                          this->differentialOperator,
+//                                                          this->boundaryCondition,
+//                                                          this->rightHandSide,
+//                                                          u,
+//                                                          matrix,
+//                                                          b );
+//   /*matrix.print( cout );
+//   cout << endl << b << endl;
+//   cout << endl << u << endl;
+//   abort();*/
+}
+
+#endif /* TNLMEANCURVATUREFLOWPROBLEM_IMPL_H_ */
diff --git a/examples/mean-curvature-flow/tnlOperatorQ.h b/examples/mean-curvature-flow/tnlOperatorQ.h
new file mode 100644
index 0000000000000000000000000000000000000000..5887515c6076c013c0fa7c35ff76c124ba938b5e
--- /dev/null
+++ b/examples/mean-curvature-flow/tnlOperatorQ.h
@@ -0,0 +1,135 @@
+#ifndef TNLOPERATORQ_H
+#define	TNLOPERATORQ_H
+
+#include <core/vectors/tnlVector.h>
+#include <mesh/tnlGrid.h>
+#include "tnlBackwardFiniteDifference.h"
+#include "tnlForwardFiniteDifference.h"
+
+template< typename Mesh,
+          typename Real = typename Mesh::RealType,
+          typename Index = typename Mesh::IndexType >
+class tnlOperatorQ
+{
+ 
+};
+
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+class tnlOperatorQ< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index >
+{
+   public: 
+   
+   typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;
+   typedef typename MeshType::CoordinatesType CoordinatesType;
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+
+   static tnlString getType();
+
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif   
+   Real getValue( const MeshType& mesh,
+                   const IndexType cellIndex,
+                   const Vector& u) const 
+   {}
+   
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif   
+   Real getValueStriped( const MeshType& mesh,
+                   const IndexType cellIndex,
+                   const Vector& u) const 
+   {}
+};
+
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+class tnlOperatorQ< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >
+{
+   public: 
+   
+   typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;
+   typedef typename MeshType::CoordinatesType CoordinatesType;
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+
+   static tnlString getType(); 
+
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif   
+   Real getValue( const MeshType& mesh,
+                   const IndexType cellIndex,
+                   const Vector& u) const;
+   
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif   
+   Real getValueStriped( const MeshType& mesh,
+                   const IndexType cellIndex,
+                   const Vector& u) const;
+   
+   private:
+   
+   tnlForwardFiniteDifference<MeshType> fDifference;
+   tnlBackwardFiniteDifference<MeshType> bDifference;
+};
+
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+class tnlOperatorQ< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >
+{
+   public: 
+   
+   typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;
+   typedef typename MeshType::CoordinatesType CoordinatesType;
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+
+   static tnlString getType();
+
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif   
+   Real getValue( const MeshType& mesh,
+                   const IndexType cellIndex,
+                   const Vector& u) const 
+   {}
+   
+   template< typename Vector >
+#ifdef HAVE_CUDA
+   __device__ __host__
+#endif   
+   Real getValueStriped( const MeshType& mesh,
+                   const IndexType cellIndex,
+                   const Vector& u) const 
+   {}
+};
+
+
+#include "tnlOperatorQ_impl.h"
+
+
+#endif	/* TNLOPERATORQ_H */
diff --git a/examples/mean-curvature-flow/tnlOperatorQ_impl.h b/examples/mean-curvature-flow/tnlOperatorQ_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..dd60e10a32f1c3a6be6e5368f45c0b0e6c31e256
--- /dev/null
+++ b/examples/mean-curvature-flow/tnlOperatorQ_impl.h
@@ -0,0 +1,93 @@
+
+#ifndef TTNLOPERATORQ_IMPL_H
+#define	TTNLOPERATORQ_IMPL_H
+
+#include "tnlOperatorQ.h"
+#include <mesh/tnlGrid.h>
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+tnlString
+tnlOperatorQ< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+   return tnlString( "tnlOperatorQ< " ) +
+          MeshType::getType() + ", " +
+          ::getType< Real >() + ", " +
+          ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+tnlString
+tnlOperatorQ< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+   return tnlString( "tnlOperatorQ< " ) +
+          MeshType::getType() + ", " +
+          ::getType< Real >() + ", " +
+          ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlOperatorQ< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getValue( const MeshType& mesh,
+          const IndexType cellIndex,
+          const Vector& u) const
+{
+   return sqrt( 1.0 + (fDifference.getValueX(mesh,cellIndex,u))*(fDifference.getValueX(mesh,cellIndex,u))
+          + (fDifference.getValueY(mesh,cellIndex,u))*(fDifference.getValueY(mesh,cellIndex,u)) );
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlOperatorQ< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getValueStriped( const MeshType& mesh,
+          const IndexType cellIndex,
+          const Vector& u) const
+{
+   return sqrt( 1.0 + 0.5*( (fDifference.getValueX(mesh,cellIndex,u))*(fDifference.getValueX(mesh,cellIndex,u))
+          + (fDifference.getValueY(mesh,cellIndex,u))*(fDifference.getValueY(mesh,cellIndex,u)) 
+          + (bDifference.getValueX(mesh,cellIndex,u))*(bDifference.getValueX(mesh,cellIndex,u))
+          + (bDifference.getValueY(mesh,cellIndex,u))*(bDifference.getValueY(mesh,cellIndex,u)) ) );
+}
+
+template< typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Real,
+          typename Index >
+tnlString
+tnlOperatorQ< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+   return tnlString( "tnlOperatorQ< " ) +
+          MeshType::getType() + ", " +
+          ::getType< Real >() + ", " +
+          ::getType< Index >() + " >";
+}
+
+#endif	/* TTNLOPERATORQ_IMPL_H */
diff --git a/examples/mean-curvature/CMakeLists.txt b/examples/mean-curvature/CMakeLists.txt
deleted file mode 100755
index 79273b47a83aeddd9da7491153445b23b234fcf9..0000000000000000000000000000000000000000
--- a/examples/mean-curvature/CMakeLists.txt
+++ /dev/null
@@ -1,23 +0,0 @@
-set( tnl_heat_equation_SOURCES     
-     tnl-heat-equation.cpp
-     tnl-heat-equation-eoc.cpp )
-               
-IF( BUILD_CUDA )
-   CUDA_ADD_EXECUTABLE(tnl-heat-equation${debugExt} tnl-heat-equation.cu)
-   CUDA_ADD_EXECUTABLE(tnl-heat-equation-eoc-test${debugExt} tnl-heat-equation-eoc.cu)
-ELSE(  BUILD_CUDA )               
-   ADD_EXECUTABLE(tnl-heat-equation${debugExt} tnl-heat-equation.cpp)     
-   ADD_EXECUTABLE(tnl-heat-equation-eoc-test${debugExt} tnl-heat-equation-eoc.cpp)   
-ENDIF( BUILD_CUDA )
-
-target_link_libraries (tnl-heat-equation${debugExt} tnl${debugExt}-${tnlVersion} )
-target_link_libraries (tnl-heat-equation-eoc-test${debugExt} tnl${debugExt}-${tnlVersion} )
-
-INSTALL( TARGETS tnl-heat-equation${debugExt}
-                 tnl-heat-equation-eoc-test${debugExt}
-         RUNTIME DESTINATION bin
-         PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE )
-        
-INSTALL( FILES tnl-run-heat-equation-eoc-test
-               tnl-run-heat-equation
-         DESTINATION share/tnl-${tnlVersion}/examples/heat-equation )