diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index b0d5d7f7d3f32843b3981c157aee52e5bff13fd5..80e1cf9833e94a871c812d04cc1ede21ab3ebaa3 100755
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1,3 +1,3 @@
 add_subdirectory( heat-equation )
 add_subdirectory( navier-stokes )
-#add_subdirectory( mean-curvature-flow )
+add_subdirectory( mean-curvature-flow )
diff --git a/src/operators/diffusion/tnlExactNonlinearDiffusion.h b/src/operators/diffusion/tnlExactNonlinearDiffusion.h
index fbc31a35880866e2c45ea280ae74bf4127b0d36c..1522343b231f2735e99d54d463191c218f9e65ab 100644
--- a/src/operators/diffusion/tnlExactNonlinearDiffusion.h
+++ b/src/operators/diffusion/tnlExactNonlinearDiffusion.h
@@ -18,14 +18,14 @@
 #ifndef TNLEXACTNONLINEARDIFFUSION_H_
 #define TNLEXACTNONLINEARDIFFUSION_H_
 
-#include <functions/tnlFunction.h>
+#include <functions/tnlDomain.h>
 
 template< typename OperatorQ, int Dimensions >
 class tnlExactNonlinearDiffusion
 {};
 
 template< typename OperatorQ >
-class tnlExactNonlinearDiffusion< OperatorQ, 1 > : public tnlFunction< 1, SpaceDomain >
+class tnlExactNonlinearDiffusion< OperatorQ, 1 > : public tnlDomain< 1, SpaceDomain >
 {
    public:
 
@@ -40,13 +40,13 @@ class tnlExactNonlinearDiffusion< OperatorQ, 1 > : public tnlFunction< 1, SpaceD
 #endif      
       
       __cuda_callable__
-      static Real operator()( const Function& function,
-                            const Vertex& v,
-                            const Real& time = 0.0 );
+      Real operator()( const Function& function,
+                       const Vertex& v,
+                       const Real& time = 0.0 ) const;
 };
 
 template< typename OperatorQ >
-class tnlExactNonlinearDiffusion< OperatorQ, 2 > : public tnlFunction< 2, SpaceDomain >
+class tnlExactNonlinearDiffusion< OperatorQ, 2 > : public tnlDomain< 2, SpaceDomain >
 {
    public:
 
@@ -61,13 +61,13 @@ class tnlExactNonlinearDiffusion< OperatorQ, 2 > : public tnlFunction< 2, SpaceD
 #endif 
 
       __cuda_callable__
-      static Real operator()( const Function& function,
-                            const Vertex& v,
-                            const Real& time = 0.0 );
+      Real operator()( const Function& function,
+                       const Vertex& v,
+                       const Real& time = 0.0 ) const;
 };
 
 template< typename OperatorQ >
-class tnlExactNonlinearDiffusion< OperatorQ, 3 > : public tnlFunction< 3, SpaceDomain >
+class tnlExactNonlinearDiffusion< OperatorQ, 3 > : public tnlDomain< 3, SpaceDomain >
 {
    public:
 
@@ -82,9 +82,9 @@ class tnlExactNonlinearDiffusion< OperatorQ, 3 > : public tnlFunction< 3, SpaceD
 #endif 
 
       __cuda_callable__
-      static Real operator()( const Function& function,
-                            const Vertex& v,
-                            const Real& time = 0.0 );
+      Real operator()( const Function& function,
+                       const Vertex& v,
+                       const Real& time = 0.0 ) const;
 };
 
 #include "tnlExactNonlinearDiffusion_impl.h"
