diff --git a/examples/mean-curvature-flow/tnl-mean-curvature-flow-eoc.h b/examples/mean-curvature-flow/tnl-mean-curvature-flow-eoc.h
index af19067c8505e4350008092d46cbfb4f8759dcfb..483bff18cfbff76bd82db2eebfbf8c6cee6f05c7 100644
--- a/examples/mean-curvature-flow/tnl-mean-curvature-flow-eoc.h
+++ b/examples/mean-curvature-flow/tnl-mean-curvature-flow-eoc.h
@@ -19,9 +19,9 @@
 #define TNL_MEAN_CURVATURE_FLOW_EOC_H_
 
 #include <solvers/tnlSolver.h>
-#include <solvers/tnlFastBuildConfig.h>
-#include <solvers/tnlConfigTags.h>
-#include <functions/tnlTestFunction.h>
+#include <solvers/tnlFastBuildConfigTag.h>
+#include <solvers/tnlBuildConfigTags.h>
+#include <functors/tnlTestFunction.h>
 #include <operators/tnlAnalyticDirichletBoundaryConditions.h>
 #include <operators/tnlAnalyticNeumannBoundaryConditions.h>
 #include <problems/tnlMeanCurvatureFlowEocRhs.h>
diff --git a/src/functions/initial_conditions/CMakeLists.txt b/src/functors/initial_conditions/CMakeLists.txt
similarity index 100%
rename from src/functions/initial_conditions/CMakeLists.txt
rename to src/functors/initial_conditions/CMakeLists.txt
diff --git a/src/functions/initial_conditions/CMakeLists.txt~ b/src/functors/initial_conditions/CMakeLists.txt~
similarity index 100%
rename from src/functions/initial_conditions/CMakeLists.txt~
rename to src/functors/initial_conditions/CMakeLists.txt~
diff --git a/src/functions/initial_conditions/level_set_functions/CMakeLists.txt b/src/functors/initial_conditions/level_set_functions/CMakeLists.txt
similarity index 100%
rename from src/functions/initial_conditions/level_set_functions/CMakeLists.txt
rename to src/functors/initial_conditions/level_set_functions/CMakeLists.txt
diff --git a/src/functions/initial_conditions/level_set_functions/CMakeLists.txt~ b/src/functors/initial_conditions/level_set_functions/CMakeLists.txt~
similarity index 100%
rename from src/functions/initial_conditions/level_set_functions/CMakeLists.txt~
rename to src/functors/initial_conditions/level_set_functions/CMakeLists.txt~
diff --git a/src/functions/initial_conditions/level_set_functions/tnlBlobFunction.h b/src/functors/initial_conditions/level_set_functions/tnlBlobFunction.h
similarity index 97%
rename from src/functions/initial_conditions/level_set_functions/tnlBlobFunction.h
rename to src/functors/initial_conditions/level_set_functions/tnlBlobFunction.h
index fe7df502d494fe443470426bae927b42932072c2..1fde719d23b71f6178e77fe0f638ab2838c74eae 100644
--- a/src/functions/initial_conditions/level_set_functions/tnlBlobFunction.h
+++ b/src/functors/initial_conditions/level_set_functions/tnlBlobFunction.h
@@ -20,7 +20,7 @@
 
 #include <config/tnlParameterContainer.h>
 #include <core/vectors/tnlStaticVector.h>
-#include <functions/tnlFunctionType.h>
+#include <functors/tnlFunctionType.h>
 
 template< typename Real >
 class tnlBlobFunctionBase
@@ -154,7 +154,7 @@ class tnlFunctionType< tnlBlobFunction< FunctionDimensions, Real > >
 };
 
 
-#include <functions/initial_conditions/level_set_functions/tnlBlobFunction_impl.h>
+#include <functors/initial_conditions/level_set_functions/tnlBlobFunction_impl.h>
 
 
 #endif /* TNLBLOBFUNCTION_H_ */
diff --git a/src/functions/initial_conditions/level_set_functions/tnlBlobFunction.h~ b/src/functors/initial_conditions/level_set_functions/tnlBlobFunction.h~
similarity index 95%
rename from src/functions/initial_conditions/level_set_functions/tnlBlobFunction.h~
rename to src/functors/initial_conditions/level_set_functions/tnlBlobFunction.h~
index 0790387e279b85efce31161e311a44baf62b9946..c085289494fb4d25b5e0abe78e707d266f833be2 100644
--- a/src/functions/initial_conditions/level_set_functions/tnlBlobFunction.h~
+++ b/src/functors/initial_conditions/level_set_functions/tnlBlobFunction.h~
@@ -20,7 +20,7 @@
 
 #include <config/tnlParameterContainer.h>
 #include <core/vectors/tnlStaticVector.h>
-#include <functions/tnlFunctionType.h>
+#include <functors/tnlFunctionType.h>
 
 template< typename Real >
 class tnlBlobFunctionBase
@@ -72,10 +72,6 @@ class tnlBlobFunction< 1, Real > : public tnlBlobFunctionBase< Real >
 #endif
       RealType getValue( const Vertex& v,
                          const Real& time = 0.0 ) const;
-                         
-   protected:
-
-      RealType height;
 };
 
 template< typename Real >
@@ -107,10 +103,6 @@ class tnlBlobFunction< 2, Real > : public tnlBlobFunctionBase< Real >
 #endif
       RealType getValue( const Vertex& v,
                          const Real& time = 0.0 ) const;
-                         
-   protected:
-
-      RealType height;
 };
 
 template< typename Real >
@@ -142,10 +134,6 @@ class tnlBlobFunction< 3, Real > : public tnlBlobFunctionBase< Real >
 #endif
       RealType getValue( const Vertex& v,
                          const Real& time = 0.0 ) const;
-
-   protected:
-
-      RealType height;
 };
 
 template< int Dimensions,
diff --git a/src/functions/initial_conditions/level_set_functions/tnlBlobFunction_impl.h b/src/functors/initial_conditions/level_set_functions/tnlBlobFunction_impl.h
similarity index 97%
rename from src/functions/initial_conditions/level_set_functions/tnlBlobFunction_impl.h
rename to src/functors/initial_conditions/level_set_functions/tnlBlobFunction_impl.h
index 6e05320396a2f903146eb82e44730ca8c2a349b7..284acbddb87732dbd7848c27b1e59db4c421477b 100644
--- a/src/functions/initial_conditions/level_set_functions/tnlBlobFunction_impl.h
+++ b/src/functors/initial_conditions/level_set_functions/tnlBlobFunction_impl.h
@@ -18,7 +18,7 @@
 #ifndef TNLBLOBFUNCTION_IMPL_H_
 #define TNLBLOBFUNCTION_IMPL_H_
 
-#include <functions/initial_conditions/level_set_functions/tnlBlobFunction.h>
+#include <functors/initial_conditions/level_set_functions/tnlBlobFunction.h>
 
 template< typename Real >
 bool
