diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 157374d4c7c1224768332810c40103f6aea4646c..9670a1ed73767ac132db7268e6169e9da17038fe 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -60,12 +60,11 @@ stages:
                 -DWITH_EXAMPLES=${WITH_EXAMPLES}
                 -DWITH_TOOLS=${WITH_TOOLS}
                 -DWITH_PYTHON=${WITH_PYTHON}
-#        - make
-#        - make test
+        # "install" implies the "all" target
 #        - make install
-        - ninja ${NINJAFLAGS}
+#        - make test
+        - ninja ${NINJAFLAGS} install
         - ninja test
-        - ninja install
         - popd
     variables:
         <<: *default_cmake_flags
@@ -73,15 +72,23 @@ stages:
 
 # Cuda builds are specified first because they take more time than host-only builds,
 # which can be allocated on hosts whitout GPUs.
+# Similarly, release builds are launched first to avoid the tail effect (they take
+# significantly more time than debug builds).
 
-cuda_base_Debug:
+cuda_full_Release:
     <<: *build_template
     tags:
+        - openmp
         - gpu
     variables:
         <<: *default_cmake_flags
+        WITH_OPENMP: "yes"
         WITH_CUDA: "yes"
-        BUILD_TYPE: Debug
+        BUILD_TYPE: Release
+        WITH_BENCHMARKS: "yes"
+        WITH_EXAMPLES: "yes"
+        WITH_TOOLS: "yes"
+        WITH_PYTHON: "yes"
 
 cuda_base_Release:
     <<: *build_template
@@ -92,19 +99,6 @@ cuda_base_Release:
         WITH_CUDA: "yes"
         BUILD_TYPE: Release
 
-cuda_mpi_Debug:
-    <<: *build_template
-    tags:
-        - openmp
-        - gpu
-        - mpi
-    variables:
-        <<: *default_cmake_flags
-        WITH_OPENMP: "yes"
-        WITH_CUDA: "yes"
-        WITH_MPI: "yes"
-        BUILD_TYPE: Debug
-
 cuda_mpi_Release:
     <<: *build_template
     tags:
@@ -133,20 +127,28 @@ cuda_full_Debug:
         WITH_TOOLS: "yes"
         WITH_PYTHON: "yes"
 
-cuda_full_Release:
+cuda_base_Debug:
+    <<: *build_template
+    tags:
+        - gpu
+    variables:
+        <<: *default_cmake_flags
+        WITH_CUDA: "yes"
+        BUILD_TYPE: Debug
+
+cuda_mpi_Debug:
     <<: *build_template
     tags:
         - openmp
         - gpu
+        - mpi
     variables:
         <<: *default_cmake_flags
         WITH_OPENMP: "yes"
         WITH_CUDA: "yes"