diff --git a/src/operators/diffusion/tnlExactNonlinearDiffusion_impl.h b/src/operators/diffusion/tnlExactNonlinearDiffusion_impl.h
index 64c9d7013327aeca39b8ad2ca902149228f98727..1444f124505e48bee71dfd82132a7819582eb71f 100644
--- a/src/operators/diffusion/tnlExactNonlinearDiffusion_impl.h
+++ b/src/operators/diffusion/tnlExactNonlinearDiffusion_impl.h
@@ -33,7 +33,7 @@ Real
 tnlExactNonlinearDiffusion< OperatorQ, 1 >::
 operator()( const Function& function,
           const Vertex& v,
-          const Real& time )
+          const Real& time ) const
 {
    return function.template getPartialDerivative< 2, 0, 0 >( v, time ) - function.template getPartialDerivative< 1, 0, 0 >( v, time ) 
           * OperatorQ::template getPartialDerivative<1, 0, 0>(function, v, time ) / OperatorQ::template getPartialDerivative<0, 0, 0>(function, v, time );
@@ -54,7 +54,7 @@ Real
 tnlExactNonlinearDiffusion< OperatorQ, 2 >::
 operator()( const Function& function,
           const Vertex& v,
-          const Real& time )
+          const Real& time ) const
 {
    return  function.template getPartialDerivative< 2, 0, 0 >( v, time ) +  function.template getPartialDerivative< 0, 2, 0 >( v, time )
            -( OperatorQ::template getPartialDerivative<1, 0, 0> (function, v, time) * function.template getPartialDerivative< 1, 0, 0 >( v, time ) 
@@ -77,7 +77,7 @@ Real
 tnlExactNonlinearDiffusion< OperatorQ, 3 >::
 operator()( const Function& function,
           const Vertex& v,
-          const Real& time )
+          const Real& time ) const
 {
    return  function.template getPartialDerivative< 2, 0, 0 >( v, time ) +  function.template getPartialDerivative< 0, 2, 0 >( v, time )
            +  function.template getPartialDerivative< 0, 0, 2 >( v, time )
diff --git a/src/operators/diffusion/tnlNonlinearDiffusion_impl.h b/src/operators/diffusion/tnlNonlinearDiffusion_impl.h
index 8b4a05cbe90d949ce265f208cec72bf7e4c84197..2f2bbec8806a702d8e3e8b9dfe4c7a9a450dab00 100644
--- a/src/operators/diffusion/tnlNonlinearDiffusion_impl.h
+++ b/src/operators/diffusion/tnlNonlinearDiffusion_impl.h
@@ -38,7 +38,7 @@ getValue( const MeshEntity& entity,
           const Vector& u,
           const Real& time ) const
 {
-    return nonlinearDiffusionOperator.getValue( entity, u, time );
+    return nonlinearDiffusionOperator.getValue( u, entity, time );
 }
 
 template< typename MeshReal,
@@ -115,7 +115,7 @@ getValue( const MeshEntity& entity,
           const Vector& u,
           const Real& time ) const
 {
-    return nonlinearDiffusionOperator.getValue( entity, u, time );
+    return nonlinearDiffusionOperator.getValue( u, entity, time );
 }
        
 template< typename MeshReal,
@@ -192,7 +192,7 @@ getValue( const MeshEntity& entity,
           const Vector& u,
           const Real& time ) const
 {
-    return nonlinearDiffusionOperator.getValue( entity, u, time );
+    return nonlinearDiffusionOperator.getValue( u, entity, time );
 }
 
 template< typename MeshReal,
diff --git a/src/operators/operator-Q/tnlExactOperatorQ.h b/src/operators/operator-Q/tnlExactOperatorQ.h
index 3ad7bedb4ce92444d4ac6b600d975c3faf6216a2..5aacd5c796a9b18846f5c1be16bb21ab635c7e9f 100644
--- a/src/operators/operator-Q/tnlExactOperatorQ.h
+++ b/src/operators/operator-Q/tnlExactOperatorQ.h
@@ -4,14 +4,14 @@
 #include <core/vectors/tnlVector.h>
 #include <core/vectors/tnlSharedVector.h>
 #include <mesh/tnlGrid.h>
-#include <functions/tnlFunction.h>
+#include <functions/tnlDomain.h>
 
 template< int Dimensions >
 class tnlExactOperatorQ
 {};
 
 template<>
-class tnlExactOperatorQ< 1 >
+class tnlExactOperatorQ< 1 > : public tnlDomain< 1, SpaceDomain >
 {
    public:
 
@@ -24,12 +24,11 @@ class tnlExactOperatorQ< 1 >
 #else   
       template< int XDiffOrder = 0, int YDiffOrder = 0, int ZDiffOrder = 0, 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, const Real& eps = 1.0 );
+      __cuda_callable__
+      static Real getPartialDerivative( const Function& function,
+                                        const Vertex& v,
+                                        const Real& time = 0.0,
+                                        const Real& eps = 1.0 );
       
 };
 
@@ -47,12 +46,12 @@ class tnlExactOperatorQ< 2 >
 #else   
       template< int XDiffOrder = 0, int YDiffOrder = 0, int ZDiffOrder = 0, 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, const Real& eps = 1.0 );