diff --git a/src/functions/initial_conditions/level_set_functions/tnlBlobFunction_impl.h~ b/src/functors/initial_conditions/level_set_functions/tnlBlobFunction_impl.h~
similarity index 100%
rename from src/functions/initial_conditions/level_set_functions/tnlBlobFunction_impl.h~
rename to src/functors/initial_conditions/level_set_functions/tnlBlobFunction_impl.h~
diff --git a/src/functions/initial_conditions/level_set_functions/tnlPseudoSquareFunction.h b/src/functors/initial_conditions/level_set_functions/tnlPseudoSquareFunction.h
similarity index 97%
rename from src/functions/initial_conditions/level_set_functions/tnlPseudoSquareFunction.h
rename to src/functors/initial_conditions/level_set_functions/tnlPseudoSquareFunction.h
index 2f343be12c3a6ac6514e9bdd86df4327f7d616b7..af48877c1ab6148e7989c51738329eb5b847285c 100644
--- a/src/functions/initial_conditions/level_set_functions/tnlPseudoSquareFunction.h
+++ b/src/functors/initial_conditions/level_set_functions/tnlPseudoSquareFunction.h
@@ -20,7 +20,7 @@
 
 #include <config/tnlParameterContainer.h>
 #include <core/vectors/tnlStaticVector.h>
-#include <functions/tnlFunctionType.h>
+#include <functors/tnlFunctionType.h>
 
 template< typename Real >
 class tnlPseudoSquareFunctionBase
@@ -154,7 +154,7 @@ class tnlFunctionType< tnlPseudoSquareFunction< FunctionDimensions, Real > >
 };
 
 
-#include <functions/initial_conditions/level_set_functions/tnlPseudoSquareFunction_impl.h>
+#include <functors/initial_conditions/level_set_functions/tnlPseudoSquareFunction_impl.h>
 
 
 #endif /* TNLPSEUDOSQUAREFUNCTION_H_ */
diff --git a/src/functions/initial_conditions/level_set_functions/tnlPseudoSquareFunction.h~ b/src/functors/initial_conditions/level_set_functions/tnlPseudoSquareFunction.h~
similarity index 95%
rename from src/functions/initial_conditions/level_set_functions/tnlPseudoSquareFunction.h~
rename to src/functors/initial_conditions/level_set_functions/tnlPseudoSquareFunction.h~
index 59a8e2d90dc4aff46cf750efaf15333cfc6011e7..5c2350e9bbbbd2cd01834d99ffc4589b85b21913 100644
--- a/src/functions/initial_conditions/level_set_functions/tnlPseudoSquareFunction.h~
+++ b/src/functors/initial_conditions/level_set_functions/tnlPseudoSquareFunction.h~
@@ -20,7 +20,7 @@
 
 #include <config/tnlParameterContainer.h>
 #include <core/vectors/tnlStaticVector.h>
-#include <functions/tnlFunctionType.h>
+#include <functors/tnlFunctionType.h>
 
 template< typename Real >
 class tnlPseudoSquareFunctionBase
@@ -72,10 +72,6 @@ class tnlPseudoSquareFunction< 1, Real > : public tnlPseudoSquareFunctionBase< R
 #endif
       RealType getValue( const Vertex& v,
                          const Real& time = 0.0 ) const;
-                         
-   protected:
-
-      RealType height;
 };
 
 template< typename Real >
@@ -107,10 +103,6 @@ class tnlPseudoSquareFunction< 2, Real > : public tnlPseudoSquareFunctionBase< R
 #endif
       RealType getValue( const Vertex& v,
                          const Real& time = 0.0 ) const;
-                         
-   protected:
-
-      RealType height;
 };
 
 template< typename Real >
@@ -142,10 +134,6 @@ class tnlPseudoSquareFunction< 3, Real > : public tnlPseudoSquareFunctionBase< R
 #endif
       RealType getValue( const Vertex& v,
                          const Real& time = 0.0 ) const;
-
-   protected:
-
-      RealType height;
 };
 
 template< int Dimensions,
diff --git a/src/functions/initial_conditions/level_set_functions/tnlPseudoSquareFunction_impl.h~ b/src/functors/initial_conditions/level_set_functions/tnlPseudoSquareFunction_impl.h
similarity index 96%
rename from src/functions/initial_conditions/level_set_functions/tnlPseudoSquareFunction_impl.h~
rename to src/functors/initial_conditions/level_set_functions/tnlPseudoSquareFunction_impl.h
index 6fbcae784ed379871873741edeb8a102ed7809f5..2e02db5e5df616e56c8396eb9c3073cb6f7f47df 100644
--- a/src/functions/initial_conditions/level_set_functions/tnlPseudoSquareFunction_impl.h~
+++ b/src/functors/initial_conditions/level_set_functions/tnlPseudoSquareFunction_impl.h
@@ -18,11 +18,11 @@
 #ifndef TNLPSEUDOSQUAREFUNCTION_IMPL_H_
 #define TNLPSEUDOSQUAREFUNCTION_IMPL_H_
 
-#include <functions/initial_conditions/level_set_functions/tnlPseudoSquareFunction.h>
+#include <functors/initial_conditions/level_set_functions/tnlPseudoSquareFunction.h>
 
 template< typename Real >
 bool
-tnlPseudoSquareFunction< Real >::
+tnlPseudoSquareFunctionBase< Real >::
 setup( const tnlParameterContainer& parameters,
        const tnlString& prefix )
 {
diff --git a/src/functions/initial_conditions/level_set_functions/tnlPseudoSquareFunction_impl.h b/src/functors/initial_conditions/level_set_functions/tnlPseudoSquareFunction_impl.h~
similarity index 100%
rename from src/functions/initial_conditions/level_set_functions/tnlPseudoSquareFunction_impl.h
rename to src/functors/initial_conditions/level_set_functions/tnlPseudoSquareFunction_impl.h~
diff --git a/src/functions/initial_conditions/tnlCylinderFunction.h b/src/functors/initial_conditions/tnlCylinderFunction.h
similarity index 97%
rename from src/functions/initial_conditions/tnlCylinderFunction.h
rename to src/functors/initial_conditions/tnlCylinderFunction.h
index 4d28b4c3f6c2538e4a9deee073d7e94baeefc71b..dd5d41064624c3a2863d84051aae1f2927a7d495 100644
--- a/src/functions/initial_conditions/tnlCylinderFunction.h
+++ b/src/functors/initial_conditions/tnlCylinderFunction.h
@@ -20,7 +20,7 @@
 
 #include <config/tnlParameterContainer.h>
 #include <core/vectors/tnlStaticVector.h>
-#include <functions/tnlFunctionType.h>
+#include <functors/tnlFunctionType.h>
 
 template< typename Real >
 class tnlCylinderFunctionBase
@@ -158,7 +158,7 @@ class tnlFunctionType< tnlCylinderFunction< FunctionDimensions, Real > >
 };
 
 
-#include <functions/initial_conditions/tnlCylinderFunction_impl.h>
+#include <functors/initial_conditions/tnlCylinderFunction_impl.h>
 
 
 #endif /* TNLEXPBUMPFUNCTION_H_ */
diff --git a/src/functions/initial_conditions/tnlCylinderFunction.h~ b/src/functors/initial_conditions/tnlCylinderFunction.h~
similarity index 94%
rename from src/functions/initial_conditions/tnlCylinderFunction.h~
rename to src/functors/initial_conditions/tnlCylinderFunction.h~
index 440e53495a5e2e77ef5790cabdb369cc7a2d8d33..ce90d3f447ebc63929ca07dbdd0c03193d159a8d 100644
--- a/src/functions/initial_conditions/tnlCylinderFunction.h~
+++ b/src/functors/initial_conditions/tnlCylinderFunction.h~
@@ -20,7 +20,7 @@
 
 #include <config/tnlParameterContainer.h>
 #include <core/vectors/tnlStaticVector.h>
