From d294aca01c1e72c1f2abf80259906d58aaeae50e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ond=C5=99ej=20Sz=C3=A9kely?= <ondra.szekely@gmail.com>
Date: Fri, 2 May 2014 01:46:23 +0200
Subject: [PATCH] Fixing errors.

---
 examples/heat-equation/CMakeLists.txt         |   2 +-
 examples/heat-equation/Makefile               |   8 +-
 examples/heat-equation/heat-equation-conf.h   |   2 +-
 .../heat-equation/heatEquationSolver_impl.h   |   4 +-
 examples/heat-equation/main.cpp               |   6 +-
 examples/heat-equation/run-init-mesh          |   2 +-
 examples/heat-equation/run-simple-solver      |  21 --
 examples/heat-equation/simple-solver.cfg.desc |  42 ----
 examples/heat-equation/simpleProblemSetter.h  |  63 -----
 .../heat-equation/simpleProblemSetter_impl.h  |  79 ------
 examples/heat-equation/simpleProblemSolver.h  |  81 -------
 .../heat-equation/simpleProblemSolver_impl.h  | 227 ------------------
 examples/heat-equation/tnlAnalyticSolution.h  |  25 +-
 .../heat-equation/tnlAnalyticSolution_impl.h  | 101 +++++++-
 examples/heat-equation/tnlLinearDiffusion.h   |  25 +-
 .../heat-equation/tnlLinearDiffusion_impl.h   |  22 +-
 examples/heat-equation/tnlTimeFunction.h      |  50 +++-
 examples/heat-equation/tnlTimeFunction_impl.h |  61 ++++-
 18 files changed, 230 insertions(+), 591 deletions(-)
 delete mode 100644 examples/heat-equation/run-simple-solver
 delete mode 100644 examples/heat-equation/simple-solver.cfg.desc
 delete mode 100644 examples/heat-equation/simpleProblemSetter.h
 delete mode 100644 examples/heat-equation/simpleProblemSetter_impl.h
 delete mode 100644 examples/heat-equation/simpleProblemSolver.h
 delete mode 100644 examples/heat-equation/simpleProblemSolver_impl.h

diff --git a/examples/heat-equation/CMakeLists.txt b/examples/heat-equation/CMakeLists.txt
index 7b4491c17e..8a5f18cc47 100755
--- a/examples/heat-equation/CMakeLists.txt
+++ b/examples/heat-equation/CMakeLists.txt
@@ -8,7 +8,7 @@ INSTALL( FILES Makefile
                tnlNeumannBoundaryConditions_impl.h
                tnlNeumannBoundaryConditions.h
                tnlLinearDiffusion.h
-	       tnlLinearDiffusion_impl.h
+               tnlLinearDiffusion_impl.h
                tnlDirichletBoundaryConditions_impl.h
                tnlDirichletBoundaryConditions.h
                tnlAnalyticSolution_impl.h
diff --git a/examples/heat-equation/Makefile b/examples/heat-equation/Makefile
index 259230476c..5918c5e8d3 100644
--- a/examples/heat-equation/Makefile
+++ b/examples/heat-equation/Makefile
@@ -2,7 +2,7 @@ TNL_VERSION=0.1
 TNL_INSTALL_DIR=${HOME}/local/lib
 TNL_INCLUDE_DIR=${HOME}/local/include/tnl-${TNL_VERSION}
 
-TARGET = simple-solver
+TARGET = heat-equation
 CONFIG_FILE = $(TARGET).cfg.desc
 INSTALL_DIR = ${HOME}/local
 CXX = g++
@@ -26,11 +26,11 @@ dist: $(DIST)
 
 install: $(TARGET)
 	cp $(TARGET) $(INSTALL_DIR)/bin
-	cp $(CONFIG_FILE) $(INSTALL_DIR)/share/tnl-0.1/examples/simple-solver
+	cp $(CONFIG_FILE) $(INSTALL_DIR)/share/tnl-0.1/examples/heat-equation
 
 uninstall: $(TARGET)
 	rm -f $(INSTALL_DIR)/bin/$(TARGET) 
-	rm -f $(CONFIG_FILE) $(INSTALL_DIR)/share/tnl-0.1/examples/simple-solver
+	rm -f $(CONFIG_FILE) $(INSTALL_DIR)/share/tnl-0.1/examples/heat-equation
 
 $(TARGET): $(TARGET)-conf.h $(OBJECTS)
 	$(CXX) -o $(TARGET) $(OBJECTS) $(LD_FLAGS)
@@ -39,5 +39,5 @@ $(TARGET): $(TARGET)-conf.h $(OBJECTS)
 	$(CXX) -c -o $@ $(CXX_FLAGS) $<
 
 $(TARGET)-conf.h:
-	echo "#define CONFIG_FILE \"${INSTALL_DIR}/share/tnl-0.1/examples/simple-solver/${CONFIG_FILE}\" " > $(TARGET)-conf.h 
+	echo "#define CONFIG_FILE \"${INSTALL_DIR}/share/tnl-0.1/examples/heat-equation/${CONFIG_FILE}\" " > $(TARGET)-conf.h 
 
diff --git a/examples/heat-equation/heat-equation-conf.h b/examples/heat-equation/heat-equation-conf.h
index 97a4216ca2..020c206e26 100644
--- a/examples/heat-equation/heat-equation-conf.h
+++ b/examples/heat-equation/heat-equation-conf.h
@@ -1 +1 @@
-#define CONFIG_FILE "/home/ondra/local/share/tnl-0.1/examples/simple-solver/heat-equation.cfg.desc" 
+#define CONFIG_FILE "/home/ondra/local/share/tnl-0.1/examples/heat-equation/heat-equation.cfg.desc" 
diff --git a/examples/heat-equation/heatEquationSolver_impl.h b/examples/heat-equation/heatEquationSolver_impl.h
index 55c751ac88..9657be6104 100644
--- a/examples/heat-equation/heatEquationSolver_impl.h
+++ b/examples/heat-equation/heatEquationSolver_impl.h
@@ -137,7 +137,7 @@ bool heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide,TimeFunc
       if( ! this -> numericalLaplace. save( fileName ) )
          return false;
       
-      exit(0);
+      EXIT_SUCCESS;
    }
    
    return true;
