diff --git a/src/TNL/Assert.h b/src/TNL/Assert.h
index f8f60200d0bdfdc4910334174bb7b9ecd3a63977..cb117501fbe72b8b845ff1463a5c02f9c31009fd 100644
--- a/src/TNL/Assert.h
+++ b/src/TNL/Assert.h
@@ -10,26 +10,246 @@
 
 #pragma once
 
+#include <TNL/Devices/CudaCallable.h>
+
 /****
  * Debugging assert
  */
 
 #ifndef NDEBUG
 
+#include <sstream>
 #include <iostream>
-#include <stdlib.h>
-#include <assert.h>
+#include <stdio.h>
+
+namespace TNL {
+namespace Assert {
+
+inline void
+printDiagnosticsHost( const char* assertion,
+                      const char* message,
+                      const char* file,
+                      const char* function,
+                      int line,
+                      const char* diagnostics )
+{
+   std::cerr << "Assertion '" << assertion << "' failed !!!\n"
+             << "Message: " << message << "\n"
+             << "File: " << file << "\n"
+             << "Function: " << function << "\n"
+             << "Line: " << line << "\n"
+             << "Diagnostics:\n" << diagnostics << std::endl;
+}
+
+__cuda_callable__
+inline void
+printDiagnosticsCuda( const char* assertion,
+                      const char* message,
+                      const char* file,
+                      const char* function,
+                      int line,
+                      const char* diagnostics )
+{
+   printf( "Assertion '%s' failed !!!\n"
+           "Message: %s\n"
+           "File: %s\n"
+           "Function: %s\n"
+           "Line: %d\n"
+           "Diagnostics: %s\n",
+           assertion, message, file, function, line, diagnostics );
+}
 
+__cuda_callable__
+inline void
+fatalFailure()
+{
+#ifdef __CUDA_ARCH__
+   // https://devtalk.nvidia.com/default/topic/509584/how-to-cancel-a-running-cuda-kernel-/
+   // TODO: it is reported as "illegal instruction", but that leads to an abort as well...
+   asm("trap;");
+#else
+   throw EXIT_FAILURE;
 #endif
+}
 
-#ifndef NDEBUG   
-   
+template< typename T >
+std::string
+printToString( const T& value )
+{
+   ::std::stringstream ss;
+   ss << value;
+   return ss.str();
+}
+
+template<>
+inline std::string
+printToString( const bool& value )
+{
+   if( value ) return "true";
+   else return "false";
+}
+
+template< typename T1, typename T2 >
+__cuda_callable__ void
+cmpHelperOpFailure( const char* assertion,
+                    const char* message,
+                    const char* file,
+                    const char* function,
+                    int line,
+                    const char* lhs_expression,
+                    const char* rhs_expression,
+                    const T1& lhs_value,
+                    const T2& rhs_value,
+                    const char* op )
+{
+#ifdef __CUDA_ARCH__
+   // diagnostics is not supported - we don't have the machinery
+   // to construct the dynamic error message
+   printDiagnosticsCuda( assertion, message, file, function, line,
+                         "Not supported in CUDA kernels." );
+#else
+   std::stringstream str;
+   if( std::string(op) == "==" ) {
+      str << "      Expected: " << lhs_expression;
+      if( printToString(lhs_value) != lhs_expression ) {
+         str << "\n      Which is: " << lhs_value;
+      }
+      str << "\nTo be equal to: " << rhs_expression;
+      if( printToString(rhs_value) != rhs_expression ) {
+         str << "\n      Which is: " << rhs_value;
+      }
+      str << std::endl;
+   }
+   else {
+      str << "Expected: (" << lhs_expression << ") " << op << " (" << rhs_expression << "), "
+          << "actual: " << lhs_value << " vs " << rhs_value << std::endl;
+   }
+   printDiagnosticsHost( assertion, message, file, function, line,
+                         str.str().c_str() );
+#endif
+   fatalFailure();
+}
+
+template< typename T1, typename T2 >
+__cuda_callable__ void
+cmpHelperTrue( const char* assertion,
+               const char* message,
+               const char* file,
+               const char* function,
+               int line,
+               const char* expr1,
+               const char* expr2,
+               const T1& val1,
+               const T2& val2 )
+{
+   // explicit cast is necessary, because T1::operator! might not be defined
+   if( ! (bool) val1 )
+      ::TNL::Assert::cmpHelperOpFailure( assertion, message, file, function, line,
+                                         expr1, "true", val1, true, "==" );
+}
+
+template< typename T1, typename T2 >
+__cuda_callable__ void
+cmpHelperFalse( const char* assertion,
+                const char* message,
+                const char* file,
+                const char* function,
+                int line,
+                const char* expr1,
+                const char* expr2,
+                const T1& val1,
+                const T2& val2 )
+{
+   if( val1 )
+      ::TNL::Assert::cmpHelperOpFailure( assertion, message, file, function, line,
+                                         expr1, "false", val1, false, "==" );
+}
+
+// A macro for implementing the helper functions needed to implement
+// TNL_ASSERT_??. It is here just to avoid copy-and-paste of similar code.
+#define TNL_IMPL_CMP_HELPER_( op_name, op ) \
+template< typename T1, typename T2 > \
+__cuda_callable__ void \
+cmpHelper##op_name( const char* assertion, \
+                    const char* message, \
+                    const char* file, \
+                    const char* function, \
+                    int line, \
+                    const char* expr1, \
+                    const char* expr2, \
+                    const T1& val1, \
+                    const T2& val2 ) \
+{\
+   if( ! ( (val1) op (val2) ) ) \
+      ::TNL::Assert::cmpHelperOpFailure( assertion, message, file, function, line, \
+                                         expr1, expr2, val1, val2, #op );\
+}
+
+// Implements the helper function for TNL_ASSERT_EQ
+TNL_IMPL_CMP_HELPER_( EQ, == );
+// Implements the helper function for TNL_ASSERT_NE
+TNL_IMPL_CMP_HELPER_( NE, != );
+// Implements the helper function for TNL_ASSERT_LE
+TNL_IMPL_CMP_HELPER_( LE, <= );
+// Implements the helper function for TNL_ASSERT_LT
+TNL_IMPL_CMP_HELPER_( LT, < );
+// Implements the helper function for TNL_ASSERT_GE
+TNL_IMPL_CMP_HELPER_( GE, >= );
+// Implements the helper function for TNL_ASSERT_GT
+TNL_IMPL_CMP_HELPER_( GT, > );
+
+#undef TNL_IMPL_CMP_HELPER_
+
+} // namespace Assert
+} // namespace TNL
+
+// Internal macro wrapping the __PRETTY_FUNCTION__ "magic".
 #if defined( __NVCC__ ) && ( __CUDACC_VER__ < 80000 )
-    #define TNL_PRETTY_FUNCTION "(not known in CUDA 7.5 or older)"
+    #define __TNL_PRETTY_FUNCTION "(not known in CUDA 7.5 or older)"
 #else
-    #define TNL_PRETTY_FUNCTION __PRETTY_FUNCTION__
+    #define __TNL_PRETTY_FUNCTION __PRETTY_FUNCTION__
 #endif
 
+// Internal macro to compose the string representing the assertion.
+// We can't do it easily at runtime, because we have to support assertions
+// in CUDA kernels, which can't use std::string objects. Instead, we do it
+// at compile time - adjacent strings are joined at the language level.
+#define __TNL_JOIN_STRINGS( val1, op, val2 ) \
+   __STRING( val1 ) " " __STRING( op ) " " __STRING( val2 )
+
+// Internal macro to pass all the arguments to the specified cmpHelperOP
+#define __TNL_ASSERT_PRED2( pred, op, val1, val2, msg ) \
+   pred( __TNL_JOIN_STRINGS( val1, op, val2 ), \
+         msg, __FILE__, __TNL_PRETTY_FUNCTION, __LINE__, \
+         #val1, #val2, val1, val2 )
+   
+// Main definitions of the TNL_ASSERT_* macros
+// unary
+#define TNL_ASSERT_TRUE( val, msg ) \
+   __TNL_ASSERT_PRED2( ::TNL::Assert::cmpHelperTrue, ==, val, true, msg )
+#define TNL_ASSERT_FALSE( val, msg ) \
+   __TNL_ASSERT_PRED2( ::TNL::Assert::cmpHelperFalse, ==, val, false, msg )
+// binary
+#define TNL_ASSERT_EQ( val1, val2, msg ) \
+   __TNL_ASSERT_PRED2( ::TNL::Assert::cmpHelperEQ, ==, val1, val2, msg )
+#define TNL_ASSERT_NE( val1, val2, msg ) \
+   __TNL_ASSERT_PRED2( ::TNL::Assert::cmpHelperNE, !=, val1, val2, msg )
+#define TNL_ASSERT_LE( val1, val2, msg ) \
+   __TNL_ASSERT_PRED2( ::TNL::Assert::cmpHelperLE, <=, val1, val2, msg )
+#define TNL_ASSERT_LT( val1, val2, msg ) \
+   __TNL_ASSERT_PRED2( ::TNL::Assert::cmpHelperLT, <,  val1, val2, msg )
+#define TNL_ASSERT_GE( val1, val2, msg ) \
+   __TNL_ASSERT_PRED2( ::TNL::Assert::cmpHelperGE, >=, val1, val2, msg )
+#define TNL_ASSERT_GT( val1, val2, msg ) \
+   __TNL_ASSERT_PRED2( ::TNL::Assert::cmpHelperGT, >,  val1, val2, msg )
+
+
+
+
+/****
+ * Original assert macro with custom command for diagnostic.
+ */
+
 // __CUDA_ARCH__ is defined by the compiler only for code executed on GPU
 #ifdef __CUDA_ARCH__
 #define TNL_ASSERT( ___tnl__assert_condition, ___tnl__assert_command )                                     \
@@ -39,7 +259,7 @@
            __STRING( ___tnl__assert_condition ),                                                           \
            __FILE__,                                                                                       \
            __LINE__ );                                                                                     \
-                                                                                                           \
+   asm("trap;");                                                                                           \
    }
 
 #else // __CUDA_ARCH__