-#include <functions/tnlFunctionType.h>
+#include <functors/tnlFunctionType.h>
 
 template< typename Real >
 class tnlCylinderFunctionBase
@@ -38,7 +38,7 @@ class tnlCylinderFunctionBase
 
    protected:
 
-      RealType amplitude, sigma;
+      RealType diameter;
 };
 
 template< int Dimensions,
@@ -142,9 +142,9 @@ class tnlCylinderFunction< 3, Real > : public tnlCylinderFunctionBase< Real >
 
 template< int Dimensions,
           typename Real >
-ostream& operator << ( ostream& str, const tnlExpBumpFunction< Dimensions, Real >& f )
+ostream& operator << ( ostream& str, const tnlCylinderFunction< Dimensions, Real >& f )
 {
-   str << "Cylinder.;
+   str << "Cylinder function.";
    return str;
 }
 
@@ -158,7 +158,7 @@ class tnlFunctionType< tnlCylinderFunction< FunctionDimensions, Real > >
 };
 
 
-#include <functions/tnlCylinderFunction_impl.h>
+#include <functions/initial_conditions/tnlCylinderFunction_impl.h>
 
 
 #endif /* TNLEXPBUMPFUNCTION_H_ */
diff --git a/src/functions/initial_conditions/tnlCylinderFunction_impl.h~ b/src/functors/initial_conditions/tnlCylinderFunction_impl.h
similarity index 61%
rename from src/functions/initial_conditions/tnlCylinderFunction_impl.h~
rename to src/functors/initial_conditions/tnlCylinderFunction_impl.h
index a2335af6473c86b41162ea5658ddfaec9e6e2d35..c5c35359d044ecc19dd79747a00de2f59fe277bb 100644
--- a/src/functions/initial_conditions/tnlCylinderFunction_impl.h~
+++ b/src/functors/initial_conditions/tnlCylinderFunction_impl.h
@@ -18,7 +18,7 @@
 #ifndef TNLCYLINDERFUNCTION_IMPL_H_
 #define TNLCYLINDERFUNCTION_IMPL_H_
 
-#include <functions/tnlCylinderFunction.h>
+#include <functors/initial_conditions/tnlCylinderFunction.h>
 
 template< typename Real >
 bool
@@ -149,25 +149,7 @@ getValue( const Vertex& v,
    const RealType& y = v.y();
    const RealType& z = v.z();
    if( XDiffOrder == 0 && YDiffOrder == 0 && ZDiffOrder == 0 )
-      return this->amplitude * exp( ( -x*x - y*y -z*z ) / ( this->sigma * this->sigma ) );
-   if( XDiffOrder == 1 && YDiffOrder == 0 && ZDiffOrder == 0 )
-      return -2.0 * x / ( this->sigma * this->sigma ) * this->amplitude * exp( ( -x*x - y*y -z*z ) / ( this->sigma * this->sigma ) );
-   if( XDiffOrder == 2 && YDiffOrder == 0 && ZDiffOrder == 0 )
-      return -2.0 / ( this->sigma * this->sigma ) * this->amplitude * exp( ( -x*x - y*y -z*z ) / ( this->sigma * this->sigma ) ) + 4.0 * x * x / ( this->sigma * this->sigma * this->sigma * this->sigma ) * this->amplitude * exp( ( -x*x - y*y -z*z ) / ( this->sigma * this->sigma ) );
-   if( XDiffOrder == 0 && YDiffOrder == 1 && ZDiffOrder == 0 )
-      return -2.0 * y / ( this->sigma * this->sigma ) * this->amplitude * exp( ( -x*x - y*y -z*z ) / ( this->sigma * this->sigma ) );
-   if( XDiffOrder == 0 && YDiffOrder == 2 && ZDiffOrder == 0 )
-      return -2.0 / ( this->sigma * this->sigma ) * this->amplitude * exp( ( -x*x - y*y -z*z ) / ( this->sigma * this->sigma ) ) + 4.0 * y * y / ( this->sigma * this->sigma * this->sigma * this->sigma ) * this->amplitude * exp( ( -x*x - y*y -z*z ) / ( this->sigma * this->sigma ) );
-   if( XDiffOrder == 0 && YDiffOrder == 0 && ZDiffOrder == 1 )
-      return -2.0 * z / ( this->sigma * this->sigma ) * this->amplitude * exp( ( -x*x - y*y -z*z ) / ( this->sigma * this->sigma ) );
-   if( XDiffOrder == 0 && YDiffOrder == 0 && ZDiffOrder == 2 )
-      return -2.0 / ( this->sigma * this->sigma ) * this->amplitude * exp( ( -x*x - y*y -z*z ) / ( this->sigma * this->sigma ) ) + 4.0 * z * z / ( this->sigma * this->sigma * this->sigma * this->sigma ) * this->amplitude * exp( ( -x*x - y*y -z*z ) / ( this->sigma * this->sigma ) );
-   if( XDiffOrder == 1 && YDiffOrder == 1 && ZDiffOrder == 0 )
-      return 4.0 * x * y / ( ( this->sigma * this->sigma ) * ( this->sigma * this->sigma ) ) * this->amplitude * exp( ( -x*x - y*y -z*z ) / ( this->sigma * this->sigma ) );
-   if( XDiffOrder == 1 && YDiffOrder == 0 && ZDiffOrder == 1 )
-      return 4.0 * x * z / ( ( this->sigma * this->sigma ) * ( this->sigma * this->sigma ) ) * this->amplitude * exp( ( -x*x - y*y -z*z ) / ( this->sigma * this->sigma ) );
-   if( XDiffOrder == 0 && YDiffOrder == 1 && ZDiffOrder == 1 )
-      return 4.0 * y * z / ( ( this->sigma * this->sigma ) * ( this->sigma * this->sigma ) ) * this->amplitude * exp( ( -x*x - y*y -z*z ) / ( this->sigma * this->sigma ) );
+      return ( ( x*x + y*y + z*z - this->diameter ) < 0 ) - ( ( x*x + y*y + z*z - this->diameter ) > 0 ) + 1;
    return 0.0;
 }
 
diff --git a/src/functions/initial_conditions/tnlCylinderFunction_impl.h b/src/functors/initial_conditions/tnlCylinderFunction_impl.h~
similarity index 100%
rename from src/functions/initial_conditions/tnlCylinderFunction_impl.h
rename to src/functors/initial_conditions/tnlCylinderFunction_impl.h~
diff --git a/src/functions/initial_conditions/tnlFlowerpotFunction.h b/src/functors/initial_conditions/tnlFlowerpotFunction.h
similarity index 97%
rename from src/functions/initial_conditions/tnlFlowerpotFunction.h
rename to src/functors/initial_conditions/tnlFlowerpotFunction.h
index 36a5adaf0892035387b40d86f940784aef01b686..56fbb2ee726702c4cc4a5c4598b7cce1c27dd9dc 100644
--- a/src/functions/initial_conditions/tnlFlowerpotFunction.h
+++ b/src/functors/initial_conditions/tnlFlowerpotFunction.h
@@ -20,7 +20,7 @@
 
 #include <config/tnlParameterContainer.h>
 #include <core/vectors/tnlStaticVector.h>
