From 151756a4b42420864d256202eb1e23400ab4cec9 Mon Sep 17 00:00:00 2001
From: Tomas Oberhuber <tomas.oberhuber@fjfi.cvut.cz>
Date: Fri, 25 Dec 2015 15:30:27 +0100
Subject: [PATCH] Implementing mesh function. Refactoring the Dirichlet and the
 Neumann boundary conditions.

---
 .../heat-equation/tnl-heat-equation-eoc.h     |   6 +-
 examples/heat-equation/tnl-heat-equation.h    |  20 +-
 src/CMakeLists.txt                            |   6 +-
 src/core/vectors/tnlSharedVector.h            |   5 +-
 src/core/vectors/tnlVector.h                  |   5 +-
 src/{functors => functions}/CMakeLists.txt    |   9 +-
 .../tnlConstantFunction.h                     |   6 +-
 .../tnlConstantFunction_impl.h                |   0
 .../tnlExpBumpFunction.h                      |   6 +-
 .../tnlExpBumpFunction_impl.h                 |   2 +-
 src/{functors => functions}/tnlFunction.h     |  16 +-
 .../tnlFunctionAdapter.h                      | 100 +++--
 .../tnlFunctionDiscretizer.h                  |   2 +-
 .../tnlFunctionDiscretizer_impl.h             |   0
 .../tnlFunctionEnumerator.h                   |  10 +-
 .../tnlFunctionEnumerator_impl.h              |   2 +-
 src/functions/tnlMeshFunction.h               |  91 ++++
 src/functions/tnlMeshFunction_impl.h          | 176 ++++++++
 .../tnlSinBumpsFunction.h                     |   6 +-
 .../tnlSinBumpsFunction_impl.h                |   2 +-
 .../tnlSinWaveFunction.h                      |   6 +-
 .../tnlSinWaveFunction_impl.h                 |   2 +-
 src/{functors => functions}/tnlTestFunction.h |   5 +-
 .../tnlTestFunction_impl.cpp                  |   2 +-
 .../tnlTestFunction_impl.cu                   |   0
 .../tnlTestFunction_impl.h                    |   8 +-
 src/functors/tnlFunctorAdapter.h              | 278 -------------
 src/functors/tnlMeshFunction.h                |  72 ----
 src/functors/tnlMeshFunction_impl.h           | 115 -----
 src/matrices/tnlMatrixSetter.h                |  12 +-
 src/mesh/grids/tnlTraverser_Grid1D_impl.h     |   6 +-
 src/mesh/grids/tnlTraverser_Grid2D_impl.h     |  10 +-
 src/mesh/grids/tnlTraverser_Grid3D_impl.h     |  85 ++--
 src/mesh/tnlDummyMesh.h                       |   5 +-
 src/operators/CMakeLists.txt                  |   4 -
 .../diffusion/tnlExactLinearDiffusion.h       |   6 +-
 src/operators/diffusion/tnlLinearDiffusion.h  |   3 -
 .../diffusion/tnlLinearDiffusion_impl.h       |  38 +-
 .../tnlAnalyticDirichletBoundaryConditions.h  | 114 -----
 ...AnalyticDirichletBoundaryConditions_impl.h | 159 -------
 .../tnlAnalyticNeumannBoundaryConditions.h    | 227 ----------
 ...nlAnalyticNeumannBoundaryConditions_impl.h | 393 ------------------
 .../tnlDirichletBoundaryConditions.h          |  40 +-
 .../tnlDirichletBoundaryConditions_impl.h     |  73 ++--
 src/operators/tnlExactOperatorEvaluator.h     |  21 +-
 src/operators/tnlNeumannBoundaryConditions.h  |  43 +-
 .../tnlNeumannBoundaryConditions_impl.h       | 161 ++++---
 src/problems/tnlHeatEquationEocProblem_impl.h |   3 +
 src/problems/tnlHeatEquationEocRhs.h          |   7 +-
 src/solvers/pde/tnlExplicitUpdater.h          |  28 +-
 src/solvers/pde/tnlLinearSystemAssembler.h    |  22 +-
 .../diffusion/tnlLinearDiffusionTest.cpp      |   2 +-
 .../operators/tnlPDEOperatorEocTestSetter.h   |   2 +-
 tests/unit-tests/tnlApproximationError.h      |   8 +-
 tests/unit-tests/tnlApproximationError_impl.h |   2 +-
 tools/src/tnl-init.cpp                        |   2 +-
 tools/src/tnl-init.h                          |   4 +-
 57 files changed, 677 insertions(+), 1761 deletions(-)
 rename src/{functors => functions}/CMakeLists.txt (84%)
 rename src/{functors => functions}/tnlConstantFunction.h (94%)
 rename src/{functors => functions}/tnlConstantFunction_impl.h (100%)
 rename src/{functors => functions}/tnlExpBumpFunction.h (96%)
 rename src/{functors => functions}/tnlExpBumpFunction_impl.h (99%)
 rename src/{functors => functions}/tnlFunction.h (83%)
 rename src/{functors => functions}/tnlFunctionAdapter.h (55%)
 rename src/{functors => functions}/tnlFunctionDiscretizer.h (97%)
 rename src/{functors => functions}/tnlFunctionDiscretizer_impl.h (100%)
 rename src/{functors => functions}/tnlFunctionEnumerator.h (96%)
 rename src/{functors => functions}/tnlFunctionEnumerator_impl.h (99%)
 create mode 100644 src/functions/tnlMeshFunction.h
 create mode 100644 src/functions/tnlMeshFunction_impl.h
 rename src/{functors => functions}/tnlSinBumpsFunction.h (97%)
 rename src/{functors => functions}/tnlSinBumpsFunction_impl.h (99%)
 rename src/{functors => functions}/tnlSinWaveFunction.h (95%)
 rename src/{functors => functions}/tnlSinWaveFunction_impl.h (99%)
 rename src/{functors => functions}/tnlTestFunction.h (96%)
 rename src/{functors => functions}/tnlTestFunction_impl.cpp (97%)
 rename src/{functors => functions}/tnlTestFunction_impl.cu (100%)
 rename src/{functors => functions}/tnlTestFunction_impl.h (98%)
 delete mode 100644 src/functors/tnlFunctorAdapter.h
 delete mode 100644 src/functors/tnlMeshFunction.h
 delete mode 100644 src/functors/tnlMeshFunction_impl.h
 delete mode 100644 src/operators/tnlAnalyticDirichletBoundaryConditions.h
 delete mode 100644 src/operators/tnlAnalyticDirichletBoundaryConditions_impl.h
 delete mode 100644 src/operators/tnlAnalyticNeumannBoundaryConditions.h
 delete mode 100644 src/operators/tnlAnalyticNeumannBoundaryConditions_impl.h

diff --git a/examples/heat-equation/tnl-heat-equation-eoc.h b/examples/heat-equation/tnl-heat-equation-eoc.h
index a5b2bdccb6..04a9ade3a9 100644
--- a/examples/heat-equation/tnl-heat-equation-eoc.h
+++ b/examples/heat-equation/tnl-heat-equation-eoc.h
@@ -21,12 +21,12 @@
 #include <solvers/tnlSolver.h>
 #include <solvers/tnlFastBuildConfigTag.h>
 #include <solvers/tnlBuildConfigTags.h>
-#include <functors/tnlTestFunction.h>
+#include <functions/tnlTestFunction.h>
 #include <operators/diffusion/tnlLinearDiffusion.h>
 #include <operators/diffusion/tnlExactLinearDiffusion.h>
-#include <operators/tnlAnalyticDirichletBoundaryConditions.h>
 #include <problems/tnlHeatEquationEocRhs.h>
 #include <problems/tnlHeatEquationEocProblem.h>
+#include <operators/tnlDirichletBoundaryConditions.h>
 
 //typedef tnlDefaultBuildMeshConfig BuildConfig;
 typedef tnlFastBuildConfig BuildConfig;
@@ -67,7 +67,7 @@ class heatEquationSetter
       typedef tnlTestFunction< MeshType::meshDimensions, Real, Device > TestFunction;
       typedef tnlHeatEquationEocRhs< ExactOperator, TestFunction > RightHandSide;
       typedef tnlStaticVector < MeshType::meshDimensions, Real > Vertex;
-      typedef tnlAnalyticDirichletBoundaryConditions< MeshType, TestFunction, Real, Index > BoundaryConditions;
+      typedef tnlDirichletBoundaryConditions< MeshType, TestFunction, Real, Index > BoundaryConditions;
       typedef tnlHeatEquationEocProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Solver;
       SolverStarter solverStarter;
       return solverStarter.template run< Solver >( parameters );
diff --git a/examples/heat-equation/tnl-heat-equation.h b/examples/heat-equation/tnl-heat-equation.h
index 209ab2d577..02ef5c4431 100644
--- a/examples/heat-equation/tnl-heat-equation.h
+++ b/examples/heat-equation/tnl-heat-equation.h
@@ -22,12 +22,12 @@
 #include <solvers/tnlFastBuildConfigTag.h>
 #include <solvers/tnlBuildConfigTags.h>
 #include <operators/diffusion/tnlLinearDiffusion.h>
-#include <operators/tnlAnalyticDirichletBoundaryConditions.h>
 #include <operators/tnlDirichletBoundaryConditions.h>
-#include <operators/tnlAnalyticNeumannBoundaryConditions.h>
 #include <operators/tnlNeumannBoundaryConditions.h>
-#include <functors/tnlConstantFunction.h>
+#include <functions/tnlConstantFunction.h>
+#include <functions/tnlMeshFunction.h>
 #include <problems/tnlHeatEquationProblem.h>
+#include <mesh/tnlGrid.h>
 
 //typedef tnlDefaultBuildMeshConfig BuildConfig;
 typedef tnlFastBuildConfig BuildConfig;
@@ -43,6 +43,10 @@ class heatEquationConfig
             config.addEntryEnum< tnlString >( "dirichlet" );
             config.addEntryEnum< tnlString >( "neumann" );
 
+         typedef tnlGrid< 1, double, tnlHost, int > Mesh;
+         typedef tnlMeshFunction< Mesh > MeshFunction;
+         tnlDirichletBoundaryConditions< Mesh, MeshFunction >::configSetup( config );
+         tnlDirichletBoundaryConditions< Mesh, tnlConstantFunction< 1 > >::configSetup( config );
          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< double >( "right-hand-side-constant", "This sets a constant value for the right-hand side.", 0.0 );
@@ -79,25 +83,25 @@ class heatEquationSetter
          typedef tnlConstantFunction< Dimensions, Real > ConstantFunction;
          if( boundaryConditionsType == "dirichlet" )
          {
-            typedef tnlAnalyticDirichletBoundaryConditions< MeshType, ConstantFunction, Real, Index > BoundaryConditions;
+            typedef tnlDirichletBoundaryConditions< MeshType, ConstantFunction > BoundaryConditions;
             typedef tnlHeatEquationProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
             SolverStarter solverStarter;
             return solverStarter.template run< Problem >( parameters );
          }
-         typedef tnlAnalyticNeumannBoundaryConditions< MeshType, ConstantFunction, Real, Index > BoundaryConditions;
+         typedef tnlNeumannBoundaryConditions< MeshType, ConstantFunction, Real, Index > BoundaryConditions;
          typedef tnlHeatEquationProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
          SolverStarter solverStarter;
          return solverStarter.template run< Problem >( parameters );
       }
-      typedef tnlVector< Real, Device, Index > VectorType;
+      typedef tnlMeshFunction< MeshType > MeshFunction;
       if( boundaryConditionsType == "dirichlet" )
       {
-         typedef tnlDirichletBoundaryConditions< MeshType, VectorType, Real, Index > BoundaryConditions;
+         typedef tnlDirichletBoundaryConditions< MeshType, MeshFunction > BoundaryConditions;
          typedef tnlHeatEquationProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
          SolverStarter solverStarter;
          return solverStarter.template run< Problem >( parameters );
       }
-      typedef tnlNeumannBoundaryConditions< MeshType, VectorType, Real, Index > BoundaryConditions;
+      typedef tnlNeumannBoundaryConditions< MeshType, MeshFunction, Real, Index > BoundaryConditions;
       typedef tnlHeatEquationProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
       SolverStarter solverStarter;
       return solverStarter.template run< Problem >( parameters );
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 6851150184..96884212e8 100755
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,5 +1,5 @@
 INCLUDE_DIRECTORIES( config )
-ADD_SUBDIRECTORY( functors )
+ADD_SUBDIRECTORY( functions )
 ADD_SUBDIRECTORY( config )
 ADD_SUBDIRECTORY( core )
 ADD_SUBDIRECTORY( debug )
@@ -10,7 +10,7 @@ ADD_SUBDIRECTORY( problems )
 ADD_SUBDIRECTORY( solvers )
 ADD_SUBDIRECTORY( legacy )
 
