diff --git a/src/TNL/Algorithms/detail/DistributedScan.h b/src/TNL/Algorithms/detail/DistributedScan.h
index 4ed71db7ea0a148cf1ac69ff627f095aa6b903a2..94e1d11a750cf2de56c8b0dc38017bccf6733728 100644
--- a/src/TNL/Algorithms/detail/DistributedScan.h
+++ b/src/TNL/Algorithms/detail/DistributedScan.h
@@ -34,7 +34,7 @@ struct DistributedScan
       using ValueType = typename OutputDistributedArray::ValueType;
       using DeviceType = typename OutputDistributedArray::DeviceType;
 
-      const auto communicator = input.getCommunicator();
+      const auto& communicator = input.getCommunicator();
       if( communicator != MPI_COMM_NULL ) {
          // adjust begin and end for the local range
          const auto localRange = input.getLocalRange();
@@ -49,7 +49,7 @@ struct DistributedScan
          const ValueType local_result = block_results.getElement( block_results.getSize() - 1 );
 
          // exchange local results between ranks
-         const int nproc = MPI::GetSize( communicator );
+         const int nproc = communicator.size();
          std::unique_ptr< ValueType[] > dataForScatter{ new ValueType[ nproc ] };
          for( int i = 0; i < nproc; i++ )
             dataForScatter[ i ] = local_result;
@@ -62,7 +62,7 @@ struct DistributedScan
             rank_results, rank_results, 0, nproc, 0, reduction, identity );
 
          // perform the second phase, using the per-block and per-rank results
-         const int rank = MPI::GetRank( communicator );
+         const int rank = communicator.rank();
          Scan< DeviceType, Type, PhaseType >::performSecondPhase(
             inputLocalView, outputLocalView, block_results, begin, end, begin, reduction, identity, rank_results[ rank ] );
       }
diff --git a/src/TNL/Containers/DistributedArray.h b/src/TNL/Containers/DistributedArray.h
index b05fdf0322d067deebcbb455bb6703754d602650..1f40b9c4edb64d7b916082ddd6cc8d92ac33adb9 100644
--- a/src/TNL/Containers/DistributedArray.h
+++ b/src/TNL/Containers/DistributedArray.h
@@ -78,11 +78,11 @@ public:
    DistributedArray( LocalRangeType localRange,
                      Index ghosts,
                      Index globalSize,
-                     MPI_Comm communicator,
+                     const MPI::Comm& communicator,
                      const AllocatorType& allocator = AllocatorType() );
 
    void
-   setDistribution( LocalRangeType localRange, Index ghosts, Index globalSize, MPI_Comm communicator );
+   setDistribution( LocalRangeType localRange, Index ghosts, Index globalSize, const MPI::Comm& communicator );
 
    const LocalRangeType&
    getLocalRange() const;
@@ -90,7 +90,7 @@ public:
    IndexType
    getGhosts() const;
 
-   MPI_Comm
+   const MPI::Comm&
    getCommunicator() const;
 
    AllocatorType
diff --git a/src/TNL/Containers/DistributedArray.hpp b/src/TNL/Containers/DistributedArray.hpp
index 367633d8fc9040c994f8d9584aa4b558709c845b..4dc65636e17fcdb1d2878bd2c5efa6e02b909be0 100644
--- a/src/TNL/Containers/DistributedArray.hpp
+++ b/src/TNL/Containers/DistributedArray.hpp
@@ -45,7 +45,7 @@ template< typename Value, typename Device, typename Index, typename Allocator >
 DistributedArray< Value, Device, Index, Allocator >::DistributedArray( LocalRangeType localRange,
                                                                        IndexType ghosts,
                                                                        IndexType globalSize,
-                                                                       MPI_Comm communicator,
+                                                                       const MPI::Comm& communicator,
                                                                        const Allocator& allocator )
 : localData( allocator )
 {
@@ -57,7 +57,7 @@ void
 DistributedArray< Value, Device, Index, Allocator >::setDistribution( LocalRangeType localRange,
                                                                       IndexType ghosts,
                                                                       IndexType globalSize,
-                                                                      MPI_Comm communicator )
+                                                                      const MPI::Comm& communicator )
 {
    TNL_ASSERT_LE( localRange.getEnd(), globalSize, "end of the local range is outside of the global range" );
    if( communicator != MPI_COMM_NULL )
@@ -80,7 +80,7 @@ DistributedArray< Value, Device, Index, Allocator >::getGhosts() const
 }
 
 template< typename Value, typename Device, typename Index, typename Allocator >