-#include <functions/tnlFunctionType.h>
+#include <functors/tnlFunctionType.h>
 
 template< typename Real >
 class tnlFlowerpotFunctionBase
@@ -158,7 +158,7 @@ class tnlFunctionType< tnlFlowerpotFunction< FunctionDimensions, Real > >
 };
 
 
-#include <functions/initial_conditions/tnlFlowerpotFunction_impl.h>
+#include <functors/initial_conditions/tnlFlowerpotFunction_impl.h>
 
 
 #endif /* TNLFLOWERPOTFUNCTION_H_ */
diff --git a/src/functors/initial_conditions/tnlFlowerpotFunction.h~ b/src/functors/initial_conditions/tnlFlowerpotFunction.h~
new file mode 100644
index 0000000000000000000000000000000000000000..04fc344a6d61112bb15b6b807edc90ed60871de6
--- /dev/null
+++ b/src/functors/initial_conditions/tnlFlowerpotFunction.h~
@@ -0,0 +1,164 @@
+/***************************************************************************
+                          tnlExpBumpFunction.h  -  description
+                             -------------------
+    begin                : Dec 5, 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 TNLFLOWERPOTFUNCTION_H_
+#define TNLFLOWERPOTFUNCTION_H_
+
+#include <config/tnlParameterContainer.h>
+#include <core/vectors/tnlStaticVector.h>
+#include <functors/tnlFunctionType.h>
+
+template< typename Real >
+class tnlFlowerpotFunctionBase
+{
+   public:
+
+      typedef Real RealType;
+
+      bool setup( const tnlParameterContainer& parameters,
+                 const tnlString& prefix = "" );
+
+      void setDiameter( const RealType& sigma );
+
+      const RealType& getDiameter() const;
+
+   protected:
+
+      RealType diameter;
+};
+
+template< int Dimensions,
+          typename Real >
+class tnlFlowerpotFunction
+{
+};
+
+template< typename Real >
+class tnlFlowerpotFunction< 1, Real > : public tnlFlowerpotFunctionBase< Real >
+{
+   public:
+
+      enum { Dimensions = 1 };
+      typedef Real RealType;
+      typedef tnlStaticVector< Dimensions, Real > VertexType;
+
+      static tnlString getType();
+
+      tnlFlowerpotFunction();
+
+#ifdef HAVE_NOT_CXX11
+      template< int XDiffOrder,
+                int YDiffOrder,
+                int ZDiffOrder,
+                typename Vertex >
+#else
+      template< int XDiffOrder = 0,
+                int YDiffOrder = 0,
+                int ZDiffOrder = 0,
+                typename Vertex = VertexType >
+#endif   
+#ifdef HAVE_CUDA
+      __device__ __host__
+#endif
+      RealType getValue( const Vertex& v,
+                         const Real& time = 0.0 ) const;
+};
+
+template< typename Real >
+class tnlFlowerpotFunction< 2, Real > : public tnlFlowerpotFunctionBase< Real >
+{
+   public:
+
+      enum { Dimensions = 2 };
+      typedef Real RealType;
+      typedef tnlStaticVector< Dimensions, Real > VertexType;
+
+      static tnlString getType();
+
+      tnlFlowerpotFunction();
+
+#ifdef HAVE_NOT_CXX11
+      template< int XDiffOrder,
+                int YDiffOrder,
+                int ZDiffOrder,
+                typename Vertex >
+#else
+      template< int XDiffOrder = 0,
+                int YDiffOrder = 0,
+                int ZDiffOrder = 0,
+                typename Vertex = VertexType >
+#endif
+#ifdef HAVE_CUDA
+      __device__ __host__
+#endif
+      RealType getValue( const Vertex& v,
+                         const Real& time = 0.0 ) const;
+};
+
+template< typename Real >
+class tnlFlowerpotFunction< 3, Real > : public tnlFlowerpotFunctionBase< Real >
+{
+   public:
+
+      enum { Dimensions = 3 };
+      typedef Real RealType;
+      typedef tnlStaticVector< Dimensions, Real > VertexType;
+
+      static tnlString getType();
+
+      tnlFlowerpotFunction();
+
+#ifdef HAVE_NOT_CXX11
+      template< int XDiffOrder,
+                int YDiffOrder,
+                int ZDiffOrder,
+                typename Vertex >
+#else
+      template< int XDiffOrder = 0,
+                int YDiffOrder = 0,
+                int ZDiffOrder = 0,
+                typename Vertex = VertexType >
+#endif   
+#ifdef HAVE_CUDA
+      __device__ __host__
+#endif
+      RealType getValue( const Vertex& v,
+                         const Real& time = 0.0 ) const;
+};
+
+template< int Dimensions,
+          typename Real >
+ostream& operator << ( ostream& str, const tnlFlowerpotFunction< Dimensions, Real >& f )
+{
+   str << "Flowerpot function.";
+   return str;
+}
+
+template< int FunctionDimensions,
+          typename Real >
+class tnlFunctionType< tnlFlowerpotFunction< FunctionDimensions, Real > >
+{
+   public:
+
+      enum { Type = tnlAnalyticFunction };
+};
+
+
+#include <functions/initial_conditions/tnlFlowerpotFunction_impl.h>
+
+
+#endif /* TNLFLOWERPOTFUNCTION_H_ */
diff --git a/src/functors/initial_conditions/tnlFlowerpotFunction_impl.h b/src/functors/initial_conditions/tnlFlowerpotFunction_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..ba3ffdc662ac1cda78d31095fbc1dd16e5f97c09
--- /dev/null
+++ b/src/functors/initial_conditions/tnlFlowerpotFunction_impl.h
@@ -0,0 +1,157 @@
+/***************************************************************************
+                          tnlExpBumpFunction_impl.h  -  description
+                             -------------------
+    begin                : Dec 5, 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 TNLFLOWERPOTFUNCTION_IMPL_H_
+#define TNLFLOWERPOTFUNCTION_IMPL_H_
+
+#include <functors/initial_conditions/tnlFlowerpotFunction.h>
+
+template< typename Real >
+bool
+tnlFlowerpotFunctionBase< Real >::
+setup( const tnlParameterContainer& parameters,
+       const tnlString& prefix )
+{
+   this->diameter = parameters.getParameter< double >( prefix + "diameter" );
+   return true;
+}
+
+template< typename Real >
+void tnlFlowerpotFunctionBase< Real >::setDiameter( const Real& sigma )
+{
+   this->diameter = diameter;
+}
+
+template< typename Real >
+const Real& tnlFlowerpotFunctionBase< Real >::getDiameter() const
+{
+   return this->diameter;
+}
+
+/***
+ * 1D
+ */
+
+template< typename Real >
+tnlString
+tnlFlowerpotFunction< 1, Real >::getType()
+{
+   return "tnlFlowerpotFunction< 1, " + ::getType< Real >() + tnlString( " >" );
+}
+
+template< typename Real >
+tnlFlowerpotFunction< 1, Real >::tnlFlowerpotFunction()
+{
+}
+
+template< typename Real >
+   template< int XDiffOrder, 
+             int YDiffOrder,
+             int ZDiffOrder,
+             typename Vertex >
+#ifdef HAVE_CUDA
+      __device__ __host__
+#endif
+Real
+tnlFlowerpotFunction< 1, Real >::getValue( const Vertex& v,
+                                         const Real& time ) const
+{
+   const RealType& x = v.x();
+   if( YDiffOrder != 0 || ZDiffOrder != 0 )
+      return 0.0;
+   if( XDiffOrder == 0 )
+      return sin( M_PI * tanh( 5 * ( x * x - this->diameter ) ) );
+   return 0.0;
+}
+
+/****
+ * 2D
+ */
+
+template< typename Real >
+tnlString
+tnlFlowerpotFunction< 2, Real >::getType()
+{
+   return tnlString( "tnlFlowerpotFunction< 2, " ) + ::getType< Real >() + " >";
+}
+
+template< typename Real >
+tnlFlowerpotFunction< 2, Real >::tnlFlowerpotFunction()
+{
+}
+
+template< typename Real >
+   template< int XDiffOrder,
+             int YDiffOrder,
+             int ZDiffOrder,
+             typename Vertex >
+#ifdef HAVE_CUDA
+      __device__ __host__
+#endif
+Real
+tnlFlowerpotFunction< 2, Real >::
+getValue( const Vertex& v,
+          const Real& time ) const
+{
+   const RealType& x = v.x();
+   const RealType& y = v.y();
+   if( ZDiffOrder != 0 )
+      return 0.0;
+   if( XDiffOrder == 0 && YDiffOrder == 0 )
+      return sin( M_PI * tanh( 5 * ( x * x + y * y - this->diameter ) ) );
+   return 0.0;
+}
+
+/****
+ * 3D
+ */
+
+template< typename Real >
+tnlString
+tnlFlowerpotFunction< 3, Real >::getType()
+{
+   return tnlString( "tnlFlowerpotFunction< 3, " ) + ::getType< Real >() + " >";
+}
+
+template< typename Real >
+tnlFlowerpotFunction< 3, Real >::tnlFlowerpotFunction()
+{
+}
+
+template< typename Real >
+   template< int XDiffOrder,
+             int YDiffOrder,
+             int ZDiffOrder,
+             typename Vertex >
+#ifdef HAVE_CUDA
+      __device__ __host__
+#endif
+Real
+tnlFlowerpotFunction< 3, Real >::
+getValue( const Vertex& v,
+          const Real& time ) const
+{
+   const RealType& x = v.x();
+   const RealType& y = v.y();
+   const RealType& z = v.z();
+   if( XDiffOrder == 0 && YDiffOrder == 0 && ZDiffOrder == 0 )
+      return sin( M_PI * tanh( 5 * ( x * x + y * y + z * z - 0.25 ) ) );
+   return 0.0;
+}
+
+
+#endif /* TNLFLOWERPOTFUNCTION_IMPL_H_ */
diff --git a/src/functions/initial_conditions/tnlFlowerpotFunction_impl.h b/src/functors/initial_conditions/tnlFlowerpotFunction_impl.h~
similarity index 100%
rename from src/functions/initial_conditions/tnlFlowerpotFunction_impl.h
rename to src/functors/initial_conditions/tnlFlowerpotFunction_impl.h~
diff --git a/src/functions/initial_conditions/tnlTwinsFunction.h b/src/functors/initial_conditions/tnlTwinsFunction.h
similarity index 97%
rename from src/functions/initial_conditions/tnlTwinsFunction.h
rename to src/functors/initial_conditions/tnlTwinsFunction.h
index f65d4ad2c1f7879e1293d7149d3abf028380129a..485c429e4eef0acdcbf980d613b0de68eece4a2e 100644
--- a/src/functions/initial_conditions/tnlTwinsFunction.h
+++ b/src/functors/initial_conditions/tnlTwinsFunction.h
@@ -20,7 +20,7 @@
 
 #include <config/tnlParameterContainer.h>
 #include <core/vectors/tnlStaticVector.h>