@@ -168,6 +168,8 @@ void heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide,TimeFunc
 
    if( DeviceType :: getDevice() == tnlHostDevice )
    {
+      boundaryCondition.applyBoundaryTimeDerivation(mesh, _fu, time, timeFunction, analyticSpaceFunction);
+      
       boundaryCondition.applyBoundaryConditions(mesh, _u, time, timeFunction, analyticSpaceFunction);
       
       diffusion.getExplicitRHS(mesh,_u,_fu);  
diff --git a/examples/heat-equation/main.cpp b/examples/heat-equation/main.cpp
index c36c7ac0d4..756220ea9d 100644
--- a/examples/heat-equation/main.cpp
+++ b/examples/heat-equation/main.cpp
@@ -15,13 +15,13 @@
  *                                                                         *
  ***************************************************************************/
 
-#include "simple-solver-conf.h"
-#include "simpleProblemSetter.h"
+#include "heat-equation-conf.h"
+#include "heatEquationSetter.h"
 #include <solvers/tnlSolver.h>
 
 int main( int argc, char* argv[] )
 {
-   tnlSolver< simpleProblemSetter > solver;
+   tnlSolver< heatEquationSetter > solver;
    if( ! solver. run( CONFIG_FILE, argc, argv ) )
       return EXIT_FAILURE;
    return EXIT_SUCCESS;
diff --git a/examples/heat-equation/run-init-mesh b/examples/heat-equation/run-init-mesh
index 24ffb2f14e..1de910039e 100644
--- a/examples/heat-equation/run-init-mesh
+++ b/examples/heat-equation/run-init-mesh
@@ -1,6 +1,6 @@
 #!/bin/bash
 
-dofSize=32
+dofSize=4
 dimension=1
 
 analyticFunction="sin-wave"
diff --git a/examples/heat-equation/run-simple-solver b/examples/heat-equation/run-simple-solver
deleted file mode 100644
index ee7d7234e6..0000000000
--- a/examples/heat-equation/run-simple-solver
+++ /dev/null
@@ -1,21 +0,0 @@
-#!/bin/bash
-
-tnl-grid-setup --dimensions 2 \
-               --origin-x 0.0 \
-               --origin-y 0.0 \
-               --proportions-x 1.0 \
-               --proportions-y 1.0 \
-               --size-x 100 \
-               --size-y 100
-               
-tnl-discrete --function sin-waves \
-             --output-file u-ini.tnl               
-
-simple-solver --dimensions 2 \
-              --time-discretisation explicit \
-              --discrete-solver merson \
-              --snapshot-period 0.01 \
-              --final-time 1.0
-              
-tnl-view --mesh mesh.tnl *tnl              
-              
\ No newline at end of file
diff --git a/examples/heat-equation/simple-solver.cfg.desc b/examples/heat-equation/simple-solver.cfg.desc
deleted file mode 100644
index 6999b1df04..0000000000
--- a/examples/heat-equation/simple-solver.cfg.desc
+++ /dev/null
@@ -1,42 +0,0 @@
-group Problem
-{
-   integer dimensions(!)                  [Number of dimenions the problem is defined in.];
-   string problem-name( "simple" )        [This defines particular problem.];
-   string log-file("simple-problem.log")  [File name for the log.];
-},[The problem parameters.];
-
-group Discretisation
-{
-   string mesh( "mesh.tnl" )            [A file containing the numerical mesh.];
-   string initial-condition("u-ini.tnl");
-   string time-discretisation(!)     [Time discretisation for the time dependent problems. Can be explicit, semi-implicit or fully-implicit.];
-},[The numerical scheme set-up.];
-
-group Solver
-{
-   string real-type( "double" )                     [Precision of the floating point arithmetics. Can be float, double, long-double.];
-   string index-type( "int" )                       [Indexing type for arrays, vectors, matrices etc. Can be int or long-int.];
-   string device( "host" )                          [Device to use for the computations. Can be host or cuda.];
-   string discrete-solver(!)                        [The main solver for the discretised problem.];
-   real merson-adaptivity(1.0e-4)                   [Parameter controling adaptivity of the Runge-Kutta-Merson solver.];
-   string time-function("time-independent")         [Time function.];
-   string analytic-space-function(!)                [Analytic space solution.];
-   real amplitude(1.0)                              [Amplitude of boundary and analytic function.];
-   real wave-length(1.0)                            [Wave length of 1D boundary and 1D analytic function.];
-   real wave-length-x(1.0)                          [Wave length in x-axis of boundary and analytic function.];
-   real wave-length-y(1.0)                          [Wave length in y-axis of boundary and analytic function.];
-   real wave-length-z(1.0)                          [Wave length in z-axis of boundary and analytic function.];
-   real waves-number(0.0)                           [Number of waves];
-   real waves-number-x(0.0)                         [Number of waves in x-axis];
-   real waves-number-y(0.0)                         [Number of waves in y-axis];
-   real waves-number-z(0.0)                         [Number of waves in z-axis];
-   real phase(0.0)                                  [Phase of 1D boundary and 1D analytic function.];
-   real phase-x(0.0)                                [Phase in x-axis of boundary and analytic function.];
-   real phase-y(0.0)                                [Phase in y-axis of boundary and analytic function.];
-   real phase-z(0.0)                                [Phase in z-axis of boundary and analytic function.];
-   real sigma(1.0)                                  [Sigma.];
-   real initial-tau(1.0e-5)                         [Initial time step for the time dependent PDE solver. It can be changed if an adaptive solver is used.];
-   real snapshot-period(!)                          [Time period for writing a state of the time dependent problem to a file.];
-   real final-time(!)                               [Stop time of the time dependent simulation.];
-   integer verbose(1)                               [Set the verbose mode. The higher number the more messages are generated.];
-},[The solver set-up.];
diff --git a/examples/heat-equation/simpleProblemSetter.h b/examples/heat-equation/simpleProblemSetter.h
deleted file mode 100644
index 22484bc815..0000000000
--- a/examples/heat-equation/simpleProblemSetter.h
+++ /dev/null
@@ -1,63 +0,0 @@
-/***************************************************************************
-                          simpleProblemSetter.h  -  description
-                             -------------------
-    begin                : Feb 23, 2013
-    copyright            : (C) 2013 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/***************************************************************************
- *                                                                         *
- *   This program is free software; you can redistribute it and/or modify  *
- *   it under the terms of the GNU General Public License as published by  *
- *   the Free Software Foundation; either version 2 of the License, or     *
- *   (at your option) any later version.                                   *
- *                                                                         *
- ***************************************************************************/
-
-#ifndef SIMPLEPROBLEMTYPESSETTER_H_
-#define SIMPLEPROBLEMTYPESSETTER_H_
-
-#include <config/tnlParameterContainer.h>
-#include <mesh/tnlGrid.h>
-#include "simpleProblemSolver.h"
-#include <generators/functions/tnlSinWaveFunction.h>
-#include <generators/functions/tnlExpBumpFunction.h>
-#include <generators/functions/tnlSinBumpsFunction.h>
-#include "tnlTimeFunction.h"
-#include "tnlDirichletBoundaryConditions.h"
-#include "tnlLinearDiffusion.h"
-#include "tnlNeumannBoundaryConditions.h"
-#include "tnlZeroRightHandSide.h"
-#include "tnlRightHandSide.h"
-
-   
-template< typename MeshType,
-          typename SolverStarter >
-class simpleProblemSetter
-{
-   public:
- 
-   typedef TimeFunction<MeshType,TimeFunctionBase::TimeIndependent> TimeIndependent;
-   typedef TimeFunction<MeshType,TimeFunctionBase::Linear> Linear;
-   typedef TimeFunction<MeshType,TimeFunctionBase::Quadratic> Quadratic;
-   typedef TimeFunction<MeshType,TimeFunctionBase::Cosinus> Cosinus;
-   typedef typename MeshType::RealType RealType;
-   typedef tnlStaticVector<MeshType::Dimensions, RealType> Vertex;
-
-     
-   template< typename RealType, typename DeviceType, typename IndexType, typename TimeFunction>
-   static bool setAnalyticSpaceFunction (const tnlParameterContainer& parameters);  
-    
-   template< typename RealType, typename DeviceType, typename IndexType>
-   static bool setTimeFunction (const tnlParameterContainer& parameters);
-      
-   template< typename RealType,
-             typename DeviceType,
-             typename IndexType >
-   static bool run( const tnlParameterContainer& parameters );
-};
-
-#include "simpleProblemSetter_impl.h"
-
-#endif /* SIMPLEPROBLEMSETTER_H_ */
diff --git a/examples/heat-equation/simpleProblemSetter_impl.h b/examples/heat-equation/simpleProblemSetter_impl.h
deleted file mode 100644
index 6ccc30e28b..0000000000
--- a/examples/heat-equation/simpleProblemSetter_impl.h
+++ /dev/null
@@ -1,79 +0,0 @@
-/***************************************************************************
-                          simpleProblemSetter_impl.h  -  description
-                             -------------------
-    begin                : Mar 10, 2013
-    copyright            : (C) 2013 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/***************************************************************************
- *                                                                         *
- *   This program is free software; you can redistribute it and/or modify  *
- *   it under the terms of the GNU General Public License as published by  *
- *   the Free Software Foundation; either version 2 of the License, or     *
- *   (at your option) any later version.                                   *
- *                                                                         *
- ***************************************************************************/
-
-#ifndef SIMPLEPROBLEMSETTER_IMPL_H_
-#define SIMPLEPROBLEMSETTER_IMPL_H_
- 
-#include "simpleProblemSetter.h"
-
-template< typename MeshType, typename SolverStarter >
-template< typename RealType, typename DeviceType, typename IndexType, typename TimeFunction>
-bool simpleProblemSetter< MeshType, SolverStarter > ::setAnalyticSpaceFunction (const tnlParameterContainer& parameters)
-{
-   SolverStarter solverStarter;
-   
-   const tnlString& analyticSpaceFunctionParameter = parameters.GetParameter<tnlString>("analytic-space-function");
-   
-   if (analyticSpaceFunctionParameter == "sin-wave")
-      return solverStarter.run< simpleProblemSolver< MeshType,tnlLinearDiffusion<MeshType>,
-                                 tnlDirichletBoundaryConditions<MeshType>,tnlRightHandSide<MeshType>,
-                                 TimeFunction,tnlSinWaveFunction<MeshType::Dimensions,Vertex,DeviceType>>>(parameters);
-   if (analyticSpaceFunctionParameter == "sin-bumps")
-     return solverStarter.run< simpleProblemSolver< MeshType,tnlLinearDiffusion<MeshType>,
-                               tnlDirichletBoundaryConditions<MeshType>,tnlRightHandSide<MeshType>,
-                               TimeFunction, tnlSinBumpsFunction<MeshType::Dimensions,Vertex,DeviceType>>>(parameters);
-   if (analyticSpaceFunctionParameter == "exp-bump")
-      return solverStarter.run< simpleProblemSolver< MeshType,tnlLinearDiffusion<MeshType>,
-                                tnlDirichletBoundaryConditions<MeshType>,tnlRightHandSide<MeshType>,
-                                TimeFunction, tnlExpBumpFunction<MeshType::Dimensions,Vertex,DeviceType>>>(parameters);
-   
-   cerr<<"Unknown analytic-space-function parameter: "<<analyticSpaceFunctionParameter<<". ";
-   return 0;
-}
-
-template< typename MeshType, typename SolverStarter >
-template< typename RealType, typename DeviceType, typename IndexType>
-bool simpleProblemSetter< MeshType, SolverStarter > ::setTimeFunction (const tnlParameterContainer& parameters)
-{
-   const tnlString& timeFunctionParameter = parameters.GetParameter<tnlString>("time-function");
-   
-   if (timeFunctionParameter == "time-independent")
-      return setAnalyticSpaceFunction<RealType, DeviceType, IndexType, TimeIndependent>(parameters);
-   if (timeFunctionParameter == "linear")
-      return setAnalyticSpaceFunction<RealType, DeviceType, IndexType, Linear>(parameters);
-   if (timeFunctionParameter == "quadratic")
-      return setAnalyticSpaceFunction<RealType, DeviceType, IndexType, Quadratic>(parameters);
-   if (timeFunctionParameter == "cosinus")
-      return setAnalyticSpaceFunction<RealType, DeviceType, IndexType, Cosinus>(parameters);
-   
-   cerr<<"Unknown time-function parameter: "<<timeFunctionParameter<<". ";
-   return 0;
-}
-
-
-template< typename MeshType,
-          typename SolverStarter >
-   template< typename RealType,
-             typename DeviceType,
-             typename IndexType >
-bool simpleProblemSetter< MeshType, SolverStarter > :: run( const tnlParameterContainer& parameters )
-{
-   return setTimeFunction<RealType, DeviceType, IndexType>(parameters);
-}
-
-
-#endif /* SIMPLEPROBLEMSETTER_IMPL_H_ */
diff --git a/examples/heat-equation/simpleProblemSolver.h b/examples/heat-equation/simpleProblemSolver.h
deleted file mode 100644
index 6babb18618..0000000000
--- a/examples/heat-equation/simpleProblemSolver.h
+++ /dev/null
@@ -1,81 +0,0 @@
-/***************************************************************************
-                          simpleProblemSolver.h  -  description
-                             -------------------
-    begin                : Feb 23, 2013
-    copyright            : (C) 2013 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/***************************************************************************
- *                                                                         *
- *   This program is free software; you can redistribute it and/or modify  *
- *   it under the terms of the GNU General Public License as published by  *
- *   the Free Software Foundation; either version 2 of the License, or     *
- *   (at your option) any later version.                                   *
- *                                                                         *
- ***************************************************************************/
-
-#ifndef SIMPLEPROBLEMSOLVER_H_
-#define SIMPLEPROBLEMSOLVER_H_
-
-#include <matrices/tnlCSRMatrix.h>
-#include <solvers/preconditioners/tnlDummyPreconditioner.h>
-#include <solvers/tnlSolverMonitor.h>
-#include <core/tnlLogger.h>
-#include <core/vectors/tnlVector.h>
-#include <core/vectors/tnlSharedVector.h>
-#include "simpleProblemSolver.h"
-#include "tnlAnalyticSolution.h"
-
-
-template< typename Mesh, typename Diffusion, typename BoundaryCondition, typename RightHandSide, typename TimeFunction, typename AnalyticSpaceFunction>
-class simpleProblemSolver
-{
-   public:
-
-   typedef typename Mesh :: RealType RealType;
-   typedef typename Mesh :: DeviceType DeviceType;
-   typedef typename Mesh :: IndexType IndexType;
-   typedef Mesh MeshType;
-   typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
-   typedef tnlCSRMatrix< RealType, DeviceType, IndexType > DiscreteSolverMatrixType;
-   typedef tnlDummyPreconditioner< RealType, DeviceType, IndexType > DiscreteSolverPreconditioner;
-
-   static tnlString getTypeStatic();
-
-   tnlString getPrologHeader() const;
-
-   void writeProlog( tnlLogger& logger,
-                     const tnlParameterContainer& parameters ) const;
-
-   bool init( const tnlParameterContainer& parameters );
-
-   bool setInitialCondition( const tnlParameterContainer& parameters );
-
-   bool makeSnapshot( const RealType& time, const IndexType& step );
-
-   DofVectorType& getDofVector();
-
-   void GetExplicitRHS( const RealType& time,
-                        const RealType& tau,
-                        DofVectorType& _u,
-                        DofVectorType& _fu );
-
-   tnlSolverMonitor< RealType, IndexType >* getSolverMonitor();
-   
-   protected:
-
-   DofVectorType dofVector,dofVector2;
-   tnlSharedVector< RealType, DeviceType, IndexType > u,v;
-   MeshType mesh;
-   AnalyticSpaceFunction analyticSpaceFunction;
-   TimeFunction timeFunction;
-   AnalyticSolution<MeshType> analyticSolution;
-   BoundaryCondition boundaryCondition;
-   Diffusion diffusion;
-   RightHandSide RHS;
-};
-
-#include "simpleProblemSolver_impl.h"
-
-#endif /* SIMPLEPROBLEM_H_ */
diff --git a/examples/heat-equation/simpleProblemSolver_impl.h b/examples/heat-equation/simpleProblemSolver_impl.h
deleted file mode 100644
index 17060502bc..0000000000
--- a/examples/heat-equation/simpleProblemSolver_impl.h
+++ /dev/null
@@ -1,227 +0,0 @@
-/***************************************************************************
-                          simpleProblemSolver_impl.h  -  description
-                             -------------------
-    begin                : Mar 10, 2013
-    copyright            : (C) 2013 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/***************************************************************************
- *                                                                         *
- *   This program is free software; you can redistribute it and/or modify  *
- *   it under the terms of the GNU General Public License as published by  *
- *   the Free Software Foundation; either version 2 of the License, or     *
- *   (at your option) any later version.                                   *
- *                                                                         *
- ***************************************************************************/
-
-#ifndef SIMPLEPROBLEMSOLVER_IMPL_H_
-#define SIMPLEPROBLEMSOLVER_IMPL_H_
-
-#include <core/mfilename.h>
-#include "simpleProblemSolver.h"
-
-template< typename Mesh, typename Diffusion, typename BoundaryCondition, typename RightHandSide, 
-          typename TimeFunction, typename AnalyticSpaceFunction>
-tnlString simpleProblemSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide,TimeFunction,AnalyticSpaceFunction> 
-::getTypeStatic()
-{
-   /****
-    * Replace 'simpleProblemSolver' by the name of your solver.
-    */
-   return tnlString( "simpleProblemSolver< " ) + Mesh :: getTypeStatic() + " >";
-}
-
-template< typename Mesh, typename Diffusion, typename BoundaryCondition, typename RightHandSide,
-          typename TimeFunction, typename AnalyticSpaceFunction>
-tnlString simpleProblemSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide,TimeFunction,AnalyticSpaceFunction>
-:: getPrologHeader() const
-{
-   /****
-    * Replace 'Simple Problem' by the your desired title in the log table.
-    */
-   return tnlString( "Simple Problem" );
-}
-
-template< typename Mesh, typename Diffusion, typename BoundaryCondition, typename RightHandSide,
-          typename TimeFunction, typename AnalyticSpaceFunction>
-void simpleProblemSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide,TimeFunction,AnalyticSpaceFunction>
-:: writeProlog( tnlLogger& logger, const tnlParameterContainer& parameters ) const
-{
-   /****
-    * In prolog, write all input parameters which define the numerical simulation.
-    * Use methods:
-    *
-    *    logger. writeParameters< Type >( "Label:", "name", parameters );
-    *
-    *  or
-    *
-    *    logger. writeParameter< Type >( "Label:", value );
-    *
-    *  See tnlLogger.h for more details.
-    */
-
-   logger. WriteParameter< tnlString >( "Problem name:", "problem-name", parameters );
-   logger. WriteParameter< int >( "Simple parameter:", 1 );
-}
-
-template< typename Mesh, typename Diffusion, typename BoundaryCondition, typename RightHandSide,
-          typename TimeFunction, typename AnalyticSpaceFunction>
-bool simpleProblemSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide,TimeFunction,AnalyticSpaceFunction>
-:: init( const tnlParameterContainer& parameters )
-{
-   /****
-    * Set-up your solver here. It means:
-    * 1. Read input parameters and model coefficients like these
-    */
-       
-   analyticSpaceFunction.init(parameters);
-
-   /****
-    * 2. Set-up geometry of the problem domain using some mesh like tnlGrid.
-    * Implement additional template specializations of the method initMesh
-    * if necessary.
-    */
-   const tnlString& meshFile = parameters.GetParameter< tnlString >( "mesh" );
-   if( ! this->mesh.load( meshFile ) )
-   {
-      cerr << "I am not able to load the mesh from the file " << meshFile << "." << endl;
-      return false;
-   }
-   
-   /****
-    * 3. Set-up DOFs and supporting grid functions
-    */
-   const IndexType& dofs = this->mesh.getDofs();
-   dofVector. setSize(dofs);
-   dofVector2. setSize(dofs);
-
-   /****
-    * You may use tnlSharedVector if you need to split the dofVector into more
-    * grid functions like the following example:
-    */
-   this -> u. bind( & dofVector. getData()[ 0 * dofs ], dofs );
-   this -> v. bind( & dofVector2. getData()[ 0 * dofs ], dofs );
-
-   /****
-    * You may now treat u and v as usual vectors and indirectly work with this->dofVector.
-    */
-
-   return true;
-}
-
-template< typename Mesh, typename Diffusion, typename BoundaryCondition, typename RightHandSide,
-          typename TimeFunction, typename AnalyticSpaceFunction>
-bool simpleProblemSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide,TimeFunction,AnalyticSpaceFunction>
-:: setInitialCondition( const tnlParameterContainer& parameters )
-{
-   /****
-    * Set the initial condition here. Manipulate only this -> dofVector.
-    */
-   const tnlString& initialConditionFile = parameters.GetParameter< tnlString >( "initial-condition" );
-   if( ! this->u.load( initialConditionFile ) )
-   {
-      cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." << endl;
-      return false;
-   }
-   
-   boundaryCondition.applyBoundaryConditions(mesh,u,0.0,timeFunction,analyticSpaceFunction);
-   
-   return true;
-}
-
-template< typename Mesh, typename Diffusion, typename BoundaryCondition, typename RightHandSide,
-          typename TimeFunction, typename AnalyticSpaceFunction>
-bool simpleProblemSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide,TimeFunction,AnalyticSpaceFunction> 
-:: makeSnapshot( const RealType& time, const IndexType& step )
-{
-   /****
-    * Use this method to write state of the solver to file(s).
-    * All data are stored in this -> dofVector. You may use
-    * supporting vectors and bind them with the dofVector as before.
-    */
-   cout << endl << "Writing output at time " << time << " step " << step << "." << endl;
-   
-   /****
-    * Now write them to files.
-    */
-   
-   analyticSolution.compute(mesh,time,v,u,timeFunction,analyticSpaceFunction);
-   
-   tnlString fileName;
-   FileNameBaseNumberEnding( "u-", step, 5, ".tnl", fileName );
-   if( ! this -> u. save( fileName ) )
-      return false;
-   
-   FileNameBaseNumberEnding( "v-", step, 5, ".tnl", fileName );
-   if( ! this -> v. save( fileName ) )
-      return false;
-
-   if(time == 1.0)
-      return true;
-   
-   return true;
-}
-
-template< typename Mesh, typename Diffusion, typename BoundaryCondition, typename RightHandSide,
-          typename TimeFunction, typename AnalyticSpaceFunction>
-typename simpleProblemSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide,TimeFunction,AnalyticSpaceFunction>
-:: DofVectorType& simpleProblemSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide,TimeFunction,AnalyticSpaceFunction> 
-:: getDofVector()
-{
-   /****
-    * You do not need to change this usually.
-    */
-   return dofVector;
-}
-
-template< typename Mesh, typename Diffusion, typename BoundaryCondition, typename RightHandSide,
-          typename TimeFunction, typename AnalyticSpaceFunction>
-void simpleProblemSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide,TimeFunction,AnalyticSpaceFunction> 
-:: GetExplicitRHS( const RealType& time, const RealType& tau, DofVectorType& _u, DofVectorType& _fu )
-{
-   /****
-    * If you use an explicit solver like tnlEulerSolver or tnlMersonSolver, you
-    * need to implement this method. Compute the right-hand side of
-    *
-    *   d/dt u(x) = fu( x, u )
-    *
-    * You may use supporting vectors again if you need.
-    */
-
-   if( DeviceType :: getDevice() == tnlHostDevice )
-   {
-      
-      /****
-       *  Write the host solver here.
-       */
-      boundaryCondition.applyBoundaryConditions(mesh, _u, time, timeFunction, analyticSpaceFunction);
-      
-      diffusion.getExplicitRHS(mesh,time,tau,_u,_fu);  
-      
-      boundaryCondition.applyBoundaryTimeDerivation(mesh, _fu, time, timeFunction, analyticSpaceFunction);
-      
-      RHS.applyRHSValues(mesh, time, _fu, timeFunction, analyticSpaceFunction);
-      
-   }
-#ifdef HAVE_CUDA
-   if( DeviceType :: getDevice() == tnlCudaDevice )
-   {
-      /****
-       * Write the CUDA solver here.
-       */
-   }
-#endif
-}
-
-template< typename Mesh, typename Diffusion, typename BoundaryCondition, typename RightHandSide, 
-          typename TimeFunction, typename AnalyticSpaceFunction>
-tnlSolverMonitor< typename simpleProblemSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide,TimeFunction,AnalyticSpaceFunction>:: RealType,
-                  typename simpleProblemSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide,TimeFunction,AnalyticSpaceFunction> :: IndexType >*
-                  simpleProblemSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide,TimeFunction,AnalyticSpaceFunction> 
-::  getSolverMonitor()
-{
-   return 0;
-}
-
-#endif /* SIMPLEPROBLEM_IMPL_H_ */
diff --git a/examples/heat-equation/tnlAnalyticSolution.h b/examples/heat-equation/tnlAnalyticSolution.h
index ef6dc3f2ac..2b1ae652bc 100644
--- a/examples/heat-equation/tnlAnalyticSolution.h
+++ b/examples/heat-equation/tnlAnalyticSolution.h
@@ -22,8 +22,12 @@ class AnalyticSolution<tnlGrid<1,Real,Device,Index,tnlIdenticalGridGeometry>>
    typedef tnlStaticVector< 1, RealType > VertexType;
    
    template< typename TimeFunction, typename AnalyticSpaceFunction>
-   void compute(const MeshType& mesh, const RealType& time, SharedVector& output, SharedVector& numerical, 
-                TimeFunction timeFunction, AnalyticSpaceFunction anlayticSpaceFunction);
+   void computeAnalyticSolution(const MeshType& mesh, const RealType& time, SharedVector& output, 
+                const TimeFunction timeFunction, const AnalyticSpaceFunction analyticSpaceFunction);
+   
+   template< typename TimeFunction, typename AnalyticSpaceFunction>
+   void computeLaplace(const MeshType& mesh, const RealType& time, DofVectorType& output,
+                const TimeFunction timeFunction, const AnalyticSpaceFunction analyticSpaceFunction);
 };
 
 template<typename Real, typename Device, typename Index>
@@ -41,8 +45,12 @@ class AnalyticSolution<tnlGrid<2,Real,Device,Index,tnlIdenticalGridGeometry>>
    typedef tnlStaticVector< 2, RealType > VertexType;
    
    template< typename TimeFunction, typename AnalyticSpaceFunction>
-   void compute(const MeshType& mesh, const RealType& time, SharedVector& output, SharedVector& numerical,
-                TimeFunction timeFunction, AnalyticSpaceFunction anlayticSpaceFunction);
+   void computeAnalyticSolution(const MeshType& mesh, const RealType& time, SharedVector& output,
+                const TimeFunction timeFunction, const AnalyticSpaceFunction analyticSpaceFunction);
+
+   template< typename TimeFunction, typename AnalyticSpaceFunction>
+   void computeLaplace(const MeshType& mesh, const RealType& time, DofVectorType& output,
+                const TimeFunction timeFunction, const AnalyticSpaceFunction analyticSpaceFunction);
 };
 
 template<typename Real, typename Device, typename Index>