@@ -48,7 +268,7 @@
    {                                                                                                    \
    std::cerr << "Assertion '" << __STRING( ___tnl__assert_condition ) << "' failed !!!" << std::endl    \
              << "File: " << __FILE__ << std::endl                                                       \
-             << "Function: " << TNL_PRETTY_FUNCTION << std::endl                                        \
+             << "Function: " << __TNL_PRETTY_FUNCTION << std::endl                                      \
              << "Line: " << __LINE__ << std::endl                                                       \
              << "Diagnostics: ";                                                                        \
         ___tnl__assert_command;                                                                         \
@@ -57,5 +277,16 @@
 #endif // __CUDA_ARCH__
 
 #else /* #ifndef NDEBUG */
+
+// empty macros for optimized build
+#define TNL_ASSERT_TRUE( val, msg )
+#define TNL_ASSERT_FALSE( val, msg )
+#define TNL_ASSERT_EQ( val1, val2, msg )
+#define TNL_ASSERT_NE( val1, val2, msg )
+#define TNL_ASSERT_LE( val1, val2, msg )
+#define TNL_ASSERT_LT( val1, val2, msg )
+#define TNL_ASSERT_GE( val1, val2, msg )
+#define TNL_ASSERT_GT( val1, val2, msg )
 #define TNL_ASSERT( ___tnl__assert_condition, ___tnl__assert_command )
+
 #endif /* #ifndef NDEBUG */
diff --git a/src/TNL/Containers/Algorithms/ArrayOperationsCuda_impl.h b/src/TNL/Containers/Algorithms/ArrayOperationsCuda_impl.h
index 174053098c23495863f48ec3b844611e5b292315..ce2f3bedcbba060ba31afca49db17a6f0e1b40ea 100644
--- a/src/TNL/Containers/Algorithms/ArrayOperationsCuda_impl.h
+++ b/src/TNL/Containers/Algorithms/ArrayOperationsCuda_impl.h
@@ -49,7 +49,7 @@ bool
 ArrayOperations< Devices::Cuda >::
 freeMemory( Element* data )
 {
-   TNL_ASSERT( data, );
+   TNL_ASSERT_TRUE( data, "Attempted to free a nullptr." );
 #ifdef HAVE_CUDA
    TNL_CHECK_CUDA_DEVICE;
    cudaFree( data );
@@ -65,7 +65,7 @@ ArrayOperations< Devices::Cuda >::
 setMemoryElement( Element* data,
                   const Element& value )
 {
-   TNL_ASSERT( data, );
+   TNL_ASSERT_TRUE( data, "Attempted to set data through a nullptr." );
    ArrayOperations< Devices::Cuda >::setMemory( data, value, 1 );
 }
 
@@ -74,7 +74,7 @@ Element
 ArrayOperations< Devices::Cuda >::
 getMemoryElement( const Element* data )
 {
-   TNL_ASSERT( data, );
+   TNL_ASSERT_TRUE( data, "Attempted to get data through a nullptr." );
    Element result;
    ArrayOperations< Devices::Host, Devices::Cuda >::copyMemory< Element, Element, int >( &result, data, 1 );
    return result;
@@ -85,7 +85,7 @@ Element&
 ArrayOperations< Devices::Cuda >::
 getArrayElementReference( Element* data, const Index i )
 {
-   TNL_ASSERT( data, );
+   TNL_ASSERT_TRUE( data, "Attempted to access data through a nullptr." );
    return data[ i ];
 }
 
@@ -94,7 +94,7 @@ const
 Element& ArrayOperations< Devices::Cuda >::
 getArrayElementReference( const Element* data, const Index i )
 {
-   TNL_ASSERT( data, );
+   TNL_ASSERT_TRUE( data, "Attempted to access data through a nullptr." );
    return data[ i ];
 }
 
@@ -123,7 +123,7 @@ setMemory( Element* data,
            const Element& value,
            const Index size )
 {
-   TNL_ASSERT( data, );
+   TNL_ASSERT_TRUE( data, "Attempted to set data through a nullptr." );
 #ifdef HAVE_CUDA
    dim3 blockSize( 0 ), gridSize( 0 );
    blockSize. x = 256;
@@ -164,8 +164,8 @@ copyMemory( DestinationElement* destination,
             const SourceElement* source,
             const Index size )
 {
-   TNL_ASSERT( destination, );
-   TNL_ASSERT( source, );
+   TNL_ASSERT_TRUE( destination, "Attempted to copy data to a nullptr." );
+   TNL_ASSERT_TRUE( source, "Attempted to copy data from a nullptr." );
 #ifdef HAVE_CUDA
    if( std::is_same< DestinationElement, SourceElement >::value )
    {
@@ -198,8 +198,8 @@ compareMemory( const Element1* destination,
                const Element2* source,
                const Index size )
 {
-   TNL_ASSERT( destination, );
-   TNL_ASSERT( source, );
+   TNL_ASSERT_TRUE( destination, "Attempted to compare data through a nullptr." );
+   TNL_ASSERT_TRUE( source, "Attempted to compare data through a nullptr." );
    //TODO: The parallel reduction on the CUDA device with different element types is needed.
    bool result;
    Algorithms::tnlParallelReductionEqualities< Element1, Index > reductionEqualities;
@@ -220,8 +220,8 @@ copyMemory( DestinationElement* destination,
             const SourceElement* source,
             const Index size )
 {
-   TNL_ASSERT( destination, );
-   TNL_ASSERT( source, );
+   TNL_ASSERT_TRUE( destination, "Attempted to copy data to a nullptr." );
+   TNL_ASSERT_TRUE( source, "Attempted to copy data from a nullptr." );
 #ifdef HAVE_CUDA
    if( std::is_same< DestinationElement, SourceElement >::value )
    {
@@ -276,9 +276,9 @@ compareMemory( const Element1* destination,
    /***
     * Here, destination is on host and source is on CUDA device.
     */
-   TNL_ASSERT( destination, );
-   TNL_ASSERT( source, );
-   TNL_ASSERT( size >= 0, std::cerr << "size = " << size );
+   TNL_ASSERT_TRUE( destination, "Attempted to compare data through a nullptr." );
+   TNL_ASSERT_TRUE( source, "Attempted to compare data through a nullptr." );
+   TNL_ASSERT_GE( size, 0, "Array size must be non-negative." );
 #ifdef HAVE_CUDA
    Element2* host_buffer = new Element2[ Devices::Cuda::getGPUTransferBufferSize() ];
    Index compared( 0 );
@@ -320,9 +320,9 @@ copyMemory( DestinationElement* destination,
             const SourceElement* source,
             const Index size )
 {
-   TNL_ASSERT( destination, );
-   TNL_ASSERT( source, );
-   TNL_ASSERT( size >= 0, std::cerr << "size = " << size );
+   TNL_ASSERT_TRUE( destination, "Attempted to copy data to a nullptr." );
+   TNL_ASSERT_TRUE( source, "Attempted to copy data from a nullptr." );
+   TNL_ASSERT_GE( size, 0, "Array size must be non-negative." );
 #ifdef HAVE_CUDA
    if( std::is_same< DestinationElement, SourceElement >::value )
    {
@@ -373,9 +373,9 @@ compareMemory( const Element1* hostData,
                const Element2* deviceData,
                const Index size )
 {
-   TNL_ASSERT( hostData, );
-   TNL_ASSERT( deviceData, );
-   TNL_ASSERT( size >= 0, std::cerr << "size = " << size );
+   TNL_ASSERT_TRUE( hostData, "Attempted to compare data through a nullptr." );
+   TNL_ASSERT_TRUE( deviceData, "Attempted to compare data through a nullptr." );
+   TNL_ASSERT_GE( size, 0, "Array size must be non-negative." );
    return ArrayOperations< Devices::Host, Devices::Cuda >::compareMemory( deviceData, hostData, size );
 }
 
diff --git a/src/TNL/Containers/Algorithms/ArrayOperationsHost_impl.h b/src/TNL/Containers/Algorithms/ArrayOperationsHost_impl.h
index f023de9e8b83b9a3eefa4a4b2143a534d42eed59..ff53ef94d0533c64d6bd633dac3f80df871496bf 100644
--- a/src/TNL/Containers/Algorithms/ArrayOperationsHost_impl.h
+++ b/src/TNL/Containers/Algorithms/ArrayOperationsHost_impl.h
@@ -30,7 +30,8 @@ allocateMemory( Element*& data,
    // According to the standard, new either throws, or returns non-nullptr.
    // Some (old) compilers don't comply:
    // https://stackoverflow.com/questions/550451/will-new-return-null-in-any-case
-   TNL_ASSERT( data, );
+   TNL_ASSERT_TRUE( data, "Operator 'new' returned a nullptr. This should never happen - there is "
+                          "either a bug or the compiler does not comply to the standard." );
    return true;
 }
 
diff --git a/src/TNL/Containers/Algorithms/Multireduction_impl.h b/src/TNL/Containers/Algorithms/Multireduction_impl.h
index b0edd1e03027448815e42b451d9450e37cf3038c..98d8c4a4b67e7a57adfc416c012ce5a0c5ef4208 100644
--- a/src/TNL/Containers/Algorithms/Multireduction_impl.h
+++ b/src/TNL/Containers/Algorithms/Multireduction_impl.h
@@ -60,8 +60,8 @@ reduce( Operation& operation,
         typename Operation::ResultType* hostResult )
 {
 #ifdef HAVE_CUDA
-   TNL_ASSERT( n > 0, );
-   TNL_ASSERT( size <= ldInput1, );
+   TNL_ASSERT_GT( n, 0, "The number of datasets must be positive." );
+   TNL_ASSERT_LE( size, ldInput1, "The size of the input cannot exceed its leading dimension." );
 
    typedef typename Operation::IndexType IndexType;
    typedef typename Operation::RealType RealType;
@@ -171,8 +171,8 @@ reduce( Operation& operation,
         const typename Operation::RealType* input2,
         typename Operation::ResultType* result )
 {
-   TNL_ASSERT( n > 0, );
-   TNL_ASSERT( size <= ldInput1, );
+   TNL_ASSERT_GT( n, 0, "The number of datasets must be positive." );
+   TNL_ASSERT_LE( size, ldInput1, "The size of the input cannot exceed its leading dimension." );
 
    typedef typename Operation::IndexType IndexType;
    typedef typename Operation::RealType RealType;
diff --git a/src/TNL/Containers/Algorithms/VectorOperationsCuda_impl.h b/src/TNL/Containers/Algorithms/VectorOperationsCuda_impl.h
index 3e3241a7a22fa1956a43331fc207e85c11a13d58..e1d5ce50b7c7cb8d4e74c9e02628bd61ccc06d18 100644
--- a/src/TNL/Containers/Algorithms/VectorOperationsCuda_impl.h
+++ b/src/TNL/Containers/Algorithms/VectorOperationsCuda_impl.h
@@ -49,7 +49,7 @@ getVectorMax( const Vector& v )
    typedef typename Vector::RealType Real;
    typedef typename Vector::IndexType Index;
 
-   TNL_ASSERT( v.getSize() > 0, );
+   TNL_ASSERT_GT( v.getSize(), 0, "Vector size must be positive." );
 
    Real result( 0 );
    Algorithms::tnlParallelReductionMax< Real, Index > operation;
@@ -69,7 +69,7 @@ getVectorMin( const Vector& v )
    typedef typename Vector::RealType Real;
    typedef typename Vector::IndexType Index;
 
-   TNL_ASSERT( v.getSize() > 0, );
+   TNL_ASSERT_GT( v.getSize(), 0, "Vector size must be positive." );
 
    Real result( 0 );
    Algorithms::tnlParallelReductionMin< Real, Index > operation;
@@ -89,7 +89,7 @@ getVectorAbsMax( const Vector& v )
    typedef typename Vector::RealType Real;
    typedef typename Vector::IndexType Index;
 
-   TNL_ASSERT( v.getSize() > 0, );
+   TNL_ASSERT_GT( v.getSize(), 0, "Vector size must be positive." );
 
    Real result( 0 );
    Algorithms::tnlParallelReductionAbsMax< Real, Index > operation;
@@ -109,7 +109,7 @@ getVectorAbsMin( const Vector& v )
    typedef typename Vector::RealType Real;
    typedef typename Vector::IndexType Index;
 
-   TNL_ASSERT( v.getSize() > 0, );
+   TNL_ASSERT_GT( v.getSize(), 0, "Vector size must be positive." );
 
    Real result( 0 );
    Algorithms::tnlParallelReductionAbsMin< Real, Index > operation;
@@ -129,7 +129,7 @@ getVectorL1Norm( const Vector& v )
    typedef typename Vector::RealType Real;
    typedef typename Vector::IndexType Index;
 
-   TNL_ASSERT( v.getSize() > 0, );
+   TNL_ASSERT_GT( v.getSize(), 0, "Vector size must be positive." );
 
    Real result( 0 );
    Algorithms::tnlParallelReductionAbsSum< Real, Index > operation;
@@ -149,7 +149,7 @@ getVectorL2Norm( const Vector& v )
    typedef typename Vector::RealType Real;
    typedef typename Vector::IndexType Index;
 
-   TNL_ASSERT( v.getSize() > 0, );
+   TNL_ASSERT_GT( v.getSize(), 0, "Vector size must be positive." );
 
    Real result( 0 );
    Algorithms::tnlParallelReductionL2Norm< Real, Index > operation;
@@ -171,9 +171,8 @@ getVectorLpNorm( const Vector& v,
    typedef typename Vector::RealType Real;
    typedef typename Vector::IndexType Index;
 
-   TNL_ASSERT( v.getSize() > 0, );
-   TNL_ASSERT( p > 0.0,
-              std::cerr << " p = " << p );
+   TNL_ASSERT_GT( v.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_GE( p, 1.0, "Parameter of the L^p norm must be at least 1.0." );
  
    if( p == 1 )
       return getVectorL1Norm( v );
@@ -198,7 +197,7 @@ getVectorSum( const Vector& v )
    typedef typename Vector::RealType Real;
    typedef typename Vector::IndexType Index;
 
-   TNL_ASSERT( v.getSize() > 0, );
+   TNL_ASSERT_GT( v.getSize(), 0, "Vector size must be positive." );
 
    Real result( 0 );
    Algorithms::tnlParallelReductionSum< Real, Index > operation;
@@ -219,8 +218,8 @@ getVectorDifferenceMax( const Vector1& v1,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( v1.getSize() > 0, );
-   TNL_ASSERT( v1.getSize() == v2.getSize(), );
+   TNL_ASSERT_GT( v1.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( v1.getSize(), v2.getSize(), "The vector sizes must be the same." );
 
    Real result( 0 );
    Algorithms::tnlParallelReductionDiffMax< Real, Index > operation;
@@ -241,8 +240,8 @@ getVectorDifferenceMin( const Vector1& v1,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( v1.getSize() > 0, );
-   TNL_ASSERT( v1.getSize() == v2.getSize(), );
+   TNL_ASSERT_GT( v1.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( v1.getSize(), v2.getSize(), "The vector sizes must be the same." );
 
    Real result( 0 );
    Algorithms::tnlParallelReductionDiffMin< Real, Index > operation;
@@ -264,8 +263,8 @@ getVectorDifferenceAbsMax( const Vector1& v1,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( v1.getSize() > 0, );
-   TNL_ASSERT( v1.getSize() == v2.getSize(), );
+   TNL_ASSERT_GT( v1.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( v1.getSize(), v2.getSize(), "The vector sizes must be the same." );
 
    Real result( 0 );
    Algorithms::tnlParallelReductionDiffAbsMax< Real, Index > operation;
@@ -286,8 +285,8 @@ getVectorDifferenceAbsMin( const Vector1& v1,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( v1.getSize() > 0, );
-   TNL_ASSERT( v1.getSize() == v2.getSize(), );
+   TNL_ASSERT_GT( v1.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( v1.getSize(), v2.getSize(), "The vector sizes must be the same." );
 
    Real result( 0 );
    Algorithms::tnlParallelReductionDiffAbsMin< Real, Index > operation;
@@ -308,8 +307,8 @@ getVectorDifferenceL1Norm( const Vector1& v1,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( v1.getSize() > 0, );
-   TNL_ASSERT( v1.getSize() == v2.getSize(), );
+   TNL_ASSERT_GT( v1.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( v1.getSize(), v2.getSize(), "The vector sizes must be the same." );
 
    Real result( 0 );
    Algorithms::tnlParallelReductionDiffAbsSum< Real, Index > operation;
@@ -330,8 +329,8 @@ getVectorDifferenceL2Norm( const Vector1& v1,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( v1.getSize() > 0, );
-   TNL_ASSERT( v1.getSize() == v2.getSize(), );
+   TNL_ASSERT_GT( v1.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( v1.getSize(), v2.getSize(), "The vector sizes must be the same." );
 
    Real result( 0 );
    Algorithms::tnlParallelReductionDiffL2Norm< Real, Index > operation;
@@ -354,10 +353,9 @@ getVectorDifferenceLpNorm( const Vector1& v1,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( p > 0.0,
-              std::cerr << " p = " << p );
-   TNL_ASSERT( v1.getSize() > 0, );
-   TNL_ASSERT( v1.getSize() == v2.getSize(), );
+   TNL_ASSERT_GT( v1.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( v1.getSize(), v2.getSize(), "The vector sizes must be the same." );
+   TNL_ASSERT_GE( p, 1.0, "Parameter of the L^p norm must be at least 1.0." );
 
    Real result( 0 );
    Algorithms::tnlParallelReductionDiffLpNorm< Real, Index > operation;
@@ -379,8 +377,8 @@ getVectorDifferenceSum( const Vector1& v1,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( v1.getSize() > 0, );
-   TNL_ASSERT( v1.getSize() == v2.getSize(), );
+   TNL_ASSERT_GT( v1.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( v1.getSize(), v2.getSize(), "The vector sizes must be the same." );
 
    Real result( 0 );
    Algorithms::tnlParallelReductionDiffSum< Real, Index > operation;
@@ -418,7 +416,7 @@ vectorScalarMultiplication( Vector& v,
    typedef typename Vector::RealType Real;
    typedef typename Vector::IndexType Index;
 
-   TNL_ASSERT( v.getSize() > 0, );
+   TNL_ASSERT_GT( v.getSize(), 0, "Vector size must be positive." );
 
 #ifdef HAVE_CUDA
    dim3 blockSize( 0 ), gridSize( 0 );
@@ -445,8 +443,8 @@ getScalarProduct( const Vector1& v1,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( v1.getSize() > 0, );
-   TNL_ASSERT( v1.getSize() == v2.getSize(), );
+   TNL_ASSERT_GT( v1.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( v1.getSize(), v2.getSize(), "The vector sizes must be the same." );
 
    Real result( 0 );
 /*#if defined HAVE_CUBLAS && defined HAVE_CUDA
@@ -502,10 +500,8 @@ addVector( Vector1& y,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( y.getSize() > 0, );
-   TNL_ASSERT( y.getSize() == x.getSize(), );
-   TNL_ASSERT( y.getData() != 0, );
-   TNL_ASSERT( x.getData() != 0, );
+   TNL_ASSERT_GT( x.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( x.getSize(), y.getSize(), "The vector sizes must be the same." );
 
 #ifdef HAVE_CUDA
    dim3 blockSize( 0 ), gridSize( 0 );
@@ -573,12 +569,9 @@ addVectors( Vector1& v,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( v.getSize() > 0, );
-   TNL_ASSERT( v.getSize() == v1.getSize(), );
-   TNL_ASSERT( v.getSize() == v2.getSize(), );
-   TNL_ASSERT( v.getData() != 0, );
-   TNL_ASSERT( v1.getData() != 0, );
-   TNL_ASSERT( v2.getData() != 0, );
+   TNL_ASSERT_GT( v.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( v.getSize(), v1.getSize(), "The vector sizes must be the same." );
+   TNL_ASSERT_EQ( v.getSize(), v2.getSize(), "The vector sizes must be the same." );
 
 #ifdef HAVE_CUDA
    dim3 blockSize( 0 ), gridSize( 0 );
diff --git a/src/TNL/Containers/Algorithms/VectorOperationsHost_impl.h b/src/TNL/Containers/Algorithms/VectorOperationsHost_impl.h
index aab2a84769c1bf745351ec1e141df1eb740f04ce..492153c729e23a6022f441aee75ddcb0f4a5dabf 100644
--- a/src/TNL/Containers/Algorithms/VectorOperationsHost_impl.h
+++ b/src/TNL/Containers/Algorithms/VectorOperationsHost_impl.h
@@ -48,7 +48,7 @@ getVectorMax( const Vector& v )
    typedef typename Vector::RealType Real;
    typedef typename Vector::IndexType Index;
 
-   TNL_ASSERT( v.getSize() > 0, );
+   TNL_ASSERT_GT( v.getSize(), 0, "Vector size must be positive." );
 
    Real result = v.getElement( 0 );
    const Index n = v.getSize();
@@ -68,7 +68,7 @@ getVectorMin( const Vector& v )
    typedef typename Vector::RealType Real;
    typedef typename Vector::IndexType Index;
 
-   TNL_ASSERT( v.getSize() > 0, );
+   TNL_ASSERT_GT( v.getSize(), 0, "Vector size must be positive." );
 
    Real result = v.getElement( 0 );
    const Index n = v.getSize();
@@ -88,7 +88,7 @@ getVectorAbsMax( const Vector& v )
    typedef typename Vector::RealType Real;
    typedef typename Vector::IndexType Index;
 
-   TNL_ASSERT( v.getSize() > 0, );
+   TNL_ASSERT_GT( v.getSize(), 0, "Vector size must be positive." );
 
    Real result = std::fabs( v.getElement( 0 ) );
    const Index n = v.getSize();
@@ -109,7 +109,7 @@ getVectorAbsMin( const Vector& v )
    typedef typename Vector::RealType Real;
    typedef typename Vector::IndexType Index;
 
-   TNL_ASSERT( v.getSize() > 0, );
+   TNL_ASSERT_GT( v.getSize(), 0, "Vector size must be positive." );
 
    Real result = std::fabs( v.getElement( 0 ) );
    const Index n = v.getSize();
@@ -129,7 +129,7 @@ getVectorL1Norm( const Vector& v )
    typedef typename Vector::RealType Real;
    typedef typename Vector::IndexType Index;
 
-   TNL_ASSERT( v.getSize() > 0, );
+   TNL_ASSERT_GT( v.getSize(), 0, "Vector size must be positive." );
 
    Real result( 0.0 );
    const Index n = v.getSize();
@@ -149,7 +149,7 @@ getVectorL2Norm( const Vector& v )
    typedef typename Vector::RealType Real;
    typedef typename Vector::IndexType Index;
 
-   TNL_ASSERT( v.getSize() > 0, );
+   TNL_ASSERT_GT( v.getSize(), 0, "Vector size must be positive." );
 
    const Index n = v.getSize();
 
@@ -212,9 +212,8 @@ getVectorLpNorm( const Vector& v,
    typedef typename Vector::RealType Real;
    typedef typename Vector::IndexType Index;
 
-   TNL_ASSERT( v.getSize() > 0, );
-   TNL_ASSERT( p > 0.0,
-              std::cerr << " p = " << p );
+   TNL_ASSERT_GT( v.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_GE( p, 1.0, "Parameter of the L^p norm must be at least 1.0." );
 
    if( p == 1.0 )
       return getVectorL1Norm( v );
@@ -239,7 +238,7 @@ getVectorSum( const Vector& v )
    typedef typename Vector::RealType Real;
    typedef typename Vector::IndexType Index;
 
-   TNL_ASSERT( v.getSize() > 0, );
+   TNL_ASSERT_GT( v.getSize(), 0, "Vector size must be positive." );
 
    Real result( 0.0 );
    const Index n = v.getSize();
@@ -260,8 +259,8 @@ getVectorDifferenceMax( const Vector1& v1,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( v1.getSize() > 0, );
-   TNL_ASSERT( v1.getSize() == v2.getSize(), );
+   TNL_ASSERT_GT( v1.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( v1.getSize(), v2.getSize(), "The vector sizes must be the same." );
 
    Real result = v1.getElement( 0 ) - v2.getElement( 0 );
    const Index n = v1.getSize();
@@ -282,8 +281,8 @@ getVectorDifferenceMin( const Vector1& v1,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( v1.getSize() > 0, );
-   TNL_ASSERT( v1.getSize() == v2.getSize(), );
+   TNL_ASSERT_GT( v1.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( v1.getSize(), v2.getSize(), "The vector sizes must be the same." );
 
    Real result = v1.getElement( 0 ) - v2.getElement( 0 );
    const Index n = v1.getSize();
@@ -304,8 +303,8 @@ getVectorDifferenceAbsMax( const Vector1& v1,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( v1.getSize() > 0, );
-   TNL_ASSERT( v1.getSize() == v2.getSize(), );
+   TNL_ASSERT_GT( v1.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( v1.getSize(), v2.getSize(), "The vector sizes must be the same." );
 
    Real result = std::fabs( v1.getElement( 0 ) - v2.getElement( 0 ) );
    const Index n = v1.getSize();
@@ -326,8 +325,8 @@ getVectorDifferenceAbsMin( const Vector1& v1,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( v1.getSize() > 0, );
-   TNL_ASSERT( v1.getSize() == v2.getSize(), );
+   TNL_ASSERT_GT( v1.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( v1.getSize(), v2.getSize(), "The vector sizes must be the same." );
 
    Real result = std::fabs( v1[ 0 ] - v2[ 0 ] );
    const Index n = v1.getSize();
@@ -348,8 +347,8 @@ getVectorDifferenceL1Norm( const Vector1& v1,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( v1.getSize() > 0, );
-   TNL_ASSERT( v1.getSize() == v2.getSize(), );
+   TNL_ASSERT_GT( v1.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( v1.getSize(), v2.getSize(), "The vector sizes must be the same." );
 
    Real result( 0.0 );
    const Index n = v1.getSize();
@@ -370,8 +369,8 @@ getVectorDifferenceL2Norm( const Vector1& v1,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( v1.getSize() > 0, );
-   TNL_ASSERT( v1.getSize() == v2.getSize(), );
+   TNL_ASSERT_GT( v1.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( v1.getSize(), v2.getSize(), "The vector sizes must be the same." );
 
    Real result( 0.0 );
    const Index n = v1.getSize();
@@ -397,10 +396,9 @@ getVectorDifferenceLpNorm( const Vector1& v1,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( p > 0.0,
-              std::cerr << " p = " << p );
-   TNL_ASSERT( v1.getSize() > 0, );
-   TNL_ASSERT( v1.getSize() == v2.getSize(), );
+   TNL_ASSERT_GT( v1.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( v1.getSize(), v2.getSize(), "The vector sizes must be the same." );
+   TNL_ASSERT_GE( p, 1.0, "Parameter of the L^p norm must be at least 1.0." );
 
    if( p == 1.0 )
       return getVectorDifferenceL1Norm( v1, v2 );
@@ -426,8 +424,8 @@ getVectorDifferenceSum( const Vector1& v1,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( v1.getSize() > 0, );
-   TNL_ASSERT( v1.getSize() == v2.getSize(), );
+   TNL_ASSERT_GT( v1.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( v1.getSize(), v2.getSize(), "The vector sizes must be the same." );
 
    Real result( 0.0 );
    const Index n = v1.getSize();
@@ -449,7 +447,7 @@ vectorScalarMultiplication( Vector& v,
    typedef typename Vector::RealType Real;
    typedef typename Vector::IndexType Index;
 
-   TNL_ASSERT( v.getSize() > 0, );
+   TNL_ASSERT_GT( v.getSize(), 0, "Vector size must be positive." );
 
    const Index n = v.getSize();
 #ifdef HAVE_OPENMP
@@ -469,8 +467,8 @@ getScalarProduct( const Vector1& v1,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( v1.getSize() > 0, );
-   TNL_ASSERT( v1.getSize() == v2.getSize(), );
+   TNL_ASSERT_GT( v1.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( v1.getSize(), v2.getSize(), "The vector sizes must be the same." );
    const Index n = v1.getSize();
 
 #ifdef OPTIMIZED_VECTOR_HOST_OPERATIONS
@@ -533,8 +531,8 @@ addVector( Vector1& y,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( x.getSize() > 0, );
-   TNL_ASSERT( x.getSize() == y.getSize(), );
+   TNL_ASSERT_GT( x.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( x.getSize(), y.getSize(), "The vector sizes must be the same." );
 
    const Index n = y.getSize();
 
@@ -603,9 +601,9 @@ addVectors( Vector1& v,
    typedef typename Vector1::RealType Real;
    typedef typename Vector1::IndexType Index;
 
-   TNL_ASSERT( v.getSize() > 0, );
-   TNL_ASSERT( v.getSize() == v1.getSize(), );
-   TNL_ASSERT( v.getSize() == v2.getSize(), );
+   TNL_ASSERT_GT( v.getSize(), 0, "Vector size must be positive." );
+   TNL_ASSERT_EQ( v.getSize(), v1.getSize(), "The vector sizes must be the same." );
+   TNL_ASSERT_EQ( v.getSize(), v2.getSize(), "The vector sizes must be the same." );
  
    const Index n = v.getSize();
    if( thisMultiplicator == 1.0 )
diff --git a/src/TNL/Containers/Array_impl.h b/src/TNL/Containers/Array_impl.h
index e1d8dd351f75d32843f7e71b8ae22b78d99606b3..674b32f0c07d7f981ed8fac213a8cca27151fab6 100644
--- a/src/TNL/Containers/Array_impl.h
+++ b/src/TNL/Containers/Array_impl.h
@@ -74,10 +74,9 @@ Array( Array< Element, Device, Index >& array,
   allocationPointer( array.allocationPointer ),
   referenceCounter( 0 )
 {
-   TNL_ASSERT( begin < array.getSize(),
-              std::cerr << " begin = " << begin << " array.getSize() = " << array.getSize() );
-   TNL_ASSERT( begin + size  < array.getSize(),
-              std::cerr << " begin = " << begin << " size = " << size <<  " array.getSize() = " << array.getSize() );
+   TNL_ASSERT_LT( begin, array.getSize(), "Begin of array is out of bounds." );
+   TNL_ASSERT_LE( begin + size, array.getSize(), "End of array is out of bounds." );
+
    if( ! this->size )
       this->size = array.getSize() - begin;
    if( array.allocationPointer )
@@ -169,9 +168,8 @@ void
 Array< Element, Device, Index >::
 setSize( const Index size )
 {
-   TNL_ASSERT( size >= 0,
-              std::cerr << "You try to set size of Array to negative value."
-                        << "New size: " << size << std::endl );
+   TNL_ASSERT_GE( size, 0, "Array size must be non-negative." );
+
    if( this->size == size && allocationPointer && ! referenceCounter )
       return;
    this->releaseData();
@@ -182,8 +180,8 @@ setSize( const Index size )
       Algorithms::ArrayOperations< Device >::allocateMemory( this->allocationPointer, size );
       this->data = this->allocationPointer;
       this->size = size;
-      TNL_ASSERT( this->allocationPointer,
-                  std::cerr << "This should never happen - allocator did not throw on an error." << std::endl; );
+      TNL_ASSERT_TRUE( this->allocationPointer,
+                       "This should never happen - allocator did not throw on an error." );
    }
 }
 
@@ -225,10 +223,8 @@ bind( const ArrayT& array,
    static_assert( std::is_same< Element, typename ArrayT::ElementType >::value, "ElementType of both arrays must be the same." );
    static_assert( std::is_same< Device, typename ArrayT::DeviceType >::value, "DeviceType of both arrays must be the same." );
    static_assert( std::is_same< Index, typename ArrayT::IndexType >::value, "IndexType of both arrays must be the same." );
-   TNL_ASSERT( begin <= array.getSize(),
-              std::cerr << " begin = " << begin << " array.getSize() = " << array.getSize() );
-   TNL_ASSERT( begin + size  <= array.getSize(),
-              std::cerr << " begin = " << begin << " size = " << size <<  " array.getSize() = " << array.getSize() );
+   TNL_ASSERT_LT( begin, array.getSize(), "Begin of array is out of bounds." );
+   TNL_ASSERT_LE( begin + size, array.getSize(), "End of array is out of bounds." );
  
    this->releaseData();
    if( size )
@@ -307,10 +303,8 @@ void
 Array< Element, Device, Index >::
 setElement( const Index& i, const Element& x )
 {
-   TNL_ASSERT( 0 <= i && i < this->getSize(),
-              std::cerr << "Wrong index for setElement method in Array "
-                        << " index is " << i
-                        << " and array size is " << this->getSize() );
+   TNL_ASSERT_GE( i, 0, "Element index must be non-negative." );
+   TNL_ASSERT_LT( i, this->getSize(), "Element index is out of bounds." );
    return Algorithms::ArrayOperations< Device > :: setMemoryElement( &( this->data[ i ] ), x );
 }
 
@@ -321,10 +315,8 @@ Element
 Array< Element, Device, Index >::
 getElement( const Index& i ) const
 {
-   TNL_ASSERT( 0 <= i && i < this->getSize(),
-              std::cerr << "Wrong index for getElement method in Array "
-                        << " index is " << i
-                        << " and array size is " << this->getSize() );
+   TNL_ASSERT_GE( i, 0, "Element index must be non-negative." );
+   TNL_ASSERT_LT( i, this->getSize(), "Element index is out of bounds." );
    return Algorithms::ArrayOperations< Device >::getMemoryElement( & ( this->data[ i ] ) );
 }
 
@@ -336,10 +328,8 @@ inline Element&
 Array< Element, Device, Index >::
 operator[] ( const Index& i )
 {
-   TNL_ASSERT( 0 <= i && i < this->getSize(),
-              std::cerr << "Wrong index for operator[] in Array "
-                        << " index is " << i
-                        << " and array size is " << this->getSize() );
+   TNL_ASSERT_GE( i, 0, "Element index must be non-negative." );
+   TNL_ASSERT_LT( i, this->getSize(), "Element index is out of bounds." );
    return this->data[ i ];
 }
 
@@ -351,10 +341,8 @@ inline const Element&
 Array< Element, Device, Index >::
 operator[] ( const Index& i ) const
 {
-   TNL_ASSERT( 0 <= i && i < this->getSize(),
-              std::cerr << "Wrong index for operator[] in Array "
-                        << " index is " << i
-                        << " and array size is " << this->getSize() );
+   TNL_ASSERT_GE( i, 0, "Element index must be non-negative." );
+   TNL_ASSERT_LT( i, this->getSize(), "Element index is out of bounds." );
    return this->data[ i ];
 }
 
@@ -365,9 +353,7 @@ Array< Element, Device, Index >&
 Array< Element, Device, Index >::
 operator = ( const Array< Element, Device, Index >& array )
 {
-   TNL_ASSERT( array. getSize() == this->getSize(),
-              std::cerr << "Source size: " << array. getSize() << std::endl
-                        << "Target size: " << this->getSize() << std::endl );
+   TNL_ASSERT_EQ( array.getSize(), this->getSize(), "Array sizes must be the same." );
    if( this->getSize() > 0 )
       Algorithms::ArrayOperations< Device >::
          template copyMemory< Element,
@@ -387,9 +373,7 @@ Array< Element, Device, Index >&
 Array< Element, Device, Index >::
 operator = ( const ArrayT& array )
 {
-   TNL_ASSERT( array. getSize() == this->getSize(),
-              std::cerr << "Source size: " << array. getSize() << std::endl
-                        << "Target size: " << this->getSize() << std::endl );
+   TNL_ASSERT_EQ( array.getSize(), this->getSize(), "Array sizes must be the same." );
    if( this->getSize() > 0 )
       Algorithms::ArrayOperations< Device, typename ArrayT::DeviceType >::
          template copyMemory< Element,
@@ -409,7 +393,7 @@ bool
 Array< Element, Device, Index >::
 operator == ( const ArrayT& array ) const
 {
-   if( array. getSize() != this -> getSize() )
+   if( array.getSize() != this->getSize() )
       return false;
    if( this->getSize() == 0 )
       return true;
@@ -437,7 +421,7 @@ template< typename Element,
           typename Index >
 void Array< Element, Device, Index > :: setValue( const Element& e )
 {
-   TNL_ASSERT( this->getData(),);
+   TNL_ASSERT_TRUE( this->getData(), "Attempted to set a value of an empty array." );
    Algorithms::ArrayOperations< Device >::setMemory( this->getData(), e, this->getSize() );
 }
 
diff --git a/src/TNL/Containers/StaticArray1D_impl.h b/src/TNL/Containers/StaticArray1D_impl.h
index cfbabafac53aaab24b32cadefb1ae85ce7fd0347..22146311f77309e4817f4a1e6e970a279f4ad70c 100644
--- a/src/TNL/Containers/StaticArray1D_impl.h
+++ b/src/TNL/Containers/StaticArray1D_impl.h
@@ -77,8 +77,8 @@ template< typename Element >
 __cuda_callable__
 inline const Element& StaticArray< 1, Element >::operator[]( int i ) const
 {
-   TNL_ASSERT( i >= 0 && i < size,
-            std::cerr << "i = " << i << " size = " << size << std::endl; );
+   TNL_ASSERT_GE( i, 0, "Element index must be non-negative." );
+   TNL_ASSERT_LT( i, size, "Element index is out of bounds." );
    return data[ i ];
 }
 
@@ -86,8 +86,8 @@ template< typename Element >
 __cuda_callable__
 inline Element& StaticArray< 1, Element >::operator[]( int i )
 {
-   TNL_ASSERT( i >= 0 && i < size,
-            std::cerr << "i = " << i << " size = " << size << std::endl; );
+   TNL_ASSERT_GE( i, 0, "Element index must be non-negative." );
+   TNL_ASSERT_LT( i, size, "Element index is out of bounds." );
    return data[ i ];
 }
 
diff --git a/src/TNL/Containers/StaticArray2D_impl.h b/src/TNL/Containers/StaticArray2D_impl.h
index a0301a4bbf13f19abe1f78577de2c45cc2ece4c7..b780903e3e7586bdd78a96a37dcec4b70ee376e8 100644
--- a/src/TNL/Containers/StaticArray2D_impl.h
+++ b/src/TNL/Containers/StaticArray2D_impl.h
@@ -89,8 +89,8 @@ template< typename Element >
 __cuda_callable__
 inline const Element& StaticArray< 2, Element >::operator[]( int i ) const
 {
-   TNL_ASSERT( i >= 0 && i < size,
-            std::cerr << "i = " << i << " size = " << size << std::endl; );
+   TNL_ASSERT_GE( i, 0, "Element index must be non-negative." );
+   TNL_ASSERT_LT( i, size, "Element index is out of bounds." );
    return data[ i ];
 }
 
@@ -98,8 +98,8 @@ template< typename Element >
 __cuda_callable__
 inline Element& StaticArray< 2, Element >::operator[]( int i )
 {
-   TNL_ASSERT( i >= 0 && i < size,
-            std::cerr << "i = " << i << " size = " << size << std::endl; );
+   TNL_ASSERT_GE( i, 0, "Element index must be non-negative." );
+   TNL_ASSERT_LT( i, size, "Element index is out of bounds." );
    return data[ i ];
 }
 
diff --git a/src/TNL/Containers/StaticArray3D_impl.h b/src/TNL/Containers/StaticArray3D_impl.h
index a246b2e3e6e2e1a26bca5c7c922a0b0371fb3dc1..7eea0909e7f8e771d391bc34b96e104775127024 100644
--- a/src/TNL/Containers/StaticArray3D_impl.h
+++ b/src/TNL/Containers/StaticArray3D_impl.h
@@ -92,8 +92,8 @@ template< typename Element >
 __cuda_callable__
 inline const Element& StaticArray< 3, Element >::operator[]( int i ) const
 {
-   TNL_ASSERT( i >= 0 && i < size,
-            std::cerr << "i = " << i << " size = " << size << std::endl; );
+   TNL_ASSERT_GE( i, 0, "Element index must be non-negative." );
+   TNL_ASSERT_LT( i, size, "Element index is out of bounds." );
    return data[ i ];
 }
 
@@ -101,8 +101,8 @@ template< typename Element >
 __cuda_callable__
 inline Element& StaticArray< 3, Element >::operator[]( int i )
 {
-   TNL_ASSERT( i >= 0 && i < size,
-            std::cerr << "i = " << i << " size = " << size << std::endl; );
+   TNL_ASSERT_GE( i, 0, "Element index must be non-negative." );
+   TNL_ASSERT_LT( i, size, "Element index is out of bounds." );
    return data[ i ];
 }
 
diff --git a/src/TNL/Containers/StaticArray_impl.h b/src/TNL/Containers/StaticArray_impl.h
index 8e639351abffeadd605031efe69cc0c8092dbb0e..60e9257dea5526bb1cd58a74a28d641ff79d9de8 100644
--- a/src/TNL/Containers/StaticArray_impl.h
+++ b/src/TNL/Containers/StaticArray_impl.h
@@ -81,8 +81,8 @@ template< int Size, typename Element >
 __cuda_callable__
 inline const Element& StaticArray< Size, Element >::operator[]( int i ) const
 {
-   TNL_ASSERT( i >= 0 && i < size,
-            std::cerr << "i = " << i << " size = " << size << std::endl; );
+   TNL_ASSERT_GE( i, 0, "Element index must be non-negative." );
+   TNL_ASSERT_LT( i, size, "Element index is out of bounds." );
    return data[ i ];
 }
 
@@ -90,8 +90,8 @@ template< int Size, typename Element >
 __cuda_callable__
 inline Element& StaticArray< Size, Element >::operator[]( int i )
 {
-   TNL_ASSERT( i >= 0 && i < size,
-            std::cerr << "i = " << i << " size = " << size << std::endl; );
+   TNL_ASSERT_GE( i, 0, "Element index must be non-negative." );
+   TNL_ASSERT_LT( i, size, "Element index is out of bounds." );
    return data[ i ];
 }
 
diff --git a/src/TNL/File_impl.h b/src/TNL/File_impl.h
index 5565b683869bb6f6fc3d0e6483e87b62489f0faf..67609929cb0da39343e85678765f3354c42cd98f 100644
--- a/src/TNL/File_impl.h
+++ b/src/TNL/File_impl.h
@@ -10,6 +10,7 @@
 
 #pragma once 
 
+#include <TNL/File.h>
 #include <TNL/Exceptions/CudaSupportMissing.h>
 
 namespace TNL {
@@ -31,8 +32,7 @@ template< typename Type, typename Device, typename Index >
 bool File :: read( Type* buffer,
                    const Index& _elements )
 {
-   TNL_ASSERT( _elements >= 0,
-           std::cerr << " elements = " << _elements << std::endl; );
+   TNL_ASSERT_GE( _elements, 0, "Number of elements to read must be non-negative." );
 
    // convert _elements from Index to size_t, which is *unsigned* type
    // (expected by fread etc)
@@ -127,8 +127,7 @@ template< class Type, typename Device, typename Index >
 bool File :: write( const Type* buffer,
                     const Index _elements )
 {
-   TNL_ASSERT( _elements >= 0,
-           std::cerr << " elements = " << _elements << std::endl; );
+   TNL_ASSERT_GE( _elements, 0, "Number of elements to write must be non-negative." );
 
    // convert _elements from Index to size_t, which is *unsigned* type
    // (expected by fread etc)
@@ -220,5 +219,3 @@ bool File :: write( const Type* buffer,
 };
 
 } // namespace TNL
-
-
diff --git a/src/TNL/Functions/Analytic/SinWaveSDF_impl.h b/src/TNL/Functions/Analytic/SinWaveSDF_impl.h
index d7ac240901502c9d3fcfb391da969954dc40609d..d104ca32210896e9103dbd4205670da4202c1940 100644
--- a/src/TNL/Functions/Analytic/SinWaveSDF_impl.h
+++ b/src/TNL/Functions/Analytic/SinWaveSDF_impl.h
@@ -13,8 +13,8 @@
 #include <TNL/Functions/Analytic/SinWaveSDF.h>
 
 namespace TNL {
-   namespace Functions {
-      namespace Analytic {
+namespace Functions {
+namespace Analytic {
 
 template< int dimensions, typename Real >
 SinWaveSDFBase< dimensions, Real >::SinWaveSDFBase()
@@ -115,7 +115,7 @@ getPartialDerivative( const PointType& v,
    const RealType distance = ::sqrt( x * x ) + this->phase * this->waveLength / (2.0*M_PI);
    if( XDiffOrder == 0 )
       return this->sinWaveFunctionSDF( distance );
-   TNL_ASSERT( false, std::cerr << "TODO: implement this" );
+   TNL_ASSERT_TRUE( false, "TODO: implement this" );
    return 0.0;
 }
 
@@ -138,7 +138,7 @@ getPartialDerivative( const PointType& v,
    const RealType distance  = ::sqrt( x * x + y * y ) + this->phase * this->waveLength / (2.0*M_PI);
    if( XDiffOrder == 0 && YDiffOrder == 0)
       return this->sinWaveFunctionSDF( distance );
-   TNL_ASSERT( false, std::cerr << "TODO: implement this" );
+   TNL_ASSERT_TRUE( false, "TODO: implement this" );
    return 0.0;
 }
 
@@ -158,10 +158,10 @@ getPartialDerivative( const PointType& v,
    const RealType distance  = ::sqrt( x * x +  y * y + z * z ) +  this->phase * this->waveLength / (2.0*M_PI);
    if( XDiffOrder == 0 && YDiffOrder == 0 && ZDiffOrder == 0 )
       return this->sinWaveFunctionSDF( distance );
-   TNL_ASSERT( false, std::cerr << "TODO: implement this" );
+   TNL_ASSERT_TRUE( false, "TODO: implement this" );
    return 0.0;
 }
 
-      } // namespace Analytic
-   } // namespace Functions
+} // namespace Analytic
+} // namespace Functions
 } // namespace TNL
diff --git a/src/TNL/Functions/Analytic/VectorNorm.h b/src/TNL/Functions/Analytic/VectorNorm.h
index 447b31973815ddbdee4ff99c09c518d7fd75b753..b54513705980249db0037044b00dabeca045bce6 100644
--- a/src/TNL/Functions/Analytic/VectorNorm.h
+++ b/src/TNL/Functions/Analytic/VectorNorm.h
@@ -211,7 +211,7 @@ class VectorNorm< 2, Real > : public VectorNormBase< 2, Real >
             return ( std::pow( std::pow( TNL::abs( x ), this->power ) * this->anisotropy.x() + 
                                std::pow( TNL::abs( y ), this->power ) * this->anisotropy.y(), 1.0 / this-> power ) - this->radius ) * this->multiplicator;
          }
-         TNL_ASSERT( false, std::cerr << "Not implemented yet." << std::endl );
+         TNL_ASSERT_TRUE( false, "Not implemented yet." );
          return 0.0;
       }
  
@@ -262,7 +262,7 @@ class VectorNorm< 3, Real > : public VectorNormBase< 3, Real >
                                std::pow( TNL::abs( y ), this->power ) * this->anisotropy.y() +
                                std::pow( TNL::abs( z ), this->power ) * this->anisotropy.z(), 1.0 / this-> power ) - this->radius ) * this->multiplicator;
          }
-         TNL_ASSERT( false, std::cerr << "Not implemented yet." << std::endl );
+         TNL_ASSERT_TRUE( false, "Not implemented yet." );
          return 0.0;
       }
       
diff --git a/src/TNL/Functions/MeshFunctionNormGetter.h b/src/TNL/Functions/MeshFunctionNormGetter.h
index 6f7e127c173b0dc0b67f51f800ced2f9d333cc6a..50e39c6de767fa79828d7e4bc9c488356c1d5be6 100644
--- a/src/TNL/Functions/MeshFunctionNormGetter.h
+++ b/src/TNL/Functions/MeshFunctionNormGetter.h
@@ -134,7 +134,7 @@ class MeshFunctionNormGetter< MeshFunction< Meshes::Grid< Dimension, MeshReal, D
          }
          if( EntityDimension > 0 )
          {
-            TNL_ASSERT( false, std::cerr << "Not implemented yet." << std::endl );
+            TNL_ASSERT_TRUE( false, "Not implemented yet." );
          }
  
          if( p == 1.0 )
diff --git a/src/TNL/Functions/OperatorFunction.h b/src/TNL/Functions/OperatorFunction.h
index 4ae8907daaee25d9c8d3305e0c447f44691f1420..c5dbde42363c27f2d58d53018d7efb41096b76ad 100644
--- a/src/TNL/Functions/OperatorFunction.h
+++ b/src/TNL/Functions/OperatorFunction.h
@@ -87,13 +87,13 @@ class OperatorFunction< Operator, MeshFunctionT, void, true, IsAnalytic >
  
       const MeshType& getMesh() const
       {
-         TNL_ASSERT( this->preimageFunction, std::cerr << "The preimage function was not set." << std::endl );
+         TNL_ASSERT_TRUE( this->preimageFunction, "The preimage function was not set." );
          return this->preimageFunction->getMesh();
       };
       
       const MeshPointer& getMeshPointer() const
       { 
-         TNL_ASSERT( this->preimageFunction, std::cerr << "The preimage function was not set." << std::endl );
+         TNL_ASSERT_TRUE( this->preimageFunction, "The preimage function was not set." );
          return this->preimageFunction->getMeshPointer(); 
       };
 
@@ -114,7 +114,7 @@ class OperatorFunction< Operator, MeshFunctionT, void, true, IsAnalytic >
          const MeshEntity& meshEntity,
          const RealType& time = 0.0 ) const
       {
-         TNL_ASSERT( this->preimageFunction, std::cerr << "The preimage function was not set." << std::endl );
+         TNL_ASSERT_TRUE( this->preimageFunction, "The preimage function was not set." );
          return operator_( *preimageFunction, meshEntity, time );
       }
  
@@ -298,13 +298,13 @@ class OperatorFunction< Operator, Function, BoundaryConditions, false, IsAnalyti
  
       const PreimageFunctionType& getPreimageFunction() const
       {
-         TNL_ASSERT( this->preimageFunction, );
+         TNL_ASSERT_TRUE( this->preimageFunction, "The preimage function was not set." );
          return *this->preimageFunction;
       };
  
       PreimageFunctionType& getPreimageFunction()
       {
-         TNL_ASSERT( this->preimageFunction, );
+         TNL_ASSERT_TRUE( this->preimageFunction, "The preimage function was not set." );
          return *this->preimageFunction;
       };
  
diff --git a/src/TNL/Meshes/GridDetails/Grid1D_impl.h b/src/TNL/Meshes/GridDetails/Grid1D_impl.h
index 3c388fff332f6a8717f4e8e063e1c16e63a5da1e..6fe0a1557a16ead07740dc543aac56fd0f12baef 100644
--- a/src/TNL/Meshes/GridDetails/Grid1D_impl.h
+++ b/src/TNL/Meshes/GridDetails/Grid1D_impl.h
@@ -14,6 +14,7 @@
 #include <iomanip>
 #include <TNL/String.h>
 #include <TNL/Assert.h>
+#include <TNL/Logger.h>
 #include <TNL/Meshes/GridDetails/GnuplotWriter.h>
 #include <TNL/Meshes/GridDetails/GridEntityGetter_impl.h>
 #include <TNL/Meshes/GridDetails/NeighbourGridEntityGetter1D_impl.h>
@@ -90,7 +91,7 @@ template< typename Real,
           typename Index  >
 void Grid< 1, Real, Device, Index >::setDimensions( const Index xSize )
 {
-   TNL_ASSERT( xSize > 0, std::cerr << "xSize = " << xSize );
+   TNL_ASSERT_GT( xSize, 0, "Grid size must be positive." );
    this->dimensions.x() = xSize;
    this->numberOfCells = xSize;
    this->numberOfVertices = xSize + 1;
@@ -231,8 +232,7 @@ const Real&
 Grid< 1, Real, Device, Index >::
 getSpaceStepsProducts() const
 {
-   TNL_ASSERT( xPow >= -2 && xPow <= 2,
-              std::cerr << " xPow = " << xPow );
+   static_assert( xPow >= -2 && xPow <= 2, "unsupported value of xPow" );
    return this->spaceStepsProducts[ xPow + 2 ];
 }
 
diff --git a/src/TNL/Meshes/GridDetails/Grid2D_impl.h b/src/TNL/Meshes/GridDetails/Grid2D_impl.h
index 90a1ae99bbc931afd6aa38ac4cb838b0ce4b18a5..459748e93c5c7cb5189eb4de26e28c5d8c29a0bc 100644
--- a/src/TNL/Meshes/GridDetails/Grid2D_impl.h
+++ b/src/TNL/Meshes/GridDetails/Grid2D_impl.h
@@ -12,7 +12,9 @@
 
 #include <fstream>
 #include <iomanip>
+#include <TNL/String.h>
 #include <TNL/Assert.h>
+#include <TNL/Logger.h>
 #include <TNL/Meshes/GridDetails/GnuplotWriter.h>
 #include <TNL/Meshes/GridDetails/GridEntityGetter_impl.h>
 #include <TNL/Meshes/GridDetails/NeighbourGridEntityGetter2D_impl.h>
@@ -135,8 +137,8 @@ template< typename Real,
           typename Index >
 void Grid< 2, Real, Device, Index > :: setDimensions( const Index xSize, const Index ySize )
 {
-   TNL_ASSERT( xSize > 0, std::cerr << "xSize = " << xSize );
-   TNL_ASSERT( ySize > 0, std::cerr << "ySize = " << ySize );
+   TNL_ASSERT_GT( xSize, 0, "Grid size must be positive." );
+   TNL_ASSERT_GT( ySize, 0, "Grid size must be positive." );
 
    this->dimensions.x() = xSize;
    this->dimensions.y() = ySize;
@@ -284,11 +286,8 @@ const Real&
 Grid< 2, Real, Device, Index >::
 getSpaceStepsProducts() const
 {
-   TNL_ASSERT( xPow >= -2 && xPow <= 2,
-              std::cerr << " xPow = " << xPow );
-   TNL_ASSERT( yPow >= -2 && yPow <= 2,
-              std::cerr << " yPow = " << yPow );
-
+   static_assert( xPow >= -2 && xPow <= 2, "unsupported value of xPow" );
+   static_assert( yPow >= -2 && yPow <= 2, "unsupported value of yPow" );
    return this->spaceStepsProducts[ xPow + 2 ][ yPow + 2 ];
 }
 
diff --git a/src/TNL/Meshes/GridDetails/Grid3D_impl.h b/src/TNL/Meshes/GridDetails/Grid3D_impl.h
index ee3e10d4b106ffa119f2059ac2a62e892b0ab87e..40e8aa820e5e998432f4c5c16d06d4f0a4ab9608 100644
--- a/src/TNL/Meshes/GridDetails/Grid3D_impl.h
+++ b/src/TNL/Meshes/GridDetails/Grid3D_impl.h
@@ -10,8 +10,12 @@
 
 #pragma once
 
+#include <fstream>
 #include <iomanip>
+#include <TNL/String.h>
 #include <TNL/Assert.h>
+#include <TNL/Logger.h>
+#include <TNL/Meshes/GridDetails/GnuplotWriter.h>
 #include <TNL/Meshes/GridDetails/GridEntityGetter_impl.h>
 #include <TNL/Meshes/GridDetails/NeighbourGridEntityGetter3D_impl.h>
 #include <TNL/Meshes/GridDetails/Grid3D.h>
@@ -164,9 +168,9 @@ template< typename Real,
           typename Index >
 void Grid< 3, Real, Device, Index > :: setDimensions( const Index xSize, const Index ySize, const Index zSize )
 {
-   TNL_ASSERT( xSize > 0, std::cerr << "xSize = " << xSize );
-   TNL_ASSERT( ySize > 0, std::cerr << "ySize = " << ySize );
-   TNL_ASSERT( zSize > 0, std::cerr << "zSize = " << zSize );
+   TNL_ASSERT_GT( xSize, 0, "Grid size must be positive." );
+   TNL_ASSERT_GT( ySize, 0, "Grid size must be positive." );
+   TNL_ASSERT_GT( zSize, 0, "Grid size must be positive." );
 
    this->dimensions.x() = xSize;
    this->dimensions.y() = ySize;
@@ -331,13 +335,9 @@ const Real&
 Grid< 3, Real, Device, Index >::
 getSpaceStepsProducts() const
 {
-   TNL_ASSERT( xPow >= -2 && xPow <= 2,
-              std::cerr << " xPow = " << xPow );
-   TNL_ASSERT( yPow >= -2 && yPow <= 2,
-              std::cerr << " yPow = " << yPow );
-   TNL_ASSERT( zPow >= -2 && zPow <= 2,
-              std::cerr << " zPow = " << zPow );
-
+   static_assert( xPow >= -2 && xPow <= 2, "unsupported value of xPow" );
+   static_assert( yPow >= -2 && yPow <= 2, "unsupported value of yPow" );
+   static_assert( zPow >= -2 && zPow <= 2, "unsupported value of zPow" );
    return this->spaceStepsProducts[ xPow + 2 ][ yPow + 2 ][ zPow + 2 ];
 }
 
@@ -511,7 +511,9 @@ template< typename Real,
 bool Grid< 3, Real, Device, Index >::writeMesh( const String& fileName,
                                                    const String& format ) const
 {
-   TNL_ASSERT( false, std::cerr << "TODO: FIX THIS"); // TODO: FIX THIS
+   /*****
+    * TODO: implement this
+    */
    return true;
 }
 
diff --git a/src/TNL/Meshes/GridDetails/GridEntityGetter_impl.h b/src/TNL/Meshes/GridDetails/GridEntityGetter_impl.h
index 15ad36fac490dd33f351566bde7cbd0af1228044..49bc86956d659edff903b15927a1f4753c80fab2 100644
--- a/src/TNL/Meshes/GridDetails/GridEntityGetter_impl.h
+++ b/src/TNL/Meshes/GridDetails/GridEntityGetter_impl.h
@@ -44,10 +44,8 @@ class GridEntityGetter<
       static GridEntity getEntity( const GridType& grid,
                                    const IndexType& index )
       {
-         TNL_ASSERT( index >= 0 && index < grid.template getEntitiesCount< GridEntity >(),
-              std::cerr << " index = " << index
-                   << " grid.getEntitiesCount<>() = " << grid.template getEntitiesCount< GridEntity >()
-                   << " entityDimension = " << entityDimension );
+         TNL_ASSERT_GE( index, 0, "Index must be non-negative." );
+         TNL_ASSERT_LT( index, grid.template getEntitiesCount< GridEntity >(), "Index is out of bounds." );
          return GridEntity
             ( grid,
               CoordinatesType( index ),
@@ -59,11 +57,8 @@ class GridEntityGetter<
       static IndexType getEntityIndex( const GridType& grid,
                                        const GridEntity& entity )
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0 ) &&
-                    entity.getCoordinates() < grid.getDimensions() + CoordinatesType( 1 - entityDimension ),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " grid.getDimensions() = " << grid.getDimensions()
-                   << " EntityDimension = " << entityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), grid.getDimensions() + CoordinatesType( 1 - entityDimension ), "wrong coordinates" );
          return entity.getCoordinates().x();
       }
 };
@@ -90,10 +85,8 @@ class GridEntityGetter< Meshes::Grid< 2, Real, Device, Index >, GridEntity, 2 >
       static GridEntity getEntity( const GridType& grid,
                                    const IndexType& index )
       {
-         TNL_ASSERT( index >= 0 && index < grid.template getEntitiesCount< GridEntity >(),
-           std::cerr << " index = " << index
-                << " grid.getEntitiesCount<>() = " << grid.template getEntitiesCount< GridEntity >()
-                << " entityDimension = " << entityDimension );
+         TNL_ASSERT_GE( index, 0, "Index must be non-negative." );
+         TNL_ASSERT_LT( index, grid.template getEntitiesCount< GridEntity >(), "Index is out of bounds." );
 
          const CoordinatesType dimensions = grid.getDimensions();
 
@@ -109,10 +102,8 @@ class GridEntityGetter< Meshes::Grid< 2, Real, Device, Index >, GridEntity, 2 >
       static IndexType getEntityIndex( const GridType& grid,
                                        const GridEntity& entity )
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0 ) &&
-                    entity.getCoordinates() < grid.getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " grid.getDimensions() = " << grid.getDimensions() );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), grid.getDimensions(), "wrong coordinates" );
 
          //const CoordinatesType coordinates = entity.getCoordinates();
          //const CoordinatesType dimensions = grid.getDimensions();
@@ -142,10 +133,8 @@ class GridEntityGetter< Meshes::Grid< 2, Real, Device, Index >, GridEntity, 1 >
       static GridEntity getEntity( const GridType& grid,
                                    const IndexType& index )
       {
-         TNL_ASSERT( index >= 0 && index < grid.template getEntitiesCount< GridEntity >(),
-           std::cerr << " index = " << index
-                << " grid.getEntitiesCount<>() = " << grid.template getEntitiesCount< GridEntity >()
-                << " entityDimension = " << entityDimension );
+         TNL_ASSERT_GE( index, 0, "Index must be non-negative." );
+         TNL_ASSERT_LT( index, grid.template getEntitiesCount< GridEntity >(), "Index is out of bounds." );
  
          const CoordinatesType dimensions = grid.getDimensions();
 
@@ -171,11 +160,8 @@ class GridEntityGetter< Meshes::Grid< 2, Real, Device, Index >, GridEntity, 1 >
       static IndexType getEntityIndex( const GridType& grid,
                                        const GridEntity& entity )
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0 ) &&
-                    entity.getCoordinates() < grid.getDimensions() + abs( entity.getOrientation() ),
-                 std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                      << " grid.getDimensions() = " << grid.getDimensions()
-                      << " abs( entity.getOrientation() ) = " << abs( entity.getOrientation() ) );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), grid.getDimensions() + abs( entity.getOrientation() ), "wrong coordinates" );
  
          const CoordinatesType coordinates = entity.getCoordinates();
          const CoordinatesType dimensions = grid.getDimensions();
@@ -205,10 +191,8 @@ class GridEntityGetter< Meshes::Grid< 2, Real, Device, Index >, GridEntity, 0 >
       static GridEntity getEntity( const GridType& grid,
                                    const IndexType& index )
       {
-         TNL_ASSERT( index >= 0 && index < grid.template getEntitiesCount< GridEntity >(),
-           std::cerr << " index = " << index
-                << " grid.getEntitiesCount<>() = " << grid.template getEntitiesCount< GridEntity >()
-                << " entityDimension = " << entityDimension );
+         TNL_ASSERT_GE( index, 0, "Index must be non-negative." );
+         TNL_ASSERT_LT( index, grid.template getEntitiesCount< GridEntity >(), "Index is out of bounds." );
 
          const CoordinatesType dimensions = grid.getDimensions();
 
@@ -225,9 +209,8 @@ class GridEntityGetter< Meshes::Grid< 2, Real, Device, Index >, GridEntity, 0 >
       static IndexType getEntityIndex( const GridType& grid,
                                        const GridEntity& entity )
       {
-         TNL_ASSERT( entity.getCoordinates() >= 0 && entity.getCoordinates() <= grid.getDimensions(),
-            std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                 << " grid.getDimensions() = " << grid.getDimensions() );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LE( entity.getCoordinates(), grid.getDimensions(), "wrong coordinates" );
  
          const CoordinatesType coordinates = entity.getCoordinates();
          const CoordinatesType dimensions = grid.getDimensions();
@@ -258,10 +241,8 @@ class GridEntityGetter< Meshes::Grid< 3, Real, Device, Index >, GridEntity, 3 >
       static GridEntity getEntity( const GridType& grid,
                                    const IndexType& index )
       {
-         TNL_ASSERT( index >= 0 && index < grid.template getEntitiesCount< GridEntity >(),
-           std::cerr << " index = " << index
-                << " grid.getEntitiesCount<>() = " << grid.template getEntitiesCount< GridEntity >()
-                << " entityDimension = " << entityDimension );
+         TNL_ASSERT_GE( index, 0, "Index must be non-negative." );
+         TNL_ASSERT_LT( index, grid.template getEntitiesCount< GridEntity >(), "Index is out of bounds." );
 
          const CoordinatesType dimensions = grid.getDimensions();
 
@@ -278,10 +259,8 @@ class GridEntityGetter< Meshes::Grid< 3, Real, Device, Index >, GridEntity, 3 >
       static IndexType getEntityIndex( const GridType& grid,
                                        const GridEntity& entity )
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0, 0 ) &&
-                    entity.getCoordinates() < grid.getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " grid.getDimensions() = " << grid.getDimensions() );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), grid.getDimensions(), "wrong coordinates" );
 
          const CoordinatesType coordinates = entity.getCoordinates();
          const CoordinatesType dimensions = grid.getDimensions();
@@ -310,10 +289,8 @@ class GridEntityGetter< Meshes::Grid< 3, Real, Device, Index >, GridEntity, 2 >
       static GridEntity getEntity( const GridType& grid,
                                    const IndexType& index )
       {
-         TNL_ASSERT( index >= 0 && index < grid.template getEntitiesCount< GridEntity >(),
-           std::cerr << " index = " << index
-                << " grid.getEntitiesCount<>() = " << grid.template getEntitiesCount< GridEntity >()
-                << " entityDimension = " << entityDimension );
+         TNL_ASSERT_GE( index, 0, "Index must be non-negative." );
+         TNL_ASSERT_LT( index, grid.template getEntitiesCount< GridEntity >(), "Index is out of bounds." );
 
          const CoordinatesType dimensions = grid.getDimensions();
  
@@ -354,11 +331,8 @@ class GridEntityGetter< Meshes::Grid< 3, Real, Device, Index >, GridEntity, 2 >
       static IndexType getEntityIndex( const GridType& grid,
                                        const GridEntity& entity )
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0, 0 ) &&
-                    entity.getCoordinates() < grid.getDimensions() + abs( entity.getOrientation() ),
-                 std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                      << " grid.getDimensions() = " << grid.getDimensions()
-                      << " abs( entity.getOrientation() ) = " << abs( entity.getOrientation() ) );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), grid.getDimensions() + abs( entity.getOrientation() ), "wrong coordinates" );
  
          const CoordinatesType coordinates = entity.getCoordinates();
          const CoordinatesType dimensions = grid.getDimensions();
@@ -400,10 +374,8 @@ class GridEntityGetter< Meshes::Grid< 3, Real, Device, Index >, GridEntity, 1 >
       static GridEntity getEntity( const GridType& grid,
                                    const IndexType& index )
       {
-         TNL_ASSERT( index >= 0 && index < grid.template getEntitiesCount< GridEntity >(),
-           std::cerr << " index = " << index
-                << " grid.getEntitiesCount<>() = " << grid.template getEntitiesCount< GridEntity >()
-                << " entityDimension = " << entityDimension );
+         TNL_ASSERT_GE( index, 0, "Index must be non-negative." );
+         TNL_ASSERT_LT( index, grid.template getEntitiesCount< GridEntity >(), "Index is out of bounds." );
  
          const CoordinatesType dimensions = grid.getDimensions();
 
@@ -447,12 +419,10 @@ class GridEntityGetter< Meshes::Grid< 3, Real, Device, Index >, GridEntity, 1 >
       static IndexType getEntityIndex( const GridType& grid,
                                        const GridEntity& entity )
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0, 0 ) &&
-                    entity.getCoordinates() < grid.getDimensions() +
-                       CoordinatesType( 1, 1, 1 ) - entity.getBasis(),
-            std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                 << " grid.getDimensions() = " << grid.getDimensions()
-                 << " CoordinatesType( 1, 1, 1 ) - entity.getBasis() = " << CoordinatesType( 1, 1, 1 ) - entity.getBasis() );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(),
+                        grid.getDimensions() + CoordinatesType( 1, 1, 1 ) - entity.getBasis(),
+                        "wrong coordinates" );
  
          const CoordinatesType coordinates = entity.getCoordinates();
          const CoordinatesType dimensions = grid.getDimensions();
@@ -490,10 +460,8 @@ class GridEntityGetter< Meshes::Grid< 3, Real, Device, Index >, GridEntity, 0 >
       static GridEntity getEntity( const GridType& grid,
                                    const IndexType& index )
       {
-         TNL_ASSERT( index >= 0 && index < grid.template getEntitiesCount< GridEntity >(),
-           std::cerr << " index = " << index
-                << " grid.getEntitiesCount<>() = " << grid.template getEntitiesCount< GridEntity >()
-                << " entityDimension = " << entityDimension );
+         TNL_ASSERT_GE( index, 0, "Index must be non-negative." );
+         TNL_ASSERT_LT( index, grid.template getEntitiesCount< GridEntity >(), "Index is out of bounds." );
 
          const CoordinatesType dimensions = grid.getDimensions();
  
@@ -512,9 +480,8 @@ class GridEntityGetter< Meshes::Grid< 3, Real, Device, Index >, GridEntity, 0 >
       static IndexType getEntityIndex( const GridType& grid,
                                        const GridEntity& entity )
       {
-         TNL_ASSERT( entity.getCoordinates() >= 0 && entity.getCoordinates() <= grid.getDimensions(),
-            std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                 << " grid.getDimensions() = " << grid.getDimensions() );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LE( entity.getCoordinates(), grid.getDimensions(), "wrong coordinates" );
  
          const CoordinatesType coordinates = entity.getCoordinates();
          const CoordinatesType dimensions = grid.getDimensions();
@@ -527,4 +494,3 @@ class GridEntityGetter< Meshes::Grid< 3, Real, Device, Index >, GridEntity, 0 >
 
 } // namespace Meshes
 } // namespace TNL
-
diff --git a/src/TNL/Meshes/GridDetails/GridEntity_impl.h b/src/TNL/Meshes/GridDetails/GridEntity_impl.h
index e484e0185890ab81df433a402a916d834ce4c824..b7828f12f2681d307d7d32303d8d2e414a1a2195 100644
--- a/src/TNL/Meshes/GridDetails/GridEntity_impl.h
+++ b/src/TNL/Meshes/GridDetails/GridEntity_impl.h
@@ -139,13 +139,11 @@ getIndex() const
 {
    typedef Meshes::Grid< Dimension, Real, Device, Index > GridType;
    typedef typename GridType::template EntityType< EntityDimension > EntityType;
-   TNL_ASSERT( this->entityIndex >= 0 &&
-              this-> entityIndex < grid.template getEntitiesCount< EntityType >(),
-              std::cerr << "this->entityIndex = " << this->entityIndex
-                   << " grid.template getEntitiesCount< EntityDimension >() = " << grid.template getEntitiesCount< EntityType >() );
-   TNL_ASSERT( this->entityIndex == grid.getEntityIndex( *this ),
-              std::cerr << "this->entityIndex = " << this->entityIndex
-                   << " grid.getEntityIndex( *this ) = " << grid.getEntityIndex( *this ) );
+   TNL_ASSERT_GE( this->entityIndex, 0, "Entity index is not non-negative." );
+   TNL_ASSERT_LT( this->entityIndex, grid.template getEntitiesCount< EntityDimension >(),
+                  "Entity index is out of bounds." );
+   TNL_ASSERT_EQ( this->entityIndex, grid.getEntityIndex( *this ),
+                  "Wrong value of stored index." );
    return this->entityIndex;
 }
 
@@ -387,13 +385,11 @@ Index
 GridEntity< Meshes::Grid< Dimension, Real, Device, Index >, Dimension, Config >::
 getIndex() const
 {
-   TNL_ASSERT( this->entityIndex >= 0 &&
-              this-> entityIndex < grid.template getEntitiesCount< ThisType >(),
-              std::cerr << "this->entityIndex = " << this->entityIndex
-                   << " grid.template getEntitiesCount< Dimension >() = " << grid.template getEntitiesCount< ThisType >() );
-   TNL_ASSERT( this->entityIndex == grid.getEntityIndex( *this ),
-              std::cerr << "this->index = " << this->entityIndex
-                   << " grid.getEntityIndex( *this ) = " << grid.getEntityIndex( *this ) );
+   TNL_ASSERT_GE( this->entityIndex, 0, "Entity index is not non-negative." );
+   TNL_ASSERT_LT( this->entityIndex, grid.template getEntitiesCount< Dimension >(),
+                  "Entity index is out of bounds." );
+   TNL_ASSERT_EQ( this->entityIndex, grid.getEntityIndex( *this ),
+                  "Wrong value of stored index." );
    return this->entityIndex;
 }
 
@@ -605,13 +601,11 @@ getIndex() const
 {
    typedef Meshes::Grid< Dimension, Real, Device, Index > GridType;
    typedef typename GridType::Vertex Vertex;
-   TNL_ASSERT( this->entityIndex >= 0 &&
-              this-> entityIndex < grid.template getEntitiesCount< Vertex >(),
-              std::cerr << "this->entityIndex = " << this->entityIndex
-                   << " grid.template getEntitiesCount< 0 >() = " << grid.template getEntitiesCount< Vertex >() );
-   TNL_ASSERT( this->entityIndex == grid.getEntityIndex( *this ),
-              std::cerr << "this->entityIndex = " << this->entityIndex
-                   << " grid.getEntityIndex( *this ) = " << grid.getEntityIndex( *this ) );
+   TNL_ASSERT_GE( this->entityIndex, 0, "Entity index is not non-negative." );
+   TNL_ASSERT_LT( this->entityIndex, grid.template getEntitiesCount< 0 >(),
+                  "Entity index is out of bounds." );
+   TNL_ASSERT_EQ( this->entityIndex, grid.getEntityIndex( *this ),
+                  "Wrong value of stored index." );
    return this->entityIndex;
 }
 
diff --git a/src/TNL/Meshes/GridDetails/NeighbourGridEntityGetter1D_impl.h b/src/TNL/Meshes/GridDetails/NeighbourGridEntityGetter1D_impl.h
index a8e91bcefc73ef9e4b5eaaf41bf1bf4f8e6f0595..dc6430e1a8976aa5a9a2164d9f01eeeef2052d1d 100644
--- a/src/TNL/Meshes/GridDetails/NeighbourGridEntityGetter1D_impl.h
+++ b/src/TNL/Meshes/GridDetails/NeighbourGridEntityGetter1D_impl.h
@@ -56,11 +56,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       NeighbourGridEntityType getEntity() const
       {
-         TNL_ASSERT( this->entity.getCoordinates() >= CoordinatesType( 0 ) &&
-                    this->entity.getCoordinates() < this->entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << this->entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() + CoordinatesType( step ) >= CoordinatesType( 0 ) &&
                     entity.getCoordinates() + CoordinatesType( step ) < entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
@@ -73,11 +70,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       IndexType getEntityIndex() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0 ) &&
-                    entity.getCoordinates() < entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() + CoordinatesType( step ) >= CoordinatesType( 0 ) &&
                     entity.getCoordinates() + CoordinatesType( step ) < entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
@@ -136,11 +130,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       NeighbourGridEntityType getEntity() const
       {
-         TNL_ASSERT( this->entity.getCoordinates() >= CoordinatesType( 0 ) &&
-                    this->entity.getCoordinates() < this->entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << this->entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() + CoordinatesType( step ) >= CoordinatesType( 0 ) &&
                     entity.getCoordinates() + CoordinatesType( step ) < entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
@@ -153,11 +144,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       IndexType getEntityIndex() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0 ) &&
-                    entity.getCoordinates() < entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() + CoordinatesType( step ) >= CoordinatesType( 0 ) &&
                     entity.getCoordinates() + CoordinatesType( step ) < entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
@@ -237,11 +225,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       NeighbourGridEntityType getEntity() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0 ) &&
-                    entity.getCoordinates() < entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates().x() + step + ( step < 0 ) >= CoordinatesType( 0 ) &&
                     entity.getCoordinates().x() + step + ( step < 0 ) <= entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
@@ -254,11 +239,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       IndexType getEntityIndex() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0 ) &&
-                    entity.getCoordinates() < entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates().x() + step + ( step < 0 ) >= CoordinatesType( 0 ).x() &&
                     entity.getCoordinates().x() + step + ( step < 0 ) <= entity.getMesh().getDimensions().x(),
               std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
@@ -318,11 +300,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       NeighbourGridEntityType getEntity() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0 ) &&
-                    entity.getCoordinates() <= entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0 ), "wrong coordinates" );
+         TNL_ASSERT_LE( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates().x() + step - ( step > 0 ) >= CoordinatesType( 0 ) &&
                     entity.getCoordinates().x() + step - ( step > 0 ) < entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
@@ -335,11 +314,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       IndexType getEntityIndex() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0 ) &&
-                    entity.getCoordinates() <= entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0 ), "wrong coordinates" );
+         TNL_ASSERT_LE( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates().x() + step - ( step > 0 ) >= 0 &&
                     entity.getCoordinates().x() + step - ( step > 0 ) < entity.getMesh().getDimensions().x(),
               std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
@@ -394,11 +370,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       NeighbourGridEntityType getEntity() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0 ) &&
-                    entity.getCoordinates() <= entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0 ), "wrong coordinates" );
+         TNL_ASSERT_LE( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates().x() + step - ( step > 0 ) >= CoordinatesType( 0 ) &&
                     entity.getCoordinates().x() + step - ( step > 0 ) < entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
@@ -411,11 +384,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       IndexType getEntityIndex() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0 ) &&
-                    entity.getCoordinates() <= entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0 ), "wrong coordinates" );
+         TNL_ASSERT_LE( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates().x() + step - ( step > 0 ) >= CoordinatesType( 0 ) &&
                     entity.getCoordinates().x() + step - ( step > 0 ) < entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
@@ -470,11 +440,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       NeighbourGridEntityType getEntity() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0 ) &&
-                    entity.getCoordinates() <= entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0 ), "wrong coordinates" );
+         TNL_ASSERT_LE( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates().x() + step >= CoordinatesType( 0 ) &&
                     entity.getCoordinates().x() + step <= entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
@@ -487,11 +454,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       IndexType getEntityIndex() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0 ) &&
-                    entity.getCoordinates() <= entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0 ), "wrong coordinates" );
+         TNL_ASSERT_LE( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates().x() + step >= CoordinatesType( 0 ) &&
                     entity.getCoordinates().x() + step <= entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
@@ -512,4 +476,3 @@ class NeighbourGridEntityGetter<
 
 } // namespace Meshes
 } // namespace TNL
-
diff --git a/src/TNL/Meshes/GridDetails/NeighbourGridEntityGetter2D_impl.h b/src/TNL/Meshes/GridDetails/NeighbourGridEntityGetter2D_impl.h
index 2ccb4568fa73c30356d306f1a0bf98075481d3aa..219042e879f4a175c1830b0dc5accfc0794ecf1d 100644
--- a/src/TNL/Meshes/GridDetails/NeighbourGridEntityGetter2D_impl.h
+++ b/src/TNL/Meshes/GridDetails/NeighbourGridEntityGetter2D_impl.h
@@ -56,11 +56,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       NeighbourGridEntityType getEntity() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0 ) &&
-                    entity.getCoordinates() < entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() + CoordinatesType( stepX, stepY ) >= CoordinatesType( 0, 0 ) &&
                     entity.getCoordinates() + CoordinatesType( stepX, stepY ) < entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates()  + CoordinatesType( stepX, stepY ) = " << entity.getCoordinates()  + CoordinatesType( stepX, stepY )
@@ -75,11 +72,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       IndexType getEntityIndex() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0 ) &&
-                    entity.getCoordinates() < entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() + CoordinatesType( stepX, stepY ) >= CoordinatesType( 0, 0 ) &&
                     entity.getCoordinates() + CoordinatesType( stepX, stepY ) < entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates()  + CoordinatesType( stepX, stepY ) = " << entity.getCoordinates()  + CoordinatesType( stepX, stepY )
@@ -140,11 +134,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       NeighbourGridEntityType getEntity() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0 ) &&
-                    entity.getCoordinates() < entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() + CoordinatesType( stepX, stepY ) >= CoordinatesType( 0, 0 ) &&
                     entity.getCoordinates() + CoordinatesType( stepX, stepY ) < entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates()  + CoordinatesType( stepX, stepY ) = " << entity.getCoordinates()  + CoordinatesType( stepX, stepY )
@@ -159,11 +150,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       IndexType getEntityIndex() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0 ) &&
-                    entity.getCoordinates() < entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() + CoordinatesType( stepX, stepY ) >= CoordinatesType( 0, 0 ) &&
                     entity.getCoordinates() + CoordinatesType( stepX, stepY ) < entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates()  + CoordinatesType( stepX, stepY ) = " << entity.getCoordinates()  + CoordinatesType( stepX, stepY )
@@ -269,13 +257,9 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       NeighbourGridEntityType getEntity() const
       {
-         TNL_ASSERT( ! stepX + ! stepY == 1,
-                    std::cerr << "Only one of the steps can be non-zero: stepX = " << stepX << " stepY = " << stepY );
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0 ) &&
-                    entity.getCoordinates() < entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         static_assert( ! stepX + ! stepY == 1, "Only one of the steps can be non-zero." );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() +
                        CoordinatesType( stepX + ( stepX < 0 ),
                                         stepY + ( stepY < 0 ) ) >= CoordinatesType( 0, 0 ) &&
@@ -350,11 +334,8 @@ class NeighbourGridEntityGetter<
       {
          TNL_ASSERT( stepX != 0 && stepY != 0,
                     std::cerr << " stepX = " << stepX << " stepY = " << stepY );
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0 ) &&
-                    entity.getCoordinates() < entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() +
                        CoordinatesType( stepX + ( stepX < 0 ), stepY + ( stepY < 0 ) ) >= CoordinatesType( 0, 0 ) &&
                     entity.getCoordinates() +
@@ -429,11 +410,8 @@ class NeighbourGridEntityGetter<
                     ( ( !! stepY ) == ( !! entity.getOrientation().y() ) ),
                     std::cerr << "( stepX, stepY ) cannot be perpendicular to entity coordinates: stepX = " << stepX << " stepY = " << stepY
                          << " entity.getOrientation() = " << entity.getOrientation() );*/
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0 ) &&
-                    entity.getCoordinates() < entity.getMesh().getDimensions() + entity.getOrientation(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() + entity.getOrientation() = " << entity.getMesh().getDimensions() + entity.getOrientation()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() +
                        CoordinatesType( stepX - ( stepX > 0 ) * ( entity.getOrientation().x() != 0.0 ),
                                         stepY - ( stepY > 0 ) * ( entity.getOrientation().y() != 0.0 ) ) >= CoordinatesType( 0, 0 ) &&
@@ -502,11 +480,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       NeighbourGridEntityType getEntity() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0 ) &&
-                    entity.getCoordinates() <= entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LE( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() + CoordinatesType( stepX, stepY ) >= CoordinatesType( 0, 0 ) &&
                     entity.getCoordinates() + CoordinatesType( stepX, stepY ) <= entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates()  + CoordinatesType( stepX, stepY ) = " << entity.getCoordinates()  + CoordinatesType( stepX, stepY )
@@ -521,11 +496,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       IndexType getEntityIndex() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0 ) &&
-                    entity.getCoordinates() <= entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LE( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() + CoordinatesType( stepX, stepY ) >= CoordinatesType( 0, 0 ) &&
                     entity.getCoordinates() + CoordinatesType( stepX, stepY ) <= entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates()  + CoordinatesType( stepX, stepY ) = " << entity.getCoordinates()  + CoordinatesType( stepX, stepY )
diff --git a/src/TNL/Meshes/GridDetails/NeighbourGridEntityGetter3D_impl.h b/src/TNL/Meshes/GridDetails/NeighbourGridEntityGetter3D_impl.h
index 682850bb1b3a442f7bf6a227418859496c690375..b6128e00d456fe3ecde1684f15c9636503ac22b2 100644
--- a/src/TNL/Meshes/GridDetails/NeighbourGridEntityGetter3D_impl.h
+++ b/src/TNL/Meshes/GridDetails/NeighbourGridEntityGetter3D_impl.h
@@ -56,11 +56,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       NeighbourGridEntityType getEntity() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0, 0 ) &&
-                    entity.getCoordinates() < entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() + CoordinatesType( stepX, stepY ) >= CoordinatesType( 0, 0, 0 ) &&
                     entity.getCoordinates() + CoordinatesType( stepX, stepY ) < entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates()  + CoordinatesType( stepX, stepY ) = " << entity.getCoordinates()  + CoordinatesType( stepX, stepY )
