diff --git a/src/TNL/Experimental/Hamilton-Jacobi/CMakeLists.txt b/src/TNL/Experimental/Hamilton-Jacobi/CMakeLists.txt
index 538e5d99ccaacd6892be2a378da2f19349ec3fa4..3c15fe33154a4d5ce4dd2dff3ae524b768e21797 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/CMakeLists.txt
+++ b/src/TNL/Experimental/Hamilton-Jacobi/CMakeLists.txt
@@ -1,3 +1,3 @@
 ADD_SUBDIRECTORY( Operators )
-#ADD_SUBDIRECTORY( Solvers )
+ADD_SUBDIRECTORY( Solvers )
 
diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/CMakeLists.txt b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/CMakeLists.txt
index 04348b8b5baa4bdbd772b67feaf9da3b32aed53d..176c9859e3c47a85a1e79356daff54251b34ec02 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/CMakeLists.txt
+++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/CMakeLists.txt
@@ -7,14 +7,14 @@ set( tnl_hamilton_jacobi_SOURCES
      hamiltonJacobiProblemConfig.h
      tnl-direct-eikonal-solver.h )
                
-ADD_EXECUTABLE(tnl-hamilton-jacobi${debugExt} main.cpp)
-target_link_libraries (tnl-hamilton-jacobi${debugExt} tnl${debugExt}-${tnlVersion} )
+#ADD_EXECUTABLE(tnl-hamilton-jacobi${debugExt} main.cpp)
+#target_link_libraries (tnl-hamilton-jacobi${debugExt} tnl${debugExt}-${tnlVersion} )
 
 ADD_EXECUTABLE(tnl-direct-eikonal-solver${debugExt} tnl-direct-eikonal-solver.cpp )
 target_link_libraries (tnl-direct-eikonal-solver${debugExt} tnl${debugExt}-${tnlVersion} )
 
 