-#include <functions/tnlFunctionType.h>
+#include <functors/tnlFunctionType.h>
 
 template< typename Real >
 class tnlTwinsFunctionBase
@@ -150,7 +150,7 @@ class tnlFunctionType< tnlTwinsFunction< FunctionDimensions, Real > >
 };
 
 
-#include <functions/initial_conditions/tnlTwinsFunction_impl.h>
+#include <functors/initial_conditions/tnlTwinsFunction_impl.h>
 
 
 #endif /* TNLTWINSFUNCTION_H_ */
diff --git a/src/functors/initial_conditions/tnlTwinsFunction.h~ b/src/functors/initial_conditions/tnlTwinsFunction.h~
new file mode 100644
index 0000000000000000000000000000000000000000..e91c1f0e998689a657417f96d7907cb50d7c9716
--- /dev/null
+++ b/src/functors/initial_conditions/tnlTwinsFunction.h~
@@ -0,0 +1,156 @@
+/***************************************************************************
+                          tnlExpBumpFunction.h  -  description
+                             -------------------
+    begin                : Dec 5, 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 TNLTWINSFUNCTION_H_
+#define TNLTWINSFUNCTION_H_
+
+#include <config/tnlParameterContainer.h>
+#include <core/vectors/tnlStaticVector.h>
+#include <functors/tnlFunctionType.h>
+
+template< typename Real >
+class tnlTwinsFunctionBase
+{
+   public:
+
+      typedef Real RealType;
+
+      bool setup( const tnlParameterContainer& parameters,
+                 const tnlString& prefix = "" );
+};
+
+template< int Dimensions,
+          typename Real >
+class tnlTwinsFunction
+{
+};
+
+template< typename Real >
+class tnlTwinsFunction< 1, Real > : public tnlTwinsFunctionBase< Real >
+{
+   public:
+
+      enum { Dimensions = 1 };
+      typedef Real RealType;
+      typedef tnlStaticVector< Dimensions, Real > VertexType;
+
+      static tnlString getType();
+
+      tnlTwinsFunction();
+
+#ifdef HAVE_NOT_CXX11
+      template< int XDiffOrder,
+                int YDiffOrder,
+                int ZDiffOrder,
+                typename Vertex >
+#else
+      template< int XDiffOrder = 0,
+                int YDiffOrder = 0,
+                int ZDiffOrder = 0,
+                typename Vertex = VertexType >
+#endif   
+#ifdef HAVE_CUDA
+      __device__ __host__
+#endif
+      RealType getValue( const Vertex& v,
+                         const Real& time = 0.0 ) const;
+};
+
+template< typename Real >
+class tnlTwinsFunction< 2, Real > : public tnlTwinsFunctionBase< Real >
+{
+   public:
+
+      enum { Dimensions = 2 };
+      typedef Real RealType;
+      typedef tnlStaticVector< Dimensions, Real > VertexType;
+
+      static tnlString getType();
+
+      tnlTwinsFunction();
+
+#ifdef HAVE_NOT_CXX11
+      template< int XDiffOrder,
+                int YDiffOrder,
+                int ZDiffOrder,
+                typename Vertex >
+#else
+      template< int XDiffOrder = 0,
+                int YDiffOrder = 0,
+                int ZDiffOrder = 0,
+                typename Vertex = VertexType >
+#endif
+#ifdef HAVE_CUDA
+      __device__ __host__
+#endif
+      RealType getValue( const Vertex& v,
+                         const Real& time = 0.0 ) const;
+};
+
+template< typename Real >
+class tnlTwinsFunction< 3, Real > : public tnlTwinsFunctionBase< Real >
+{
+   public:
+
+      enum { Dimensions = 3 };
+      typedef Real RealType;
+      typedef tnlStaticVector< Dimensions, Real > VertexType;
+
+      static tnlString getType();
+
+      tnlTwinsFunction();
+
+#ifdef HAVE_NOT_CXX11
+      template< int XDiffOrder,
+                int YDiffOrder,
+                int ZDiffOrder,
+                typename Vertex >
+#else
+      template< int XDiffOrder = 0,
+                int YDiffOrder = 0,
+                int ZDiffOrder = 0,
+                typename Vertex = VertexType >
+#endif   
+#ifdef HAVE_CUDA
+      __device__ __host__
+#endif
+      RealType getValue( const Vertex& v,
+                         const Real& time = 0.0 ) const;
+};
+
+template< int Dimensions,
+          typename Real >
+ostream& operator << ( ostream& str, const tnlTwinsFunction< Dimensions, Real >& f )
+{
+   str << "Twins function.";
+   return str;
+}
+
+template< int FunctionDimensions,
+          typename Real >
+class tnlFunctionType< tnlTwinsFunction< FunctionDimensions, Real > >
+{
+   public:
+
+      enum { Type = tnlAnalyticFunction };
+};
+
+
+#include <functions/initial_conditions/tnlTwinsFunction_impl.h>
+
+
+#endif /* TNLTWINSFUNCTION_H_ */
diff --git a/src/functors/initial_conditions/tnlTwinsFunction_impl.h b/src/functors/initial_conditions/tnlTwinsFunction_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..910cab6ccba6265a66679613549512bf87a2f7a7
--- /dev/null
+++ b/src/functors/initial_conditions/tnlTwinsFunction_impl.h
@@ -0,0 +1,145 @@
+/***************************************************************************
+                          tnlExpBumpFunction_impl.h  -  description
+                             -------------------
+    begin                : Dec 5, 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 TNLTWINSFUNCTION_IMPL_H_
+#define TNLTWINSFUNCTION_IMPL_H_
+
+#include <functors/initial_conditions/tnlTwinsFunction.h>
+
+template< typename Real >
+bool
+tnlTwinsFunctionBase< Real >::
+setup( const tnlParameterContainer& parameters,
+       const tnlString& prefix )
+{
+   return true;
+}
+
+
+/***
+ * 1D
+ */
+
+template< typename Real >
+tnlString
+tnlTwinsFunction< 1, Real >::getType()
+{
+   return "tnlTwinsFunction< 1, " + ::getType< Real >() + tnlString( " >" );
+}
+
+template< typename Real >
+tnlTwinsFunction< 1, Real >::tnlTwinsFunction()
+{
+}
+
+template< typename Real >
+   template< int XDiffOrder, 
+             int YDiffOrder,
+             int ZDiffOrder,
+             typename Vertex >
+#ifdef HAVE_CUDA
+      __device__ __host__
+#endif
+Real
+tnlTwinsFunction< 1, Real >::getValue( const Vertex& v,
+                                         const Real& time ) const
+{
+   const RealType& x = v.x();
+   if( YDiffOrder != 0 || ZDiffOrder != 0 )
+      return 0.0;
+   if( XDiffOrder == 0 )
+      return 0.0;
+   return 0.0;
+}
+
+/****
+ * 2D
+ */
+
+template< typename Real >
+tnlString
+tnlTwinsFunction< 2, Real >::getType()
+{
+   return tnlString( "tnlTwinsFunction< 2, " ) + ::getType< Real >() + " >";
+}
+
+template< typename Real >
+tnlTwinsFunction< 2, Real >::tnlTwinsFunction()
+{
+}
+
+template< typename Real >
+   template< int XDiffOrder,
+             int YDiffOrder,
+             int ZDiffOrder,
+             typename Vertex >
+#ifdef HAVE_CUDA
+      __device__ __host__
+#endif
+Real
+tnlTwinsFunction< 2, Real >::
+getValue( const Vertex& v,
+          const Real& time ) const
+{
+   const RealType& x = v.x();
+   const RealType& y = v.y();
+   if( ZDiffOrder != 0 )
+      return 0.0;
+   if( XDiffOrder == 0 && YDiffOrder == 0 )
+      return -0.5 * sin( M_PI * x) * sin( M_PI * x) * ( 1 - ( y - 2 ) * ( y - 2 ) ) * ( 1 - tanh ( 10 * ( sqrt( x * x + y * y ) - 0.6 ) ) );
+   return 0.0;
+}
+
+/****
+ * 3D
+ */
+
+template< typename Real >
+tnlString
+tnlTwinsFunction< 3, Real >::getType()
+{
+   return tnlString( "tnlTwinsFunction< 3, " ) + ::getType< Real >() + " >";
+}
+
+template< typename Real >
+tnlTwinsFunction< 3, Real >::tnlTwinsFunction()
+{
+}
+
+template< typename Real >
+   template< int XDiffOrder,
+             int YDiffOrder,
+             int ZDiffOrder,
+             typename Vertex >
+#ifdef HAVE_CUDA
+      __device__ __host__
+#endif
+Real
+tnlTwinsFunction< 3, Real >::
+getValue( const Vertex& v,
+          const Real& time ) const
+{
+   const RealType& x = v.x();
+   const RealType& y = v.y();
+   const RealType& z = v.z();
+   if( XDiffOrder == 0 && YDiffOrder == 0 && ZDiffOrder == 0 )
+      return 0.0;
+   return 0.0;
+}
+
+
+#endif /* TNLTWINSFUNCTION_IMPL_H_ */
diff --git a/src/functions/initial_conditions/tnlTwinsFunction_impl.h b/src/functors/initial_conditions/tnlTwinsFunction_impl.h~
similarity index 100%
rename from src/functions/initial_conditions/tnlTwinsFunction_impl.h
rename to src/functors/initial_conditions/tnlTwinsFunction_impl.h~
diff --git a/src/functors/tnlTestFunction_impl.h b/src/functors/tnlTestFunction_impl.h
index 23a05de8f82964f242a6b4d1e814a0027b4769f4..4e7ad770ba309970ee55aaebca9a6482d753ed1c 100644
--- a/src/functors/tnlTestFunction_impl.h
+++ b/src/functors/tnlTestFunction_impl.h
@@ -23,11 +23,11 @@
 #include <functors/tnlExpBumpFunction.h>
 #include <functors/tnlSinBumpsFunction.h>
 #include <functors/tnlSinWaveFunction.h>