@@ -60,11 +68,14 @@ class AnalyticSolution<tnlGrid<3,Real,Device,Index,tnlIdenticalGridGeometry>>
    typedef tnlStaticVector< 3, RealType > VertexType;
    
    template< typename TimeFunction, typename AnalyticSpaceFunction>
-   void compute(const MeshType& mesh, const RealType& time, SharedVector& output, SharedVector& numerical,
-                TimeFunction timeFunction, AnalyticSpaceFunction anlayticSpaceFunction);
+   void computeAnalyticSolution(const MeshType& mesh, const RealType& time, SharedVector& output,
+                const TimeFunction timeFunction, const AnalyticSpaceFunction analyticSpaceFunction);
+
+   template< typename TimeFunction, typename AnalyticSpaceFunction>
+   void computeLaplace(const MeshType& mesh, const RealType& time, DofVectorType& output,
+                const TimeFunction timeFunction, const AnalyticSpaceFunction analyticSpaceFunction);
 };
 
 #include "tnlAnalyticSolution_impl.h"
 
 #endif	/* TNLANALYTICSOLUTION_H */
-
diff --git a/examples/heat-equation/tnlAnalyticSolution_impl.h b/examples/heat-equation/tnlAnalyticSolution_impl.h
index 252c007fc3..6d63ac9345 100644
--- a/examples/heat-equation/tnlAnalyticSolution_impl.h
+++ b/examples/heat-equation/tnlAnalyticSolution_impl.h
@@ -6,8 +6,8 @@
 template<typename Real, typename Device, typename Index>
 template< typename TimeFunction, typename AnalyticSpaceFunction>
 void AnalyticSolution<tnlGrid<1,Real,Device,Index,tnlIdenticalGridGeometry>>::