-INSTALL( TARGETS tnl-hamilton-jacobi${debugExt} 
+INSTALL( TARGETS #tnl-hamilton-jacobi${debugExt} 
                  tnl-direct-eikonal-solver${debugExt}
          RUNTIME DESTINATION bin
          PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE )
diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblemConfig.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblemConfig.h
index 4e84c7d98966b3e35d86e263850a7df14d85ff13..52a441165aba70c6641edb9590b0a6164375f2f1 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblemConfig.h
+++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblemConfig.h
@@ -17,7 +17,7 @@
 
 #pragma once 
 
-#include <config/tnlConfigDescription.h>
+#include <TNL/Config/ConfigDescription.h>
 
 template< typename ConfigTag >
 class HamiltonJacobiProblemConfig
diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnl-direct-eikonal-solver.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnl-direct-eikonal-solver.h
index b453d0e6a909c9145e07d674267167e9e8dda213..225ce04b9851b93dd00d37644a7c042f563a3b9d 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnl-direct-eikonal-solver.h
+++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnl-direct-eikonal-solver.h
@@ -13,23 +13,24 @@
 
 #pragma once
 
-#include <solvers/tnlSolver.h>
-#include <solvers/tnlFastBuildConfigTag.h>
-#include <solvers/tnlBuildConfigTags.h>
-#include <functions/tnlConstantFunction.h>
-#include <functions/tnlMeshFunction.h>
-//#include <problems/tnlHeatEquationProblem.h>
-#include <mesh/tnlGrid.h>
+#include <TNL/Solvers/Solver.h>
+#include <TNL/Solvers/FastBuildConfigTag.h>
+#include <TNL/Solvers/BuildConfigTags.h>
+#include <TNL/Functions/Analytic/Constant.h>
+#include <TNL/Functions/MeshFunction.h>
+#include <TNL/Meshes/Grid.h>
 #include "tnlDirectEikonalProblem.h"
 
+using namespace TNL;
+
 //typedef tnlDefaultBuildMeshConfig BuildConfig;
-typedef tnlFastBuildConfig BuildConfig;
+typedef Solvers::FastBuildConfig BuildConfig;
 
 template< typename MeshConfig >
 class tnlDirectEikonalSolverConfig
 {
    public:
-      static void configSetup( tnlConfigDescription& config )
+      static void configSetup( Config::ConfigDescription& config )
       {
          config.addDelimiter( "Direct eikonal equation solver settings:" );
          config.addRequiredEntry< String >( "input-file", "Input file." );
@@ -50,12 +51,12 @@ class tnlDirectEikonalSolverSetter
    typedef Device DeviceType;
    typedef Index IndexType;
 
-   typedef Containers::StaticVector< MeshType::meshDimensions, Real > Point;
+   typedef Containers::StaticVector< MeshType::getMeshDimension(), Real > Point;
 
    static bool run( const Config::ParameterContainer& parameters )
    {
-      enum { Dimensions = MeshType::meshDimensions };
-      typedef tnlConstantFunction< Dimensions, Real > Anisotropy;
+      static const int Dimension = MeshType::getMeshDimension();
+      typedef Functions::Analytic::Constant< Dimension, Real > Anisotropy;
       typedef tnlDirectEikonalProblem< MeshType, Anisotropy > Problem;
       SolverStarter solverStarter;
       return solverStarter.template run< Problem >( parameters );
@@ -64,7 +65,7 @@ class tnlDirectEikonalSolverSetter
 
 int main( int argc, char* argv[] )
 {
-   if( ! tnlSolver< tnlDirectEikonalSolverSetter, tnlDirectEikonalSolverConfig, BuildConfig >::run( argc, argv ) )
+   if( ! Solvers::Solver< tnlDirectEikonalSolverSetter, tnlDirectEikonalSolverConfig, BuildConfig >::run( argc, argv ) )
       return EXIT_FAILURE;
    return EXIT_SUCCESS;
 }
diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h
index 409a2d93a1d7eb719e3b4f93b6eb9b7da53f72e5..c6385f64a9668272bee961d0bb5dc02b182655cd 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h
+++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h
@@ -7,8 +7,10 @@
 
 #pragma once 
 
-#include <mesh/tnlGrid.h>
-#include <functions/tnlMeshFunction.h>
+#include <TNL/Meshes/Grid.h>
+#include <TNL/Functions/MeshFunction.h>
+
+using namespace TNL;
 
 template< typename Mesh >
 class tnlDirectEikonalMethodsBase
@@ -18,16 +20,16 @@ class tnlDirectEikonalMethodsBase
 template< typename Real,
           typename Device,
           typename Index >
-class tnlDirectEikonalMethodsBase< tnlGrid< 1, Real, Device, Index > >
+class tnlDirectEikonalMethodsBase< Meshes::Grid< 1, Real, Device, Index > >
 {
    public:
       
-      typedef tnlGrid< 1, Real, Device, Index > MeshType;
+      typedef Meshes::Grid< 1, Real, Device, Index > MeshType;
       typedef Real RealType;
       typedef Device DevcieType;
       typedef Index IndexType;
-      typedef tnlMeshFunction< MeshType > MeshFunctionType;
-      typedef tnlMeshFunction< MeshType, 1, bool > InterfaceMapType;
+      typedef Functions::MeshFunction< MeshType > MeshFunctionType;
+      typedef Functions::MeshFunction< MeshType, 1, bool > InterfaceMapType;
       
       void initInterface( const MeshFunctionType& input,
                           MeshFunctionType& output,
@@ -43,16 +45,16 @@ class tnlDirectEikonalMethodsBase< tnlGrid< 1, Real, Device, Index > >
 template< typename Real,
           typename Device,
           typename Index >
-class tnlDirectEikonalMethodsBase< tnlGrid< 2, Real, Device, Index > >
+class tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > >
 {
    public:
       
-      typedef tnlGrid< 2, Real, Device, Index > MeshType;
+      typedef Meshes::Grid< 2, Real, Device, Index > MeshType;
       typedef Real RealType;
       typedef Device DevcieType;
       typedef Index IndexType;
-      typedef tnlMeshFunction< MeshType > MeshFunctionType;
-      typedef tnlMeshFunction< MeshType, 2, bool > InterfaceMapType;
+      typedef Functions::MeshFunction< MeshType > MeshFunctionType;
+      typedef Functions::MeshFunction< MeshType, 2, bool > InterfaceMapType;
 
       void initInterface( const MeshFunctionType& input,
                           MeshFunctionType& output,
@@ -66,16 +68,16 @@ class tnlDirectEikonalMethodsBase< tnlGrid< 2, Real, Device, Index > >
 template< typename Real,
           typename Device,
           typename Index >
-class tnlDirectEikonalMethodsBase< tnlGrid< 3, Real, Device, Index > >
+class tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > >
 {
    public:
       
-      typedef tnlGrid< 3, Real, Device, Index > MeshType;
+      typedef Meshes::Grid< 3, Real, Device, Index > MeshType;
       typedef Real RealType;
       typedef Device DevcieType;
       typedef Index IndexType;
-      typedef tnlMeshFunction< MeshType > MeshFunctionType;
-      typedef tnlMeshFunction< MeshType, 3, bool > InterfaceMapType;
+      typedef Functions::MeshFunction< MeshType > MeshFunctionType;
+      typedef Functions::MeshFunction< MeshType, 3, bool > InterfaceMapType;
 
       void initInterface( const MeshFunctionType& input,
                           MeshFunctionType& output,
diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h
index 88ed4ddf4f6cbf10d45c89ed1555b564b73de3ac..904253116e7e45a16f1cd10988bba28f4dbc5994 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h
+++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h
@@ -7,14 +7,13 @@
 
 #pragma once
 
-#include <core/tnlTypeInfo.h>
-#include <functions/tnlFunctions.h>
+#include <TNL/TypeInfo.h>
 
 template< typename Real,
           typename Device,
           typename Index >
 void
-tnlDirectEikonalMethodsBase< tnlGrid< 1, Real, Device, Index > >::
+tnlDirectEikonalMethodsBase< Meshes::Grid< 1, Real, Device, Index > >::
 initInterface( const MeshFunctionType& input,
                MeshFunctionType& output,
                InterfaceMapType& interfaceMap  )
@@ -66,8 +65,8 @@ initInterface( const MeshFunctionType& input,
          }
       }
       output[ cell.getIndex() ] =
-      c > 0 ? tnlTypeInfo< RealType >::getMaxValue() :
-             -tnlTypeInfo< RealType >::getMaxValue();
+      c > 0 ? TypeInfo< RealType >::getMaxValue() :
+             -TypeInfo< RealType >::getMaxValue();
       interfaceMap[ cell.getIndex() ] = false;
    }
 }
@@ -77,7 +76,7 @@ template< typename Real,
           typename Index >
    template< typename MeshEntity >
 void
-tnlDirectEikonalMethodsBase< tnlGrid< 1, Real, Device, Index > >::
+tnlDirectEikonalMethodsBase< Meshes::Grid< 1, Real, Device, Index > >::
 updateCell( MeshFunctionType& u,
             const MeshEntity& cell )
 {
@@ -88,7 +87,7 @@ template< typename Real,
           typename Device,
           typename Index >
 void
-tnlDirectEikonalMethodsBase< tnlGrid< 2, Real, Device, Index > >::
+tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > >::
 initInterface( const MeshFunctionType& input,
                MeshFunctionType& output,
                InterfaceMapType& interfaceMap  )
@@ -124,8 +123,8 @@ initInterface( const MeshFunctionType& input,
             }
          }
          output[ cell.getIndex() ] =
-            c > 0 ? tnlTypeInfo< RealType >::getMaxValue() :
-                   -tnlTypeInfo< RealType >::getMaxValue();  
+            c > 0 ? TypeInfo< RealType >::getMaxValue() :
+                   -TypeInfo< RealType >::getMaxValue();  
          interfaceMap[ cell.getIndex() ] = false;
       }
 }
@@ -135,7 +134,7 @@ template< typename Real,
           typename Index >
    template< typename MeshEntity >
 void
-tnlDirectEikonalMethodsBase< tnlGrid< 2, Real, Device, Index > >::
+tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > >::
 updateCell( MeshFunctionType& u,
             const MeshEntity& cell )
 {
@@ -152,7 +151,7 @@ updateCell( MeshFunctionType& u,
       a = u[ neighborEntities.template getEntityIndex< -1,  0 >() ];
    else
    {
-      a = ArgAbsMin( u[ neighborEntities.template getEntityIndex< -1,  0 >() ],
+      a = argAbsMin( u[ neighborEntities.template getEntityIndex< -1,  0 >() ],
                      u[ neighborEntities.template getEntityIndex<  1,  0 >() ] );
    }
 
@@ -162,18 +161,18 @@ updateCell( MeshFunctionType& u,
       b = u[ neighborEntities.template getEntityIndex< 0,  -1 >() ];
    else
    {
-      b = ArgAbsMin( u[ neighborEntities.template getEntityIndex< 0,  -1 >() ],
+      b = argAbsMin( u[ neighborEntities.template getEntityIndex< 0,  -1 >() ],
                      u[ neighborEntities.template getEntityIndex< 0,   1 >() ] );
    }
 
-   if( fabs( a ) == tnlTypeInfo< Real >::getMaxValue() && 
-       fabs( b ) == tnlTypeInfo< Real >::getMaxValue() )
+   if( fabs( a ) == TypeInfo< Real >::getMaxValue() && 
+       fabs( b ) == TypeInfo< Real >::getMaxValue() )
       return;
-   if( fabs( a ) == tnlTypeInfo< Real >::getMaxValue() ||
-       fabs( b ) == tnlTypeInfo< Real >::getMaxValue() ||
+   if( fabs( a ) == TypeInfo< Real >::getMaxValue() ||
+       fabs( b ) == TypeInfo< Real >::getMaxValue() ||
        fabs( a - b ) >= h )
    {
-      tmp = ArgAbsMin( a, b ) + sign( value ) * h;
+      tmp = argAbsMin( a, b ) + sign( value ) * h;
       /*   std::cerr << "a = " << a << " b = " << b << " h = " << h 
              << " ArgAbsMin( a, b ) = " << ArgAbsMin( a, b ) << " sign( value ) = " << sign( value )
              << " sign( value ) * h = " << sign( value ) * h
@@ -190,7 +189,7 @@ updateCell( MeshFunctionType& u,
    else
       tmp = 0.5 * ( a + b + sign( value ) * sqrt( 2.0 * h * h - ( a - b ) * ( a - b ) ) );
 
-   u[ cell.getIndex() ] = ArgAbsMin( value, tmp );
+   u[ cell.getIndex() ] = argAbsMin( value, tmp );
    //std::cerr << ArgAbsMin( value, tmp ) << " ";   
 }
 
@@ -199,7 +198,7 @@ template< typename Real,
           typename Device,
           typename Index >
 void
-tnlDirectEikonalMethodsBase< tnlGrid< 3, Real, Device, Index > >::
+tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > >::
 initInterface( const MeshFunctionType& input,
                MeshFunctionType& output,
                InterfaceMapType& interfaceMap  )
@@ -240,8 +239,8 @@ initInterface( const MeshFunctionType& input,
                }
             }
             output[ cell.getIndex() ] =
-               c > 0 ? tnlTypeInfo< RealType >::getMaxValue() :
-                      -tnlTypeInfo< RealType >::getMaxValue();
+               c > 0 ? TypeInfo< RealType >::getMaxValue() :
+                      -TypeInfo< RealType >::getMaxValue();
             interfaceMap[ cell.getIndex() ] = false;
          }
 }
@@ -251,7 +250,7 @@ template< typename Real,
           typename Index >
    template< typename MeshEntity >
 void
-tnlDirectEikonalMethodsBase< tnlGrid< 3, Real, Device, Index > >::
+tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > >::
 updateCell( MeshFunctionType& u,
             const MeshEntity& cell )
 {
diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem.h
index f4baba1528165a9ca95bf18b68eaebe09d0c7848..ca3b8368ccf7a78650ba5fc985f89a52877a9a70 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem.h
+++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem.h
@@ -13,8 +13,9 @@
 
 #pragma once
 
-#include <problems/tnlPDEProblem.h>
-#include <functions/tnlMeshFunction.h>
+#include <TNL/Problems/PDEProblem.h>
+#include <TNL/Functions/MeshFunction.h>
+#include <TNL/SharedPointer.h>
 #include "tnlFastSweepingMethod.h"
 
 template< typename Mesh,
@@ -22,49 +23,55 @@ template< typename Mesh,
           typename Real = typename Mesh::RealType,
           typename Index = typename Mesh::IndexType >
 class tnlDirectEikonalProblem
-   : public tnlPDEProblem< Mesh,
-                           TimeIndependentProblem,
-                           Real,
-                           typename Mesh::DeviceType,
-                           Index  >
+   : public Problems::PDEProblem< Mesh,
+                                  Real,
+                                  typename Mesh::DeviceType,
+                                  Index  >
 {
    public:
    
       typedef Real RealType;
       typedef typename Mesh::DeviceType DeviceType;
       typedef Index IndexType;
-      typedef tnlMeshFunction< Mesh > MeshFunctionType;
-      typedef tnlPDEProblem< Mesh, TimeIndependentProblem, RealType, DeviceType, IndexType > BaseType;
-      typedef Anisotropy AnisotropyType;
+      typedef Functions::MeshFunction< Mesh > MeshFunctionType;
+      typedef Problems::PDEProblem< Mesh, RealType, DeviceType, IndexType > BaseType;
+      using AnisotropyType = Anisotropy;
 
       using typename BaseType::MeshType;
       using typename BaseType::DofVectorType;
       using typename BaseType::MeshDependentDataType;
+      using MeshPointer = SharedPointer< MeshType >;
+      using DofVectorPointer = SharedPointer< DofVectorType >;
+      using MeshDependentDataPointer = SharedPointer< MeshDependentDataType >;
+      
+      static constexpr bool isTimeDependent() { return false; };
 
       static String getTypeStatic();
 
       String getPrologHeader() const;
 
-      void writeProlog( tnlLogger& logger,
+      void writeProlog( Logger& logger,
                         const Config::ParameterContainer& parameters ) const;
       
-      bool writeEpilog( tnlLogger& logger );
+      bool writeEpilog( Logger& logger );
 
 
-      bool setup( const Config::ParameterContainer& parameters );
+      bool setup( const MeshPointer& mesh,
+                  const Config::ParameterContainer& parameters,
+                  const String& prefix );
 
-      IndexType getDofs( const MeshType& mesh ) const;
+      IndexType getDofs( const MeshPointer& mesh ) const;
 
-      void bindDofs( const MeshType& mesh,
-                     const DofVectorType& dofs );
+      void bindDofs( const MeshPointer& mesh,
+                     const DofVectorPointer& dofs );
       
-      bool setInitialData( const Config::ParameterContainer& parameters,
-                           const MeshType& mesh,
-                           DofVectorType& dofs,
-                           MeshDependentDataType& meshdependentData );
+      bool setInitialCondition( const Config::ParameterContainer& parameters,
+                                const MeshPointer& mesh,
+                                DofVectorPointer& dofs,
+                                MeshDependentDataPointer& meshdependentData );
 
-      bool solve( const MeshType& mesh,
-                  DofVectorType& dosf );
+      bool solve( const MeshPointer& mesh,
+                  DofVectorPointer& dosf );
 
 
       protected:
diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem_impl.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem_impl.h
index 9dcf2fb5fa40b49462bed8ced05b4ebf605c38d6..b10c6949e6ffea16eec6e7a4129186c7d0d817b5 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem_impl.h
+++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem_impl.h
@@ -21,7 +21,11 @@ String
 tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::
 getTypeStatic()
 {
-   
+   return String( "DirectEikonalProblem< " + 
+                  Mesh::getTypeStatic() + ", " +
+                  Anisotropy::getTypeStatic() + ", " +
+                  Real::getTypeStatic() + ", " +
+                  Index::getTypeStatic() + " >" );
 }
 
 template< typename Mesh,
@@ -41,7 +45,7 @@ template< typename Mesh,
           typename Index >
 void
 tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::
-writeProlog( tnlLogger& logger,
+writeProlog( Logger& logger,
              const Config::ParameterContainer& parameters ) const
 {
    
@@ -53,7 +57,7 @@ template< typename Mesh,
           typename Index >
 bool
 tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::
-writeEpilog( tnlLogger& logger )
+writeEpilog( Logger& logger )
 {
    return true;
 }
@@ -64,7 +68,9 @@ template< typename Mesh,
           typename Index >
 bool
 tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::
-setup( const Config::ParameterContainer& parameters )
+setup( const MeshPointer& mesh,
+       const Config::ParameterContainer& parameters,
+       const String& prefix )
 {
    return true;
 }
@@ -75,9 +81,9 @@ template< typename Mesh,
           typename Index >
 Index
 tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::
-getDofs( const MeshType& mesh ) const
+getDofs( const MeshPointer& mesh ) const
 {
-   return mesh.template getEntitiesCount< typename MeshType::Cell >();
+   return mesh->template getEntitiesCount< typename MeshType::Cell >();
 }
 
 template< typename Mesh,
@@ -86,8 +92,8 @@ template< typename Mesh,
           typename Index >
 void
 tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::
-bindDofs( const MeshType& mesh,
-          const DofVectorType& dofs )
+bindDofs( const MeshPointer& mesh,
+          const DofVectorPointer& dofs )
 {
    this->u.bind( mesh, dofs );
 }
@@ -98,10 +104,10 @@ template< typename Mesh,
           typename Index >
 bool
 tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::
-setInitialData( const Config::ParameterContainer& parameters,
-                const MeshType& mesh,
-                DofVectorType& dofs,
-                MeshDependentDataType& meshdependentData )
+setInitialCondition( const Config::ParameterContainer& parameters,
+                     const MeshPointer& mesh,
+                     DofVectorPointer& dofs,
+                     MeshDependentDataPointer& meshdependentData )
 {
    String inputFile = parameters.getParameter< String >( "input-file" );
    this->initialData.setMesh( mesh );
@@ -117,9 +123,10 @@ template< typename Mesh,
           typename Index >
 bool
 tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::
-solve( const MeshType& mesh,
-       DofVectorType& dofs )
+solve( const MeshPointer& mesh,
+       DofVectorPointer& dofs )
 {
-   tnlFastSweepingMethod< MeshType, AnisotropyType > fsm;
+   FastSweepingMethod< MeshType, AnisotropyType > fsm;
    fsm.solve( mesh, anisotropy, initialData );
+   return true;
 }
diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod.h
index 59ddc2a690b5e6aecb62d5cf4c0997b7c567ff82..17d139b0ec8b01c0dbf87188968be67ee5fa5cca 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod.h
+++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod.h
@@ -1,25 +1,23 @@
-/*
- * To change this license header, choose License Headers in Project Properties.
- * To change this template file, choose Tools | Templates
- * and open the template in the editor.
- */
+/***************************************************************************
+                          FastSweepingMethod.h  -  description
+                             -------------------
+    begin                : Jul 14, 2016
+    copyright            : (C) 2017 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
 
-/* 
- * File:   tnlFastSweepingMethod.h
- * Author: oberhuber
- *
- * Created on July 14, 2016, 10:04 AM
- */
+/* See Copyright Notice in tnl/Copyright */
 
 #pragma once
 
-#include <mesh/tnlGrid.h>
-#include <functions/tnlConstantFunction.h>
+#include <TNL/Meshes/Grid.h>
+#include <TNL/Functions/Analytic/Constant.h>
+#include <TNL/SharedPointer.h>
 #include "tnlDirectEikonalMethodsBase.h"
 
 template< typename Mesh,
-          typename Anisotropy = tnlConstantFunction< Mesh::getMeshDimension(), typename Mesh::RealType > >
-class tnlFastSweepingMethod
+          typename Anisotropy = Functions::Analytic::Constant< Mesh::getMeshDimension(), typename Mesh::RealType > >
+class FastSweepingMethod
 {   
 };
 
@@ -27,30 +25,31 @@ template< typename Real,
           typename Device,
           typename Index,
           typename Anisotropy >
-class tnlFastSweepingMethod< tnlGrid< 1, Real, Device, Index >, Anisotropy >
-   : public tnlDirectEikonalMethodsBase< tnlGrid< 1, Real, Device, Index > >
+class FastSweepingMethod< Meshes::Grid< 1, Real, Device, Index >, Anisotropy >
+   : public tnlDirectEikonalMethodsBase< Meshes::Grid< 1, Real, Device, Index > >
 {
    static_assert(  std::is_same< Device, TNL::Devices::Host >::value, "The fast sweeping method works only on CPU." );
    
    public:
       
-      typedef tnlGrid< 1, Real, Device, Index > MeshType;
+      typedef Meshes::Grid< 1, Real, Device, Index > MeshType;
       typedef Real RealType;
       typedef TNL::Devices::Host DeviceType;
       typedef Index IndexType;
       typedef Anisotropy AnisotropyType;
-      typedef tnlDirectEikonalMethodsBase< tnlGrid< 1, Real, Device, Index > > BaseType;
+      typedef tnlDirectEikonalMethodsBase< Meshes::Grid< 1, Real, Device, Index > > BaseType;
+      using MeshPointer = SharedPointer< MeshType >;
       
       using typename BaseType::InterfaceMapType;
       using typename BaseType::MeshFunctionType;
       
-      tnlFastSweepingMethod();
+      FastSweepingMethod();
       
       const IndexType& getMaxIterations() const;
       
       void setMaxIterations( const IndexType& maxIterations );
       
-      void solve( const MeshType& mesh,
+      void solve( const MeshPointer& mesh,
                   const AnisotropyType& anisotropy,
                   MeshFunctionType& u );
       
@@ -64,30 +63,31 @@ template< typename Real,
           typename Device,
           typename Index,
           typename Anisotropy >
-class tnlFastSweepingMethod< tnlGrid< 2, Real, Device, Index >, Anisotropy >
-   : public tnlDirectEikonalMethodsBase< tnlGrid< 2, Real, Device, Index > >
+class FastSweepingMethod< Meshes::Grid< 2, Real, Device, Index >, Anisotropy >
+   : public tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > >
 {
    static_assert(  std::is_same< Device, TNL::Devices::Host >::value, "The fast sweeping method works only on CPU." );
    
    public:
       
-      typedef tnlGrid< 2, Real, Device, Index > MeshType;
+      typedef Meshes::Grid< 2, Real, Device, Index > MeshType;
       typedef Real RealType;
       typedef TNL::Devices::Host DeviceType;
       typedef Index IndexType;
       typedef Anisotropy AnisotropyType;
-      typedef tnlDirectEikonalMethodsBase< tnlGrid< 2, Real, Device, Index > > BaseType;
+      typedef tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > > BaseType;
+      using MeshPointer = SharedPointer< MeshType >;
 
       using typename BaseType::InterfaceMapType;
       using typename BaseType::MeshFunctionType;
 
-      tnlFastSweepingMethod();
+      FastSweepingMethod();
       
       const IndexType& getMaxIterations() const;
       
       void setMaxIterations( const IndexType& maxIterations );
       
-      void solve( const MeshType& mesh,
+      void solve( const MeshPointer& mesh,
                   const AnisotropyType& anisotropy,
                   MeshFunctionType& u );
       
@@ -101,30 +101,31 @@ template< typename Real,
           typename Device,
           typename Index,
           typename Anisotropy >
-class tnlFastSweepingMethod< tnlGrid< 3, Real, Device, Index >, Anisotropy >
-   : public tnlDirectEikonalMethodsBase< tnlGrid< 3, Real, Device, Index > >
+class FastSweepingMethod< Meshes::Grid< 3, Real, Device, Index >, Anisotropy >
+   : public tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > >
 {
    static_assert(  std::is_same< Device, TNL::Devices::Host >::value, "The fast sweeping method works only on CPU." );
    
    public:
       
-      typedef tnlGrid< 3, Real, Device, Index > MeshType;
+      typedef Meshes::Grid< 3, Real, Device, Index > MeshType;
       typedef Real RealType;
       typedef TNL::Devices::Host DeviceType;
       typedef Index IndexType;
       typedef Anisotropy AnisotropyType;
-      typedef tnlDirectEikonalMethodsBase< tnlGrid< 3, Real, Device, Index > > BaseType;
+      typedef tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > > BaseType;
+      using MeshPointer = SharedPointer< MeshType >;
       
       using typename BaseType::InterfaceMapType;
       using typename BaseType::MeshFunctionType;
       
-      tnlFastSweepingMethod();
+      FastSweepingMethod();
       
       const IndexType& getMaxIterations() const;
       
       void setMaxIterations( const IndexType& maxIterations );
       
-      void solve( const MeshType& mesh,
+      void solve( const MeshPointer& mesh,
                   const AnisotropyType& anisotropy,
                   MeshFunctionType& u );
       
diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod1D_impl.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod1D_impl.h
index 050bc6075454c3d23685f25a6e9d607b36c88d73..874945eafd7f0cb1eff4667eba6ef026ba026634 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod1D_impl.h
+++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod1D_impl.h
@@ -19,8 +19,8 @@ template< typename Real,
           typename Device,
           typename Index,
           typename Anisotropy >
-tnlFastSweepingMethod< tnlGrid< 1, Real, Device, Index >, Anisotropy >::
-tnlFastSweepingMethod()
+FastSweepingMethod< Meshes::Grid< 1, Real, Device, Index >, Anisotropy >::
+FastSweepingMethod()
 : maxIterations( 1 )
 {
    
@@ -31,7 +31,7 @@ template< typename Real,
           typename Index,
           typename Anisotropy >
 const Index&
-tnlFastSweepingMethod< tnlGrid< 1, Real, Device, Index >, Anisotropy >::
+FastSweepingMethod< Meshes::Grid< 1, Real, Device, Index >, Anisotropy >::
 getMaxIterations() const
 {
    
@@ -42,7 +42,7 @@ template< typename Real,
           typename Index,
           typename Anisotropy >
 void
-tnlFastSweepingMethod< tnlGrid< 1, Real, Device, Index >, Anisotropy >::
+FastSweepingMethod< Meshes::Grid< 1, Real, Device, Index >, Anisotropy >::
 setMaxIterations( const IndexType& maxIterations )
 {
    
@@ -53,8 +53,8 @@ template< typename Real,
           typename Index,
           typename Anisotropy >
 void
-tnlFastSweepingMethod< tnlGrid< 1, Real, Device, Index >, Anisotropy >::
-solve( const MeshType& mesh,
+FastSweepingMethod< Meshes::Grid< 1, Real, Device, Index >, Anisotropy >::
+solve( const MeshPointer& mesh,
        const AnisotropyType& anisotropy,
        MeshFunctionType& u )
 {
diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod2D_impl.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod2D_impl.h
index cc2cfe9401c54aa7bff5e88f4c9d9ecf40f919cd..13e5c824912ef387c47ddaca1c5eb98acfe30fe6 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod2D_impl.h
+++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod2D_impl.h
@@ -19,8 +19,8 @@ template< typename Real,
           typename Device,
           typename Index,
           typename Anisotropy >
-tnlFastSweepingMethod< tnlGrid< 2, Real, Device, Index >, Anisotropy >::
-tnlFastSweepingMethod()
+FastSweepingMethod< Meshes::Grid< 2, Real, Device, Index >, Anisotropy >::
+FastSweepingMethod()
 : maxIterations( 1 )
 {
    
@@ -31,7 +31,7 @@ template< typename Real,
           typename Index,
           typename Anisotropy >
 const Index&
-tnlFastSweepingMethod< tnlGrid< 2, Real, Device, Index >, Anisotropy >::
+FastSweepingMethod< Meshes::Grid< 2, Real, Device, Index >, Anisotropy >::
 getMaxIterations() const
 {
    
@@ -42,7 +42,7 @@ template< typename Real,
           typename Index,
           typename Anisotropy >
 void
-tnlFastSweepingMethod< tnlGrid< 2, Real, Device, Index >, Anisotropy >::
+FastSweepingMethod< Meshes::Grid< 2, Real, Device, Index >, Anisotropy >::
 setMaxIterations( const IndexType& maxIterations )
 {
    
@@ -53,8 +53,8 @@ template< typename Real,
           typename Index,
           typename Anisotropy >
 void
-tnlFastSweepingMethod< tnlGrid< 2, Real, Device, Index >, Anisotropy >::
-solve( const MeshType& mesh,
+FastSweepingMethod< Meshes::Grid< 2, Real, Device, Index >, Anisotropy >::
+solve( const MeshPointer& mesh,
        const AnisotropyType& anisotropy,
        MeshFunctionType& u )
 {
@@ -66,18 +66,18 @@ solve( const MeshType& mesh,
    BaseType::initInterface( u, aux, interfaceMap );
    //aux.save( "aux-ini.tnl" );
 
-   typename MeshType::Cell cell( mesh );
+   typename MeshType::Cell cell( *mesh );
    
    IndexType iteration( 0 );
    while( iteration < this->maxIterations )
    {
 
       for( cell.getCoordinates().y() = 0;
-           cell.getCoordinates().y() < mesh.getDimensions().y();
+           cell.getCoordinates().y() < mesh->getDimensions().y();
            cell.getCoordinates().y()++ )
       {
          for( cell.getCoordinates().x() = 0;
-              cell.getCoordinates().x() < mesh.getDimensions().x();
+              cell.getCoordinates().x() < mesh->getDimensions().x();
               cell.getCoordinates().x()++ )
             {
                cell.refresh();
@@ -88,10 +88,10 @@ solve( const MeshType& mesh,
       //aux.save( "aux-1.tnl" );
 
       for( cell.getCoordinates().y() = 0;
-           cell.getCoordinates().y() < mesh.getDimensions().y();
+           cell.getCoordinates().y() < mesh->getDimensions().y();
            cell.getCoordinates().y()++ )
       {
-         for( cell.getCoordinates().x() = mesh.getDimensions().x() - 1;
+         for( cell.getCoordinates().x() = mesh->getDimensions().x() - 1;
               cell.getCoordinates().x() >= 0 ;
               cell.getCoordinates().x()-- )		
             {
@@ -103,12 +103,12 @@ solve( const MeshType& mesh,
       }
       //aux.save( "aux-2.tnl" );
 
-      for( cell.getCoordinates().y() = mesh.getDimensions().y() - 1;
+      for( cell.getCoordinates().y() = mesh->getDimensions().y() - 1;
            cell.getCoordinates().y() >= 0 ;
            cell.getCoordinates().y()-- )
          {
          for( cell.getCoordinates().x() = 0;
-              cell.getCoordinates().x() < mesh.getDimensions().x();
+              cell.getCoordinates().x() < mesh->getDimensions().x();
               cell.getCoordinates().x()++ )
             {
                //std::cerr << "3 -> ";
@@ -120,11 +120,11 @@ solve( const MeshType& mesh,
       //aux.save( "aux-3.tnl" );
 
 
-      for( cell.getCoordinates().y() = mesh.getDimensions().y() - 1;
+      for( cell.getCoordinates().y() = mesh->getDimensions().y() - 1;
            cell.getCoordinates().y() >= 0;
            cell.getCoordinates().y()-- )
          {
-         for( cell.getCoordinates().x() = mesh.getDimensions().x() - 1;
+         for( cell.getCoordinates().x() = mesh->getDimensions().x() - 1;
               cell.getCoordinates().x() >= 0 ;
               cell.getCoordinates().x()-- )		
             {
diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod3D_impl.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod3D_impl.h
index f252537148b9968b1a84421da7ca563169525ae7..6fad87c7ea63f50a13e7b9e38c65994998e2b441 100644
--- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod3D_impl.h
+++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod3D_impl.h
@@ -19,8 +19,8 @@ template< typename Real,
           typename Device,
           typename Index,
           typename Anisotropy >
-tnlFastSweepingMethod< tnlGrid< 3, Real, Device, Index >, Anisotropy >::
-tnlFastSweepingMethod()
+FastSweepingMethod< Meshes::Grid< 3, Real, Device, Index >, Anisotropy >::
+FastSweepingMethod()
 : maxIterations( 1 )
 {
    
@@ -31,7 +31,7 @@ template< typename Real,
           typename Index,
           typename Anisotropy >
 const Index&
-tnlFastSweepingMethod< tnlGrid< 3, Real, Device, Index >, Anisotropy >::
+FastSweepingMethod< Meshes::Grid< 3, Real, Device, Index >, Anisotropy >::
 getMaxIterations() const
 {
    
@@ -42,7 +42,7 @@ template< typename Real,
           typename Index,
           typename Anisotropy >
 void
-tnlFastSweepingMethod< tnlGrid< 3, Real, Device, Index >, Anisotropy >::
+FastSweepingMethod< Meshes::Grid< 3, Real, Device, Index >, Anisotropy >::
 setMaxIterations( const IndexType& maxIterations )
 {
    
@@ -53,8 +53,8 @@ template< typename Real,
           typename Index,
           typename Anisotropy >
 void
-tnlFastSweepingMethod< tnlGrid< 3, Real, Device, Index >, Anisotropy >::
-solve( const MeshType& mesh,
+FastSweepingMethod< Meshes::Grid< 3, Real, Device, Index >, Anisotropy >::
+solve( const MeshPointer& mesh,
        const AnisotropyType& anisotropy,
        MeshFunctionType& u )
 {
diff --git a/src/TNL/Functions/MeshFunction_impl.h b/src/TNL/Functions/MeshFunction_impl.h
index 7164170f1ab51a20bcdeaeb9011cb97bb97b4590..8704fc6868b364c444d086916aa870556f41ce83 100644
--- a/src/TNL/Functions/MeshFunction_impl.h
+++ b/src/TNL/Functions/MeshFunction_impl.h
@@ -312,7 +312,7 @@ MeshFunction< Mesh, MeshEntityDimension, Real >::
 getValue( const EntityType& meshEntity ) const
 {
    static_assert( EntityType::getEntityDimension() == MeshEntityDimension, "Calling with wrong EntityType -- entity dimensions do not match." );
-   return this->data.getValue( meshEntity.getIndex() );
+   return this->data.getElement( meshEntity.getIndex() );
 }
 
 template< typename Mesh,
@@ -325,7 +325,7 @@ setValue( const EntityType& meshEntity,
           const RealType& value )
 {
    static_assert( EntityType::getEntityDimension() == MeshEntityDimension, "Calling with wrong EntityType -- entity dimensions do not match." );
-   this->data.setValue( meshEntity.getIndex(), value );
+   this->data.setElement( meshEntity.getIndex(), value );
 }
 
 template< typename Mesh,
diff --git a/src/TNL/Math.h b/src/TNL/Math.h
index e6993069e55cf9a999df56a06588ca13ac5231e3..c546a58f78f94cac4ac882bf91c26c0bbc214236 100644
--- a/src/TNL/Math.h
+++ b/src/TNL/Math.h
@@ -75,7 +75,6 @@ ResultType max( const T1& a, const T2& b )
 #endif
 }
 
-
 /***
  * This function returns absolute value of given number.
  */
@@ -92,6 +91,45 @@ T abs( const T& n )
 #endif
 }
 
+/***
+ * This function returns argument of minimum of two numbers.
+ */
+template< typename T1, typename T2, typename ResultType = larger_type< T1, T2 > >
+__cuda_callable__ inline
+ResultType argMin( const T1& a, const T2& b )
+{
+   return ( a < b ) ?  a : b;
+}
+
+/***
+ * This function returns argument of maximum of two numbers.
+ */
+template< typename T1, typename T2, typename ResultType = larger_type< T1, T2 > >
+__cuda_callable__
+ResultType argMax( const T1& a, const T2& b )
+{
+   return ( a > b ) ?  a : b;   
+}
+
+/***
+ * This function returns argument of minimum of absolute values of two numbers.
+ */
+template< typename T1, typename T2, typename ResultType = larger_type< T1, T2 > >
+__cuda_callable__ inline
+ResultType argAbsMin( const T1& a, const T2& b )
+{
+   return ( TNL::abs( a ) < TNL::abs( b ) ) ?  a : b;
+}
+
+/***
+ * This function returns argument of maximum of absolute values of two numbers.
+ */
+template< typename T1, typename T2, typename ResultType = larger_type< T1, T2 > >
+__cuda_callable__
+ResultType argAbsMax( const T1& a, const T2& b )
+{
+   return ( TNL::abs( a ) > TNL::abs( b ) ) ?  a : b;   
+}
 
 template< typename T1, typename T2, typename ResultType = larger_type< T1, T2 > >
 __cuda_callable__ inline
diff --git a/src/TNL/Meshes/GridDetails/GridTraverser_impl.h b/src/TNL/Meshes/GridDetails/GridTraverser_impl.h
index 4443e161c3d3296f7a7924d927f00a3581248cec..b61837ae9e8673038c1ac263b0c94f37913dbdd5 100644
--- a/src/TNL/Meshes/GridDetails/GridTraverser_impl.h
+++ b/src/TNL/Meshes/GridDetails/GridTraverser_impl.h
@@ -381,23 +381,6 @@ GridTraverser2D(
    coordinates.x() = begin.x() + Devices::Cuda::getGlobalThreadIdx_x( gridIdx );
    coordinates.y() = begin.y() + Devices::Cuda::getGlobalThreadIdx_y( gridIdx );
    
-   /*if( processOnlyBoundaryEntities && 
-      ( GridEntity::getMeshDimension() == 2 || GridEntity::getMeshDimension() == 0 ) )
-   {
-      if( coordinates.x() == begin.x() || coordinates.x() == end.x() ||
-          coordinates.y() == begin.y() || coordinates.y() == end.y() )
-      {
-         GridEntity entity( *grid, coordinates, gridEntityParameters... );
-         entity.refresh();
-         EntitiesProcessor::processEntity
-         ( *grid,
-           *userData,
-           entity );
-      }
-      return;
-   }*/
-   
-   
    if( coordinates <= end )
    {
       GridEntity entity( *grid, coordinates, gridEntityParameters... );
@@ -504,7 +487,7 @@ processEntities(
 {
 #ifdef HAVE_CUDA
    if( processOnlyBoundaryEntities && 
-      ( GridEntity::getMeshDimension() == 2 || GridEntity::getMeshDimension() == 0 ) )
+       ( GridEntity::getEntityDimension() == 2 || GridEntity::getEntityDimension() == 0 ) )
    {
       dim3 cudaBlockSize( 256 );
       dim3 cudaBlocksCountAlongX, cudaGridsCountAlongX,
@@ -998,7 +981,7 @@ processEntities(
 {
 #ifdef HAVE_CUDA   
    if( processOnlyBoundaryEntities && 
-      ( GridEntity::getMeshDimension() == 3 || GridEntity::getMeshDimension() == 0 ) )
+       ( GridEntity::getEntityDimension() == 3 || GridEntity::getEntityDimension() == 0 ) )
    {
       dim3 cudaBlockSize( 16, 16 );
       const IndexType entitiesAlongX = end.x() - begin.x() + 1;
diff --git a/src/TNL/ParallelFor.h b/src/TNL/ParallelFor.h
index b1e81ca97bf66576875aa1d0dcb34ecc0298d5df..7002ecafe0fd7cc97271bf2fd885f3bfb50ac331 100644
--- a/src/TNL/ParallelFor.h
+++ b/src/TNL/ParallelFor.h
@@ -76,6 +76,8 @@ struct ParallelFor< Devices::Cuda >
 
          Devices::Cuda::synchronizeDevice();
          ParallelForKernel<<< gridSize, blockSize >>>( start, end, f, args... );
+         cudaDeviceSynchronize();
+         TNL_CHECK_CUDA_DEVICE;
       }
 #endif
    }
diff --git a/src/TNL/Problems/PDEProblem.h b/src/TNL/Problems/PDEProblem.h
index 2da2f9981f1bc12673096d278c20a4b3b1b7eedc..f8535261b95e4748dcd6e0d1236ec3474e4ee02b 100644
--- a/src/TNL/Problems/PDEProblem.h
+++ b/src/TNL/Problems/PDEProblem.h
@@ -13,6 +13,7 @@
 #include <TNL/Problems/Problem.h>
 #include <TNL/SharedPointer.h>
 #include <TNL/Matrices/SlicedEllpack.h>
+#include <TNL/Solvers/PDE/TimeDependentPDESolver.h>
 
 namespace TNL {
 namespace Problems {
@@ -38,11 +39,13 @@ class PDEProblem : public Problem< Real, Device, Index >
       typedef Containers::Vector< RealType, DeviceType, IndexType > MeshDependentDataType;
       typedef SharedPointer< MeshDependentDataType, DeviceType > MeshDependentDataPointer;
 
+      static constexpr bool isTimeDependent() { return true; };
+      
       /****
        * This means that the time stepper will be set from the command line arguments.
        */
       typedef void TimeStepper;
-
+      
       static String getTypeStatic();
 
       String getPrologHeader() const;
diff --git a/src/TNL/Solvers/DummyProblem.h b/src/TNL/Solvers/DummyProblem.h
index 6c036bb2dfccc7fd95654c0c577cd96368136df0..b34206865ef9948639590f193707094750aa429c 100644
--- a/src/TNL/Solvers/DummyProblem.h
+++ b/src/TNL/Solvers/DummyProblem.h
@@ -30,7 +30,12 @@ class DummyProblem
       typedef Containers::Vector< Real, Device, Index > DofVectorType;
       typedef Meshes::Grid< 1, Real, Device, Index > MeshType;
       typedef DofVectorType MeshDependentDataType;
+      
+      static constexpr bool isTimeDependent(){ return true; };      
 };
 
+class DummySolver
+{};
+
 } // namespace Solvers
 } // namespace TNL
diff --git a/src/TNL/Solvers/PDE/CMakeLists.txt b/src/TNL/Solvers/PDE/CMakeLists.txt
index 346e2acab6284faf682c3dfa53212ef2d5f1864d..283b1d4b3d389cae7e1905e1629fb2f32c89f7af 100644
--- a/src/TNL/Solvers/PDE/CMakeLists.txt
+++ b/src/TNL/Solvers/PDE/CMakeLists.txt
@@ -5,6 +5,9 @@ SET( headers BackwardTimeDiscretisation.h
              ExplicitUpdater.h
              LinearSystemAssembler.h
              NoTimeDiscretisation.h
+             PDESolver.h
+             PDESolver_impl.h
+             PDESolverTypeResolver.h
              TimeDependentPDESolver.h
              TimeDependentPDESolver_impl.h
              TimeIndependentPDESolver.h
diff --git a/src/TNL/Solvers/PDE/ExplicitTimeStepper.h b/src/TNL/Solvers/PDE/ExplicitTimeStepper.h
index 381803f2f9d2d8a2bfed04b3022a2bb5898fbb05..a4ae5927f9eb5ec68e33ca0d41d0198db25c957d 100644
--- a/src/TNL/Solvers/PDE/ExplicitTimeStepper.h
+++ b/src/TNL/Solvers/PDE/ExplicitTimeStepper.h
@@ -39,6 +39,8 @@ class ExplicitTimeStepper
    typedef SharedPointer< DofVectorType, DeviceType > DofVectorPointer;
    typedef SharedPointer< MeshDependentDataType, DeviceType > MeshDependentDataPointer;
    typedef IterativeSolverMonitor< RealType, IndexType > SolverMonitorType;
+   
+   static_assert( ProblemType::isTimeDependent(), "The problem is not time dependent." );
 
    ExplicitTimeStepper();
 
@@ -73,7 +75,7 @@ class ExplicitTimeStepper
                         DofVectorPointer& _u,
                         DofVectorPointer& _fu );
    
-   bool writeEpilog( Logger& logger );
+   bool writeEpilog( Logger& logger ) const;
 
    protected:
 
diff --git a/src/TNL/Solvers/PDE/ExplicitTimeStepper_impl.h b/src/TNL/Solvers/PDE/ExplicitTimeStepper_impl.h
index a2450e38d2acdd331da85a59ef3ca639c022660d..c2d6a03e4f22db6f3e2da1e9a1d7e0202ec82eec 100644
--- a/src/TNL/Solvers/PDE/ExplicitTimeStepper_impl.h
+++ b/src/TNL/Solvers/PDE/ExplicitTimeStepper_impl.h
@@ -196,7 +196,7 @@ template< typename Problem,
           template < typename OdeProblem > class OdeSolver >
 bool
 ExplicitTimeStepper< Problem, OdeSolver >::
-writeEpilog( Logger& logger )
+writeEpilog( Logger& logger ) const
 {
    logger.writeParameter< long long int >( "Iterations count:", this->allIterations );
    logger.writeParameter< const char* >( "Pre-iterate time:", "" );
diff --git a/src/TNL/Solvers/PDE/PDESolver.h b/src/TNL/Solvers/PDE/PDESolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..595b57c026abe9c7aab82283f6ffc8dff5a4a5ec
--- /dev/null
+++ b/src/TNL/Solvers/PDE/PDESolver.h
@@ -0,0 +1,69 @@
+/***************************************************************************
+                          PDESolver.h  -  description
+                             -------------------
+    begin                : Nov 11, 2017
+    copyright            : (C) 2017 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#pragma once
+
+#include <TNL/Object.h>
+#include <TNL/Timer.h>
+#include <TNL/Logger.h>
+#include <TNL/Config/ConfigDescription.h>
+#include <TNL/Config/ParameterContainer.h>
+#include <TNL/Solvers/IterativeSolverMonitor.h>
+
+namespace TNL {
+namespace Solvers {
+namespace PDE { 
+   
+template< typename Real,
+          typename Index >
+class PDESolver : public Object
+{
+   public:
+      using RealType = Real;
+      using IndexType = Index;
+      using SolverMonitorType = IterativeSolverMonitor< RealType, IndexType >;
+      
+      
+      PDESolver();
+
+      static void configSetup( Config::ConfigDescription& config,
+                               const String& prefix = "" );
+
+      bool setup( const Config::ParameterContainer& parameters,
+                 const String& prefix = "" );
+
+      bool writeProlog( Logger& logger,
+                        const Config::ParameterContainer& parameters );
+      
+      void setIoTimer( Timer& ioTimer);
+
+      void setComputeTimer( Timer& computeTimer );
+      
+      void setTotalTimer( Timer& totalTimer );
+      
+      void setSolverMonitor( SolverMonitorType& solverMonitor );
+      
+      SolverMonitorType& getSolverMonitor();
+
+      bool writeEpilog( Logger& logger ) const;      
+      
+   protected:
+
+      Timer *ioTimer, *computeTimer, *totalTimer;
+      
+      SolverMonitorType *solverMonitorPointer;
+};
+ 
+} // namespace PDE
+} // namespace Solvers
+} // namespace TNL
+
+#include <TNL/Solvers/PDE/PDESolver_impl.h>
+
diff --git a/src/TNL/Solvers/PDE/PDESolverTypeResolver.h b/src/TNL/Solvers/PDE/PDESolverTypeResolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..bf82fce66445ddefe5f7a2476c2eba60e6626eda
--- /dev/null
+++ b/src/TNL/Solvers/PDE/PDESolverTypeResolver.h
@@ -0,0 +1,52 @@
+/***************************************************************************
+                          PDESolverTypeResolver.h  -  description
+                             -------------------
+    begin                : Nov 10, 2017
+    copyright            : (C) 2017 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#pragma once
+
+#include <TNL/Solvers/PDE/TimeDependentPDESolver.h>
+#include <TNL/Solvers/PDE/TimeIndependentPDESolver.h>
+
+namespace TNL {
+namespace Solvers {
+namespace PDE { 
+   
+template< typename Problem,
+          typename DiscreteSolver,
+          typename TimeStepper,
+          bool TimeDependent = Problem::isTimeDependent() >
+class PDESolverTypeResolver
+{
+};
+  
+template< typename Problem,
+          typename DiscreteSolver,
+          typename TimeStepper >
+class PDESolverTypeResolver< Problem, DiscreteSolver, TimeStepper, true >
+{
+   public:
+      
+      using SolverType = TimeDependentPDESolver< Problem, DiscreteSolver, TimeStepper >;
+};
+
+template< typename Problem,
+          typename DiscreteSolver,
+          typename TimeStepper >
+class PDESolverTypeResolver< Problem, DiscreteSolver, TimeStepper, false >
+{
+   public:
+      
+      using SolverType = TimeIndependentPDESolver< Problem, DiscreteSolver >;
+};
+   
+ 
+} // namespace PDE
+} // namespace Solvers
+} // namespace TNL
+
diff --git a/src/TNL/Solvers/PDE/PDESolver_impl.h b/src/TNL/Solvers/PDE/PDESolver_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..16d6fc43ce457d02ac2bd813e2523be3abb2a712
--- /dev/null
+++ b/src/TNL/Solvers/PDE/PDESolver_impl.h
@@ -0,0 +1,121 @@
+/***************************************************************************
+                          PDESolver_impl.h  -  description
+                             -------------------
+    begin                : Nov 11, 2017
+    copyright            : (C) 2017 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#pragma once
+
+#include <TNL/Solvers/PDE/PDESolver.h>
+
+namespace TNL {
+namespace Solvers {
+namespace PDE { 
+
+template< typename Real,
+          typename Index >   
+PDESolver< Real, Index >::PDESolver()   
+: ioTimer( 0 ),
+  computeTimer( 0 ),
+  totalTimer( 0 ),
+  solverMonitorPointer( 0 )
+{
+}
+   
+template< typename Real,
+          typename Index >
+void
+PDESolver< Real, Index >::
+configSetup( Config::ConfigDescription& config,
+             const String& prefix )
+{
+}
+
+template< typename Real,
+          typename Index >
+bool
+PDESolver< Real, Index >::
+setup( const Config::ParameterContainer& parameters,
+       const String& prefix )
+{
+   return true;
+}
+
+template< typename Real,
+          typename Index >
+void
+PDESolver< Real, Index >::
+setSolverMonitor( typename PDESolver< Real, Index >::SolverMonitorType& solverMonitor )
+{
+   this->solverMonitorPointer = &solverMonitor;
+}
+
+template< typename Real,
+          typename Index >
+typename PDESolver< Real, Index >::SolverMonitorType&
+PDESolver< Real, Index >::
+getSolverMonitor()
+{
+   return *this->solverMonitorPointer;
+}
+
+template< typename Real,
+          typename Index >
+bool
+PDESolver< Real, Index >::
+writeProlog( Logger& logger,
+             const Config::ParameterContainer& parameters )
+{
+   logger.writeParameter< String >( "Real type:", "real-type", parameters, 0 );
+   logger.writeParameter< String >( "Index type:", "index-type", parameters, 0 );
+   logger.writeParameter< String >( "Device:", "device", parameters, 0 );
+   if( parameters.getParameter< String >( "device" ) == "host" )
+   {
+      if( Devices::Host::isOMPEnabled() )
+      {
+         logger.writeParameter< String >( "OMP enabled:", "yes", 1 );
+         logger.writeParameter< int >( "OMP threads:", Devices::Host::getMaxThreadsCount(), 1 );
+      }
+      else
+         logger.writeParameter< String >( "OMP enabled:", "no", 1 );
+   }
+   logger.writeSeparator();
+   logger.writeSystemInformation( parameters );
+   logger.writeSeparator();
+   logger.writeCurrentTime( "Started at:" );
+   logger.writeSeparator();
+   return true;
+}
+
+template< typename Real,
+          typename Index >
+void PDESolver< Real, Index >::
+setIoTimer( Timer& ioTimer )
+{
+   this->ioTimer = &ioTimer;
+}
+
+template< typename Real,
+          typename Index >
+void PDESolver< Real, Index >::
+setComputeTimer( Timer& computeTimer )
+{
+   this->computeTimer = &computeTimer;
+}
+
+template< typename Real,
+          typename Index >
+void PDESolver< Real, Index >::
+setTotalTimer( Timer& totalTimer )
+{
+   this->totalTimer = &totalTimer;
+}  
+   
+} // namespace PDE
+} // namespace Solvers
+} // namespace TNL
+   
\ No newline at end of file
diff --git a/src/TNL/Solvers/PDE/SemiImplicitTimeStepper.h b/src/TNL/Solvers/PDE/SemiImplicitTimeStepper.h
index b40ac527c01e831269e405c109fe60be2856f503..9d8b4c376cb82f531bcda6301d83acbedc3734f3 100644
--- a/src/TNL/Solvers/PDE/SemiImplicitTimeStepper.h
+++ b/src/TNL/Solvers/PDE/SemiImplicitTimeStepper.h
@@ -72,7 +72,7 @@ class SemiImplicitTimeStepper
                DofVectorPointer& dofVectorPointer,
                MeshDependentDataPointer& meshDependentData );
  
-   bool writeEpilog( Logger& logger );
+   bool writeEpilog( Logger& logger ) const;
 
    protected:
 
diff --git a/src/TNL/Solvers/PDE/SemiImplicitTimeStepper_impl.h b/src/TNL/Solvers/PDE/SemiImplicitTimeStepper_impl.h
index 1c324f8f944d43ae73cea4036ae252768112bc32..ba5fc1ae437ad75cec22c3301adec71183259dc5 100644
--- a/src/TNL/Solvers/PDE/SemiImplicitTimeStepper_impl.h
+++ b/src/TNL/Solvers/PDE/SemiImplicitTimeStepper_impl.h
@@ -245,7 +245,7 @@ template< typename Problem,
           typename LinearSystemSolver >
 bool
 SemiImplicitTimeStepper< Problem, LinearSystemSolver >::
-writeEpilog( Logger& logger )
+writeEpilog( Logger& logger ) const
 {
    logger.writeParameter< long long int >( "Iterations count:", this->allIterations );
    logger.writeParameter< const char* >( "Pre-iterate time:", "" );
diff --git a/src/TNL/Solvers/PDE/TimeDependentPDESolver.h b/src/TNL/Solvers/PDE/TimeDependentPDESolver.h
index c6a1baa174b64b985e9590372eb913eda784cd47..8b400c994e1b05d560e41c1123aaf3d73c702753 100644
--- a/src/TNL/Solvers/PDE/TimeDependentPDESolver.h
+++ b/src/TNL/Solvers/PDE/TimeDependentPDESolver.h
@@ -10,33 +10,38 @@
 
 #pragma once
 
-#include <TNL/Object.h>
 #include <TNL/Config/ConfigDescription.h>
 #include <TNL/Config/ParameterContainer.h>
 #include <TNL/Logger.h>
-#include <TNL/Timer.h>
 #include <TNL/SharedPointer.h>
+#include <TNL/Solvers/PDE/PDESolver.h>
 
 namespace TNL {
 namespace Solvers {
 namespace PDE {   
 
 template< typename Problem,
+          typename DiscreteSolver,
           typename TimeStepper >
-class TimeDependentPDESolver : public Object
+class TimeDependentPDESolver : public PDESolver< typename Problem::RealType, 
+                                                 typename Problem::IndexType >
 {
    public:
 
-      typedef typename TimeStepper::RealType RealType;
-      typedef typename TimeStepper::DeviceType DeviceType;
-      typedef typename TimeStepper::IndexType IndexType;
-      typedef Problem ProblemType;
+      using RealType = typename Problem::RealType;
+      using DeviceType = typename Problem::DeviceType;
+      using IndexType = typename Problem::IndexType;
+      using BaseType = PDESolver< RealType, IndexType >;
+      using ProblemType = Problem;
       typedef typename ProblemType::MeshType MeshType;
       typedef typename ProblemType::DofVectorType DofVectorType;
       typedef typename ProblemType::MeshDependentDataType MeshDependentDataType;
       typedef SharedPointer< MeshType, DeviceType > MeshPointer;
       typedef SharedPointer< DofVectorType, DeviceType > DofVectorPointer;
       typedef SharedPointer< MeshDependentDataType, DeviceType > MeshDependentDataPointer;
+      typedef IterativeSolverMonitor< typename Problem::RealType, typename Problem::IndexType > SolverMonitorType;
+      
+      static_assert( ProblemType::isTimeDependent(), "The problem is not time dependent." );
 
       TimeDependentPDESolver();
 
@@ -49,8 +54,6 @@ class TimeDependentPDESolver : public Object
       bool writeProlog( Logger& logger,
                         const Config::ParameterContainer& parameters );
 
-      void setTimeStepper( TimeStepper& timeStepper );
-
       void setProblem( ProblemType& problem );
 
       void setInitialTime( const RealType& initialT );
@@ -73,10 +76,6 @@ class TimeDependentPDESolver : public Object
 
       const RealType& getSnapshotPeriod() const;
 
-      void setIoTimer( Timer& ioTimer);
-
-      void setComputeTimer( Timer& computeTimer );
-
       bool solve();
 
       bool writeEpilog( Logger& logger ) const;
@@ -89,13 +88,13 @@ class TimeDependentPDESolver : public Object
 
       MeshDependentDataPointer meshDependentDataPointer;
 
-      TimeStepper* timeStepper;
-
-      RealType initialTime, finalTime, snapshotPeriod, timeStep, timeStepOrder;
-
+      TimeStepper timeStepper;
+      
+      DiscreteSolver discreteSolver;
+      
       ProblemType* problem;
 
-      Timer *ioTimer, *computeTimer;
+      RealType initialTime, finalTime, snapshotPeriod, timeStep, timeStepOrder;
 };
 
 } // namespace PDE
diff --git a/src/TNL/Solvers/PDE/TimeDependentPDESolver_impl.h b/src/TNL/Solvers/PDE/TimeDependentPDESolver_impl.h
index 6d20972d28c9f06cdf8da05e90567ac50f656e67..991789e8abf1898c5af82653de4875db1d5153e3 100644
--- a/src/TNL/Solvers/PDE/TimeDependentPDESolver_impl.h
+++ b/src/TNL/Solvers/PDE/TimeDependentPDESolver_impl.h
@@ -17,28 +17,28 @@ namespace Solvers {
 namespace PDE {   
 
 template< typename Problem,
+          typename DiscreteSolver,
           typename TimeStepper >
-TimeDependentPDESolver< Problem, TimeStepper >::
+TimeDependentPDESolver< Problem, DiscreteSolver, TimeStepper >::
 TimeDependentPDESolver()
-: timeStepper( 0 ),
+: problem( 0 ),
   initialTime( 0.0 ),
   finalTime( 0.0 ),
   snapshotPeriod( 0.0 ),
   timeStep( 1.0 ),
-  timeStepOrder( 0.0 ),
-  problem( 0 ),
-  ioTimer( 0 ),
-  computeTimer( 0 )
+  timeStepOrder( 0.0 )
 {
 }
 
 template< typename Problem,
+          typename DiscreteSolver,   
           typename TimeStepper >
 void
-TimeDependentPDESolver< Problem, TimeStepper >::
+TimeDependentPDESolver< Problem, DiscreteSolver, TimeStepper >::
 configSetup( Config::ConfigDescription& config,
              const String& prefix )
 {
+   BaseType::configSetup( config, prefix );
    config.addEntry< String >( prefix + "initial-condition", "File name with the initial condition.", "init.tnl" );
    config.addRequiredEntry< double >( prefix + "final-time", "Stop time of the time dependent problem." );
    config.addEntry< double >( prefix + "initial-time", "Initial time of the time dependent problem.", 0 );
@@ -48,12 +48,16 @@ configSetup( Config::ConfigDescription& config,
 }
 
 template< typename Problem,
+          typename DiscreteSolver,   
           typename TimeStepper >
 bool
-TimeDependentPDESolver< Problem, TimeStepper >::
+TimeDependentPDESolver< Problem, DiscreteSolver, TimeStepper >::
 setup( const Config::ParameterContainer& parameters,
        const String& prefix )
 {
+   
+   BaseType::setup( parameters, prefix );
+   
    /****
     * Load the mesh from the mesh file
     */
@@ -108,16 +112,31 @@ setup( const Config::ParameterContainer& parameters,
    this->setSnapshotPeriod( parameters.getParameter< double >( "snapshot-period" ) );
    this->setTimeStep( parameters.getParameter< double >( "time-step") );
    this->setTimeStepOrder( parameters.getParameter< double >( "time-step-order" ) );
+
+   /****
+    * Set-up the discrete solver
+    */
+   if( ! this->discreteSolver.setup( parameters ) )
+      return false;
+   
+   /****
+    * Set-up the time stepper
+    */
+   if( ! this->timeStepper.setup( parameters ) )
+      return false;
+   this->timeStepper.setSolver( this->discreteSolver );
+   this->timeStepper.setSolverMonitor( *this->solverMonitorPointer );      
    return true;
 }
 
 template< typename Problem,
+          typename DiscreteSolver,   
           typename TimeStepper >
 bool
-TimeDependentPDESolver< Problem, TimeStepper >::
+TimeDependentPDESolver< Problem, DiscreteSolver, TimeStepper >::
 writeProlog( Logger& logger,
              const Config::ParameterContainer& parameters )
-{
+{   
    logger.writeHeader( problem->getPrologHeader() );
    problem->writeProlog( logger, parameters );
    logger.writeSeparator();
@@ -145,68 +164,44 @@ writeProlog( Logger& logger,
    logger.writeParameter< int >( "Maximal number of iterations:", "max-iterations", parameters );
    logger.writeParameter< int >( "Minimal number of iterations:", "min-iterations", parameters );
    logger.writeSeparator();
-   logger.writeParameter< String >( "Real type:", "real-type", parameters, 0 );
-   logger.writeParameter< String >( "Index type:", "index-type", parameters, 0 );
-   logger.writeParameter< String >( "Device:", "device", parameters, 0 );
-   if( parameters.getParameter< String >( "device" ) == "host" )
-   {
-      if( Devices::Host::isOMPEnabled() )
-      {
-         logger.writeParameter< String >( "OMP enabled:", "yes", 1 );
-         logger.writeParameter< int >( "OMP threads:", Devices::Host::getMaxThreadsCount(), 1 );
-      }
-      else
-         logger.writeParameter< String >( "OMP enabled:", "no", 1 );
-   }
-   logger.writeSeparator();
-   logger.writeSystemInformation( parameters );
-   logger.writeSeparator();
-   logger.writeCurrentTime( "Started at:" );
-   logger.writeSeparator();
-   return true;
+   return BaseType::writeProlog( logger, parameters );
 }
 
 template< typename Problem,
+          typename DiscreteSolver,   
           typename TimeStepper >
 void
-TimeDependentPDESolver< Problem, TimeStepper >::
-setTimeStepper( TimeStepper& timeStepper )
-{
-   this->timeStepper = &timeStepper;
-}
-
-template< typename Problem,
-          typename TimeStepper >
-void
-TimeDependentPDESolver< Problem, TimeStepper >::
+TimeDependentPDESolver< Problem, DiscreteSolver, TimeStepper >::
 setProblem( ProblemType& problem )
 {
    this->problem = &problem;
 }
 
 template< typename Problem,
+          typename DiscreteSolver,   
           typename TimeStepper >
 void
-TimeDependentPDESolver< Problem, TimeStepper >::
+TimeDependentPDESolver< Problem, DiscreteSolver, TimeStepper >::
 setInitialTime( const RealType& initialTime )
 {
    this->initialTime = initialTime;
 }
 
 template< typename Problem,
+          typename DiscreteSolver,   
           typename TimeStepper >
-const typename TimeStepper :: RealType&
-TimeDependentPDESolver< Problem, TimeStepper >::
+const typename Problem::RealType&
+TimeDependentPDESolver< Problem, DiscreteSolver, TimeStepper >::
 getInitialTime() const
 {
    return this->initialTime;
 }
 
-
 template< typename Problem,
+          typename DiscreteSolver,   
           typename TimeStepper >
 bool
-TimeDependentPDESolver< Problem, TimeStepper >::
+TimeDependentPDESolver< Problem, DiscreteSolver, TimeStepper >::
 setFinalTime( const RealType& finalTime )
 {
    if( finalTime <= this->initialTime )
@@ -219,18 +214,20 @@ setFinalTime( const RealType& finalTime )
 }
 
 template< typename Problem,
+          typename DiscreteSolver,   
           typename TimeStepper >
-const typename TimeStepper :: RealType&
-TimeDependentPDESolver< Problem, TimeStepper >::
+const typename Problem::RealType&
+TimeDependentPDESolver< Problem, DiscreteSolver, TimeStepper >::
 getFinalTime() const
 {
    return this->finalTime;
 }
 
 template< typename Problem,
+          typename DiscreteSolver,   
           typename TimeStepper >
 bool
-TimeDependentPDESolver< Problem, TimeStepper >::
+TimeDependentPDESolver< Problem, DiscreteSolver, TimeStepper >::
 setSnapshotPeriod( const RealType& period )
 {
    if( period <= 0 )
@@ -243,18 +240,20 @@ setSnapshotPeriod( const RealType& period )
 }
 
 template< typename Problem,
+          typename DiscreteSolver,   
           typename TimeStepper >
-const typename TimeStepper::RealType&
-TimeDependentPDESolver< Problem, TimeStepper >::
+const typename Problem::RealType&
+TimeDependentPDESolver< Problem, DiscreteSolver, TimeStepper >::
 getSnapshotPeriod() const
 {
    return this->snapshotPeriod;
 }
 
 template< typename Problem,
+          typename DiscreteSolver,   
           typename TimeStepper >
 bool
-TimeDependentPDESolver< Problem, TimeStepper >::
+TimeDependentPDESolver< Problem, DiscreteSolver, TimeStepper >::
 setTimeStep( const RealType& timeStep )
 {
    if( timeStep <= 0 )
@@ -267,18 +266,20 @@ setTimeStep( const RealType& timeStep )
 }
  
 template< typename Problem,
+          typename DiscreteSolver,   
           typename TimeStepper >
-const typename TimeStepper::RealType&
-TimeDependentPDESolver< Problem, TimeStepper >::
+const typename Problem::RealType&
+TimeDependentPDESolver< Problem, DiscreteSolver, TimeStepper >::
 getTimeStep() const
 {
    return this->timeStep;
 }
 
 template< typename Problem,
+          typename DiscreteSolver,   
           typename TimeStepper >
 bool
-TimeDependentPDESolver< Problem, TimeStepper >::
+TimeDependentPDESolver< Problem, DiscreteSolver, TimeStepper >::
 setTimeStepOrder( const RealType& timeStepOrder )
 {
    if( timeStepOrder < 0 )
@@ -291,32 +292,22 @@ setTimeStepOrder( const RealType& timeStepOrder )
 }
 
 template< typename Problem,
+          typename DiscreteSolver,   
           typename TimeStepper >
-const typename TimeStepper::RealType&
-TimeDependentPDESolver< Problem, TimeStepper >::
+const typename Problem::RealType&
+TimeDependentPDESolver< Problem, DiscreteSolver, TimeStepper >::
 getTimeStepOrder() const
 {
    return this->timeStepOrder;
 }
 
-template< typename Problem, typename TimeStepper >
-void TimeDependentPDESolver< Problem, TimeStepper > :: setIoTimer( Timer& ioTimer )
-{
-   this->ioTimer = &ioTimer;
-}
-
-template< typename Problem, typename TimeStepper >
-void TimeDependentPDESolver< Problem, TimeStepper > :: setComputeTimer( Timer& computeTimer )
-{
-   this->computeTimer = &computeTimer;
-}
-
-template< typename Problem, typename TimeStepper >
+template< typename Problem,
+          typename DiscreteSolver,   
+          typename TimeStepper >
 bool
-TimeDependentPDESolver< Problem, TimeStepper >::
+TimeDependentPDESolver< Problem, DiscreteSolver, TimeStepper >::
 solve()
 {
-   TNL_ASSERT_TRUE( timeStepper, "No time stepper was set in PDESolver." );
    TNL_ASSERT_TRUE( problem, "No problem was set in PDESolver." );
 
    if( snapshotPeriod == 0 )
@@ -343,14 +334,14 @@ solve()
    /****
     * Initialize the time stepper
     */
-   this->timeStepper->setProblem( * ( this->problem ) );
-   this->timeStepper->init( this->meshPointer );
-   this->timeStepper->setTimeStep( this->timeStep * std::pow( this->meshPointer.getData().getSmallestSpaceStep(), this->timeStepOrder ) );
+   this->timeStepper.setProblem( * ( this->problem ) );
+   this->timeStepper.init( this->meshPointer );
+   this->timeStepper.setTimeStep( this->timeStep * std::pow( this->meshPointer.getData().getSmallestSpaceStep(), this->timeStepOrder ) );
    while( step < allSteps )
    {
       RealType tau = min( this->snapshotPeriod,
                           this->finalTime - t );
-      if( ! this->timeStepper->solve( t, t + tau, this->meshPointer, this->dofsPointer, this->meshDependentDataPointer ) )
+      if( ! this->timeStepper.solve( t, t + tau, this->meshPointer, this->dofsPointer, this->meshDependentDataPointer ) )
          return false;
       step ++;
       t += tau;
@@ -366,15 +357,20 @@ solve()
       this->computeTimer->start();
    }
    this->computeTimer->stop();
+   
+   this->solverMonitorPointer->stopMainLoop();
+   
    return true;
 }
 
-template< typename Problem, typename TimeStepper >
+template< typename Problem,
+          typename DiscreteSolver,   
+          typename TimeStepper >
 bool
-TimeDependentPDESolver< Problem, TimeStepper >::
+TimeDependentPDESolver< Problem, DiscreteSolver, TimeStepper >::
 writeEpilog( Logger& logger ) const
 {
-   return ( this->timeStepper->writeEpilog( logger ) &&
+   return ( this->timeStepper.writeEpilog( logger ) &&
       this->problem->writeEpilog( logger ) );
 }
 
diff --git a/src/TNL/Solvers/PDE/TimeIndependentPDESolver.h b/src/TNL/Solvers/PDE/TimeIndependentPDESolver.h
index 41b6a47b5135fa74822804b5551f0dabbc81a3d4..996dab830698fbbb825b33fdd828d1c1495d252a 100644
--- a/src/TNL/Solvers/PDE/TimeIndependentPDESolver.h
+++ b/src/TNL/Solvers/PDE/TimeIndependentPDESolver.h
@@ -17,14 +17,21 @@
 
 #pragma once
 
-#include <core/tnlObject.h>
-#include <config/tnlConfigDescription.h>
+#include <TNL/Object.h>
+#include <TNL/Logger.h>
+#include <TNL/Config/ConfigDescription.h>
 #include <TNL/Config/ParameterContainer.h>
-#include <solvers/tnlSolverMonitor.h>
-#include <core/tnlLogger.h>
+#include <TNL/Solvers/PDE/PDESolver.h>
 
-template< typename Problem >
-class tnlTimeIndependentPDESolver : public tnlObject
+
+namespace TNL {
+namespace Solvers {   
+namespace PDE {
+
+template< typename Problem,
+          typename DiscreteSolver >
+class TimeIndependentPDESolver : public PDESolver< typename Problem::RealType,
+                                                   typename Problem::IndexType >
 {
    public:
 
@@ -34,42 +41,46 @@ class tnlTimeIndependentPDESolver : public tnlObject
       typedef typename ProblemType::IndexType IndexType;
       typedef typename ProblemType::MeshType MeshType;
       typedef typename ProblemType::DofVectorType DofVectorType;
-      typedef typename ProblemType::MeshDependentDataType MeshDependentDataType;      
+      typedef typename ProblemType::MeshDependentDataType MeshDependentDataType; 
+      typedef SharedPointer< MeshType, DeviceType > MeshPointer;
+      typedef SharedPointer< DofVectorType, DeviceType > DofVectorPointer;
+      typedef SharedPointer< MeshDependentDataType, DeviceType > MeshDependentDataPointer;
 
-      tnlTimeIndependentPDESolver();
+      TimeIndependentPDESolver();
 
-      static void configSetup( tnlConfigDescription& config,
+      static void configSetup( Config::ConfigDescription& config,
                                const String& prefix = "" );
 
       bool setup( const Config::ParameterContainer& parameters,
                   const String& prefix = "" );
 
-      bool writeProlog( tnlLogger& logger,
+      bool writeProlog( Logger& logger,
                         const Config::ParameterContainer& parameters );
 
 
       void setProblem( ProblemType& problem );
 
-      void setComputeTimer( tnlTimer& computeTimer );
-      
-      void setIoTimer( tnlTimer& ioTimer );
-
       bool solve();
 
-      bool writeEpilog( tnlLogger& logger ) const;
+      bool writeEpilog( Logger& logger ) const;
 
    protected:
 
-      MeshType mesh;
+      MeshPointer mesh;
 
-      DofVectorType dofs;
+      DofVectorPointer dofs;
 
-      MeshDependentDataType meshDependentData;
+      MeshDependentDataPointer meshDependentData;
+      
+      DiscreteSolver discreteSolver;
 
       ProblemType* problem;
-
-      tnlTimer *computeTimer;
 };
 
-#include <solvers/pde/tnlTimeIndependentPDESolver_impl.h>
+
+} // namespace PDE
+} // namespace Solvers
+} // namespace TNL
+
+#include <TNL/Solvers/PDE/TimeIndependentPDESolver_impl.h>
 
diff --git a/src/TNL/Solvers/PDE/TimeIndependentPDESolver_impl.h b/src/TNL/Solvers/PDE/TimeIndependentPDESolver_impl.h
index 957f710300bbc373b184dbd5736028377b68a7ba..9190fd8ed8a8c34a223a218215704d1540391ca0 100644
--- a/src/TNL/Solvers/PDE/TimeIndependentPDESolver_impl.h
+++ b/src/TNL/Solvers/PDE/TimeIndependentPDESolver_impl.h
@@ -1,5 +1,5 @@
 /***************************************************************************
-                          tnlTimeIndependentPDESolver_impl.h  -  description
+                          TimeIndependentPDESolver_impl.h  -  description
                              -------------------
     begin                : Jan 15, 2013
     copyright            : (C) 2013 by Tomas Oberhuber
@@ -17,27 +17,34 @@
 
 #pragma once 
 
-#include <solvers/pde/tnlTimeIndependentPDESolver.h>
+#include <TNL/Solvers/PDE/TimeIndependentPDESolver.h>
 
-template< typename Problem >
-tnlTimeIndependentPDESolver< Problem >::
-tnlTimeIndependentPDESolver()
-: problem( 0 ),
-  computeTimer( 0 )
+namespace TNL {
+namespace Solvers {
+namespace PDE {   
+
+
+template< typename Problem,
+          typename DiscreteSolver >
+TimeIndependentPDESolver< Problem, DiscreteSolver >::
+TimeIndependentPDESolver()
+: problem( 0 )
 {
 }
 
-template< typename Problem >
+template< typename Problem,
+          typename DiscreteSolver >
 void
-tnlTimeIndependentPDESolver< Problem >::
-configSetup( tnlConfigDescription& config,
+TimeIndependentPDESolver< Problem, DiscreteSolver >::
+configSetup( Config::ConfigDescription& config,
              const String& prefix )
 {
 }
 
-template< typename Problem >
+template< typename Problem,
+          typename DiscreteSolver >
 bool
-tnlTimeIndependentPDESolver< Problem >::
+TimeIndependentPDESolver< Problem, DiscreteSolver >::
 setup( const Config::ParameterContainer& parameters,
        const String& prefix )
 {
@@ -45,22 +52,22 @@ setup( const Config::ParameterContainer& parameters,
     * Load the mesh from the mesh file
     */
    const String& meshFile = parameters.getParameter< String >( "mesh" );
-   cout << "Loading a mesh from the file " << meshFile << "...";
-   if( ! this->mesh.load( meshFile ) )
+   std::cout << "Loading a mesh from the file " << meshFile << "...";
+   if( ! this->mesh->load( meshFile ) )
    {
-      cerr << endl;
-      cerr << "I am not able to load the mesh from the file " << meshFile << "." << endl;
-      cerr << " You may create it with tools like tnl-grid-setup or tnl-mesh-convert." << endl;
+      std::cerr << std::endl;
+      std::cerr << "I am not able to load the mesh from the file " << meshFile << "." << std::endl;
+      std::cerr << " You may create it with tools like tnl-grid-setup or tnl-mesh-convert." << std::endl;
       return false;
    }
-   cout << " [ OK ] " << endl;
+   std::cout << " [ OK ] " << std::endl;
 
    /****
     * Set DOFs (degrees of freedom)
     */
    TNL_ASSERT_GT( problem->getDofs( this->mesh ), 0, "number of DOFs must be positive" );
-   this->dofs.setSize( problem->getDofs( this->mesh ) );
-   this->dofs.setValue( 0.0 );
+   this->dofs->setSize( problem->getDofs( this->mesh ) );
+   this->dofs->setValue( 0.0 );
    this->problem->bindDofs( this->mesh, this->dofs );   
    
    /****
@@ -72,30 +79,29 @@ setup( const Config::ParameterContainer& parameters,
    /***
     * Set-up the initial condition
     */
-   cout << "Setting up the initial condition ... ";
+   std::cout << "Setting up the initial condition ... ";
    typedef typename Problem :: DofVectorType DofVectorType;
-   if( ! this->problem->setInitialData( parameters, this->mesh, this->dofs, this->meshDependentData ) )
+   if( ! this->problem->setInitialCondition( parameters, this->mesh, this->dofs, this->meshDependentData ) )
       return false;
-   cout << " [ OK ]" << endl;
+   std::cout << " [ OK ]" << std::endl;
    
    return true;
 }
 
-template< typename Problem >
+template< typename Problem,
+          typename DiscreteSolver >
 bool
-tnlTimeIndependentPDESolver< Problem >::
-writeProlog( tnlLogger& logger,
+TimeIndependentPDESolver< Problem, DiscreteSolver >::
+writeProlog( Logger& logger,
              const Config::ParameterContainer& parameters )
 {
    logger.writeHeader( problem->getPrologHeader() );
    problem->writeProlog( logger, parameters );
    logger.writeSeparator();
-   mesh.writeProlog( logger );
+   mesh->writeProlog( logger );
    logger.writeSeparator();
    const String& solverName = parameters. getParameter< String >( "discrete-solver" );
    logger.writeParameter< String >( "Discrete solver:", "discrete-solver", parameters );
-   if( solverName == "merson" )
-      logger.writeParameter< double >( "Adaptivity:", "merson-adaptivity", parameters, 1 );
    if( solverName == "sor" )
       logger.writeParameter< double >( "Omega:", "sor-omega", parameters, 1 );
    if( solverName == "gmres" )
@@ -116,29 +122,19 @@ writeProlog( tnlLogger& logger,
    return true;
 }
 
-template< typename Problem >
+template< typename Problem,
+          typename DiscreteSolver >
 void
-tnlTimeIndependentPDESolver< Problem >::
+TimeIndependentPDESolver< Problem, DiscreteSolver >::
 setProblem( ProblemType& problem )
 {
    this->problem = &problem;
 }
 
-template< typename Problem >
-void tnlTimeIndependentPDESolver< Problem > :: setIoTimer( tnlTimer& ioTimer )
-{
-  // this->ioTimer = &ioTimer;
-}
-
-template< typename Problem >
-void tnlTimeIndependentPDESolver< Problem > :: setComputeTimer( tnlTimer& computeTimer )
-{
-   this->computeTimer = &computeTimer;
-}
-
-template< typename Problem >
+template< typename Problem,
+          typename DiscreteSolver >
 bool
-tnlTimeIndependentPDESolver< Problem >::
+TimeIndependentPDESolver< Problem, DiscreteSolver >::
 solve()
 {
    TNL_ASSERT_TRUE( problem, "No problem was set in tnlPDESolver." );
@@ -154,10 +150,15 @@ solve()
    return true;
 }
 
-template< typename Problem >
+template< typename Problem,
+          typename DiscreteSolver >
 bool
-tnlTimeIndependentPDESolver< Problem >::
-writeEpilog( tnlLogger& logger ) const
+TimeIndependentPDESolver< Problem, DiscreteSolver >::
+writeEpilog( Logger& logger ) const
 {
    return this->problem->writeEpilog( logger );
 }
+
+} // namespace PDE
+} // namespace Solvers
+} // namespace TNL
diff --git a/src/TNL/Solvers/SolverConfig_impl.h b/src/TNL/Solvers/SolverConfig_impl.h
index e66397310414edbecc68d7a1090e2186b952f851..69b895b50164449e8aa255ea52800d314e3e0da5 100644
--- a/src/TNL/Solvers/SolverConfig_impl.h
+++ b/src/TNL/Solvers/SolverConfig_impl.h
@@ -89,7 +89,8 @@ bool SolverConfig< ConfigTag, ProblemConfig >::configSetup( Config::ConfigDescri
     */
    config.addDelimiter( " === Time discretisation parameters ==== " );
    typedef PDE::ExplicitTimeStepper< DummyProblemType, ODE::Euler > ExplicitTimeStepper;
-   PDE::TimeDependentPDESolver< DummyProblemType, ExplicitTimeStepper >::configSetup( config );
+   typedef Solvers::DummySolver DiscreteSolver;
+   PDE::TimeDependentPDESolver< DummyProblemType, DiscreteSolver, ExplicitTimeStepper >::configSetup( config );
    ExplicitTimeStepper::configSetup( config );
    if( ConfigTagTimeDiscretisation< ConfigTag, ExplicitTimeDiscretisationTag >::enabled ||
        ConfigTagTimeDiscretisation< ConfigTag, SemiImplicitTimeDiscretisationTag >::enabled ||
diff --git a/src/TNL/Solvers/SolverStarter_impl.h b/src/TNL/Solvers/SolverStarter_impl.h
index a8fae47d95980bf6f1414c6adeb38c1413a82fa0..a7e810e1700db1c80045a7f30b97f387f6306d2e 100644
--- a/src/TNL/Solvers/SolverStarter_impl.h
+++ b/src/TNL/Solvers/SolverStarter_impl.h
@@ -32,14 +32,21 @@
 #include <TNL/Solvers/PDE/ExplicitTimeStepper.h>
 #include <TNL/Solvers/PDE/SemiImplicitTimeStepper.h>
 #include <TNL/Solvers/PDE/TimeDependentPDESolver.h>
+#include <TNL/Solvers/PDE/PDESolverTypeResolver.h>
 
 namespace TNL {
 namespace Solvers {   
 
+template< typename Problem,
+          typename ConfigTag,
+          bool TimeDependent = Problem::isTimeDependent() >
+class TimeDependencyResolver
+{};
+   
 template< typename Problem,
           typename ConfigTag,
           typename TimeStepper = typename Problem::TimeStepper >
-class tnlUserDefinedTimeDiscretisationSetter;
+class UserDefinedTimeDiscretisationSetter;
 
 template< typename Problem,
           typename TimeDiscretisation,
@@ -58,7 +65,7 @@ template< typename Problem,
           template<typename, typename, typename> class Preconditioner,
           typename ConfigTag,
           bool enabled = ConfigTagSemiImplicitSolver< ConfigTag, SemiImplicitSolver >::enabled >
-class SolverStarterSemiImplicitSolverSetter{};
+class SolverStarterLinearSolverSetter{};
 
 template< typename Problem,
           typename SemiImplicitSolverTag,
@@ -83,13 +90,39 @@ bool SolverStarter< ConfigTag > :: run( const Config::ParameterContainer& parame
        ! Devices::Cuda::setup( parameters ) )
       return false;
    Problem problem;
-   return tnlUserDefinedTimeDiscretisationSetter< Problem, ConfigTag >::run( problem, parameters );
+   //return UserDefinedTimeDiscretisationSetter< Problem, ConfigTag >::run( problem, parameters );
+   return TimeDependencyResolver< Problem, ConfigTag >::run( problem, parameters );
 }
 
+template< typename Problem,
+          typename ConfigTag>
+class TimeDependencyResolver< Problem, ConfigTag, true >
+{
+   public:
+      static bool run( Problem& problem,
+                       const Config::ParameterContainer& parameters )
+      {
+         return UserDefinedTimeDiscretisationSetter< Problem, ConfigTag >::run( problem, parameters );
+      }
+};
+
+template< typename Problem,
+          typename ConfigTag>
+class TimeDependencyResolver< Problem, ConfigTag, false >
+{
+   public:
+      static bool run( Problem& problem,
+                       const Config::ParameterContainer& parameters )
+      {
+         // TODO: This should be improved - at least rename to LinearSolverSetter
+         return SolverStarterTimeDiscretisationSetter< Problem, SemiImplicitTimeDiscretisationTag, ConfigTag, true >::run( problem, parameters );   
+      }
+};
+
 template< typename Problem,
           typename ConfigTag,
           typename TimeStepper >
-class tnlUserDefinedTimeDiscretisationSetter
+class UserDefinedTimeDiscretisationSetter
 {
    public:
       static bool run( Problem& problem,
@@ -108,7 +141,7 @@ class tnlUserDefinedTimeDiscretisationSetter
 
 template< typename Problem,
           typename ConfigTag >
-class tnlUserDefinedTimeDiscretisationSetter< Problem, ConfigTag, void >
+class UserDefinedTimeDiscretisationSetter< Problem, ConfigTag, void >
 {
    public:
       static bool run( Problem& problem,
@@ -291,11 +324,11 @@ class SolverStarterPreconditionerSetter
          const String& preconditioner = parameters.getParameter< String>( "preconditioner" );
 
          if( preconditioner == "none" )
-            return SolverStarterSemiImplicitSolverSetter< Problem, SemiImplicitSolverTag, Linear::Preconditioners::Dummy, ConfigTag >::run( problem, parameters );
+            return SolverStarterLinearSolverSetter< Problem, SemiImplicitSolverTag, Linear::Preconditioners::Dummy, ConfigTag >::run( problem, parameters );
          if( preconditioner == "diagonal" )
-            return SolverStarterSemiImplicitSolverSetter< Problem, SemiImplicitSolverTag, Linear::Preconditioners::Diagonal, ConfigTag >::run( problem, parameters );
+            return SolverStarterLinearSolverSetter< Problem, SemiImplicitSolverTag, Linear::Preconditioners::Diagonal, ConfigTag >::run( problem, parameters );
          if( preconditioner == "ilu0" )
-            return SolverStarterSemiImplicitSolverSetter< Problem, SemiImplicitSolverTag, Linear::Preconditioners::ILU0, ConfigTag >::run( problem, parameters );
+            return SolverStarterLinearSolverSetter< Problem, SemiImplicitSolverTag, Linear::Preconditioners::ILU0, ConfigTag >::run( problem, parameters );
 
          std::cerr << "Unknown preconditioner " << preconditioner << ". It can be only: none, diagonal, ilu0." << std::endl;
          return false;
@@ -306,7 +339,7 @@ template< typename Problem,
           typename SemiImplicitSolverTag,
           template<typename, typename, typename> class Preconditioner,
           typename ConfigTag >
-class SolverStarterSemiImplicitSolverSetter< Problem, SemiImplicitSolverTag, Preconditioner, ConfigTag, false >
+class SolverStarterLinearSolverSetter< Problem, SemiImplicitSolverTag, Preconditioner, ConfigTag, false >
 {
    public:
       static bool run( Problem& problem,
@@ -321,7 +354,7 @@ template< typename Problem,
           typename SemiImplicitSolverTag,
           template<typename, typename, typename> class Preconditioner,
           typename ConfigTag >
-class SolverStarterSemiImplicitSolverSetter< Problem, SemiImplicitSolverTag, Preconditioner, ConfigTag, true >
+class SolverStarterLinearSolverSetter< Problem, SemiImplicitSolverTag, Preconditioner, ConfigTag, true >
 {
    public:
       static bool run( Problem& problem,
@@ -349,22 +382,11 @@ bool SolverStarter< ConfigTag > :: runPDESolver( Problem& problem,
 {
    this->totalTimer.reset();
    this->totalTimer.start();
+   
+   using SolverMonitorType = IterativeSolverMonitor< typename Problem::RealType,
+                                                     typename Problem::IndexType >;
+   SolverMonitorType solverMonitor, *solverMonitorPointer( &solverMonitor );
  
-   /****
-    * Set-up the discrete solver
-    */
-   DiscreteSolver discreteSolver;
-   if( ! discreteSolver.setup( parameters ) )
-      return false;
-
-   /****
-    * Set-up the time stepper
-    */
-   TimeStepper timeStepper;
-   if( ! timeStepper.setup( parameters ) )
-      return false;
-   timeStepper.setSolver( discreteSolver );
-
    /****
     * Open the log file
     */
@@ -378,14 +400,25 @@ bool SolverStarter< ConfigTag > :: runPDESolver( Problem& problem,
    /****
     * Set-up the PDE solver
     */
-   PDE::TimeDependentPDESolver< Problem, TimeStepper > solver;
+   //PDE::TimeDependentPDESolver< Problem, TimeStepper > solver;
+   typename PDE::PDESolverTypeResolver< Problem, DiscreteSolver, TimeStepper >::SolverType solver;
+   solver.setComputeTimer( this->computeTimer );
+   solver.setIoTimer( this->ioTimer );
+   solver.setTotalTimer( this->totalTimer );
+
+   if( problem.getSolverMonitor() )
+      solverMonitorPointer = ( SolverMonitorType* ) problem.getSolverMonitor();
+   solverMonitorPointer->setVerbose( parameters.getParameter< int >( "verbose" ) );
+   solverMonitorPointer->setTimer( this->totalTimer );
+   solver.setSolverMonitor( *solverMonitorPointer );
+   
    // catching exceptions ala gtest:
    // https://github.com/google/googletest/blob/59c795ce08be0c8b225bc894f8da6c7954ea5c14/googletest/src/gtest.cc#L2409-L2431
    const int catch_exceptions = parameters.getParameter< bool >( "catch-exceptions" );
    if( catch_exceptions ) {
       try {
          solver.setProblem( problem );
-         solver.setTimeStepper( timeStepper );
+         //solver.setTimeStepper( timeStepper ); // TODO: BETTER FIX: This does not make sense for time independent problem
          if( ! solver.setup( parameters ) )
             return false;
       }
@@ -402,7 +435,7 @@ bool SolverStarter< ConfigTag > :: runPDESolver( Problem& problem,
    }
    else {
       solver.setProblem( problem );
-      solver.setTimeStepper( timeStepper );
+      //solver.setTimeStepper( timeStepper );
       if( ! solver.setup( parameters ) )
          return false;
    }
@@ -419,27 +452,16 @@ bool SolverStarter< ConfigTag > :: runPDESolver( Problem& problem,
    Logger logger( logWidth, logFile );
    solver.writeProlog( logger, parameters  );
 
-   /****
-    * Set-up solver monitor and launch the main loop.
-    */
-   typedef IterativeSolverMonitor< typename Problem::RealType, typename Problem::IndexType > SolverMonitorType;
-   SolverMonitorType _tempSolverMonitor;
-   SolverMonitorType* solverMonitorPointer = &_tempSolverMonitor;
-   if( problem.getSolverMonitor() )
-      solverMonitorPointer = ( SolverMonitorType* ) problem.getSolverMonitor();
-
-   timeStepper.setSolverMonitor( *solverMonitorPointer );
-   solverMonitorPointer->setVerbose( verbose );
-   solverMonitorPointer->setTimer( this->totalTimer );
-   SolverMonitorThread t( *solverMonitorPointer );
-
    /****
     * Set-up timers
     */
    this->computeTimer.reset();
    this->ioTimer.reset();
-   solver.setComputeTimer( this->computeTimer );
-   solver.setIoTimer( this->ioTimer );
+   
+   /****
+    * Create solver monitor thread
+    */
+   SolverMonitorThread t( solver.getSolverMonitor() );
 
    /****
     * Start the solver
@@ -466,7 +488,6 @@ bool SolverStarter< ConfigTag > :: runPDESolver( Problem& problem,
       returnCode = solver.solve();
    }
 
-   solverMonitorPointer->stopMainLoop();
    if( ! returnCode ) {
       if( verbose )
          std::cerr << std::endl << "The solver did not converge. " << std::endl;
diff --git a/src/TNL/Timer.cpp b/src/TNL/Timer.cpp
index dc561202e7f6f7b759cc599edba7e37d2a1b7362..2d11949edc05706fd5d6200faf84e7f263bed921 100644
--- a/src/TNL/Timer.cpp
+++ b/src/TNL/Timer.cpp
@@ -46,99 +46,71 @@ void Timer::stop()
 
    if( ! this->stopState )
    {
-      /****
-       * Real time
-       */
-#ifdef HAVE_TIME
-      struct timeval tp;
-      int rtn = gettimeofday( &tp, NULL );
-      this->totalRealTime += ( double ) tp. tv_sec + 1.0e-6 * ( double ) tp. tv_usec - this->initialRealTime;
-#endif
- 
-      /****
-       * CPU time
-       */
-#ifdef HAVE_SYS_RESOURCE_H
-      rusage initUsage;
-      getrusage(  RUSAGE_SELF, &initUsage );
-      this->totalCPUTime += initUsage. ru_utime. tv_sec + 1.0e-6 * ( double ) initUsage. ru_utime. tv_usec - this->initialCPUTime;
-#endif
- 
-      /****
-       * CPU cycles
-       */
-      this->totalCPUCycles += this->rdtsc() - this->initialCPUCycles;
+      this->totalRealTime += this->readRealTime() - this->initialRealTime;
+      this->totalCPUTime += this->readCPUTime() - this->initialCPUTime;
+      this->totalCPUCycles += this->readCPUCycles() - this->initialCPUCycles;
       this->stopState = true;
    }
 }
 
 void Timer::start()
 {
-   /****
-    * Real time
-    */
-#ifdef HAVE_TIME
-   struct timeval tp;
-   int rtn = gettimeofday( &tp, NULL );
-   this->initialRealTime = ( double ) tp. tv_sec + 1.0e-6 * ( double ) tp. tv_usec;
-#endif
-
-   /****
-    * CPU Time
-    */
-#ifdef HAVE_SYS_RESOURCE_H
-   rusage initUsage;
-   getrusage( RUSAGE_SELF, &initUsage );
-   this->initialCPUTime = initUsage. ru_utime. tv_sec + 1.0e-6 * ( double ) initUsage. ru_utime. tv_usec;
-#endif
- 
-   /****
-    * CPU cycles
-    */
-   this->initialCPUCycles = this->rdtsc();
- 
+   this->initialRealTime = this->readRealTime();
+   this->initialCPUTime = this->readCPUTime();
+   this->initialCPUCycles = this->readCPUCycles(); 
    this->stopState = false;
 }
 
-double Timer::getRealTime()
+double Timer::getRealTime() const
 {
-#ifdef HAVE_TIME
    if( ! this->stopState )
-   {
-      this->stop();
-      this->start();
-   }
+    return this->readRealTime() - this->initialRealTime;
    return this->totalRealTime;
+}
+
+double Timer::getCPUTime() const
+{
+   if( ! this->stopState )
+    return this->readCPUTime() - this->initialCPUTime;
+   return this->totalCPUTime;
+}
+
+unsigned long long int Timer::getCPUCycles() const
+{
+   if( ! this->stopState )
+    return this->readCPUCycles() - this->initialCPUCycles;
+   return this->totalCPUCycles;
+}
+
+double Timer::readRealTime() const
+{
+#ifdef HAVE_TIME
+   struct timeval tp;
+   int rtn = gettimeofday( &tp, NULL );
+   return ( double ) tp. tv_sec + 1.0e-6 * ( double ) tp. tv_usec;
 #else
    return -1;
 #endif
 }
 
-double Timer::getCPUTime()
+double Timer::readCPUTime() const
 {
 #ifdef HAVE_SYS_RESOURCE_H
-   if( ! this->stopState )
-   {
-      this->stop();
-      this->start();
-   }
-   return this->totalCPUTime;
+   rusage initUsage;
+   getrusage( RUSAGE_SELF, &initUsage );
+   return initUsage. ru_utime. tv_sec + 1.0e-6 * ( double ) initUsage. ru_utime. tv_usec;
 #else
    return -1;
 #endif
 }
 
-unsigned long long int Timer::getCPUCycles()
+unsigned long long int Timer::readCPUCycles() const
 {
-   if( ! this->stopState )
-   {
-      this->stop();
-      this->start();
-   }
-   return this->totalCPUCycles;
+   return this->rdtsc();
 }
 
-bool Timer::writeLog( Logger& logger, int logLevel )
+
+bool Timer::writeLog( Logger& logger, int logLevel ) const
 {
    logger.writeParameter< double                 >( "Real time:",  this->getRealTime(),  logLevel );
    logger.writeParameter< double                 >( "CPU time:",   this->getCPUTime(),   logLevel );
diff --git a/src/TNL/Timer.h b/src/TNL/Timer.h
index 5019ab46f3bf2ff0dcef2bed65d96ac9af848457..f885bdd443bbeeec0af8d496a7daf9c2193daa27 100644
--- a/src/TNL/Timer.h
+++ b/src/TNL/Timer.h
@@ -27,15 +27,22 @@ class Timer
 
       void start();
 
-      double getRealTime();
+      double getRealTime() const;
 
-      double getCPUTime();
+      double getCPUTime() const;
 
-      unsigned long long int getCPUCycles();
+      unsigned long long int getCPUCycles() const;
  
-      bool writeLog( Logger& logger, int logLevel = 0 );
+      bool writeLog( Logger& logger, int logLevel = 0 ) const;
  
    protected:
+      
+      double readRealTime() const;
+
+      double readCPUTime() const;
+
+      unsigned long long int readCPUCycles() const;
+      
 
    double initialRealTime, totalRealTime,
           initialCPUTime, totalCPUTime;
@@ -44,7 +51,7 @@ class Timer
  
    bool stopState;
  
-   inline unsigned long long rdtsc()
+   inline unsigned long long rdtsc() const
    {
      unsigned hi, lo;
      __asm__ __volatile__ ("rdtsc" : "=a"(lo), "=d"(hi));
diff --git a/src/TNL/TypeInfo.h b/src/TNL/TypeInfo.h
new file mode 100644
index 0000000000000000000000000000000000000000..514fbacadebeea27e18dca4669db499b4e3430c7
--- /dev/null
+++ b/src/TNL/TypeInfo.h
@@ -0,0 +1,46 @@
+/***************************************************************************
+                          TypeInfo.h  -  description
+                             -------------------
+    begin                : Nov 10, 2017
+    copyright            : (C) 2017 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#pragma once
+
+#include <cfloat>
+#include <TNL/Devices/Cuda.h>
+
+namespace TNL 
+{
+
+template< typename Type >
+class TypeInfo
+{
+};
+
+template<>
+class TypeInfo< float >
+{
+   public:
+      
+      static __cuda_callable__ float getMaxValue() { return FLT_MAX; }
+      
+      static __cuda_callable__ float getMinValue() { return FLT_MIN; }   
+};
+
+template<>
+class TypeInfo< double >
+{
+   public:
+      
+      static __cuda_callable__ double getMaxValue() { return DBL_MAX; }
+      
+      static __cuda_callable__ double getMinValue() { return DBL_MIN; }
+   
+   
+};
+
+} // namespace TNL
\ No newline at end of file