-#include <functions/initial_conditions/tnlCylinderFunction.h>
-#include <functions/initial_conditions/tnlFlowerpotFunction.h>
-#include <functions/initial_conditions/tnlTwinsFunction.h>
-#include <functions/initial_conditions/level_set_functions/tnlBlobFunction.h>
-#include <functions/initial_conditions/level_set_functions/tnlPseudoSquareFunction.h>
+#include <functors/initial_conditions/tnlCylinderFunction.h>
+#include <functors/initial_conditions/tnlFlowerpotFunction.h>
+#include <functors/initial_conditions/tnlTwinsFunction.h>
+#include <functors/initial_conditions/level_set_functions/tnlBlobFunction.h>
+#include <functors/initial_conditions/level_set_functions/tnlPseudoSquareFunction.h>
 
 template< int FunctionDimensions,
           typename Real,
diff --git a/src/operators/diffusion/tnlExactNonlinearDiffusion.h b/src/operators/diffusion/tnlExactNonlinearDiffusion.h
index 07da1dc96c14344f73be5b16e480a7f2bc4b2395..76a4f8f57f286c2d8cd30c3967c03c5c52f2afcc 100644
--- a/src/operators/diffusion/tnlExactNonlinearDiffusion.h
+++ b/src/operators/diffusion/tnlExactNonlinearDiffusion.h
@@ -18,7 +18,7 @@
 #ifndef TNLEXACTNONLINEARDIFFUSION_H_
 #define TNLEXACTNONLINEARDIFFUSION_H_
 