-        BUILD_TYPE: Release
-        WITH_BENCHMARKS: "yes"
-        WITH_EXAMPLES: "yes"
-        WITH_TOOLS: "yes"
-        WITH_PYTHON: "yes"
+        WITH_MPI: "yes"
+        BUILD_TYPE: Debug
+
 
 default_base_Debug:
     <<: *build_template
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
new file mode 100644
index 0000000000000000000000000000000000000000..21c479a2cdf65ff588bbc82077baf11c1825cb22
--- /dev/null
+++ b/CONTRIBUTING.md
@@ -0,0 +1,79 @@
+## How to configure git
+
+It is important to [configure your git username and email address](
+https://mmg-gitlab.fjfi.cvut.cz/gitlab/help/gitlab-basics/start-using-git.md#add-your-git-username-and-set-your-email),
+since every git commit will use this information to identify you as the author:
+
+    git config --global user.name "John Doe"
+    git config --global user.email "john.doe@example.com"
+
+In the TNL project, this username should be your real name (given name and
+family name) and the email address should be the email address used in the
+Gitlab profile. You should use the same configuration on all computers where
+you make commits. If you have made some commits with a different email address
+in the past, you can also add secondary email addresses to the Gitlab profile.
+
+## How to write good commit messages
+
+Begin with a short summary line a.k.a. message subject:
+
+- Use up to 50 characters; this is the git official preference.
+- Finish without a sentence-ending period.
+
+Continue with a longer description a.k.a. message body:
+
+- Add a blank line after the summary line, then write as much as you want.
+- Use up to 72 characters per line for typical text for word wrap.
+- Use as many characters as needed for atypical text, such as URLs, terminal
+  output, formatted messages, etc.
+- Include any kind of notes, links, examples, etc. as you want.
+
+See [5 Useful Tips For A Better Commit Message](
+https://robots.thoughtbot.com/5-useful-tips-for-a-better-commit-message) and
+[How to Write a Git Commit Message](https://chris.beams.io/posts/git-commit/)
+for examples and reasoning.
+
+## How to split commits
+
+Extensive changes should be split into multiple commits so that only related
+changes are made in each commit. All changes made in a commit should be
+described in its commit message. If describing all changes would not result in
+a good commit message, you should probably make multiple separate commits.
+
+Multiple small commits are better than one big commit, because later they can be
+easily [squashed](https://git-scm.com/book/en/v2/Git-Tools-Rewriting-History)
+together, whereas splitting a big commit into logical parts is significantly
+more difficult.
+
+> __Tip:__ Use `git add -p` to stage only part of your working tree for the next
+> commit. See [git add -p: The most powerful git feature you're not using yet](
+https://johnkary.net/blog/git-add-p-the-most-powerful-git-feature-youre-not-using-yet/).
+
+## Rebase-based workflow
+
+The development of new features should follow the rebase-based workflow:
+
+- Create a new _feature_ branch based on the main branch (`develop`).
+- Make commits with your work in the feature branch.
+- When there are other new commits in the `develop` branch, you should do a
+  [rebase](https://git-scm.com/book/en/v2/Git-Branching-Rebasing) to rewind your
+  commits on top of the current develop branch.
+    - If there are conflicts, you will need to resolve them manually. Hence, it
+      is a good practice to rebase as often as possible (generally as soon as
+      the changes appear in the `develop` branch) and to split commits into
+      logical parts.
+- When your work is ready for review, you can open a [merge request](
+  https://mmg-gitlab.fjfi.cvut.cz/gitlab/tnl/tnl-dev/merge_requests/new).
+    - If the branch is not ready for merge yet, prepend `[WIP]` to the merge
+      request title to indicate _work in progress_ and to prevent a premature
+      merge.
+    - This is also a good time to squash small commits (e.g. typos, forgotten
+      changes or trivial corrections) with relevant bigger commits to make the
+      review easier.
+- When reviewed, the feature branch can be merged into the develop branch.
+
+The main advantages of this workflow are linear history, clear commits and
+reduction of merge conflicts. See [A rebase-based workflow](
+https://brokenco.de/2010/04/02/a-rebase-based-workflow.html) and
+[Why is rebase-then-merge better than just merge](https://stackoverflow.com/a/457988)
+([complement](https://stackoverflow.com/a/804178)) for reference.
\ No newline at end of file
diff --git a/src/TNL/CMakeLists.txt b/src/TNL/CMakeLists.txt
index 306bd82a3c5c6633c893beafa90aa768f0e83bee..5d2bf5b15bb9999a70bed7efcb3718bad92c5322 100644
--- a/src/TNL/CMakeLists.txt
+++ b/src/TNL/CMakeLists.txt
@@ -1,5 +1,6 @@
 ADD_SUBDIRECTORY( Config )
 ADD_SUBDIRECTORY( Containers )
+ADD_SUBDIRECTORY( DistributedContainers )
 ADD_SUBDIRECTORY( Communicators )
 ADD_SUBDIRECTORY( Debugging )
 ADD_SUBDIRECTORY( Devices )
diff --git a/src/TNL/Containers/List_impl.h b/src/TNL/Containers/List_impl.h
index 136f4cc986a35063a8be9f52ab4fa4e705b6c273..e67be136cd19e28db2f3f2da1183e8b2d7bf3236 100644
--- a/src/TNL/Containers/List_impl.h
+++ b/src/TNL/Containers/List_impl.h
@@ -38,7 +38,7 @@ List< T >::~List()
 template< typename T >
 String List< T >::getType()
 {
-   return String( "List< " ) + TNL::getType< T >() +  String( " >" );
+   return String( "Containers::List< " ) + TNL::getType< T >() +  String( " >" );
 }
 
 template< typename T >
diff --git a/src/TNL/Devices/Cuda.cpp b/src/TNL/Devices/Cuda.cpp
index 1bc85e3c829ff4a2f44518f4016c7ef90a447b27..9dc2ff9d00b035154195d4ef81dc39a94a36345a 100644
--- a/src/TNL/Devices/Cuda.cpp
+++ b/src/TNL/Devices/Cuda.cpp
@@ -23,7 +23,7 @@ Timer Cuda::smartPointersSynchronizationTimer;
 
 String Cuda::getDeviceType()
 {
-   return String( "Cuda" );
+   return String( "Devices::Cuda" );
 }
 
 int Cuda::getNumberOfBlocks( const int threads,
diff --git a/src/TNL/Devices/MIC.h b/src/TNL/Devices/MIC.h
index 776b7c36fe7a44fd0063142708fd0789f350c45c..5b6a9a4f8c40adf1b26bfdc687254db1d1241872 100644
--- a/src/TNL/Devices/MIC.h
+++ b/src/TNL/Devices/MIC.h
@@ -71,7 +71,7 @@ class MIC
    
         static String getDeviceType()
         {
-            return String( "MIC" );
+            return String( "Devices::MIC" );
         };
         
 #ifdef HAVE_MIC  
diff --git a/src/TNL/DistributedContainers/CMakeLists.txt b/src/TNL/DistributedContainers/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..7aa8b2cde74b43c23641157e906c6632260c405e
--- /dev/null
+++ b/src/TNL/DistributedContainers/CMakeLists.txt
@@ -0,0 +1,16 @@
+SET( headers DistributedArray.h
+             DistributedArray_impl.h
+             DistributedArrayView.h
+             DistributedArrayView_impl.h
+             DistributedMatrix.h
+             DistributedMatrix_impl.h
+             DistributedSpMV.h
+             DistributedVector.h
+             DistributedVector_impl.h
+             DistributedVectorView.h
+             DistributedVectorView_impl.h
+             Partitioner.h
+             Subrange.h
+    )
+
+INSTALL( FILES ${headers} DESTINATION ${TNL_TARGET_INCLUDE_DIRECTORY}/DistributedContainers )
diff --git a/src/TNL/Solvers/Linear/CMakeLists.txt b/src/TNL/Solvers/Linear/CMakeLists.txt
index 2e33af0b8e081c2374990cbe6ff1eb70a0a897d0..f5815d34f7f161d1082779eadffe565249a6a9ed 100644
--- a/src/TNL/Solvers/Linear/CMakeLists.txt
+++ b/src/TNL/Solvers/Linear/CMakeLists.txt
@@ -16,6 +16,7 @@ SET( headers BICGStab.h
              SOR_impl.h
              TFQMR.h
              TFQMR_impl.h
+             Traits.h
              UmfpackWrapper.h
              UmfpackWrapper_impl.h
    )
diff --git a/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h b/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h
index 2efc01001b9bb5438f2c6613409065c03bf79164..8d294bd8ae73e4a3e2e41cc7e4ede262a80740da 100644
--- a/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h
+++ b/src/TNL/Solvers/Linear/Preconditioners/Preconditioner.h
@@ -15,8 +15,7 @@
 #include <TNL/Containers/VectorView.h>
 #include <TNL/Pointers/SharedPointer.h>
 #include <TNL/Config/ParameterContainer.h>
-
-#include "../Traits.h"
+#include <TNL/Solvers/Linear/Traits.h>
 
 namespace TNL {
 namespace Solvers {
diff --git a/src/Tools/CMakeLists.txt b/src/Tools/CMakeLists.txt
index 79aa6243008c41f37856c1a0762e45098be52b6e..3dc914c9d0ed8e3feac4c463cd3ab23e05a32833 100644
--- a/src/Tools/CMakeLists.txt
+++ b/src/Tools/CMakeLists.txt
@@ -31,9 +31,6 @@ target_link_libraries (tnl-dicom-reader tnl ${DCMTK_LIBRARIES} )
 ADD_EXECUTABLE(tnl-lattice-init tnl-lattice-init.cpp )
 target_link_libraries (tnl-lattice-init tnl )
 
-ADD_EXECUTABLE( tnl-functions-benchmark functions-benchmark.cpp )
-target_link_libraries( tnl-functions-benchmark tnl )
-
 IF( BUILD_CUDA )
    CUDA_ADD_EXECUTABLE( tnl-cuda-arch tnl-cuda-arch.cu )
    INSTALL( TARGETS tnl-cuda-arch
@@ -55,9 +52,7 @@ INSTALL( TARGETS tnl-init
 INSTALL( FILES ${PROJECT_TOOLS_PATH}/tnl-bindir
                ${PROJECT_TOOLS_PATH}/tnl-compile
                ${PROJECT_TOOLS_PATH}/tnl-link
-               tnl-time-series2png
                tnl-err2eoc
-               tnl-eoc-test-log
                tnl-log-to-html.py
          DESTINATION bin
          PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE )
diff --git a/src/Tools/functions-benchmark.cpp b/src/Tools/functions-benchmark.cpp
deleted file mode 100644
index a47dcced31b33c61052a68331d89e9df746fe552..0000000000000000000000000000000000000000
--- a/src/Tools/functions-benchmark.cpp
+++ /dev/null
@@ -1,38 +0,0 @@
-/***************************************************************************
-                          functions-benchmark.cpp  -  description
-                             -------------------
-    begin                : Jul 4, 2010
-    copyright            : (C) 2010 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/* See Copyright Notice in tnl/Copyright */
-
-#include "functions-benchmark.h"
-
-int main( int argc, char* argv[] )
-{
-   const long int loops = 1 << 24;
-
-   std::cout << "Runnning benchmarks in single precision on CPU ... " << std::endl;
-   benchmarkAddition< float >( loops );
-   benchmarkMultiplication< float >( loops );
-   benchmarkDivision< float >( loops );
-   benchmarkSqrt< float >( loops );
-   benchmarkSin< float >( loops );
-   benchmarkExp< float >( loops );
-   benchmarkPow< float >( loops );
-
-   std::cout << "Runnning benchmarks in double precision on CPU ... " << std::endl;
-   benchmarkAddition< double >( loops );
-   benchmarkMultiplication< float >( loops );
-   benchmarkDivision< double >( loops );
-   benchmarkSqrt< double >( loops );
-   benchmarkSin< double >( loops );
-   benchmarkExp< double >( loops );
-   benchmarkPow< double >( loops );
-
-
-
-   return EXIT_SUCCESS;
-}
diff --git a/src/Tools/functions-benchmark.h b/src/Tools/functions-benchmark.h
deleted file mode 100644
index 2f04f6fbeb8577e9af131c76a0f6bf8c8b100b68..0000000000000000000000000000000000000000
--- a/src/Tools/functions-benchmark.h
+++ /dev/null
@@ -1,183 +0,0 @@
-/***************************************************************************
-                          functions-benchmark.h  -  description
-                             -------------------
-    begin                : Jul 4, 2010
-    copyright            : (C) 2010 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/* See Copyright Notice in tnl/Copyright */
-
-#ifndef FUNCTIONSBENCHMARK_H_
-#define FUNCTIONSBENCHMARK_H_
-
-#include <iostream>
-#include <math.h>
-
-#include <TNL/Timer.h>
-
-using namespace TNL;
-
-template< typename REAL > void benchmarkAddition( long int loops )
-{
-   std::cout << "Benchmarking addition on CPU ( " << loops << " loops ) ... " << std::flush;
-   Timer timer;
-   timer.start();
-
-   REAL a1 = 1.2;
-   REAL a2 = 1.2;
-   REAL a3 = 1.2;
-   REAL a4 = 1.2;
-   for( long int i = 0; i < loops; i ++ )
-   {
-      a1 += REAL( 0.1 );
-      a2 += REAL( 0.1 );
-      a3 += REAL( 0.1 );
-      a4 += REAL( 0.1 );
-   }
-
-   double cpu_time = timer.getCPUTime();
-   std::cout << " ( " << a1 + a2 + a3 + a4 << " ) " <<  cpu_time << "secs. " << 4.0 * ( ( double ) loops ) / cpu_time * 1.0e-9 << " GFLOPS." << std::endl;
-}
-
-template< typename REAL > void benchmarkMultiplication( const long int loops )
-{
-   std::cout << "Benchmarking multiplication on CPU ( " << loops << " loops ) ... " << std::flush;
-   Timer timer;
-   timer.start();
-
-   REAL a1 = 1.0e9;
-   REAL a2 = 1.0e9;
-   REAL a3 = 1.0e9;
-   REAL a4 = 1.0e9;
-   for( long int i = 0; i < loops; i ++ )
-   {
-      {
-         a1 *= REAL( 0.99 );
-         a2 *= REAL( 0.99 );
-         a3 *= REAL( 0.99 );
-         a4 *= REAL( 0.99 );
-         if( a1 < REAL( 0.01 ) ) a1 = a2 = a3 = a4 = REAL( 1.0e9 );
-      }
-   }
-
-   double cpu_time = timer.getCPUTime();
-   std::cout << " ( " << a1 * a2 * a3 * a4 << " ) " <<  cpu_time << "secs. " << 4.0 * ( ( double ) loops ) / cpu_time * 1.0e-9 << " GFLOPS." << std::endl;
-}
-
-template< typename REAL > void benchmarkDivision( long int loops )
-{
-   std::cout << "Benchmarking division on CPU ( " << loops << " loops ) ... " << std::flush;
-   Timer timer;
-   timer.start();
-
-   REAL a1( 1.0e9 );
-   REAL a2( 1.0e9 );
-   REAL a3( 1.0e9 );
-   REAL a4( 1.0e9 );
-   for( long int i = 0; i < loops; i ++ )
-   {
-      a1 /= REAL( 1.1 );
-      a2 /= REAL( 1.1 );
-      a3 /= REAL( 1.1 );
-      a4 /= REAL( 1.1 );
-      if( a1 < REAL( 0.01 ) ) a1 = a2 = a3 = a4 = REAL( 1.0e9 );
-   }
-
-   double cpu_time = timer.getCPUTime();
-   std::cout << " ( " << a1 / a2 / a3 / a4 << " ) " << cpu_time << "secs. " << 4.0 * ( ( double ) loops / 2 ) / cpu_time * 1.0e-9 << " GFLOPS." << std::endl;
-}
-
-template< typename REAL > void benchmarkSqrt( long int loops )
-{
-   std::cout << "Benchmarking sqrt on CPU ( " << loops << " loops ) ... " << std::flush;
-   Timer timer;
-   timer.start();
-
-   REAL a1( 1.0e9 );
-   REAL a2( 1.0e9 );
-   REAL a3( 1.0e9 );
-   REAL a4( 1.0e9 );
-   for( long int i = 0; i < loops; i ++ )
-   {
-      a1 = ::sqrt( a1 );
-      a2 = ::sqrt( a2 );
-      a3 = ::sqrt( a3 );
-      a4 = ::sqrt( a4 );
-      if( a1 < REAL( 100.0 ) ) a1 = a2 = a3 = a4 = REAL( 1.0e9 );
-   }
-
-   double cpu_time = timer.getCPUTime();
-   std::cout << " ( " << a1 + a2 + a3 + a4 << " ) " << cpu_time << "secs. " << 4.0 * ( ( double ) loops / 2 ) / cpu_time * 1.0e-9 << " GFLOPS." << std::endl;
-}
-
-template< typename REAL > void benchmarkSin( long int loops )
-{
-   std::cout << "Benchmarking sin on CPU ( " << loops << " loops ) ... " << std::flush;
-   Timer timer;
-   timer.start();
-
-   REAL a1( 1.0e9 );
-   REAL a2( 1.0e9 );
-   REAL a3( 1.0e9 );
-   REAL a4( 1.0e9 );
-   for( long int i = 0; i < loops; i ++ )
-   {
-      a1 = ::sin( a1 );
-      a2 = ::sin( a2 );
-      a3 = ::sin( a3 );
-      a4 = ::sin( a4 );
-   }
-
-   double cpu_time = timer.getCPUTime();
-   std::cout << " ( " << a1 + a2 + a3 + a4 << " ) " << cpu_time << "secs. " << 4.0 * ( ( double ) loops ) / cpu_time * 1.0e-9 << " GFLOPS." << std::endl;
-}
-
-template< typename REAL > void benchmarkExp( long int loops )
-{
-   std::cout << "Benchmarking exp on CPU ( " << loops << " loops ) ... " << std::flush;
-   Timer timer;
-   timer.start();
-
-   REAL a1( 1.1 );
-   REAL a2( 1.1 );
-   REAL a3( 1.1 );
-   REAL a4( 1.1 );
-   for( long int i = 0; i < loops; i ++ )
-   {
-      a1 = exp( a1 );
-      a2 = exp( a2 );
-      a3 = exp( a3 );
-      a4 = exp( a4 );
-      if( a1 > REAL( 1.0e9 ) ) a1 = a2 = a3 = a4 = REAL( 1.1 );
-   }
-
-   double cpu_time = timer.getCPUTime();
-   std::cout << " ( " << a1 + a2 + a3 + a4 << " ) " << cpu_time << "secs. " << 4.0 * ( ( double ) loops) / cpu_time * 1.0e-9 << " GFLOPS." << std::endl;
-}
-
-template< typename REAL > void benchmarkPow( long int loops )
-{
-   std::cout << "Benchmarking pow on CPU ( " << loops << " loops ) ... " << std::flush;
-   Timer timer;
-   timer.start();
-
-   REAL a1( 1.0e9 );
-   REAL a2( 1.0e9 );
-   REAL a3( 1.0e9 );
-   REAL a4( 1.0e9 );
-   for( long int i = 0; i < loops; i ++ )
-   {
-      a1 = ::pow( a1, REAL( 0.9 ) );
-      a2 = ::pow( a2, REAL( 0.9 ) );
-      a3 = ::pow( a3, REAL( 0.9 ) );
-      a4 = ::pow( a4, REAL( 0.9 ) );
-      if( a1 < REAL( 1.0 ) ) a1 = a2 = a3 = a4 = REAL( 1.0e9 );
-   }
-
-   double cpu_time = timer.getCPUTime();
-   std::cout << " ( " << a1 + a2 + a3 + a4 << " ) " << cpu_time << "secs. " << 4.0 * ( ( double ) loops) / cpu_time * 1.0e-9 << " GFLOPS." << std::endl;
-}
-
-
-#endif /* FUNCTIONSBENCHMARK_H_ */
diff --git a/src/Tools/polygonizer.cpp b/src/Tools/polygonizer.cpp
deleted file mode 100644
index c80a72b286c7048c5b13a507d108bf4217c07022..0000000000000000000000000000000000000000
--- a/src/Tools/polygonizer.cpp
+++ /dev/null
@@ -1,489 +0,0 @@
-/***************************************************************************
-                          polygonizer.cpp  -  description
-                             -------------------
-    begin                : Mon Feb 11 2002
-    copyright            : (C) 2002 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/* See Copyright Notice in tnl/Copyright */
-
-#include <string.h>
-#include "polygonizer.h"
-#include "../solids/reciever.h"
-#include "../solids/solid.h"
-#include "hash.h"
-
-POLYGONIZER polygonizer;
-
-static int corner1[ 12 ] =
-   { LBN, LTN, LBN, LBF, RBN, RTN, RBN, RBF, LBN, LBF, LTN, LTF };
-
-static int corner2[ 12 ] =
-   { LBF, LTF, LTN, LTF, RBF, RTF, RTN, RTF, RBN, RBF, RTN, RTF };
-
-// these are fields of corners for edges in order
-//   LB,  LT,  LN,  LF,  RB,  RT,  RN,  RF,  BN,  BF,  TN,  TF
-
-static int left_face[ 12 ] =
-   { B, L, L, F, R, T, N, R, N, B, T, F };
-// face on left when going corner1 to corner2
-
-static int right_face[ 12 ] =
-   { L, T, N, L, B, R, R, F, B, F, N, T };
-// face on right when going coner1 to corner2
-
-//---------------------------------------------------------------------------
-CUBE :: CUBE( int _i, int _j, int _k )
-      : i( _i ), j( _j ), k( _k )
-{
-   bzero( corners, 8 * sizeof( CORNER* ) );
-}
-//---------------------------------------------------------------------------
-POLYGONIZER :: POLYGONIZER()
-{
-   edge_hash = new ( list< EDGE* >* )[ 2 * hash_size ];
-   bzero( edge_hash, 2 * hash_size * sizeof( list< EDGE* >* ) );
-   corner_hash = new ( list< CORNER* >* )[ 2 * hash_size ];
-   bzero( corner_hash, 2 * hash_size * sizeof( list< CORNER* >* ) );
-   center_hash = new ( list< CENTER* >* )[ 2 * hash_size ];
-   bzero( center_hash, 2 * hash_size * sizeof( list< CENTER* >* ) );
-   Make_Cube_Table();
-}
-//---------------------------------------------------------------------------
-POLYGONIZER :: ~POLYGONIZER()
-{
-   assert( edge_hash );
-   delete[] edge_hash;
-   assert( corner_hash );
-   delete[] corner_hash;
-   assert( center_hash );
-   delete[] center_hash;
-}
-//---------------------------------------------------------------------------
-int POLYGONIZER :: Parametric( const SOLID* sld, RECIEVER* reciever,
-                               const unsigned int slices, const unsigned int stacks)
-{
-   assert( slices && stacks );
-   const double x_step = 1.0 / slices;
-   const double y_step = 1.0 / stacks;
-   for( unsigned int i = 0; i < stacks; i ++ )
-      for( unsigned int j = 0; j < slices; j ++ )
-      {
-         double x = j * x_step;
-         double y = i * y_step;
-         VECTOR n1, n2, n3, n4;
-         VECTOR u1 = sld -> Surface_Point( x, y, &n1 );
-         VECTOR u2 = sld -> Surface_Point( x + x_step, y, &n2 );
-         VECTOR u3 = sld -> Surface_Point( x + x_step, y + y_step, &n3 );
-         VECTOR u4 = sld -> Surface_Point( x, y + y_step, &n4 );
-
-         VERTEX v1( u1, n1 );
-         VERTEX v2( u2, n2 );
-         VERTEX v3( u3, n3 );
-         VERTEX v4( u4, n4 );
-
-         /*TRIANGLE* t1 = new TRIANGLE( new VERTEX( v1 ),
-                                      new VERTEX( v4 ),
-                                      new VERTEX( v2 ) ) ;
-         TRIANGLE* t2 = new TRIANGLE( new VERTEX( v3 ),
-                                      new VERTEX( v2 ),
-                                      new VERTEX( v4 ) ) ;*/
-
-          if( ! reciever -> Insert_Triangle( v1, v4, v2 ) ) return 0;
-          if( ! reciever -> Insert_Triangle( v3, v2, v4 ) ) return 0;
-      }
-   return 1;
-}
-//---------------------------------------------------------------------------
-int POLYGONIZER :: Init( const VECTOR& pos, const SOLID* sld, const double& cb_sz )
-{
-   cube_size = cb_sz;
-   solid = sld;
-   VECTOR in( pos ), out( pos );
-   if( ! Find_Point( 1, in, 1.0, solid ) ||
-       ! Find_Point( 0, out, 1.0, solid ) ) return 0;
-   CUBE *c = NULL;
-   if( cube_stack. empty() ) // this is the initial cube
-   {
-      position = Compute_Surface_Point( in, out, solid -> Continuous_Function( in ) );
-      c = new CUBE( 0, 0, 0 );
-   }
-   else
-   {
-      VECTOR _pos = Compute_Surface_Point( in, out, solid -> Continuous_Function( in ) ) - position;
-      _pos *= 1.0 / cube_size;
-      int i( ( int ) _pos. x ),
-          j( ( int ) _pos. y ),
-          k( ( int ) _pos. z );
-      if( Set_Center( i, j, k ) ) return 1; //we have already found this cube
-      c = new CUBE( i, j, k );
-   }
-   for( int i = 0; i < 8; i ++ ) c -> corners[ i ] =
-      Set_Corner( ( i >> 2 ) & 1, ( i >> 1 ) & 1, i & 1 );
-
-   cube_stack. push( c );
-   return 1;
-}
-//---------------------------------------------------------------------------
-int POLYGONIZER :: Implicit( const SOLID* sld, RECIEVER* _reciever,
-                             const VECTOR& ps, const VECTOR& cr,
-                             const double& cb_size, int mode )
-{
-   solid = sld;
-   reciever = _reciever;
-   cube_size = cb_size;
-   bound_ps = ps;
-   bound_cr = cr;
-   if( cube_stack. empty() && ! Init( solid -> Position(), solid, cb_size ) )
-   {
-      std::cerr << "Can not find initial points for polygonization" << std::endl;
-      return 0;
-   }
-   std::cout << "Starting polygonizer ... " << std::endl;
-   /*VECTOR in( sld -> Position() ), out( sld -> Position() );
-   if( ! Find_Point( 1, in, 1.0 ) ||
-       ! Find_Point( 0, out, 1.0 ) )
-      {
-         std::cerr << "Can not find initial points for polygonization" << std::endl;
-         return 0;
-      }
-   std::cout << "Starting polygonizer ... " << std::endl;
-   position = Compute_Surface_Point( in, out, solid -> Continuous_Function( in ) );
-
-   CUBE *c = new CUBE( 0, 0, 0 );
-   for( int i = 0; i < 8; i ++ ) c -> corners[ i ] =
-      Set_Corner( ( i >> 2 ) & 1, ( i >> 1 ) & 1, i & 1 );
-
-   cube_stack. push( c );*/
-   int jkl;
-   CUBE* c;
-   while( ! cube_stack. empty() )
-   {
-      c = cube_stack. top();
-      std::cout << "Current stack size is " << cube_stack. size() << " cubes. " << '\r' << std::flush;
-      cube_stack. pop();
-      jkl = cube_stack. size();
-      if( ! Triangulate_Cube( c ) )
-      {
-         Free_Memory();
-         return 0;
-      }
-      Test_Face( c -> i - 1, c -> j, c -> k, c, L, LBN, LBF, LTN, LTF );
-      std::cout << "Current stack size is " << cube_stack. size() << " cubes. " << '\r' << std::flush;
-      Test_Face( c -> i + 1, c -> j, c -> k, c, R, RBN, RBF, RTN, RTF );
-      std::cout << "Current stack size is " << cube_stack. size() << " cubes. " << '\r' << std::flush;
-      Test_Face( c -> i, c -> j - 1, c -> k, c, B, LBN, LBF, RBN, RBF );
-      std::cout << "Current stack size is " << cube_stack. size() << " cubes. " << '\r' << std::flush;
-      Test_Face( c -> i, c -> j + 1, c -> k, c, T, LTN, LTF, RTN, RTF );
-      std::cout << "Current stack size is " << cube_stack. size() << " cubes. " << '\r' << std::flush;
-      Test_Face( c -> i, c -> j, c -> k - 1, c, N, LBN, LTN, RBN, RTN );
-      std::cout << "Current stack size is " << cube_stack. size() << " cubes. " << '\r' << std::flush;
-      Test_Face( c -> i, c -> j, c -> k + 1, c, F, LBF, LTF, RBF, RTF );
-      std::cout << "Current stack size is " << cube_stack. size() << " cubes. " << '\r' << std::flush;
-   }
-   std::cout << std::endl;
-   Free_Memory();
-   return 1;
-}
-//---------------------------------------------------------------------------
-int POLYGONIZER :: Triangulate_Cube( CUBE* cube )
-{
-   #ifdef DBG_POLYGONIZER
-   std::cout << "Cube polygonize: " << cube -> i << " " << cube -> j << " " << cube -> k << std::endl;;
-   #endif
-   int index = 0;
-   // this is a index of cube in cube_table, it will be count by
-   // the signs of function in corners
-   for( int i = 0; i < 8; i ++ )
-      if( cube -> corners[ i ] -> value > 0.0 ) index += ( 1 << i );
-   for( size_t i = 0; i < cube_table[ index ]. size(); i ++ )
-   {
-      vector< int >& edges_vector = ( cube_table[ index ] )[ i ];
-      VERTEX *v1( 0 ), *v2( 0 ), *v3( 0 );
-      for( size_t j = 0; j < edges_vector. size(); j ++ )
-      {
-         CORNER* c1 = cube -> corners[ corner1[ edges_vector[ j ] ] ];
-         CORNER* c2 = cube -> corners[ corner2[ edges_vector[ j ] ] ];
-         VERTEX* c = Get_Vertex( c1, c2 );
-         if( j == 0 ) v1 = c;
-         if( j == 1 ) v2 = c;
-         if( j > 1 )
-         {
-            v3 = c;
-            if( ! reciever -> Insert_Triangle( *v1, *v3, *v2 ) ) return 0;
-            v2 = v3;
-         }
-      }
-   }
-   return 1;
-}
-//---------------------------------------------------------------------------
-VERTEX* POLYGONIZER :: Get_Vertex( CORNER* c1, CORNER* c2 )
-{
-   VERTEX* vertex;
-   if( Check_Edge( c1, c2, vertex ) ) return vertex;
-   // the vertex has been already computed
-   VECTOR p = Compute_Surface_Point( c1 -> position, c2 -> position, c1 -> value );
-   vertex = new VERTEX( p, solid -> Normal( p ) );
-   Set_Edge( c1, c2, vertex );
-   return vertex;
-}
-//---------------------------------------------------------------------------
-VECTOR POLYGONIZER :: Compute_Surface_Point( const VECTOR& v1, const VECTOR& v2, const double& val )
-{
-   VECTOR pos, neg, p;
-   if( val < 0.0 )
-   {
-      pos = v2; neg = v1;
-   }
-   else
-   {
-      pos = v1; neg = v2;
-   }
-   int i = 0;
-   while( i ++ < max_iter )
-   {
-      p = 0.5 * ( pos + neg );
-      if( solid -> Continuous_Function( p ) > 0 ) pos = p;
-      else neg = p;
-   }
-   return p;
-}
-//---------------------------------------------------------------------------
-void POLYGONIZER :: Set_Edge( CORNER* c1, CORNER* c2, VERTEX* v )
-{
-   int& i1 = c1 -> i;
-   int& j1 = c1 -> j;
-   int& k1 = c1 -> k;
-   int& i2 = c2 -> i;
-   int& j2 = c2 -> j;
-   int& k2 = c2 -> k;
-   if( i1 > i2 || ( i1 == i2 && ( j1 > j2 || ( j1 == j2  && k1 > k2 ) ) ) )
-   {
-      CORNER* c = c1; c1 = c2; c2 = c;
-   }
-   unsigned int index = hash_index( i1, j1, k1 ) + hash_index( i2, j2, k2 );
-   if( ! edge_hash[ index ] ) edge_hash[ index ] = new list< EDGE* >;
-   edge_hash[ index ] -> push_back( new EDGE( c1, c2, v ) );
-}
-//---------------------------------------------------------------------------
-int POLYGONIZER :: Check_Edge( CORNER* c1, CORNER* c2, VERTEX*& vertex )
-{
-   int& i1 = c1 -> i;
-   int& j1 = c1 -> j;
-   int& k1 = c1 -> k;
-   int& i2 = c2 -> i;
-   int& j2 = c2 -> j;
-   int& k2 = c2 -> k;
-   if( i1 > i2 || ( i1 == i2 && ( j1 > j2 || ( j1 == j2  && k1 > k2 ) ) ) )
-   {
-      CORNER* c = c1; c1 = c2; c2 = c;
-   }
-   list< EDGE* >* el = edge_hash[ hash_index( i1, j1, k1 ) + hash_index( i2, j2, k2 ) ];
-   if( ! el ) return 0;
-   list< EDGE* > :: iterator it = el -> begin();
-   while( it != el -> end() )
-   {
-      EDGE* e = * it ++;
-      if( e -> corners[ 0 ] == c1 && e -> corners[ 1 ] == c2 )
-      {
-         vertex = e -> vertex;
-         return 1;
-      }
-   }
-   return 0;
-}
-//---------------------------------------------------------------------------
-CORNER* POLYGONIZER :: Set_Corner( int i, int j, int k )
-{
-   #ifdef DBG_POLYGONIZER
-      std::cout << "Set Corner index: " << i << " " << j << " " << k << std::endl;
-   #endif
-   int index = hash_index( i, j, k );
-   list< CORNER* >*& cl = corner_hash[ index ];
-   if( ! cl ) cl = new list< CORNER* >;
-   list< CORNER* > :: iterator it = cl -> begin();
-   CORNER* cr;
-   while( it != cl -> end() )
-   {
-      cr = * it ++;
-      if( ( cr -> i == i ) && ( cr -> j == j ) && ( cr -> k == k ) ) return cr;
-   }
-   VECTOR tmp = ( VECTOR( i, j, k ) - VECTOR( 0.5 ) );
-   VECTOR p = position + cube_size * tmp;
-   double val = solid -> Continuous_Function( p );
-   #ifdef DBG_POLYGONIZER
-   std::cout << "Set Corer: value " << val << std::endl;
-   #endif
-   cr = new CORNER( i, j, k, p, val );
-   cl -> push_back( cr );
-   return cr;
-}
-//---------------------------------------------------------------------------
-int POLYGONIZER :: Set_Center( int i, int j, int k )
-{
-   #ifdef DBG_POLYGONIZER
-   std::cout << "Set Center: " << i << " " << j << " " << k << std::endl;
-   #endif
-   list< CENTER* >*& cl = center_hash[ hash_index( i, j, k ) ];
-   if( ! cl ) cl = new list< CENTER* >;
-   list< CENTER* > :: iterator it = cl -> begin();
-   while( it != cl -> end() )
-   {
-      CENTER* cr = * it ++;
-      if( cr -> i == i && cr -> j == j && cr -> k == k ) return 1;
-   }
-   cl -> push_back( new CENTER( i, j, k ) );
-   return 0;
-}
-//---------------------------------------------------------------------------
-void POLYGONIZER :: Make_Cube_Table()
-{
-   int done_edges[ 12 ], positive[ 8 ];
-   // we have 256 different possibilities to sign all 8 corners of a cube
-   // they are describe here in cube_table
-   for( int i = 0; i < 256; i ++ )
-   {
-      for( int e = 0; e < 12; e ++ ) done_edges[ e ] = 0;
-      for( int c = 0; c < 8; c ++ ) positive[ c ] = ( i >> c ) & 1;
-      // each bit in 'i' represent one corner of cube and sign of function value
-      // '1' means positive value, '0' means negative value
-      for( int e = 0; e < 12; e ++ )
-      {
-         if( ! done_edges[ e ] &&
-             // we didn't process this edge yet ...
-             ( positive[ corner1[ e ] ] != positive[ corner2[ e ] ] ) )
-             // ... and corners of this edge have different sign
-            {
-               int start = e;
-               int edge = e;
-               int face = positive[ corner1[ e ] ] ? right_face[ e ] : left_face[ e ];
-               // get face that is to right of edge from positive to negatve corner
-               vector< int > temp_vector;
-               while( 1 )
-               {
-                  edge = Next_Clockwise_Edge( edge, face );
-                  done_edges[ edge ] = 1;
-                  if( positive[ corner1[ edge ] ] != positive[ corner2[ edge ] ] )
-                  {
-                     temp_vector. push_back( edge );
-                     if( edge == start ) break;
-                     face = Other_Face( edge, face );
-                  }
-               }
-               cube_table[ i ]. push_back( temp_vector );
-            }
-      }
-   }
-}
-//---------------------------------------------------------------------------
-int POLYGONIZER :: Next_Clockwise_Edge( int edge, int face )
-{
-   switch( edge )
-   {
-      case LB: return ( face == L ) ? LF : BN;
-      case LT: return ( face == L ) ? LN : TF;
-      case LN: return ( face == L ) ? LB : TN;
-      case LF: return ( face == L ) ? LT : BF;
-      case RB: return ( face == R ) ? RN : BF;
-      case RT: return ( face == R ) ? RF : TN;
-      case RN: return ( face == R ) ? RT : BN;
-      case RF: return ( face == R ) ? RB : TF;
-      case BN: return ( face == B ) ? RB : LN;
-      case BF: return ( face == B ) ? LB : RF;
-      case TN: return ( face == T ) ? LT : RN;
-      case TF: return ( face == T ) ? RT : LF;
-   }
-   return 0; // this is just for avoiding compiler warning
-}
-//---------------------------------------------------------------------------
-int POLYGONIZER :: Other_Face( int edge, int face )
-{
-   int other_face = left_face[ edge ];
-   return face == other_face ? right_face[ edge ] : other_face;
-}
-//---------------------------------------------------------------------------
-void POLYGONIZER :: Test_Face( int i, int j, int k, CUBE* old, int face,
-                               int c1, int c2, int c3, int c4 )
-{
-   static int face_bit[ 6 ] = { 2, 2, 1, 1, 0, 0 };
-   int bit = face_bit[ face ];
-   int pos = old -> corners[ c1 ] -> value > 0.0 ? 1 : 0;
-   // test id no surface crossing, cube out of bounds, or already visited
-   if( ( old -> corners[ c2 ] -> value > 0.0 ) == pos &&
-       ( old -> corners[ c3 ] -> value > 0.0 ) == pos &&
-       ( old -> corners[ c4 ] -> value > 0.0 ) == pos ) return;
-   // test bounds
-   if( bound_ps < bound_cr )
-   {
-      VECTOR p( position + cube_size * VECTOR( ( double ) i, ( double ) j, ( double ) k ) );
-      if( ! ( p >= bound_ps && p <= bound_cr ) ) return;
-   }
-   else
-      if( abs( i ) > 50 || abs( j ) > 50 || abs( k ) > 50 ) return;
-   if( Set_Center( i, j, k ) ) return;
-   CUBE* new_cube = new CUBE( i, j, k );
-   // given face of old cube is the same as the opposite face of new_cube
-   new_cube -> corners[ flip_bit( c1, bit ) ] = old -> corners[ c1 ];
-   new_cube -> corners[ flip_bit( c2, bit ) ] = old -> corners[ c2 ];
-   new_cube -> corners[ flip_bit( c3, bit ) ] = old -> corners[ c3 ];
-   new_cube -> corners[ flip_bit( c4, bit ) ] = old -> corners[ c4 ];
-   #ifdef DBG_POLYGONIZER
-   std::cout << "Test Face indexes: " << i << " " << j << " " << k << std::endl;
-   #endif
-   for( int n = 0; n < 8; n ++ )
-      if( ! new_cube -> corners[ n ] )
-         new_cube -> corners[ n ] = Set_Corner( i + ( ( n >> 2 ) & 1 ),
-                                                j + ( ( n >> 1 ) & 1 ),
-                                                k + ( n & 1 ) );
-   cube_stack. push( new_cube );
-}
-//---------------------------------------------------------------------------
-int POLYGONIZER :: Find_Point( int sign, VECTOR& point, double size, const SOLID* sld )
-{
-   VECTOR init = point;
-   VECTOR ps( -0.5 );
-   VECTOR cr(  0.5 );
-   for( int i = 0; i < 10000; i ++ )
-   {
-      point = init + size * Random_Vector( ps, cr );
-      //cout << point << " - " << sld -> Continuous_Function( point ) << std::endl;
-      if( sign == ( sld -> Sign_Function( point ) == 1 ) )
-         return 1;
-      size *= 1.005;
-   }
-   return 0;
-}
-//---------------------------------------------------------------------------
-void POLYGONIZER :: Free_Memory()
-{
-   // clear hashes
-   for( size_t i = 0; i < 2 * hash_size; i ++ )
-   {
-      if( edge_hash[ i ] )
-      {
-         list< EDGE* > :: iterator it = edge_hash[ i ] -> begin();
-         while( it != edge_hash[ i ] -> end() ) delete * it ++;
-         delete edge_hash[ i ];
-         edge_hash[ i ] = NULL;
-
-      }
-
-      if( corner_hash[ i ] )
-      {
-         list< CORNER* > :: iterator it = corner_hash[ i ] -> begin();
-         while( it != corner_hash[ i ] -> end() ) delete * it ++;
-         delete corner_hash[ i ];
-         corner_hash[ i ] = NULL;
-      }
-      if( center_hash[ i ] )
-      {
-         list< CENTER* > :: iterator it = center_hash[ i ] -> begin();
-         while( it != center_hash[ i ] -> end() ) delete * it ++;
-         delete center_hash[ i ];
-         center_hash[ i ] = NULL;
-      }
-   }
-}
-
diff --git a/src/Tools/polygonizer.h b/src/Tools/polygonizer.h
deleted file mode 100644
index d0d9e6df5f16528533e8a4eef519f4657a5275fb..0000000000000000000000000000000000000000
--- a/src/Tools/polygonizer.h
+++ /dev/null
@@ -1,194 +0,0 @@
-/***************************************************************************
-                          polygonizer.h  -  description
-                             -------------------
-    begin                : Mon Feb 11 2002
-    copyright            : (C) 2002 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/* See Copyright Notice in tnl/Copyright */
-
-#ifndef POLYGONIZER_H
-#define POLYGONIZER_H
-
-
-/**
-  *@author Tomas Oberhuber
-  *
-  * this class is rewritten polygonizer by Jules Bloomenthal, Xerox PARC.
-  * Copyright of original code (c) Xerox Corporation, 1991.
-  * See 'implicit.eecs.wsu.edu'.
-  *
-  */
-
-#include <list.h>
-#include <stack.h>
-#include <vector.h>
-#include "vctr.h"
-#include "vertex.h"
-
-//using namespace :: std;
-
-class SOLID;
-class RECIEVER;
-
-enum { poly_parametric, poly_implicit_cube, poly_implicit_tetrahedron };
-
-const int max_iter = 10;
-inline int flip_bit( int i, int bit )
-   { return i ^ 1 << bit; };
-
-
-enum { L = 0,  // left direction   - x
-       R,      // right direction  + x
-       B,      // bottom direction - y
-       T,      // top direction    + y
-       N,      // near direction   - z
-       F       // far direction    + z
-       };
-enum { LB = 0, // left bottom edge
-       LT,     // left top edge
-       LN,     // left near edge
-       LF,     // left far edge
-       RB,     // right bottom edge
-       RT,     // right top edge
-       RN,     // right near edge
-       RF,     // right far edge
-       BN,     // bottom near edge
-       BF,     // bottom far edge
-       TN,     // top near edge
-       TF      // top far edge
-       };
-
-enum { LBN = 0,   // left bottom near corner
-       LBF,       // left bottom far corner
-       LTN,       // left top near corner
-       LTF,       // left top far corner
-       RBN,       // right bottom near corner
-       RBF,       // right bottom far corner
-       RTN,       // right top near corner
-       RTF        // right top far corner
-       };
-
-struct CORNER
-{
-   int i, j, k;
-   VECTOR position;
-   double value;
-   CORNER( int _i, int _j, int _k, const VECTOR& p, const double& val )
-      :i( _i ), j( _j ), k( _k ), position( p ), value( val ){};
-};
-
-struct EDGE
-{
-   CORNER* corners[ 2 ];
-   VERTEX* vertex;
-   EDGE( CORNER* c1, CORNER* c2, VERTEX* v )
-      : vertex( v ){ corners[ 0 ] = c1; corners[ 1 ] = c2;};
-   ~EDGE() { assert( vertex ); delete vertex; };
-};
-
-struct CUBE
-{
-   int i, j, k;
-   // cube lattice
-   CORNER* corners[ 8 ];
-   CUBE( int, int, int );
-};
-
-struct CENTER
-{
-   int i, j, k;
-   CENTER( int _i, int _j, int _k )
-      : i( _i ), j( _j ), k( _k ){};
-};
-
-class POLYGONIZER
-{
-   public:
-   POLYGONIZER();
-   ~POLYGONIZER();
-
-   //! Creates initial cubes for implicite polygonizer
-   /*! We need to call this method several times in case of
-       non-connected surfaces. It si at least once called automaticly
-       if it was not called before starting the process of tesselation.
-    */
-   int Init( const VECTOR& pos, const SOLID* sld, const double& cb_sz );
-
-   //! Method for implicite surfaces
-   /*! start polygonizer with given SOLID, triangles are inserted into TRIANGLES
-    two const VECTOR&s define the bound for polygonizations
-    const double& is size of elementar cube
-    int is poligonization mode - polygonize cube or tetrahedrons
-    */
-   int Implicit( const SOLID*, RECIEVER*, const VECTOR&, const VECTOR&, const double&, int );
- 
-   int Parametric( const SOLID*, RECIEVER*, const unsigned int, const unsigned int );
-   protected:
- 
-   //! This solid will be polygonized using method SOLID :: Solid_Function( const VECTOR& );
-   const SOLID* solid;
-
-   //! Reciver of emitted triangles
-   RECIEVER* reciever;
-
-   //! Position of the first cube
-   /*! Thus, it is also position of the origin for cubes indexing
-    */
-   VECTOR position;
- 
-   //! Bounds for polygonizer
-   /*! All coubes outside culled.
-    */
-   VECTOR bound_ps, bound_cr;
-
-   //! Size of the cube and
-   double cube_size;
-
-   stack< CUBE* > cube_stack;
-   vector< vector< int > > cube_table[ 256 ];
-   list< EDGE* >** edge_hash;
-   // hash field of edges lists
-   list< CORNER* >** corner_hash;
-   // another hash for corners
-   list< CENTER* >** center_hash;
-   // yet another hash for cube centers
-
-   int Triangulate_Cube( CUBE* );
-   // triangulate the cube directly, without decomposition
-   VERTEX* Get_Vertex( CORNER*, CORNER* );
-   // return vertex for given edge using edge hash table ( Get_Edge method )
-   // both corners values are presumed of different sign
-   VECTOR Compute_Surface_Point( const VECTOR&, const VECTOR&, const double& );
-   // compute vertex on given edge using solid function of given solid
-   void Set_Edge( CORNER*, CORNER*, VERTEX* );
-   // insert edge to hash table
-   int Check_Edge( CORNER*, CORNER*, VERTEX*& );
-   // get vertex on edge from hash table
-   // return 1 if edge is in hash table
-   //        0 if edge is not in hash table
-   CORNER* Set_Corner( int, int, int );
-   // return corner with given lattice location
-   // set and cache its function value
-   int Set_Center( int, int, int );
-   // set ( i, j, k ) entry to center_hash table
-   // return 1 if already set; otherwise set and return 0
-   void Make_Cube_Table();
-   // create cube_table
-   int Next_Clockwise_Edge( int, int );
-   // return next clockwise edge from given edge around given face
-   int Other_Face( int, int );
-   // return face adjoining edge is not the given face
-   void Test_Face( int, int, int, CUBE*, int, int, int, int, int );
-   // given cube at lattice ( i, j, k ) and four corners of face,
-   // if surface crosses face, compute other four corners of adjacent cube
-   // and add new cube to cube stack
-   int Find_Point( int, VECTOR&, double, const SOLID* sld );
-   // find point with given value sign
-   void Free_Memory();
-};
-
-extern POLYGONIZER polygonizer;
-
-#endif
diff --git a/src/Tools/tnl-compile.in b/src/Tools/tnl-compile.in
index 1a813fe2f41d1bb5ed6a9a60e5baff3cfeb038ce..b106c613c368681841bcf36f0388657d553da942 100644
--- a/src/Tools/tnl-compile.in
+++ b/src/Tools/tnl-compile.in
@@ -12,5 +12,5 @@ do
     esac
 done
 
-echo -I@CMAKE_INSTALL_PREFIX@/@TNL_TARGET_INCLUDE_DIRECTORY@ ${CUDA_FLAGS} ${CXX_STD_FLAGS} ${DEBUG_FLAGS}
+echo -I@CMAKE_INSTALL_PREFIX@/include ${CUDA_FLAGS} ${CXX_STD_FLAGS} ${DEBUG_FLAGS}
 
diff --git a/src/Tools/tnl-eoc-test-log b/src/Tools/tnl-eoc-test-log
deleted file mode 100644
index c89d141eb32d707a1e1ade87c0692442d6a3e2cb..0000000000000000000000000000000000000000
--- a/src/Tools/tnl-eoc-test-log
+++ /dev/null
@@ -1,4 +0,0 @@
-#!/usr/bin/env python
-
-import sys, string, math
-
diff --git a/src/Tools/tnl-quickstart/build-config-tag.h.in b/src/Tools/tnl-quickstart/build-config-tag.h.in
index c72d8c1b9dd30c104ff3530cef8c715c04a0e7ee..f507baf061cbd1d76c187cba42d54dcc3c765cdb 100644
--- a/src/Tools/tnl-quickstart/build-config-tag.h.in
+++ b/src/Tools/tnl-quickstart/build-config-tag.h.in
@@ -20,31 +20,21 @@ template<> struct ConfigTagIndex< {problemBaseName}BuildConfigTag, short int >{{
 template<> struct ConfigTagIndex< {problemBaseName}BuildConfigTag, long int >{{ enum {{ enabled = false }}; }};
     
 /****
- * With how many dimensions may have the problem to be solved...
- */    
-template< int Dimension > struct ConfigTagDimension< {problemBaseName}BuildConfigTag, Dimension >{{ enum {{ enabled = ( Dimension == 1 ) }}; }};
-
-/****
- * Use of Meshes::Grid is enabled for allowed dimensions and Real, Device and Index types.
+ * The mesh type will be resolved by the Solver by default.
+ * (The detailed mesh configuration is in TNL/Meshes/BuildConfigTags.h)
  */
-template< int Dimension, typename Real, typename Device, typename Index >
-   struct ConfigTagMesh< {problemBaseName}BuildConfigTag, Meshes::Grid< Dimension, Real, Device, Index > >
-      {{ enum {{ enabled = ConfigTagDimension< {problemBaseName}BuildConfigTag, Dimension >::enabled  &&
-                         ConfigTagReal< {problemBaseName}BuildConfigTag, Real >::enabled &&
-                         ConfigTagDevice< {problemBaseName}BuildConfigTag, Device >::enabled &&
-                         ConfigTagIndex< {problemBaseName}BuildConfigTag, Index >::enabled }}; }};
+template<> struct ConfigTagMeshResolve< {problemBaseName}BuildConfigTag >{{ enum {{ enabled = true }}; }};
 
 /****
- * Please, chose your preferred time discretization  here.
+ * All time discretisations (explicit, semi-impicit and implicit ) are
+ * enabled by default.
  */
-template<> struct ConfigTagTimeDiscretisation< {problemBaseName}BuildConfigTag, ExplicitTimeDiscretisationTag >{{ enum {{ enabled = true }}; }};
-template<> struct ConfigTagTimeDiscretisation< {problemBaseName}BuildConfigTag, SemiImplicitTimeDiscretisationTag >{{ enum {{ enabled = true }}; }};
-template<> struct ConfigTagTimeDiscretisation< {problemBaseName}BuildConfigTag, ImplicitTimeDiscretisationTag >{{ enum {{ enabled = false }}; }};
+template< typename TimeDiscretisation > struct ConfigTagTimeDiscretisation< {problemBaseName}BuildConfigTag, TimeDiscretisation >{{ enum {{ enabled = true }}; }};
 
 /****
- * Only the Runge-Kutta-Merson solver is enabled by default.
+ * All explicit solvers are enabled by default
  */
-template<> struct ConfigTagExplicitSolver< {problemBaseName}BuildConfigTag, ExplicitEulerSolverTag >{{ enum {{ enabled = false }}; }};
+template< typename ExplicitSolver > struct ConfigTagExplicitSolver< {problemBaseName}BuildConfigTag, ExplicitSolver >{{ enum {{ enabled = true }}; }};
 
 }} // namespace Solvers
-}} // namespace TNL
\ No newline at end of file
+}} // namespace TNL
diff --git a/src/Tools/tnl-quickstart/main.h.in b/src/Tools/tnl-quickstart/main.h.in
index f836967c83ccc88fed8898a0c344a9ea08e40906..305f59e1a21757db59c91a2ca80d8e9772c11f72 100644
--- a/src/Tools/tnl-quickstart/main.h.in
+++ b/src/Tools/tnl-quickstart/main.h.in
@@ -47,7 +47,8 @@ template< typename Real,
           typename Index,
           typename MeshType,
           typename ConfigTag,
-          typename SolverStarter >
+          typename SolverStarter,
+          typename Communicator >
 class {problemBaseName}Setter
 {{
    public:
@@ -75,12 +76,12 @@ class {problemBaseName}Setter
              if( boundaryConditionsType == "dirichlet" )
              {{
                 typedef Operators::DirichletBoundaryConditions< MeshType, ConstantFunction, MeshType::getMeshDimension(), Real, Index > BoundaryConditions;
-                typedef {problemBaseName}Problem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
+                typedef {problemBaseName}Problem< MeshType, Communicator, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
                 SolverStarter solverStarter;
                 return solverStarter.template run< Problem >( parameters );
              }}
              typedef Operators::NeumannBoundaryConditions< MeshType, ConstantFunction, Real, Index > BoundaryConditions;
-             typedef {problemBaseName}Problem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
+             typedef {problemBaseName}Problem< MeshType, Communicator, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
              SolverStarter solverStarter;
              return solverStarter.template run< Problem >( parameters );
           }}
@@ -88,12 +89,12 @@ class {problemBaseName}Setter
           if( boundaryConditionsType == "dirichlet" )
           {{    
              typedef Operators::DirichletBoundaryConditions< MeshType, MeshFunction, MeshType::getMeshDimension(), Real, Index > BoundaryConditions;
-             typedef {problemBaseName}Problem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
+             typedef {problemBaseName}Problem< MeshType, Communicator, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
              SolverStarter solverStarter;
              return solverStarter.template run< Problem >( parameters );
           }}
           typedef Operators::NeumannBoundaryConditions< MeshType, MeshFunction, Real, Index > BoundaryConditions;
-          typedef {problemBaseName}Problem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
+          typedef {problemBaseName}Problem< MeshType, Communicator, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
           SolverStarter solverStarter;
           return solverStarter.template run< Problem >( parameters );
       }}
diff --git a/src/Tools/tnl-quickstart/operator-grid-specialization_impl.h.in b/src/Tools/tnl-quickstart/operator-grid-specialization_impl.h.in
index 3e7d0d670f5de48659c146bbf7b1fb62287e4a38..da4da6d635d10d1681d689972d5e93695e47b4dd 100644
--- a/src/Tools/tnl-quickstart/operator-grid-specialization_impl.h.in
+++ b/src/Tools/tnl-quickstart/operator-grid-specialization_impl.h.in
@@ -34,7 +34,7 @@ operator()( const MeshFunction& u,
     * The following example is the Laplace operator approximated 
     * by the Finite difference method.
     */  
-   static_assert( MeshEntity::entityDimension == {meshDimension}, "Wrong mesh entity dimensions." );
+   static_assert( MeshEntity::getEntityDimension() == {meshDimension}, "Wrong mesh entity dimensions." );
    static_assert( MeshFunction::getEntitiesDimension() == {meshDimension}, "Wrong preimage function" );
    const typename MeshEntity::template NeighborEntities< {meshDimension} >& neighborEntities = entity.getNeighborEntities(); 
 
@@ -84,7 +84,7 @@ setMatrixElements( const PreimageFunction& u,
                    Matrix& matrix,
                    Vector& b ) const
 {{
-   static_assert( MeshEntity::entityDimension == {meshDimension}, "Wrong mesh entity dimensions." );
+   static_assert( MeshEntity::getEntityDimension() == {meshDimension}, "Wrong mesh entity dimensions." );
    static_assert( PreimageFunction::getEntitiesDimension() == {meshDimension}, "Wrong preimage function" );
 
    /****
diff --git a/src/Tools/tnl-quickstart/problem.h.in b/src/Tools/tnl-quickstart/problem.h.in
index dd53ffe43f77e58a5dd49cdc51306820a0b07249..9006f7cf7c5f1e39e5b59333fe888e5d56b479d5 100644
--- a/src/Tools/tnl-quickstart/problem.h.in
+++ b/src/Tools/tnl-quickstart/problem.h.in
@@ -8,33 +8,35 @@
 
 
 template< typename Mesh,
+          typename Communicator,
           typename BoundaryCondition,
           typename RightHandSide,
           typename DifferentialOperator >
 class {problemBaseName}Problem:
    public TNL::Problems::PDEProblem< Mesh,
-                         typename DifferentialOperator::RealType,
-                         typename Mesh::DeviceType,
-                         typename DifferentialOperator::IndexType >
+                                     Communicator,
+                                     typename DifferentialOperator::RealType,
+                                     typename Mesh::DeviceType,
+                                     typename Mesh::IndexType >
 {{
    public:
 
       typedef typename DifferentialOperator::RealType RealType;
       typedef typename Mesh::DeviceType DeviceType;
-      typedef typename DifferentialOperator::IndexType IndexType;
-      typedef TNL::Functions::MeshFunction< Mesh > MeshFunctionType;
-      typedef TNL::Problems::PDEProblem< Mesh, RealType, DeviceType, IndexType > BaseType;
-      typedef TNL::SharedPointer< DifferentialOperator > DifferentialOperatorPointer;
-      typedef TNL::SharedPointer< BoundaryCondition > BoundaryConditionPointer;
-      typedef TNL::SharedPointer< RightHandSide, DeviceType > RightHandSidePointer;
+      typedef typename Mesh::IndexType IndexType;
+      typedef TNL::Functions::MeshFunction< Mesh, Mesh::getMeshDimension(), RealType > MeshFunctionType;
+      typedef TNL::Problems::PDEProblem< Mesh, Communicator, RealType, DeviceType, IndexType > BaseType;
+      typedef TNL::Pointers::SharedPointer< MeshFunctionType > MeshFunctionPointer;
+      typedef TNL::Pointers::SharedPointer< DifferentialOperator > DifferentialOperatorPointer;
+      typedef TNL::Pointers::SharedPointer< BoundaryCondition > BoundaryConditionPointer;
+      typedef TNL::Pointers::SharedPointer< RightHandSide, DeviceType > RightHandSidePointer;
 
       using typename BaseType::MeshType;
       using typename BaseType::MeshPointer;
       using typename BaseType::DofVectorType;
       using typename BaseType::DofVectorPointer;
-      using typename BaseType::MeshDependentDataType;
-      using typename BaseType::MeshDependentDataPointer;
 
+      using CommunicatorType = Communicator;
 
       static TNL::String getTypeStatic();
 
@@ -43,51 +45,42 @@ class {problemBaseName}Problem:
       void writeProlog( TNL::Logger& logger,
                         const TNL::Config::ParameterContainer& parameters ) const;
 
-      bool setup( const MeshPointer& meshPointer,
-                  const TNL::Config::ParameterContainer& parameters,
+      bool setup( const TNL::Config::ParameterContainer& parameters,
                   const TNL::String& prefix );
 
 
       bool setInitialCondition( const TNL::Config::ParameterContainer& parameters,
-                                const MeshPointer& mesh,
-                                DofVectorPointer& dofs,
-                                MeshDependentDataPointer& meshDependentData );
+                                DofVectorPointer& dofs );
 
       template< typename MatrixPointer >
-      bool setupLinearSystem( const MeshPointer& mesh,
-                              MatrixPointer& matrixPointer );
+      bool setupLinearSystem( MatrixPointer& matrixPointer );
 
       bool makeSnapshot( const RealType& time,
                          const IndexType& step,
-                         const MeshPointer& mesh,
-                         DofVectorPointer& dofs,
-                         MeshDependentDataPointer& meshDependentData );
+                         DofVectorPointer& dofs );
 
-      IndexType getDofs( const MeshPointer& mesh ) const;
+      IndexType getDofs() const;
 
-      void bindDofs( const MeshPointer& mesh,
-                     DofVectorPointer& dofs );
+      void bindDofs( DofVectorPointer& dofs );
 
       void getExplicitUpdate( const RealType& time,
                               const RealType& tau,
-                              const MeshPointer& mesh,
                               DofVectorPointer& _u,
-                              DofVectorPointer& _fu,
-                              MeshDependentDataPointer& meshDependentData );
+                              DofVectorPointer& _fu );
 
       template< typename MatrixPointer >
       void assemblyLinearSystem( const RealType& time,
                                  const RealType& tau,
-                                 const MeshPointer& mesh,
                                  DofVectorPointer& dofs,
                                  MatrixPointer& matrixPointer,
-                                 DofVectorPointer& rightHandSide,
-                                 MeshDependentDataPointer& meshDependentData );
+                                 DofVectorPointer& rightHandSide );
 
    protected:
     
       DifferentialOperatorPointer differentialOperator;
+
       BoundaryConditionPointer boundaryCondition;
+
       RightHandSidePointer rightHandSide;
    
       TNL::Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, DifferentialOperator, BoundaryCondition, RightHandSide > explicitUpdater;
diff --git a/src/Tools/tnl-quickstart/problem_impl.h.in b/src/Tools/tnl-quickstart/problem_impl.h.in
index d0b1162383eea26cf35282ecb78b9b5143f4abdb..f196ebcec1922b51539ca2f5794ba8b8324be368 100644
--- a/src/Tools/tnl-quickstart/problem_impl.h.in
+++ b/src/Tools/tnl-quickstart/problem_impl.h.in
@@ -8,33 +8,36 @@
 #include <TNL/Solvers/PDE/BackwardTimeDiscretisation.h>
 
 template< typename Mesh,
+          typename Communicator,
           typename BoundaryCondition,
           typename RightHandSide,
           typename DifferentialOperator >
 TNL::String
-{problemBaseName}Problem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+{problemBaseName}Problem< Mesh, Communicator, BoundaryCondition, RightHandSide, DifferentialOperator >::
 getTypeStatic()
 {{
    return TNL::String( "{problemBaseName}Problem< " ) + Mesh :: getTypeStatic() + " >";
 }}
 
 template< typename Mesh,
+          typename Communicator,
           typename BoundaryCondition,
           typename RightHandSide,
           typename DifferentialOperator >
 TNL::String 
-{problemBaseName}Problem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+{problemBaseName}Problem< Mesh, Communicator, BoundaryCondition, RightHandSide, DifferentialOperator >::
 getPrologHeader() const
 {{    
    return TNL::String( "{problemName}" );
 }}
 
 template< typename Mesh,
+          typename Communicator,
           typename BoundaryCondition,
           typename RightHandSide,
           typename DifferentialOperator >
 void
-{problemBaseName}Problem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::    
+{problemBaseName}Problem< Mesh, Communicator, BoundaryCondition, RightHandSide, DifferentialOperator >::    
 writeProlog( TNL::Logger& logger, const TNL::Config::ParameterContainer& parameters ) const
 {{
    /****
@@ -44,60 +47,60 @@ writeProlog( TNL::Logger& logger, const TNL::Config::ParameterContainer& paramet
 }}
 
 template< typename Mesh,
+          typename Communicator,
           typename BoundaryCondition,
           typename RightHandSide,
           typename DifferentialOperator >
 bool
-{problemBaseName}Problem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::        
-setup( const MeshPointer& meshPointer,
-       const TNL::Config::ParameterContainer& parameters,
+{problemBaseName}Problem< Mesh, Communicator, BoundaryCondition, RightHandSide, DifferentialOperator >::        
+setup( const TNL::Config::ParameterContainer& parameters,
        const TNL::String& prefix )
 {{
-   if( ! this->boundaryCondition->setup( meshPointer, parameters, "boundary-conditions-" ) ||
+   if( ! this->boundaryCondition->setup( this->getMesh(), parameters, "boundary-conditions-" ) ||
        ! this->rightHandSide->setup( parameters, "right-hand-side-" ) )
       return false;
    return true;
 }}
 
 template< typename Mesh,
+          typename Communicator,
           typename BoundaryCondition,
           typename RightHandSide,
           typename DifferentialOperator >
-typename {problemBaseName}Problem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::IndexType
-    {problemBaseName}Problem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::            
-getDofs( const MeshPointer& mesh ) const
+typename {problemBaseName}Problem< Mesh, Communicator, BoundaryCondition, RightHandSide, DifferentialOperator >::IndexType
+    {problemBaseName}Problem< Mesh, Communicator, BoundaryCondition, RightHandSide, DifferentialOperator >::            
+getDofs() const
 {{
    /****
     * Return number of  DOFs (degrees of freedom) i.e. number
     * of unknowns to be resolved by the main solver.
     */
-   return mesh->template getEntitiesCount< typename MeshType::Cell >();
+   return this->getMesh()->template getEntitiesCount< typename MeshType::Cell >();
 }}
 
 template< typename Mesh,
+          typename Communicator,
           typename BoundaryCondition,
           typename RightHandSide,
           typename DifferentialOperator >
 void
-{problemBaseName}Problem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::            
-bindDofs( const MeshPointer& mesh,
-          DofVectorPointer& dofVector )    
+{problemBaseName}Problem< Mesh, Communicator, BoundaryCondition, RightHandSide, DifferentialOperator >::            
+bindDofs( DofVectorPointer& dofVector )    
 {{
 }}
 
 template< typename Mesh,
+          typename Communicator,
           typename BoundaryCondition,
           typename RightHandSide,
           typename DifferentialOperator >
 bool
-{problemBaseName}Problem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::            
+{problemBaseName}Problem< Mesh, Communicator, BoundaryCondition, RightHandSide, DifferentialOperator >::            
 setInitialCondition( const TNL::Config::ParameterContainer& parameters,    
-                     const MeshPointer& mesh,
-                     DofVectorPointer& dofs,
-                     MeshDependentDataPointer& meshDependentData )
+                     DofVectorPointer& dofs )
 {{
    const TNL::String& initialConditionFile = parameters.getParameter< TNL::String >( "initial-condition" );
-   TNL::Functions::MeshFunction< Mesh > u( mesh, dofs );
+   MeshFunctionType u( this->getMesh(), dofs );
    if( ! u.boundLoad( initialConditionFile ) )
    {{
       std::cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." << std::endl;
@@ -107,45 +110,42 @@ setInitialCondition( const TNL::Config::ParameterContainer& parameters,
 }}
 
 template< typename Mesh,
+          typename Communicator,
           typename BoundaryCondition,
           typename RightHandSide,
           typename DifferentialOperator >
    template< typename MatrixPointer >
 bool
-{problemBaseName}Problem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::                
-setupLinearSystem( const MeshPointer& meshPointer,
-                   MatrixPointer& matrixPointer )
+{problemBaseName}Problem< Mesh, Communicator, BoundaryCondition, RightHandSide, DifferentialOperator >::                
+setupLinearSystem( MatrixPointer& matrixPointer )
 {{
-   const IndexType dofs = this->getDofs( meshPointer );
+   const IndexType dofs = this->getDofs();
    typedef typename MatrixPointer::ObjectType::CompressedRowLengthsVector CompressedRowLengthsVectorType;
-   TNL::SharedPointer< CompressedRowLengthsVectorType > rowLengthsPointer;
-   if( ! rowLengthsPointer->setSize( dofs ) )
-      return false;
+   TNL::Pointers::SharedPointer< CompressedRowLengthsVectorType > rowLengthsPointer;
+   rowLengthsPointer->setSize( dofs );
    TNL::Matrices::MatrixSetter< MeshType, DifferentialOperator, BoundaryCondition, CompressedRowLengthsVectorType > matrixSetter;
-   matrixSetter.template getCompressedRowLengths< typename Mesh::Cell >( meshPointer,
+   matrixSetter.template getCompressedRowLengths< typename Mesh::Cell >( this->getMesh(),
                                                                          differentialOperator,
                                                                          boundaryCondition,
                                                                          rowLengthsPointer );
    matrixPointer->setDimensions( dofs, dofs );
-   if( ! matrixPointer->setCompressedRowLengths( *rowLengthsPointer ) )
-      return false;
+   matrixPointer->setCompressedRowLengths( *rowLengthsPointer );
    return true;
 }}
 
 template< typename Mesh,
+          typename Communicator,
           typename BoundaryCondition,
           typename RightHandSide,
           typename DifferentialOperator >
 bool
-{problemBaseName}Problem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::                    
+{problemBaseName}Problem< Mesh, Communicator, BoundaryCondition, RightHandSide, DifferentialOperator >::                    
 makeSnapshot( const RealType& time,
               const IndexType& step,
-              const MeshPointer& mesh,
-              DofVectorPointer& dofs,
-              MeshDependentDataPointer& meshDependentData )
+              DofVectorPointer& dofs )
 {{
    std::cout << std::endl << "Writing output at time " << time << " step " << step << "." << std::endl;
-   this->bindDofs( mesh, dofs );
+   this->bindDofs( dofs );
    TNL::FileName fileName;
    fileName.setFileNameBase( "u-" );
    fileName.setExtension( "tnl" );
@@ -156,17 +156,16 @@ makeSnapshot( const RealType& time,
 }}
 
 template< typename Mesh,
+          typename Communicator,
           typename BoundaryCondition,
           typename RightHandSide,
           typename DifferentialOperator >
 void
-{problemBaseName}Problem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::                    
+{problemBaseName}Problem< Mesh, Communicator, BoundaryCondition, RightHandSide, DifferentialOperator >::                    
 getExplicitUpdate( const RealType& time,
                    const RealType& tau,
-                   const MeshPointer& mesh,
                    DofVectorPointer& _u,
-                   DofVectorPointer& _fu,
-                   MeshDependentDataPointer& meshDependentData )
+                   DofVectorPointer& _fu )
 {{
    /****
     * If you use an explicit solver like EulerSolver or MersonSolver, you
@@ -177,28 +176,27 @@ getExplicitUpdate( const RealType& time,
     * You may use supporting mesh dependent data if you need.
     */
 
-   TNL::SharedPointer< MeshFunctionType > uPointer( mesh, _u ); 
-   TNL::SharedPointer< MeshFunctionType > fuPointer( mesh, _fu ); 
-   this->explicitUpdater.setDifferentialOperator( this->differentialOperator ),
-   this->explicitUpdater.setBoundaryConditions( this->boundaryCondition ),
-   this->explicitUpdater.setRightHandSide( this->rightHandSide ),
-   this->explicitUpdater.template update< typename Mesh::Cell >( time, tau, mesh, uPointer, fuPointer );
+   TNL::Pointers::SharedPointer< MeshFunctionType > uPointer( this->getMesh(), _u ); 
+   TNL::Pointers::SharedPointer< MeshFunctionType > fuPointer( this->getMesh(), _fu ); 
+   this->explicitUpdater.setDifferentialOperator( this->differentialOperator );
+   this->explicitUpdater.setBoundaryConditions( this->boundaryCondition );
+   this->explicitUpdater.setRightHandSide( this->rightHandSide );
+   this->explicitUpdater.template update< typename Mesh::Cell, CommunicatorType >( time, tau, this->getMesh(), uPointer, fuPointer );
 }}
 
 template< typename Mesh,
+          typename Communicator,
           typename BoundaryCondition,
           typename RightHandSide,
           typename DifferentialOperator >
    template< typename MatrixPointer >
 void
-{problemBaseName}Problem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::                
+{problemBaseName}Problem< Mesh, Communicator, BoundaryCondition, RightHandSide, DifferentialOperator >::                
 assemblyLinearSystem( const RealType& time,
                       const RealType& tau,
-                      const MeshPointer& mesh,
                       DofVectorPointer& _u,
                       MatrixPointer& matrixPointer,
-                      DofVectorPointer& b,
-                      MeshDependentDataPointer& meshDependentData )
+                      DofVectorPointer& b )
 {{
    /****
     * If you implement a (semi-)implicit solver, this method is supposed
@@ -206,14 +204,14 @@ assemblyLinearSystem( const RealType& time,
     * You may use supporting mesh dependent data if you need.
     */
 
-   TNL::SharedPointer< TNL::Functions::MeshFunction< Mesh > > uPointer( mesh, _u );
+   MeshFunctionPointer uPointer( this->getMesh(), _u );
    this->systemAssembler.setDifferentialOperator( this->differentialOperator );
    this->systemAssembler.setBoundaryConditions( this->boundaryCondition );
    this->systemAssembler.setRightHandSide( this->rightHandSide );
    this->systemAssembler.template assembly< typename Mesh::Cell, typename MatrixPointer::ObjectType >( 
       time,
       tau,
-      mesh,
+      this->getMesh(),
       uPointer,
       matrixPointer,
       b );
diff --git a/src/Tools/tnl-quickstart/rhs.h.in b/src/Tools/tnl-quickstart/rhs.h.in
index ffd1b78813dbfcd2f17809fcea69f2d999e86b89..07a97832d03f62e45912a3f229b52f71a31db11f 100644
--- a/src/Tools/tnl-quickstart/rhs.h.in
+++ b/src/Tools/tnl-quickstart/rhs.h.in
@@ -4,7 +4,7 @@
 
 template< typename Mesh, typename Real >
 class {problemBaseName}Rhs
-  : public TNL::Functions::Domain< Mesh::meshDimension, TNL::Functions::MeshDomain >
+  : public TNL::Functions::Domain< Mesh::getMeshDimension(), TNL::Functions::MeshDomain >
 {{
    public:
 
diff --git a/src/Tools/tnl-time-series2png b/src/Tools/tnl-time-series2png
deleted file mode 100644
index 8d59b386b8f15536fb7d360e63d1b00c2d02f332..0000000000000000000000000000000000000000
--- a/src/Tools/tnl-time-series2png
+++ /dev/null
@@ -1,36 +0,0 @@
-#!/usr/bin/env python
-
-class timeSerie:
-   def __init__( self, index, label ):
-      self.index = index
-      self.label = label
-
-inputFile = ""
-outputFile = ""
-verbose = "false"
-timeSeries = []
-
-i = 0
-while i < len( arguments ):
-   if arguments[ i ] == "--input-file":
-      refinement = arguments[ i + 1 ]
-      i = i + 2
-      continue
-   if arguments[ i ] == "--output-file":
-      output_file_name = arguments[ i + 1 ]
-      i = i + 2
-      continue
-   if arguments[ i ] == "--time-serie":
-      index = arguments[ i + 1 ]
-      label = arguments[ i + 2 ]
-      timeSeries.append( timeSerie( index, label ) )
-      i = i + 3
-      continue
-   if arguments[ i ] == "--verbose":
-       verbose = float( arguments[ i + 1 ] )
-       i = i +2
-       continue       
-   input_files. append( arguments[ i ] )
-   i = i + 1
-   
-