+
+      __cuda_callable__
+      static Real getPartialDerivative( const Function& function,
+                                        const Vertex& v,
+                                        const Real& time = 0.0,
+                                        const Real& eps = 1.0 );
 };
 
 template<>
@@ -69,12 +68,12 @@ class tnlExactOperatorQ< 3 >
 #else   
       template< int XDiffOrder = 0, int YDiffOrder = 0, int ZDiffOrder = 0, 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, const Real& eps = 1.0 );
+
+      __cuda_callable__
+      static Real getPartialDerivative( const Function& function,
+                                        const Vertex& v,
+                                        const Real& time = 0.0,
+                                        const Real& eps = 1.0 );
 };
 
 #include <operators/operator-Q/tnlExactOperatorQ_impl.h>
diff --git a/src/operators/operator-Q/tnlExactOperatorQ_impl.h b/src/operators/operator-Q/tnlExactOperatorQ_impl.h
index 19fdbd917aea0fc0184b25678e7f54175550ff3b..72d72943fe1195e2de6a62a9b400fbe7d40b1ac0 100644
--- a/src/operators/operator-Q/tnlExactOperatorQ_impl.h
+++ b/src/operators/operator-Q/tnlExactOperatorQ_impl.h
@@ -28,21 +28,19 @@ getType()
 }
 
 template< int XDiffOrder, int YDiffOrder, int ZDiffOrder, typename Function, typename Vertex, typename Real >
-#ifdef HAVE_CUDA
-   __device__ __host__
-#endif
+__cuda_callable__
 Real
 tnlExactOperatorQ< 1 >::
-getValue( const Function& function,
-          const Vertex& v,
-          const Real& time, const Real& eps )
+getPartialDerivative( const Function& function,
+                      const Vertex& v,
+                      const Real& time, const Real& eps )
 {
    if( YDiffOrder != 0 || ZDiffOrder != 0 )
         return 0.0;
    if (XDiffOrder == 0)
-        return sqrt(eps * eps + function.template getValue< 1, 0, 0 >( v, time ) * function.template getValue< 1, 0, 0 >( v, time ) );
+        return sqrt(eps * eps + function.template getPartialDerivative< 1, 0, 0 >( v, time ) * function.template getPartialDerivative< 1, 0, 0 >( v, time ) );
    if (XDiffOrder == 1)
-        return (function.template getValue< 1, 0, 0 >( v, time ) * function.template getValue< 2, 0, 0 >( v, time ) ) / getValue( function, v, time, eps);
+        return (function.template getPartialDerivative< 1, 0, 0 >( v, time ) * function.template getPartialDerivative< 2, 0, 0 >( v, time ) ) / getPartialDerivative( function, v, time, eps);
    return 0;
 }
 
@@ -54,26 +52,25 @@ getType()
 }
 
 template< int XDiffOrder, int YDiffOrder, int ZDiffOrder, typename Function, typename Vertex, typename Real >
-#ifdef HAVE_CUDA
-   __device__ __host__
-#endif
+
+__cuda_callable__
 Real
 tnlExactOperatorQ< 2 >::