-#include <functions/tnlFunctionType.h>
+#include <functors/tnlFunctionType.h>
 
 template< typename OperatorQ, int Dimensions >
 class tnlExactNonlinearDiffusion
diff --git a/src/operators/operator-Q/tnlExactOperatorQ.h b/src/operators/operator-Q/tnlExactOperatorQ.h
index e9e731aa759e3534f8be94228d1ed3d018c158ed..e319024f3f89440d67673172d7d4b46b1af21f95 100644
--- a/src/operators/operator-Q/tnlExactOperatorQ.h
+++ b/src/operators/operator-Q/tnlExactOperatorQ.h
@@ -4,7 +4,7 @@
 #include <core/vectors/tnlVector.h>
 #include <core/vectors/tnlSharedVector.h>
 #include <mesh/tnlGrid.h>
-#include <functions/tnlFunctionType.h>
+#include <functors/tnlFunctionType.h>
 
 template< int Dimensions >
 class tnlExactOperatorQ
diff --git a/src/operators/operator-Q/tnlOneSideDiffOperatorQ.h b/src/operators/operator-Q/tnlOneSideDiffOperatorQ.h
index dc91446d82aeb7058f06a832e9f3e4fe3775f8b6..954152b9a163c9a382752cd5a80145f61e275945 100644
--- a/src/operators/operator-Q/tnlOneSideDiffOperatorQ.h
+++ b/src/operators/operator-Q/tnlOneSideDiffOperatorQ.h
@@ -62,7 +62,7 @@ class tnlOneSideDiffOperatorQ< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, I
           const Vector& u,
           const Real& time )const;
           
-   void setEps(const Real& eps);
+   bool setEps(const Real& eps);
       
    private:
    
@@ -117,7 +117,7 @@ class tnlOneSideDiffOperatorQ< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, I
           const Vector& u,
           const Real& time )const;
         
-   void setEps(const Real& eps);
+   bool setEps(const Real& eps);
    
    private:
    
@@ -172,7 +172,7 @@ class tnlOneSideDiffOperatorQ< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, I
           const Vector& u,
           const Real& time ) const;
         
-   void setEps(const Real& eps);
+   bool setEps(const Real& eps);
    
    private:
    
@@ -225,7 +225,7 @@ class tnlOneSideDiffOperatorQ< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, I
           const Vector& u,
           const Real& time ) const;
           
-   void setEps(const Real& eps);
+   bool setEps(const Real& eps);
    
    private:
    
@@ -282,7 +282,7 @@ class tnlOneSideDiffOperatorQ< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, I
           const Vector& u,
           const Real& time )const;
           
-   void setEps(const Real& eps);
+   bool setEps(const Real& eps);
    
    private:
    
@@ -338,7 +338,7 @@ class tnlOneSideDiffOperatorQ< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, I
           const Vector& u,
           const Real& time )const;
           
-   void setEps(const Real& eps);
+   bool setEps(const Real& eps);
    
    private:
    
diff --git a/src/operators/operator-Q/tnlOneSideDiffOperatorQ_impl.h b/src/operators/operator-Q/tnlOneSideDiffOperatorQ_impl.h
index ba76eec3ebe86c4e45aca482d9418a75c3ed2f92..eba301d2baaceaa9353358dfe5541379226e06b8 100644
--- a/src/operators/operator-Q/tnlOneSideDiffOperatorQ_impl.h
+++ b/src/operators/operator-Q/tnlOneSideDiffOperatorQ_impl.h
@@ -40,9 +40,11 @@ template< typename MeshReal,
           typename MeshIndex,
           typename Real,
           typename Index >
-void tnlOneSideDiffOperatorQ< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 0 >::setEps( const Real& eps )
+bool tnlOneSideDiffOperatorQ< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 0 >::setEps( const Real& eps )
 {
   this->eps = eps;
+  
+  return true;
 }
 
 template< typename MeshReal,
@@ -50,9 +52,11 @@ template< typename MeshReal,
           typename MeshIndex,
           typename Real,
           typename Index >
-void tnlOneSideDiffOperatorQ< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 1 >::setEps( const Real& eps )
+bool tnlOneSideDiffOperatorQ< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index, 1 >::setEps( const Real& eps )
 {
   this->eps = eps;
+  
+  return true;
 }
 
 template< typename MeshReal,
@@ -215,9 +219,11 @@ template< typename MeshReal,
           typename MeshIndex,
           typename Real,
           typename Index >
-void tnlOneSideDiffOperatorQ< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index, 1 >::setEps( const Real& eps )
+bool tnlOneSideDiffOperatorQ< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index, 1 >::setEps( const Real& eps )
 {
   this->eps = eps;
+  
+  return true;
 }
 
 template< typename MeshReal,
@@ -225,9 +231,11 @@ template< typename MeshReal,
           typename MeshIndex,
           typename Real,
           typename Index >
-void tnlOneSideDiffOperatorQ< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index, 0 >::setEps( const Real& eps )
+bool tnlOneSideDiffOperatorQ< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index, 0 >::setEps( const Real& eps )
 {
   this->eps = eps;
+  
+  return true;
 }
 
 template< typename MeshReal,
@@ -397,9 +405,11 @@ template< typename MeshReal,
           typename MeshIndex,
           typename Real,
           typename Index >
-void tnlOneSideDiffOperatorQ< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index, 1 >::setEps( const Real& eps )
+bool tnlOneSideDiffOperatorQ< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index, 1 >::setEps( const Real& eps )
 {
   this->eps = eps;
+  
+  return true;
 }
 
 template< typename MeshReal,
@@ -407,9 +417,11 @@ template< typename MeshReal,
           typename MeshIndex,
           typename Real,
           typename Index >
-void tnlOneSideDiffOperatorQ< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index, 0 >::setEps( const Real& eps )
+bool tnlOneSideDiffOperatorQ< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index, 0 >::setEps( const Real& eps )
 {
   this->eps = eps;
+  
+  return true;
 }
 
 template< typename MeshReal,
diff --git a/src/problems/tnlMeanCurvatureFlowEocProblem_impl.h b/src/problems/tnlMeanCurvatureFlowEocProblem_impl.h
index dd7e1fb23036a5f3cde940c888df476d6a36e731..50c4ddcd7dfbd4554cba0da08262b41b67e06087 100644
--- a/src/problems/tnlMeanCurvatureFlowEocProblem_impl.h
+++ b/src/problems/tnlMeanCurvatureFlowEocProblem_impl.h
@@ -38,8 +38,10 @@ tnlMeanCurvatureFlowEocProblem< Mesh, BoundaryCondition, RightHandSide, Differen
 setup( const tnlParameterContainer& parameters )
 {
    if( ! this->boundaryCondition.setup( parameters ) ||
-       ! this->rightHandSide.setup( parameters ) )
+       ! this->rightHandSide.setup( parameters ) || 
+       ! this->differentialOperator.nonlinearDiffusionOperator.operatorQ.setEps(parameters.getParameter< double >("eps")) )
       return false;