-MPI_Comm
+const MPI::Comm&
 DistributedArray< Value, Device, Index, Allocator >::getCommunicator() const
 {
    return view.getCommunicator();
diff --git a/src/TNL/Containers/DistributedArrayView.h b/src/TNL/Containers/DistributedArrayView.h
index 0c4a616b2852895c80bba3fe6b8a0eefc14c4655..ef09ce3f800b9b429150694e2a2fb29f3dcbf39f 100644
--- a/src/TNL/Containers/DistributedArrayView.h
+++ b/src/TNL/Containers/DistributedArrayView.h
@@ -44,9 +44,9 @@ public:
    DistributedArrayView( const LocalRangeType& localRange,
                          IndexType ghosts,
                          IndexType globalSize,
-                         MPI_Comm communicator,
+                         MPI::Comm communicator,
                          LocalViewType localData )
-   : localRange( localRange ), ghosts( ghosts ), globalSize( globalSize ), communicator( communicator ), localData( localData )
+   : localRange( localRange ), ghosts( ghosts ), globalSize( globalSize ), communicator( std::move( communicator ) ), localData( localData )
    {
       TNL_ASSERT_EQ( localData.getSize(),
                      localRange.getSize() + ghosts,
@@ -71,7 +71,7 @@ public:
    bind( const LocalRangeType& localRange,
          IndexType ghosts,
          IndexType globalSize,
-         MPI_Comm communicator,
+         const MPI::Comm& communicator,
          LocalViewType localData );
 
    // Note that you can also bind directly to DistributedArray and other types implicitly
@@ -91,7 +91,7 @@ public:
    IndexType
    getGhosts() const;
 
-   MPI_Comm
+   const MPI::Comm&
    getCommunicator() const;
 
    LocalViewType
@@ -256,7 +256,7 @@ protected:
    LocalRangeType localRange;
    IndexType ghosts = 0;
    IndexType globalSize = 0;
-   MPI_Comm communicator = MPI_COMM_NULL;
+   MPI::Comm communicator = MPI_COMM_NULL;
    LocalViewType localData;
 
    std::shared_ptr< SynchronizerType > synchronizer = nullptr;
diff --git a/src/TNL/Containers/DistributedArrayView.hpp b/src/TNL/Containers/DistributedArrayView.hpp
index f260c005f8a7e0880df1ffda00ca936359a68759..7b502a5c3c729f8116a66f369c97e779bacd97f5 100644
--- a/src/TNL/Containers/DistributedArrayView.hpp
+++ b/src/TNL/Containers/DistributedArrayView.hpp
@@ -38,7 +38,7 @@ void
 DistributedArrayView< Value, Device, Index >::bind( const LocalRangeType& localRange,
                                                     IndexType ghosts,
                                                     IndexType globalSize,
-                                                    MPI_Comm communicator,
+                                                    const MPI::Comm& communicator,
                                                     LocalViewType localData )
 {
    TNL_ASSERT_EQ( localData.getSize(),
@@ -93,7 +93,7 @@ DistributedArrayView< Value, Device, Index >::getGhosts() const
 }
 
 template< typename Value, typename Device, typename Index >
-MPI_Comm
+const MPI::Comm&
 DistributedArrayView< Value, Device, Index >::getCommunicator() const
 {
    return communicator;
diff --git a/src/TNL/Containers/DistributedNDArray.h b/src/TNL/Containers/DistributedNDArray.h
index 4f9e99c26b70f04ebc937d2026bf8de204df0827..1725be01e827e9c745adc201d858531a5c9669f2 100644
--- a/src/TNL/Containers/DistributedNDArray.h
+++ b/src/TNL/Containers/DistributedNDArray.h
@@ -83,7 +83,7 @@ public:
       return localArray.getAllocator();
    }
 
-   MPI_Comm
+   const MPI::Comm&
    getCommunicator() const
    {
       return communicator;
@@ -389,7 +389,7 @@ public:
 
    template< std::size_t level >
    void
-   setDistribution( IndexType begin, IndexType end, MPI_Comm communicator = MPI_COMM_WORLD )
+   setDistribution( IndexType begin, IndexType end, const MPI::Comm& communicator = MPI_COMM_WORLD )
    {
       static_assert( SizesHolderType::template getStaticSize< level >() == 0,
                      "NDArray cannot be distributed in static dimensions." );
@@ -472,7 +472,7 @@ public:
 
 protected:
    NDArray localArray;
-   MPI_Comm communicator = MPI_COMM_NULL;
+   MPI::Comm communicator = MPI_COMM_NULL;
    SizesHolderType globalSizes;
    // static sizes should have different type: localBegin is always 0, localEnd is always the full size
    LocalBeginsType localBegins;
diff --git a/src/TNL/Containers/DistributedNDArraySynchronizer.h b/src/TNL/Containers/DistributedNDArraySynchronizer.h
index 6391c302d5285e7f0187db341ac9b9829361e087..7e97b5dbfcc147cb67366a3f6a41a58b76dd3af4 100644
--- a/src/TNL/Containers/DistributedNDArraySynchronizer.h
+++ b/src/TNL/Containers/DistributedNDArraySynchronizer.h
@@ -345,7 +345,7 @@ protected:
 
       // issue all send and receive async operations
       RequestsVector requests;
-      const MPI_Comm communicator = array_view.getCommunicator();
+      const MPI::Comm& communicator = array_view.getCommunicator();
       Algorithms::staticFor< std::size_t, 0, DistributedNDArray::getDimension() >(
          [ & ]( auto dim )
          {
@@ -419,9 +419,9 @@ protected:
       dim_buffers.right_recv_offsets.template setSize< dim >( localEnds.template getSize< dim >() );
 
       // set default neighbor IDs
-      const MPI_Comm communicator = array_view.getCommunicator();
-      const int rank = MPI::GetRank( communicator );
-      const int nproc = MPI::GetSize( communicator );
+      const MPI::Comm& communicator = array_view.getCommunicator();
+      const int rank = communicator.rank();
+      const int nproc = communicator.size();
       if( dim_buffers.left_neighbor < 0 )
          dim_buffers.left_neighbor = ( rank + nproc - 1 ) % nproc;
       if( dim_buffers.right_neighbor < 0 )
@@ -485,7 +485,7 @@ protected:
    static void
    sendHelper( Buffers& buffers,
                RequestsVector& requests,
-               MPI_Comm communicator,
+               const MPI::Comm& communicator,
                int tag_from_left,
                int tag_to_left,
                int tag_from_right,
diff --git a/src/TNL/Containers/DistributedNDArrayView.h b/src/TNL/Containers/DistributedNDArrayView.h
index f77d4d8d9e9d4037eb12e3f50f057c7c8a9e9e5e..3bba5d7da746b89c35d5a00e9a1da81ea05770db 100644
--- a/src/TNL/Containers/DistributedNDArrayView.h
+++ b/src/TNL/Containers/DistributedNDArrayView.h
@@ -40,8 +40,8 @@ public:
                            SizesHolderType globalSizes,
                            LocalBeginsType localBegins,
                            SizesHolderType localEnds,
-                           MPI_Comm communicator )
-   : localView( localView ), communicator( communicator ), globalSizes( globalSizes ), localBegins( localBegins ),
+                           MPI::Comm communicator )
+   : localView( localView ), communicator( std::move( communicator ) ), globalSizes( globalSizes ), localBegins( localBegins ),
      localEnds( localEnds )
    {}
 
@@ -113,7 +113,7 @@ public:
       return NDArrayView::getDimension();
    }
 
-   MPI_Comm
+   const MPI::Comm&
    getCommunicator() const
    {
       return communicator;
@@ -405,7 +405,7 @@ public:
 
 protected:
    NDArrayView localView;
-   MPI_Comm communicator = MPI_COMM_NULL;
+   MPI::Comm communicator = MPI_COMM_NULL;
    SizesHolderType globalSizes;
    // static sizes should have different type: localBegin is always 0, localEnd is always the full size
    LocalBeginsType localBegins;
diff --git a/src/TNL/Containers/Expressions/DistributedExpressionTemplates.h b/src/TNL/Containers/Expressions/DistributedExpressionTemplates.h
index 1946c02f7ba7140d09b7ed7b40505755a61b373c..30003e1cdbcc2dfd37c2c8fd1860d67f4527a947 100644
--- a/src/TNL/Containers/Expressions/DistributedExpressionTemplates.h
+++ b/src/TNL/Containers/Expressions/DistributedExpressionTemplates.h
@@ -109,7 +109,7 @@ struct DistributedBinaryExpressionTemplate< T1, T2, Operation, VectorExpressionV
       return op1.getGhosts();
    }
 
-   MPI_Comm
+   const MPI::Comm&
    getCommunicator() const
    {
       return op1.getCommunicator();
@@ -201,7 +201,7 @@ struct DistributedBinaryExpressionTemplate< T1, T2, Operation, VectorExpressionV
       return op1.getGhosts();
    }
 
-   MPI_Comm
+   const MPI::Comm&
    getCommunicator() const
    {
       return op1.getCommunicator();
@@ -292,7 +292,7 @@ struct DistributedBinaryExpressionTemplate< T1, T2, Operation, ArithmeticVariabl
       return op2.getGhosts();
    }
 
-   MPI_Comm
+   const MPI::Comm&
    getCommunicator() const
    {
       return op2.getCommunicator();
@@ -385,7 +385,7 @@ struct DistributedUnaryExpressionTemplate
       return operand.getGhosts();
    }
 
-   MPI_Comm
+   const MPI::Comm&
    getCommunicator() const
    {
       return operand.getCommunicator();
diff --git a/src/TNL/Containers/Expressions/DistributedVerticalOperations.h b/src/TNL/Containers/Expressions/DistributedVerticalOperations.h
index a1d6a24ad0128105be658bdff4cf6fb94ecbaff0..b2409fd2161b5c3c64b2d1bad0bf719265fec308 100644
--- a/src/TNL/Containers/Expressions/DistributedVerticalOperations.h
+++ b/src/TNL/Containers/Expressions/DistributedVerticalOperations.h
@@ -6,9 +6,11 @@
 
 #pragma once
 
+#include <memory>
+
+#include <TNL/MPI/Comm.h>
 #include <TNL/MPI/Wrappers.h>
 #include <TNL/Algorithms/reduce.h>
-#include <memory>
 
 namespace TNL {
 namespace Containers {
@@ -42,7 +44,7 @@ DistributedExpressionArgMin( const Expression& expression )
    static_assert( std::numeric_limits< RealType >::is_specialized,
                   "std::numeric_limits is not specialized for the reduction's real type" );
    ResultType result( -1, std::numeric_limits< RealType >::max() );
-   const MPI_Comm communicator = expression.getCommunicator();
+   const MPI::Comm& communicator = expression.getCommunicator();
    if( communicator != MPI_COMM_NULL ) {
       // compute local argMin
       ResultType localResult = Algorithms::reduceWithArgument( expression.getConstLocalView(), TNL::MinWithArg{} );
@@ -101,7 +103,7 @@ DistributedExpressionArgMax( const Expression& expression )
    static_assert( std::numeric_limits< RealType >::is_specialized,
                   "std::numeric_limits is not specialized for the reduction's real type" );
    ResultType result( -1, std::numeric_limits< RealType >::lowest() );
-   const MPI_Comm communicator = expression.getCommunicator();
+   const MPI::Comm& communicator = expression.getCommunicator();
    if( communicator != MPI_COMM_NULL ) {
       // compute local argMax
       ResultType localResult = Algorithms::reduceWithArgument( expression.getConstLocalView(), TNL::MaxWithArg{} );
diff --git a/src/TNL/Containers/Partitioner.h b/src/TNL/Containers/Partitioner.h
index e7240a72cf8a3180b158f31867a66aae8402197d..5533bae4043ef442653e8eafc884804797d215b7 100644
--- a/src/TNL/Containers/Partitioner.h
+++ b/src/TNL/Containers/Partitioner.h
@@ -8,6 +8,7 @@
 
 #pragma once
 
+#include <utility>
 #include <vector>
 
 #include "Subrange.h"
@@ -25,11 +26,11 @@ public:
    using SubrangeType = Subrange< Index >;
 
    static SubrangeType
-   splitRange( Index globalSize, MPI_Comm communicator )
+   splitRange( Index globalSize, const MPI::Comm& communicator )
    {
       if( communicator != MPI_COMM_NULL ) {
-         const int rank = MPI::GetRank( communicator );
-         const int partitions = MPI::GetSize( communicator );
+         const int rank = communicator.rank();
+         const int partitions = communicator.size();
          const Index begin = TNL::min( globalSize, rank * globalSize / partitions );
          const Index end = TNL::min( globalSize, ( rank + 1 ) * globalSize / partitions );
          return SubrangeType( begin, end );
@@ -76,7 +77,7 @@ public:
 
       SubrangeType localRange;
       int overlaps;
-      MPI_Comm communicator;
+      MPI::Comm communicator;
 
    public:
       using ByteArrayView = typename Base::ByteArrayView;
@@ -91,24 +92,24 @@ public:
 
       ArraySynchronizer() = delete;
 
-      ArraySynchronizer( SubrangeType localRange, int overlaps, MPI_Comm communicator )
-      : localRange( localRange ), overlaps( overlaps ), communicator( communicator )
+      ArraySynchronizer( SubrangeType localRange, int overlaps, MPI::Comm communicator )
+      : localRange( localRange ), overlaps( overlaps ), communicator( std::move( communicator ) )
       {}
 
-      virtual void
+      void
       synchronizeByteArray( ByteArrayView array, int bytesPerValue ) override
       {
          auto requests = synchronizeByteArrayAsyncWorker( array, bytesPerValue );
          MPI::Waitall( requests.data(), requests.size() );
       }
 
-      virtual RequestsVector
+      RequestsVector
       synchronizeByteArrayAsyncWorker( ByteArrayView array, int bytesPerValue ) override
       {
          TNL_ASSERT_EQ( array.getSize(), bytesPerValue * ( localRange.getSize() + 2 * overlaps ), "unexpected array size" );
 
-         const int rank = MPI::GetRank( communicator );
-         const int nproc = MPI::GetSize( communicator );
+         const int rank = communicator.rank();
+         const int nproc = communicator.size();
          const int left = ( rank > 0 ) ? rank - 1 : nproc - 1;
          const int right = ( rank < nproc - 1 ) ? rank + 1 : 0;
 
diff --git a/src/TNL/Functions/MeshFunctionIO.h b/src/TNL/Functions/MeshFunctionIO.h
index a93fc2d32660ca20e41e12dc5086e30dd1743f76..58bb603bd4f652511206126c8f8869a7e308f37c 100644
--- a/src/TNL/Functions/MeshFunctionIO.h
+++ b/src/TNL/Functions/MeshFunctionIO.h
@@ -235,9 +235,9 @@ writeDistributedMeshFunction(
    }
 
    if( format == "pvti" ) {
-      const MPI_Comm communicator = distributedMesh.getCommunicator();
+      const MPI::Comm& communicator = distributedMesh.getCommunicator();
       std::ofstream file;
-      if( TNL::MPI::GetRank( communicator ) == 0 )
+      if( communicator.rank() == 0 )
          file.open( fileName );
 
       using PVTI = Meshes::Writers::PVTIWriter< typename MeshFunction::MeshType >;
diff --git a/src/TNL/Matrices/DistributedMatrix.h b/src/TNL/Matrices/DistributedMatrix.h
index ccc901d9db519ac47d82699772df1a3ea08993bc..9cadc3a4dca3835bdb38770d9a6d83fcc1478b3c 100644
--- a/src/TNL/Matrices/DistributedMatrix.h
+++ b/src/TNL/Matrices/DistributedMatrix.h
@@ -42,15 +42,15 @@ public:
 
    DistributedMatrix( DistributedMatrix& ) = default;
 
-   DistributedMatrix( LocalRangeType localRowRange, IndexType rows, IndexType columns, MPI_Comm communicator );
+   DistributedMatrix( LocalRangeType localRowRange, IndexType rows, IndexType columns, const MPI::Comm& communicator );
 
    void
-   setDistribution( LocalRangeType localRowRange, IndexType rows, IndexType columns, MPI_Comm communicator );
+   setDistribution( LocalRangeType localRowRange, IndexType rows, IndexType columns, const MPI::Comm& communicator );
 
    const LocalRangeType&
    getLocalRowRange() const;
 
-   MPI_Comm
+   const MPI::Comm&
    getCommunicator() const;
 
    const Matrix&
@@ -212,7 +212,7 @@ public:
 protected:
    LocalRangeType localRowRange;
    IndexType rows = 0;  // global rows count
-   MPI_Comm communicator = MPI_COMM_NULL;
+   MPI::Comm communicator = MPI_COMM_NULL;
    Matrix localMatrix;
 };
 
diff --git a/src/TNL/Matrices/DistributedMatrix_impl.h b/src/TNL/Matrices/DistributedMatrix_impl.h
index f48124946b4d88b4fb70db4a3321f1fd2a1ba32f..5cea862837ea9bef10dcf9716ce674bc55805b07 100644
--- a/src/TNL/Matrices/DistributedMatrix_impl.h
+++ b/src/TNL/Matrices/DistributedMatrix_impl.h
@@ -17,7 +17,7 @@ template< typename Matrix >
 DistributedMatrix< Matrix >::DistributedMatrix( LocalRangeType localRowRange,
                                                 IndexType rows,
                                                 IndexType columns,
-                                                MPI_Comm communicator )
+                                                const MPI::Comm& communicator )
 {
    setDistribution( localRowRange, rows, columns, communicator );
 }
@@ -27,12 +27,12 @@ void
 DistributedMatrix< Matrix >::setDistribution( LocalRangeType localRowRange,
                                               IndexType rows,
                                               IndexType columns,
-                                              MPI_Comm communicator )
+                                              const MPI::Comm& communicator )
 {
    this->localRowRange = localRowRange;
    this->rows = rows;
    this->communicator = communicator;
-   if( communicator != MPI_COMM_NULL )
+   if( this->communicator != MPI_COMM_NULL )
       localMatrix.setDimensions( localRowRange.getSize(), columns );
 }
 
@@ -44,7 +44,7 @@ DistributedMatrix< Matrix >::getLocalRowRange() const
 }
 
 template< typename Matrix >
-MPI_Comm
+const MPI::Comm&
 DistributedMatrix< Matrix >::getCommunicator() const
 {
    return communicator;
diff --git a/src/TNL/Meshes/DistributedMeshes/DistributedGrid.h b/src/TNL/Meshes/DistributedMeshes/DistributedGrid.h
index 2c4becf7e53c5e7fc62c14177226801c32e0dfab..5c3f96f578f35e9e6ab722b66f317f96f105dbf1 100644
--- a/src/TNL/Meshes/DistributedMeshes/DistributedGrid.h
+++ b/src/TNL/Meshes/DistributedMeshes/DistributedGrid.h
@@ -100,7 +100,7 @@ public:
    getSubdomainCoordinates() const;
 
    void
-   setCommunicator( MPI::Comm&& communicator );
+   setCommunicator( const MPI::Comm& communicator );
 
    const MPI::Comm&
    getCommunicator() const;
diff --git a/src/TNL/Meshes/DistributedMeshes/DistributedGrid.hpp b/src/TNL/Meshes/DistributedMeshes/DistributedGrid.hpp
index 33cf5d3b1e82464c0985a05816cca7e7d5c7adbf..01d4b2dd98b4123f58721a9ab5ff96e4fc68a036 100644
--- a/src/TNL/Meshes/DistributedMeshes/DistributedGrid.hpp
+++ b/src/TNL/Meshes/DistributedMeshes/DistributedGrid.hpp
@@ -237,9 +237,9 @@ DistributedMesh< Grid< Dimension, Real, Device, Index > >::getEntitiesCount() co
 
 template< int Dimension, typename Real, typename Device, typename Index >
 void
-DistributedMesh< Grid< Dimension, Real, Device, Index > >::setCommunicator( MPI::Comm&& communicator )
+DistributedMesh< Grid< Dimension, Real, Device, Index > >::setCommunicator( const MPI::Comm& communicator )
 {
-   this->communicator = std::move( communicator );
+   this->communicator = communicator;
 }
 
 template< int Dimension, typename Real, typename Device, typename Index >
@@ -337,7 +337,7 @@ DistributedMesh< Grid< Dimension, Real, Device, Index > >::SetupByCut(
    }
 
    // create new communicator with used nodes
-   const MPI_Comm oldCommunicator = inputDistributedGrid.getCommunicator();
+   const MPI::Comm oldCommunicator = inputDistributedGrid.getCommunicator();
    if( isInCut ) {
       this->isSet = true;
 
@@ -383,7 +383,7 @@ DistributedMesh< Grid< Dimension, Real, Device, Index > >::SetupByCut(
       // TODO: set interiorBegin, interiorEnd
 
       const int newRank = getRankOfProcCoord( this->subdomainCoordinates );
-      this->communicator = MPI::Comm::split( oldCommunicator, 1, newRank );
+      this->communicator = oldCommunicator.split( 1, newRank );
 
       setupNeighbors();
 
@@ -396,7 +396,7 @@ DistributedMesh< Grid< Dimension, Real, Device, Index > >::SetupByCut(
       return true;
    }
    else {
-      this->communicator = MPI::Comm::split( oldCommunicator, MPI_UNDEFINED, 0 );
+      this->communicator = oldCommunicator.split( MPI_UNDEFINED, 0 );
       return false;
    }
 }
diff --git a/src/TNL/Meshes/DistributedMeshes/DistributedGridSynchronizer.h b/src/TNL/Meshes/DistributedMeshes/DistributedGridSynchronizer.h
index 2f258271f166d9477c601a3cf095222dec546316..ad67a8a25236c97a2c96691a6f52e75ea715d5b6 100644
--- a/src/TNL/Meshes/DistributedMeshes/DistributedGridSynchronizer.h
+++ b/src/TNL/Meshes/DistributedMeshes/DistributedGridSynchronizer.h
@@ -145,7 +145,7 @@ public:
 
       // async send and receive
       std::unique_ptr< MPI_Request[] > requests{ new MPI_Request[ 2 * this->getNeighborsCount() ] };
-      MPI_Comm communicator = distributedGrid->getCommunicator();
+      const MPI::Comm& communicator = distributedGrid->getCommunicator();
       int requestsCount( 0 );
 
       // send everything, recieve everything
diff --git a/src/TNL/Meshes/DistributedMeshes/DistributedMesh.h b/src/TNL/Meshes/DistributedMeshes/DistributedMesh.h
index 6958288730530f768cf3b202b4513677d4906a3f..2c9c04a2376ac70b1497ddb2f4d703acfb54ee21 100644
--- a/src/TNL/Meshes/DistributedMeshes/DistributedMesh.h
+++ b/src/TNL/Meshes/DistributedMeshes/DistributedMesh.h
@@ -95,9 +95,9 @@ public:
     * Methods specific to the distributed mesh
     */
    void
-   setCommunicator( MPI::Comm&& communicator )
+   setCommunicator( const MPI::Comm& communicator )
    {
-      this->communicator = std::move( communicator );
+      this->communicator = communicator;
    }
 
    const MPI::Comm&
diff --git a/src/TNL/Meshes/DistributedMeshes/DistributedMeshSynchronizer.h b/src/TNL/Meshes/DistributedMeshes/DistributedMeshSynchronizer.h
index 45b84a9b1ed6799840a13b998adfe98e4ff8768a..04ac62ad96f7ddacce9e8e9ad17e4dfeaf6e71ce 100644
--- a/src/TNL/Meshes/DistributedMeshes/DistributedMeshSynchronizer.h
+++ b/src/TNL/Meshes/DistributedMeshes/DistributedMeshSynchronizer.h
@@ -12,6 +12,7 @@
 #include <TNL/Containers/ByteArraySynchronizer.h>
 #include <TNL/Containers/Vector.h>
 #include <TNL/Matrices/DenseMatrix.h>
+#include <TNL/MPI/Comm.h>
 #include <TNL/MPI/Wrappers.h>
 
 namespace TNL {
@@ -58,8 +59,8 @@ public:
                      "Global indices are not allocated properly." );
 
       communicator = mesh.getCommunicator();
-      const int rank = MPI::GetRank( communicator );
-      const int nproc = MPI::GetSize( communicator );
+      const int rank = communicator.rank();
+      const int nproc = communicator.size();
 
       // exchange the global index offsets so that each rank can determine the
       // owner of every entity by its global index
@@ -207,8 +208,8 @@ public:
                      bytesPerValue * ghostOffsets[ ghostOffsets.getSize() - 1 ],
                      "The array does not have the expected size." );
 
-      const int rank = MPI::GetRank( communicator );
-      const int nproc = MPI::GetSize( communicator );
+      const int rank = communicator.rank();
+      const int nproc = communicator.size();
 
       // allocate send buffers (setSize does nothing if the array size is already correct)
       sendBuffers.setSize( bytesPerValue * ghostNeighborOffsets[ nproc ] );
@@ -271,8 +272,8 @@ public:
    {
       TNL_ASSERT_EQ( pattern.getRows(), ghostOffsets[ ghostOffsets.getSize() - 1 ], "invalid sparse pattern matrix" );
 
-      const int rank = MPI::GetRank( communicator );
-      const int nproc = MPI::GetSize( communicator );
+      const int rank = communicator.rank();
+      const int nproc = communicator.size();
 
       // buffer for asynchronous communication requests
       RequestsVector requests;
@@ -461,7 +462,7 @@ public:
 
 protected:
    // communicator taken from the distributed mesh
-   MPI_Comm communicator;
+   MPI::Comm communicator;
 
    /**
     * Global offsets: array of size nproc where the i-th value is the lowest
diff --git a/src/TNL/Meshes/Readers/PVTIReader.h b/src/TNL/Meshes/Readers/PVTIReader.h
index 44e2a3ecc07cf9a6cb674ffa32e806cec8c59e8f..233ea00196506a5ff33f65806beeb6e97d7cfa7b 100644
--- a/src/TNL/Meshes/Readers/PVTIReader.h
+++ b/src/TNL/Meshes/Readers/PVTIReader.h
@@ -10,7 +10,7 @@
 
 #include <experimental/filesystem>
 
-#include <TNL/MPI/Wrappers.h>
+#include <TNL/MPI/Comm.h>
 #include <TNL/MPI/Utils.h>
 #include <TNL/Meshes/Readers/VTIReader.h>
 #include <TNL/Meshes/MeshDetails/layers/EntityTags/Traits.h>
@@ -137,14 +137,14 @@ class PVTIReader : public XMLVTK
          throw MeshReaderError( "PVTIReader", "the file does not contain any <Piece> element." );
 
       // check that the number of pieces matches the number of MPI ranks
-      const int nproc = MPI::GetSize( communicator );
+      const int nproc = communicator.size();
       if( (int) pieceSources.size() != nproc )
          throw MeshReaderError( "PVTIReader",
                                 "the number of subdomains does not match the number of MPI ranks ("
                                    + std::to_string( pieceSources.size() ) + " vs " + std::to_string( nproc ) + ")." );
 
       // read the local piece source
-      const int rank = MPI::GetRank( communicator );
+      const int rank = communicator.rank();
       localReader.setFileName( pieceSources[ rank ] );
       localReader.detectMesh();
 
@@ -174,8 +174,8 @@ class PVTIReader : public XMLVTK
 public:
    PVTIReader() = default;
 
-   PVTIReader( const std::string& fileName, MPI_Comm communicator = MPI_COMM_WORLD )
-   : XMLVTK( fileName ), communicator( communicator )
+   PVTIReader( const std::string& fileName, MPI::Comm communicator = MPI_COMM_WORLD )
+   : XMLVTK( fileName ), communicator( std::move( communicator ) )
    {}
 
    void
@@ -249,7 +249,7 @@ public:
       localReader.loadMesh( localMesh );
       if( localMesh != mesh.getLocalMesh() ) {
          std::stringstream msg;
-         msg << "The grid from the " << MPI::GetRank( communicator )
+         msg << "The grid from the " << communicator.rank()
              << "-th subdomain .vti file does not match the local grid of the DistributedGrid."
              << "\n- Grid from the .vti file:\n"
              << localMesh << "\n- Local grid from the DistributedGrid:\n"
@@ -323,7 +323,7 @@ public:
    }
 
 protected:
-   MPI_Comm communicator;
+   MPI::Comm communicator;
 
    int ghostLevels = 0;
    int minCommonVertices = 0;
diff --git a/src/TNL/Meshes/Readers/PVTUReader.h b/src/TNL/Meshes/Readers/PVTUReader.h
index f3c72d3e19a8cd263d4775486c235ab44564ec6f..bf2b3236109c2e3a937c09218189ae3cb560582f 100644
--- a/src/TNL/Meshes/Readers/PVTUReader.h
+++ b/src/TNL/Meshes/Readers/PVTUReader.h
@@ -9,8 +9,9 @@
 #pragma once
 
 #include <experimental/filesystem>
+#include <utility>
 
-#include <TNL/MPI/Wrappers.h>
+#include <TNL/MPI/Comm.h>
 #include <TNL/MPI/Utils.h>
 #include <TNL/Meshes/Readers/VTUReader.h>
 #include <TNL/Meshes/MeshDetails/layers/EntityTags/Traits.h>
@@ -65,14 +66,14 @@ class PVTUReader : public XMLVTK
          throw MeshReaderError( "PVTUReader", "the file does not contain any <Piece> element." );
 
       // check that the number of pieces matches the number of MPI ranks
-      const int nproc = MPI::GetSize( communicator );
+      const int nproc = communicator.size();
       if( (int) pieceSources.size() != nproc )
          throw MeshReaderError( "PVTUReader",
                                 "the number of subdomains does not match the number of MPI ranks ("
                                    + std::to_string( pieceSources.size() ) + " vs " + std::to_string( nproc ) + ")." );
 
       // read the local piece source
-      const int rank = MPI::GetRank( communicator );
+      const int rank = communicator.rank();
       localReader.setFileName( pieceSources[ rank ] );
       localReader.detectMesh();
 
@@ -101,8 +102,8 @@ class PVTUReader : public XMLVTK
 public:
    PVTUReader() = default;
 
-   PVTUReader( const std::string& fileName, MPI_Comm communicator = MPI_COMM_WORLD )
-   : XMLVTK( fileName ), communicator( communicator )
+   PVTUReader( const std::string& fileName, MPI::Comm communicator = MPI_COMM_WORLD )
+   : XMLVTK( fileName ), communicator( std::move( communicator ) )
    {}
 
    void
@@ -228,7 +229,7 @@ public:
       if( minCount == 0 ) {
          // split the communicator, remove the ranks which did not get a subdomain
          const int color = ( pointsCount > 0 && cellsCount > 0 ) ? 0 : MPI_UNDEFINED;
-         MPI::Comm subCommunicator = MPI::Comm::split( communicator, color, 0 );
+         MPI::Comm subCommunicator = communicator.split( color, 0 );
 
          // set the communicator
          mesh.setCommunicator( std::move( subCommunicator ) );
@@ -262,7 +263,7 @@ public:
    }
 
 protected:
-   MPI_Comm communicator;
+   MPI::Comm communicator;
 
    int ghostLevels = 0;
    int minCommonVertices = 0;
diff --git a/src/TNL/Meshes/Writers/PVTIWriter.hpp b/src/TNL/Meshes/Writers/PVTIWriter.hpp
index d4136364d50e888a8a90e44aa70d5f69e8fc7397..f0c2ad2a6156d5b95af4ea78f53717d6c8fe19d7 100644
--- a/src/TNL/Meshes/Writers/PVTIWriter.hpp
+++ b/src/TNL/Meshes/Writers/PVTIWriter.hpp
@@ -200,13 +200,13 @@ std::string
 PVTIWriter< Grid >::addPiece( const std::string& mainFileName,
                               const DistributedMeshes::DistributedMesh< Grid >& distributedMesh )
 {
-   const MPI_Comm communicator = distributedMesh.getCommunicator();
+   const MPI::Comm& communicator = distributedMesh.getCommunicator();
    const typename Grid::CoordinatesType& globalBegin = distributedMesh.getGlobalBegin() - distributedMesh.getLowerOverlap();
    const typename Grid::CoordinatesType& globalEnd =
       globalBegin + distributedMesh.getLocalSize() + distributedMesh.getUpperOverlap();
 
    // exchange globalBegin and globalEnd among the ranks
-   const int nproc = MPI::GetSize( communicator );
+   const int nproc = communicator.size();
    std::unique_ptr< typename Grid::CoordinatesType[] > beginsForScatter{ new typename Grid::CoordinatesType[ nproc ] };
    std::unique_ptr< typename Grid::CoordinatesType[] > endsForScatter{ new typename Grid::CoordinatesType[ nproc ] };
    for( int i = 0; i < nproc; i++ ) {
@@ -231,9 +231,9 @@ PVTIWriter< Grid >::addPiece( const std::string& mainFileName,
 
    // add pieces for all ranks, return the source for the current rank
    std::string source;
-   for( int i = 0; i < MPI::GetSize( communicator ); i++ ) {
+   for( int i = 0; i < communicator.size(); i++ ) {
       const std::string s = addPiece( mainFileName, i, globalBegins[ i ], globalEnds[ i ] );
-      if( i == MPI::GetRank( communicator ) )
+      if( i == communicator.rank() )
          source = s;
    }
    return source;
diff --git a/src/TNL/Solvers/Linear/Utils/Traits.h b/src/TNL/Solvers/Linear/Utils/Traits.h
index a94298036e980fc12c6efc95f48d3f3728bb07c5..86bdc225cae8492df1c2b046ae74d889820daed5 100644
--- a/src/TNL/Solvers/Linear/Utils/Traits.h
+++ b/src/TNL/Solvers/Linear/Utils/Traits.h
@@ -8,7 +8,7 @@
 
 #pragma once
 
-#include <TNL/MPI/Wrappers.h>
+#include <TNL/MPI/Comm.h>
 #include <TNL/Containers/Vector.h>
 #include <TNL/Containers/VectorView.h>
 #include <TNL/Containers/DistributedVector.h>
@@ -50,7 +50,7 @@ struct Traits
       return v;
    }
 
-   static MPI_Comm
+   static MPI::Comm
    getCommunicator( const Matrix& m )
    {
       return MPI_COMM_WORLD;
@@ -98,7 +98,7 @@ struct Traits< Matrices::DistributedMatrix< Matrix > >
       return v.getLocalView();
    }
 
-   static MPI_Comm
+   static const MPI::Comm&
    getCommunicator( const Matrices::DistributedMatrix< Matrix >& m )
    {
       return m.getCommunicator();