-compute(const MeshType& mesh, const RealType& time, SharedVector& output, 
-        SharedVector& numerical, TimeFunction timeFunction, AnalyticSpaceFunction analyticSpaceFunction)
+computeAnalyticSolution(const MeshType& mesh, const RealType& time, SharedVector& output, 
+        const TimeFunction timeFunction, const AnalyticSpaceFunction analyticSpaceFunction)
 {  
       RealType timeFunctionValue = timeFunction.getTimeValue(time);
       
@@ -23,11 +23,34 @@ compute(const MeshType& mesh, const RealType& time, SharedVector& output,
       }
 }
 
+template<typename Real, typename Device, typename Index>
+template< typename TimeFunction, typename AnalyticSpaceFunction>
+void AnalyticSolution<tnlGrid<1,Real,Device,Index,tnlIdenticalGridGeometry>>::
+computeLaplace(const MeshType& mesh, const RealType& time, DofVectorType& output,
+                const TimeFunction timeFunction, const AnalyticSpaceFunction analyticSpaceFunction)
+{        
+      RealType timeFunctionValue = timeFunction.getTimeValue(time);   
+   
+      VertexType vertex;
+      CoordinatesType coordinates;
+
+      #ifdef HAVE_OPENMP
+         #pragma omp parallel for private(coordinates,vertex)
+      #endif
+      for(IndexType i=1; i<(output.getSize()-1); i++)
+      {         
+         mesh.getElementCoordinates(i,coordinates);
+         mesh.getElementCenter(coordinates,vertex);
+
+      output[i] = timeFunctionValue*analyticSpaceFunction.template getF<2,0,0>(vertex);
+   }  
+}
+
 template<typename Real, typename Device, typename Index>
 template< typename TimeFunction, typename AnalyticSpaceFunction>
 void AnalyticSolution<tnlGrid<2,Real,Device,Index,tnlIdenticalGridGeometry>>::