@@ -75,11 +72,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       IndexType getEntityIndex() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0, 0 ) &&
-                    entity.getCoordinates() < entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() + CoordinatesType( stepX, stepY, stepZ ) >= CoordinatesType( 0, 0, 0 ) &&
                     entity.getCoordinates() + CoordinatesType( stepX, stepY, stepZ ) < entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates()  + CoordinatesType( stepX, stepY, stepZ ) = "
@@ -142,11 +136,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       NeighbourGridEntityType getEntity() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0, 0 ) &&
-                    entity.getCoordinates() < entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() + CoordinatesType( stepX, stepY, stepZ ) >= CoordinatesType( 0, 0, 0 ) &&
                     entity.getCoordinates() + CoordinatesType( stepX, stepY, stepZ ) < entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates()  + CoordinatesType( stepX, stepY ) = " << entity.getCoordinates()  + CoordinatesType( stepX, stepY, stepZ )
@@ -161,11 +152,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       IndexType getEntityIndex() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0, 0 ) &&
-                    entity.getCoordinates() < entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() + CoordinatesType( stepX, stepY, stepZ ) >= CoordinatesType( 0, 0, 0 ) &&
                     entity.getCoordinates() + CoordinatesType( stepX, stepY, stepZ ) < entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates()  + CoordinatesType( stepX, stepY, stepZ ) = "