-getValue( const Function& function,
-          const Vertex& v,
-          const Real& time, const Real& eps )
+getPartialDerivative( const Function& function,
+                      const Vertex& v,
+                      const Real& time, const Real& eps )
 {
    if( ZDiffOrder != 0 )
         return 0.0;
    if (XDiffOrder == 0 && YDiffOrder == 0 )
-        return sqrt(eps * eps + function.template getValue< 1, 0, 0 >( v, time ) * function.template getValue< 1, 0, 0 >( v, time ) 
-                + function.template getValue< 0, 1, 0 >( v, time ) * function.template getValue< 0, 1, 0 >( v, time ) );
+        return sqrt(eps * eps + function.template getPartialDerivative< 1, 0, 0 >( v, time ) * function.template getPartialDerivative< 1, 0, 0 >( v, time ) 
+                + function.template getPartialDerivative< 0, 1, 0 >( v, time ) * function.template getPartialDerivative< 0, 1, 0 >( v, time ) );
    if (XDiffOrder == 1 && YDiffOrder == 0 )
-        return (function.template getValue< 1, 0, 0 >( v, time ) * function.template getValue< 2, 0, 0 >( v, time ) + 
-                function.template getValue< 0, 1, 0 >( v, time ) * function.template getValue< 1, 1, 0 >( v, time )) / getValue( function, v, time, eps);
+        return (function.template getPartialDerivative< 1, 0, 0 >( v, time ) * function.template getPartialDerivative< 2, 0, 0 >( v, time ) + 
+                function.template getPartialDerivative< 0, 1, 0 >( v, time ) * function.template getPartialDerivative< 1, 1, 0 >( v, time )) / getPartialDerivative( function, v, time, eps);
    if (XDiffOrder == 0 && YDiffOrder == 1 )
-        return (function.template getValue< 0, 1, 0 >( v, time ) * function.template getValue< 0, 2, 0 >( v, time ) + 
-                function.template getValue< 1, 0, 0 >( v, time ) * function.template getValue< 1, 1, 0 >( v, time )) / getValue( function, v, time, eps);
+        return (function.template getPartialDerivative< 0, 1, 0 >( v, time ) * function.template getPartialDerivative< 0, 2, 0 >( v, time ) + 
+                function.template getPartialDerivative< 1, 0, 0 >( v, time ) * function.template getPartialDerivative< 1, 1, 0 >( v, time )) / getPartialDerivative( function, v, time, eps);
    return 0;
 }
 
@@ -85,31 +82,30 @@ getType()
 }
 
 template< int XDiffOrder, int YDiffOrder, int ZDiffOrder, typename Function, typename Vertex, typename Real >
-#ifdef HAVE_CUDA
-   __device__ __host__
-#endif
+
+__cuda_callable__
 Real
 tnlExactOperatorQ< 3 >::
-getValue( const Function& function,
-          const Vertex& v,
-          const Real& time, const Real& eps )
+getPartialDerivative( const Function& function,
+                      const Vertex& v,
+                      const Real& time, const Real& eps )
 {
    if ( XDiffOrder == 0 && YDiffOrder == 0  && ZDiffOrder == 0 )
-        return sqrt(eps * eps + function.template getValue< 1, 0, 0 >( v, time ) * function.template getValue< 1, 0, 0 >( v, time ) 
-                + function.template getValue< 0, 1, 0 >( v, time ) * function.template getValue< 0, 1, 0 >( v, time )
-                + function.template getValue< 0, 0, 1 >( v, time ) * function.template getValue< 0, 0, 1 >( v, time ) );
+        return sqrt(eps * eps + function.template getPartialDerivative< 1, 0, 0 >( v, time ) * function.template getPartialDerivative< 1, 0, 0 >( v, time ) 
+                + function.template getPartialDerivative< 0, 1, 0 >( v, time ) * function.template getPartialDerivative< 0, 1, 0 >( v, time )
+                + function.template getPartialDerivative< 0, 0, 1 >( v, time ) * function.template getPartialDerivative< 0, 0, 1 >( v, time ) );
    if (XDiffOrder == 1 && YDiffOrder == 0 && ZDiffOrder == 0 )
-        return (function.template getValue< 1, 0, 0 >( v, time ) * function.template getValue< 2, 0, 0 >( v, time ) + 
-                function.template getValue< 0, 1, 0 >( v, time ) * function.template getValue< 1, 1, 0 >( v, time ) + 
-                function.template getValue< 0, 0, 1 >( v, time ) * function.template getValue< 1, 0, 1 >( v, time )) / getValue( function, v, time, eps);
+        return (function.template getPartialDerivative< 1, 0, 0 >( v, time ) * function.template getPartialDerivative< 2, 0, 0 >( v, time ) + 
+                function.template getPartialDerivative< 0, 1, 0 >( v, time ) * function.template getPartialDerivative< 1, 1, 0 >( v, time ) + 
+                function.template getPartialDerivative< 0, 0, 1 >( v, time ) * function.template getPartialDerivative< 1, 0, 1 >( v, time )) / getPartialDerivative( function, v, time, eps);
    if (XDiffOrder == 0 && YDiffOrder == 1 && ZDiffOrder == 0 )
-        return (function.template getValue< 1, 0, 0 >( v, time ) * function.template getValue< 1, 1, 0 >( v, time ) + 
-                function.template getValue< 0, 1, 0 >( v, time ) * function.template getValue< 0, 2, 0 >( v, time ) + 
-                function.template getValue< 0, 0, 1 >( v, time ) * function.template getValue< 0, 1, 1 >( v, time )) / getValue( function, v, time, eps);
+        return (function.template getPartialDerivative< 1, 0, 0 >( v, time ) * function.template getPartialDerivative< 1, 1, 0 >( v, time ) + 
+                function.template getPartialDerivative< 0, 1, 0 >( v, time ) * function.template getPartialDerivative< 0, 2, 0 >( v, time ) + 
+                function.template getPartialDerivative< 0, 0, 1 >( v, time ) * function.template getPartialDerivative< 0, 1, 1 >( v, time )) / getPartialDerivative( function, v, time, eps);
    if (XDiffOrder == 0 && YDiffOrder == 0 && ZDiffOrder == 1 )
-        return (function.template getValue< 1, 0, 0 >( v, time ) * function.template getValue< 1, 0, 1 >( v, time ) + 
-                function.template getValue< 0, 1, 0 >( v, time ) * function.template getValue< 0, 1, 1 >( v, time ) + 
-                function.template getValue< 0, 0, 1 >( v, time ) * function.template getValue< 0, 0, 2 >( v, time )) / getValue( function, v, time, eps);
+        return (function.template getPartialDerivative< 1, 0, 0 >( v, time ) * function.template getPartialDerivative< 1, 0, 1 >( v, time ) + 
+                function.template getPartialDerivative< 0, 1, 0 >( v, time ) * function.template getPartialDerivative< 0, 1, 1 >( v, time ) + 
+                function.template getPartialDerivative< 0, 0, 1 >( v, time ) * function.template getPartialDerivative< 0, 0, 2 >( v, time )) / getPartialDerivative( function, v, time, eps);
    return 0;
 }
 