-compute(const MeshType& mesh, const RealType& time, SharedVector& output, 
-        SharedVector& numerical, TimeFunction timeFunction, AnalyticSpaceFunction analyticSpaceFunction)
+computeAnalyticSolution(const MeshType& mesh, const RealType& time, SharedVector& output, 
+        const TimeFunction timeFunction, const AnalyticSpaceFunction analyticSpaceFunction)
 { 
    RealType timeFunctionValue = timeFunction.getTimeValue(time);
    
@@ -47,11 +70,42 @@ compute(const MeshType& mesh, const RealType& time, SharedVector& output,
    
 }
 
+template<typename Real, typename Device, typename Index>
+template< typename TimeFunction, typename AnalyticSpaceFunction>
+void AnalyticSolution<tnlGrid<2,Real,Device,Index,tnlIdenticalGridGeometry>>::
+computeLaplace(const MeshType& mesh, const RealType& time, DofVectorType& output,
+                const TimeFunction timeFunction, const AnalyticSpaceFunction analyticSpaceFunction)
+{
+   RealType timeFunctionValue = timeFunction.getTimeValue(time);
+   
+   VertexType vertex;
+   CoordinatesType coordinates;
+   CoordinatesType dimensions = mesh.getDimensions();
+
+   #ifdef HAVE_OPENMP
+    #pragma omp parallel for private(coordinates,vertex)
+   #endif
+   for(IndexType i=1; i<(dimensions.x()-1); i++)
+   {
+      for(IndexType j=1; j<(dimensions.y()-1); j++)
+      {  
+         coordinates.x()=i;
+         coordinates.y()=j;
+         
+         mesh.getElementCenter(coordinates,vertex);
+         
+         output[j*mesh.getDimensions().x()+i] = 
+                    timeFunctionValue*(analyticSpaceFunction.template getF<2,0,0>(vertex)+
+                    analyticSpaceFunction.template getF<0,2,0>(vertex));
+      } 
+   }
+}
+   
 template<typename Real, typename Device, typename Index>
 template< typename TimeFunction, typename AnalyticSpaceFunction>
 void AnalyticSolution<tnlGrid<3,Real,Device,Index,tnlIdenticalGridGeometry>>::
-compute(const MeshType& mesh, const RealType& time, SharedVector& output,
-        SharedVector& numerical, TimeFunction timeFunction, AnalyticSpaceFunction analyticSpaceFunction)
+computeAnalyticSolution(const MeshType& mesh, const RealType& time, SharedVector& output,
+        const TimeFunction timeFunction, const AnalyticSpaceFunction analyticSpaceFunction)
 {
    RealType timeFunctionValue = timeFunction.getTimeValue(time);
    
@@ -70,6 +124,39 @@ compute(const MeshType& mesh, const RealType& time, SharedVector& output,
    }   
 }
 
+template<typename Real, typename Device, typename Index>
+template< typename TimeFunction, typename AnalyticSpaceFunction>
+void AnalyticSolution<tnlGrid<3,Real,Device,Index,tnlIdenticalGridGeometry>>::
+computeLaplace(const MeshType& mesh, const RealType& time, DofVectorType& output,
+                const TimeFunction timeFunction, const AnalyticSpaceFunction analyticSpaceFunction)
+{ 
+   RealType timeFunctionValue = timeFunction.getTimeValue(time);
+      
+   VertexType vertex;
+   CoordinatesType coordinates;
+   CoordinatesType dimensions = mesh.getDimensions();
 
+   #ifdef HAVE_OPENMP
+    #pragma omp parallel for private(coordinates,vertex)
+   #endif
+   for(IndexType i=1; i<(dimensions.x()-1); i++)
+   {
+      for(IndexType j=1; j<(dimensions.y()-1); j++)
+      {  
+         for(IndexType k=1; k<(dimensions.y()-1); k++)
+         {
+            coordinates.x()=i;
+            coordinates.y()=j;
+            coordinates.z()=k;
+            
+            mesh.getElementCenter(coordinates,vertex);
+         
+            output[mesh.getElementIndex(i,j,k)] = 
+                       timeFunctionValue*(analyticSpaceFunction.template getF<2,0,0>(vertex)+
+                       analyticSpaceFunction.template getF<0,2,0>(vertex)+
+                       analyticSpaceFunction.template getF<0,0,2>(vertex));
+         }
+      } 
+   }
+}
 #endif	/* TNLANALYTICSOLUTION_IMPL_H */
-
diff --git a/examples/heat-equation/tnlLinearDiffusion.h b/examples/heat-equation/tnlLinearDiffusion.h
index 3605dc0045..f10c9ae1ee 100644
--- a/examples/heat-equation/tnlLinearDiffusion.h
+++ b/examples/heat-equation/tnlLinearDiffusion.h
@@ -21,16 +21,12 @@ class tnlLinearDiffusion<tnlGrid<1,Real,Device,Index,tnlIdenticalGridGeometry>>
    typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
 
    void getExplicitRHS( const MeshType& mesh,
-                        const RealType& time,
-                        const RealType& tau,
                         const CoordinatesType& coordinates,
                         DofVectorType& _u,
                         DofVectorType& _fu
                         );
 
    void getExplicitRHS( const MeshType& mesh,
-                        const RealType& time,
-                        const RealType& tau,
                         DofVectorType& _u,
                         DofVectorType& _fu);
    
@@ -49,16 +45,13 @@ class tnlLinearDiffusion<tnlGrid<2,Real,Device,Index,tnlIdenticalGridGeometry>>
    typedef typename MeshType::IndexType IndexType;
    typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
 
-   void getExplicitRHS(MeshType& mesh,
-                        const RealType& time,
-                        const RealType& tau,
+   void getExplicitRHS( const MeshType& mesh,
                         const CoordinatesType& coordinates,
                         DofVectorType& _u,
-                        DofVectorType& _fu);
- 
-   void getExplicitRHS(MeshType& mesh,
-                        const RealType& time,
-                        const RealType& tau,
+                        DofVectorType& _fu
+                        );
+
+   void getExplicitRHS( const MeshType& mesh,
                         DofVectorType& _u,
                         DofVectorType& _fu);
    
@@ -78,15 +71,12 @@ class tnlLinearDiffusion<tnlGrid<3,Real,Device,Index,tnlIdenticalGridGeometry>>
    typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
 
    void getExplicitRHS( const MeshType& mesh,
-                        const RealType& time,
-                        const RealType& tau,
                         const CoordinatesType& coordinates,
                         DofVectorType& _u,
-                        DofVectorType& _fu);
+                        DofVectorType& _fu
+                        );
 
    void getExplicitRHS( const MeshType& mesh,
-                        const RealType& time,
-                        const RealType& tau,
                         DofVectorType& _u,
                         DofVectorType& _fu);
    
@@ -97,4 +87,3 @@ class tnlLinearDiffusion<tnlGrid<3,Real,Device,Index,tnlIdenticalGridGeometry>>
 
 
 #endif	/* TNLLINEARDIFFUSION_H */
-
diff --git a/examples/heat-equation/tnlLinearDiffusion_impl.h b/examples/heat-equation/tnlLinearDiffusion_impl.h
index 128415e438..618930ed2f 100644
--- a/examples/heat-equation/tnlLinearDiffusion_impl.h
+++ b/examples/heat-equation/tnlLinearDiffusion_impl.h
@@ -6,8 +6,8 @@
 
 template<typename Real, typename Device, typename Index>
 void tnlLinearDiffusion<tnlGrid<1,Real,Device,Index,tnlIdenticalGridGeometry>> ::
-getExplicitRHS( const MeshType& mesh, const RealType& time, const RealType& tau,
-                const CoordinatesType& coordinates, DofVectorType& _u, DofVectorType& _fu)
+getExplicitRHS( const MeshType& mesh, const CoordinatesType& coordinates, 
+               DofVectorType& _u, DofVectorType& _fu)
 {
    if((coordinates.x() <= 0) || (coordinates.x() >= (mesh.getDimensions().x()-1)))
    {
@@ -23,8 +23,7 @@ getExplicitRHS( const MeshType& mesh, const RealType& time, const RealType& tau,
    
 template<typename Real, typename Device, typename Index>
 void tnlLinearDiffusion<tnlGrid<1,Real,Device,Index,tnlIdenticalGridGeometry>> ::
-getExplicitRHS( const MeshType& mesh, const RealType& time, const RealType& tau,
-                DofVectorType& _u, DofVectorType& _fu)
+getExplicitRHS( const MeshType& mesh, DofVectorType& _u, DofVectorType& _fu)
 {
    RealType stepXSquare = mesh.getParametricStep().x()*mesh.getParametricStep().x();
    
@@ -42,8 +41,8 @@ getExplicitRHS( const MeshType& mesh, const RealType& time, const RealType& tau,
 
 template<typename Real, typename Device, typename Index>
 void tnlLinearDiffusion<tnlGrid<2,Real,Device,Index,tnlIdenticalGridGeometry>> :: 
-getExplicitRHS(MeshType& mesh, const RealType& time, const RealType& tau,
-                const CoordinatesType& coordinates, DofVectorType& _u, DofVectorType& _fu)
+getExplicitRHS(const MeshType& mesh, const CoordinatesType& coordinates, 
+               DofVectorType& _u, DofVectorType& _fu)
 {
    if((coordinates.x() <= 0) || (coordinates.x() >= (mesh.getDimensions().x()-1)) || 
       (coordinates.y() <= 0 ) || (coordinates.y() >= (mesh.getDimensions().y()-1)))
@@ -63,8 +62,7 @@ getExplicitRHS(MeshType& mesh, const RealType& time, const RealType& tau,
    
 template<typename Real, typename Device, typename Index>
 void tnlLinearDiffusion<tnlGrid<2,Real,Device,Index,tnlIdenticalGridGeometry>> :: 
-getExplicitRHS(MeshType& mesh, const RealType& time, const RealType& tau,
-                DofVectorType& _u, DofVectorType& _fu)
+getExplicitRHS(const MeshType& mesh, DofVectorType& _u, DofVectorType& _fu)
 {
    RealType stepXSquare = mesh.getParametricStep().x()*mesh.getParametricStep().x();
    RealType stepYSquare = mesh.getParametricStep().y()*mesh.getParametricStep().y(); 
@@ -87,8 +85,8 @@ getExplicitRHS(MeshType& mesh, const RealType& time, const RealType& tau,
 
 template<typename Real, typename Device, typename Index>
 void tnlLinearDiffusion<tnlGrid<3,Real,Device,Index,tnlIdenticalGridGeometry>> :: 
-getExplicitRHS( const MeshType& mesh, const RealType& time, const RealType& tau,
-                const CoordinatesType& coordinates, DofVectorType& _u, DofVectorType& _fu)
+getExplicitRHS( const MeshType& mesh, const CoordinatesType& coordinates, 
+               DofVectorType& _u, DofVectorType& _fu)
 {
    if((coordinates.x() <= 0) || (coordinates.x() >= (mesh.getDimensions().x()-1)) || 
       (coordinates.y() <= 0 ) || (coordinates.y() >= (mesh.getDimensions().y()-1)) || 
@@ -118,8 +116,7 @@ getExplicitRHS( const MeshType& mesh, const RealType& time, const RealType& tau,
    
 template<typename Real, typename Device, typename Index>
 void tnlLinearDiffusion<tnlGrid<3,Real,Device,Index,tnlIdenticalGridGeometry>> :: 
-getExplicitRHS( const MeshType& mesh, const RealType& time, const RealType& tau,
-                DofVectorType& _u, DofVectorType& _fu)
+getExplicitRHS( const MeshType& mesh, DofVectorType& _u, DofVectorType& _fu)
 {
    RealType stepXSquare = mesh.getParametricStep().x()*mesh.getParametricStep().x();
    RealType stepYSquare = mesh.getParametricStep().y()*mesh.getParametricStep().y();   
@@ -150,4 +147,3 @@ getExplicitRHS( const MeshType& mesh, const RealType& time, const RealType& tau,
 
 
 #endif	/* TNLLINEARDIFFUSION_IMP_H */
-
diff --git a/examples/heat-equation/tnlTimeFunction.h b/examples/heat-equation/tnlTimeFunction.h
index ee8c5a90a3..060054a0d6 100644
--- a/examples/heat-equation/tnlTimeFunction.h
+++ b/examples/heat-equation/tnlTimeFunction.h
@@ -19,11 +19,17 @@ template<int dim, typename Real, typename Device, typename Index>
 class TimeFunction<tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>,TimeFunctionBase::TimeIndependent>: public TimeFunctionBase
 {
    public:
+   enum setTimeFunction{TimeIndependent};
    typedef tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry> MeshType;
-   typedef typename MeshType::RealType RealType; 
+   typedef typename MeshType::RealType RealType;
+   typedef typename MeshType::DeviceType DeviceType;
+   typedef typename MeshType::IndexType IndexType;
+   typedef tnlSharedVector< RealType, DeviceType, IndexType > SharedVectorType;
+   typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
    
-   RealType getTimeValue(const RealType& time);
-   RealType getDerivation(const RealType& time);
+   RealType getTimeValue(const RealType& time) const;
+   RealType getDerivation(const RealType& time) const;
+   void applyInitTimeValues(SharedVectorType& u) const;
 };
 
 template<>
@@ -31,11 +37,17 @@ template<int dim, typename Real, typename Device, typename Index>
 class TimeFunction<tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>,TimeFunctionBase::Linear>: public TimeFunctionBase
 {
    public:
+   enum setTimeFunction{Linear};
    typedef tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry> MeshType;
-   typedef typename MeshType::RealType RealType; 
+   typedef typename MeshType::RealType RealType;
+   typedef typename MeshType::DeviceType DeviceType;
+   typedef typename MeshType::IndexType IndexType;
+   typedef tnlSharedVector< RealType, DeviceType, IndexType > SharedVectorType;
+   typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType; 
    
-   RealType getTimeValue(const RealType& time);
-   RealType getDerivation(const RealType& time);
+   RealType getTimeValue(const RealType& time) const;
+   RealType getDerivation(const RealType& time) const;
+   void applyInitTimeValues(SharedVectorType& u) const;
 };
 
 template<>
@@ -43,11 +55,17 @@ template<int dim, typename Real, typename Device, typename Index>
 class TimeFunction<tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>,TimeFunctionBase::Quadratic>: public TimeFunctionBase
 {
    public:
+   enum setTimeFunction{Quadratic};
    typedef tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry> MeshType;
-   typedef typename MeshType::RealType RealType; 
+   typedef typename MeshType::RealType RealType;
+   typedef typename MeshType::DeviceType DeviceType;
+   typedef typename MeshType::IndexType IndexType;
+   typedef tnlSharedVector< RealType, DeviceType, IndexType > SharedVectorType;
+   typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
    
-   RealType getTimeValue(const RealType& time);
-   RealType getDerivation(const RealType& time);
+   RealType getTimeValue(const RealType& time) const;
+   RealType getDerivation(const RealType& time) const;
+   void applyInitTimeValues(SharedVectorType& u) const;
 };
 
 template<>
@@ -55,13 +73,19 @@ template<int dim, typename Real, typename Device, typename Index>
 class TimeFunction<tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>,TimeFunctionBase::Cosinus>: public TimeFunctionBase
 {
    public:
+   enum setTimeFunction{Cosinus};
    typedef tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry> MeshType;
-   typedef typename MeshType::RealType RealType; 
+   typedef typename MeshType::RealType RealType;
+   typedef typename MeshType::DeviceType DeviceType;
+   typedef typename MeshType::IndexType IndexType;
+   typedef tnlSharedVector< RealType, DeviceType, IndexType > SharedVectorType;
+   typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
    
-   RealType getTimeValue(const RealType& time);
-   RealType getDerivation(const RealType& time);
+   RealType getTimeValue(const RealType& time) const;
+   RealType getDerivation(const RealType& time) const;
+   void applyInitTimeValues(SharedVectorType& u) const;
 };
 
 #include "tnlTimeFunction_impl.h"
 
-#endif	/* TNLTIMEFUNCTION_H */
+#endif	/* TNLTIMEFUNCTION_H */
\ No newline at end of file
diff --git a/examples/heat-equation/tnlTimeFunction_impl.h b/examples/heat-equation/tnlTimeFunction_impl.h
index fb4115254b..570bd799c5 100644
--- a/examples/heat-equation/tnlTimeFunction_impl.h
+++ b/examples/heat-equation/tnlTimeFunction_impl.h
@@ -7,7 +7,7 @@ template<>
 template<int dim, typename Real, typename Device, typename Index>
 typename tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>::RealType 
 TimeFunction<tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>,TimeFunctionBase::TimeIndependent>::
-getTimeValue(const RealType& time)
+getTimeValue(const RealType& time) const
 {
       return 1.0;
 }
@@ -16,7 +16,7 @@ template<>
 template<int dim, typename Real, typename Device, typename Index>
 typename tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>::RealType 
 TimeFunction<tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>,TimeFunctionBase::Linear>::
-getTimeValue(const RealType& time)
+getTimeValue(const RealType& time) const
 {
       return time;
 }
@@ -25,7 +25,7 @@ template<>
 template<int dim, typename Real, typename Device, typename Index>
 typename tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>::RealType 
 TimeFunction<tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>,TimeFunctionBase::Quadratic>::
-getTimeValue(const RealType& time)
+getTimeValue(const RealType& time) const
 {
       return time*time;
 }
@@ -34,7 +34,7 @@ template<>
 template<int dim, typename Real, typename Device, typename Index>
 typename tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>::RealType 
 TimeFunction<tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>,TimeFunctionBase::Cosinus>::
-getTimeValue(const RealType& time)
+getTimeValue(const RealType& time) const
 {
       return cos(time);
 }
@@ -43,7 +43,7 @@ template<>
 template<int dim, typename Real, typename Device, typename Index>
 typename tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>::RealType 
 TimeFunction<tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>,TimeFunctionBase::TimeIndependent>::
-getDerivation(const RealType& time)
+getDerivation(const RealType& time) const
 {
       return 0.0;
 }
@@ -52,7 +52,7 @@ template<>
 template<int dim, typename Real, typename Device, typename Index>
 typename tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>::RealType 
 TimeFunction<tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>,TimeFunctionBase::Linear>::
-getDerivation(const RealType& time)
+getDerivation(const RealType& time) const
 {
       return 1.0;
 }
@@ -61,7 +61,7 @@ template<>
 template<int dim, typename Real, typename Device, typename Index>
 typename tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>::RealType 
 TimeFunction<tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>,TimeFunctionBase::Quadratic>::
-getDerivation(const RealType& time)
+getDerivation(const RealType& time) const
 {
       return 2.0*time;
 }
@@ -70,10 +70,53 @@ template<>
 template<int dim, typename Real, typename Device, typename Index>
 typename tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>::RealType 
 TimeFunction<tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>,TimeFunctionBase::Cosinus>::
-getDerivation(const RealType& time)
+getDerivation(const RealType& time) const
 {
       return -sin(time);
 }
 
-#endif	/* TNLTIMEFUNCTION_IMPL_H */
+template<>
+template<int dim, typename Real, typename Device, typename Index>
+void TimeFunction<tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>,TimeFunctionBase::TimeIndependent>::
+applyInitTimeValues(SharedVectorType& u) const
+{
+   
+}
+
+template<>
+template<int dim, typename Real, typename Device, typename Index>
+void TimeFunction<tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>,TimeFunctionBase::Linear>::
+applyInitTimeValues(SharedVectorType& u) const
+{
+   #ifdef HAVE_OPENMP
+      #pragma omp parallel for
+   #endif
+   for(IndexType i=1; i<(u.getSize()-1); i++)
+   {         
+     u[i] = 0;
+   }
+}
+
+template<>
+template<int dim, typename Real, typename Device, typename Index>
+void TimeFunction<tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>,TimeFunctionBase::Quadratic>::
+applyInitTimeValues(SharedVectorType& u) const
+{
+   #ifdef HAVE_OPENMP
+      #pragma omp parallel for
+   #endif
+   for(IndexType i=1; i<(u.getSize()-1); i++)
+   {         
+     u[i] = 0;
+   }
+}
 
+template<>
+template<int dim, typename Real, typename Device, typename Index>
+void TimeFunction<tnlGrid<dim,Real,Device,Index,tnlIdenticalGridGeometry>,TimeFunctionBase::Cosinus>::
+applyInitTimeValues(SharedVectorType& u) const
+{
+
+}
+
+#endif	/* TNLTIMEFUNCTION_IMPL_H */
-- 
GitLab