+   
    return true;
 }
 
diff --git a/src/problems/tnlMeanCurvatureFlowEocRhs.h b/src/problems/tnlMeanCurvatureFlowEocRhs.h
index 60910448ee9da96d8fe413b8f90b943b39c5a761..0f97a4323c07c679b471ff65262f8eb2d2038d3b 100644
--- a/src/problems/tnlMeanCurvatureFlowEocRhs.h
+++ b/src/problems/tnlMeanCurvatureFlowEocRhs.h
@@ -18,7 +18,7 @@
 #ifndef TNLMEANCURVATUREFLOWEOCRHS_H_
 #define TNLMEANCURVATUREFLOWEOCRHS_H_
 
-#include <functions/tnlFunctionType.h>
+#include <functors/tnlFunctionType.h>
 
 template< typename ExactOperator,
           typename TestFunction >
diff --git a/src/problems/tnlMeanCurvatureFlowProblem.h b/src/problems/tnlMeanCurvatureFlowProblem.h
index 2522aac23e3687e16f1c575efc0bf309139735f8..3541cccd8257d82baca804ea21389140ccebaaed 100644
--- a/src/problems/tnlMeanCurvatureFlowProblem.h
+++ b/src/problems/tnlMeanCurvatureFlowProblem.h
@@ -46,6 +46,7 @@ class tnlMeanCurvatureFlowProblem : public tnlPDEProblem< Mesh,
 
       using typename BaseType::MeshType;
       using typename BaseType::DofVectorType;
+      using typename BaseType::MeshDependentDataType;
 
       static tnlString getTypeStatic();
 
@@ -59,7 +60,7 @@ class tnlMeanCurvatureFlowProblem : public tnlPDEProblem< Mesh,
       bool setInitialCondition( const tnlParameterContainer& parameters,
                                 const MeshType& mesh,
                                 DofVectorType& dofs,
-                                DofVectorType& auxDofs );
+                                MeshDependentDataType& meshDependentData );
 
       template< typename Matrix >
       bool setupLinearSystem( const MeshType& mesh,
@@ -69,7 +70,7 @@ class tnlMeanCurvatureFlowProblem : public tnlPDEProblem< Mesh,
                          const IndexType& step,
                          const MeshType& mesh,
                          DofVectorType& dofs,
-                         DofVectorType& auxDofs );
+                         MeshDependentDataType& meshDependentData );
 
       IndexType getDofs( const MeshType& mesh ) const;
 
@@ -80,14 +81,15 @@ class tnlMeanCurvatureFlowProblem : public tnlPDEProblem< Mesh,
                            const RealType& tau,
                            const MeshType& mesh,
                            DofVectorType& _u,
-                           DofVectorType& _fu );
+                           DofVectorType& _fu,
+			   MeshDependentDataType& meshDependentData );
 
       template< typename Matrix >
       void assemblyLinearSystem( const RealType& time,
                                  const RealType& tau,
                                  const MeshType& mesh,
                                  DofVectorType& dofs,
-                                 DofVectorType& auxDofs,
+                                 MeshDependentDataType& meshDependentData,
                                  Matrix& matrix,
                                  DofVectorType& rightHandSide );
 
diff --git a/src/problems/tnlMeanCurvatureFlowProblem_impl.h b/src/problems/tnlMeanCurvatureFlowProblem_impl.h
index aff7a2581d27325d063448c2a312eb6a55991d07..e941919bc4aa11052ee5556361801b7a4103bcad 100644
--- a/src/problems/tnlMeanCurvatureFlowProblem_impl.h
+++ b/src/problems/tnlMeanCurvatureFlowProblem_impl.h
@@ -24,6 +24,7 @@
 #include <core/tnlLogger.h>
 #include <solvers/pde/tnlExplicitUpdater.h>
 #include <solvers/pde/tnlLinearSystemAssembler.h>
+#include <solvers/pde/tnlBackwardTimeDiscretisation.h>
 
 #include "tnlMeanCurvatureFlowProblem.h"
 
@@ -69,11 +70,10 @@ tnlMeanCurvatureFlowProblem< Mesh, BoundaryCondition, RightHandSide, Differentia
 setup( const tnlParameterContainer& parameters )
 {
    if( ! this->boundaryCondition.setup( parameters, "boundary-conditions-" ) ||
-       ! this->rightHandSide.setup( parameters, "right-hand-side-" ) )
+       ! this->rightHandSide.setup( parameters, "right-hand-side-" ) ||
+       ! this->differentialOperator.nonlinearDiffusionOperator.operatorQ.setEps( parameters.getParameter< double >( "eps" ) ) )
       return false;
    
-   differentialOperator.nonlinearDiffusionOperator.operatorQ.setEps(parameters.getParameter< double >("eps"));
-   
    return true;
 }
 
@@ -115,7 +115,7 @@ tnlMeanCurvatureFlowProblem< Mesh, BoundaryCondition, RightHandSide, Differentia
 setInitialCondition( const tnlParameterContainer& parameters,
                      const MeshType& mesh,
                      DofVectorType& dofs,
-                     DofVectorType& auxiliaryDofs )
+                     MeshDependentDataType& meshDependentData )
 {
    this->bindDofs( mesh, dofs );
    const tnlString& initialConditionFile = parameters.getParameter< tnlString >( "initial-condition" );
@@ -163,7 +163,7 @@ makeSnapshot( const RealType& time,
               const IndexType& step,
               const MeshType& mesh,
               DofVectorType& dofs,
-              DofVectorType& auxiliaryDofs )
+              MeshDependentDataType& meshDependentData )
 {
    cout << endl << "Writing output at time " << time << " step " << step << "." << endl;
 
@@ -186,7 +186,8 @@ getExplicitRHS( const RealType& time,
                 const RealType& tau,
                 const MeshType& mesh,
                 DofVectorType& u,
-                DofVectorType& fu )
+                DofVectorType& fu,
+		MeshDependentDataType& meshDependentData )
 {
    /****
     * If you use an explicit solver like tnlEulerSolver or tnlMersonSolver, you
@@ -228,11 +229,17 @@ assemblyLinearSystem( const RealType& time,
                       const RealType& tau,
                       const MeshType& mesh,
                       DofVectorType& u,
-                      DofVectorType& auxDofs,
+                      MeshDependentDataType& meshDependentData,
                       Matrix& matrix,
                       DofVectorType& b )
 {
-   tnlLinearSystemAssembler< Mesh, DofVectorType, DifferentialOperator, BoundaryCondition, RightHandSide, MatrixType > systemAssembler;
+   tnlLinearSystemAssembler< Mesh, 
+			     DofVectorType, 
+			     DifferentialOperator, 
+			     BoundaryCondition, 
+			     RightHandSide, 
+			     tnlBackwardTimeDiscretisation, 
+			     MatrixType > systemAssembler;
    systemAssembler.template assembly< Mesh::Dimensions >( time,
                                                           tau,
                                                           mesh,