diff --git a/src/problems/tnlMeanCurvatureFlowEocRhs.h b/src/problems/tnlMeanCurvatureFlowEocRhs.h
index afd883975b97e1cd71bd5f101f75d1d35090f890..2f1d3521e5a578e52041fd1839e67e34d2b86ac9 100644
--- a/src/problems/tnlMeanCurvatureFlowEocRhs.h
+++ b/src/problems/tnlMeanCurvatureFlowEocRhs.h
@@ -18,12 +18,12 @@
 #ifndef TNLMEANCURVATUREFLOWEOCRHS_H_
 #define TNLMEANCURVATUREFLOWEOCRHS_H_
 
-#include <functions/tnlFunction.h>
+#include <functions/tnlDomain.h>
 
 template< typename ExactOperator,
           typename TestFunction,
           int Dimensions >
-class tnlMeanCurvatureFlowEocRhs : public tnlFunction< Dimensions, SpaceDomain >
+class tnlMeanCurvatureFlowEocRhs : public tnlDomain< Dimensions, SpaceDomain >
 {
    public:
 
@@ -43,11 +43,11 @@ class tnlMeanCurvatureFlowEocRhs : public tnlFunction< Dimensions, SpaceDomain >
       template< typename Vertex,
                 typename Real >
       __cuda_callable__
-      Real getValue( const Vertex& vertex,
-                     const Real& time ) const
+      Real operator()( const Vertex& vertex,
+                       const Real& time ) const
       {
          return testFunction.getTimeDerivative( vertex, time )
-                - exactOperator.getValue( testFunction, vertex, time );
+                - exactOperator( testFunction, vertex, time );
       };
 
    protected:
diff --git a/src/problems/tnlMeanCurvatureFlowProblem_impl.h b/src/problems/tnlMeanCurvatureFlowProblem_impl.h
index 95d6f3883a25bcb71fcc66ddd6f18c9cf9353035..7f2ad74af530f7032b262f0e40e5869ce0dc2b49 100644
--- a/src/problems/tnlMeanCurvatureFlowProblem_impl.h
+++ b/src/problems/tnlMeanCurvatureFlowProblem_impl.h
@@ -218,11 +218,10 @@ getExplicitRHS( const RealType& time,
       u,
       fu );
    
-   tnlBoundaryConditionsSetter< Mesh, MeshFunctionType, BoundaryCondition > boundaryConditionsSetter;
+   tnlBoundaryConditionsSetter< MeshFunctionType, BoundaryCondition > boundaryConditionsSetter;
    boundaryConditionsSetter.template apply< typename Mesh::Cell >(
-      time,
-      mesh,
       this->boundaryCondition,
+      time,
       u );
 
    /*cout << "u = " << u << endl;