@@ -290,15 +278,9 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       NeighbourGridEntityType getEntity() const
       {
-         TNL_ASSERT( ! stepX + ! stepY + ! stepZ == 2,
-                    std::cerr << "Only one of the steps can be non-zero: stepX = " << stepX
-                         << " stepY = " << stepY
-                         << " stepZ = " << stepZ );
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0, 0 ) &&
-                    entity.getCoordinates() < entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         static_assert( ! stepX + ! stepY + ! stepZ == 2, "Only one of the steps can be non-zero." );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() +
                        CoordinatesType( stepX + ( stepX < 0 ),
                                         stepY + ( stepY < 0 ),
@@ -379,15 +361,9 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       NeighbourGridEntityType getEntity() const
       {
-         TNL_ASSERT( ! stepX + ! stepY + ! stepZ == 2,
-                    std::cerr << "Only one of the steps can be non-zero: stepX = " << stepX
-                         << " stepY = " << stepY
-                         << " stepZ = " << stepZ );
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0, 0 ) &&
-                    entity.getCoordinates() < entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         static_assert( ! stepX + ! stepY + ! stepZ == 2, "Only one of the steps can be non-zero." );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() +
                        CoordinatesType( stepX + ( stepX < 0 ),
                                         stepY + ( stepY < 0 ),
@@ -469,15 +445,9 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       NeighbourGridEntityType getEntity() const
       {
-         TNL_ASSERT( ! stepX + ! stepY + ! stepZ == 1,
-                    std::cerr << "Exactly two of the steps must be non-zero: stepX = " << stepX
-                         << " stepY = " << stepY
-                         << " stepZ = " << stepZ );
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0, 0 ) &&
-                    entity.getCoordinates() < entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         static_assert( ! stepX + ! stepY + ! stepZ == 1, "Exactly two of the steps must be non-zero." );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() +
                        CoordinatesType( stepX + ( stepX < 0 ),
                                         stepY + ( stepY < 0 ),
@@ -558,11 +528,8 @@ class NeighbourGridEntityGetter<
                     std::cerr << " stepX = " << stepX
                          << " stepY = " << stepY
                          << " stepZ = " << stepZ );
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0, 0 ) &&
-                    entity.getCoordinates() < entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() +
                        CoordinatesType( stepX + ( stepX < 0 ),
                                         stepY + ( stepY < 0 ),
@@ -643,11 +610,8 @@ class NeighbourGridEntityGetter<
                     std::cerr << "( stepX, stepY, stepZ ) cannot be perpendicular to entity coordinates: stepX = " << stepX
                          << " stepY = " << stepY << " stepZ = " << stepZ
                          << " entity.getOrientation() = " << entity.getOrientation() );*/
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0, 0 ) &&
-                    entity.getCoordinates() < entity.getMesh().getDimensions() + entity.getOrientation(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() + entity.getOrientation() = " << entity.getMesh().getDimensions() + entity.getOrientation()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LT( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() +
                        CoordinatesType( stepX - ( stepX > 0 ) * ( entity.getOrientation().x() != 0.0 ),
                                         stepY - ( stepY > 0 ) * ( entity.getOrientation().y() != 0.0 ),
@@ -723,11 +687,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       NeighbourGridEntityType getEntity() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0, 0 ) &&
-                    entity.getCoordinates() <= entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LE( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() + CoordinatesType( stepX, stepY, stepZ ) >= CoordinatesType( 0, 0, 0 ) &&
                     entity.getCoordinates() + CoordinatesType( stepX, stepY, stepZ ) <= entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates()  + CoordinatesType( stepX, stepY, stepZ ) = "
@@ -743,11 +704,8 @@ class NeighbourGridEntityGetter<
       __cuda_callable__ inline
       IndexType getEntityIndex() const
       {
-         TNL_ASSERT( entity.getCoordinates() >= CoordinatesType( 0, 0, 0 ) &&
-                    entity.getCoordinates() <= entity.getMesh().getDimensions(),
-              std::cerr << "entity.getCoordinates() = " << entity.getCoordinates()
-                   << " entity.getMesh().getDimensions() = " << entity.getMesh().getDimensions()
-                   << " EntityDimension = " << EntityDimension );
+         TNL_ASSERT_GE( entity.getCoordinates(), CoordinatesType( 0, 0, 0 ), "wrong coordinates" );
+         TNL_ASSERT_LE( entity.getCoordinates(), entity.getMesh().getDimensions(), "wrong coordinates" );
          TNL_ASSERT( entity.getCoordinates() + CoordinatesType( stepX, stepY, stepZ ) >= CoordinatesType( 0, 0, 0 ) &&
                     entity.getCoordinates() + CoordinatesType( stepX, stepY, stepZ ) <= entity.getMesh().getDimensions(),
               std::cerr << "entity.getCoordinates()  + CoordinatesType( stepX, stepY, stepZ ) = "
@@ -770,4 +728,3 @@ class NeighbourGridEntityGetter<
 
 } // namespace Meshes
 } // namespace TNL
-
diff --git a/src/TNL/Operators/euler/fvm/LaxFridrichs_impl.h b/src/TNL/Operators/euler/fvm/LaxFridrichs_impl.h
index d0c650b0a1bf3d4f1d9a90952a7548f2b9700506..bb5784c1fbdbd5d47e4c31b24c56fd4e4ec3dabb 100644
--- a/src/TNL/Operators/euler/fvm/LaxFridrichs_impl.h
+++ b/src/TNL/Operators/euler/fvm/LaxFridrichs_impl.h
@@ -147,8 +147,8 @@ void LaxFridrichs< Meshes::Grid< 2, Real, Device, Index, GridGeometry >, Pressur
                                                                                                               RealType& rho_u2_t,
                                                                                                               const RealType& tau ) const
 {
-   TNL_ASSERT( mesh, std::cerr << "No mesh has been binded with the Lax-Fridrichs scheme." );
-   TNL_ASSERT( pressureGradient, std::cerr << "No pressure gradient was set in the the Lax-Fridrichs scheme." )
+   TNL_ASSERT_TRUE( mesh, "No mesh has been binded with the Lax-Fridrichs scheme." );
+   TNL_ASSERT_TRUE( pressureGradient, "No pressure gradient was set in the the Lax-Fridrichs scheme." )
 
    const IndexType& c = centralVolume;
    const IndexType e = this->mesh -> getElementNeighbour( centralVolume,  1,  0 );
@@ -398,8 +398,8 @@ void LaxFridrichs< Meshes::Grid< 2, Real, Device, Index, tnlIdenticalGridGeometr
                                                                                                                           RealType& rho_u2_t,
                                                                                                                           const RealType& tau ) const
 {
-   TNL_ASSERT( mesh, std::cerr << "No mesh has been binded with the Lax-Fridrichs scheme." );
-   TNL_ASSERT( pressureGradient, std::cerr << "No pressure gradient was set in the the Lax-Fridrichs scheme." )
+   TNL_ASSERT_TRUE( mesh, "No mesh has been binded with the Lax-Fridrichs scheme." );
+   TNL_ASSERT_TRUE( pressureGradient, "No pressure gradient was set in the the Lax-Fridrichs scheme." )
 
    const IndexType& xSize = this->mesh -> getDimensions(). x();
    const IndexType& ySize = this->mesh -> getDimensions(). y();
@@ -456,8 +456,8 @@ void LaxFridrichs< Meshes::Grid< 2, Real, Device, Index, tnlIdenticalGridGeometr
                                                                                                                           RealType& e_t,
                                                                                                                           const RealType& tau ) const
 {
-   TNL_ASSERT( mesh, std::cerr << "No mesh has been binded with the Lax-Fridrichs scheme." );
-   TNL_ASSERT( pressureGradient, std::cerr << "No pressure gradient was set in the the Lax-Fridrichs scheme." )
+   TNL_ASSERT_TRUE( mesh, "No mesh has been binded with the Lax-Fridrichs scheme." );
+   TNL_ASSERT_TRUE( pressureGradient, "No pressure gradient was set in the the Lax-Fridrichs scheme." )
 
    const IndexType& xSize = this->mesh -> getDimensions(). x();
    const IndexType& ySize = this->mesh -> getDimensions(). y();
diff --git a/src/TNL/Problems/cfd/navier-stokes/NavierStokesSolver_impl.h b/src/TNL/Problems/cfd/navier-stokes/NavierStokesSolver_impl.h
index 50f2138dd01edf6112ecd6d151d1776aa2910edf..229ff4667b4e6d0c84e1ba27dd5464f40dd94b7f 100644
--- a/src/TNL/Problems/cfd/navier-stokes/NavierStokesSolver_impl.h
+++ b/src/TNL/Problems/cfd/navier-stokes/NavierStokesSolver_impl.h
@@ -345,10 +345,10 @@ void NavierStokesSolver< AdvectionScheme,
                                                             SolverVectorType& u,
                                                             SolverVectorType& fu )
 {
-   TNL_ASSERT( this->advection, );
-   TNL_ASSERT( this->u1Viscosity, );
-   TNL_ASSERT( this->u2Viscosity, );
-   TNL_ASSERT( this->boundaryConditions, );
+   TNL_ASSERT_TRUE( this->advection, "advection scheme was not set" );
+   TNL_ASSERT_TRUE( this->u1Viscosity, "diffusion scheme was not set" );
+   TNL_ASSERT_TRUE( this->u2Viscosity, "diffusion scheme was not set" );
+   TNL_ASSERT_TRUE( this->boundaryConditions, "boundary conditions were not set" );
 
    SharedVector< RealType, DeviceType, IndexType > dofs_rho, dofs_rho_u1, dofs_rho_u2, dofs_e,
                                                       rho_t, rho_u1_t, rho_u2_t, e_t;
diff --git a/src/TNL/Solvers/Linear/CWYGMRES_impl.h b/src/TNL/Solvers/Linear/CWYGMRES_impl.h
index 4a9187775e9de0c63323f14c129707705775a864..01f966faa5b1df7298bdd603db893b79210ec310 100644
--- a/src/TNL/Solvers/Linear/CWYGMRES_impl.h
+++ b/src/TNL/Solvers/Linear/CWYGMRES_impl.h
@@ -125,7 +125,7 @@ bool
 CWYGMRES< Matrix, Preconditioner >::
 solve( const Vector& b, Vector& x )
 {
-   TNL_ASSERT( matrix, std::cerr << "No matrix was set in CWYGMRES. Call setMatrix() before solve()." << std::endl );
+   TNL_ASSERT_TRUE( matrix, "No matrix was set in CWYGMRES. Call setMatrix() before solve()." );
    if( restarting_min <= 0 || restarting_max <= 0 || restarting_min > restarting_max )
    {
       std::cerr << "Wrong value for the GMRES restarting parameters: r_min = " << restarting_min
diff --git a/src/TNL/Solvers/Linear/GMRES_impl.h b/src/TNL/Solvers/Linear/GMRES_impl.h
index 5b0b034761d92b3713c1697efe62968e1b0a6fe3..d3a20175926d9023228152f91a24ac04154fc80c 100644
--- a/src/TNL/Solvers/Linear/GMRES_impl.h
+++ b/src/TNL/Solvers/Linear/GMRES_impl.h
@@ -116,7 +116,7 @@ bool
 GMRES< Matrix, Preconditioner >::
 solve( const Vector& b, Vector& x )
 {
-   TNL_ASSERT( matrix, std::cerr << "No matrix was set in GMRES. Call setMatrix() before solve()." << std::endl );
+   TNL_ASSERT_TRUE( matrix, "No matrix was set in GMRES. Call setMatrix() before solve()." );
    if( restarting_min <= 0 || restarting_max <= 0 || restarting_min > restarting_max )
    {
       std::cerr << "Wrong value for the GMRES restarting parameters: r_min = " << restarting_min
diff --git a/src/TNL/Solvers/Linear/Preconditioners/Diagonal_impl.h b/src/TNL/Solvers/Linear/Preconditioners/Diagonal_impl.h
index 3a0f877dacf71e58244a9cd237e97978fb4ec71e..0c77239443c054c4429416c5c654ab53d8419a9e 100644
--- a/src/TNL/Solvers/Linear/Preconditioners/Diagonal_impl.h
+++ b/src/TNL/Solvers/Linear/Preconditioners/Diagonal_impl.h
@@ -57,7 +57,8 @@ update( const MatrixPointer& matrix )
 {
 //  std::cout << getType() << "->setMatrix()" << std::endl;
 
-   TNL_ASSERT( matrix->getRows() > 0 && matrix->getRows() == matrix->getColumns(), );
+   TNL_ASSERT_GT( matrix->getRows(), 0, "empty matrix" );
+   TNL_ASSERT_EQ( matrix->getRows(), matrix->getColumns(), "matrix must be square" );
 
    if( diagonal.getSize() != matrix->getRows() )
       diagonal.setSize( matrix->getRows() );
diff --git a/src/TNL/Solvers/Linear/Preconditioners/Dummy.h b/src/TNL/Solvers/Linear/Preconditioners/Dummy.h
index 5629568c78b0bf1571c6a2d25a46d4bdc0f3199e..2a7283b22bd31f881b2d8a320b09ddc50b3bbe8a 100644
--- a/src/TNL/Solvers/Linear/Preconditioners/Dummy.h
+++ b/src/TNL/Solvers/Linear/Preconditioners/Dummy.h
@@ -29,8 +29,7 @@ class Dummy
    template< typename Vector1, typename Vector2 >
    bool solve( const Vector1& b, Vector2& x ) const
    {
-      TNL_ASSERT( false,
-              std::cerr << "The solve() method of a dummy preconditioner should not be called." << std::endl; );
+      TNL_ASSERT_TRUE( false, "The solve() method of a dummy preconditioner should not be called." );
       return true;
    }
 
diff --git a/src/TNL/Solvers/Linear/Preconditioners/ILU0_impl.h b/src/TNL/Solvers/Linear/Preconditioners/ILU0_impl.h
index b7c773d40cc98b87e58323a7ab5c64d91cc81544..8ac5991e88b590fa1fd4662e0db29c1f2ae1b22d 100644
--- a/src/TNL/Solvers/Linear/Preconditioners/ILU0_impl.h
+++ b/src/TNL/Solvers/Linear/Preconditioners/ILU0_impl.h
@@ -14,7 +14,8 @@ void
 ILU0< Real, Devices::Host, Index >::
 update( const MatrixPointer& matrixPointer )
 {
-   TNL_ASSERT( matrixPointer->getRows() > 0 && matrixPointer->getRows() == matrixPointer->getColumns(), );
+   TNL_ASSERT_GT( matrixPointer->getRows(), 0, "empty matrix" );
+   TNL_ASSERT_EQ( matrixPointer->getRows(), matrixPointer->getColumns(), "matrix must be square" );
 
    const IndexType N = matrixPointer->getRows();
 
@@ -95,7 +96,8 @@ bool
 ILU0< Real, Devices::Host, Index >::
 solve( const Vector1& b, Vector2& x ) const
 {
-   TNL_ASSERT( b.getSize() == L.getRows() && x.getSize() == L.getRows(), );
+   TNL_ASSERT_EQ( b.getSize(), L.getRows(), "wrong size of the right hand side" );
+   TNL_ASSERT_EQ( x.getSize(), L.getRows(), "wrong size of the solution vector" );
 
    const IndexType N = x.getSize();
 
diff --git a/src/TNL/Solvers/Linear/UmfpackWrapper_impl.h b/src/TNL/Solvers/Linear/UmfpackWrapper_impl.h
index a07143a573900b4369cae396e43ea4b1a0797cc4..20d1b036899951e90d3dbefb94a4887348da01bd 100644
--- a/src/TNL/Solvers/Linear/UmfpackWrapper_impl.h
+++ b/src/TNL/Solvers/Linear/UmfpackWrapper_impl.h
@@ -63,8 +63,9 @@ bool UmfpackWrapper< Matrices::CSR< double, Devices::Host, int >, Preconditioner
 solve( const Vector& b,
        Vector& x )
 {
-    TNL_ASSERT( matrix->getRows() == matrix->getColumns(), );
-    TNL_ASSERT( matrix->getColumns() == x.getSize() && matrix->getColumns() == b.getSize(), );
+    TNL_ASSERT_EQ( matrix->getRows(), matrix->getColumns(), "matrix must be square" );
+    TNL_ASSERT_EQ( matrix->getColumns(), x.getSize(), "wrong size of the solution vector" );
+    TNL_ASSERT_EQ( matrix->getColumns(), b.getSize(), "wrong size of the right hand side" );
 
     const IndexType size = matrix -> getRows();
 
diff --git a/src/TNL/Solvers/PDE/ExplicitTimeStepper_impl.h b/src/TNL/Solvers/PDE/ExplicitTimeStepper_impl.h
index dff5f2e6c53b79ee348391d702fb4d53e32ef102..a2450e38d2acdd331da85a59ef3ca639c022660d 100644
--- a/src/TNL/Solvers/PDE/ExplicitTimeStepper_impl.h
+++ b/src/TNL/Solvers/PDE/ExplicitTimeStepper_impl.h
@@ -122,7 +122,7 @@ solve( const RealType& time,
        DofVectorPointer& dofVector,
        MeshDependentDataPointer& meshDependentData )
 {
-   TNL_ASSERT( this->odeSolver, );
+   TNL_ASSERT_TRUE( this->odeSolver, "ODE solver was not set" );
    mainTimer.start();
    this->odeSolver->setTau( this->timeStep );
    this->odeSolver->setProblem( * this );
diff --git a/src/TNL/Solvers/PDE/ExplicitUpdater.h b/src/TNL/Solvers/PDE/ExplicitUpdater.h
index cc4056b70eeb082154dedbf2e2c6318a8949a77a..ff66a0c6de9995eca234532bc0e7017fd49405fb 100644
--- a/src/TNL/Solvers/PDE/ExplicitUpdater.h
+++ b/src/TNL/Solvers/PDE/ExplicitUpdater.h
@@ -122,12 +122,12 @@ class ExplicitUpdater
                                                  typename MeshFunction::IndexType > >::value != true,
             "Error: I am getting Vector instead of MeshFunction or similar object. You might forget to bind DofVector into MeshFunction in you method getExplicitUpdate."  );
             
-         TNL_ASSERT( this->userDataPointer->differentialOperator, 
-            std::cerr << "The differential operator is not correctly set-up. Use method setDifferentialOperator() to do it." << std::endl );
-         TNL_ASSERT( this->userDataPointer->boundaryConditions, 
-            std::cerr << "The boundary conditions are not correctly set-up. Use method setBoundaryCondtions() to do it." << std::endl );
-         TNL_ASSERT( this->userDataPointer->rightHandSide, 
-            std::cerr << "The right-hand side is not correctly set-up. Use method setRightHandSide() to do it." << std::endl );
+         TNL_ASSERT_TRUE( this->userDataPointer->differentialOperator,
+                          "The differential operator is not correctly set-up. Use method setDifferentialOperator() to do it." );
+         TNL_ASSERT_TRUE( this->userDataPointer->boundaryConditions, 
+                          "The boundary conditions are not correctly set-up. Use method setBoundaryCondtions() to do it." );
+         TNL_ASSERT_TRUE( this->userDataPointer->rightHandSide, 
+                          "The right-hand side is not correctly set-up. Use method setRightHandSide() to do it." );
          
          
          this->userDataPointer->time = time;
diff --git a/src/TNL/Solvers/PDE/LinearSystemAssembler.h b/src/TNL/Solvers/PDE/LinearSystemAssembler.h
index 6b96f0b66cda28c970a38153e5ec637892746e88..d97dd15ef530ec33833d016d4ed321a83d711df6 100644
--- a/src/TNL/Solvers/PDE/LinearSystemAssembler.h
+++ b/src/TNL/Solvers/PDE/LinearSystemAssembler.h
@@ -114,7 +114,7 @@ class LinearSystemAssembler
       "Error: I am getting Vector instead of MeshFunction or similar object. You might forget to bind DofVector into MeshFunction in you method getExplicitUpdate."  );
 
       const IndexType maxRowLength = matrixPointer.template getData< Devices::Host >().getMaxRowLength();
-      TNL_ASSERT( maxRowLength > 0, );
+      TNL_ASSERT_GT( maxRowLength, 0, "maximum row length must be positive" );
       this->userDataPointer->time = time;
       this->userDataPointer->tau = tau;
       this->userDataPointer->u = &uPointer.template getData< DeviceType >();
diff --git a/src/TNL/Solvers/PDE/SemiImplicitTimeStepper_impl.h b/src/TNL/Solvers/PDE/SemiImplicitTimeStepper_impl.h
index 4b34d1df69d9453f931da492452e1500fa454c13..e9acc7556301878693b07ae46f72c58d21663ddf 100644
--- a/src/TNL/Solvers/PDE/SemiImplicitTimeStepper_impl.h
+++ b/src/TNL/Solvers/PDE/SemiImplicitTimeStepper_impl.h
@@ -153,7 +153,7 @@ solve( const RealType& time,
        DofVectorPointer& dofVector,
        MeshDependentDataPointer& meshDependentData )
 {
-   TNL_ASSERT( this->problem != 0, );
+   TNL_ASSERT_TRUE( this->problem, "problem was not set" );
    RealType t = time;
    this->linearSystemSolver->setMatrix( this->matrix );
    PreconditionerPointer preconditioner;
diff --git a/src/TNL/Solvers/PDE/TimeDependentPDESolver_impl.h b/src/TNL/Solvers/PDE/TimeDependentPDESolver_impl.h
index 6d63e9083d188eec064a83afd91ac2694d525e3c..6d20972d28c9f06cdf8da05e90567ac50f656e67 100644
--- a/src/TNL/Solvers/PDE/TimeDependentPDESolver_impl.h
+++ b/src/TNL/Solvers/PDE/TimeDependentPDESolver_impl.h
@@ -80,7 +80,7 @@ setup( const Config::ParameterContainer& parameters,
    /****
     * Set DOFs (degrees of freedom)
     */
-   TNL_ASSERT( problem->getDofs( this->meshPointer ) != 0, );
+   TNL_ASSERT_GT( problem->getDofs( this->meshPointer ), 0, "number of DOFs must be positive" );
    this->dofsPointer->setSize( problem->getDofs( this->meshPointer ) );
    this->dofsPointer->setValue( 0.0 );
    this->problem->bindDofs( this->meshPointer, this->dofsPointer );
@@ -316,10 +316,8 @@ bool
 TimeDependentPDESolver< Problem, TimeStepper >::
 solve()
 {
-   TNL_ASSERT( timeStepper != 0,
-              std::cerr << "No time stepper was set in PDESolver." );
-   TNL_ASSERT( problem != 0,
-              std::cerr << "No problem was set in PDESolver." );
+   TNL_ASSERT_TRUE( timeStepper, "No time stepper was set in PDESolver." );
+   TNL_ASSERT_TRUE( problem, "No problem was set in PDESolver." );
 
    if( snapshotPeriod == 0 )
    {
diff --git a/src/TNL/Solvers/PDE/TimeIndependentPDESolver_impl.h b/src/TNL/Solvers/PDE/TimeIndependentPDESolver_impl.h
index 4cf63b4ac44d59c99a445d81cfffa75cdc450eef..957f710300bbc373b184dbd5736028377b68a7ba 100644
--- a/src/TNL/Solvers/PDE/TimeIndependentPDESolver_impl.h
+++ b/src/TNL/Solvers/PDE/TimeIndependentPDESolver_impl.h
@@ -58,7 +58,7 @@ setup( const Config::ParameterContainer& parameters,
    /****
     * Set DOFs (degrees of freedom)
     */
-   TNL_ASSERT( problem->getDofs( this->mesh ) != 0, );
+   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->problem->bindDofs( this->mesh, this->dofs );   
@@ -141,8 +141,7 @@ bool
 tnlTimeIndependentPDESolver< Problem >::
 solve()
 {
-   TNL_ASSERT( problem != 0,
-              cerr << "No problem was set in tnlPDESolver." );
+   TNL_ASSERT_TRUE( problem, "No problem was set in tnlPDESolver." );
 
    this->computeTimer->reset();
    this->computeTimer->start();