-set( tnl_SOURCES ${tnl_functors_SOURCES}
+set( tnl_SOURCES ${tnl_functions_SOURCES}
                  ${tnl_config_SOURCES}
                  ${tnl_core_SOURCES}
                  ${tnl_legacy_SOURCES}
@@ -21,7 +21,7 @@ set( tnl_SOURCES ${tnl_functors_SOURCES}
                  ${tnl_problems_SOURCES}
                   )
 
-set( tnl_CUDA__SOURCES ${tnl_functors_CUDA__SOURCES}
+set( tnl_CUDA__SOURCES ${tnl_functions_CUDA__SOURCES}
                        ${tnl_config_CUDA__SOURCES}
                        ${tnl_core_CUDA__SOURCES}
                        ${tnl_legacy_CUDA__SOURCES}
diff --git a/src/core/vectors/tnlSharedVector.h b/src/core/vectors/tnlSharedVector.h
index d68f4c8a5d..25eeed52b4 100644
--- a/src/core/vectors/tnlSharedVector.h
+++ b/src/core/vectors/tnlSharedVector.h
@@ -20,7 +20,7 @@
 
 #include <core/arrays/tnlSharedArray.h>
 #include <core/vectors/tnlVector.h>
-#include <functors/tnlFunction.h>
+#include <functions/tnlFunction.h>
 
 class tnlHost;
 
@@ -37,9 +37,6 @@ class tnlSharedVector : public tnlSharedArray< Real, Device, Index >
    typedef tnlSharedVector< Real, tnlHost, Index > HostType;
    typedef tnlSharedVector< Real, tnlCuda, Index > CudaType;
 
-   //static constexpr tnlFunctionType getFunctionType() { return tnlDiscreteFunction; }
-   enum { functionType = tnlDiscreteFunction };
-
    tnlSharedVector();
 
    tnlSharedVector( Real* data,
diff --git a/src/core/vectors/tnlVector.h b/src/core/vectors/tnlVector.h
index 5108b18355..84b8a8159a 100644
--- a/src/core/vectors/tnlVector.h
+++ b/src/core/vectors/tnlVector.h
@@ -19,7 +19,7 @@
 #define TNLVECTOR_H_
 
 #include <core/arrays/tnlArray.h>
-#include <functors/tnlFunction.h>
+#include <functions/tnlFunction.h>
 
 class tnlHost;
 
@@ -36,9 +36,6 @@ class tnlVector : public tnlArray< Real, Device, Index >
    typedef tnlVector< Real, tnlHost, Index > HostType;
    typedef tnlVector< Real, tnlCuda, Index > CudaType;
 
-   //static constexpr tnlFunctionType getFunctionType() { return tnlDiscreteFunction; }
-   enum { functionType = tnlDiscreteFunction };
-
    tnlVector();
 
    tnlVector( const Index size );
diff --git a/src/functors/CMakeLists.txt b/src/functions/CMakeLists.txt
similarity index 84%
rename from src/functors/CMakeLists.txt
rename to src/functions/CMakeLists.txt
index cb2b36f11f..16b7e41311 100755
--- a/src/functors/CMakeLists.txt
+++ b/src/functions/CMakeLists.txt
@@ -2,7 +2,6 @@ SET( headers tnlFunctionDiscretizer.h
              tnlFunctionDiscretizer_impl.h
              tnlFunctionEnumerator.h
              tnlFunctionEnumerator_impl.h
-             tnlFunctorAdapter.h
              tnlFunctionAdapter.h
              tnlConstantFunction.h
              tnlConstantFunction_impl.h
@@ -16,19 +15,19 @@ SET( headers tnlFunctionDiscretizer.h
              tnlFunction.h
              tnlTestFunction_impl.h )
 
-SET( CURRENT_DIR ${CMAKE_SOURCE_DIR}/src/functors )
+SET( CURRENT_DIR ${CMAKE_SOURCE_DIR}/src/functions )
 set( common_SOURCES
      ${CURRENT_DIR}/tnlTestFunction_impl.cpp )       
 
 IF( BUILD_CUDA )
-   set( tnl_functors_CUDA__SOURCES
+   set( tnl_functions_CUDA__SOURCES
         ${common_SOURCES} 
         ${CURRENT_DIR}/tnlTestFunction_impl.cu
         PARENT_SCOPE )
 ENDIF()    
 
-set( tnl_functors_SOURCES     
+set( tnl_functions_SOURCES     
      ${common_SOURCES}
      PARENT_SCOPE )
         
-INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/functors )
\ No newline at end of file
+INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/functions )
\ No newline at end of file
diff --git a/src/functors/tnlConstantFunction.h b/src/functions/tnlConstantFunction.h
similarity index 94%
rename from src/functors/tnlConstantFunction.h
rename to src/functions/tnlConstantFunction.h
index 1d8480b233..c023e75531 100644
--- a/src/functors/tnlConstantFunction.h
+++ b/src/functions/tnlConstantFunction.h
@@ -20,11 +20,11 @@
 
 #include <iostream>
 #include <core/vectors/tnlStaticVector.h>
-#include <functors/tnlFunction.h>
+#include <functions/tnlFunction.h>
 
 template< int dimensions,
           typename Real = double >
-class tnlConstantFunction : public tnlFunction< dimensions, tnlAnalyticConstantFunction >
+class tnlConstantFunction : public tnlFunction< dimensions, AnalyticConstantFunction >
 {
    public:
       
@@ -85,6 +85,6 @@ std::ostream& operator << ( std::ostream& str, const tnlConstantFunction< dimens
    return str;
 }
 
-#include <functors/tnlConstantFunction_impl.h>
+#include <functions/tnlConstantFunction_impl.h>
 
 #endif /* TNLCONSTANTFUNCTION_H_ */
diff --git a/src/functors/tnlConstantFunction_impl.h b/src/functions/tnlConstantFunction_impl.h
similarity index 100%
rename from src/functors/tnlConstantFunction_impl.h
rename to src/functions/tnlConstantFunction_impl.h
diff --git a/src/functors/tnlExpBumpFunction.h b/src/functions/tnlExpBumpFunction.h
similarity index 96%
rename from src/functors/tnlExpBumpFunction.h
rename to src/functions/tnlExpBumpFunction.h
index 3387344b02..c6f3ffa851 100644
--- a/src/functors/tnlExpBumpFunction.h
+++ b/src/functions/tnlExpBumpFunction.h
@@ -20,11 +20,11 @@
 
 #include <config/tnlParameterContainer.h>
 #include <core/vectors/tnlStaticVector.h>
-#include <functors/tnlFunction.h>
+#include <functions/tnlFunction.h>
 
 template< int dimensions,
           typename Real >
-class tnlExpBumpFunctionBase : public tnlFunction< dimensions, tnlAnalyticFunction >
+class tnlExpBumpFunctionBase : public tnlFunction< dimensions, AnalyticFunction >
 {
    public:
      
@@ -137,7 +137,7 @@ ostream& operator << ( ostream& str, const tnlExpBumpFunction< Dimensions, Real
    return str;
 }
 
-#include <functors/tnlExpBumpFunction_impl.h>
+#include <functions/tnlExpBumpFunction_impl.h>
 
 
 #endif /* TNLEXPBUMPFUNCTION_H_ */
diff --git a/src/functors/tnlExpBumpFunction_impl.h b/src/functions/tnlExpBumpFunction_impl.h
similarity index 99%
rename from src/functors/tnlExpBumpFunction_impl.h
rename to src/functions/tnlExpBumpFunction_impl.h
index 914eac9a9b..3990902b86 100644
--- a/src/functors/tnlExpBumpFunction_impl.h
+++ b/src/functions/tnlExpBumpFunction_impl.h
@@ -18,7 +18,7 @@
 #ifndef TNLEXPBUMPFUNCTION_IMPL_H_
 #define TNLEXPBUMPFUNCTION_IMPL_H_
 
-#include <functors/tnlExpBumpFunction.h>
+#include <functions/tnlExpBumpFunction.h>
 
 template< int dimensions, typename Real >
 bool
diff --git a/src/functors/tnlFunction.h b/src/functions/tnlFunction.h
similarity index 83%
rename from src/functors/tnlFunction.h
rename to src/functions/tnlFunction.h
index da56210367..4ddaf58ac6 100644
--- a/src/functors/tnlFunction.h
+++ b/src/functions/tnlFunction.h
@@ -21,23 +21,21 @@
 
 #include <core/vectors/tnlStaticVector.h>
 
-enum tnlFunctionType { tnlGeneralFunction, 
-                       tnlDiscreteFunction,
-                       tnlAnalyticFunction,
-                       tnlAnalyticConstantFunction };
+enum tnlFunctionType { GeneralFunction, 
+                       MeshFunction,
+                       AnalyticFunction,
+                       AnalyticConstantFunction };
 
 template< int Dimensions,
-          tnlFunctionType FunctionType = tnlGeneralFunction >
+          tnlFunctionType FunctionType = GeneralFunction >
 class tnlFunction
 {
    public:
       
       static const int dimensions = Dimensions;
-      // TODO: restore constexpr when CUDA allows it
-      //static constexpr int getDimensions() { return Dimensions; }
+      static constexpr int getDimensions() { return Dimensions; }
       
-      //static constexpr tnlFunctionType getFunctionType() { return FunctionType; }
-      enum { functionType = FunctionType };
+      static constexpr tnlFunctionType getFunctionType() { return FunctionType; }
 };
 
 
diff --git a/src/functors/tnlFunctionAdapter.h b/src/functions/tnlFunctionAdapter.h
similarity index 55%
rename from src/functors/tnlFunctionAdapter.h
rename to src/functions/tnlFunctionAdapter.h
index b531c23fed..351c9f58fa 100644
--- a/src/functors/tnlFunctionAdapter.h
+++ b/src/functions/tnlFunctionAdapter.h
@@ -30,17 +30,17 @@
 template< typename Mesh,
           typename Function,
           //tnlFunctionType functionType = Function::functionType >
-          int functionType = Function::functionType >
+          int functionType = Function::getFunctionType() >
 class tnlFunctionAdapter
 {
 };
 
 /***
- * Specialization for analytic functions
+ * Specialization for general functions
  */
 template< typename Mesh,
           typename Function >
-class tnlFunctionAdapter< Mesh, Function, tnlAnalyticFunction >
+class tnlFunctionAdapter< Mesh, Function, GeneralFunction >
 {
    public:
       
@@ -50,74 +50,86 @@ class tnlFunctionAdapter< Mesh, Function, tnlAnalyticFunction >
       typedef typename MeshType::IndexType     IndexType;      
       typedef typename FunctionType::VertexType VertexType;
       
-      template< int EntityDimensions >
+      template< typename EntityType >
+      __cuda_callable__ inline
+      static RealType getValue( const FunctionType& function,
+                                const EntityType& meshEntity,
+                                const RealType& time )
+      {         
+         return function.getValue( meshEntity, time );
+      }
+};
+
+/***
+ * Specialization for mesh functions
+ */
+template< typename Mesh,
+          typename Function >
+class tnlFunctionAdapter< Mesh, Function, MeshFunction >
+{
+   public:
+      
+      typedef Function FunctionType;
+      typedef Mesh MeshType;
+      typedef typename FunctionType::RealType  RealType;
+      typedef typename MeshType::IndexType     IndexType;      
+      
+      template< typename EntityType >
       __cuda_callable__ inline
-      static RealType getValue( const MeshType& mesh,
-                                const FunctionType& function,
-                                const IndexType meshEntityIndex,
+      static RealType getValue( const FunctionType& function,
+                                const EntityType& meshEntity,
                                 const RealType& time )
       {         
-         return function.getValue( mesh.template getEntityCenter< EntityDimensions >( meshEntityIndex ), time );
+         return function( meshEntity );
       }
 };
 
 /***
- * Specialization for analytic functions and grids.
+ * Specialization for analytic functions
  */
-template< int Dimensions,
-          typename Real,
-          typename Device,
-          typename Index,
+template< typename Mesh,
           typename Function >
-class tnlFunctionAdapter< tnlGrid< Dimensions, Real, Device, Index >, Function, tnlAnalyticFunction >
+class tnlFunctionAdapter< Mesh, Function, AnalyticFunction >
 {
    public:
-     
-      typedef Function FunctionType; 
-      typedef tnlGrid< Dimensions, Real, Device, Index > MeshType;
-      typedef typename FunctionType::RealType RealType;
-      typedef typename MeshType::IndexType IndexType;
-      typedef typename MeshType::CoordinatesType CoordinatesType;
+      
+      typedef Function FunctionType;
+      typedef Mesh MeshType;
+      typedef typename FunctionType::RealType  RealType;
+      typedef typename MeshType::IndexType     IndexType;      
+      typedef typename FunctionType::VertexType VertexType;
       
       template< typename EntityType >
       __cuda_callable__ inline
-      static RealType getValue( const MeshType& mesh,
-                                const FunctionType& function,
-                                const IndexType meshEntytiIndex,
-                                const EntityType& entity,
+      static RealType getValue( const FunctionType& function,
+                                const EntityType& meshEntity,
                                 const RealType& time )
-      {
-         return function.getValue( entity.getCenter(), time );
-         //return 0.0;
+      {         
+         return function.getValue( meshEntity.getCenter(), time );
       }
 };
 
 /***
- * Specialization for constant functions and grids.
+ * Specialization for constant analytic functions
  */
-template< int Dimensions,
-          typename Real,
-          typename Device,
-          typename Index,
+template< typename Mesh,
           typename Function >
-class tnlFunctionAdapter< tnlGrid< Dimensions, Real, Device, Index >, Function, tnlAnalyticConstantFunction >
+class tnlFunctionAdapter< Mesh, Function, AnalyticConstantFunction >
 {
    public:
-     
-      typedef Function FunctionType; 
-      typedef tnlGrid< Dimensions, Real, Device, Index > MeshType;
-      typedef typename FunctionType::RealType RealType;
-      typedef typename MeshType::IndexType IndexType;
-      typedef typename MeshType::CoordinatesType CoordinatesType;
+      
+      typedef Function FunctionType;
+      typedef Mesh MeshType;
+      typedef typename FunctionType::RealType  RealType;
+      typedef typename MeshType::IndexType     IndexType;      
+      typedef typename FunctionType::VertexType VertexType;
       
       template< typename EntityType >
       __cuda_callable__ inline
-      static RealType getValue( const MeshType& mesh,
-                                const FunctionType& function,
-                                const IndexType meshEntytiIndex,
-                                const EntityType& entity,
+      static RealType getValue( const FunctionType& function,
+                                const EntityType& meshEntity,
                                 const RealType& time )
-      {
+      {         
          return function.getValue( time );
       }
 };
diff --git a/src/functors/tnlFunctionDiscretizer.h b/src/functions/tnlFunctionDiscretizer.h
similarity index 97%
rename from src/functors/tnlFunctionDiscretizer.h
rename to src/functions/tnlFunctionDiscretizer.h
index d5b7ed757a..bceebb6463 100644
--- a/src/functors/tnlFunctionDiscretizer.h
+++ b/src/functions/tnlFunctionDiscretizer.h
@@ -48,6 +48,6 @@ class tnlFunctionDiscretizer
    
 };
 
-#include <functors/tnlFunctionDiscretizer_impl.h>
+#include <functions/tnlFunctionDiscretizer_impl.h>
 
 #endif /* TNLFUNCTIONDISCRETIZER_H_ */
diff --git a/src/functors/tnlFunctionDiscretizer_impl.h b/src/functions/tnlFunctionDiscretizer_impl.h
similarity index 100%
rename from src/functors/tnlFunctionDiscretizer_impl.h
rename to src/functions/tnlFunctionDiscretizer_impl.h
diff --git a/src/functors/tnlFunctionEnumerator.h b/src/functions/tnlFunctionEnumerator.h
similarity index 96%
rename from src/functors/tnlFunctionEnumerator.h
rename to src/functions/tnlFunctionEnumerator.h
index 2e25baf478..57e2e39b5a 100644
--- a/src/functors/tnlFunctionEnumerator.h
+++ b/src/functions/tnlFunctionEnumerator.h
@@ -17,7 +17,7 @@
 #ifndef SRC_FUNCTIONS_TNLFUNCTIONENUMERATOR_H_
 #define SRC_FUNCTIONS_TNLFUNCTIONENUMERATOR_H_
 
-#include <functors/tnlFunctorAdapter.h>
+#include <functions/tnlFunctionAdapter.h>
 
 template< typename Function,
           typename DofVector >
@@ -85,7 +85,7 @@ class tnlFunctionEnumerator
                                        TraverserUserData& userData,
                                        const IndexType index )
             {
-               typedef tnlFunctorAdapter< MeshType, Function > FunctionAdapter;
+               typedef tnlFunctionAdapter< MeshType, Function > FunctionAdapter;
                if( ! *userData.dofVectorCoefficient  )
                   ( *userData.u )[ index ] =
                      ( *userData.functionCoefficient ) * FunctionAdapter::getValue( mesh,
@@ -148,7 +148,7 @@ class tnlFunctionEnumerator< tnlGrid< Dimensions, Real, Device, Index >,
                                      const CoordinatesType& coordinates )
             {
                //printf( "Enumerator::processCell mesh =%p \n", &mesh );
-               typedef tnlFunctorAdapter< MeshType, Function > FunctionAdapter;
+               typedef tnlFunctionAdapter< MeshType, Function > FunctionAdapter;
                if( ! ( *userData.dofVectorCoefficient ) )
                   ( *userData.u )[ index ] =
                      ( *userData.functionCoefficient ) * FunctionAdapter::getValue( mesh,
@@ -175,7 +175,7 @@ class tnlFunctionEnumerator< tnlGrid< Dimensions, Real, Device, Index >,
                                      const IndexType index,
                                      const CoordinatesType& coordinates )
             {
-               typedef tnlFunctorAdapter< MeshType, Function > FunctionAdapter;
+               typedef tnlFunctionAdapter< MeshType, Function > FunctionAdapter;
                if( ! ( *userData.dofVectorCoefficient ) )
                   ( *userData.u )[ index ] =
                      ( *userData.functionCoefficient ) * FunctionAdapter::getValue( mesh,
@@ -196,7 +196,7 @@ class tnlFunctionEnumerator< tnlGrid< Dimensions, Real, Device, Index >,
 
 };
 
-#include <functors/tnlFunctionEnumerator_impl.h>
+#include <functions/tnlFunctionEnumerator_impl.h>
 
 
 
diff --git a/src/functors/tnlFunctionEnumerator_impl.h b/src/functions/tnlFunctionEnumerator_impl.h
similarity index 99%
rename from src/functors/tnlFunctionEnumerator_impl.h
rename to src/functions/tnlFunctionEnumerator_impl.h
index 4e3d824b9d..fb15ab2253 100644
--- a/src/functors/tnlFunctionEnumerator_impl.h
+++ b/src/functions/tnlFunctionEnumerator_impl.h
@@ -17,7 +17,7 @@
 #ifndef SRC_FUNCTIONS_TNLFUNCTIONENUMERATOR_IMPL_H_
 #define SRC_FUNCTIONS_TNLFUNCTIONENUMERATOR_IMPL_H_
 
-#include <functors/tnlFunctionEnumerator.h>
+#include <functions/tnlFunctionEnumerator.h>
 #include <mesh/tnlTraverser_Grid1D.h>
 #include <mesh/tnlTraverser_Grid2D.h>
 #include <mesh/tnlTraverser_Grid3D.h>
diff --git a/src/functions/tnlMeshFunction.h b/src/functions/tnlMeshFunction.h
new file mode 100644
index 0000000000..382a6b6bf0
--- /dev/null
+++ b/src/functions/tnlMeshFunction.h
@@ -0,0 +1,91 @@
+/***************************************************************************
+                          tnlMeshFunction.h  -  description
+                             -------------------
+    begin                : Nov 8, 2015
+    copyright            : (C) 2015 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 <functions/tnlFunction.h>
+
+#ifndef TNLMESHFUNCTION_H
+#define TNLMESHFUNCTION_H
+
+template< typename Mesh,
+          int MeshEntityDimensions = Mesh::meshDimensions,
+          typename Real = typename Mesh::RealType >
+class tnlMeshFunction : public tnlFunction< Mesh::meshDimensions,
+                                            MeshFunction >
+{
+   //static_assert( Mesh::DeviceType::DeviceType == Vector::DeviceType::DeviceType,
+   //               "Both mesh and vector of a mesh function must reside on the same device.");
+   public:
+      
+      typedef Mesh MeshType;      
+      typedef typename MeshType::DeviceType DeviceType;
+      typedef typename MeshType::IndexType IndexType;
+      typedef Real RealType;
+      typedef tnlVector< RealType, DeviceType, IndexType > VectorType;
+      
+      static constexpr int getMeshEntityDimensions() { return  MeshEntityDimensions; };
+      
+      tnlMeshFunction();
+      
+      template< typename Vector >
+      tnlMeshFunction( const MeshType& mesh,
+                       Vector& data,
+                       const IndexType& offset = 0 );
+      
+      static void configSetup( tnlConfigDescription& config,
+                               const tnlString& prefix = "" );
+
+      bool setup( const tnlParameterContainer& parameters,
+                  const tnlString& prefix = "" );      
+      
+      // TODO: implement bind tnlVector ( using shared pointers )
+      /*template< typename Vector >
+      bool bind( const MeshType& mesh,
+                 Vector& data,
+                 const IndexType& offset = 0 );*/
+      
+      void setMesh( const MeshType& mesh ) const;      
+      
+      const MeshType& getMesh() const;      
+      
+      const VectorType& getData() const;      
+      
+      VectorType& getData();      
+      
+      template< typename EntityType >
+      RealType getValue( const EntityType& meshEntity ) const;
+      
+      template< typename EntityType >
+      void setValue( const EntityType& meshEntity,
+                     const RealType& value );
+      
+      template< typename EntityType >
+      RealType& operator()( const EntityType& meshEntityIndex );
+      
+      template< typename EntityType >
+      const RealType& operator()( const EntityType& meshEntityIndex ) const;
+            
+   protected:
+      
+      const MeshType* mesh;
+      
+      VectorType data;     
+};
+
+#include <functions/tnlMeshFunction_impl.h>
+
+#endif	/* TNLMESHFUNCTION_H */
+
diff --git a/src/functions/tnlMeshFunction_impl.h b/src/functions/tnlMeshFunction_impl.h
new file mode 100644
index 0000000000..fc0d396cc9
--- /dev/null
+++ b/src/functions/tnlMeshFunction_impl.h
@@ -0,0 +1,176 @@
+/***************************************************************************
+                          tnlMeshFunction_impl.h  -  description
+                             -------------------
+    begin                : Nov 8, 2015
+    copyright            : (C) 2015 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 <core/tnlAssert.h>
+
+#ifndef TNLMESHFUNCTION_IMPL_H
+#define	TNLMESHFUNCTION_IMPL_H
+
+template< typename Mesh,
+          int MeshEntityDimensions,
+          typename Real >
+tnlMeshFunction< Mesh, MeshEntityDimensions, Real >::
+tnlMeshFunction()
+: mesh( 0 )
+{
+}
+
+template< typename Mesh,
+          int MeshEntityDimensions,
+          typename Real >
+   template< typename Vector >
+tnlMeshFunction< Mesh, MeshEntityDimensions, Real >::
+tnlMeshFunction( const MeshType& mesh,
+                 Vector& data,
+                 const IndexType& offset )
+{
+   //this->bind( mesh, data, offset );   
+}
+
+template< typename Mesh,
+          int MeshEntityDimensions,
+          typename Real >
+void
+tnlMeshFunction< Mesh, MeshEntityDimensions, Real >::
+configSetup( tnlConfigDescription& config,
+             const tnlString& prefix )
+{
+   config.addEntry< tnlString >( prefix + "file", "Dataset for the mesh function." );   
+}
+
+template< typename Mesh,
+          int MeshEntityDimensions,
+          typename Real >
+bool
+tnlMeshFunction< Mesh, MeshEntityDimensions, Real >::
+setup( const tnlParameterContainer& parameters,
+       const tnlString& prefix )
+{
+   if( parameters.checkParameter( prefix + "file" ) )
+   {
+      tnlString fileName = parameters.getParameter< tnlString >( prefix + "file" );
+      if( ! this->data.load( fileName ) )
+         return false;
+   }
+   return true;
+}
+
+/*template< typename Mesh,
+          int MeshEntityDimensions,
+          typename Real >
+   template< typename Vector >
+bool
+tnlMeshFunction< Mesh, MeshEntityDimensions, Real >::
+bind( const MeshType& mesh,
+      Vector& data,
+      const IndexType& offset )
+{
+   this->mesh = &mesh;
+   return this->data.bind( data, offset, mesh.template getEntitiesCount< MeshEntity >() );      
+}*/
+
+template< typename Mesh,
+          int MeshEntityDimensions,
+          typename Real >
+void
+tnlMeshFunction< Mesh, MeshEntityDimensions, Real >::
+setMesh( const MeshType& mesh ) const
+{
+   this->mesh = &mesh;
+}
+
+template< typename Mesh,
+          int MeshEntityDimensions,
+          typename Real >
+const typename tnlMeshFunction< Mesh, MeshEntityDimensions, Real >::MeshType& 
+tnlMeshFunction< Mesh, MeshEntityDimensions, Real >::
+getMesh() const
+{
+   return this->mesh;
+}
+
+template< typename Mesh,
+          int MeshEntityDimensions,
+          typename Real >
+const typename tnlMeshFunction< Mesh, MeshEntityDimensions, Real >::VectorType& 
+tnlMeshFunction< Mesh, MeshEntityDimensions, Real >::
+getData() const
+{
+   return this->data;
+}
+
+template< typename Mesh,
+          int MeshEntityDimensions,
+          typename Real >
+typename tnlMeshFunction< Mesh, MeshEntityDimensions, Real >::VectorType& 
+tnlMeshFunction< Mesh, MeshEntityDimensions, Real >::
+getData()
+{
+   return this->data;
+}
+
+template< typename Mesh,
+          int MeshEntityDimensions,
+          typename Real >
+   template< typename EntityType >          
+typename tnlMeshFunction< Mesh, MeshEntityDimensions, Real >::RealType
+tnlMeshFunction< Mesh, MeshEntityDimensions, Real >::
+getValue( const EntityType& meshEntity ) const
+{
+   static_assert( EntityType::entityDimensions == MeshEntityDimensions, "Calling with wrong EntityType -- entity dimensions do not match." );
+   return this->data.getValue( meshEntity.getIndex() );
+}
+
+template< typename Mesh,
+          int MeshEntityDimensions,
+          typename Real >
+   template< typename EntityType >
+void 
+tnlMeshFunction< Mesh, MeshEntityDimensions, Real >::
+setValue( const EntityType& meshEntity,
+          const RealType& value )
+{
+   static_assert( EntityType::entityDimensions == MeshEntityDimensions, "Calling with wrong EntityType -- entity dimensions do not match." );
+   this->data.setValue( meshEntity.getIndex(), value );
+}
+
+template< typename Mesh,
+          int MeshEntityDimensions,
+          typename Real >
+   template< typename EntityType >
+typename tnlMeshFunction< Mesh, MeshEntityDimensions, Real >::RealType& 
+tnlMeshFunction< Mesh, MeshEntityDimensions, Real >::
+operator()( const EntityType& meshEntity )
+{
+   static_assert( EntityType::entityDimensions == MeshEntityDimensions, "Calling with wrong EntityType -- entity dimensions do not match." );
+   return this->data[ meshEntity.getIndex() ];
+}
+
+template< typename Mesh,
+          int MeshEntityDimensions,
+          typename Real >
+   template< typename EntityType >
+const typename tnlMeshFunction< Mesh, MeshEntityDimensions, Real >::RealType& 
+tnlMeshFunction< Mesh, MeshEntityDimensions, Real >::
+operator()( const EntityType& meshEntity ) const
+{
+   static_assert( EntityType::entityDimensions == MeshEntityDimensions, "Calling with wrong EntityType -- entity dimensions do not match." );
+   return this->data[ meshEntity.getIndex() ];
+}
+
+#endif	/* TNLMESHFUNCTION_IMPL_H */
+
diff --git a/src/functors/tnlSinBumpsFunction.h b/src/functions/tnlSinBumpsFunction.h
similarity index 97%
rename from src/functors/tnlSinBumpsFunction.h
rename to src/functions/tnlSinBumpsFunction.h
index f307d2d237..4af21e0792 100644
--- a/src/functors/tnlSinBumpsFunction.h
+++ b/src/functions/tnlSinBumpsFunction.h
@@ -20,10 +20,10 @@
 
 #include <config/tnlParameterContainer.h>
 #include <core/vectors/tnlStaticVector.h>
-#include <functors/tnlFunction.h>
+#include <functions/tnlFunction.h>
 
 template< typename Vertex >
-class tnlSinBumpsFunctionBase : public tnlFunction< Vertex::size, tnlAnalyticFunction >
+class tnlSinBumpsFunctionBase : public tnlFunction< Vertex::size, AnalyticFunction >
 {
    public:
       
@@ -148,7 +148,7 @@ ostream& operator << ( ostream& str, const tnlSinBumpsFunction< Dimensions, Real
    return str;
 }
 
-#include <functors/tnlSinBumpsFunction_impl.h>
+#include <functions/tnlSinBumpsFunction_impl.h>
 
 
 #endif /* TNLSINBUMPSFUNCTION_H_ */
diff --git a/src/functors/tnlSinBumpsFunction_impl.h b/src/functions/tnlSinBumpsFunction_impl.h
similarity index 99%
rename from src/functors/tnlSinBumpsFunction_impl.h
rename to src/functions/tnlSinBumpsFunction_impl.h
index e9bd9dcf62..308ff76f97 100644
--- a/src/functors/tnlSinBumpsFunction_impl.h
+++ b/src/functions/tnlSinBumpsFunction_impl.h
@@ -18,7 +18,7 @@
 #ifndef TNLSINBUMPSFUNCTION_IMPL_H_
 #define TNLSINBUMPSFUNCTION_IMPL_H_
 
-#include <functors/tnlSinBumpsFunction.h>
+#include <functions/tnlSinBumpsFunction.h>
 
 template< typename Vertex >
 void tnlSinBumpsFunctionBase< Vertex >::setWaveLength( const Vertex& waveLength )
diff --git a/src/functors/tnlSinWaveFunction.h b/src/functions/tnlSinWaveFunction.h
similarity index 95%
rename from src/functors/tnlSinWaveFunction.h
rename to src/functions/tnlSinWaveFunction.h
index 5a215362d5..adb5bc47a9 100644
--- a/src/functors/tnlSinWaveFunction.h
+++ b/src/functions/tnlSinWaveFunction.h
@@ -20,11 +20,11 @@
 
 #include <config/tnlParameterContainer.h>
 #include <core/vectors/tnlStaticVector.h>
-#include <functors/tnlFunction.h>
+#include <functions/tnlFunction.h>
 
 template< int dimensions,
           typename Real = double >
-class tnlSinWaveFunctionBase : public tnlFunction< dimensions, tnlAnalyticFunction >
+class tnlSinWaveFunctionBase : public tnlFunction< dimensions, AnalyticFunction >
 {
    public:
       
@@ -134,6 +134,6 @@ ostream& operator << ( ostream& str, const tnlSinWaveFunction< Dimensions, Real
    return str;
 }
 
-#include <functors/tnlSinWaveFunction_impl.h>
+#include <functions/tnlSinWaveFunction_impl.h>
 
 #endif /* TNLSINWAVEFUNCTION_H_ */
diff --git a/src/functors/tnlSinWaveFunction_impl.h b/src/functions/tnlSinWaveFunction_impl.h
similarity index 99%
rename from src/functors/tnlSinWaveFunction_impl.h
rename to src/functions/tnlSinWaveFunction_impl.h
index a4f73c49aa..14e7976f2b 100644
--- a/src/functors/tnlSinWaveFunction_impl.h
+++ b/src/functions/tnlSinWaveFunction_impl.h
@@ -18,7 +18,7 @@
 #ifndef TNLSINWAVEFUNCTION_IMPL_H_
 #define TNLSINWAVEFUNCTION_IMPL_H_
 
-#include <functors/tnlSinWaveFunction.h>
+#include <functions/tnlSinWaveFunction.h>
 
 template< int dimensions, typename Real >
 tnlSinWaveFunctionBase< dimensions, Real >::tnlSinWaveFunctionBase()
diff --git a/src/functors/tnlTestFunction.h b/src/functions/tnlTestFunction.h
similarity index 96%
rename from src/functors/tnlTestFunction.h
rename to src/functions/tnlTestFunction.h
index 90359c3ded..6684bd23a5 100644
--- a/src/functors/tnlTestFunction.h
+++ b/src/functions/tnlTestFunction.h
@@ -22,11 +22,12 @@
 #include <core/vectors/tnlStaticVector.h>
 #include <config/tnlConfigDescription.h>
 #include <config/tnlParameterContainer.h>
+#include <functions/tnlFunction.h>
 
 template< int FunctionDimensions,
           typename Real = double,
           typename Device = tnlHost >
-class tnlTestFunction
+class tnlTestFunction : public tnlFunction< FunctionDimensions, AnalyticFunction >
 {
    protected:
 
@@ -140,6 +141,6 @@ ostream& operator << ( ostream& str, const tnlTestFunction< FunctionDimensions,
    return f.print( str );
 }
 
-#include <functors/tnlTestFunction_impl.h>
+#include <functions/tnlTestFunction_impl.h>
 
 #endif /* TNLTESTFUNCTION_H_ */
diff --git a/src/functors/tnlTestFunction_impl.cpp b/src/functions/tnlTestFunction_impl.cpp
similarity index 97%
rename from src/functors/tnlTestFunction_impl.cpp
rename to src/functions/tnlTestFunction_impl.cpp
index 948a11afa1..c8baf09f22 100644
--- a/src/functors/tnlTestFunction_impl.cpp
+++ b/src/functions/tnlTestFunction_impl.cpp
@@ -18,7 +18,7 @@
 
 #ifdef TEMPLATE_EXPLICIT_INSTANTIATION
 
-#include <functors/tnlTestFunction.h>
+#include <functions/tnlTestFunction.h>
 
 template class tnlTestFunction< 1, float, tnlHost >;
 template class tnlTestFunction< 2, float, tnlHost >;
diff --git a/src/functors/tnlTestFunction_impl.cu b/src/functions/tnlTestFunction_impl.cu
similarity index 100%
rename from src/functors/tnlTestFunction_impl.cu
rename to src/functions/tnlTestFunction_impl.cu
diff --git a/src/functors/tnlTestFunction_impl.h b/src/functions/tnlTestFunction_impl.h
similarity index 98%
rename from src/functors/tnlTestFunction_impl.h
rename to src/functions/tnlTestFunction_impl.h
index 04c393d4a5..5609a46a1f 100644
--- a/src/functors/tnlTestFunction_impl.h
+++ b/src/functions/tnlTestFunction_impl.h
@@ -19,10 +19,10 @@
 #define TNLTESTFUNCTION_IMPL_H_
 
 #include <core/tnlCuda.h>
-#include <functors/tnlConstantFunction.h>
-#include <functors/tnlExpBumpFunction.h>
-#include <functors/tnlSinBumpsFunction.h>
-#include <functors/tnlSinWaveFunction.h>
+#include <functions/tnlConstantFunction.h>
+#include <functions/tnlExpBumpFunction.h>
+#include <functions/tnlSinBumpsFunction.h>
+#include <functions/tnlSinWaveFunction.h>
 
 template< int FunctionDimensions,
           typename Real,
diff --git a/src/functors/tnlFunctorAdapter.h b/src/functors/tnlFunctorAdapter.h
deleted file mode 100644
index 4369754f45..0000000000
--- a/src/functors/tnlFunctorAdapter.h
+++ /dev/null
@@ -1,278 +0,0 @@
-/***************************************************************************
-                          tnlFunctorAdapter.h  -  description
-                             -------------------
-    begin                : Nov 28, 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 tnlFunctorAdapter_H_
-#define tnlFunctorAdapter_H_
-
-#include <functors/tnlConstantFunction.h>
-#include <functors/tnlFunction.h>
-
-template< typename Mesh,
-          typename Function,
-          int FunctionType = Function::functionType >
-class tnlFunctorAdapter
-{
-};
-
-/****
- * General implementation:
- * - it passes both mesh entity center and mesh entity index
- */
-template< typename Mesh,
-          typename Function >
-class tnlFunctorAdapter< Mesh, Function, tnlGeneralFunction >
-{
-   public:
-
-      typedef Mesh MeshType;
-      typedef Function FunctionType;
-      typedef typename FunctionType::RealType RealType;
-      typedef typename MeshType::IndexType IndexType;
-
-      template< int MeshEntityDimension >
-      __cuda_callable__ 
-      static RealType getValue( const MeshType& mesh,
-                                const FunctionType& function,
-                                const IndexType index,
-                                const RealType& time = 0.0 )
-      {
-         return function.getValue( mesh, //.template getEntityCenter< MeshEntityDimension >,
-                                   index,
-                                   time );
-      }
-};
-
-/****
- * General implementation with specialization for grid functions
- *  - it takes grid coordinates for faster entity center evaluation
- */
-template< int Dimensions,
-          typename Real,
-          typename Device,
-          typename Index,
-          typename Function >
-class tnlFunctorAdapter< tnlGrid< Dimensions, Real, Device, Index >, Function, tnlGeneralFunction >
-{
-         public:
-
-      typedef tnlGrid< Dimensions, Real, Device, Index > MeshType;
-      typedef Function FunctionType;
-      typedef Real RealType;
-      typedef Index IndexType;
-      typedef typename MeshType::VertexType VertexType;
-      typedef typename MeshType::CoordinatesType CoordinatesType;
-
-      __cuda_callable__
-      static RealType getValue( const MeshType& mesh,
-                                const FunctionType& function,
-                                const IndexType index,
-                                const CoordinatesType& coordinates,
-                                const RealType& time = 0.0 )
-      {
-         return function.getValue( mesh, //.template getCellCenter< VertexType >( coordinates ),
-                                   index,
-                                   time );
-      }
-};
-
-/****
- * Specialization for discrete functions:
- * - it passes only the mesh entity index
- */
-template< typename Mesh,
-          typename Function >
-class tnlFunctorAdapter< Mesh, Function, tnlDiscreteFunction >
-{
-   public:
-
-      typedef Mesh MeshType;
-      typedef Function FunctionType;
-      typedef typename FunctionType::RealType RealType;
-      typedef typename MeshType::IndexType IndexType;
-
-      template< int MeshEntityDimension >
-      __cuda_callable__
-      static RealType getValue( const MeshType& mesh,
-                                const FunctionType& function,
-                                const IndexType index,
-                                const RealType& time = 0.0 )
-      {
-         return function.getValue( index,
-                                   time );
-      }
-};
-
-/****
- * Specialization for discrete functions:
- * - it passes only the mesh entity index
- */
-template< int Dimensions,
-          typename Real,
-          typename Device,
-          typename Index,
-          typename Function >
-class tnlFunctorAdapter< tnlGrid< Dimensions, Real, Device, Index >, Function, tnlDiscreteFunction >
-{
-   public:
-
-      typedef tnlGrid< Dimensions, Real, Device, Index > MeshType;
-      typedef Function FunctionType;
-      typedef typename FunctionType::RealType RealType;
-      typedef typename MeshType::IndexType IndexType;
-      typedef typename MeshType::VertexType VertexType;
-      typedef typename MeshType::CoordinatesType CoordinatesType;
-
-      //template< int MeshEntityDimension >
-      __cuda_callable__
-      static RealType getValue( const MeshType& mesh,
-                                const FunctionType& function,
-                                const IndexType index,
-                                const CoordinatesType& coordinates,
-                                const RealType& time = 0.0 )
-      {
-         return function.getValue( index,
-                                   time );
-      }
-};
-
-
-/****
- * Specialization for analytic functions:
- * - it does not pass the mesh entity index
- */
-template< typename Mesh,
-          typename Function >
-class tnlFunctorAdapter< Mesh, Function, tnlAnalyticFunction >
-{
-   public:
-
-      typedef Mesh MeshType;
-      typedef Function FunctionType;
-      typedef typename FunctionType::RealType RealType;
-      typedef typename MeshType::IndexType IndexType;
-
-      template< int MeshEntityDimension >
-      __cuda_callable__
-      static RealType getValue( const MeshType& mesh,
-                                const FunctionType& function,
-                                const IndexType index,
-                                const RealType& time = 0.0 )
-      {
-         return function.getValue( mesh.template getEntityCenter< MeshEntityDimension >,
-                                   time );
-      }
-};
-
-/****
- * Specialization for analytic grid functions:
- * - it does not pass the mesh entity index
- */
-template< int Dimensions,
-          typename Real,
-          typename Device,
-          typename Index,
-          typename Function >
-class tnlFunctorAdapter< tnlGrid< Dimensions, Real, Device, Index >, Function, tnlAnalyticFunction >
-{
-         public:
-
-      typedef tnlGrid< Dimensions, Real, Device, Index > MeshType;
-      typedef Function FunctionType;
-      typedef Real RealType;
-      typedef Index IndexType;
-      typedef typename MeshType::VertexType VertexType;
-      typedef typename MeshType::CoordinatesType CoordinatesType;
-
-      __cuda_callable__
-      static RealType getValue( const MeshType& mesh,
-                                const FunctionType& function,
-                                const IndexType index,
-                                const CoordinatesType& coordinates,
-                                const RealType& time = 0.0 )
-      {
-         return function.getValue( mesh.template getCellCenter< VertexType >( coordinates ),
-                                   time );
-      }
-};
-
-// TODO: Fix the specializations for the constant function.
-#ifdef UNDEF
-/****
- * Specialization for constant function
- *  - it does not ask the mesh for the mesh entity center
- */
-template< typename Mesh,
-          int FunctionDimensions,
-          typename Real >
-class tnlFunctorAdapter< Mesh, tnlConstantFunction< FunctionDimensions, Real >, tnlAnalyticFunction >
-{
-   public:
-
-      typedef Mesh MeshType;
-      typedef tnlConstantFunction< FunctionDimensions, Real > FunctionType;
-      typedef typename FunctionType::RealType RealType;
-      typedef typename MeshType::IndexType IndexType;
-      typedef typename MeshType::VertexType VertexType;
-
-      template< int MeshEntityDimension >
-      __cuda_callable__
-      static RealType getValue( const MeshType& mesh,
-                                const FunctionType& function,
-                                const IndexType index,
-                                const RealType& time = 0.0 )
-      {
-         VertexType v;
-         return function.getValue( v, time );
-      }
-};
-
-/****
- * Specialization for grids and constant function
- *  - it takes grid coordinates for faster entity center evaluation
- */
-template< int Dimensions,
-          typename Real,
-          typename Device,
-          typename Index >
-class tnlFunctorAdapter< tnlGrid< Dimensions, Real, Device, Index >,
-                          tnlConstantFunction< Dimensions, Real >,
-                          tnlAnalyticFunction >
-{
-   public:
-
-      typedef tnlGrid< Dimensions, Real, Device, Index > MeshType;
-      typedef tnlConstantFunction< Dimensions, Real > FunctionType;
-      typedef Real RealType;
-      typedef Index IndexType;
-      typedef typename MeshType::VertexType VertexType;
-      typedef typename MeshType::CoordinatesType CoordinatesType;
-
-      __cuda_callable__
-      static RealType getValue( const MeshType& mesh,
-                                const FunctionType& function,
-                                const IndexType index,
-                                const CoordinatesType& coordinates,
-                                const RealType& time = 0.0 )
-      {
-         VertexType v;
-         return function.getValue( v, time );
-      }
-};
-
-#endif /* UNDEF */
-
-#endif /* tnlFunctorAdapter_H_ */
diff --git a/src/functors/tnlMeshFunction.h b/src/functors/tnlMeshFunction.h
deleted file mode 100644
index 7aa538a04b..0000000000
--- a/src/functors/tnlMeshFunction.h
+++ /dev/null
@@ -1,72 +0,0 @@
-/***************************************************************************
-                          tnlMeshFunction.h  -  description
-                             -------------------
-    begin                : Nov 8, 2015
-    copyright            : (C) 2015 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 <functors/tnlFunction.h>
-
-#ifndef TNLMESHFUNCTION_H
-#define TNLMESHFUNCTION_H
-
-template< typename Mesh,
-          typename Real = typename Mesh::RealType,
-          int MeshEntitiesDimensions_ = Mesh::Dimensions >
-class tnlMeshFunction : public tnlFunction< Mesh::Dimensions,
-                                            tnlDiscreteFunction >
-{
-   static_assert( Mesh::DeviceType::DeviceType == Vector::DeviceType::DeviceType,
-                  "Both mesh and vector of a mesh function must reside on the same device.");
-   public:
-      
-      typedef Mesh MeshType;      
-      typedef typename MeshType::DeviceType DeviceType;
-      typedef typename MeshType::IndexType IndexType;
-      typedef Real RealType;
-      typedef tnlSharedVector< RealType, DeviceType, IndexType > VectorType;
-      
-      static constexpr int MeshEntitiesDimensions = MeshEntitiesDimensions_;
-      
-      tnlMeshFunction();
-      
-      tnlMeshFunction( const MeshType* mesh,
-                       RealType* data );
-      
-      void bind( const MeshType* mesh,
-                 RealType* data );
-      
-      const MeshType& getMesh() const;
-      
-      const VectorType& getData() const;
-      
-      const RealType& getValue( const IndexType meshEntityIndex );
-      
-      void setValue( const IndexType meshEntityIndex,
-                     const RealType& value );
-      
-      RealType& operator()( const IndexType meshEntityIndex );
-      
-      const RealType& operator()( const IndexType meshEntityIndex ) const;
-      
-   protected:
-      
-      const MeshType* mesh;
-      
-      tnlSharedVector<  > data;
-      
-};
-
-
-#endif	/* TNLMESHFUNCTION_H */
-
diff --git a/src/functors/tnlMeshFunction_impl.h b/src/functors/tnlMeshFunction_impl.h
deleted file mode 100644
index 72af176271..0000000000
--- a/src/functors/tnlMeshFunction_impl.h
+++ /dev/null
@@ -1,115 +0,0 @@
-/***************************************************************************
-                          tnlMeshFunction_impl.h  -  description
-                             -------------------
-    begin                : Nov 8, 2015
-    copyright            : (C) 2015 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 <core/tnlAssert.h>
-
-#ifndef TNLMESHFUNCTION_IMPL_H
-#define	TNLMESHFUNCTION_IMPL_H
-
-template< typename Mesh,
-          typename Vector,
-          int MeshEntitiesDimensions >
-tnlMeshFunction< Mesh, Vector, MeshEntitiesDimensions >::
-tnlMeshFunction()
-: mesh( 0 )
-{
-}
-      
-template< typename Mesh,
-          typename Vector,
-          int MeshEntitiesDimensions >
-tnlMeshFunction< Mesh, Vector, MeshEntitiesDimensions >::
-tnlMeshFunction( const MeshType* mesh )
-{
-   this->setMesh( mesh );
-}
-
-template< typename Mesh,
-          typename Vector,
-          int MeshEntitiesDimensions >
-void 
-tnlMeshFunction< Mesh, Vector, MeshEntitiesDimensions >::
-setMesh( const MeshType* mesh )
-{
-   tnlAssert( mesh != NULL, );
-   this->mesh = mesh;
-   this->data.setSize( this->mesh->template getNumberOfEntities< MeshEntitiesDimensions >() );
-}
-
-template< typename Mesh,
-          typename Vector,
-          int MeshEntitiesDimensions >
-const typename tnlMeshFunction< Mesh, Vector, MeshEntitiesDimensions >::MeshType& 
-tnlMeshFunction< Mesh, Vector, MeshEntitiesDimensions >::
-getMesh() const
-{
-   return * this->mesh;
-}
-
-template< typename Mesh,
-          typename Vector,
-          int MeshEntitiesDimensions >
-const typename tnlMeshFunction< Mesh, Vector, MeshEntitiesDimensions >::VectorType& 
-tnlMeshFunction< Mesh, Vector, MeshEntitiesDimensions >::
-getData() const
-{
-   return this->data;
-}
-
-template< typename Mesh,
-          typename Vector,
-          int MeshEntitiesDimensions >
-const typename tnlMeshFunction< Mesh, Vector, MeshEntitiesDimensions >::RealType& 
-tnlMeshFunction< Mesh, Vector, MeshEntitiesDimensions >::
-getValue( const IndexType meshEntityIndex )
-{
-   return this->data.getValue( meshEntityIndex );
-}
-
-template< typename Mesh,
-          typename Vector,
-          int MeshEntitiesDimensions >
-void 
-tnlMeshFunction< Mesh, Vector, MeshEntitiesDimensions >::
-setValue( const IndexType meshEntityIndex,
-          const RealType& value )
-{
-   this->data.setValue( meshEntityIndex, value );
-}
-
-template< typename Mesh,
-          typename Vector,
-          int MeshEntitiesDimensions >
-typename tnlMeshFunction< Mesh, Vector, MeshEntitiesDimensions >::RealType& 
-tnlMeshFunction< Mesh, Vector, MeshEntitiesDimensions >::
-operator[]( const IndexType meshEntityIndex )
-{
-   return this->data[ meshEntityIndex ];
-}
-
-template< typename Mesh,
-          typename Vector,
-          int MeshEntitiesDimensions >
-const typename tnlMeshFunction< Mesh, Vector, MeshEntitiesDimensions >::RealType& 
-tnlMeshFunction< Mesh, Vector, MeshEntitiesDimensions >::
-operator[]( const IndexType meshEntityIndex ) const
-{
-   return this->data[ meshEntityIndex ];
-}
-
-#endif	/* TNLMESHFUNCTION_IMPL_H */
-
diff --git a/src/matrices/tnlMatrixSetter.h b/src/matrices/tnlMatrixSetter.h
index 3eba3c5675..516c86a54e 100644
--- a/src/matrices/tnlMatrixSetter.h
+++ b/src/matrices/tnlMatrixSetter.h
@@ -130,12 +130,11 @@ class tnlMatrixSetter< tnlGrid< Dimensions, Real, Device, Index >,
          template< typename EntityType >
          __cuda_callable__
          static void processEntity( const MeshType& mesh,
-                                    TraversalUserData& userData,
-                                    const IndexType index,
+                                    TraversalUserData& userData,                                    
                                     const EntityType& entity )
          {
-            ( *userData.rowLengths )[ index ] =
-                     userData.boundaryConditions->getLinearSystemRowLength( mesh, index, entity );
+            ( *userData.rowLengths )[ entity.getIndex() ] =
+                     userData.boundaryConditions->getLinearSystemRowLength( mesh, entity.getIndex(), entity );
          }
 
          
@@ -172,11 +171,10 @@ class tnlMatrixSetter< tnlGrid< Dimensions, Real, Device, Index >,
          __cuda_callable__
          static void processEntity( const MeshType& mesh,
                                     TraversalUserData& userData,
-                                    const IndexType index,
                                     const EntityType& entity )
          {
-            ( *userData.rowLengths )[ index ] =
-                     userData.differentialOperator->getLinearSystemRowLength( mesh, index, entity );
+            ( *userData.rowLengths )[ entity.getIndex() ] =
+                     userData.differentialOperator->getLinearSystemRowLength( mesh, entity.getIndex(), entity );
          }
 
 
diff --git a/src/mesh/grids/tnlTraverser_Grid1D_impl.h b/src/mesh/grids/tnlTraverser_Grid1D_impl.h
index fedb59d610..7989ecf82d 100644
--- a/src/mesh/grids/tnlTraverser_Grid1D_impl.h
+++ b/src/mesh/grids/tnlTraverser_Grid1D_impl.h
@@ -38,10 +38,10 @@ processBoundaryEntities( const GridType& grid,
    const IndexType& xSize = grid.getDimensions().x();   
    coordinates.x() = 0;
    entity.refresh();
-   EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+   EntitiesProcessor::processEntity( grid, userData, entity );
    coordinates.x() = xSize - 1;
    entity.refresh();
-   EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+   EntitiesProcessor::processEntity( grid, userData, entity );
 }
 
 template< typename Real,
@@ -66,7 +66,7 @@ processInteriorEntities( const GridType& grid,
         cell.getCoordinates().x()++ )
    {
       cell.refresh();
-      EntitiesProcessor::processEntity( grid, userData, grid.getEntityIndex( cell ), cell );
+      EntitiesProcessor::processEntity( grid, userData, cell );
    }
 }
 
diff --git a/src/mesh/grids/tnlTraverser_Grid2D_impl.h b/src/mesh/grids/tnlTraverser_Grid2D_impl.h
index 7946e871b3..d1ce44d46d 100644
--- a/src/mesh/grids/tnlTraverser_Grid2D_impl.h
+++ b/src/mesh/grids/tnlTraverser_Grid2D_impl.h
@@ -43,19 +43,19 @@ processBoundaryEntities( const GridType& grid,
    {
       coordinates.y() = 0;
       entity.refresh();
-      EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+      EntitiesProcessor::processEntity( grid, userData, entity );
       coordinates.y() = ySize - 1;
       entity.refresh();
-      EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+      EntitiesProcessor::processEntity( grid, userData, entity );
    }
    for( coordinates.y() = 1; coordinates.y() < ySize - 1; coordinates.y() ++ )
    {
       coordinates.x() = 0;
       entity.refresh();
-      EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+      EntitiesProcessor::processEntity( grid, userData, entity );
       coordinates.x() = xSize - 1;
       entity.refresh();
-      EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+      EntitiesProcessor::processEntity( grid, userData, entity );
    }
 }
 
@@ -87,7 +87,7 @@ processInteriorEntities( const GridType& grid,
       for( coordinates.x() = 1; coordinates.x() < xSize - 1; coordinates.x() ++ )
       {
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
       }
 }
 
diff --git a/src/mesh/grids/tnlTraverser_Grid3D_impl.h b/src/mesh/grids/tnlTraverser_Grid3D_impl.h
index 6103c800a5..8bc9f1c7f3 100644
--- a/src/mesh/grids/tnlTraverser_Grid3D_impl.h
+++ b/src/mesh/grids/tnlTraverser_Grid3D_impl.h
@@ -44,10 +44,10 @@ processBoundaryEntities( const GridType& grid,
       {
          coordinates.z() = 0;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
          coordinates.z() = zSize - 1;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
       }
 
    for( coordinates.z() = 0; coordinates.z() < zSize; coordinates.z() ++ )
@@ -55,10 +55,10 @@ processBoundaryEntities( const GridType& grid,
       {
          coordinates.y() = 0;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
          coordinates.y() = ySize - 1;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
       }
 
    for( coordinates.z() = 0; coordinates.z() < zSize; coordinates.z() ++ )
@@ -66,10 +66,10 @@ processBoundaryEntities( const GridType& grid,
       {
          coordinates.x() = 0;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
          coordinates.x() = xSize - 1;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
       }
 }
 template< typename Real,
@@ -102,7 +102,7 @@ processInteriorEntities( const GridType& grid,
          for( coordinates.x() = 1; coordinates.x() < xSize - 1; coordinates.x() ++ )
          {
             entity.refresh();
-            EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+            EntitiesProcessor::processEntity( grid, userData, entity );
          }
 }
 
@@ -135,10 +135,10 @@ processBoundaryEntities( const GridType& grid,
       {
          coordinates.x() = 0;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
          coordinates.x() = xSize;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
       }
 
    entity.setOrientation( typename EntityType::EntityOrientationType( 0, 1, 0 ) );
@@ -147,10 +147,10 @@ processBoundaryEntities( const GridType& grid,
       {
          coordinates.y() = 0;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
          coordinates.y() = ySize;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
       }
   
    entity.setOrientation( typename EntityType::EntityOrientationType( 0, 0, 1 ) );
@@ -159,10 +159,10 @@ processBoundaryEntities( const GridType& grid,
       {         
          coordinates.z() = 0;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
          coordinates.z() = zSize;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
       }     
 }
 
@@ -199,7 +199,7 @@ processInteriorEntities( const GridType& grid,
          for( coordinates.x() = 1; coordinates.x() < xSize; coordinates.x() ++ )
          {
             entity.refresh();         
-            EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+            EntitiesProcessor::processEntity( grid, userData, entity );
          }
 
    entity.setOrientation( typename EntityType::EntityOrientationType( 0, 1, 0 ) );
@@ -208,7 +208,7 @@ processInteriorEntities( const GridType& grid,
          for( coordinates.x() = 0; coordinates.x() < xSize; coordinates.x() ++ )
          {
             entity.refresh();
-            EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+            EntitiesProcessor::processEntity( grid, userData, entity );
          }
             
    
@@ -218,7 +218,7 @@ processInteriorEntities( const GridType& grid,
          for( coordinates.x() = 0; coordinates.x() < xSize; coordinates.x() ++ )
          {
             entity.refresh();
-            EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+            EntitiesProcessor::processEntity( grid, userData, entity );
          }
 }
 
@@ -252,10 +252,10 @@ processBoundaryEntities( const GridType& grid,
       {
          coordinates.z() = 0;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );        
+         EntitiesProcessor::processEntity( grid, userData, entity );        
          coordinates.z() = zSize;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
       }
    
    entity.setBasis( typename EntityType::EntityBasisType( 0, 1, 0 ) );
@@ -264,10 +264,10 @@ processBoundaryEntities( const GridType& grid,
       {
          coordinates.z() = 0;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
          coordinates.z() = zSize;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
       }
    
    entity.setBasis( typename EntityType::EntityBasisType( 1, 0, 0 ) );
@@ -276,10 +276,10 @@ processBoundaryEntities( const GridType& grid,
       {
          coordinates.y() = 0;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
          coordinates.y() = ySize;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
       }
    
    entity.setBasis( typename EntityType::EntityBasisType( 0, 0, 1 ) );
@@ -288,10 +288,10 @@ processBoundaryEntities( const GridType& grid,
       {
          coordinates.y() = 0;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
          coordinates.y() = ySize;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
       }
 
    entity.setBasis( typename EntityType::EntityBasisType( 0, 1, 0 ) );
@@ -300,10 +300,10 @@ processBoundaryEntities( const GridType& grid,
       {
          coordinates.x() = 0;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
          coordinates.x() = xSize;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
       }
    
    entity.setBasis( typename EntityType::EntityBasisType( 0, 0, 1 ) );
@@ -312,10 +312,10 @@ processBoundaryEntities( const GridType& grid,
       {
          coordinates.x() = 0;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, grid.getEdgeIndex( entity ), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
          coordinates.x() = xSize;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, grid.getEdgeIndex( entity ), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
       }
 }
 
@@ -352,7 +352,7 @@ processInteriorEntities( const GridType& grid,
          for( coordinates.x() = 0; coordinates.x() < xSize; coordinates.x() ++ )
          {
             entity.refresh();
-            EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+            EntitiesProcessor::processEntity( grid, userData, entity );
          }
 
    entity.setBasis( typename EntityType::EntityBasisType( 0, 1, 0 ) );
@@ -361,7 +361,7 @@ processInteriorEntities( const GridType& grid,
          for( coordinates.x() = 1; coordinates.x() < xSize; coordinates.x() ++ )
          {
             entity.refresh();
-            EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+            EntitiesProcessor::processEntity( grid, userData, entity );
          }
    
    entity.setBasis( typename EntityType::EntityBasisType( 1, 0, 0 ) );
@@ -370,7 +370,7 @@ processInteriorEntities( const GridType& grid,
          for( coordinates.x() = 1; coordinates.x() < xSize; coordinates.x() ++ )
          {
             entity.refresh();
-            EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+            EntitiesProcessor::processEntity( grid, userData, entity );
          }
 }
 
@@ -404,10 +404,10 @@ processBoundaryEntities( const GridType& grid,
       {
          coordinates.z() = 0;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
          coordinates.z() = zSize;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
       }
 
    for( coordinates.z() = 0; coordinates.z() <= zSize; coordinates.z() ++ )
@@ -415,10 +415,10 @@ processBoundaryEntities( const GridType& grid,
       {
          coordinates.y() = 0;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
          coordinates.y() = ySize;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
       }
 
    for( coordinates.z() = 0; coordinates.z() <= zSize; coordinates.z() ++ )
@@ -426,10 +426,10 @@ processBoundaryEntities( const GridType& grid,
       {
          coordinates.x() = 0;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
          coordinates.x() = xSize;
          entity.refresh();
-         EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+         EntitiesProcessor::processEntity( grid, userData, entity );
       }
 }
 
@@ -464,7 +464,7 @@ processInteriorEntities( const GridType& grid,
          for( coordinates.x() = 1; coordinates.x() < xSize; coordinates.x() ++ )
          {
             entity.refresh();
-            EntitiesProcessor::processEntity( grid, userData, entity.getIndex(), entity );
+            EntitiesProcessor::processEntity( grid, userData, entity );
          }
 }
 
@@ -511,7 +511,6 @@ __global__ void tnlTraverserGrid3DCells( const tnlGrid< 3, Real, tnlCuda, Index
          EntitiesProcessor::template processEntity
             ( *grid,
               *userData,
-              grid->getEntityIndex( entity ),
               entity );
       }
    }
@@ -535,6 +534,7 @@ __global__ void tnlTraverserGrid3DFaces( const tnlGrid< 3, Real, tnlCuda, Index
    typedef tnlGrid< 3, Real, tnlCuda, Index > GridType;
    const int FaceDimensions = GridType::meshDimensions - 1;
    typename GridType::template GridEntity< FaceDimensions > entity( *grid );
+   entity.setOrientation( nx, ny, nz );
    typedef typename GridType::CoordinatesType CoordinatesType;
    CoordinatesType& coordinates = entity.getCoordinates();
 
@@ -556,9 +556,7 @@ __global__ void tnlTraverserGrid3DFaces( const tnlGrid< 3, Real, tnlCuda, Index
          EntitiesProcessor::processEntity
             ( *grid,
               *userData,
-              grid->geEntityIndex( entity ),
-              entity,
-              nx, ny, nz );
+              entity );
       }
    }
 }
@@ -581,6 +579,7 @@ __global__ void tnlTraverserGrid3DEdges( const tnlGrid< 3, Real, tnlCuda, Index
    typedef tnlGrid< 3, Real, tnlCuda, Index > GridType;
    const int EdgeDimensions = 1;
    typename GridType::template GridEntity< EdgeDimensions > entity( *grid );
+   entity.setBasis( dx, dy, dz );
    typedef typename GridType::CoordinatesType CoordinatesType;
    CoordinatesType& coordinates = entity.getCoordinates();
 
@@ -602,7 +601,6 @@ __global__ void tnlTraverserGrid3DEdges( const tnlGrid< 3, Real, tnlCuda, Index
          EntitiesProcessor::processEntity
             ( *grid,
               *userData,
-              grid->getEntityIndex( entity ),
               entity );
       }
    }
@@ -644,7 +642,6 @@ __global__ void tnlTraverserGrid3DVertices( const tnlGrid< 3, Real, tnlCuda, Ind
          EntitiesProcessor::template processEntity
          ( *grid,
            *userData,
-           grid->getEntityIndex( entity ),
            entity );
       }
    }
diff --git a/src/mesh/tnlDummyMesh.h b/src/mesh/tnlDummyMesh.h
index 9935628b02..4c2d3b2ddd 100644
--- a/src/mesh/tnlDummyMesh.h
+++ b/src/mesh/tnlDummyMesh.h
@@ -18,7 +18,9 @@
 #ifndef TNLDUMMYMESH_H_
 #define TNLDUMMYMESH_H_
 
-template< typename Real, typename Device, typename Index >
+template< typename Real = double,
+          typename Device = tnlHost,
+          typename Index = int >
 class tnlDummyMesh
 {
    public:
@@ -27,6 +29,7 @@ class tnlDummyMesh
    typedef Device DeviceType;
    typedef Index IndexType;
 
+   static const int meshDimensions = 1;
 
    const Real& getParametricStep(){ return 0.0; }
 
diff --git a/src/operators/CMakeLists.txt b/src/operators/CMakeLists.txt
index 6ddc8f02ab..f03203db39 100755
--- a/src/operators/CMakeLists.txt
+++ b/src/operators/CMakeLists.txt
@@ -6,12 +6,8 @@ SET( headers tnlFiniteDifferences.h
              tnlFiniteDifferences_impl.h
              tnlDirichletBoundaryConditions.h
              tnlDirichletBoundaryConditions_impl.h
-             tnlAnalyticDirichletBoundaryConditions.h
-             tnlAnalyticDirichletBoundaryConditions_impl.h
              tnlNeumannBoundaryConditions.h
              tnlNeumannBoundaryConditions_impl.h
-             tnlAnalyticNeumannBoundaryConditions.h
-             tnlAnalyticNeumannBoundaryConditions_impl.h
              tnlExactOperatorEvaluator.h
              tnlOperatorEnumerator.h
              tnlOperatorEnumerator_impl.h )
diff --git a/src/operators/diffusion/tnlExactLinearDiffusion.h b/src/operators/diffusion/tnlExactLinearDiffusion.h
index 006167f71f..2f8665e243 100644
--- a/src/operators/diffusion/tnlExactLinearDiffusion.h
+++ b/src/operators/diffusion/tnlExactLinearDiffusion.h
@@ -18,14 +18,14 @@
 #ifndef TNLEXACTLINEARDIFFUSION_H_
 #define TNLEXACTLINEARDIFFUSION_H_
 
-#include <functors/tnlFunction.h>
+#include <functions/tnlFunction.h>
 
 template< int Dimensions >
 class tnlExactLinearDiffusion
 {};
 
 template<>
-class tnlExactLinearDiffusion< 1 > : public tnlFunction< 1, tnlAnalyticFunction >
+class tnlExactLinearDiffusion< 1 > : public tnlFunction< 1, AnalyticFunction >
 {
    public:
 
@@ -39,7 +39,7 @@ class tnlExactLinearDiffusion< 1 > : public tnlFunction< 1, tnlAnalyticFunction
 };
 
 template<>
-class tnlExactLinearDiffusion< 2 > : public tnlFunction< 2, tnlAnalyticFunction >
+class tnlExactLinearDiffusion< 2 > : public tnlFunction< 2, AnalyticFunction >
 {
    public:
       
diff --git a/src/operators/diffusion/tnlLinearDiffusion.h b/src/operators/diffusion/tnlLinearDiffusion.h
index 2cf89fb77f..c416c92243 100644
--- a/src/operators/diffusion/tnlLinearDiffusion.h
+++ b/src/operators/diffusion/tnlLinearDiffusion.h
@@ -52,7 +52,6 @@ class tnlLinearDiffusion< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index
       template< typename Vector >
       __cuda_callable__
       inline Real getValue( const MeshType& mesh,
-                            const IndexType cellIndex,
                             const CellType& cell,
                             const Vector& u,
                             const RealType& time ) const;
@@ -98,7 +97,6 @@ class tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index
                 typename EntityType >
       __cuda_callable__
       inline Real getValue( const MeshType& mesh,
-                            const IndexType cellIndex,
                             const EntityType& entity,
                             const Vector& u,
                             const Real& time ) const;
@@ -146,7 +144,6 @@ class tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index
                 typename EntityType >
       __cuda_callable__
       inline Real getValue( const MeshType& mesh,
-                            const IndexType cellIndex,
                             const EntityType& entity,
                             const Vector& u,
                             const Real& time ) const;
diff --git a/src/operators/diffusion/tnlLinearDiffusion_impl.h b/src/operators/diffusion/tnlLinearDiffusion_impl.h
index 86e5f1675c..298e1b01df 100644
--- a/src/operators/diffusion/tnlLinearDiffusion_impl.h
+++ b/src/operators/diffusion/tnlLinearDiffusion_impl.h
@@ -47,7 +47,6 @@ inline
 Real
 tnlLinearDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
 getValue( const MeshType& mesh,
-          const IndexType cellIndex,
           const CellType& cell,
           const Vector& u,
           const Real& time ) const
@@ -55,7 +54,7 @@ getValue( const MeshType& mesh,
    auto neighbourEntities = cell.getNeighbourEntities();
    const RealType& hxSquareInverse = mesh.template getSpaceStepsProducts< - 2 >();
    return ( u[ neighbourEntities.template getEntityIndex< -1 >() ]
-            - 2.0 * u[ cellIndex ]
+            - 2.0 * u[ cell.getIndex() ]
             + u[ neighbourEntities.template getEntityIndex< 1 >() ] ) * hxSquareInverse;
 }
 
@@ -147,7 +146,6 @@ inline
 Real
 tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
 getValue( const MeshType& mesh,
-          const IndexType cellIndex,
           const EntityType& entity,
           const Vector& u,
           const Real& time ) const
@@ -155,12 +153,11 @@ getValue( const MeshType& mesh,
    auto neighbourEntities = entity.getNeighbourEntities();
    const RealType& hxSquareInverse = mesh.template getSpaceStepsProducts< -2, 0 >();
    const RealType& hySquareInverse = mesh.template getSpaceStepsProducts< 0, -2 >();
-   return ( u[ neighbourEntities.template getEntityIndex< -1, 0 >() ]
-            - 2.0 * u[ cellIndex ]
-            + u[ neighbourEntities.template getEntityIndex< 1, 0 >() ] ) * hxSquareInverse +
-           ( u[ neighbourEntities.template getEntityIndex< 0, -1 >() ]
-             - 2.0 * u[ cellIndex ]
-             + u[ neighbourEntities.template getEntityIndex< 0, 1 >() ] ) * hySquareInverse;
+   return ( u[ neighbourEntities.template getEntityIndex< -1,  0 >() ]
+          + u[ neighbourEntities.template getEntityIndex<  1,  0 >() ] ) * hxSquareInverse +
+          ( u[ neighbourEntities.template getEntityIndex<  0, -1 >() ]
+          + u[ neighbourEntities.template getEntityIndex<  0,  1 >() ] ) * hySquareInverse
+          - 2.0 * u[ entity.getIndex() ] * ( hxSquareInverse + hySquareInverse );
 }
 
 template< typename MeshReal,
@@ -223,24 +220,21 @@ inline
 Real
 tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
 getValue( const MeshType& mesh,
-          const IndexType cellIndex,
           const EntityType& entity,
           const Vector& u,
           const Real& time ) const
 {
    auto neighbourEntities = entity.getNeighbourEntities();
-   const RealType& hxSquareInverse = mesh.template getSpaceStepsProducts< -2, 0, 0 >();
-   const RealType& hySquareInverse = mesh.template getSpaceStepsProducts< 0, -2, 0 >();
-   const RealType& hzSquareInverse = mesh.template getSpaceStepsProducts< 0, 0, -2 >();
-   return (   u[ neighbourEntities.template getEntityIndex< -1, 0, 0 >() ]
-            - 2.0 * u[ cellIndex ]
-            + u[ neighbourEntities.template getEntityIndex< 1, 0, 0 >() ] ) * hxSquareInverse +
-          (   u[ neighbourEntities.template getEntityIndex< 0, -1, 0 >() ]
-            - 2.0 * u[ cellIndex ]
-            + u[ neighbourEntities.template getEntityIndex< 0, 1, 0 >() ] ) * hySquareInverse +
-          (   u[ neighbourEntities.template getEntityIndex< 0, 0, -1 >() ]
-            - 2.0 * u[ cellIndex ]
-            + u[ neighbourEntities.template getEntityIndex< 0, 0, 1 >() ] ) * hzSquareInverse;
+   const RealType& hxSquareInverse = mesh.template getSpaceStepsProducts< -2,  0,  0 >();
+   const RealType& hySquareInverse = mesh.template getSpaceStepsProducts<  0, -2,  0 >();
+   const RealType& hzSquareInverse = mesh.template getSpaceStepsProducts<  0,  0, -2 >();
+   return (   u[ neighbourEntities.template getEntityIndex< -1,  0,  0 >() ]
+            + u[ neighbourEntities.template getEntityIndex<  1,  0,  0 >() ] ) * hxSquareInverse +
+          (   u[ neighbourEntities.template getEntityIndex<  0, -1,  0 >() ]
+            + u[ neighbourEntities.template getEntityIndex<  0,  1,  0 >() ] ) * hySquareInverse +
+          (   u[ neighbourEntities.template getEntityIndex<  0,  0, -1 >() ]
+            + u[ neighbourEntities.template getEntityIndex<  0,  0,  1 >() ] ) * hzSquareInverse
+         - 2.0 * u[ entity.getIndex() ] * ( hxSquareInverse + hySquareInverse + hzSquareInverse );
 }
 
 template< typename MeshReal,
diff --git a/src/operators/tnlAnalyticDirichletBoundaryConditions.h b/src/operators/tnlAnalyticDirichletBoundaryConditions.h
deleted file mode 100644
index 5564c0cb88..0000000000
--- a/src/operators/tnlAnalyticDirichletBoundaryConditions.h
+++ /dev/null
@@ -1,114 +0,0 @@
-/***************************************************************************
-           tnlAnalyticDirichletBoundaryConditions.h  -  description
-                             -------------------
-    begin                : Nov 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 tnlAnalyticDirichletBoundaryConditions_H
-#define	tnlAnalyticDirichletBoundaryConditions_H
-
-#include <core/vectors/tnlStaticVector.h>
-#include <config/tnlParameterContainer.h>
-#include <functors/tnlConstantFunction.h>
-#include <core/vectors/tnlSharedVector.h>
-
-template< typename Mesh,
-          typename Function = tnlConstantFunction< Mesh::Dimensions,
-                                                   typename Mesh::RealType >,
-          typename Real = typename Mesh::RealType,
-          typename Index = typename Mesh::IndexType >
-class tnlAnalyticDirichletBoundaryConditions
-{
-   
-};
-
-template< int Dimensions,
-          typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-class tnlAnalyticDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Function, Real, Index >
-{
-   public:
-   
-      typedef tnlGrid< Dimensions, MeshReal, Device, MeshIndex > MeshType;
-      typedef Real RealType;
-      typedef Device DeviceType;
-      typedef Index IndexType;
-      typedef tnlAnalyticDirichletBoundaryConditions< MeshType, Function, Real, Index > ThisType;
-
-      typedef tnlSharedVector< RealType, DeviceType, IndexType > SharedVector;
-      typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
-      typedef tnlStaticVector< Dimensions, RealType > VertexType;
-      typedef typename MeshType::CoordinatesType CoordinatesType;
-
-      static void configSetup( tnlConfigDescription& config,
-                               const tnlString& prefix = "" );
-
-      bool setup( const tnlParameterContainer& parameters,
-                  const tnlString& prefix = "" );
-
-      void setFunction( const Function& function );
-
-      Function& getFunction();
-
-      const Function& getFunction() const;
-
-      template< typename EntityType >
-      __cuda_callable__
-      void setBoundaryConditions( const RealType& time,
-                                  const MeshType& mesh,
-                                  const IndexType index,
-                                  const EntityType& entity,
-                                  DofVectorType& u,
-                                  DofVectorType& fu ) const;
-
-      template< typename EntityType >
-      __cuda_callable__
-      Index getLinearSystemRowLength( const MeshType& mesh,
-                                      const IndexType& index,
-                                      const EntityType& entity ) const;
-
-      template< typename MatrixRow,
-                typename EntityType >
-      __cuda_callable__
-      void updateLinearSystem( const RealType& time,
-                               const MeshType& mesh,
-                               const IndexType& index,
-                               const EntityType& entity,
-                               DofVectorType& u,
-                               DofVectorType& b,
-                               MatrixRow& matrixRow ) const;
-
-   protected:
-
-      Function function;
-};
-
-template< typename Mesh,
-          typename Function, 
-          typename Real,
-          typename Index >
-ostream& operator << ( ostream& str, const tnlAnalyticDirichletBoundaryConditions< Mesh, Function, Real, Index >& bc )
-{
-   str << "Analytic Dirichlet boundary conditions: function = " << bc.getFunction();
-   return str;
-}
-
-#include <operators/tnlAnalyticDirichletBoundaryConditions_impl.h>
-
-#endif	/* tnlAnalyticDirichletBoundaryConditions_H */
diff --git a/src/operators/tnlAnalyticDirichletBoundaryConditions_impl.h b/src/operators/tnlAnalyticDirichletBoundaryConditions_impl.h
deleted file mode 100644
index 4f203d543f..0000000000
--- a/src/operators/tnlAnalyticDirichletBoundaryConditions_impl.h
+++ /dev/null
@@ -1,159 +0,0 @@
-/***************************************************************************
-           tnlAnalyticDirichletBoundaryConditions_impl.h  -  description
-                             -------------------
-    begin                : Nov 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 tnlAnalyticDirichletBoundaryConditions_IMPL_H
-#define	tnlAnalyticDirichletBoundaryConditions_IMPL_H
-
-template< int Dimensions,
-          typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-void
-tnlAnalyticDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Function, Real, Index >::
-configSetup( tnlConfigDescription& config,
-             const tnlString& prefix )
-{
-   Function::configSetup( config, prefix );
-}
-
-template< int Dimensions,
-          typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-bool
-tnlAnalyticDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Function, Real, Index >::
-setup( const tnlParameterContainer& parameters,
-       const tnlString& prefix )
-{
-   return function.setup( parameters, prefix );
-}
-
-template< int Dimensions,
-          typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-void
-tnlAnalyticDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Function, Real, Index >::
-setFunction( const Function& function )
-{
-   this->function = function;
-}
-
-template< int Dimensions,
-          typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-Function&
-tnlAnalyticDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Function, Real, Index >::
-getFunction()
-{
-   return this->function;
-}
-
-template< int Dimensions,
-          typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-const Function&
-tnlAnalyticDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Function, Real, Index >::
-getFunction() const
-{
-   return this->function;
-}
-
-template< int Dimensions,
-          typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-   template< typename EntityType >
-__cuda_callable__
-void
-tnlAnalyticDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Function, Real, Index >::
-setBoundaryConditions( const RealType& time,
-                       const MeshType& mesh,
-                       const IndexType index,
-                       const EntityType& entity,
-                       DofVectorType& u,
-                       DofVectorType& fu ) const
-{
-   fu[ index ] = 0;
-   u[ index ] = function.getValue( entity.getCenter(), time );
-}
-
-template< int Dimensions,
-          typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-   template< typename EntityType >
-__cuda_callable__
-Index
-tnlAnalyticDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Function, Real, Index >::
-getLinearSystemRowLength( const MeshType& mesh,
-                          const IndexType& index,
-                          const EntityType& entity ) const
-{
-   return 1;
-}
-
-template< int Dimensions,
-          typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-   template< typename Matrix,
-             typename EntityType >          
-__cuda_callable__
-void
-tnlAnalyticDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Function, Real, Index >::
-updateLinearSystem( const RealType& time,
-                    const MeshType& mesh,
-                    const IndexType& index,
-                    const EntityType& entity,
-                    DofVectorType& u,
-                    DofVectorType& b,
-                    Matrix& matrix ) const
-{
-   typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
-   matrixRow.setElement( 0, index, 1.0 );
-   b[ index ] = function.getValue( entity.getCenter(), time );
-}
-
-#endif	/* tnlAnalyticDirichletBoundaryConditions_IMPL_H */
-
diff --git a/src/operators/tnlAnalyticNeumannBoundaryConditions.h b/src/operators/tnlAnalyticNeumannBoundaryConditions.h
deleted file mode 100644
index 30e75f60cd..0000000000
--- a/src/operators/tnlAnalyticNeumannBoundaryConditions.h
+++ /dev/null
@@ -1,227 +0,0 @@
-/***************************************************************************
-                          tnlAnalyticNeumannBoundaryConditions.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 TNLANALYTICNEUMANNBOUNDARYCONDITIONS_H_
-#define TNLANALYTICNEUMANNBOUNDARYCONDITIONS_H_
-
-
-template< typename Mesh,
-          typename Function = tnlConstantFunction< Mesh::Dimensions,
-                                                   typename Mesh::RealType >,
-          typename Real = typename Mesh::RealType,
-          typename Index = typename Mesh::IndexType >
-class tnlAnalyticNeumannBoundaryConditions
-{
-
-};
-
-template< typename Function >
-class tnlAnalyticNeumannBoundaryConditionsBase
-{
-   public:
-
-      static void configSetup( tnlConfigDescription& config,
-                               const tnlString& prefix = "" );
-
-      bool setup( const tnlParameterContainer& parameters,
-                  const tnlString& prefix = "" );
-
-      void setFunction( const Function& function );
-
-      Function& getFunction();
-
-      const Function& getFunction() const;
-
-   protected:
-
-      Function function;
-};
-
-/****
- * 1D grid
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-class tnlAnalyticNeumannBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >
-   : public tnlAnalyticNeumannBoundaryConditionsBase< Function >
-{
-   public:
-
-
-   typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;
-   typedef Real RealType;
-   typedef Device DeviceType;
-   typedef Index IndexType;
-   typedef tnlAnalyticNeumannBoundaryConditions< MeshType, Function, Real, Index > ThisType;
-
-   typedef tnlSharedVector< RealType, DeviceType, IndexType > SharedVector;
-   typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
-   typedef tnlStaticVector< 1, RealType > VertexType;
-   typedef typename MeshType::CoordinatesType CoordinatesType;
-
-
-   template< typename EntityType >
-   __cuda_callable__
-   void setBoundaryConditions( const RealType& time,
-                               const MeshType& mesh,
-                               const IndexType index,
-                               const EntityType& entity,
-                               DofVectorType& u,
-                               DofVectorType& fu ) const;
-
-   template< typename EntityType >
-   __cuda_callable__
-   Index getLinearSystemRowLength( const MeshType& mesh,
-                                   const IndexType& index,
-                                   const EntityType& entity ) const;
-
-   template< typename MatrixRow,
-             typename EntityType >
-   __cuda_callable__
-      void updateLinearSystem( const RealType& time,
-                               const MeshType& mesh,
-                               const IndexType& index,
-                               const EntityType& entity,
-                               DofVectorType& u,
-                               DofVectorType& b,
-                               MatrixRow& matrixRow ) const;
-
-};
-
-/****
- * 2D grid
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-class tnlAnalyticNeumannBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >
-   : public tnlAnalyticNeumannBoundaryConditionsBase< Function >
-{
-   public:
-
-   typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;
-   typedef Real RealType;
-   typedef Device DeviceType;
-   typedef Index IndexType;
-
-   typedef tnlSharedVector< RealType, DeviceType, IndexType > SharedVector;
-   typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
-   typedef tnlStaticVector< 2, RealType > VertexType;
-   typedef typename MeshType::CoordinatesType CoordinatesType;
-
-   template< typename EntityType >
-   __cuda_callable__  
-   void setBoundaryConditions( const RealType& time,
-                               const MeshType& mesh,
-                               const IndexType index,
-                               const EntityType& entity,
-                               DofVectorType& u,
-                               DofVectorType& fu ) const;
-
-   template< typename EntityType >
-   __cuda_callable__
-   Index getLinearSystemRowLength( const MeshType& mesh,
-                                   const IndexType& index,
-                                   const EntityType& entity ) const;
-
-   template< typename MatrixRow,
-             typename EntityType >
-   __cuda_callable__
-      void updateLinearSystem( const RealType& time,
-                               const MeshType& mesh,
-                               const IndexType& index,
-                               const EntityType& entity,
-                               DofVectorType& u,
-                               DofVectorType& b,
-                               MatrixRow& matrixRow ) const;
-};
-
-/****
- * 3D grid
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-class tnlAnalyticNeumannBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >
-   : public tnlAnalyticNeumannBoundaryConditionsBase< Function >
-{
-   public:
-
-   typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;
-   typedef Real RealType;
-   typedef Device DeviceType;
-   typedef Index IndexType;
-
-   typedef tnlSharedVector< RealType, DeviceType, IndexType > SharedVector;
-   typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
-   typedef tnlStaticVector< 3, RealType > VertexType;
-   typedef typename MeshType::CoordinatesType CoordinatesType;
-
-
-   template< typename EntityType >
-   __cuda_callable__
-   void setBoundaryConditions( const RealType& time,
-                               const MeshType& mesh,
-                               const IndexType index,
-                               const EntityType& entity,
-                               DofVectorType& u,
-                               DofVectorType& fu ) const;
-
-   template< typename EntityType >
-   __cuda_callable__
-   Index getLinearSystemRowLength( const MeshType& mesh,
-                                   const IndexType& index,
-                                   const EntityType& entity ) const;
-
-   template< typename MatrixRow,
-             typename EntityType >
-   __cuda_callable__
-      void updateLinearSystem( const RealType& time,
-                               const MeshType& mesh,
-                               const IndexType& index,
-                               const EntityType& enity,
-                               DofVectorType& u,
-                               DofVectorType& b,
-                               MatrixRow& matrixRow ) const;
-};
-
-template< typename Mesh,
-          typename Function,
-          typename Real,
-          typename Index >
-ostream& operator << ( ostream& str, const tnlAnalyticNeumannBoundaryConditions< Mesh, Function, Real, Index >& bc )
-{
-   str << "Analytic Dirichlet boundary conditions: function = " << bc.getFunction();
-   return str;
-}
-
-#include <operators/tnlAnalyticNeumannBoundaryConditions_impl.h>
-
-#endif /* TNLANALYTICNEUMANNBOUNDARYCONDITIONS_H_ */
diff --git a/src/operators/tnlAnalyticNeumannBoundaryConditions_impl.h b/src/operators/tnlAnalyticNeumannBoundaryConditions_impl.h
deleted file mode 100644
index 298f5c1cda..0000000000
--- a/src/operators/tnlAnalyticNeumannBoundaryConditions_impl.h
+++ /dev/null
@@ -1,393 +0,0 @@
-/***************************************************************************
-                          tnlAnalyticNeumannBoundaryConditions_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 TNLANALYTICNEUMANNBOUNDARYCONDITIONS_IMPL_H_
-#define TNLANALYTICNEUMANNBOUNDARYCONDITIONS_IMPL_H_
-
-/****
- * Base
- */
-template< typename Function >
-void
-tnlAnalyticNeumannBoundaryConditionsBase< Function >::
-configSetup( tnlConfigDescription& config,
-             const tnlString& prefix )
-{
-   Function::configSetup( config, prefix );
-}
-
-template< typename Function >
-bool
-tnlAnalyticNeumannBoundaryConditionsBase< Function >::
-setup( const tnlParameterContainer& parameters,
-       const tnlString& prefix )
-{
-   return function.setup( parameters, prefix );
-}
-
-template< typename Function >
-void
-tnlAnalyticNeumannBoundaryConditionsBase< Function >::
-setFunction( const Function& function )
-{
-   this->function = function;
-}
-
-template< typename Function >
-Function&
-tnlAnalyticNeumannBoundaryConditionsBase< Function >::
-getFunction()
-{
-   return this->function;
-}
-
-template< typename Function >
-const Function&
-tnlAnalyticNeumannBoundaryConditionsBase< Function >::
-getFunction() const
-{
-   return this->function;
-}
-
-/****
- * 1D grid
- */
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-   template< typename EntityType >
-__cuda_callable__
-void
-tnlAnalyticNeumannBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >::
-setBoundaryConditions( const RealType& time,
-                       const MeshType& mesh,
-                       const IndexType index,
-                       const EntityType& entity,
-                       DofVectorType& u,
-                       DofVectorType& fu ) const
-{
-   auto neighbourEntities = entity.getNeighbourEntities();
-   fu[ index ] = 0;
-   const Real functionValue = this->function.getValue( entity.getCenter(), time );
-   if( entity.getCoordinates().x() == 0 )
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< 1 >() ] - mesh.getSpaceSteps().x() * functionValue;
-   else
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< -1 >() ] + mesh.getSpaceSteps().x() * functionValue;
-}
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-   template< typename EntityType >
-__cuda_callable__
-Index
-tnlAnalyticNeumannBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >::
-getLinearSystemRowLength( const MeshType& mesh,
-                          const IndexType& index,
-                          const EntityType& entity ) const
-{
-   return 2;
-}
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-   template< typename Matrix,
-             typename EntityType >
-__cuda_callable__
-void
-tnlAnalyticNeumannBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >::
-updateLinearSystem( const RealType& time,
-                    const MeshType& mesh,
-                    const IndexType& index,
-                    const EntityType& entity,
-                    DofVectorType& u,
-                    DofVectorType& b,
-                    Matrix& matrix ) const
-{
-   auto neighbourEntities = entity.getNeighbourEntities();
-   typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
-   const Real functionValue = this->function.getValue( entity.getCenter(), time );
-   if( entity.getCoordinates().x() == 0 )
-   {
-      matrixRow.setElement( 0, index,                            1.0 );
-      matrixRow.setElement( 1, neighbourEntities.template getEntityIndex< 1 >(), -1.0 );
-      b[ index ] = - mesh.getSpaceSteps().x() * functionValue;
-   }
-   else
-   {
-      matrixRow.setElement( 0, neighbourEntities.template getEntityIndex< -1 >(), -1.0 );
-      matrixRow.setElement( 1, index,                              1.0 );
-      b[ index ] = mesh.getSpaceSteps().x() * functionValue;
-   }
-}
-
-/****
- * 2D grid
- */
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-   template< typename EntityType >          
-__cuda_callable__
-void
-tnlAnalyticNeumannBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >::
-setBoundaryConditions( const RealType& time,
-                       const MeshType& mesh,
-                       const IndexType index,
-                       const EntityType& entity,
-                       DofVectorType& u,
-                       DofVectorType& fu ) const
-{
-   auto neighbourEntities = entity.getNeighbourEntities();
-   fu[ index ] = 0;
-   const Real functionValue = this->function.getValue( entity.getCenter(), time );
-   if( entity.getCoordinates().x() == 0 )
-   {
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< 1, 0 >() ] - mesh.getSpaceSteps().x() * functionValue;
-      return;
-   }
-   if( entity.getCoordinates().x() == mesh.getDimensions().x() - 1 )
-   {
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< -1, 0 >() ] + mesh.getSpaceSteps().x() * functionValue;
-      return;
-   }
-   if( entity.getCoordinates().y() == 0 )
-   {
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< 0, 1 >() ] - mesh.getSpaceSteps().y() * functionValue;
-      return;
-   }
-   if( entity.getCoordinates().y() == mesh.getDimensions().y() - 1 )
-   {
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< 0, -1 >() ] + mesh.getSpaceSteps().y() * functionValue;
-      return;
-   }
-}
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-   template< typename EntityType >          
-__cuda_callable__
-Index
-tnlAnalyticNeumannBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >::
-getLinearSystemRowLength( const MeshType& mesh,
-                          const IndexType& index,
-                          const EntityType& entity ) const
-{
-   return 2;
-}
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-   template< typename Matrix,
-             typename EntityType >
-__cuda_callable__
-void
-tnlAnalyticNeumannBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >::
-updateLinearSystem( const RealType& time,
-                    const MeshType& mesh,
-                    const IndexType& index,
-                    const EntityType& entity,
-                    DofVectorType& u,
-                    DofVectorType& b,
-                    Matrix& matrix ) const
-{
-   typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
-   const Real functionValue = this->function.getValue( entity.getCenter(), time );
-   auto neighbourEntities = entity.getNeighbourEntities();
-   if( entity.getCoordinates().x() == 0 )
-   {
-      matrixRow.setElement( 0, index,                                                1.0 );
-      matrixRow.setElement( 1, neighbourEntities.template getEntityIndex< 1, 0 >(), -1.0 );
-      b[ index ] = - mesh.getSpaceSteps().x() * functionValue;
-   }
-   if( entity.getCoordinates().x() == mesh.getDimensions().x() - 1 )
-   {
-      matrixRow.setElement( 0, neighbourEntities.template getEntityIndex< -1, 0 >(), -1.0 );
-      matrixRow.setElement( 1, index,                                                 1.0 );
-      b[ index ] = mesh.getSpaceSteps().x() * functionValue;
-   }
-   if( entity.getCoordinates().y() == 0 )
-   {
-      matrixRow.setElement( 0, index,                                                1.0 );
-      matrixRow.setElement( 1, neighbourEntities.template getEntityIndex< 0, 1 >(), -1.0 );
-      b[ index ] = - mesh.getSpaceSteps().y() * functionValue;
-   }
-   if( entity.getCoordinates().y() == mesh.getDimensions().y() - 1 )
-   {
-      matrixRow.setElement( 0, neighbourEntities.template getEntityIndex< 0, -1 >(), -1.0 );
-      matrixRow.setElement( 1, index,                                                 1.0 );
-      b[ index ] = mesh.getSpaceSteps().y() * functionValue;
-   }
-}
-
-/****
- * 3D grid
- */
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-   template< typename EntityType >          
-__cuda_callable__
-void
-tnlAnalyticNeumannBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >::
-setBoundaryConditions( const RealType& time,
-                       const MeshType& mesh,
-                       const IndexType index,
-                       const EntityType& entity,
-                       DofVectorType& u,
-                       DofVectorType& fu ) const
-{
-   auto neighbourEntities = entity.getNeighbourEntities();
-   fu[ index ] = 0;
-   fu[ index ] = 0;
-   const Real functionValue = this->function.getValue( entity.getCenter(), time );
-   if( entity.getCoordinates().x() == 0 )
-   {
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< 1, 0, 0 >() ] - mesh.getSpaceSteps().x() * functionValue;
-      return;
-   }
-   if( entity.getCoordinates().x() == mesh.getDimensions().x() - 1 )
-   {
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< -1, 0, 0 >() ] + mesh.getSpaceSteps().x() * functionValue;
-      return;
-   }
-   if( entity.getCoordinates().y() == 0 )
-   {
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< 0, 1, 0 >() ] - mesh.getSpaceSteps().y() * functionValue;
-      return;
-   }
-   if( entity.getCoordinates().y() == mesh.getDimensions().y() - 1 )
-   {
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< 0, -1, 0 >() ] + mesh.getSpaceSteps().y() * functionValue;
-      return;
-   }
-   if( entity.getCoordinates().z() == 0 )
-   {
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< 0, 0, 1 >() ] - mesh.getSpaceSteps().z() * functionValue;
-      return;
-   }
-   if( entity.getCoordinates().z() == mesh.getDimensions().z() - 1 )
-   {
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< 0, 0, - 1 >() ] + mesh.getSpaceSteps().z() * functionValue;
-      return;
-   }
-}
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-   template< typename EntityType >          
-__cuda_callable__
-Index
-tnlAnalyticNeumannBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >::
-getLinearSystemRowLength( const MeshType& mesh,
-                          const IndexType& index,
-                          const EntityType& entity ) const
-{
-   return 2;
-}
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Function,
-          typename Real,
-          typename Index >
-   template< typename Matrix,
-             typename EntityType >
-__cuda_callable__
-void
-tnlAnalyticNeumannBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >::
-updateLinearSystem( const RealType& time,
-                    const MeshType& mesh,
-                    const IndexType& index,
-                    const EntityType& entity,
-                    DofVectorType& u,
-                    DofVectorType& b,
-                    Matrix& matrix ) const
-{
-   auto neighbourEntities = entity.getNeighbourEntities();
-   typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
-   const Real functionValue = this->function.getValue( entity.getCenter(), time );
-   if( entity.getCoordinates().x() == 0 )
-   {
-      matrixRow.setElement( 0, index,                                                   1.0 );
-      matrixRow.setElement( 1, neighbourEntities.template getEntityIndex< 1, 0, 0 >(), -1.0 );
-      b[ index ] = - mesh.getSpaceSteps().x() * functionValue;
-   }
-   if( entity.getCoordinates().x() == mesh.getDimensions().x() - 1 )
-   {
-      matrixRow.setElement( 0, neighbourEntities.template getEntityIndex< -1, 0, 0 >(), -1.0 );
-      matrixRow.setElement( 1, index,                                                    1.0 );
-      b[ index ] = mesh.getSpaceSteps().x() * functionValue;
-   }
-   if( entity.getCoordinates().y() == 0 )
-   {
-      matrixRow.setElement( 0, index,                                                   1.0 );
-      matrixRow.setElement( 1, neighbourEntities.template getEntityIndex< 0, 1, 0 >(), -1.0 );
-      b[ index ] = - mesh.getSpaceSteps().y() * functionValue;
-   }
-   if( entity.getCoordinates().y() == mesh.getDimensions().y() - 1 )
-   {
-      matrixRow.setElement( 0, neighbourEntities.template getEntityIndex< 0, -1, 0 >(), -1.0 );
-      matrixRow.setElement( 1, index,                                                    1.0 );
-      b[ index ] = mesh.getSpaceSteps().y() * functionValue;
-   }
-   if( entity.getCoordinates().z() == 0 )
-   {
-      matrixRow.setElement( 0, index,                                                   1.0 );
-      matrixRow.setElement( 1, neighbourEntities.template getEntityIndex< 0, 0, 1 >(), -1.0 );
-      b[ index ] = - mesh.getSpaceSteps().z() * functionValue;
-   }
-   if( entity.getCoordinates().z() == mesh.getDimensions().z() - 1 )
-   {
-      matrixRow.setElement( 0, neighbourEntities.template getEntityIndex< 0, 0, -1 >(), -1.0 );
-      matrixRow.setElement( 1, index,                                                    1.0 );
-      b[ index ] = mesh.getSpaceSteps().z() * functionValue;
-   }
-}
-
-
-#endif /* TNLANALYTICNEUMANNBOUNDARYCONDITIONS_IMPL_H_ */
diff --git a/src/operators/tnlDirichletBoundaryConditions.h b/src/operators/tnlDirichletBoundaryConditions.h
index 4dc8387afe..d718eaf6ce 100644
--- a/src/operators/tnlDirichletBoundaryConditions.h
+++ b/src/operators/tnlDirichletBoundaryConditions.h
@@ -19,7 +19,7 @@
 #define TNLDIRICHLETBOUNDARYCONDITIONS_H_
 
 template< typename Mesh,
-          typename Vector,
+          typename Function,
           typename Real = typename Mesh::RealType,
           typename Index = typename Mesh::IndexType >
 class tnlDirichletBoundaryConditions
@@ -31,48 +31,50 @@ template< int Dimensions,
           typename MeshReal,
           typename Device,
           typename MeshIndex,
-          typename Vector,
+          typename Function,
           typename Real,
           typename Index >
-class tnlDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Vector, Real, Index >
+class tnlDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Function, Real, Index >
 {
    public:
 
    typedef tnlGrid< Dimensions, MeshReal, Device, MeshIndex > MeshType;
+   typedef Function FunctionType;
    typedef Real RealType;
    typedef Device DeviceType;
    typedef Index IndexType;
 
-   typedef Vector VectorType;
-   typedef tnlSharedVector< RealType, DeviceType, IndexType > SharedVector;
+   
+   //typedef tnlSharedVector< RealType, DeviceType, IndexType > SharedVector;
    typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
    typedef tnlStaticVector< Dimensions, RealType > VertexType;
-   typedef typename MeshType::CoordinatesType CoordinatesType;
+   //typedef typename MeshType::CoordinatesType CoordinatesType;
 
-   void configSetup( tnlConfigDescription& config,
-                     const tnlString& prefix );
+   static void configSetup( tnlConfigDescription& config,
+                            const tnlString& prefix = "" );
 
    bool setup( const tnlParameterContainer& parameters,
                const tnlString& prefix = "" );
 
-   Vector& getVector();
+   void setFunction( const Function& function );
+   
+   Function& getFunction();
 
-   const Vector& getVector() const;
+   const Function& getFunction() const;
 
    template< typename EntityType >
    __cuda_callable__
    void setBoundaryConditions( const RealType& time,
                                const MeshType& mesh,
-                               const IndexType index,
                                const EntityType& entity,
                                DofVectorType& u,
                                DofVectorType& fu ) const;
 
    template< typename EntityType >
    __cuda_callable__
-   Index getLinearSystemRowLength( const MeshType& mesh,
-                                   const IndexType& index,
-                                   const EntityType& entity ) const;
+   IndexType getLinearSystemRowLength( const MeshType& mesh,
+                                       const IndexType& index,
+                                       const EntityType& entity ) const;
 
    template< typename MatrixRow,
              typename EntityType >
@@ -87,14 +89,14 @@ class tnlDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, Mes
 
    protected:
 
-   Vector vector;
+   Function function;
+   
+   //static_assert( Device::DeviceType == Function::Device::DeviceType );
 };
 
 template< typename Mesh,
-          typename Function,
-          typename Real,
-          typename Index >
-ostream& operator << ( ostream& str, const tnlDirichletBoundaryConditions< Mesh, Function, Real, Index >& bc )
+          typename Function >
+ostream& operator << ( ostream& str, const tnlDirichletBoundaryConditions< Mesh, Function >& bc )
 {
    str << "Dirichlet boundary conditions: vector = " << bc.getVector();
    return str;
diff --git a/src/operators/tnlDirichletBoundaryConditions_impl.h b/src/operators/tnlDirichletBoundaryConditions_impl.h
index 475c4021d9..c91c15c775 100644
--- a/src/operators/tnlDirichletBoundaryConditions_impl.h
+++ b/src/operators/tnlDirichletBoundaryConditions_impl.h
@@ -18,68 +18,79 @@
 #ifndef TNLDIRICHLETBOUNDARYCONDITIONS_IMPL_H_
 #define TNLDIRICHLETBOUNDARYCONDITIONS_IMPL_H_
 
+#include <functions/tnlFunctionAdapter.h>
+
 template< int Dimensions,
           typename MeshReal,
           typename Device,
           typename MeshIndex,
-          typename Vector,
+          typename Function,
           typename Real,
           typename Index >
 void
-tnlDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Vector, Real, Index >::
+tnlDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Function, Real, Index >::
 configSetup( tnlConfigDescription& config,
              const tnlString& prefix )
 {
-   config.addEntry     < tnlString >( prefix + "file", "Data for the boundary conditions." );
+   Function::configSetup( config );
 }
 
 template< int Dimensions,
           typename MeshReal,
           typename Device,
           typename MeshIndex,
-          typename Vector,
+          typename Function,
           typename Real,
           typename Index >
 bool
-tnlDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Vector, Real, Index >::
+tnlDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Function, Real, Index >::
 setup( const tnlParameterContainer& parameters,
        const tnlString& prefix )
 {
-   if( parameters.checkParameter( prefix + "file" ) )
-   {
-      tnlString fileName = parameters.getParameter< tnlString >( prefix + "file" );
-      if( ! this->vector.load( fileName ) )
-         return false;
-   }
-   return true;
+   return this->function.setup( parameters );
+}
+
+template< int Dimensions,
+          typename MeshReal,
+          typename Device,
+          typename MeshIndex,
+          typename Function,
+          typename Real,
+          typename Index >
+void
+tnlDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Function, Real, Index >::
+setFunction( const Function& function )
+{
+   this->function = function;
 }
 
+
 template< int Dimensions,
           typename MeshReal,
           typename Device,
           typename MeshIndex,
-          typename Vector,
+          typename Function,
           typename Real,
           typename Index >
-Vector&
-tnlDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Vector, Real, Index >::
-getVector()
+Function&
+tnlDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Function, Real, Index >::
+getFunction()
 {
-   return this->vector;
+   return this->function;
 }
 
 template< int Dimensions,
           typename MeshReal,
           typename Device,
           typename MeshIndex,
-          typename Vector,
+          typename Function,
           typename Real,
           typename Index >
-const Vector&
-tnlDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Vector, Real, Index >::
-getVector() const
+const Function&
+tnlDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Function, Real, Index >::
+getFunction() const
 {
-   return this->vector;
+   return *this->function;
 }
 
 
@@ -87,35 +98,35 @@ template< int Dimensions,
           typename MeshReal,
           typename Device,
           typename MeshIndex,
-          typename Vector,
+          typename Function,
           typename Real,
           typename Index >
    template< typename EntityType >
 __cuda_callable__
 void
-tnlDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Vector, Real, Index >::
+tnlDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Function, Real, Index >::
 setBoundaryConditions( const RealType& time,
                        const MeshType& mesh,
-                       const IndexType index,
                        const EntityType& entity,
                        DofVectorType& u,
                        DofVectorType& fu ) const
 {
+   const IndexType& index = entity.getIndex();
    fu[ index ] = 0;
-   u[ index ] = this->vector[ index ];
+   u[ index ] = tnlFunctionAdapter< MeshType, Function >::template getValue( this->function, entity, time );
 }
 
 template< int Dimensions,
           typename MeshReal,
           typename Device,
           typename MeshIndex,
-          typename Vector,
+          typename Function,
           typename Real,
           typename Index >
    template< typename EntityType >
 __cuda_callable__
 Index
-tnlDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Vector, Real, Index >::
+tnlDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Function, Real, Index >::
 getLinearSystemRowLength( const MeshType& mesh,
                           const IndexType& index,
                           const EntityType& entity ) const
@@ -127,14 +138,14 @@ template< int Dimensions,
           typename MeshReal,
           typename Device,
           typename MeshIndex,
-          typename Vector,
+          typename Function,
           typename Real,
           typename Index >
    template< typename Matrix,
              typename EntityType >
 __cuda_callable__
 void
-tnlDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Vector, Real, Index >::
+tnlDirichletBoundaryConditions< tnlGrid< Dimensions, MeshReal, Device, MeshIndex >, Function, Real, Index >::
 updateLinearSystem( const RealType& time,
                     const MeshType& mesh,
                     const IndexType& index,
@@ -145,7 +156,7 @@ updateLinearSystem( const RealType& time,
 {
    typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
    matrixRow.setElement( 0, index, 1.0 );
-   b[ index ] = this->vector[ index ];
+   b[ index ] = tnlFunctionAdapter< MeshType, Function >::getValue( this->function, entity, time );
 }
 
 
diff --git a/src/operators/tnlExactOperatorEvaluator.h b/src/operators/tnlExactOperatorEvaluator.h
index f01376a626..c9c5d2eac9 100644
--- a/src/operators/tnlExactOperatorEvaluator.h
+++ b/src/operators/tnlExactOperatorEvaluator.h
@@ -156,33 +156,15 @@ class tnlExactOperatorEvaluator< tnlGrid< Dimensions, Real, Device, Index >, Dof
             __cuda_callable__
             static void processEntity( const MeshType& mesh,
                                        TraversalUserData& userData,
-                                       const IndexType index,
                                        const EntityType& entity )
             {
                userData.boundaryConditions.setBoundaryConditions
                   ( userData.time,
                     mesh,
-                    index,
                     entity,
                     userData.fu,
                     userData.fu );
             }
-
-            
-            __cuda_callable__
-            static void processCell( const MeshType& mesh,
-                                     TraversalUserData& userData,
-                                     const IndexType index,
-                                     const CoordinatesType& c )
-            {
-               userData.boundaryConditions.setBoundaryConditions( userData.time,
-                                                                  mesh,
-                                                                  index,
-                                                                  c,
-                                                                  userData.fu,
-                                                                  userData.fu );
-            }
-
       };
 
       class TraversalInteriorEntitiesProcessor
@@ -193,10 +175,9 @@ class tnlExactOperatorEvaluator< tnlGrid< Dimensions, Real, Device, Index >, Dof
             __cuda_callable__
             static void processEntity( const MeshType& mesh,
                                        TraversalUserData& userData,
-                                       const IndexType index,
                                        const EntityType& entity )
             {
-               userData.fu[ index ] = 
+               userData.fu[ entity.getIndex() ] = 
                   userData.differentialOperator.getValue
                      ( userData.function,
                        entity.getCenter(),                                                                           
diff --git a/src/operators/tnlNeumannBoundaryConditions.h b/src/operators/tnlNeumannBoundaryConditions.h
index 2df09c3e30..dc87909308 100644
--- a/src/operators/tnlNeumannBoundaryConditions.h
+++ b/src/operators/tnlNeumannBoundaryConditions.h
@@ -2,7 +2,7 @@
 #define	TNLNEUMANNBOUNDARYCONDITIONS_H
 
 template< typename Mesh,
-          typename Vector,
+          typename Function,
           typename Real = typename Mesh::RealType,
           typename Index = typename Mesh::IndexType >
 class tnlNeumannBoundaryConditions
@@ -13,10 +13,12 @@ class tnlNeumannBoundaryConditions
 /****
  * Base
  */
-template< typename Vector >
+template< typename Function >
 class tnlNeumannBoundaryConditionsBase
 {
    public:
+      
+      typedef Function FunctionType;
 
       static void configSetup( tnlConfigDescription& config,
                                const tnlString& prefix = "" );
@@ -24,13 +26,15 @@ class tnlNeumannBoundaryConditionsBase
       bool setup( const tnlParameterContainer& parameters,
                   const tnlString& prefix = "" );
 
-      Vector& getVector();
+      void setFunction( const FunctionType& function );
+      
+      FunctionType& getFunction();
 
-      const Vector& getVector() const;
+      const FunctionType& getFunction() const;
 
    protected:
 
-      Vector vector;
+      FunctionType function;
 
 };
 
@@ -40,11 +44,11 @@ class tnlNeumannBoundaryConditionsBase
 template< typename MeshReal,
           typename Device,
           typename MeshIndex,
-          typename Vector,
+          typename Function,
           typename Real,
           typename Index >
-class tnlNeumannBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >, Vector, Real, Index >
-   : public tnlNeumannBoundaryConditionsBase< Vector >
+class tnlNeumannBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >
+   : public tnlNeumannBoundaryConditionsBase< Function >
 {
    public:
 
@@ -53,7 +57,7 @@ class tnlNeumannBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >, V
    typedef Device DeviceType;
    typedef Index IndexType;
 
-   typedef Vector VectorType;
+   typedef Function FunctionType;
    typedef tnlSharedVector< RealType, DeviceType, IndexType > SharedVector;
    typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
    typedef tnlStaticVector< 1, RealType > VertexType;
@@ -63,7 +67,6 @@ class tnlNeumannBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >, V
    __cuda_callable__
    void setBoundaryConditions( const RealType& time,
                                const MeshType& mesh,
-                               const IndexType index,
                                const EntityType& entity,
                                DofVectorType& u,
                                DofVectorType& fu ) const;
@@ -92,11 +95,11 @@ class tnlNeumannBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >, V
 template< typename MeshReal,
           typename Device,
           typename MeshIndex,
-          typename Vector,
+          typename Function,
           typename Real,
           typename Index >
-class tnlNeumannBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, Vector, Real, Index >
-   : public tnlNeumannBoundaryConditionsBase< Vector >
+class tnlNeumannBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >
+   : public tnlNeumannBoundaryConditionsBase< Function >
 {
    public:
 
@@ -105,7 +108,7 @@ class tnlNeumannBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, V
    typedef Device DeviceType;
    typedef Index IndexType;
 
-   typedef Vector VectorType;
+   typedef Function FunctionType;
    typedef tnlSharedVector< RealType, DeviceType, IndexType > SharedVector;
    typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
    typedef tnlStaticVector< 2, RealType > VertexType;
@@ -115,7 +118,6 @@ class tnlNeumannBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, V
    __cuda_callable__
    void setBoundaryConditions( const RealType& time,
                                const MeshType& mesh,
-                               const IndexType index,
                                const EntityType& entity,
                                DofVectorType& u,
                                DofVectorType& fu ) const;
@@ -144,11 +146,11 @@ class tnlNeumannBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, V
 template< typename MeshReal,
           typename Device,
           typename MeshIndex,
-          typename Vector,
+          typename Function,
           typename Real,
           typename Index >
-class tnlNeumannBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >, Vector, Real, Index >
-   : public tnlNeumannBoundaryConditionsBase< Vector >
+class tnlNeumannBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >
+   : public tnlNeumannBoundaryConditionsBase< Function >
 {
    public:
 
@@ -157,7 +159,7 @@ class tnlNeumannBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >, V
    typedef Device DeviceType;
    typedef Index IndexType;
 
-   typedef Vector VectorType;
+   typedef Function FunctionType;
    typedef tnlSharedVector< RealType, DeviceType, IndexType > SharedVector;
    typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
    typedef tnlStaticVector< 3, RealType > VertexType;
@@ -167,7 +169,6 @@ class tnlNeumannBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >, V
    __cuda_callable__
    void setBoundaryConditions( const RealType& time,
                                const MeshType& mesh,
-                               const IndexType index,
                                const EntityType& entity,
                                DofVectorType& u,
                                DofVectorType& fu ) const;
@@ -196,7 +197,7 @@ template< typename Mesh,
           typename Index >
 ostream& operator << ( ostream& str, const tnlNeumannBoundaryConditions< Mesh, Function, Real, Index >& bc )
 {
-   str << "Dirichlet boundary conditions: vector = " << bc.getVector();
+   str << "Dirichlet boundary conditions: vector = " << bc.getFunction();
    return str;
 }
 
diff --git a/src/operators/tnlNeumannBoundaryConditions_impl.h b/src/operators/tnlNeumannBoundaryConditions_impl.h
index a7e5e3df6b..4133edda92 100644
--- a/src/operators/tnlNeumannBoundaryConditions_impl.h
+++ b/src/operators/tnlNeumannBoundaryConditions_impl.h
@@ -1,44 +1,47 @@
 #ifndef TNLNEUMANNBOUNDARYCONDITIONS_IMPL_H
 #define	TNLNEUMANNBOUNDARYCONDITIONS_IMPL_H
 
-template< typename Vector >
+template< typename Function >
 void
-tnlNeumannBoundaryConditionsBase< Vector >::
+tnlNeumannBoundaryConditionsBase< Function >::
 configSetup( tnlConfigDescription& config,
              const tnlString& prefix )
 {
-   config.addEntry     < tnlString >( prefix + "file", "Data for the boundary conditions." );
+   Function::configSetup( config );
 }
 
-template< typename Vector >
+template< typename Function >
 bool
-tnlNeumannBoundaryConditionsBase< Vector >::
+tnlNeumannBoundaryConditionsBase< Function >::
 setup( const tnlParameterContainer& parameters,
        const tnlString& prefix )
 {
-   if( parameters.checkParameter( prefix + "file" ) )
-   {
-      tnlString fileName = parameters.getParameter< tnlString >( prefix + "file" );
-      if( ! this->vector.load( fileName ) )
-         return false;
-   }
-   return true;
+   return this->function.setup( parameters );
 }
 
-template< typename Vector >
-Vector&
-tnlNeumannBoundaryConditionsBase< Vector >::
-getVector()
+template< typename Function >
+void
+tnlNeumannBoundaryConditionsBase< Function >::
+setFunction( const Function& function )
+{
+   this->function = function;
+}
+
+
+template< typename Function >
+Function&
+tnlNeumannBoundaryConditionsBase< Function >::
+getFunction()
 {
-   return this->vector;
+   return this->function;
 }
 
-template< typename Vector >
-const Vector&
-tnlNeumannBoundaryConditionsBase< Vector >::
-getVector() const
+template< typename Function >
+const Function&
+tnlNeumannBoundaryConditionsBase< Function >::
+getFunction() const
 {
-   return this->vector;
+   return this->function;
 }
 
 /****
@@ -47,38 +50,40 @@ getVector() const
 template< typename MeshReal,
           typename Device,
           typename MeshIndex,
-          typename Vector,
+          typename Function,
           typename Real,
           typename Index >
    template< typename EntityType >          
 __cuda_callable__
 void
-tnlNeumannBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >, Vector, Real, Index >::
+tnlNeumannBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >::
 setBoundaryConditions( const RealType& time,
                        const MeshType& mesh,
-                       const IndexType index,
                        const EntityType& entity,
                        DofVectorType& u,
                        DofVectorType& fu ) const
 {
    auto neighbourEntities = entity.getNeighbourEntities();
+   const IndexType& index = entity.getIndex();
    fu[ index ] = 0;
    if( entity.getCoordinates().x() == 0 )
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< 1 >() ] - mesh.getSpaceSteps().x() * this->vector[ index ];
+      u[ index ] = u[ neighbourEntities.template getEntityIndex< 1 >() ] - mesh.getSpaceSteps().x() * 
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
    else
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< -1 >() ] + mesh.getSpaceSteps().x() * this->vector[ index ];
+      u[ index ] = u[ neighbourEntities.template getEntityIndex< -1 >() ] + mesh.getSpaceSteps().x() * 
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
 }
 
 template< typename MeshReal,
           typename Device,
           typename MeshIndex,
-          typename Vector,
+          typename Function,
           typename Real,
           typename Index >
    template< typename EntityType >
 __cuda_callable__
 Index
-tnlNeumannBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >, Vector, Real, Index >::
+tnlNeumannBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >::
 getLinearSystemRowLength( const MeshType& mesh,
                           const IndexType& index,
                           const EntityType& entity ) const
@@ -89,14 +94,14 @@ getLinearSystemRowLength( const MeshType& mesh,
 template< typename MeshReal,
           typename Device,
           typename MeshIndex,
-          typename Vector,
+          typename Function,
           typename Real,
           typename Index >
    template< typename Matrix,
              typename EntityType >
 __cuda_callable__
 void
-tnlNeumannBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >, Vector, Real, Index >::
+tnlNeumannBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >::
 updateLinearSystem( const RealType& time,
                     const MeshType& mesh,
                     const IndexType& index,
@@ -111,13 +116,15 @@ updateLinearSystem( const RealType& time,
    {
       matrixRow.setElement( 0, index, 1.0 );
       matrixRow.setElement( 1, neighbourEntities.template getEntityIndex< 1 >(), -1.0 );
-      b[ index ] = - mesh.getSpaceSteps().x() * this->vector[ index];
+      b[ index ] = - mesh.getSpaceSteps().x() * 
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
    }
    else
    {
       matrixRow.setElement( 0, neighbourEntities.template getEntityIndex< -1 >(), -1.0 );
       matrixRow.setElement( 1, index, 1.0 );
-      b[ index ] = mesh.getSpaceSteps().x() * this->vector[ index ];
+      b[ index ] = mesh.getSpaceSteps().x() *
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
    }
 }
 
@@ -127,40 +134,44 @@ updateLinearSystem( const RealType& time,
 template< typename MeshReal,
           typename Device,
           typename MeshIndex,
-          typename Vector,
+          typename Function,
           typename Real,
           typename Index >
    template< typename EntityType >          
 __cuda_callable__
 void
-tnlNeumannBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, Vector, Real, Index >::
+tnlNeumannBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >::
 setBoundaryConditions( const RealType& time,
                        const MeshType& mesh,
-                       const IndexType index,
                        const EntityType& entity,
                        DofVectorType& u,
                        DofVectorType& fu ) const
 {
    auto neighbourEntities = entity.getNeighbourEntities();
+   const IndexType& index = entity.getIndex();
    fu[ index ] = 0;
    if( entity.getCoordinates().x() == 0 )
    {
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< 1, 0 >() ] - mesh.getSpaceSteps().x() * this->vector[ index ];
+      u[ index ] = u[ neighbourEntities.template getEntityIndex< 1, 0 >() ] - mesh.getSpaceSteps().x() *
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
       return;
    }
    if( entity.getCoordinates().x() == mesh.getDimensions().x() - 1 )
    {
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< -1, 0 >() ] + mesh.getSpaceSteps().x() * this->vector[ index ];
+      u[ index ] = u[ neighbourEntities.template getEntityIndex< -1, 0 >() ] + mesh.getSpaceSteps().x() *
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
       return;
    }
    if( entity.getCoordinates().y() == 0 )
    {
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< 0, 1 >() ] - mesh.getSpaceSteps().y() * this->vector[ index ];
+      u[ index ] = u[ neighbourEntities.template getEntityIndex< 0, 1 >() ] - mesh.getSpaceSteps().y() *
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
       return;
    }
    if( entity.getCoordinates().y() == mesh.getDimensions().y() - 1 )
    {
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< 0, -1 >() ] + mesh.getSpaceSteps().y() * this->vector[ index ];
+      u[ index ] = u[ neighbourEntities.template getEntityIndex< 0, -1 >() ] + mesh.getSpaceSteps().y() *
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
       return;
    }
 }
@@ -168,13 +179,13 @@ setBoundaryConditions( const RealType& time,
 template< typename MeshReal,
           typename Device,
           typename MeshIndex,
-          typename Vector,
+          typename Function,
           typename Real,
           typename Index >
    template< typename EntityType >          
 __cuda_callable__
 Index
-tnlNeumannBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, Vector, Real, Index >::
+tnlNeumannBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >::
 getLinearSystemRowLength( const MeshType& mesh,
                           const IndexType& index,
                           const EntityType& entity ) const
@@ -185,14 +196,14 @@ getLinearSystemRowLength( const MeshType& mesh,
 template< typename MeshReal,
           typename Device,
           typename MeshIndex,
-          typename Vector,
+          typename Function,
           typename Real,
           typename Index >
    template< typename Matrix,
              typename EntityType >
 __cuda_callable__
 void
-tnlNeumannBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, Vector, Real, Index >::
+tnlNeumannBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >::
 updateLinearSystem( const RealType& time,
                     const MeshType& mesh,
                     const IndexType& index,
@@ -207,25 +218,29 @@ updateLinearSystem( const RealType& time,
    {
       matrixRow.setElement( 0, index,                                                1.0 );
       matrixRow.setElement( 1, neighbourEntities.template getEntityIndex< 1, 0 >(), -1.0 );
-      b[ index ] = - mesh.getSpaceSteps().x() * this->vector[ index ];
+      b[ index ] = - mesh.getSpaceSteps().x() *
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
    }
    if( entity.getCoordinates().x() == mesh.getDimensions().x() - 1 )
    {
       matrixRow.setElement( 0, neighbourEntities.template getEntityIndex< -1, 0 >(), -1.0 );
       matrixRow.setElement( 1, index,                                                 1.0 );
-      b[ index ] = mesh.getSpaceSteps().x() * this->vector[ index ];
+      b[ index ] = mesh.getSpaceSteps().x() *
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
    }
    if( entity.getCoordinates().y() == 0 )
    {
       matrixRow.setElement( 0, index,                                                1.0 );
       matrixRow.setElement( 1, neighbourEntities.template getEntityIndex< 0, 1 >(), -1.0 );
-      b[ index ] = - mesh.getSpaceSteps().y() * this->vector[ index ];
+      b[ index ] = - mesh.getSpaceSteps().y() *
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
    }
    if( entity.getCoordinates().y() == mesh.getDimensions().y() - 1 )
    {
       matrixRow.setElement( 0, neighbourEntities.template getEntityIndex< 0, -1 >(), -1.0 );
       matrixRow.setElement( 1, index,                                                 1.0 );
-      b[ index ] = mesh.getSpaceSteps().y() * this->vector[ index ];
+      b[ index ] = mesh.getSpaceSteps().y() *
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
    }
 }
 
@@ -235,50 +250,56 @@ updateLinearSystem( const RealType& time,
 template< typename MeshReal,
           typename Device,
           typename MeshIndex,
-          typename Vector,
+          typename Function,
           typename Real,
           typename Index >
    template< typename EntityType >          
 __cuda_callable__
 void
-tnlNeumannBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >, Vector, Real, Index >::
+tnlNeumannBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >::
 setBoundaryConditions( const RealType& time,
                        const MeshType& mesh,
-                       const IndexType index,
                        const EntityType& entity,
                        DofVectorType& u,
                        DofVectorType& fu ) const
 {
    auto neighbourEntities = entity.getNeighbourEntities();
+   const IndexType& index = entity.getIndex();
    fu[ index ] = 0;
    if( entity.getCoordinates().x() == 0 )
    {
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< 1, 0, 0 >() ] - mesh.getSpaceSteps().x() * this->vector[ index ];
+      u[ index ] = u[ neighbourEntities.template getEntityIndex< 1, 0, 0 >() ] - mesh.getSpaceSteps().x() *
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
       return;
    }
    if( entity.getCoordinates().x() == mesh.getDimensions().x() - 1 )
    {
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< -1, 0, 0 >() ] + mesh.getSpaceSteps().x() * this->vector[ index ];
+      u[ index ] = u[ neighbourEntities.template getEntityIndex< -1, 0, 0 >() ] + mesh.getSpaceSteps().x() *
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
       return;
    }
    if( entity.getCoordinates().y() == 0 )
    {
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< 0, 1, 0 >() ] - mesh.getSpaceSteps().y() * this->vector[ index ];
+      u[ index ] = u[ neighbourEntities.template getEntityIndex< 0, 1, 0 >() ] - mesh.getSpaceSteps().y() *
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
       return;
    }
    if( entity.getCoordinates().y() == mesh.getDimensions().y() - 1 )
    {
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< 0, -1, 0 >() ] + mesh.getSpaceSteps().y() * this->vector[ index ];
+      u[ index ] = u[ neighbourEntities.template getEntityIndex< 0, -1, 0 >() ] + mesh.getSpaceSteps().y() *
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
       return;
    }
    if( entity.getCoordinates().z() == 0 )
    {
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< 0, 0, 1 >() ] - mesh.getSpaceSteps().z() * this->vector[ index ];
+      u[ index ] = u[ neighbourEntities.template getEntityIndex< 0, 0, 1 >() ] - mesh.getSpaceSteps().z() *
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
       return;
    }
    if( entity.getCoordinates().z() == mesh.getDimensions().z() - 1 )
    {
-      u[ index ] = u[ neighbourEntities.template getEntityIndex< 0, 0, -1 >() ] + mesh.getSpaceSteps().z() * this->vector[ index ];
+      u[ index ] = u[ neighbourEntities.template getEntityIndex< 0, 0, -1 >() ] + mesh.getSpaceSteps().z() *
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
       return;
    }
 }
@@ -286,13 +307,13 @@ setBoundaryConditions( const RealType& time,
 template< typename MeshReal,
           typename Device,
           typename MeshIndex,
-          typename Vector,
+          typename Function,
           typename Real,
           typename Index >
    template< typename EntityType >          
 __cuda_callable__
 Index
-tnlNeumannBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >, Vector, Real, Index >::
+tnlNeumannBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >::
 getLinearSystemRowLength( const MeshType& mesh,
                           const IndexType& index,
                           const EntityType& entity ) const
@@ -303,14 +324,14 @@ getLinearSystemRowLength( const MeshType& mesh,
 template< typename MeshReal,
           typename Device,
           typename MeshIndex,
-          typename Vector,
+          typename Function,
           typename Real,
           typename Index >
    template< typename Matrix,
              typename EntityType >
 __cuda_callable__
 void
-tnlNeumannBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >, Vector, Real, Index >::
+tnlNeumannBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >::
 updateLinearSystem( const RealType& time,
                     const MeshType& mesh,
                     const IndexType& index,
@@ -325,37 +346,43 @@ updateLinearSystem( const RealType& time,
    {
       matrixRow.setElement( 0, index,                                                   1.0 );
       matrixRow.setElement( 1, neighbourEntities.template getEntityIndex< 1, 0, 0 >(), -1.0 );
-      b[ index ] = - mesh.getSpaceSteps().x() * this->vector[ index ];
+      b[ index ] = - mesh.getSpaceSteps().x() *
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
    }
    if( entity.getCoordinates().x() == mesh.getDimensions().x() - 1 )
    {
       matrixRow.setElement( 0, neighbourEntities.template getEntityIndex< -1, 0, 0 >(), -1.0 );
       matrixRow.setElement( 1, index,                                                    1.0 );
-      b[ index ] = mesh.getSpaceSteps().x() * this->vector[ index ];
+      b[ index ] = mesh.getSpaceSteps().x() *
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
    }
    if( entity.getCoordinates().y() == 0 )
    {
       matrixRow.setElement( 0, index,                                                   1.0 );
       matrixRow.setElement( 1, neighbourEntities.template getEntityIndex< 0, 1, 0 >(), -1.0 );
-      b[ index ] = - mesh.getSpaceSteps().y() * this->vector[ index ];
+      b[ index ] = - mesh.getSpaceSteps().y() * 
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
    }
    if( entity.getCoordinates().y() == mesh.getDimensions().y() - 1 )
    {
       matrixRow.setElement( 0, neighbourEntities.template getEntityIndex< 0, -1, 0 >(), -1.0 );
       matrixRow.setElement( 1, index,                                                    1.0 );
-      b[ index ] = mesh.getSpaceSteps().y() * this->vector[ index ];
+      b[ index ] = mesh.getSpaceSteps().y() *
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
    }
    if( entity.getCoordinates().z() == 0 )
    {
       matrixRow.setElement( 0, index,                                                   1.0 );
       matrixRow.setElement( 1, neighbourEntities.template getEntityIndex< 0, 0, 1 >(), -1.0 );
-      b[ index ] = - mesh.getSpaceSteps().z() * this->vector[ index ];
+      b[ index ] = - mesh.getSpaceSteps().z() *
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
    }
    if( entity.getCoordinates().z() == mesh.getDimensions().z() - 1 )
    {
       matrixRow.setElement( 0, neighbourEntities.template getEntityIndex< 0, 0, -1 >(), -1.0 );
       matrixRow.setElement( 1, index,                                                    1.0 );
-      b[ index ] = mesh.getSpaceSteps().z() * this->vector[ index ];
+      b[ index ] = mesh.getSpaceSteps().z() *
+         tnlFunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
    }
 }
 
diff --git a/src/problems/tnlHeatEquationEocProblem_impl.h b/src/problems/tnlHeatEquationEocProblem_impl.h
index 0193cfd37d..fe60dbe774 100644
--- a/src/problems/tnlHeatEquationEocProblem_impl.h
+++ b/src/problems/tnlHeatEquationEocProblem_impl.h
@@ -25,6 +25,9 @@
 #ifndef TNLHEATEQUATIONEOCPROBLEM_IMPL_H_
 #define TNLHEATEQUATIONEOCPROBLEM_IMPL_H_
 
+#include "tnlHeatEquationProblem.h"
+
+
 template< typename Mesh,
           typename BoundaryCondition,
           typename RightHandSide,
diff --git a/src/problems/tnlHeatEquationEocRhs.h b/src/problems/tnlHeatEquationEocRhs.h
index d086983aa8..494962f055 100644
--- a/src/problems/tnlHeatEquationEocRhs.h
+++ b/src/problems/tnlHeatEquationEocRhs.h
@@ -24,12 +24,12 @@
 #ifndef TNLHEATEQUATIONEOCRHS_H_
 #define TNLHEATEQUATIONEOCRHS_H_
 
-#include <functors/tnlFunction.h>
+#include <functions/tnlFunction.h>
 
 template< typename ExactOperator,
           typename TestFunction >
 class tnlHeatEquationEocRhs : public tnlFunction< TestFunction::Dimensions,
-                                                  tnlAnalyticFunction >
+                                                  AnalyticFunction >
 {
    public:
 
@@ -38,8 +38,7 @@ class tnlHeatEquationEocRhs : public tnlFunction< TestFunction::Dimensions,
       typedef typename TestFunction::RealType RealType;
       typedef typename TestFunction::VertexType VertexType;
 
-      //static constexpr tnlFunctionType getFunctionType() { return tnlAnalyticFunction; }     
-      enum { functionType = tnlAnalyticFunction };
+      static constexpr tnlFunctionType getFunctionType() { return AnalyticFunction; }     
       
       bool setup( const tnlParameterContainer& parameters,
                   const tnlString& prefix = "" )
diff --git a/src/solvers/pde/tnlExplicitUpdater.h b/src/solvers/pde/tnlExplicitUpdater.h
index 68d88e5573..40606705b7 100644
--- a/src/solvers/pde/tnlExplicitUpdater.h
+++ b/src/solvers/pde/tnlExplicitUpdater.h
@@ -18,7 +18,7 @@
 #ifndef TNLEXPLICITUPDATER_H_
 #define TNLEXPLICITUPDATER_H_
 
-#include <functors/tnlFunctionAdapter.h>
+#include <functions/tnlFunctionAdapter.h>
 
 template< typename Real,
           typename DofVector,
@@ -90,11 +90,11 @@ class tnlExplicitUpdater
             __cuda_callable__
             static inline void processEntity( const MeshType& mesh,
                                               TraverserUserData& userData,
-                                              const IndexType index )
+                                              const EntityType& entity )
             {
                userData.boundaryConditions->setBoundaryConditions( *userData.time,
                                                                    mesh,
-                                                                   index,
+                                                                   entity,
                                                                    *userData.u,
                                                                    *userData.fu );
             }
@@ -109,16 +109,16 @@ class tnlExplicitUpdater
             __cuda_callable__
             static inline void processEntity( const MeshType& mesh,
                                               TraverserUserData& userData,
-                                              const IndexType index )
+                                              const EntityType& entity )
             {
-               (* userData.fu )[ index ] = userData.differentialOperator->getValue( mesh,
-                                                                                    index,
+               (* userData.fu )[ entity.getIndex() ] = userData.differentialOperator->getValue( mesh,
+                                                                                    entity,
                                                                                     *userData.u,
                                                                                     *userData.time );
                typedef tnlFunctionAdapter< MeshType, RightHandSide > FunctionAdapter;
-               ( *userData.fu )[ index ] += FunctionAdapter::getValue( mesh,
+               ( *userData.fu )[ entity.getIndex() ] += FunctionAdapter::getValue( mesh,
                                                                        *userData.rightHandSide,
-                                                                       index,
+                                                                       entity,
                                                                        *userData.time );
             }
 
@@ -170,13 +170,11 @@ class tnlExplicitUpdater< tnlGrid< Dimensions, Real, Device, Index >,
             __cuda_callable__
             static inline void processEntity( const MeshType& mesh,
                                               TraverserUserData& userData,
-                                              const IndexType index,
                                               const GridEntity& entity )
             {
                userData.boundaryConditions->setBoundaryConditions
                ( *userData.time,
                  mesh,
-                 index,
                  entity,
                  *userData.u,
                  *userData.fu );
@@ -194,25 +192,21 @@ class tnlExplicitUpdater< tnlGrid< Dimensions, Real, Device, Index >,
             __cuda_callable__
             static inline void processEntity( const MeshType& mesh,
                                               TraverserUserData& userData,
-                                              const IndexType index,
                                               const EntityType& entity )
             {
-               ( *userData.fu)[ index ] = 
+               ( *userData.fu)[ entity.getIndex() ] = 
                   userData.differentialOperator->getValue(
                      mesh,
-                     index,
                      entity,
                      *userData.u,
                      *userData.time );
 
                typedef tnlFunctionAdapter< MeshType, RightHandSide > FunctionAdapter;
-               ( * userData.fu )[ index ] += 
+               ( * userData.fu )[ entity.getIndex() ] += 
                   FunctionAdapter::getValue(
-                     mesh,
                      *userData.rightHandSide,
-                     index,
                      entity,
-                    *userData.time );
+                     *userData.time );
             }
       };
 };
diff --git a/src/solvers/pde/tnlLinearSystemAssembler.h b/src/solvers/pde/tnlLinearSystemAssembler.h
index 6f5c43f9b6..66436f5a05 100644
--- a/src/solvers/pde/tnlLinearSystemAssembler.h
+++ b/src/solvers/pde/tnlLinearSystemAssembler.h
@@ -18,7 +18,7 @@
 #ifndef TNLLINEARSYSTEMASSEMBLER_H_
 #define TNLLINEARSYSTEMASSEMBLER_H_
 
-#include <functors/tnlFunctorAdapter.h>
+#include <functions/tnlFunctionAdapter.h>
 
 template< typename Real,
           typename DofVector,
@@ -215,14 +215,13 @@ class tnlLinearSystemAssembler< tnlGrid< Dimensions, Real, Device, Index >,
          __cuda_callable__
          static void processEntity( const MeshType& mesh,
                                     TraverserUserData& userData,
-                                    const IndexType index,
                                     const EntityType& entity )
          {
-             ( *userData.b )[ index ] = 0.0;           
+             ( *userData.b )[ entity.getIndex() ] = 0.0;           
              userData.boundaryConditions->updateLinearSystem
                ( *userData.time + *userData.tau,
                  mesh,
-                 index,
+                 entity.getIndex(),
                  entity,
                  *userData.u,
                  *userData.b,
@@ -238,15 +237,14 @@ class tnlLinearSystemAssembler< tnlGrid< Dimensions, Real, Device, Index >,
          __cuda_callable__
          static void processEntity( const MeshType& mesh,
                                     TraverserUserData& userData,
-                                    const IndexType index,
                                     const EntityType& entity )
          {
-            ( *userData.b )[ index ] = 0.0;            
+            ( *userData.b )[ entity.getIndex() ] = 0.0;            
             userData.differentialOperator->updateLinearSystem
                ( *userData.time,
                  *userData.tau,
                  mesh,
-                 index,
+                 entity.getIndex(),
                  entity,
                  *userData.u,
                  *userData.b,
@@ -254,15 +252,13 @@ class tnlLinearSystemAssembler< tnlGrid< Dimensions, Real, Device, Index >,
             
             typedef tnlFunctionAdapter< MeshType, RightHandSide > FunctionAdapter;
             const RealType& rhs = FunctionAdapter::getValue
-               ( mesh,
-                 *userData.rightHandSide,
-                 index,
+               ( *userData.rightHandSide,
                  entity,
                  *userData.time );
             TimeDiscretisation::applyTimeDiscretisation( *userData.matrix,
-                                                         ( *userData.b )[ index ],
-                                                         index,
-                                                         ( *userData.u )[ index ],
+                                                         ( *userData.b )[ entity.getIndex() ],
+                                                         entity.getIndex(),
+                                                         ( *userData.u )[ entity.getIndex() ],
                                                          ( *userData.tau ),
                                                          rhs );
          }
diff --git a/tests/unit-tests/operators/diffusion/tnlLinearDiffusionTest.cpp b/tests/unit-tests/operators/diffusion/tnlLinearDiffusionTest.cpp
index e20eba5283..687e485c2a 100644
--- a/tests/unit-tests/operators/diffusion/tnlLinearDiffusionTest.cpp
+++ b/tests/unit-tests/operators/diffusion/tnlLinearDiffusionTest.cpp
@@ -25,7 +25,7 @@
 #include <operators/diffusion/tnlLinearDiffusion.h>
 #include <operators/diffusion/tnlExactLinearDiffusion.h>
 #include "../tnlPDEOperatorEocTestResult.h"
-#include <functors/tnlExpBumpFunction.h>
+#include <functions/tnlExpBumpFunction.h>
 
 template< int Dimensions,
           typename Real,
diff --git a/tests/unit-tests/operators/tnlPDEOperatorEocTestSetter.h b/tests/unit-tests/operators/tnlPDEOperatorEocTestSetter.h
index 69d1c31843..16542c201d 100644
--- a/tests/unit-tests/operators/tnlPDEOperatorEocTestSetter.h
+++ b/tests/unit-tests/operators/tnlPDEOperatorEocTestSetter.h
@@ -19,7 +19,7 @@
 #define TNLPDEOPERATOREOCTESTSETTER_H_
 
 #include <mesh/tnlGrid.h>
-#include <functors/tnlExpBumpFunction.h>
+#include <functions/tnlExpBumpFunction.h>
 
 template< typename ApproximateOperator,
           typename ExactOperator,
diff --git a/tests/unit-tests/tnlApproximationError.h b/tests/unit-tests/tnlApproximationError.h
index 13b9bc2da8..6210c50594 100644
--- a/tests/unit-tests/tnlApproximationError.h
+++ b/tests/unit-tests/tnlApproximationError.h
@@ -19,8 +19,8 @@
 #define TNLAPPROXIMATIONERROR_H_
 
 #include <mesh/tnlGrid.h>
-#include <functors/tnlConstantFunction.h>
-#include <operators/tnlAnalyticDirichletBoundaryConditions.h>
+#include <functions/tnlConstantFunction.h>
+#include <operators/tnlDirichletBoundaryConditions.h>
 #include <solvers/pde/tnlExplicitUpdater.h>
 
 class tnlExplicitApproximation
@@ -66,7 +66,7 @@ class tnlApproximationError< Mesh, ExactOperator, ApproximateOperator, Function,
       typedef typename MeshType::IndexType IndexType;
       typedef typename MeshType::VertexType VertexType;
       typedef tnlConstantFunction< MeshType::meshDimensions, RealType > ConstantFunctionType;
-      typedef tnlAnalyticDirichletBoundaryConditions< MeshType, Function  > BoundaryConditionsType;
+      typedef tnlDirichletBoundaryConditions< MeshType, Function  > BoundaryConditionsType;
 
       static void getError( const Mesh& mesh,
                             const ExactOperator& exactOperator,
@@ -91,7 +91,7 @@ class tnlApproximationError< Mesh, ExactOperator, ApproximateOperator, Function,
       typedef typename MeshType::IndexType IndexType;
       typedef typename MeshType::VertexType VertexType;
       typedef tnlConstantFunction< MeshType::meshDimensions, RealType > ConstantFunctionType;
-      typedef tnlAnalyticDirichletBoundaryConditions< MeshType, Function  > BoundaryConditionsType;
+      typedef tnlDirichletBoundaryConditions< MeshType, Function  > BoundaryConditionsType;
 
       static void getError( const Mesh& mesh,
                             const ExactOperator& exactOperator,
diff --git a/tests/unit-tests/tnlApproximationError_impl.h b/tests/unit-tests/tnlApproximationError_impl.h
index 692b75e9ba..f171e2ffb5 100644
--- a/tests/unit-tests/tnlApproximationError_impl.h
+++ b/tests/unit-tests/tnlApproximationError_impl.h
@@ -20,7 +20,7 @@
 
 #include <mesh/tnlTraverser.h>
 #include <core/vectors/tnlVector.h>
-#include <functors/tnlFunctionDiscretizer.h>
+#include <functions/tnlFunctionDiscretizer.h>
 #include <matrices/tnlCSRMatrix.h>
 #include <matrices/tnlMatrixSetter.h>
 #include <solvers/pde/tnlLinearSystemAssembler.h>
diff --git a/tools/src/tnl-init.cpp b/tools/src/tnl-init.cpp
index f71d0e41cc..d3141c87dc 100644
--- a/tools/src/tnl-init.cpp
+++ b/tools/src/tnl-init.cpp
@@ -21,7 +21,7 @@
 #include <debug/tnlDebug.h>
 #include <config/tnlConfigDescription.h>
 #include <config/tnlParameterContainer.h>
-#include <functors/tnlTestFunction.h>
+#include <functions/tnlTestFunction.h>
 #include <mesh/tnlDummyMesh.h>
 #include <mesh/tnlGrid.h>
 
diff --git a/tools/src/tnl-init.h b/tools/src/tnl-init.h
index 427c4eb102..57bca2bfd0 100644
--- a/tools/src/tnl-init.h
+++ b/tools/src/tnl-init.h
@@ -21,8 +21,8 @@
 #include <config/tnlParameterContainer.h>
 #include <core/vectors/tnlVector.h>
 #include <mesh/tnlGrid.h>
-#include <functors/tnlFunctionDiscretizer.h>
-#include <functors/tnlTestFunction.h>
+#include <functions/tnlFunctionDiscretizer.h>
+#include <functions/tnlTestFunction.h>
 #include <operators/tnlFiniteDifferences.h>
 #include <core/mfilename.h>
 
-- 
GitLab