diff --git a/Applications/ApplicationsLib/LinearSolverLibrarySetup.h b/Applications/ApplicationsLib/LinearSolverLibrarySetup.h
index 59c94cd317dd620fda5c1493ca00edba7bce144b..e6b48ce2f6ad7506099b3811ca8fd604fe6b47f3 100644
--- a/Applications/ApplicationsLib/LinearSolverLibrarySetup.h
+++ b/Applications/ApplicationsLib/LinearSolverLibrarySetup.h
@@ -19,6 +19,7 @@
 /// The default implementation is empty providing polymorphic behaviour when
 /// using this class.
 
+#include "BaseLib/MPI.h"
 #include "NumLib/DOF/GlobalMatrixProviders.h"
 
 #if defined(USE_PETSC)
@@ -28,9 +29,8 @@ namespace ApplicationsLib
 {
 struct LinearSolverLibrarySetup final
 {
-    LinearSolverLibrarySetup(int argc, char* argv[])
+    LinearSolverLibrarySetup(int argc, char* argv[]) : mpi_setup(argc, argv)
     {
-        MPI_Init(&argc, &argv);
         char help[] = "ogs6 with PETSc \n";
         PetscInitialize(&argc, &argv, nullptr, help);
         MPI_Comm_set_errhandler(PETSC_COMM_WORLD, MPI_ERRORS_RETURN);
@@ -40,8 +40,9 @@ struct LinearSolverLibrarySetup final
     {
         NumLib::cleanupGlobalMatrixProviders();
         PetscFinalize();
-        MPI_Finalize();
     }
+
+    BaseLib::MPI::Setup mpi_setup;
 };
 }    // ApplicationsLib
 #elif defined(USE_LIS)
diff --git a/Applications/ApplicationsLib/TestDefinition.cpp b/Applications/ApplicationsLib/TestDefinition.cpp
index 6c936e4eff30dcbadc6db6e0ff1a974d79899d27..43856010bb4fa0f683fce3baa255ded5ee573a38 100644
--- a/Applications/ApplicationsLib/TestDefinition.cpp
+++ b/Applications/ApplicationsLib/TestDefinition.cpp
@@ -20,9 +20,9 @@
 #include "BaseLib/ConfigTree.h"
 #include "BaseLib/Error.h"
 #include "BaseLib/FileTools.h"
-#ifdef USE_PETSC
-#include <petsc.h>
+#include "BaseLib/MPI.h"
 
+#ifdef USE_PETSC
 #include "MeshLib/IO/VtkIO/VtuInterface.h"  // For petsc file name conversion.
 #endif
 
@@ -197,16 +197,13 @@ TestDefinition::TestDefinition(BaseLib::ConfigTree const& config_tree,
                 //! \ogs_file_param{prj__test_definition__vtkdiff__file}
                 vtkdiff_config.getConfigParameter<std::string>("file");
 #ifdef USE_PETSC
-            int mpi_size;
-            MPI_Comm_size(PETSC_COMM_WORLD, &mpi_size);
-            if (mpi_size > 1)
+            BaseLib::MPI::Mpi mpi;
+            if (mpi.size > 1)
             {
-                int rank;
-                MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
                 filename =
                     MeshLib::IO::getVtuFileNameForPetscOutputWithoutExtension(
                         filename) +
-                    "_" + std::to_string(rank) + ".vtu";
+                    "_" + std::to_string(mpi.rank) + ".vtu";
             }
 #endif  // OGS_USE_PETSC
             filenames.push_back(filename);
diff --git a/Applications/Utils/FileConverter/ConvertSHPToGLI.cpp b/Applications/Utils/FileConverter/ConvertSHPToGLI.cpp
index 2fb39d512d68b563c2384c9ab34ab0f7eac93320..6023763a1b083f5aeeb3ab264eca6c76faab5ed4 100644
--- a/Applications/Utils/FileConverter/ConvertSHPToGLI.cpp
+++ b/Applications/Utils/FileConverter/ConvertSHPToGLI.cpp
@@ -14,10 +14,6 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 // STL
 #include <fstream>
 #include <vector>
@@ -25,6 +21,7 @@
 // ShapeLib
 #include <shapefil.h>
 
+#include "BaseLib/MPI.h"
 #include "GeoLib/GEOObjects.h"
 #include "GeoLib/IO/XmlIO/Qt/XmlGmlInterface.h"
 #include "GeoLib/IO/XmlIO/Qt/XmlStnInterface.h"
@@ -185,9 +182,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::string fname(shapefile_arg.getValue());
 
@@ -207,9 +202,6 @@ int main(int argc, char* argv[])
             ERR("Shape file contains {:d} polylines.", number_of_elements);
             ERR("This programme only handles only files containing points.");
             SHPClose(hSHP);
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_SUCCESS;
         }
         SHPClose(hSHP);
@@ -304,8 +296,5 @@ int main(int argc, char* argv[])
         ERR("Could not open the database file.");
     }
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/FileConverter/FEFLOW2OGS.cpp b/Applications/Utils/FileConverter/FEFLOW2OGS.cpp
index acd33c06c8998202652ef1f1fdc85b3fb216e0c2..e299c2b56d2b6cfe1c98f0b0247a0d269f2008c4 100644
--- a/Applications/Utils/FileConverter/FEFLOW2OGS.cpp
+++ b/Applications/Utils/FileConverter/FEFLOW2OGS.cpp
@@ -14,12 +14,8 @@
 // ThirdParty
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
-// BaseLib
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "BaseLib/RunTime.h"
 #include "InfoLib/GitInfo.h"
 #ifndef WIN32
@@ -62,9 +58,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     // *** read mesh
     INFO("Reading {:s}.", feflow_mesh_arg.getValue());
@@ -81,9 +75,6 @@ int main(int argc, char* argv[])
     if (mesh == nullptr)
     {
         INFO("Could not read mesh from {:s}.", feflow_mesh_arg.getValue());
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 #ifndef WIN32
@@ -99,8 +90,5 @@ int main(int argc, char* argv[])
     INFO("Writing {:s}.", ogs_mesh_fname);
     MeshLib::IO::writeMeshToFile(*mesh, ogs_mesh_fname);
     INFO("\tDone.");
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/FileConverter/GMSH2OGS.cpp b/Applications/Utils/FileConverter/GMSH2OGS.cpp
index a528076749b32809089eaf5bf03ad9a226cb3c3f..6674dbc902fdf8828b4626c5c575ca517f782fdd 100644
--- a/Applications/Utils/FileConverter/GMSH2OGS.cpp
+++ b/Applications/Utils/FileConverter/GMSH2OGS.cpp
@@ -19,12 +19,8 @@
 // ThirdParty
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
-// BaseLib
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "BaseLib/RunTime.h"
 #include "InfoLib/GitInfo.h"
 #ifndef WIN32
@@ -176,9 +172,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     // *** read mesh
     INFO("Reading {:s}.", gmsh_mesh_arg.getValue());
@@ -194,9 +188,6 @@ int main(int argc, char* argv[])
     if (mesh == nullptr)
     {
         INFO("Could not read mesh from {:s}.", gmsh_mesh_arg.getValue());
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return -1;
     }
 #ifndef WIN32
@@ -279,8 +270,5 @@ int main(int argc, char* argv[])
     MeshLib::IO::writeMeshToFile(*mesh, ogs_mesh_arg.getValue());
 
     delete mesh;
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/FileConverter/GocadSGridReader.cpp b/Applications/Utils/FileConverter/GocadSGridReader.cpp
index 0f3e6784c141519e3c80ef9a3b96d3c416148c28..a28124cc5c313cf0e8e36b78d718179614f285bb 100644
--- a/Applications/Utils/FileConverter/GocadSGridReader.cpp
+++ b/Applications/Utils/FileConverter/GocadSGridReader.cpp
@@ -12,16 +12,13 @@
 #include <spdlog/spdlog.h>
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <fstream>
 #include <sstream>
 #include <string>
 
 #include "Applications/FileIO/GocadIO/GenerateFaceSetMeshes.h"
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/Elements/Element.h"
 #include "MeshLib/IO/writeMeshToFile.h"
@@ -61,9 +58,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     // read the Gocad SGrid
     INFO("Start reading Gocad SGrid.");
@@ -82,8 +77,5 @@ int main(int argc, char* argv[])
 
     INFO("Writing mesh to '{:s}'.", mesh_output_arg.getValue());
     MeshLib::IO::writeMeshToFile(*mesh, mesh_output_arg.getValue());
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/FileConverter/GocadTSurfaceReader.cpp b/Applications/Utils/FileConverter/GocadTSurfaceReader.cpp
index dd93c5024e1aa173a98d928661f29279c39fd536..8bb1391f3e41f3e5af1d40d23f7cabe88ac010f1 100644
--- a/Applications/Utils/FileConverter/GocadTSurfaceReader.cpp
+++ b/Applications/Utils/FileConverter/GocadTSurfaceReader.cpp
@@ -9,11 +9,8 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include "Applications/FileIO/GocadIO/GocadAsciiReader.h"
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/VtkIO/VtuInterface.h"
 #include "MeshLib/Mesh.h"
@@ -65,17 +62,12 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     if (export_lines_arg.isSet() && export_surfaces_arg.isSet())
     {
         ERR("Both the 'lines-only'-flag and 'surfaces-only'-flag are set. Only "
             "one is allowed at a time.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return 2;
     }
 
@@ -94,9 +86,6 @@ int main(int argc, char* argv[])
     if (!FileIO::Gocad::GocadAsciiReader::readFile(file_name, meshes, t))
     {
         ERR("Error reading file.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return 1;
     }
     INFO("{:d} meshes found.", meshes.size());
@@ -115,8 +104,5 @@ int main(int argc, char* argv[])
         MeshLib::IO::VtuInterface vtu(mesh.get(), data_mode, compressed);
         vtu.writeToFile(dir + delim + mesh->getName() + ".vtu");
     }
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return 0;
 }
diff --git a/Applications/Utils/FileConverter/Mesh2Raster.cpp b/Applications/Utils/FileConverter/Mesh2Raster.cpp
index 4e32f111d57e707572a6fb0c5ac85e8f0fb556c4..9b090554776f0416475b75b7fa709e443a1fcd77 100644
--- a/Applications/Utils/FileConverter/Mesh2Raster.cpp
+++ b/Applications/Utils/FileConverter/Mesh2Raster.cpp
@@ -9,15 +9,12 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <filesystem>
 #include <fstream>
 #include <memory>
 #include <string>
 
+#include "BaseLib/MPI.h"
 #include "GeoLib/AABB.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/readMeshFromFile.h"
@@ -52,9 +49,7 @@ int main(int argc, char* argv[])
     cmd.add(input_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     INFO("Rasterising mesh...");
     std::unique_ptr<MeshLib::Mesh> const mesh(
@@ -62,18 +57,12 @@ int main(int argc, char* argv[])
     if (mesh == nullptr)
     {
         ERR("Error reading mesh file.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return 1;
     }
     if (mesh->getDimension() != 2)
     {
         ERR("The programme requires a mesh containing two-dimensional elements "
             "(i.e. triangles or quadrilaterals.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return 2;
     }
 
@@ -178,8 +167,5 @@ int main(int argc, char* argv[])
     }
     out.close();
     INFO("Result written to {:s}", output_name);
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return 0;
 }
diff --git a/Applications/Utils/FileConverter/Mesh2Shape.cpp b/Applications/Utils/FileConverter/Mesh2Shape.cpp
index b9ffd3f0493f43684b74c9d7fce4982f2e43840c..44bb45b383b51c4a733d30b8df1883fce431d38a 100644
--- a/Applications/Utils/FileConverter/Mesh2Shape.cpp
+++ b/Applications/Utils/FileConverter/Mesh2Shape.cpp
@@ -9,11 +9,8 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include "Applications/FileIO/SHPInterface.h"
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/readMeshFromFile.h"
 #include "MeshLib/Mesh.h"
@@ -43,22 +40,14 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::string const file_name(input_arg.getValue());
     std::unique_ptr<MeshLib::Mesh> const mesh(
         MeshLib::IO::readMeshFromFile(file_name));
     if (FileIO::SHPInterface::write2dMeshToSHP(output_arg.getValue(), *mesh))
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_SUCCESS;
     }
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_FAILURE;
 }
diff --git a/Applications/Utils/FileConverter/NetCdfConverter.cpp b/Applications/Utils/FileConverter/NetCdfConverter.cpp
index 84e3be1b7ab368f37932339a0ab59c6c7906f93c..8009dd19e5f1553d900ee6c375a3f401f872ceb6 100644
--- a/Applications/Utils/FileConverter/NetCdfConverter.cpp
+++ b/Applications/Utils/FileConverter/NetCdfConverter.cpp
@@ -9,10 +9,6 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 // STL
 #include <cctype>
 #include <iostream>
@@ -26,6 +22,7 @@
 
 #include "BaseLib/FileTools.h"
 #include "BaseLib/Logging.h"
+#include "BaseLib/MPI.h"
 #include "GeoLib/IO/AsciiRasterInterface.h"
 #include "GeoLib/Raster.h"
 #include "InfoLib/GitInfo.h"
@@ -740,18 +737,13 @@ int main(int argc, char* argv[])
     cmd.add(arg_input);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     NcFile dataset(arg_input.getValue().c_str(), NcFile::read);
 
     if (dataset.isNull())
     {
         ERR("Error opening file.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return -1;
     }
 
@@ -761,9 +753,6 @@ int main(int argc, char* argv[])
     {
         ERR("Only one output format can be specified (single-file, multi-file, "
             "or images)");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -783,9 +772,6 @@ int main(int argc, char* argv[])
     if (var.isNull())
     {
         ERR("Variable \"{:s}\" not found in file.", arg_varname.getValue());
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -797,9 +783,6 @@ int main(int argc, char* argv[])
         if (!assignDimParams(var, dim_idx_map, arg_dim_time, arg_dim1, arg_dim2,
                              arg_dim3))
         {
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
     }
@@ -856,15 +839,9 @@ int main(int argc, char* argv[])
     if (!convert(dataset, var, output_name, dim_idx_map, is_time_dep,
                  time_bounds, output, elem_type))
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
     std::cout << "Conversion finished successfully.\n";
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/FileConverter/OGS2VTK.cpp b/Applications/Utils/FileConverter/OGS2VTK.cpp
index 78f360c28831677e5674e49776b26366f2517e62..46bc603a270d4bc555575f66305ed4d88d67c131 100644
--- a/Applications/Utils/FileConverter/OGS2VTK.cpp
+++ b/Applications/Utils/FileConverter/OGS2VTK.cpp
@@ -13,13 +13,10 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <memory>
 #include <string>
 
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/VtkIO/VtuInterface.h"
 #include "MeshLib/IO/readMeshFromFile.h"
@@ -52,17 +49,12 @@ int main(int argc, char* argv[])
     cmd.add(use_ascii_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::unique_ptr<MeshLib::Mesh const> mesh(
         MeshLib::IO::readMeshFromFile(mesh_in.getValue()));
     if (!mesh)
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     INFO("Mesh read: {:d} nodes, {:d} elements.", mesh->getNumberOfNodes(),
@@ -73,8 +65,5 @@ int main(int argc, char* argv[])
 
     MeshLib::IO::writeVtu(*mesh, mesh_out.getValue(), data_mode);
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/FileConverter/PVD2XDMF.cpp b/Applications/Utils/FileConverter/PVD2XDMF.cpp
index 838aa202559ec5b4e0aefe8870fc7c8d3273bbdb..6c114ef6673ced9e09f7d28a6658af1b2f01f796 100644
--- a/Applications/Utils/FileConverter/PVD2XDMF.cpp
+++ b/Applications/Utils/FileConverter/PVD2XDMF.cpp
@@ -9,10 +9,6 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <array>
 #include <boost/property_tree/ptree.hpp>
 #include <boost/property_tree/xml_parser.hpp>
@@ -20,6 +16,7 @@
 
 #include "BaseLib/FileTools.h"
 #include "BaseLib/Logging.h"
+#include "BaseLib/MPI.h"
 #include "BaseLib/MemWatch.h"
 #include "BaseLib/RunTime.h"
 #include "BaseLib/StringTools.h"
@@ -141,9 +138,7 @@ int main(int argc, char* argv[])
     cmd.add(outdir_arg);
 
     cmd.parse(argc, argv);
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
     BaseLib::setConsoleLogLevel(log_level_arg.getValue());
 
     auto const pvd_file_dir = BaseLib::extractPath(pvd_file_arg.getValue());
@@ -237,8 +232,5 @@ int main(int argc, char* argv[])
 
         mesh_xdmf_hdf_writer->writeStep(time);
     }
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/FileConverter/Raster2ASC.cpp b/Applications/Utils/FileConverter/Raster2ASC.cpp
index de81adf7d1929dd4b8a32d491478d7f6b257d942..398d8665c11a6ef2099cf069e6b096019f20175a 100644
--- a/Applications/Utils/FileConverter/Raster2ASC.cpp
+++ b/Applications/Utils/FileConverter/Raster2ASC.cpp
@@ -11,10 +11,7 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
+#include "BaseLib/MPI.h"
 #include "GeoLib/IO/AsciiRasterInterface.h"
 #include "GeoLib/Raster.h"
 #include "InfoLib/GitInfo.h"
@@ -41,9 +38,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::unique_ptr<GeoLib::Raster> raster(
         FileIO::AsciiRasterInterface::readRaster(input_arg.getValue()));
@@ -61,8 +56,5 @@ int main(int argc, char* argv[])
     }
 
     FileIO::AsciiRasterInterface::writeRasterAsASC(*raster, output_name);
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/FileConverter/TIN2VTK.cpp b/Applications/Utils/FileConverter/TIN2VTK.cpp
index 91cc0a1962bf41ebe3f606f76122d14027b97cac..1d2b7290baabbb87ab31c838a49758dd55131d2c 100644
--- a/Applications/Utils/FileConverter/TIN2VTK.cpp
+++ b/Applications/Utils/FileConverter/TIN2VTK.cpp
@@ -9,17 +9,13 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 // STL
 #include <memory>
 #include <string>
 #include <vector>
 
-// BaseLib
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 
 // GeoLib
@@ -53,9 +49,7 @@ int main(int argc, char* argv[])
     cmd.add(outArg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     INFO("reading the TIN file...");
     const std::string tinFileName(inArg.getValue());
@@ -66,9 +60,6 @@ int main(int argc, char* argv[])
         GeoLib::IO::TINInterface::readTIN(tinFileName, point_vec));
     if (!sfc)
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     INFO("TIN read:  {:d} points, {:d} triangles", point_vec.size(),
@@ -85,8 +76,5 @@ int main(int argc, char* argv[])
     MeshLib::IO::VtuInterface writer(mesh.get());
     writer.writeToFile(outArg.getValue());
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/FileConverter/TecPlotTools.cpp b/Applications/Utils/FileConverter/TecPlotTools.cpp
index 0d70176cf1972bd2893645627e388d9a0ff5b391..bf209269396a1bb33f9ef92cd83b66116a705e63 100644
--- a/Applications/Utils/FileConverter/TecPlotTools.cpp
+++ b/Applications/Utils/FileConverter/TecPlotTools.cpp
@@ -9,16 +9,13 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <iostream>
 #include <memory>
 #include <string>
 #include <string_view>
 #include <vector>
 
+#include "BaseLib/MPI.h"
 #include "BaseLib/StringTools.h"
 #include "GeoLib/Point.h"
 #include "InfoLib/GitInfo.h"
@@ -480,25 +477,17 @@ int main(int argc, char* argv[])
     cmd.add(input_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     if (!input_arg.isSet())
     {
         ERR("No input file given. Please specify TecPlot (*.plt) file");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return -1;
     }
 
     if (convert_arg.getValue() && !output_arg.isSet())
     {
         ERR("No output file given. Please specify OGS mesh (*.vtu) file");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return -1;
     }
 
@@ -506,18 +495,12 @@ int main(int argc, char* argv[])
     if (!in.is_open())
     {
         ERR("Could not open file {:s}.", input_arg.getValue());
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return -2;
     }
 
     if (!convert_arg.isSet() && !split_arg.isSet())
     {
         INFO("Nothing to do. Use -s to split or -c to convert.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return 0;
     }
 
@@ -534,8 +517,5 @@ int main(int argc, char* argv[])
     }
 
     in.close();
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return return_val;
 }
diff --git a/Applications/Utils/FileConverter/VTK2OGS.cpp b/Applications/Utils/FileConverter/VTK2OGS.cpp
index 2ec723dc2e79bf35f90701f6a7e72441c7d50a12..80c424df8a1e37bd1a45d44f717755352e4dcc81 100644
--- a/Applications/Utils/FileConverter/VTK2OGS.cpp
+++ b/Applications/Utils/FileConverter/VTK2OGS.cpp
@@ -14,12 +14,9 @@
 // STL
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <string>
 
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/Legacy/MeshIO.h"
 #include "MeshLib/IO/VtkIO/VtuInterface.h"
@@ -47,9 +44,7 @@ int main(int argc, char* argv[])
     cmd.add(mesh_out);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     MeshLib::Mesh* mesh(
         MeshLib::IO::VtuInterface::readVTUFile(mesh_in.getValue()));
@@ -60,8 +55,5 @@ int main(int argc, char* argv[])
     meshIO.setMesh(mesh);
     BaseLib::IO::writeStringToFile(meshIO.writeToString(), mesh_out.getValue());
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/FileConverter/VTK2TIN.cpp b/Applications/Utils/FileConverter/VTK2TIN.cpp
index b0401d1b199b66373e9b1a76e1518c38e5f3f640..bde167facc2857c2d42ff72106a214945206a3b0 100644
--- a/Applications/Utils/FileConverter/VTK2TIN.cpp
+++ b/Applications/Utils/FileConverter/VTK2TIN.cpp
@@ -9,17 +9,13 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 // STL
 #include <fstream>
 #include <memory>
 #include <string>
 
-// BaseLib
 #include "BaseLib/Logging.h"
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 
 // GeoLib
@@ -56,9 +52,7 @@ int main(int argc, char* argv[])
     cmd.add(mesh_out);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
     std::unique_ptr<MeshLib::Mesh> mesh(
         MeshLib::IO::VtuInterface::readVTUFile(mesh_in.getValue()));
     INFO("Mesh read: {:d} nodes, {:d} elements.", mesh->getNumberOfNodes(),
@@ -74,8 +68,5 @@ int main(int argc, char* argv[])
             mesh_out.getValue());
     }
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/FileConverter/convertGEO.cpp b/Applications/Utils/FileConverter/convertGEO.cpp
index c7f76b7f641356e020ac363d3d3e341a8c62a6f9..3a7d859d519b1299aeb8c17471a887feae2dc2dd 100644
--- a/Applications/Utils/FileConverter/convertGEO.cpp
+++ b/Applications/Utils/FileConverter/convertGEO.cpp
@@ -10,15 +10,12 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <string>
 #include <vector>
 
 #include "Applications/FileIO/readGeometryFromFile.h"
 #include "Applications/FileIO/writeGeometryToFile.h"
+#include "BaseLib/MPI.h"
 #include "GeoLib/GEOObjects.h"
 #include "InfoLib/GitInfo.h"
 
@@ -51,9 +48,7 @@ int main(int argc, char* argv[])
     cmd.add(gmsh_path_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     GeoLib::GEOObjects geoObjects;
     FileIO::readGeometryFromFile(argInputFileName.getValue(), geoObjects,
@@ -64,8 +59,5 @@ int main(int argc, char* argv[])
     FileIO::writeGeometryToFile(geo_names[0], geoObjects,
                                 argOutputFileName.getValue());
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/FileConverter/generateMatPropsFromMatID.cpp b/Applications/Utils/FileConverter/generateMatPropsFromMatID.cpp
index 9291eacb3b7bdbfa7c53b025b73bc72d7bde5d50..88032d6354ea69f369c9deafe262a90275eb4e39 100644
--- a/Applications/Utils/FileConverter/generateMatPropsFromMatID.cpp
+++ b/Applications/Utils/FileConverter/generateMatPropsFromMatID.cpp
@@ -14,14 +14,11 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <fstream>
 #include <memory>
 
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/Elements/Element.h"
 #include "MeshLib/IO/readMeshFromFile.h"
@@ -50,9 +47,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     // read mesh
     std::unique_ptr<MeshLib::Mesh> mesh(
@@ -61,9 +56,6 @@ int main(int argc, char* argv[])
     if (!mesh)
     {
         INFO("Could not read mesh from file '{:s}'.", mesh_arg.getValue());
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -92,9 +84,6 @@ int main(int argc, char* argv[])
     else
     {
         ERR("Could not create property '{:s}' file.", new_matname);
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -106,8 +95,5 @@ int main(int argc, char* argv[])
 
     INFO("New files '{:s}' and '{:s}' written.", new_mshname, new_matname);
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/GeoTools/MoveGeometry.cpp b/Applications/Utils/GeoTools/MoveGeometry.cpp
index 33505f1d59475065fed2c5305b8595b99f50f4bc..dd14e7210d22b9be1c2e9d7a59984b4ab4254a30 100644
--- a/Applications/Utils/GeoTools/MoveGeometry.cpp
+++ b/Applications/Utils/GeoTools/MoveGeometry.cpp
@@ -14,10 +14,7 @@
 // ThirdParty
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
+#include "BaseLib/MPI.h"
 #include "GeoLib/GEOObjects.h"
 #include "GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.h"
 #include "InfoLib/GitInfo.h"
@@ -49,9 +46,7 @@ int main(int argc, char* argv[])
     cmd.add(geo_input_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     GeoLib::GEOObjects geo_objects;
     GeoLib::IO::BoostXmlGmlInterface xml(geo_objects);
@@ -59,9 +54,6 @@ int main(int argc, char* argv[])
     {
         if (!xml.readFile(geo_input_arg.getValue()))
         {
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
     }
@@ -69,9 +61,6 @@ int main(int argc, char* argv[])
     {
         ERR("Failed to read file `{:s}'.", geo_input_arg.getValue());
         ERR("{:s}", err.what());
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -106,8 +95,5 @@ int main(int argc, char* argv[])
     BaseLib::IO::writeStringToFile(xml.writeToString(),
                                    geo_output_arg.getValue());
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/GeoTools/addDataToRaster.cpp b/Applications/Utils/GeoTools/addDataToRaster.cpp
index 425bfb5d68995ded4d0f5e2ed67ef65c81a03dcc..c78dc0549cc968b4fa5c77f25eb82545c8765116 100644
--- a/Applications/Utils/GeoTools/addDataToRaster.cpp
+++ b/Applications/Utils/GeoTools/addDataToRaster.cpp
@@ -11,16 +11,13 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <algorithm>
 #include <cmath>
 #include <memory>
 #include <numbers>
 #include <numeric>
 
+#include "BaseLib/MPI.h"
 #include "GeoLib/AABB.h"
 #include "GeoLib/IO/AsciiRasterInterface.h"
 #include "GeoLib/Point.h"
@@ -153,9 +150,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::array input_points = {
         GeoLib::Point{{ll_x_arg.getValue(), ll_y_arg.getValue(), 0}},
@@ -210,8 +205,5 @@ int main(int argc, char* argv[])
 
     FileIO::AsciiRasterInterface::writeRasterAsASC(*raster,
                                                    out_raster_arg.getValue());
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/GeoTools/createRaster.cpp b/Applications/Utils/GeoTools/createRaster.cpp
index 8d36fc335aa1e31bc76ffac6b557fa1aec426ac7..b3aa8b4b8f03bd945fde2a6f0e6111cf59dbd73b 100644
--- a/Applications/Utils/GeoTools/createRaster.cpp
+++ b/Applications/Utils/GeoTools/createRaster.cpp
@@ -11,12 +11,9 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <memory>
 
+#include "BaseLib/MPI.h"
 #include "GeoLib/AABB.h"
 #include "GeoLib/IO/AsciiRasterInterface.h"
 #include "GeoLib/Point.h"
@@ -71,9 +68,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     GeoLib::RasterHeader header{
         n_cols.getValue(),
@@ -88,8 +83,5 @@ int main(int argc, char* argv[])
     FileIO::AsciiRasterInterface::writeRasterAsASC(raster,
                                                    output_arg.getValue());
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/GeoTools/generateGeometry.cpp b/Applications/Utils/GeoTools/generateGeometry.cpp
index b84b167e37f1402c53aa0ac1113f331ad4d480ed..56c054e79e08d3304fdc27a8c4b97157bdf73327 100644
--- a/Applications/Utils/GeoTools/generateGeometry.cpp
+++ b/Applications/Utils/GeoTools/generateGeometry.cpp
@@ -11,12 +11,9 @@
 // ThirdParty
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <numeric>
 
+#include "BaseLib/MPI.h"
 #include "GeoLib/GEOObjects.h"
 #include "GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.h"
 #include "GeoLib/Point.h"
@@ -241,9 +238,7 @@ int main(int argc, char* argv[])
     cmd.add(geo_output_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     auto const p0 = GeoLib::Point{x0.getValue(), y0.getValue(), z0.getValue()};
     auto const p1 = GeoLib::Point{x1.getValue(), y1.getValue(), z1.getValue()};
@@ -286,9 +281,6 @@ int main(int argc, char* argv[])
                 std::move(polyline_name.getValue()), geometry_name.getValue(),
                 geometry) == EXIT_FAILURE)
         {
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
     }
@@ -298,8 +290,5 @@ int main(int argc, char* argv[])
     BaseLib::IO::writeStringToFile(xml.writeToString(),
                                    geo_output_arg.getValue());
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/AddElementQuality.cpp b/Applications/Utils/MeshEdit/AddElementQuality.cpp
index 25753925c288c2d753e9bf3b0b122affd6ca3730..08c7087db3b853aa5f9f0b787df975fc2bcf7c4e 100644
--- a/Applications/Utils/MeshEdit/AddElementQuality.cpp
+++ b/Applications/Utils/MeshEdit/AddElementQuality.cpp
@@ -9,13 +9,10 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <array>
 #include <string>
 
+#include "BaseLib/MPI.h"
 #include "BaseLib/RunTime.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/readMeshFromFile.h"
@@ -53,9 +50,7 @@ int main(int argc, char* argv[])
     cmd.add(mesh_in_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     // read the mesh file
     BaseLib::RunTime run_time;
@@ -64,9 +59,6 @@ int main(int argc, char* argv[])
         mesh_in_arg.getValue(), true /* compute_element_neighbors */));
     if (!mesh)
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     INFO("Time for reading: {:g} s", run_time.elapsed());
@@ -82,8 +74,5 @@ int main(int argc, char* argv[])
     INFO("Writing mesh '{:s}' ... ", mesh_out_arg.getValue());
     MeshLib::IO::writeMeshToFile(*mesh, mesh_out_arg.getValue());
     INFO("done.");
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/AddFaultToVoxelGrid.cpp b/Applications/Utils/MeshEdit/AddFaultToVoxelGrid.cpp
index e27be76f84654260c39335f24a183166d856f77e..2d4e47d92a53ddd7d2881760a84f021d4cf2751c 100644
--- a/Applications/Utils/MeshEdit/AddFaultToVoxelGrid.cpp
+++ b/Applications/Utils/MeshEdit/AddFaultToVoxelGrid.cpp
@@ -15,10 +15,7 @@
 // ThirdParty
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/VtkIO/VtuInterface.h"
 #include "MeshLib/IO/readMeshFromFile.h"
@@ -64,9 +61,7 @@ int main(int argc, char* argv[])
     cmd.add(input_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::string const input_name = input_arg.getValue();
     std::string const fault_name = fault_arg.getValue();
@@ -81,17 +76,11 @@ int main(int argc, char* argv[])
     if (mesh == nullptr)
     {
         ERR("Input mesh not found...");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     auto const& mat_ids = MeshLib::materialIDs(*mesh);
     if (!mat_ids)
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         ERR("Input mesh has no material IDs");
         return EXIT_FAILURE;
     }
@@ -105,15 +94,9 @@ int main(int argc, char* argv[])
     {
         MeshLib::IO::VtuInterface vtu(mesh.get());
         vtu.writeToFile(output_name);
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         INFO("The fault was successfully added.");
         return EXIT_SUCCESS;
     }
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     ERR("No fault could be added.");
     return EXIT_FAILURE;
 }
diff --git a/Applications/Utils/MeshEdit/AddLayer.cpp b/Applications/Utils/MeshEdit/AddLayer.cpp
index 7958e404e2de021985d5a768e0f156c75b8f388b..90e2d51efb788096293673bde20d6417dba8e092 100644
--- a/Applications/Utils/MeshEdit/AddLayer.cpp
+++ b/Applications/Utils/MeshEdit/AddLayer.cpp
@@ -11,13 +11,10 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <memory>
 
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/readMeshFromFile.h"
 #include "MeshLib/IO/writeMeshToFile.h"
@@ -78,15 +75,10 @@ int main(int argc, char* argv[])
     {
         ERR("It is not possible to set both options '--copy-material-ids' and "
             "'--set-material-id'.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     INFO("Reading mesh '{:s}' ... ", mesh_arg.getValue());
     auto subsfc_mesh = std::unique_ptr<MeshLib::Mesh>(
@@ -94,9 +86,6 @@ int main(int argc, char* argv[])
     if (!subsfc_mesh)
     {
         ERR("Error reading mesh '{:s}'.", mesh_arg.getValue());
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     INFO("done.");
@@ -114,9 +103,6 @@ int main(int argc, char* argv[])
     if (!result)
     {
         ERR("Failure while adding layer.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -124,8 +110,5 @@ int main(int argc, char* argv[])
     MeshLib::IO::writeMeshToFile(*result, mesh_out_arg.getValue());
     INFO("done.");
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/CreateBoundaryConditionsAlongPolylines.cpp b/Applications/Utils/MeshEdit/CreateBoundaryConditionsAlongPolylines.cpp
index 070bcedc532acb904b30592aae6d6e0d53d2aa0f..1debc20fb54ab5c2b8ae3d198d87fa3fc9c41bb1 100644
--- a/Applications/Utils/MeshEdit/CreateBoundaryConditionsAlongPolylines.cpp
+++ b/Applications/Utils/MeshEdit/CreateBoundaryConditionsAlongPolylines.cpp
@@ -12,10 +12,6 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <fstream>
 #include <map>
 #include <string>
@@ -24,6 +20,7 @@
 #include "Applications/FileIO/readGeometryFromFile.h"
 #include "Applications/FileIO/writeGeometryToFile.h"
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "GeoLib/GEOObjects.h"
 #include "GeoLib/Point.h"
 #include "InfoLib/GitInfo.h"
@@ -183,9 +180,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     // *** read mesh
     INFO("Reading mesh '{:s}' ... ", mesh_arg.getValue());
@@ -216,9 +211,6 @@ int main(int argc, char* argv[])
     {
         ERR("Could not get vector of polylines out of geometry '{:s}'.",
             geo_name);
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -251,9 +243,6 @@ int main(int argc, char* argv[])
     if (geo_names.empty())
     {
         ERR("Did not find mesh nodes along polylines.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -313,8 +302,5 @@ int main(int argc, char* argv[])
         BaseLib::dropFileExtension(output_base_fname.getValue()));
     writeBCsAndGeometry(geometry_sets, surface_name, base_fname,
                         bc_type.getValue(), gml_arg.getValue());
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/ExtractBoundary.cpp b/Applications/Utils/MeshEdit/ExtractBoundary.cpp
index 744929e8c64338a9b513091c2ca478988bc58cd4..f643dbd7d1052d5f249414fea16971560036902c 100644
--- a/Applications/Utils/MeshEdit/ExtractBoundary.cpp
+++ b/Applications/Utils/MeshEdit/ExtractBoundary.cpp
@@ -11,16 +11,13 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <algorithm>
 #include <memory>
 #include <string>
 #include <vector>
 
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "BaseLib/StringTools.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/Elements/Element.h"
@@ -63,18 +60,13 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::unique_ptr<MeshLib::Mesh const> mesh(MeshLib::IO::readMeshFromFile(
         mesh_in.getValue(), true /* compute_element_neighbors */));
 
     if (!mesh)
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -102,8 +94,5 @@ int main(int argc, char* argv[])
         use_ascii_arg.getValue() ? vtkXMLWriter::Ascii : vtkXMLWriter::Binary;
     MeshLib::IO::writeVtu(*surface_mesh, out_fname, data_mode);
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/ExtractMaterials.cpp b/Applications/Utils/MeshEdit/ExtractMaterials.cpp
index 9edf690aae1eb59b6c4f177e7bc46fa77442d728..050cf13656a9e731ff256df2939684162887676c 100644
--- a/Applications/Utils/MeshEdit/ExtractMaterials.cpp
+++ b/Applications/Utils/MeshEdit/ExtractMaterials.cpp
@@ -12,11 +12,8 @@
 // ThirdParty
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/VtkIO/VtuInterface.h"
 #include "MeshLib/IO/readMeshFromFile.h"
@@ -73,9 +70,7 @@ int main(int argc, char* argv[])
     cmd.add(input_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::string const input_name = input_arg.getValue();
     std::string const output_name = output_arg.getValue();
@@ -87,9 +82,6 @@ int main(int argc, char* argv[])
     if (mesh == nullptr)
     {
         ERR("Error reading input mesh. Aborting...");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -97,9 +89,6 @@ int main(int argc, char* argv[])
     if (mat_ids == nullptr)
     {
         ERR("No material IDs found in mesh. Aborting...");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -107,9 +96,6 @@ int main(int argc, char* argv[])
     if (id_range.first == id_range.second)
     {
         ERR("Mesh only contains one material, no extraction required.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     int min_id, max_id;
@@ -119,9 +105,6 @@ int main(int argc, char* argv[])
         if (min_id < *id_range.first || min_id > *id_range.second)
         {
             ERR("Specified material ID does not exist.");
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
         max_id = min_id;
@@ -160,8 +143,5 @@ int main(int argc, char* argv[])
     {
         ostream.close();
     }
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/ExtractSurface.cpp b/Applications/Utils/MeshEdit/ExtractSurface.cpp
index da959739497657ca995cd29790fadf317770955e..f8a95e304de731ab126bcf5e20caee46ff0fa74c 100644
--- a/Applications/Utils/MeshEdit/ExtractSurface.cpp
+++ b/Applications/Utils/MeshEdit/ExtractSurface.cpp
@@ -11,16 +11,13 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <algorithm>
 #include <memory>
 #include <string>
 #include <vector>
 
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "BaseLib/StringTools.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/VtkIO/VtuInterface.h"
@@ -75,9 +72,7 @@ int main(int argc, char* argv[])
     cmd.add(mesh_in);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::unique_ptr<MeshLib::Mesh const> mesh(MeshLib::IO::readMeshFromFile(
         mesh_in.getValue(), true /* compute_element_neighbors */));
@@ -85,18 +80,12 @@ int main(int argc, char* argv[])
     if (!mesh)
     {
         ERR("Error reading mesh file.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
     if (mesh->getDimension() != 3)
     {
         ERR("Surfaces can currently only be extracted from 3D meshes.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -122,8 +111,5 @@ int main(int argc, char* argv[])
         use_ascii_arg.getValue() ? vtkXMLWriter::Ascii : vtkXMLWriter::Binary;
     MeshLib::IO::writeVtu(*surface_mesh, out_fname, data_mode);
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/Layers2Grid.cpp b/Applications/Utils/MeshEdit/Layers2Grid.cpp
index 6108bee99fd7ff6ee37ae6c35819c27e6b01ab51..12525d1326c185e082bc66469e391a0520b6588f 100644
--- a/Applications/Utils/MeshEdit/Layers2Grid.cpp
+++ b/Applications/Utils/MeshEdit/Layers2Grid.cpp
@@ -15,11 +15,8 @@
 // ThirdParty
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include "BaseLib/IO/readStringListFromFile.h"
+#include "BaseLib/MPI.h"
 #include "GeoLib/AABB.h"
 #include "InfoLib/GitInfo.h"
 #include "MathLib/Point3d.h"
@@ -78,18 +75,13 @@ int main(int argc, char* argv[])
     cmd.add(input_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     if ((y_arg.isSet() && !z_arg.isSet()) ||
         ((!y_arg.isSet() && z_arg.isSet())))
     {
         ERR("For equilateral cubes, only x needs to be set. For unequal "
             "cuboids, all three edge lengths (x/y/z) need to be specified.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -104,9 +96,6 @@ int main(int argc, char* argv[])
     if (layer_names.size() < 2)
     {
         ERR("At least two layers are required to create a 3D Mesh");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -124,9 +113,6 @@ int main(int argc, char* argv[])
         if (mesh == nullptr)
         {
             ERR("Input layer '{:s}' not found. Aborting...", layer);
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
         layers.emplace_back(mesh);
@@ -141,15 +127,9 @@ int main(int argc, char* argv[])
     if (mesh == nullptr)
     {
         ERR("The VoxelGrid could not be created.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     MeshLib::IO::VtuInterface vtu(mesh.get());
     vtu.writeToFile(output_name);
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/MapGeometryToMeshSurface.cpp b/Applications/Utils/MeshEdit/MapGeometryToMeshSurface.cpp
index d32118f004c8c06014c5410b6f95b5adcf4b227c..106e8305d44f62f4e78d5f96aab0ad1c81c2a6da 100644
--- a/Applications/Utils/MeshEdit/MapGeometryToMeshSurface.cpp
+++ b/Applications/Utils/MeshEdit/MapGeometryToMeshSurface.cpp
@@ -11,14 +11,11 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <algorithm>
 #include <cstdlib>
 #include <vector>
 
+#include "BaseLib/MPI.h"
 #include "GeoLib/GEOObjects.h"
 #include "GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.h"
 #include "InfoLib/GitInfo.h"
@@ -60,9 +57,7 @@ int main(int argc, char* argv[])
     cmd.add(output_geometry_fname);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     // *** read geometry
     GeoLib::GEOObjects geometries;
@@ -75,9 +70,6 @@ int main(int argc, char* argv[])
         }
         else
         {
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
     }
@@ -105,8 +97,5 @@ int main(int argc, char* argv[])
         BaseLib::IO::writeStringToFile(xml_io.writeToString(),
                                        output_geometry_fname.getValue());
     }
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/MeshMapping.cpp b/Applications/Utils/MeshEdit/MeshMapping.cpp
index 6b939c1821d47d9e96f3a97599b963b5ff1d35cc..3d3f3139fda3eec08eefe48f29609afa385644ba 100644
--- a/Applications/Utils/MeshEdit/MeshMapping.cpp
+++ b/Applications/Utils/MeshEdit/MeshMapping.cpp
@@ -15,11 +15,8 @@
 #include <memory>
 #include <string>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "GeoLib/AABB.h"
 #include "GeoLib/IO/AsciiRasterInterface.h"
 #include "GeoLib/Raster.h"
@@ -54,9 +51,7 @@ double getClosestPointElevation(MeshLib::Node const& p,
 
 int main(int argc, char* argv[])
 {
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
     TCLAP::CmdLine cmd(
         "Changes the elevation of 2D mesh nodes based on either raster data or "
         "another 2D mesh. In addition, a low pass filter can be applied to "
@@ -107,9 +102,6 @@ int main(int argc, char* argv[])
     if (mesh == nullptr)
     {
         ERR("Error reading mesh file.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -118,18 +110,12 @@ int main(int argc, char* argv[])
     {
         ERR("Nothing to do. Please choose mapping based on a raster or mesh "
             "file, or to a static value.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
     if (map_raster_arg.isSet() && map_mesh_arg.isSet())
     {
         ERR("Please select mapping based on *either* a mesh or a raster file.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -148,9 +134,6 @@ int main(int argc, char* argv[])
         if (!file_stream.good())
         {
             ERR("Opening raster file {} failed.", raster_path);
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
         file_stream.close();
@@ -168,9 +151,6 @@ int main(int argc, char* argv[])
         if (ground_truth == nullptr)
         {
             ERR("Error reading mesh file.");
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
 
@@ -243,15 +223,9 @@ int main(int argc, char* argv[])
 
     if (MeshLib::IO::writeMeshToFile(*mesh, output_arg.getValue()) != 0)
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
     INFO("Result successfully written.");
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/MoveMesh.cpp b/Applications/Utils/MeshEdit/MoveMesh.cpp
index a755a792143525429c83aad7152557609c65b3bf..65cefe2c3bb0ca7d1df876d104fc4af237a8aeff 100644
--- a/Applications/Utils/MeshEdit/MoveMesh.cpp
+++ b/Applications/Utils/MeshEdit/MoveMesh.cpp
@@ -11,11 +11,8 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "BaseLib/StringTools.h"
 #include "GeoLib/AABB.h"
 #include "InfoLib/GitInfo.h"
@@ -64,9 +61,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
     std::string fname(mesh_arg.getValue());
 
     std::unique_ptr<MeshLib::Mesh> mesh(MeshLib::IO::readMeshFromFile(fname));
@@ -74,9 +69,6 @@ int main(int argc, char* argv[])
     if (!mesh)
     {
         ERR("Could not read mesh from file '{:s}'.", fname);
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -107,8 +99,5 @@ int main(int argc, char* argv[])
 
     MeshLib::IO::writeMeshToFile(*mesh, out_fname);
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/NodeReordering.cpp b/Applications/Utils/MeshEdit/NodeReordering.cpp
index cc83457c3a48ef1528d41a41ff9ef977e780cb73..f98f006397521cd019316c530621a0b842c26afb 100644
--- a/Applications/Utils/MeshEdit/NodeReordering.cpp
+++ b/Applications/Utils/MeshEdit/NodeReordering.cpp
@@ -11,16 +11,13 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <algorithm>
 #include <array>
 #include <memory>
 #include <vector>
 
 #include "BaseLib/Algorithm.h"
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/Elements/Element.h"
 #include "MeshLib/IO/readMeshFromFile.h"
@@ -191,18 +188,13 @@ int main(int argc, char* argv[])
     cmd.add(input_mesh_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::unique_ptr<MeshLib::Mesh> mesh(
         MeshLib::IO::readMeshFromFile(input_mesh_arg.getValue()));
 
     if (!mesh)
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -228,8 +220,5 @@ int main(int argc, char* argv[])
 
     INFO("VTU file written.");
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/PVTU2VTU/PVTU2VTU.cpp b/Applications/Utils/MeshEdit/PVTU2VTU/PVTU2VTU.cpp
index 3e166e0734b297801bfa74e2ec0047bf21b02a93..753004bf2bcd49485eeb033a6948002406e485cc 100644
--- a/Applications/Utils/MeshEdit/PVTU2VTU/PVTU2VTU.cpp
+++ b/Applications/Utils/MeshEdit/PVTU2VTU/PVTU2VTU.cpp
@@ -25,6 +25,7 @@
 #include <vector>
 
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "BaseLib/RunTime.h"
 #include "GeoLib/AABB.h"
 #include "GeoLib/OctTree.h"
@@ -229,9 +230,8 @@ getMergedNodesVector(std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes)
     ranges::transform(meshes, std::back_inserter(number_of_nodes_per_partition),
                       [](auto const& mesh)
                       { return mesh->getNumberOfNodes(); });
-    std::vector<std::size_t> offsets(meshes.size() + 1, 0);
-    std::partial_sum(number_of_nodes_per_partition.begin(),
-                     number_of_nodes_per_partition.end(), offsets.begin() + 1);
+    std::vector<std::size_t> const offsets =
+        BaseLib::sizesToOffsets(number_of_nodes_per_partition);
 
     std::vector<MeshLib::Node*> all_nodes;
     all_nodes.reserve(offsets.back());
@@ -387,6 +387,8 @@ int main(int argc, char* argv[])
     cmd.add(input_arg);
     cmd.parse(argc, argv);
 
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
+
     if (BaseLib::getFileExtension(input_arg.getValue()) != ".pvtu")
     {
         OGS_FATAL("The extension of input file name {:s} is not \"pvtu\"",
diff --git a/Applications/Utils/MeshEdit/RemoveGhostData.cpp b/Applications/Utils/MeshEdit/RemoveGhostData.cpp
index 1824b8661c7e0d14d3ad1092bd15d75d210302db..1384ac736b2a08ad1a20bf397c06a00e52b8b978 100644
--- a/Applications/Utils/MeshEdit/RemoveGhostData.cpp
+++ b/Applications/Utils/MeshEdit/RemoveGhostData.cpp
@@ -15,11 +15,8 @@
 #include <vtkXMLPUnstructuredGridReader.h>
 #include <vtkXMLUnstructuredGridWriter.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include "BaseLib/Logging.h"
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/writeMeshToFile.h"
 
@@ -51,9 +48,7 @@ int main(int argc, char* argv[])
     cmd.add(input_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     vtkSmartPointer<vtkXMLPUnstructuredGridReader> reader =
         vtkSmartPointer<vtkXMLPUnstructuredGridReader>::New();
@@ -73,8 +68,5 @@ int main(int argc, char* argv[])
     writer->SetFileName(output_arg.getValue().c_str());
     writer->Write();
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/ReorderMesh.cpp b/Applications/Utils/MeshEdit/ReorderMesh.cpp
index 220571ce22a940bd520706af5f09d3698836c9ff..5de0c9879c880ae1ac40eb9be2f1ce603e7233fe 100644
--- a/Applications/Utils/MeshEdit/ReorderMesh.cpp
+++ b/Applications/Utils/MeshEdit/ReorderMesh.cpp
@@ -12,6 +12,7 @@
 #include <memory>
 #include <string>
 
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/Elements/Element.h"
 #include "MeshLib/IO/readMeshFromFile.h"
@@ -132,6 +133,8 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
+
     const std::string filename(mesh_in_arg.getValue());
 
     // read the mesh file
diff --git a/Applications/Utils/MeshEdit/ResetPropertiesInPolygonalRegion.cpp b/Applications/Utils/MeshEdit/ResetPropertiesInPolygonalRegion.cpp
index 76631d3ea973f19b88d0cf9de818b6bed11a343f..ada2f29b6dd0f78adf46abfac93e1dc0d0fc7ed0 100644
--- a/Applications/Utils/MeshEdit/ResetPropertiesInPolygonalRegion.cpp
+++ b/Applications/Utils/MeshEdit/ResetPropertiesInPolygonalRegion.cpp
@@ -16,11 +16,8 @@
 #include <cstdlib>
 #include <vector>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include "Applications/FileIO/readGeometryFromFile.h"
+#include "BaseLib/MPI.h"
 #include "GeoLib/GEOObjects.h"
 #include "GeoLib/Polygon.h"
 #include "InfoLib/GitInfo.h"
@@ -99,9 +96,7 @@ int main(int argc, char* argv[])
     cmd.add(gmsh_path_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     // *** read geometry
     GeoLib::GEOObjects geometries;
@@ -117,9 +112,6 @@ int main(int argc, char* argv[])
     {
         ERR("Could not get vector of polylines out of geometry '{:s}'.",
             geo_name);
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -129,9 +121,6 @@ int main(int argc, char* argv[])
     if (!ply)
     {
         ERR("Polyline '{:s}' not found.", polygon_name_arg.getValue());
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -140,9 +129,6 @@ int main(int argc, char* argv[])
     {
         ERR("Polyline '{:s}' is not closed, i.e. does not describe a region.",
             polygon_name_arg.getValue());
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -154,9 +140,6 @@ int main(int argc, char* argv[])
     if (!mesh)
     {
         // error message written already by readMeshFromFile
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     std::string const& property_name(property_name_arg.getValue());
@@ -186,8 +169,5 @@ int main(int argc, char* argv[])
 
     MeshLib::IO::writeMeshToFile(*mesh, mesh_out.getValue());
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/Vtu2Grid.cpp b/Applications/Utils/MeshEdit/Vtu2Grid.cpp
index 7429e1d9f72e7f45db1944012b1f0dc217891ec1..55a8bc6f5821b1e41b2c64c20b7deba7a3b17ca3 100644
--- a/Applications/Utils/MeshEdit/Vtu2Grid.cpp
+++ b/Applications/Utils/MeshEdit/Vtu2Grid.cpp
@@ -10,6 +10,7 @@
 #include <tclap/CmdLine.h>
 #include <vtkXMLUnstructuredGridReader.h>
 
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/writeMeshToFile.h"
 #include "MeshLib/Mesh.h"
@@ -18,9 +19,6 @@
 #include "MeshToolsLib/MeshGenerators/MeshGenerator.h"
 #include "MeshToolsLib/MeshGenerators/VoxelGridFromMesh.h"
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
 int main(int argc, char* argv[])
 {
     TCLAP::CmdLine cmd(
@@ -64,18 +62,13 @@ int main(int argc, char* argv[])
     cmd.add(input_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     if ((y_arg.isSet() && !z_arg.isSet()) ||
         ((!y_arg.isSet() && z_arg.isSet())))
     {
         ERR("For equilateral cubes, only x needs to be set. For unequal "
             "cuboids, all three edge lengths (x/y/z) need to be specified.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return -1;
     }
     using namespace MeshToolsLib::MeshGenerator;
@@ -88,9 +81,6 @@ int main(int argc, char* argv[])
     {
         ERR("A cellsize ({},{},{}) is not allowed to be <= 0", x_size, y_size,
             z_size);
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return -1;
     }
 
@@ -113,9 +103,6 @@ int main(int argc, char* argv[])
     {
         ERR("The range ({},{},{}) is not allowed to be < 0", ranges[0],
             ranges[1], ranges[2]);
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return -1;
     }
     std::array<std::size_t, 3> const dims =
@@ -134,9 +121,6 @@ int main(int argc, char* argv[])
 
     if (!VoxelGridFromMesh::removeUnusedGridCells(mesh, grid))
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -144,14 +128,8 @@ int main(int argc, char* argv[])
 
     if (MeshLib::IO::writeMeshToFile(*grid, output_arg.getValue()) != 0)
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/appendLinesAlongPolyline.cpp b/Applications/Utils/MeshEdit/appendLinesAlongPolyline.cpp
index e7ab944bd296bf74ca5e078b603b0ced23f90f28..69e3baa152283e077226a0e0912bb58b7b7fe714 100644
--- a/Applications/Utils/MeshEdit/appendLinesAlongPolyline.cpp
+++ b/Applications/Utils/MeshEdit/appendLinesAlongPolyline.cpp
@@ -11,12 +11,9 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include "Applications/FileIO/readGeometryFromFile.h"
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "GeoLib/GEOObjects.h"
 #include "GeoLib/PolylineVec.h"
 #include "InfoLib/GitInfo.h"
@@ -58,9 +55,7 @@ int main(int argc, char* argv[])
     // parse arguments
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     // read GEO objects
     GeoLib::GEOObjects geo_objs;
@@ -71,9 +66,6 @@ int main(int argc, char* argv[])
     if (geo_names.empty())
     {
         ERR("No geometries found.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     const GeoLib::PolylineVec* ply_vec(
@@ -81,9 +73,6 @@ int main(int argc, char* argv[])
     if (!ply_vec)
     {
         ERR("Could not find polylines in geometry '{:s}'.", geo_names.front());
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -93,9 +82,6 @@ int main(int argc, char* argv[])
     if (!mesh)
     {
         ERR("Mesh file '{:s}' not found", mesh_in.getValue());
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     INFO("Mesh read: {:d} nodes, {:d} elements.", mesh->getNumberOfNodes(),
@@ -109,8 +95,5 @@ int main(int argc, char* argv[])
 
     MeshLib::IO::writeMeshToFile(*new_mesh, mesh_out.getValue());
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/checkMesh.cpp b/Applications/Utils/MeshEdit/checkMesh.cpp
index 4795937702cf6b75c7a8bc3db00bce38059f22a6..2d31816fd60fb5270ae2dd56f01827c5cb4d66d1 100644
--- a/Applications/Utils/MeshEdit/checkMesh.cpp
+++ b/Applications/Utils/MeshEdit/checkMesh.cpp
@@ -7,18 +7,14 @@
  *              http://www.opengeosys.org/project/license
  */
 
-#include <tclap/CmdLine.h>
-
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <spdlog/fmt/bundled/ranges.h>
+#include <tclap/CmdLine.h>
 
 #include <array>
 #include <string>
 
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "BaseLib/MemWatch.h"
 #include "BaseLib/RunTime.h"
 #include "BaseLib/StringTools.h"
@@ -52,9 +48,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     // read the mesh file
     BaseLib::MemWatch mem_watch;
@@ -65,9 +59,6 @@ int main(int argc, char* argv[])
         mesh_arg.getValue(), true /* compute_element_neighbors */));
     if (!mesh)
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -106,8 +97,5 @@ int main(int argc, char* argv[])
         // Remark: MeshValidation can modify the original mesh
         MeshToolsLib::MeshInformation::writeMeshValidationResults(*mesh);
     }
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/convertToLinearMesh.cpp b/Applications/Utils/MeshEdit/convertToLinearMesh.cpp
index 2eeea03f024b7bb8b6cde9bb0179ebe0203291c3..61569eaaba66bfe1250318530b3ac76511ed0c0f 100644
--- a/Applications/Utils/MeshEdit/convertToLinearMesh.cpp
+++ b/Applications/Utils/MeshEdit/convertToLinearMesh.cpp
@@ -11,13 +11,10 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <memory>
 #include <string>
 
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/readMeshFromFile.h"
 #include "MeshLib/IO/writeMeshToFile.h"
@@ -42,25 +39,17 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::unique_ptr<MeshLib::Mesh> mesh(
         MeshLib::IO::readMeshFromFile(input_arg.getValue()));
     if (!mesh)
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     if (!mesh->hasNonlinearElement())
     {
         ERR("The input mesh is linear. Exit.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -71,8 +60,5 @@ int main(int argc, char* argv[])
     INFO("Save the new mesh into a file");
     MeshLib::IO::writeMeshToFile(*new_mesh, output_arg.getValue());
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/createLayeredMeshFromRasters.cpp b/Applications/Utils/MeshEdit/createLayeredMeshFromRasters.cpp
index 91ab223b4bc1145f9cb1d90d8f828b7a48a9d860..56805e54f9470572249e6783ae31f8a7c213d731 100644
--- a/Applications/Utils/MeshEdit/createLayeredMeshFromRasters.cpp
+++ b/Applications/Utils/MeshEdit/createLayeredMeshFromRasters.cpp
@@ -18,12 +18,9 @@
 #include <string>
 #include <vector>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include "BaseLib/FileTools.h"
 #include "BaseLib/IO/readStringListFromFile.h"
+#include "BaseLib/MPI.h"
 #include "GeoLib/IO/AsciiRasterInterface.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/VtkIO/VtuInterface.h"
@@ -79,9 +76,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     if (min_thickness_arg.isSet())
     {
@@ -89,9 +84,6 @@ int main(int argc, char* argv[])
         if (min_thickness < 0)
         {
             ERR("Minimum layer thickness must be non-negative value.");
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
     }
@@ -102,17 +94,11 @@ int main(int argc, char* argv[])
     if (!sfc_mesh)
     {
         ERR("Error reading mesh '{:s}'.", mesh_arg.getValue());
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     if (sfc_mesh->getDimension() != 2)
     {
         ERR("Input mesh must be a 2D mesh.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     INFO("done.");
@@ -122,9 +108,6 @@ int main(int argc, char* argv[])
     if (raster_paths.size() < 2)
     {
         ERR("At least two raster files needed to create 3D mesh.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     std::reverse(raster_paths.begin(), raster_paths.end());
@@ -134,18 +117,12 @@ int main(int argc, char* argv[])
     {
         if (!mapper.createLayers(*sfc_mesh, *rasters, min_thickness))
         {
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
     }
     else
     {
         ERR("Reading raster files.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -160,9 +137,6 @@ int main(int argc, char* argv[])
     if (result_mesh == nullptr)
     {
         ERR("Mapper returned empty result for 'SubsurfaceMesh'.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -172,8 +146,5 @@ int main(int argc, char* argv[])
     MeshLib::IO::writeVtu(*result_mesh, output_name, data_mode);
     INFO("done.");
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/createQuadraticMesh.cpp b/Applications/Utils/MeshEdit/createQuadraticMesh.cpp
index 735c37537c458dd400be9c7355755e55e5fbf9b5..6758cf37146db4e312c87ef78190a21000d745d2 100644
--- a/Applications/Utils/MeshEdit/createQuadraticMesh.cpp
+++ b/Applications/Utils/MeshEdit/createQuadraticMesh.cpp
@@ -9,13 +9,10 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <memory>
 #include <string>
 
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/readMeshFromFile.h"
 #include "MeshLib/IO/writeMeshToFile.h"
@@ -44,17 +41,12 @@ int main(int argc, char* argv[])
     cmd.add(add_centre_node_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::unique_ptr<MeshLib::Mesh> mesh(
         MeshLib::IO::readMeshFromFile(input_arg.getValue()));
     if (!mesh)
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -65,8 +57,5 @@ int main(int argc, char* argv[])
     INFO("Save the new mesh into a file");
     MeshLib::IO::writeMeshToFile(*new_mesh, output_arg.getValue());
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/createTetgenSmeshFromRasters.cpp b/Applications/Utils/MeshEdit/createTetgenSmeshFromRasters.cpp
index 838eb9bde55d8c5fc218a334ea36f0550136e614..4b2a881ff699154c1c2e7cf62e54c5ae57549c25 100644
--- a/Applications/Utils/MeshEdit/createTetgenSmeshFromRasters.cpp
+++ b/Applications/Utils/MeshEdit/createTetgenSmeshFromRasters.cpp
@@ -18,13 +18,10 @@
 #include <string>
 #include <vector>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include "Applications/FileIO/TetGenInterface.h"
 #include "BaseLib/FileTools.h"
 #include "BaseLib/IO/readStringListFromFile.h"
+#include "BaseLib/MPI.h"
 #include "GeoLib/IO/AsciiRasterInterface.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/VtkIO/VtuInterface.h"
@@ -85,9 +82,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     if (min_thickness_arg.isSet())
     {
@@ -95,9 +90,6 @@ int main(int argc, char* argv[])
         if (min_thickness < 0)
         {
             ERR("Minimum layer thickness must be non-negative value.");
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
     }
@@ -108,17 +100,11 @@ int main(int argc, char* argv[])
     if (!sfc_mesh)
     {
         ERR("Error reading mesh '{:s}'.", mesh_arg.getValue());
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     if (sfc_mesh->getDimension() != 2)
     {
         ERR("Input mesh must be a 2D mesh.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     INFO("done.");
@@ -128,9 +114,6 @@ int main(int argc, char* argv[])
     if (raster_paths.size() < 2)
     {
         ERR("At least two raster files needed to create 3D mesh.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     std::reverse(raster_paths.begin(), raster_paths.end());
@@ -148,18 +131,12 @@ int main(int argc, char* argv[])
         if (!lv.createLayers(*sfc_mesh, *rasters, min_thickness))
         {
             ERR("Creating the layers was erroneous.");
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
     }
     else
     {
         ERR("The raster files are not accessible.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     std::unique_ptr<MeshLib::Mesh> tg_mesh =
@@ -174,14 +151,8 @@ int main(int argc, char* argv[])
     else
     {
         ERR("The tetgen-smesh could not be created.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/editMaterialID.cpp b/Applications/Utils/MeshEdit/editMaterialID.cpp
index 179997796907f5644b346a5f4ae24b3c69e9ff75..61d0892de4625b0e7a1319510e2716f47d1705e0 100644
--- a/Applications/Utils/MeshEdit/editMaterialID.cpp
+++ b/Applications/Utils/MeshEdit/editMaterialID.cpp
@@ -9,12 +9,9 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <memory>
 
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/Elements/Element.h"
 #include "MeshLib/IO/readMeshFromFile.h"
@@ -68,24 +65,16 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     if (!replaceArg.isSet() && !condenseArg.isSet() && !specifyArg.isSet())
     {
         INFO("Please select editing mode: -r or -c or -s");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return 0;
     }
     if (replaceArg.isSet() && condenseArg.isSet())
     {
         INFO("Please select only one editing mode: -r or -c or -s");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return 0;
     }
     if (replaceArg.isSet())
@@ -95,9 +84,6 @@ int main(int argc, char* argv[])
             INFO(
                 "current and new material IDs must be provided for "
                 "replacement");
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return 0;
         }
     }
@@ -108,9 +94,6 @@ int main(int argc, char* argv[])
             INFO(
                 "element type and new material IDs must be provided to specify "
                 "elements");
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return 0;
         }
     }
@@ -170,8 +153,5 @@ int main(int argc, char* argv[])
     }
     MeshLib::IO::writeMeshToFile(*mesh, mesh_out.getValue());
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/ipDataToPointCloud.cpp b/Applications/Utils/MeshEdit/ipDataToPointCloud.cpp
index 0004f85b724c6b3cadb827bea47016696e9703e8..d1fda1f1f9d93b4c06492c3ccfef13a36fe45f86 100644
--- a/Applications/Utils/MeshEdit/ipDataToPointCloud.cpp
+++ b/Applications/Utils/MeshEdit/ipDataToPointCloud.cpp
@@ -11,6 +11,7 @@
 
 #include <unordered_map>
 
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/readMeshFromFile.h"
 #include "MeshLib/IO/writeMeshToFile.h"
@@ -251,6 +252,8 @@ int main(int argc, char** argv)
 
     cmd.parse(argc, argv);
 
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
+
     std::unique_ptr<MeshLib::Mesh const> mesh_in(
         MeshLib::IO::readMeshFromFile(arg_in_file.getValue()));
 
diff --git a/Applications/Utils/MeshEdit/queryMesh.cpp b/Applications/Utils/MeshEdit/queryMesh.cpp
index fd713e8fe74b52b847d04e6e968859c0cd3adb8d..c4a3cb9d3e04517fe983acbd366d98f68aac0441 100644
--- a/Applications/Utils/MeshEdit/queryMesh.cpp
+++ b/Applications/Utils/MeshEdit/queryMesh.cpp
@@ -9,16 +9,13 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <array>
 #include <memory>
 #include <sstream>
 #include <string>
 
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "BaseLib/StringTools.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/Elements/Element.h"
@@ -52,9 +49,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     const std::string filename(mesh_arg.getValue());
 
@@ -64,9 +59,6 @@ int main(int argc, char* argv[])
             filename, true /* compute_element_neighbors */));
     if (!mesh)
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -155,8 +147,5 @@ int main(int argc, char* argv[])
         out << std::endl;
         INFO("{:s}", out.str());
     }
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/removeMeshElements.cpp b/Applications/Utils/MeshEdit/removeMeshElements.cpp
index 456fe2635f19fa2855c9a339e242044514c81591..bd20620b226dfff18a94a8fe2bcff5c0f1f92dec 100644
--- a/Applications/Utils/MeshEdit/removeMeshElements.cpp
+++ b/Applications/Utils/MeshEdit/removeMeshElements.cpp
@@ -13,12 +13,9 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <memory>
 
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/Elements/Element.h"
 #include "MeshLib/IO/readMeshFromFile.h"
@@ -182,17 +179,12 @@ int main(int argc, char* argv[])
     cmd.add(mesh_in);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::unique_ptr<MeshLib::Mesh const> mesh(
         MeshLib::IO::readMeshFromFile(mesh_in.getValue()));
     if (mesh == nullptr)
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -229,9 +221,6 @@ int main(int argc, char* argv[])
             !property_name_arg.isSet())
         {
             ERR("Specify a property name for the value/range selected.");
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
 
@@ -241,9 +230,6 @@ int main(int argc, char* argv[])
         {
             ERR("Specify a value or range ('-min-value' and '-max_value') for "
                 "the property selected.");
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
 
@@ -264,9 +250,6 @@ int main(int argc, char* argv[])
             {
                 ERR("Specify if the inside or the outside of the selected "
                     "range should be removed.");
-#ifdef USE_PETSC
-                MPI_Finalize();
-#endif
                 return EXIT_FAILURE;
             }
 
@@ -299,9 +282,6 @@ int main(int argc, char* argv[])
         }
         if (aabb_error)
         {
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
 
@@ -324,17 +304,11 @@ int main(int argc, char* argv[])
 
     if (new_mesh == nullptr)
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
     // write into a file
     MeshLib::IO::writeMeshToFile(*new_mesh, mesh_out.getValue());
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/reviseMesh.cpp b/Applications/Utils/MeshEdit/reviseMesh.cpp
index ff9c0d3509ac84a03820ab764f27a845a6e23f3b..ca06a00b084831892f840d6cf2d5e34a29be6b59 100644
--- a/Applications/Utils/MeshEdit/reviseMesh.cpp
+++ b/Applications/Utils/MeshEdit/reviseMesh.cpp
@@ -9,15 +9,12 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <array>
 #include <memory>
 #include <string>
 
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "BaseLib/StringTools.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/Elements/Element.h"
@@ -59,18 +56,13 @@ int main(int argc, char* argv[])
     cmd.add(input_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     // read a mesh file
     std::unique_ptr<MeshLib::Mesh> org_mesh(
         MeshLib::IO::readMeshFromFile(input_arg.getValue()));
     if (!org_mesh)
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     INFO("Mesh read: {:d} nodes, {:d} elements.", org_mesh->getNumberOfNodes(),
@@ -92,8 +84,5 @@ int main(int argc, char* argv[])
         MeshLib::IO::writeMeshToFile(*new_mesh, output_arg.getValue());
     }
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshEdit/swapNodeCoordinateAxes.cpp b/Applications/Utils/MeshEdit/swapNodeCoordinateAxes.cpp
index dcc9558a7f1269bd8740919ad3916a803b837d48..0eaadc7013e403b19020184ca3bb47ed82f53d3e 100644
--- a/Applications/Utils/MeshEdit/swapNodeCoordinateAxes.cpp
+++ b/Applications/Utils/MeshEdit/swapNodeCoordinateAxes.cpp
@@ -9,15 +9,12 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <array>
 #include <memory>
 #include <string>
 
 #include "BaseLib/Logging.h"
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/readMeshFromFile.h"
 #include "MeshLib/IO/writeMeshToFile.h"
@@ -114,17 +111,12 @@ int main(int argc, char* argv[])
     cmd.add(new_order_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     const std::string str_order = new_order_arg.getValue();
     std::array<int, 3> new_order = {{}};
     if (!parseNewOrder(str_order, new_order))
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -132,9 +124,6 @@ int main(int argc, char* argv[])
         MeshLib::IO::readMeshFromFile(input_arg.getValue()));
     if (!mesh)
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -152,8 +141,5 @@ int main(int argc, char* argv[])
     INFO("Save the new mesh into a file");
     MeshLib::IO::writeMeshToFile(*mesh, output_arg.getValue());
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshGeoTools/AssignRasterDataToMesh.cpp b/Applications/Utils/MeshGeoTools/AssignRasterDataToMesh.cpp
index ceada5695ec3773702a1fdbe69c68812ce0a7c68..338354d0d6eb7bc62e438e33425fa29c228c6944 100644
--- a/Applications/Utils/MeshGeoTools/AssignRasterDataToMesh.cpp
+++ b/Applications/Utils/MeshGeoTools/AssignRasterDataToMesh.cpp
@@ -13,10 +13,7 @@
 // ThirdParty
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
+#include "BaseLib/MPI.h"
 #include "GeoLib/IO/AsciiRasterInterface.h"
 #include "GeoLib/Raster.h"
 #include "InfoLib/GitInfo.h"
@@ -74,9 +71,7 @@ int main(int argc, char* argv[])
     cmd.add(input_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     bool const create_cell_array(set_cells_arg.isSet());
     bool const create_node_array =
@@ -91,9 +86,6 @@ int main(int argc, char* argv[])
     if (mesh->getDimension() > 2)
     {
         ERR("Method can not be applied to 3D meshes.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -108,9 +100,6 @@ int main(int argc, char* argv[])
         if (!assigned)
         {
             ERR("Error assigning raster data to scalar node array");
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
         INFO("Created node array {:s}", array_name_arg.getValue());
@@ -124,9 +113,6 @@ int main(int argc, char* argv[])
         if (!assigned)
         {
             ERR("Error assigning raster data to scalar cell array");
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
         INFO("Created cell array {:s}", array_name_arg.getValue());
@@ -134,8 +120,5 @@ int main(int argc, char* argv[])
 
     MeshLib::IO::VtuInterface vtu(mesh.get());
     vtu.writeToFile(output_name);
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshGeoTools/IntegrateBoreholesIntoMesh.cpp b/Applications/Utils/MeshGeoTools/IntegrateBoreholesIntoMesh.cpp
index 6f0f5b701cd13be145048b1ac4f9b3b3e8b1ecd1..9c88e9215889506dd96bf3d62253eeac9552064c 100644
--- a/Applications/Utils/MeshGeoTools/IntegrateBoreholesIntoMesh.cpp
+++ b/Applications/Utils/MeshGeoTools/IntegrateBoreholesIntoMesh.cpp
@@ -16,10 +16,7 @@
 // ThirdParty
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
+#include "BaseLib/MPI.h"
 #include "GeoLib/GEOObjects.h"
 #include "GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.h"
 #include "InfoLib/GitInfo.h"
@@ -122,9 +119,7 @@ int main(int argc, char* argv[])
     cmd.add(input_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::pair<int, int> mat_limits(0, std::numeric_limits<int>::max());
     std::pair<double, double> elevation_limits(
@@ -134,9 +129,6 @@ int main(int argc, char* argv[])
     {
         ERR("If minimum MaterialID is set, maximum ID must be set, too (and "
             "vice versa).");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     if (min_id_arg.isSet() && max_id_arg.isSet())
@@ -151,18 +143,12 @@ int main(int argc, char* argv[])
     if (min_id_arg.isSet() && (mat_limits.first < 0 || mat_limits.second < 0))
     {
         ERR("Specified MaterialIDs must have non-negative values.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     if (min_elevation_arg.isSet() != max_elevation_arg.isSet())
     {
         ERR("If minimum elevation is set, maximum elevation must be set, too "
             "(and vice versa).");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     if (min_elevation_arg.isSet() && max_elevation_arg.isSet())
@@ -184,9 +170,6 @@ int main(int argc, char* argv[])
     if (!xml_io.readFile(geo_name))
     {
         ERR("Failed to read geometry file `{:s}'.", geo_name);
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     std::vector<GeoLib::Point*> const& points =
@@ -197,17 +180,11 @@ int main(int argc, char* argv[])
     if (mesh == nullptr)
     {
         ERR("Failed to read input mesh file `{:s}'.", mesh_name);
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     if (mesh->getDimension() != 3)
     {
         ERR("Method can only be applied to 3D meshes.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -216,9 +193,6 @@ int main(int argc, char* argv[])
     if (mat_ids == nullptr)
     {
         ERR("Mesh is required to have MaterialIDs");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -256,8 +230,5 @@ int main(int argc, char* argv[])
                                true /* compute_element_neighbors */, props);
     MeshLib::IO::VtuInterface vtu(&result);
     vtu.writeToFile(output_name);
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshGeoTools/Raster2Mesh.cpp b/Applications/Utils/MeshGeoTools/Raster2Mesh.cpp
index c4905205da3b5f9de2be99d36479b576b7899944..423ecc4d961e5818d6e6126e16023bd8407362e0 100644
--- a/Applications/Utils/MeshGeoTools/Raster2Mesh.cpp
+++ b/Applications/Utils/MeshGeoTools/Raster2Mesh.cpp
@@ -13,10 +13,7 @@
 // ThirdParty
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
+#include "BaseLib/MPI.h"
 #include "GeoLib/IO/AsciiRasterInterface.h"
 #include "GeoLib/Raster.h"
 #include "InfoLib/GitInfo.h"
@@ -72,9 +69,7 @@ int main(int argc, char* argv[])
     cmd.add(input_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::string const input_name = input_arg.getValue().c_str();
     std::string const output_name = output_arg.getValue().c_str();
@@ -107,16 +102,10 @@ int main(int argc, char* argv[])
     if (mesh == nullptr)
     {
         ERR("Conversion failed.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
     MeshLib::IO::VtuInterface vtu(mesh.get());
     vtu.writeToFile(output_name);
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshGeoTools/VerticalSliceFromLayers.cpp b/Applications/Utils/MeshGeoTools/VerticalSliceFromLayers.cpp
index 671cdaeeae9148c3db68de145dd1817e58845774..f9193a588cdfb36418a70beaecefd21fa5bb8183 100644
--- a/Applications/Utils/MeshGeoTools/VerticalSliceFromLayers.cpp
+++ b/Applications/Utils/MeshGeoTools/VerticalSliceFromLayers.cpp
@@ -16,16 +16,13 @@
 // ThirdParty
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <QCoreApplication>
 
 #include "Applications/FileIO/Gmsh/GMSHInterface.h"
 #include "Applications/FileIO/Gmsh/GmshReader.h"
 #include "BaseLib/FileTools.h"
 #include "BaseLib/IO/readStringListFromFile.h"
+#include "BaseLib/MPI.h"
 #include "GeoLib/AABB.h"
 #include "GeoLib/AnalyticalGeometry.h"
 #include "GeoLib/GEOObjects.h"
@@ -418,9 +415,7 @@ int main(int argc, char* argv[])
     cmd.add(input_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::string const input_name = input_arg.getValue();
     std::string const output_name = output_arg.getValue();
@@ -439,9 +434,6 @@ int main(int argc, char* argv[])
     if (layer_names.size() < 2)
     {
         ERR("At least two layers are required to extract a slice.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -453,9 +445,6 @@ int main(int argc, char* argv[])
     {
         ERR("Less than two geometries could be created from layers. Aborting "
             "extraction...");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -475,9 +464,6 @@ int main(int argc, char* argv[])
     if (mesh == nullptr)
     {
         ERR("Error generating mesh... (GMSH was unable to output mesh)");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     rotateMesh(*mesh, rot_mat, z_shift);
@@ -485,9 +471,6 @@ int main(int argc, char* argv[])
     if (new_mesh == nullptr)
     {
         ERR("Error generating mesh... (GMSH created line mesh)");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -506,8 +489,5 @@ int main(int argc, char* argv[])
     }
     MeshLib::IO::VtuInterface vtu(revised_mesh.get());
     vtu.writeToFile(output_name + ".vtu");
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshGeoTools/computeSurfaceNodeIDsInPolygonalRegion.cpp b/Applications/Utils/MeshGeoTools/computeSurfaceNodeIDsInPolygonalRegion.cpp
index 5329c0b9ad9b2bc14bbca3684241043d1ef4ee36..4ca7478f6d46974686f59d28dddc2a979503c829 100644
--- a/Applications/Utils/MeshGeoTools/computeSurfaceNodeIDsInPolygonalRegion.cpp
+++ b/Applications/Utils/MeshGeoTools/computeSurfaceNodeIDsInPolygonalRegion.cpp
@@ -12,10 +12,6 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <algorithm>
 #include <fstream>
 #include <memory>
@@ -25,6 +21,7 @@
 #include "Applications/FileIO/readGeometryFromFile.h"
 #include "BaseLib/Error.h"
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "GeoLib/GEOObjects.h"
 #include "GeoLib/Polygon.h"
 #include "InfoLib/GitInfo.h"
@@ -101,9 +98,7 @@ int main(int argc, char* argv[])
     cmd.add(gmsh_path_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::unique_ptr<MeshLib::Mesh const> mesh(
         MeshLib::IO::readMeshFromFile(mesh_in.getValue()));
@@ -194,8 +189,5 @@ int main(int argc, char* argv[])
     std::for_each(all_sfc_nodes.begin(), all_sfc_nodes.end(),
                   std::default_delete<MeshLib::Node>());
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshGeoTools/constructMeshesFromGeometry.cpp b/Applications/Utils/MeshGeoTools/constructMeshesFromGeometry.cpp
index d57974e0476a652f085eda4f5b27e36141c85477..9d99fc3e19d019c7bc758e9f668b57d0a5b3c33e 100644
--- a/Applications/Utils/MeshGeoTools/constructMeshesFromGeometry.cpp
+++ b/Applications/Utils/MeshGeoTools/constructMeshesFromGeometry.cpp
@@ -10,6 +10,7 @@
 
 #include <tclap/CmdLine.h>
 
+#include "BaseLib/MPI.h"
 #include "GeoLib/GEOObjects.h"
 #include "GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.h"
 #include "InfoLib/GitInfo.h"
@@ -19,9 +20,6 @@
 #include "MeshLib/IO/writeMeshToFile.h"
 #include "MeshLib/Mesh.h"
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
 std::unique_ptr<GeoLib::GEOObjects> readGeometry(std::string const& filename)
 {
     auto geo_objects = std::make_unique<GeoLib::GEOObjects>();
@@ -80,9 +78,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
     std::unique_ptr<MeshLib::Mesh> mesh{
         MeshLib::IO::readMeshFromFile(mesh_arg.getValue())};
 
@@ -101,9 +97,6 @@ int main(int argc, char* argv[])
         if (!m_ptr)
         {
             ERR("Could not create a mesh for each given geometry.");
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
         if (m_ptr->getNodes().empty())
@@ -117,8 +110,5 @@ int main(int argc, char* argv[])
         MeshLib::IO::writeMeshToFile(*m_ptr, m_ptr->getName() + ".vtu");
     }
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshGeoTools/createIntermediateRasters.cpp b/Applications/Utils/MeshGeoTools/createIntermediateRasters.cpp
index 2a98a25c57d1f6e2231eeba2d9dacf04b82044eb..75a7572a9641c3e81362b0b49899816277e9c509 100644
--- a/Applications/Utils/MeshGeoTools/createIntermediateRasters.cpp
+++ b/Applications/Utils/MeshGeoTools/createIntermediateRasters.cpp
@@ -9,14 +9,11 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <memory>
 #include <string>
 
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "GeoLib/IO/AsciiRasterInterface.h"
 #include "GeoLib/Raster.h"
 #include "InfoLib/GitInfo.h"
@@ -50,9 +47,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::unique_ptr<GeoLib::Raster> dem1(
         FileIO::AsciiRasterInterface::readRaster(input1_arg.getValue()));
@@ -61,9 +56,6 @@ int main(int argc, char* argv[])
 
     if (dem1 == nullptr || dem2 == nullptr)
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return 1;
     }
 
@@ -99,9 +91,6 @@ int main(int argc, char* argv[])
 
     if (errors_found)
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return 2;
     }
 
@@ -120,9 +109,6 @@ int main(int argc, char* argv[])
         if (it2 == dem2->end())
         {
             ERR("Error: File 2 is shorter than File 1.");
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return 3;
         }
         if (*it1 == h1.no_data || *it2 == h2.no_data)
@@ -147,9 +133,6 @@ int main(int argc, char* argv[])
     if (it2 != dem2->end())
     {
         ERR("Error: File 1 is shorter than File 2.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return 3;
     }
 
@@ -164,8 +147,5 @@ int main(int argc, char* argv[])
             r, basename + std::to_string(i) + ext);
         INFO("Layer {:d} written.", i + 1);
     }
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshGeoTools/geometryToGmshGeo.cpp b/Applications/Utils/MeshGeoTools/geometryToGmshGeo.cpp
index 87777c243a7af24430cd51a2b8a058458fb37ff0..87d735950a42b17ee9d2f8c377cdb1d586d0c717 100644
--- a/Applications/Utils/MeshGeoTools/geometryToGmshGeo.cpp
+++ b/Applications/Utils/MeshGeoTools/geometryToGmshGeo.cpp
@@ -12,11 +12,8 @@
 // ThirdParty
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include "Applications/FileIO/Gmsh/GMSHInterface.h"
+#include "BaseLib/MPI.h"
 #include "GeoLib/GEOObjects.h"
 #include "GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.h"
 #include "InfoLib/GitInfo.h"
@@ -70,9 +67,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     GeoLib::GEOObjects geo_objects;
     for (auto const& geometry_name : geo_input_arg.getValue())
@@ -82,9 +77,6 @@ int main(int argc, char* argv[])
         {
             if (!xml.readFile(geometry_name))
             {
-#ifdef USE_PETSC
-                MPI_Finalize();
-#endif
                 return EXIT_FAILURE;
             }
         }
@@ -92,9 +84,6 @@ int main(int argc, char* argv[])
         {
             ERR("Failed to read file '{:s}'.", geometry_name);
             ERR("{:s}", err.what());
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
         INFO("Successfully read file '{:s}'.", geometry_name);
@@ -155,13 +144,7 @@ int main(int argc, char* argv[])
                 "for better analysis of the problem.");
         }
         ERR("An error has occurred - programme will be terminated.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshGeoTools/identifySubdomains.cpp b/Applications/Utils/MeshGeoTools/identifySubdomains.cpp
index 46e44c801d5730e72c6e7bca081dec9d4fcc3254..fcfba24983f80302dcadcd15e2446d8e7fc444ab 100644
--- a/Applications/Utils/MeshGeoTools/identifySubdomains.cpp
+++ b/Applications/Utils/MeshGeoTools/identifySubdomains.cpp
@@ -10,10 +10,7 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
+#include "BaseLib/MPI.h"
 #include "BaseLib/RunTime.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshGeoToolsLib/IdentifySubdomainMesh.h"
@@ -98,9 +95,7 @@ int main(int argc, char* argv[])
     cmd.add(subdomain_meshes_filenames_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     //
     // The bulk mesh.
@@ -165,9 +160,6 @@ int main(int argc, char* argv[])
     }
     INFO("writing time: {:g} s", writing_time.elapsed());
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     INFO("Entire run time: {:g} s", run_time.elapsed());
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/ModelPreparation/ComputeNodeAreasFromSurfaceMesh.cpp b/Applications/Utils/ModelPreparation/ComputeNodeAreasFromSurfaceMesh.cpp
index c29a8bc4bbf82d1f7c42efb6c041d08f63935835..9bad7b3970af3a0b43213c6f93d088e14e957e85 100644
--- a/Applications/Utils/ModelPreparation/ComputeNodeAreasFromSurfaceMesh.cpp
+++ b/Applications/Utils/ModelPreparation/ComputeNodeAreasFromSurfaceMesh.cpp
@@ -11,10 +11,6 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <fstream>
 #include <memory>
 #include <numeric>
@@ -23,6 +19,7 @@
 
 #include "BaseLib/Error.h"
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/readMeshFromFile.h"
 #include "MeshLib/Mesh.h"
@@ -90,9 +87,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::unique_ptr<MeshLib::Mesh> surface_mesh(
         MeshLib::IO::readMeshFromFile(mesh_in.getValue()));
@@ -112,9 +107,6 @@ int main(int argc, char* argv[])
         if (!orig_node_ids)
         {
             ERR("Fatal error: could not create property.");
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
         orig_node_ids->resize(surface_mesh->getNumberOfNodes());
@@ -149,8 +141,5 @@ int main(int argc, char* argv[])
     writeToFile(id_and_area_fname, csv_fname, ids_and_areas,
                 surface_mesh->getNodes());
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/ModelPreparation/PartitionMesh/BinaryToPVTU.cpp b/Applications/Utils/ModelPreparation/PartitionMesh/BinaryToPVTU.cpp
index fd7b3baca9c08d252051cc6f872c28455a526646..81ce322f9ec6a42d72dc9d717a505db3dbf3981f 100644
--- a/Applications/Utils/ModelPreparation/PartitionMesh/BinaryToPVTU.cpp
+++ b/Applications/Utils/ModelPreparation/PartitionMesh/BinaryToPVTU.cpp
@@ -10,7 +10,6 @@
 
 */
 
-#include <mpi.h>
 #include <spdlog/spdlog.h>
 #include <tclap/CmdLine.h>
 #include <vtkMPIController.h>
@@ -18,6 +17,7 @@
 
 #include "BaseLib/CPUTime.h"
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "BaseLib/RunTime.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/VtkIO/VtuInterface.h"
@@ -72,9 +72,7 @@ int main(int argc, char* argv[])
             OGS_FATAL("spdlog logger error occurred.");
         });
 
-    // init MPI
-    MPI_Init(&argc, &argv);
-
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
     // start the timer
     BaseLib::RunTime run_timer;
     run_timer.start();
@@ -111,8 +109,5 @@ int main(int argc, char* argv[])
 
     INFO("Total runtime: {:g} s.", run_timer.elapsed());
     INFO("Total CPU time: {:g} s.", CPU_timer.elapsed());
-
-    MPI_Finalize();
-
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/ModelPreparation/PartitionMesh/PartitionMesh.cpp b/Applications/Utils/ModelPreparation/PartitionMesh/PartitionMesh.cpp
index daa2aa43cd3ea0fe8e9de0e1c61a12f2119dc7ef..3e5316313857119e63917dd33361787ea94387a8 100644
--- a/Applications/Utils/ModelPreparation/PartitionMesh/PartitionMesh.cpp
+++ b/Applications/Utils/ModelPreparation/PartitionMesh/PartitionMesh.cpp
@@ -15,12 +15,9 @@
 #include <spdlog/spdlog.h>
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include "BaseLib/CPUTime.h"
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "BaseLib/RunTime.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/readMeshFromFile.h"
@@ -102,9 +99,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     BaseLib::setConsoleLogLevel(log_level_arg.getValue());
     spdlog::set_pattern("%^%l:%$ %v");
@@ -144,9 +139,6 @@ int main(int argc, char* argv[])
         INFO("Total runtime: {:g} s.", run_timer.elapsed());
         INFO("Total CPU time: {:g} s.", CPU_timer.elapsed());
 
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_SUCCESS;
     }
 
@@ -192,9 +184,6 @@ int main(int argc, char* argv[])
         {
             INFO("Failed in system calling.");
             INFO("Return value of system call {:d} ", status);
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
     }
@@ -244,8 +233,5 @@ int main(int argc, char* argv[])
     INFO("Total runtime: {:g} s.", run_timer.elapsed());
     INFO("Total CPU time: {:g} s.", CPU_timer.elapsed());
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/ModelPreparation/convertVtkDataArrayToVtkDataArray.cpp b/Applications/Utils/ModelPreparation/convertVtkDataArrayToVtkDataArray.cpp
index dbb156149c02168c9920bfd374837e8d0e198f85..0fd23a6fe3c1ced13434a4e8591c03562737bb49 100644
--- a/Applications/Utils/ModelPreparation/convertVtkDataArrayToVtkDataArray.cpp
+++ b/Applications/Utils/ModelPreparation/convertVtkDataArrayToVtkDataArray.cpp
@@ -11,15 +11,12 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <algorithm>
 #include <cmath>
 #include <memory>
 #include <numeric>
 
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/readMeshFromFile.h"
 #include "MeshLib/IO/writeMeshToFile.h"
@@ -105,18 +102,13 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::unique_ptr<MeshLib::Mesh> mesh(
         MeshLib::IO::readMeshFromFile(mesh_arg.getValue()));
 
     if (!mesh)
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return -1;
     }
 
@@ -163,16 +155,10 @@ int main(int argc, char* argv[])
     if (!success)
     {
         ERR("{:s}", err_msg);
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return -1;
     }
 
     MeshLib::IO::writeMeshToFile(*mesh, out_mesh_arg.getValue());
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/ModelPreparation/createNeumannBc.cpp b/Applications/Utils/ModelPreparation/createNeumannBc.cpp
index 7f108478b1f70e54b7fa99ae67114c4d8dea7e28..b69f69c48613b1cd25ffd64797e59e74c9feb4d0 100644
--- a/Applications/Utils/ModelPreparation/createNeumannBc.cpp
+++ b/Applications/Utils/ModelPreparation/createNeumannBc.cpp
@@ -10,12 +10,9 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <fstream>
 
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/Elements/Element.h"
 #include "MeshLib/IO/readMeshFromFile.h"
@@ -130,9 +127,7 @@ int main(int argc, char* argv[])
     cmd.add(result_file);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     // read surface mesh
     std::unique_ptr<MeshLib::Mesh> surface_mesh(
@@ -155,9 +150,6 @@ int main(int argc, char* argv[])
     }();
     if (!node_id_pv)
     {
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
 
@@ -192,8 +184,5 @@ int main(int argc, char* argv[])
         result_out << p.first << " " << p.second << "\n";
     }
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/ModelPreparation/scaleProperty.cpp b/Applications/Utils/ModelPreparation/scaleProperty.cpp
index 03eb59846b759d122c7ccdbef0efcd9c5edd0c7d..562ee1d9212647d2ec86dc2ffadab0f53fe88d51 100644
--- a/Applications/Utils/ModelPreparation/scaleProperty.cpp
+++ b/Applications/Utils/ModelPreparation/scaleProperty.cpp
@@ -11,15 +11,12 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <algorithm>
 #include <cmath>
 #include <memory>
 #include <numeric>
 
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/readMeshFromFile.h"
 #include "MeshLib/IO/writeMeshToFile.h"
@@ -72,9 +69,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::unique_ptr<MeshLib::Mesh> mesh(
         MeshLib::IO::readMeshFromFile(mesh_arg.getValue()));
@@ -97,8 +92,5 @@ int main(int argc, char* argv[])
 
     MeshLib::IO::writeMeshToFile(*mesh, out_mesh_arg.getValue());
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/PostProcessing/Raster2PointCloud.cpp b/Applications/Utils/PostProcessing/Raster2PointCloud.cpp
index 2b22e4fde07524ff0cad6ca773068e6ea674cc42..44fa7de536fe4567393b802aeb96f573063aac21 100644
--- a/Applications/Utils/PostProcessing/Raster2PointCloud.cpp
+++ b/Applications/Utils/PostProcessing/Raster2PointCloud.cpp
@@ -13,9 +13,6 @@
 
 // ThirdParty
 #include <tclap/CmdLine.h>
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
 #include <vtkPolyDataAlgorithm.h>
 #include <vtkSmartPointer.h>
 #include <vtkXMLPolyDataWriter.h>
@@ -23,6 +20,7 @@
 #include "Applications/DataExplorer/VtkVis/VtkGeoImageSource.h"
 #include "Applications/DataExplorer/VtkVis/VtkImageDataToPointCloudFilter.h"
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "GeoLib/IO/AsciiRasterInterface.h"
 #include "InfoLib/GitInfo.h"
 
@@ -80,9 +78,7 @@ int main(int argc, char* argv[])
     cmd.add(input_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     std::string const input_name = input_arg.getValue().c_str();
     std::string const output_name = output_arg.getValue().c_str();
@@ -161,8 +157,5 @@ int main(int argc, char* argv[])
         }
     }
     std::cout << "done." << std::endl;
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/PostProcessing/postLIE.cpp b/Applications/Utils/PostProcessing/postLIE.cpp
index 86199ac4055ea3f0b34b42d7a1b655d95fef5be4..dac283371f5843970749be49b18adba01ae05e4d 100644
--- a/Applications/Utils/PostProcessing/postLIE.cpp
+++ b/Applications/Utils/PostProcessing/postLIE.cpp
@@ -9,10 +9,6 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <boost/property_tree/ptree.hpp>
 #include <boost/property_tree/xml_parser.hpp>
 #include <map>
@@ -20,6 +16,7 @@
 #include <vector>
 
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "InfoLib/GitInfo.h"
 #include "MeshLib/IO/readMeshFromFile.h"
 #include "MeshLib/IO/writeMeshToFile.h"
@@ -151,9 +148,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     auto const in_file_ext = BaseLib::getFileExtension(arg_in_file.getValue());
     if (in_file_ext == ".pvd")
@@ -170,8 +165,5 @@ int main(int argc, char* argv[])
         OGS_FATAL("The given file type ({:s}) is not supported.", in_file_ext);
     }
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/SWMMConverter/SWMMConverter.cpp b/Applications/Utils/SWMMConverter/SWMMConverter.cpp
index 7db44c2b4afa29462b28f4c25f050fc2018a8235..be95e7c61282d27ab77600e89affc6404c69bb4e 100644
--- a/Applications/Utils/SWMMConverter/SWMMConverter.cpp
+++ b/Applications/Utils/SWMMConverter/SWMMConverter.cpp
@@ -9,12 +9,9 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include "Applications/FileIO/SWMM/SWMMInterface.h"
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "BaseLib/StringTools.h"
 #include "GeoLib/GEOObjects.h"
 #include "GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.h"
@@ -199,18 +196,13 @@ int main(int argc, char* argv[])
     cmd.add(add_system_arg);
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     if (!(geo_output_arg.isSet() || mesh_output_arg.isSet() ||
           csv_output_arg.isSet()))
     {
         ERR("No output format given. Please specify OGS geometry or mesh "
             "output file.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return -1;
     }
 
@@ -219,9 +211,6 @@ int main(int argc, char* argv[])
     {
         ERR("Please specify csv output file for exporting subcatchment or "
             "system parameters.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return -1;
     }
 
@@ -240,8 +229,5 @@ int main(int argc, char* argv[])
                        add_subcatchments_arg.getValue(),
                        add_system_arg.getValue());
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return 0;
 }
diff --git a/Applications/Utils/SimpleMeshCreation/createMeshElemPropertiesFromASCRaster.cpp b/Applications/Utils/SimpleMeshCreation/createMeshElemPropertiesFromASCRaster.cpp
index 9039a44a1349f43e1b3f1167ac065e3d2c16fe83..00b7dfedf5a1e2a10c2f356a28f14b79655450ff 100644
--- a/Applications/Utils/SimpleMeshCreation/createMeshElemPropertiesFromASCRaster.cpp
+++ b/Applications/Utils/SimpleMeshCreation/createMeshElemPropertiesFromASCRaster.cpp
@@ -12,16 +12,13 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <algorithm>
 #include <cmath>
 #include <memory>
 #include <numeric>
 
 #include "BaseLib/FileTools.h"
+#include "BaseLib/MPI.h"
 #include "BaseLib/quicksort.h"
 #include "GeoLib/IO/AsciiRasterInterface.h"
 #include "GeoLib/Raster.h"
@@ -94,9 +91,7 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
 
     // read mesh
     std::unique_ptr<MeshLib::Mesh> dest_mesh(
@@ -138,8 +133,5 @@ int main(int argc, char* argv[])
         MeshLib::IO::writeMeshToFile(*dest_mesh, out_mesh_arg.getValue());
     }
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/SimpleMeshCreation/generateStructuredMesh.cpp b/Applications/Utils/SimpleMeshCreation/generateStructuredMesh.cpp
index 7e3062698a32deaf00debee231fe8011de5fe172..0392559aae06ae71e09f583e7aedac88ec4d206e 100644
--- a/Applications/Utils/SimpleMeshCreation/generateStructuredMesh.cpp
+++ b/Applications/Utils/SimpleMeshCreation/generateStructuredMesh.cpp
@@ -9,16 +9,13 @@
 
 #include <tclap/CmdLine.h>
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
 #include <algorithm>
 #include <memory>
 #include <string>
 #include <vector>
 
 #include "BaseLib/Error.h"
+#include "BaseLib/MPI.h"
 #include "BaseLib/Subdivision.h"
 #include "BaseLib/TCLAPCustomOutput.h"
 #include "InfoLib/GitInfo.h"
@@ -161,9 +158,7 @@ int main(int argc, char* argv[])
 
     // parse arguments
     cmd.parse(argc, argv);
-#ifdef USE_PETSC
-    MPI_Init(&argc, &argv);
-#endif
+    BaseLib::MPI::Setup mpi_setup(argc, argv);
     const std::string eleTypeName(eleTypeArg.getValue());
     const MeshLib::MeshElemType eleType =
         MeshLib::String2MeshElemType(eleTypeName);
@@ -194,9 +189,6 @@ int main(int argc, char* argv[])
     if (!isLengthSet)
     {
         ERR("Missing input: Length information is not provided at all.");
-#ifdef USE_PETSC
-        MPI_Finalize();
-#endif
         return EXIT_FAILURE;
     }
     for (unsigned i = 0; i < 3; i++)
@@ -206,9 +198,6 @@ int main(int argc, char* argv[])
             ERR("Missing input: Length for dimension [{:d}] is required but "
                 "missing.",
                 i);
-#ifdef USE_PETSC
-            MPI_Finalize();
-#endif
             return EXIT_FAILURE;
         }
     }
@@ -304,8 +293,5 @@ int main(int argc, char* argv[])
             *(mesh.get()), std::filesystem::path(mesh_out.getValue()));
     }
 
-#ifdef USE_PETSC
-    MPI_Finalize();
-#endif
     return EXIT_SUCCESS;
 }
diff --git a/BaseLib/Algorithm.h b/BaseLib/Algorithm.h
index 371f14622cd6648338cca23e9ea0e26c33e90965..28ee2516578e0a6217af93d3c345cc701fe5f950 100644
--- a/BaseLib/Algorithm.h
+++ b/BaseLib/Algorithm.h
@@ -16,6 +16,10 @@
 #include <optional>
 #include <range/v3/algorithm/find_if.hpp>
 #include <range/v3/range/concepts.hpp>
+#include <range/v3/range/conversion.hpp>
+#include <range/v3/view/concat.hpp>
+#include <range/v3/view/partial_sum.hpp>
+#include <range/v3/view/single.hpp>
 #include <string>
 #include <typeindex>
 #include <typeinfo>
@@ -272,6 +276,18 @@ void cleanupVectorElements(std::vector<T1*>& dependent_items, Args&&... args)
     cleanupVectorElements(std::forward<Args>(args)...);
 }
 
+/// Converts range of sizes to a vector of offsets. First offset is 0 and the
+/// resulting vector size is the size of the input range plus one.
+template <ranges::range R>
+    requires std::is_integral_v<ranges::range_value_t<R>>
+std::vector<ranges::range_value_t<R>> sizesToOffsets(R const& sizes)
+{
+    return ranges::views::concat(
+               ranges::views::single(ranges::range_value_t<R>{0}),
+               ranges::views::partial_sum(sizes)) |
+           ranges::to<std::vector<ranges::range_value_t<R>>>();
+}
+
 /// Checks if any of the elements in the given list is true.
 template <typename List>
 constexpr bool any_of(List const& values)
diff --git a/BaseLib/CMakeLists.txt b/BaseLib/CMakeLists.txt
index 3af91ae9e084128940a9071ecd01b3c597b4008a..34dea6cfceed6e229e8da5580e4971c3ffb66938 100644
--- a/BaseLib/CMakeLists.txt
+++ b/BaseLib/CMakeLists.txt
@@ -12,6 +12,7 @@ target_link_libraries(
     BaseLib
     PUBLIC Boost::algorithm
            Boost::property_tree
+           range-v3
            spdlog
            tclap
            $<$<BOOL:${MSVC}>:WinMM> # needed for timeGetTime
diff --git a/BaseLib/MPI.h b/BaseLib/MPI.h
index 3b37773e3e976ab5298ee9f83304dd86f1bfce68..649769740e615b7a0c6a4116a230fc0144e08065 100644
--- a/BaseLib/MPI.h
+++ b/BaseLib/MPI.h
@@ -5,28 +5,196 @@
  *            Distributed under a Modified BSD License.
  *              See accompanying file LICENSE.txt or
  *              http://www.opengeosys.org/project/license
- *
  */
 
 #pragma once
 
+#include <algorithm>
+
+#include "Algorithm.h"
+#include "Error.h"
+
+#ifdef USE_PETSC
+#include <mpi.h>
+#endif
+
 namespace BaseLib::MPI
 {
 
+struct Setup
+{
+    Setup(int argc, char* argv[])
+    {
+#ifdef USE_PETSC
+        MPI_Init(&argc, &argv);
+#else
+        (void)argc;
+        (void)argv;
+#endif  // USE_PETSC
+    }
+
+    ~Setup()
+    {
+#ifdef USE_PETSC
+        MPI_Finalize();
+#endif  // USE_PETSC
+    }
+};
+
 #ifdef USE_PETSC
-// Reduce operations for interprocess communications while using Petsc
-static inline int reduceMin(int val)
+struct Mpi
+{
+    Mpi(MPI_Comm const communicator = MPI_COMM_WORLD)
+        : communicator(communicator)
+    {
+        int mpi_init;
+        MPI_Initialized(&mpi_init);
+        if (mpi_init != 1)
+        {
+            OGS_FATAL("MPI is not initialized.");
+        }
+        MPI_Comm_size(communicator, &size);
+        MPI_Comm_rank(communicator, &rank);
+    }
+
+    MPI_Comm communicator;
+    int size;
+    int rank;
+};
+
+template <typename T>
+constexpr MPI_Datatype mpiType()
 {
-    int result;
-    MPI_Allreduce(&val, &result, 1, MPI_INTEGER, MPI_MIN, PETSC_COMM_WORLD);
+    using U = std::remove_const_t<T>;
+    if constexpr (std::is_same_v<U, bool>)
+    {
+        return MPI_C_BOOL;
+    }
+    if constexpr (std::is_same_v<U, char>)
+    {
+        return MPI_CHAR;
+    }
+    if constexpr (std::is_same_v<U, double>)
+    {
+        return MPI_DOUBLE;
+    }
+    if constexpr (std::is_same_v<U, float>)
+    {
+        return MPI_FLOAT;
+    }
+    if constexpr (std::is_same_v<U, int>)
+    {
+        return MPI_INT;
+    }
+    if constexpr (std::is_same_v<U, std::size_t>)
+    {
+        return MPI_UNSIGNED_LONG;
+    }
+    if constexpr (std::is_same_v<U, unsigned int>)
+    {
+        return MPI_UNSIGNED;
+    }
+}
+
+template <typename T>
+static std::vector<T> allgather(T const& value, Mpi const& mpi)
+{
+    std::vector<T> result(mpi.size);
+
+    result[mpi.rank] = value;
+
+    MPI_Allgather(&result[mpi.rank], 1, mpiType<T>(), result.data(), 1,
+                  mpiType<T>(), mpi.communicator);
+
     return result;
 }
-#else
-// Reduce operations for interprocess communications without using Petsc
-static inline int reduceMin(int val)
+
+template <typename T>
+static std::vector<T> allgather(std::vector<T> const& vector, Mpi const& mpi)
 {
-    return val;
+    std::size_t const size = vector.size();
+    // Flat in memory over all ranks;
+    std::vector<T> result(mpi.size * size);
+
+    std::copy_n(vector.begin(), size, &result[mpi.rank * size]);
+
+    MPI_Allgather(&result[mpi.rank * size], size, mpiType<T>(), result.data(),
+                  size, mpiType<T>(), mpi.communicator);
+
+    return result;
+}
+
+template <typename T>
+static T allreduce(T const& value, MPI_Op const& mpi_op, Mpi const& mpi)
+{
+    T result{};
+
+    MPI_Allreduce(&value, &result, 1, mpiType<T>(), mpi_op, mpi.communicator);
+    return result;
+}
+
+template <typename T>
+static std::vector<T> allreduce(std::vector<T> const& vector,
+                                MPI_Op const& mpi_op, Mpi const& mpi)
+{
+    std::size_t const size = vector.size();
+    std::vector<T> result(vector.size());
+
+    MPI_Allreduce(vector.data(), result.data(), size, mpiType<T>(), mpi_op,
+                  mpi.communicator);
+    return result;
+}
+
+template <typename T>
+static void allreduceInplace(std::vector<T>& vector,
+                             MPI_Op const& mpi_op,
+                             Mpi const& mpi)
+{
+    MPI_Allreduce(MPI_IN_PLACE,
+                  vector.data(),
+                  vector.size(),
+                  mpiType<T>(),
+                  mpi_op,
+                  mpi.communicator);
+}
+
+/// Allgather for variable data. Returns offsets in the receive buffer.
+/// The receive buffer is resized to accommodate the gathered data.
+template <typename T>
+static std::vector<int> allgatherv(
+    std::span<T> const send_buffer,
+    std::vector<std::remove_const_t<T>>& receive_buffer,
+    Mpi const& mpi)
+{
+    // Determine the number of elements to send
+    int const size = static_cast<int>(send_buffer.size());
+
+    // Gather sizes from all ranks
+    std::vector<int> const sizes = allgather(size, mpi);
+
+    // Compute offsets based on counts
+    std::vector<int> const offsets = BaseLib::sizesToOffsets(sizes);
+
+    // Resize receive buffer to hold all gathered data
+    receive_buffer.resize(offsets.back());
+
+    MPI_Allgatherv(send_buffer.data(), size, mpiType<T>(),
+                   receive_buffer.data(), sizes.data(), offsets.data(),
+                   mpiType<T>(), mpi.communicator);
+
+    return offsets;
 }
 #endif
 
+/// The reduction is implemented transparently for with and without MPI. In the
+/// latter case the input value is returned.
+static inline int reduceMin(int const val)
+{
+#ifdef USE_PETSC
+    return allreduce(val, MPI_MIN, Mpi{MPI_COMM_WORLD});
+#else
+    return val;
+#endif
+}
+
 }  // namespace BaseLib::MPI
diff --git a/MathLib/LinAlg/PETSc/PETScVector.cpp b/MathLib/LinAlg/PETSc/PETScVector.cpp
index bf12c6d3878ce03005e4b6b2abb060c76a1be842..a6201f116d811a352f280896be6bab9d1d3b6e8a 100644
--- a/MathLib/LinAlg/PETSc/PETScVector.cpp
+++ b/MathLib/LinAlg/PETSc/PETScVector.cpp
@@ -23,7 +23,9 @@
 #include <algorithm>
 #include <cassert>
 
+#include "BaseLib/Algorithm.h"
 #include "BaseLib/Error.h"
+#include "BaseLib/MPI.h"
 
 namespace MathLib
 {
@@ -113,33 +115,6 @@ void PETScVector::finalizeAssembly()
     VecAssemblyEnd(v_);
 }
 
-void PETScVector::gatherLocalVectors(PetscScalar local_array[],
-                                     PetscScalar global_array[]) const
-{
-    // Collect vectors from processors.
-    int size_rank;
-    MPI_Comm_size(PETSC_COMM_WORLD, &size_rank);
-
-    // number of elements to be sent for each rank
-    std::vector<PetscInt> i_cnt(size_rank);
-
-    MPI_Allgather(&size_loc_, 1, MPI_INT, &i_cnt[0], 1, MPI_INT,
-                  PETSC_COMM_WORLD);
-
-    // collect local array
-    PetscInt offset = 0;
-    // offset in the receive vector of the data from each rank
-    std::vector<PetscInt> i_disp(size_rank);
-    for (PetscInt i = 0; i < size_rank; i++)
-    {
-        i_disp[i] = offset;
-        offset += i_cnt[i];
-    }
-
-    MPI_Allgatherv(local_array, size_loc_, MPI_DOUBLE, global_array, &i_cnt[0],
-                   &i_disp[0], MPI_DOUBLE, PETSC_COMM_WORLD);
-}
-
 void PETScVector::getGlobalVector(std::vector<PetscScalar>& u) const
 {
 #ifdef TEST_MEM_PETSC
@@ -159,7 +134,8 @@ void PETScVector::getGlobalVector(std::vector<PetscScalar>& u) const
     PetscScalar* xp = nullptr;
     VecGetArray(v_, &xp);
 
-    gatherLocalVectors(xp, u.data());
+    BaseLib::MPI::Mpi mpi{PETSC_COMM_WORLD};
+    BaseLib::MPI::allgatherv(std::span(xp, size_loc_), u, mpi);
 
     // This following line may be needed late on
     //  for a communication load balance:
diff --git a/MathLib/LinAlg/PETSc/PETScVector.h b/MathLib/LinAlg/PETSc/PETScVector.h
index ed2734878fb9398b903ee37e3a4de4c1c6b69ed4..20983dacf2d3edcf875439f275e53afd505998bd 100644
--- a/MathLib/LinAlg/PETSc/PETScVector.h
+++ b/MathLib/LinAlg/PETSc/PETScVector.h
@@ -267,14 +267,6 @@ private:
     /// Map global indices of ghost entries to local indices
     mutable std::map<PetscInt, PetscInt> global_ids2local_ids_ghost_;
 
-    /*!
-          \brief  Collect local vectors
-          \param  local_array Local array
-          \param  global_array Global array
-    */
-    void gatherLocalVectors(PetscScalar local_array[],
-                            PetscScalar global_array[]) const;
-
     /*!
        Get local vector, i.e. entries in the same rank
     */
diff --git a/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp b/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp
index 3979edd184fc7dacac1dd48856b2557acf4e8202..56e83a22c75d642cdf451eeaf673168bddec0e93 100644
--- a/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp
+++ b/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp
@@ -9,10 +9,7 @@
 
 #include "BoundaryElementsAtPoint.h"
 
-#ifdef USE_PETSC
-#include <mpi.h>
-#endif
-
+#include "BaseLib/MPI.h"
 #include "GeoLib/Point.h"
 #include "MathLib/Point3d.h"
 #include "MeshGeoToolsLib/MeshNodeSearcher.h"
@@ -32,10 +29,8 @@ BoundaryElementsAtPoint::BoundaryElementsAtPoint(
 
 #ifdef USE_PETSC
     std::size_t const number_of_found_nodes_at_rank = node_ids.size();
-    std::size_t number_of_total_found_nodes = 0;
-
-    MPI_Allreduce(&number_of_found_nodes_at_rank, &number_of_total_found_nodes,
-                  1, MPI_UNSIGNED_LONG, MPI_SUM, MPI_COMM_WORLD);
+    std::size_t const number_of_total_found_nodes = BaseLib::MPI::allreduce(
+        number_of_found_nodes_at_rank, MPI_SUM, BaseLib::MPI::Mpi{});
 
     if (number_of_total_found_nodes == 0)
     {
diff --git a/MeshLib/IO/MPI_IO/NodePartitionedMeshReader.cpp b/MeshLib/IO/MPI_IO/NodePartitionedMeshReader.cpp
index ed7000c22868f2224b518e6645ac4f3033205c0e..faf10327ae865bc7c2bf1b180816d1ddc5b28df0 100644
--- a/MeshLib/IO/MPI_IO/NodePartitionedMeshReader.cpp
+++ b/MeshLib/IO/MPI_IO/NodePartitionedMeshReader.cpp
@@ -22,6 +22,7 @@
 
 #include "BaseLib/FileTools.h"
 #include "BaseLib/Logging.h"
+#include "BaseLib/MPI.h"
 #include "BaseLib/RunTime.h"
 #include "MeshLib/Elements/Elements.h"
 #include "MeshLib/MeshEnums.h"
@@ -45,12 +46,8 @@ namespace MeshLib
 {
 namespace IO
 {
-NodePartitionedMeshReader::NodePartitionedMeshReader(MPI_Comm comm)
-    : _mpi_comm(comm)
+NodePartitionedMeshReader::NodePartitionedMeshReader(MPI_Comm comm) : mpi_(comm)
 {
-    MPI_Comm_size(_mpi_comm, &_mpi_comm_size);
-    MPI_Comm_rank(_mpi_comm, &_mpi_rank);
-
     registerNodeDataMpiType();
 }
 
@@ -79,13 +76,13 @@ MeshLib::NodePartitionedMesh* NodePartitionedMeshReader::read(
 
     // Always try binary file first
     std::string const fname_new = file_name_base + "_partitioned_msh_cfg" +
-                                  std::to_string(_mpi_comm_size) + ".bin";
+                                  std::to_string(mpi_.size) + ".bin";
 
     if (!BaseLib::IsFileExisting(fname_new))  // binary file does not exist.
     {
         std::string const fname_ascii = file_name_base +
                                         "_partitioned_msh_cfg" +
-                                        std::to_string(_mpi_comm_size) + ".msh";
+                                        std::to_string(mpi_.size) + ".msh";
         if (BaseLib::IsFileExisting(fname_ascii))
         {
             ERR("Reading of ASCII mesh file {:s} is not supported since OGS "
@@ -102,7 +99,7 @@ MeshLib::NodePartitionedMesh* NodePartitionedMeshReader::read(
 
     INFO("[time] Reading the mesh took {:f} s.", timer.elapsed());
 
-    MPI_Barrier(_mpi_comm);
+    MPI_Barrier(mpi_.communicator);
 
     return mesh;
 }
@@ -124,8 +121,9 @@ bool NodePartitionedMeshReader::readDataFromFile(std::string const& filename,
     MPI_File file;
 
     char* filename_char = const_cast<char*>(filename.data());
-    int const file_status = MPI_File_open(
-        _mpi_comm, filename_char, MPI_MODE_RDONLY, MPI_INFO_NULL, &file);
+    int const file_status =
+        MPI_File_open(mpi_.communicator, filename_char, MPI_MODE_RDONLY,
+                      MPI_INFO_NULL, &file);
 
     if (file_status != 0)
     {
@@ -151,13 +149,13 @@ MeshLib::NodePartitionedMesh* NodePartitionedMeshReader::readMesh(
     //----------------------------------------------------------------------------------
     // Read headers
     const std::string fname_header = file_name_base + "_partitioned_msh_";
-    const std::string fname_num_p_ext = std::to_string(_mpi_comm_size) + ".bin";
+    const std::string fname_num_p_ext = std::to_string(mpi_.size) + ".bin";
 
     // Read the config meta data from *cfg* file into struct PartitionedMeshInfo
     // _mesh_info
     if (!readDataFromFile(
             fname_header + "cfg" + fname_num_p_ext,
-            static_cast<MPI_Offset>(static_cast<unsigned>(_mpi_rank) *
+            static_cast<MPI_Offset>(static_cast<unsigned>(mpi_.rank) *
                                     sizeof(_mesh_info)),
             MPI_LONG, _mesh_info))
         return nullptr;
@@ -228,7 +226,7 @@ void NodePartitionedMeshReader::readProperties(
 
     const std::string fname_cfg = file_name_base + "_partitioned_" + item_type +
                                   "_properties_cfg" +
-                                  std::to_string(_mpi_comm_size) + ".bin";
+                                  std::to_string(mpi_.size) + ".bin";
     std::ifstream is(fname_cfg.c_str(), std::ios::binary | std::ios::in);
     if (!is)
     {
@@ -262,21 +260,19 @@ void NodePartitionedMeshReader::readProperties(
     auto pos = is.tellg();
     auto offset =
         static_cast<long>(pos) +
-        static_cast<long>(_mpi_rank *
+        static_cast<long>(mpi_.rank *
                           sizeof(MeshLib::IO::PropertyVectorPartitionMetaData));
     is.seekg(offset);
     std::optional<MeshLib::IO::PropertyVectorPartitionMetaData> pvpmd(
         MeshLib::IO::readPropertyVectorPartitionMetaData(is));
-    bool pvpmd_read_ok = static_cast<bool>(pvpmd);
-    bool all_pvpmd_read_ok;
-    MPI_Allreduce(&pvpmd_read_ok, &all_pvpmd_read_ok, 1, MPI_C_BOOL, MPI_LOR,
-                  _mpi_comm);
+    bool const all_pvpmd_read_ok =
+        BaseLib::MPI::allreduce(static_cast<bool>(pvpmd), MPI_LOR, mpi_);
     if (!all_pvpmd_read_ok)
     {
         OGS_FATAL(
             "Error in NodePartitionedMeshReader::readProperties: "
             "Could not read the partition meta data for the mpi process {:d}",
-            _mpi_rank);
+            mpi_.rank);
     }
     DBUG("offset in the PropertyVector: {:d}", pvpmd->offset);
     DBUG("{:d} tuples in partition.", pvpmd->number_of_tuples);
@@ -284,7 +280,7 @@ void NodePartitionedMeshReader::readProperties(
 
     const std::string fname_val = file_name_base + "_partitioned_" + item_type +
                                   "_properties_val" +
-                                  std::to_string(_mpi_comm_size) + ".bin";
+                                  std::to_string(mpi_.size) + ".bin";
     is.open(fname_val.c_str(), std::ios::binary | std::ios::in);
     if (!is)
     {
@@ -407,40 +403,20 @@ MeshLib::NodePartitionedMesh* NodePartitionedMeshReader::newMesh(
     std::vector<MeshLib::Element*> const& mesh_elems,
     MeshLib::Properties const& properties) const
 {
-    std::vector<std::size_t> gathered_n_regular_base_nodes(_mpi_comm_size);
-
-    MPI_Allgather(&_mesh_info.number_of_regular_base_nodes,
-                  1,
-                  MPI_UNSIGNED_LONG,
-                  gathered_n_regular_base_nodes.data(),
-                  1,
-                  MPI_UNSIGNED_LONG,
-                  _mpi_comm);
-
-    std::vector<std::size_t> n_regular_base_nodes_at_rank;
-    n_regular_base_nodes_at_rank.push_back(0);
-    std::partial_sum(begin(gathered_n_regular_base_nodes),
-                     end(gathered_n_regular_base_nodes),
-                     back_inserter(n_regular_base_nodes_at_rank));
-
-    std::vector<std::size_t> gathered_n_regular_high_order_nodes(
-        _mpi_comm_size);
+    std::vector<std::size_t> const gathered_n_regular_base_nodes =
+        BaseLib::MPI::allgather(_mesh_info.number_of_regular_base_nodes, mpi_);
+
+    std::vector<std::size_t> n_regular_base_nodes_at_rank =
+        BaseLib::sizesToOffsets(gathered_n_regular_base_nodes);
+
     std::size_t const n_regular_high_order_nodes =
         _mesh_info.number_of_regular_nodes -
         _mesh_info.number_of_regular_base_nodes;
-    MPI_Allgather(&n_regular_high_order_nodes,
-                  1,
-                  MPI_UNSIGNED_LONG,
-                  gathered_n_regular_high_order_nodes.data(),
-                  1,
-                  MPI_UNSIGNED_LONG,
-                  _mpi_comm);
-
-    std::vector<std::size_t> n_regular_high_order_nodes_at_rank;
-    n_regular_high_order_nodes_at_rank.push_back(0);
-    std::partial_sum(begin(gathered_n_regular_high_order_nodes),
-                     end(gathered_n_regular_high_order_nodes),
-                     back_inserter(n_regular_high_order_nodes_at_rank));
+    std::vector<std::size_t> const gathered_n_regular_high_order_nodes =
+        BaseLib::MPI::allgather(n_regular_high_order_nodes, mpi_);
+
+    std::vector<std::size_t> n_regular_high_order_nodes_at_rank =
+        BaseLib::sizesToOffsets(gathered_n_regular_high_order_nodes);
 
     return new MeshLib::NodePartitionedMesh(
         mesh_name, mesh_nodes, glb_node_ids, mesh_elems, properties,
diff --git a/MeshLib/IO/MPI_IO/NodePartitionedMeshReader.h b/MeshLib/IO/MPI_IO/NodePartitionedMeshReader.h
index 9874ed6bc3fe906ea2b39cce4aa15eee61df48ab..5d531b822f7202d1f84a579f1a0e543f2476a53d 100644
--- a/MeshLib/IO/MPI_IO/NodePartitionedMeshReader.h
+++ b/MeshLib/IO/MPI_IO/NodePartitionedMeshReader.h
@@ -19,6 +19,7 @@
 #include <string>
 #include <vector>
 
+#include "BaseLib/MPI.h"
 #include "MeshLib/IO/MPI_IO/PropertyVectorMetaData.h"
 #include "MeshLib/IO/NodeData.h"
 #include "MeshLib/NodePartitionedMesh.h"
@@ -53,14 +54,8 @@ public:
     MeshLib::NodePartitionedMesh* read(const std::string& file_name_base);
 
 private:
-    /// Pointer to MPI communicator.
-    MPI_Comm _mpi_comm;
-
-    /// Number of processes in the communicator: _mpi_comm.
-    int _mpi_comm_size;
-
-    /// Rank of compute core.
-    int _mpi_rank;
+    /// Pointer to MPI communicator, the rank and the size.
+    BaseLib::MPI::Mpi mpi_;
 
     /// MPI data type for struct NodeData.
     MPI_Datatype _mpi_node_type;
@@ -205,7 +200,7 @@ private:
             OGS_FATAL(
                 "Error in NodePartitionedMeshReader::readProperties: "
                 "Could not read part {:d} of the PropertyVector.",
-                _mpi_rank);
+                mpi_.rank);
     }
 
     /// Read data for property OGS_VERSION or IntegrationPointMetaData, and
@@ -220,12 +215,12 @@ private:
             pvmd.property_name, t, pvmd.number_of_components);
 
         std::size_t const property_vector_size =
-            pvmd.number_of_tuples / _mpi_comm_size;
+            pvmd.number_of_tuples / mpi_.size;
         pv->resize(property_vector_size);
 
         // Locate the start position of the data in the file for the current
         // rank.
-        is.seekg(global_offset + property_vector_size * sizeof(T) * _mpi_rank);
+        is.seekg(global_offset + property_vector_size * sizeof(T) * mpi_.rank);
 
         // read the values
         if (!is.read(reinterpret_cast<char*>(pv->data()),
@@ -234,7 +229,7 @@ private:
             OGS_FATAL(
                 "Error in NodePartitionedMeshReader::readProperties: "
                 "Could not read part {:d} of the PropertyVector.",
-                _mpi_rank);
+                mpi_.rank);
         }
     }
 
diff --git a/MeshLib/IO/VtkIO/VtuInterface.cpp b/MeshLib/IO/VtkIO/VtuInterface.cpp
index ee28445da54fb76f66c8fdfb0f19b0ba867ec50a..5c43b0dd61fad94a5d9e53adbd4395b3a649f178 100644
--- a/MeshLib/IO/VtkIO/VtuInterface.cpp
+++ b/MeshLib/IO/VtkIO/VtuInterface.cpp
@@ -26,13 +26,9 @@
 #include <boost/algorithm/string/erase.hpp>
 
 #include "BaseLib/DisableFPE.h"
-#include "BaseLib/Logging.h"
-
-#ifdef USE_PETSC
-#include <petsc.h>
-#endif
-
 #include "BaseLib/FileTools.h"
+#include "BaseLib/Logging.h"
+#include "BaseLib/MPI.h"
 #include "MeshLib/Mesh.h"
 #include "MeshLib/Vtk/VtkMappedMeshSource.h"
 #include "VtkMeshConverter.h"
@@ -148,23 +144,15 @@ std::string getVtuFileNameForPetscOutputWithoutExtension(
 bool VtuInterface::writeToFile(std::filesystem::path const& file_path)
 {
 #ifdef USE_PETSC
-    int mpi_init;
-    MPI_Initialized(&mpi_init);
-    if (mpi_init == 1)
+    BaseLib::MPI::Mpi mpi;
+    if (mpi.size == 1)
     {
-        int mpi_size;
-        MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
-        if (mpi_size == 1)
-        {
-            return writeVTU<vtkXMLUnstructuredGridWriter>(file_path.string());
-        }
-        auto const vtu_file_name =
-            getVtuFileNameForPetscOutputWithoutExtension(file_path.string());
-        int rank;
-        MPI_Comm_rank(MPI_COMM_WORLD, &rank);
-        return writeVTU<vtkXMLPUnstructuredGridWriter>(vtu_file_name + ".pvtu",
-                                                       mpi_size, rank);
+        return writeVTU<vtkXMLUnstructuredGridWriter>(file_path.string());
     }
+    auto const vtu_file_name =
+        getVtuFileNameForPetscOutputWithoutExtension(file_path.string());
+    return writeVTU<vtkXMLPUnstructuredGridWriter>(vtu_file_name + ".pvtu",
+                                                   mpi.size, mpi.rank);
 #endif
     return writeVTU<vtkXMLUnstructuredGridWriter>(file_path.string());
 }
diff --git a/MeshLib/IO/XDMF/mpi/partition.cpp b/MeshLib/IO/XDMF/mpi/partition.cpp
index c1855dc2ac0c6a06f9d109e8f3fd327de1fac163..19255358885db07410cf8b6144ece64ad6b7fbad 100644
--- a/MeshLib/IO/XDMF/mpi/partition.cpp
+++ b/MeshLib/IO/XDMF/mpi/partition.cpp
@@ -17,7 +17,9 @@
 
 #include <numeric>
 
+#include "BaseLib/Algorithm.h"
 #include "BaseLib/Logging.h"
+#include "BaseLib/MPI.h"
 #include "MeshLib/IO/XDMF/fileIO.h"
 #include "getCommunicator.h"
 
@@ -33,37 +35,22 @@ bool isFileManager()
 PartitionInfo getPartitionInfo(std::size_t const size,
                                unsigned int const n_files)
 {
-    MPI_Comm const mpi_comm = getCommunicator(n_files).mpi_communicator;
-    int mpi_size;
-    int mpi_rank;
-    MPI_Comm_size(mpi_comm, &mpi_size);
-    MPI_Comm_rank(mpi_comm, &mpi_rank);
-
-    std::vector<std::size_t> partition_sizes;
-    partition_sizes.resize(mpi_size);
+    BaseLib::MPI::Mpi const mpi{getCommunicator(n_files).mpi_communicator};
 
-    MPI_Allgather(&size,
-                  1,
-                  MPI_UNSIGNED_LONG,
-                  partition_sizes.data(),
-                  1,
-                  MPI_UNSIGNED_LONG,
-                  mpi_comm);
+    std::vector<std::size_t> const partition_sizes =
+        BaseLib::MPI::allgather(size, mpi);
 
     // the first partition's offset is zero, offsets for subsequent
-    // partitions are the accumulated sum of all preceding size (excluding
-    // own size)
-    std::vector<std::size_t> partition_offsets(1, 0);
-    std::partial_sum(partition_sizes.begin(),
-                     partition_sizes.end(),
-                     back_inserter(partition_offsets));
+    // partitions are the accumulated sum of all preceding size.
+    std::vector<std::size_t> const partition_offsets =
+        BaseLib::sizesToOffsets(partition_sizes);
 
     // chunked
     std::size_t longest_partition =
         *max_element(partition_sizes.begin(), partition_sizes.end());
 
     // local_offset, local_length, longest_local_length, global_length
-    return {partition_offsets[mpi_rank], size, longest_partition,
+    return {partition_offsets[mpi.rank], size, longest_partition,
             partition_offsets.back()};
 }
 }  // namespace MeshLib::IO
diff --git a/MeshLib/IO/readMeshFromFile.cpp b/MeshLib/IO/readMeshFromFile.cpp
index 3545a04b858e7d8715d3c86fb392d3912a729373..515ac9738ba1af28e10045ae7a094e1a59bfe0fd 100644
--- a/MeshLib/IO/readMeshFromFile.cpp
+++ b/MeshLib/IO/readMeshFromFile.cpp
@@ -21,6 +21,7 @@
 
 #include "BaseLib/FileTools.h"
 #include "BaseLib/Logging.h"
+#include "BaseLib/MPI.h"
 #include "BaseLib/StringTools.h"
 #include "MeshLib/IO/Legacy/MeshIO.h"
 #include "MeshLib/IO/VtkIO/VtuInterface.h"
@@ -68,33 +69,27 @@ MeshLib::Mesh* readMeshFromFile(const std::string& file_name,
                                 bool const compute_element_neighbors)
 {
 #ifdef USE_PETSC
-    int mpi_init;
-    MPI_Initialized(&mpi_init);
-    if (mpi_init == 1)
+    BaseLib::MPI::Mpi mpi;
+    if (mpi.size > 1)
     {
-        int world_size;
-        MPI_Comm_size(MPI_COMM_WORLD, &world_size);
-        if (world_size > 1)
+        MeshLib::IO::NodePartitionedMeshReader read_pmesh(mpi.communicator);
+        const std::string file_name_base =
+            BaseLib::dropFileExtension(file_name);
+        return read_pmesh.read(file_name_base);
+    }
+    if (mpi.size == 1)
+    {
+        std::unique_ptr<Mesh> mesh{
+            readMeshFromFileSerial(file_name, compute_element_neighbors)};
+
+        if (!mesh)
         {
-            MeshLib::IO::NodePartitionedMeshReader read_pmesh(MPI_COMM_WORLD);
-            const std::string file_name_base =
-                BaseLib::dropFileExtension(file_name);
-            return read_pmesh.read(file_name_base);
+            return nullptr;
         }
-        if (world_size == 1)
-        {
-            std::unique_ptr<Mesh> mesh{
-                readMeshFromFileSerial(file_name, compute_element_neighbors)};
 
-            if (!mesh)
-            {
-                return nullptr;
-            }
-
-            return new MeshLib::NodePartitionedMesh(*mesh);
-        }
-        return nullptr;
+        return new MeshLib::NodePartitionedMesh(*mesh);
     }
+    return nullptr;
 #endif
     return readMeshFromFileSerial(file_name, compute_element_neighbors);
 }
diff --git a/MeshLib/Utils/transformMeshToNodePartitionedMesh.cpp b/MeshLib/Utils/transformMeshToNodePartitionedMesh.cpp
index 0c627e33cb4a112ebe458477e4e7ff6a37a00e45..4a8fed12f10f9d957e8758cac2006d985c2b612d 100644
--- a/MeshLib/Utils/transformMeshToNodePartitionedMesh.cpp
+++ b/MeshLib/Utils/transformMeshToNodePartitionedMesh.cpp
@@ -22,6 +22,7 @@
 #include <vector>
 
 #include "BaseLib/Logging.h"
+#include "BaseLib/MPI.h"
 #include "MeshLib/Elements/Element.h"
 #include "MeshLib/Mesh.h"
 #include "MeshLib/Node.h"
@@ -42,15 +43,6 @@ bool isRegularNode(
 
 namespace MeshLib
 {
-std::pair<int, int> getMPIRankAndSize(MPI_Comm const& mpi_comm)
-{
-    int mpi_comm_size;
-    MPI_Comm_size(mpi_comm, &mpi_comm_size);
-    int mpi_comm_rank;
-    MPI_Comm_rank(mpi_comm, &mpi_comm_rank);
-    return {mpi_comm_rank, mpi_comm_size};
-}
-
 std::pair<std::vector<Node*>, std::vector<Element*>> copyNodesAndElements(
     std::vector<Element*> const& input_elements)
 {
@@ -136,27 +128,20 @@ computeRegularBaseNodeGlobalNodeIDsOfSubDomainPartition(
     std::size_t const number_of_regular_nodes =
         computeNumberOfRegularNodes(bulk_mesh, subdomain_mesh);
 
-    // in the following information exchange with other ranks is required
-    MPI_Comm mpi_comm = MPI_COMM_WORLD;
-    auto const [mpi_comm_rank, mpi_comm_size] = getMPIRankAndSize(mpi_comm);
+    BaseLib::MPI::Mpi mpi;
 
     // send own number of regular nodes to all others
-    std::vector<std::size_t> gathered_number_of_regular_nodes(mpi_comm_size);
-    MPI_Allgather(&number_of_regular_nodes, 1, MPI_UNSIGNED_LONG,
-                  gathered_number_of_regular_nodes.data(), 1, MPI_UNSIGNED_LONG,
-                  mpi_comm);
+    std::vector<std::size_t> const gathered_number_of_regular_nodes =
+        BaseLib::MPI::allgather(number_of_regular_nodes, mpi);
+
     // compute the 'offset' in the global_node_ids
-    std::vector<std::size_t> numbers_of_regular_nodes_at_rank;
-    numbers_of_regular_nodes_at_rank.push_back(0);
-    std::partial_sum(begin(gathered_number_of_regular_nodes),
-                     end(gathered_number_of_regular_nodes),
-                     back_inserter(numbers_of_regular_nodes_at_rank));
+    std::vector<std::size_t> const numbers_of_regular_nodes_at_rank =
+        BaseLib::sizesToOffsets(gathered_number_of_regular_nodes);
 
     // add the offset to the partitioned-owned subdomain
     std::vector<std::size_t> subdomain_global_node_ids;
     subdomain_global_node_ids.reserve(subdomain_nodes.size());
-    auto const partition_offset =
-        numbers_of_regular_nodes_at_rank[mpi_comm_rank];
+    auto const partition_offset = numbers_of_regular_nodes_at_rank[mpi.rank];
     DBUG("[{}] partition offset: {}", subdomain_mesh->getName(),
          partition_offset);
     // set the global id for the regular base nodes
@@ -175,9 +160,7 @@ std::vector<std::size_t> computeGhostBaseNodeGlobalNodeIDsOfSubDomainPartition(
     Mesh const* subdomain_mesh,
     std::vector<std::size_t> const& global_regular_base_node_ids)
 {
-    // in the following information exchange with other ranks is required
-    MPI_Comm const mpi_comm = MPI_COMM_WORLD;
-    auto const [mpi_comm_rank, mpi_comm_size] = getMPIRankAndSize(mpi_comm);
+    BaseLib::MPI::Mpi mpi;
 
     // count regular nodes that is the offset in the mapping
     auto const local_bulk_node_ids_for_subdomain =
@@ -197,36 +180,19 @@ std::vector<std::size_t> computeGhostBaseNodeGlobalNodeIDsOfSubDomainPartition(
         }
     }
 
-    // send ids of bulk ghost nodes belonging to the subdomain mesh to all
-    // other ranks and at the same time receive from all other ranks
-    // first send the sizes to all other to are able to allocate buffer
-    auto const size = subdomain_node_id_to_bulk_node_id.size();
-    std::size_t global_number_of_subdomain_node_id_to_bulk_node_id = 0;
-    MPI_Allreduce(&size, &global_number_of_subdomain_node_id_to_bulk_node_id, 1,
-                  MPI_UNSIGNED_LONG, MPI_SUM, mpi_comm);
+    // Send ids of bulk ghost nodes belonging to the subdomain mesh to all
+    // other ranks and at the same time receive from all other ranks.
+    std::vector<std::size_t> ghost_node_ids_of_all_ranks;
+    std::vector<int> const offsets =
+        BaseLib::MPI::allgatherv(std::span{subdomain_node_id_to_bulk_node_id},
+                                 ghost_node_ids_of_all_ranks, mpi);
 
+    std::size_t const global_number_of_subdomain_node_id_to_bulk_node_id =
+        offsets.back();
     DBUG("[{}] global_number_of_subdomain_node_id_to_bulk_node_id: '{}' ",
          subdomain_mesh->getName(),
          global_number_of_subdomain_node_id_to_bulk_node_id);
 
-    std::vector<int> numbers_of_ids_at_ranks(mpi_comm_size);
-    MPI_Allgather(&size, 1, MPI_INT, numbers_of_ids_at_ranks.data(), 1, MPI_INT,
-                  mpi_comm);
-    std::vector<int> offsets;
-    offsets.push_back(0);
-    std::partial_sum(begin(numbers_of_ids_at_ranks),
-                     end(numbers_of_ids_at_ranks), back_inserter(offsets));
-    std::vector<std::size_t> ghost_node_ids_of_all_ranks(
-        global_number_of_subdomain_node_id_to_bulk_node_id);
-    MPI_Allgatherv(subdomain_node_id_to_bulk_node_id.data(), /* sendbuf */
-                   size,                                     /* sendcount */
-                   MPI_UNSIGNED_LONG,                        /* sendtype */
-                   ghost_node_ids_of_all_ranks.data(),       /* recvbuf (out) */
-                   numbers_of_ids_at_ranks.data(),           /* recvcounts */
-                   offsets.data(),                           /* displs */
-                   MPI_UNSIGNED_LONG,                        /* recvtype */
-                   mpi_comm);
-
     // construct a map for fast search of local bulk node ids
     std::map<std::size_t, std::size_t> global_to_local_bulk_node_ids;
     for (auto const id : bulk_mesh->getNodes() | MeshLib::views::ids)
@@ -244,13 +210,13 @@ std::vector<std::size_t> computeGhostBaseNodeGlobalNodeIDsOfSubDomainPartition(
         local_subdomain_node_ids[local_bulk_node_ids_for_subdomain[id]] = id;
     }
 
-    std::vector<std::size_t> local_subdomain_node_ids_of_all_ranks(
+    std::vector<std::size_t> subdomain_node_ids_of_all_ranks(
         global_number_of_subdomain_node_id_to_bulk_node_id,
         std::numeric_limits<std::size_t>::max());
     // search in all ranks within the bulk ids for the corresponding id
-    for (int rank = 0; rank < mpi_comm_size; ++rank)
+    for (int rank = 0; rank < mpi.size; ++rank)
     {
-        if (rank == mpi_comm_rank)
+        if (rank == mpi.rank)
         {
             continue;
         }
@@ -263,35 +229,24 @@ std::vector<std::size_t> computeGhostBaseNodeGlobalNodeIDsOfSubDomainPartition(
                 continue;
             }
             auto const local_bulk_node_id = it->second;
-            local_subdomain_node_ids_of_all_ranks[i] =
-                global_regular_base_node_ids
-                    [local_subdomain_node_ids.find(local_bulk_node_id)->second];
+            subdomain_node_ids_of_all_ranks[i] = global_regular_base_node_ids
+                [local_subdomain_node_ids.find(local_bulk_node_id)->second];
             DBUG(
                 "[{}] found global subdomain node id: '{}' for global bulk "
                 "node id {} ",
                 subdomain_mesh->getName(),
-                local_subdomain_node_ids_of_all_ranks[i],
+                subdomain_node_ids_of_all_ranks[i],
                 ghost_node_ids_of_all_ranks[i]);
         }
     }
 
-    // send the computed ids back
-    std::vector<std::size_t> computed_global_ids_for_subdomain_ghost_nodes(
-        global_number_of_subdomain_node_id_to_bulk_node_id);
-    MPI_Allreduce(
-        local_subdomain_node_ids_of_all_ranks.data(), /* sendbuf */
-        computed_global_ids_for_subdomain_ghost_nodes
-            .data(),                                        /* recvbuf (out) */
-        global_number_of_subdomain_node_id_to_bulk_node_id, /* sendcount */
-        MPI_UNSIGNED_LONG,                                  /* sendtype */
-        MPI_MAX,                                            /* operation */
-        mpi_comm);
+    // find maximum over all ranks
+    BaseLib::MPI::allreduceInplace(subdomain_node_ids_of_all_ranks, MPI_MAX,
+                                   mpi);
 
     std::vector<std::size_t> global_ids_for_subdomain_ghost_nodes(
-        computed_global_ids_for_subdomain_ghost_nodes.begin() +
-            offsets[mpi_comm_rank],
-        computed_global_ids_for_subdomain_ghost_nodes.begin() +
-            offsets[mpi_comm_rank + 1]);
+        subdomain_node_ids_of_all_ranks.begin() + offsets[mpi.rank],
+        subdomain_node_ids_of_all_ranks.begin() + offsets[mpi.rank + 1]);
     return global_ids_for_subdomain_ghost_nodes;
 }
 
@@ -301,51 +256,30 @@ std::vector<std::size_t> computeNumberOfRegularBaseNodesAtRank(
     auto const number_of_regular_base_nodes =
         subdomain_mesh->computeNumberOfBaseNodes();
 
-    MPI_Comm const mpi_comm = MPI_COMM_WORLD;
-    auto const [mpi_comm_rank, mpi_comm_size] = getMPIRankAndSize(mpi_comm);
+    BaseLib::MPI::Mpi mpi;
 
-    std::vector<std::size_t> gathered_number_of_regular_base_nodes(
-        mpi_comm_size);
-    MPI_Allgather(&number_of_regular_base_nodes, 1, MPI_UNSIGNED_LONG,
-                  gathered_number_of_regular_base_nodes.data(), 1,
-                  MPI_UNSIGNED_LONG, mpi_comm);
+    std::vector<std::size_t> const gathered_number_of_regular_base_nodes =
+        BaseLib::MPI::allgather(number_of_regular_base_nodes, mpi);
 
-    std::vector<std::size_t> numbers_of_regular_base_nodes_at_rank;
-    numbers_of_regular_base_nodes_at_rank.push_back(0);
-    std::partial_sum(begin(gathered_number_of_regular_base_nodes),
-                     end(gathered_number_of_regular_base_nodes),
-                     back_inserter(numbers_of_regular_base_nodes_at_rank));
-
-    return numbers_of_regular_base_nodes_at_rank;
+    return BaseLib::sizesToOffsets(gathered_number_of_regular_base_nodes);
 }
 
 // similar to the above only with regular higher order nodes
 std::vector<std::size_t> computeNumberOfRegularHigherOrderNodesAtRank(
     Mesh const* subdomain_mesh)
 {
-    // in the following information exchange with other ranks is required
-    MPI_Comm const mpi_comm = MPI_COMM_WORLD;
-    auto [mpi_comm_rank, mpi_comm_size] = getMPIRankAndSize(mpi_comm);
+    BaseLib::MPI::Mpi mpi;
 
     auto const number_of_regular_base_nodes =
         subdomain_mesh->computeNumberOfBaseNodes();
 
-    std::vector<std::size_t> gathered_number_of_regular_higher_order_nodes(
-        mpi_comm_size);
     auto const number_of_regular_higher_order_nodes =
         subdomain_mesh->getNumberOfNodes() - number_of_regular_base_nodes;
-    MPI_Allgather(&number_of_regular_higher_order_nodes, 1, MPI_UNSIGNED_LONG,
-                  gathered_number_of_regular_higher_order_nodes.data(), 1,
-                  MPI_UNSIGNED_LONG, mpi_comm);
-
-    std::vector<std::size_t> numbers_of_regular_higher_order_nodes_at_rank;
-    numbers_of_regular_higher_order_nodes_at_rank.push_back(0);
-    std::partial_sum(
-        begin(gathered_number_of_regular_higher_order_nodes),
-        end(gathered_number_of_regular_higher_order_nodes),
-        back_inserter(numbers_of_regular_higher_order_nodes_at_rank));
-
-    return numbers_of_regular_higher_order_nodes_at_rank;
+    std::vector<std::size_t> gathered_number_of_regular_higher_order_nodes =
+        BaseLib::MPI::allgather(number_of_regular_higher_order_nodes, mpi);
+
+    return BaseLib::sizesToOffsets(
+        gathered_number_of_regular_higher_order_nodes);
 }
 
 std::vector<std::size_t> computeGlobalNodeIDsOfSubDomainPartition(
@@ -382,13 +316,8 @@ std::vector<std::size_t> computeGlobalNodeIDsOfSubDomainPartition(
 unsigned long getNumberOfGlobalNodes(Mesh const* subdomain_mesh)
 {
     // sum all nodes over all partitions in number_of_global_nodes
-    unsigned long number_of_local_nodes = subdomain_mesh->getNodes().size();
-    unsigned long number_of_global_nodes = 0;
-
-    MPI_Comm mpi_comm = MPI_COMM_WORLD;
-
-    MPI_Allreduce(&number_of_local_nodes, &number_of_global_nodes, 1,
-                  MPI_UNSIGNED_LONG, MPI_SUM, mpi_comm);
+    unsigned long const number_of_global_nodes = BaseLib::MPI::allreduce(
+        subdomain_mesh->getNodes().size(), MPI_SUM, BaseLib::MPI::Mpi{});
     DBUG("[{}] number_of_global_nodes: {}'", subdomain_mesh->getName(),
          number_of_global_nodes);
     return number_of_global_nodes;
diff --git a/NumLib/DOF/DOFTableUtil.cpp b/NumLib/DOF/DOFTableUtil.cpp
index 93266fc65862929b316109aa30466fea0b2c0d24..4241919a6ecc5a3281bd9283ddb353f271096bdf 100644
--- a/NumLib/DOF/DOFTableUtil.cpp
+++ b/NumLib/DOF/DOFTableUtil.cpp
@@ -14,6 +14,7 @@
 #include <cassert>
 #include <functional>
 
+#include "BaseLib/MPI.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 namespace NumLib
 {
@@ -44,15 +45,12 @@ double norm(GlobalVector const& x, unsigned const global_component,
 double norm1(GlobalVector const& x, unsigned const global_component,
              LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
 {
-    double res =
+    double const res =
         norm(x, global_component, dof_table, mesh,
              [](double res, double value) { return res + std::abs(value); });
 
 #ifdef USE_PETSC
-    double global_result = 0.0;
-    MPI_Allreduce(&res, &global_result, 1, MPI_DOUBLE, MPI_SUM,
-                  PETSC_COMM_WORLD);
-    res = global_result;
+    return BaseLib::MPI::allreduce(res, MPI_SUM, BaseLib::MPI::Mpi{});
 #endif
     return res;
 }
@@ -60,15 +58,13 @@ double norm1(GlobalVector const& x, unsigned const global_component,
 double norm2(GlobalVector const& x, unsigned const global_component,
              LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
 {
-    double res =
+    double const res =
         norm(x, global_component, dof_table, mesh,
              [](double res, double value) { return res + value * value; });
 
 #ifdef USE_PETSC
-    double global_result = 0.0;
-    MPI_Allreduce(&res, &global_result, 1, MPI_DOUBLE, MPI_SUM,
-                  PETSC_COMM_WORLD);
-    res = global_result;
+    return std::sqrt(
+        BaseLib::MPI::allreduce(res, MPI_SUM, BaseLib::MPI::Mpi{}));
 #endif
     return std::sqrt(res);
 }
@@ -77,15 +73,12 @@ double normInfinity(GlobalVector const& x, unsigned const global_component,
                     LocalToGlobalIndexMap const& dof_table,
                     MeshLib::Mesh const& mesh)
 {
-    double res = norm(x, global_component, dof_table, mesh,
-                      [](double res, double value)
-                      { return std::max(res, std::abs(value)); });
+    double const res =
+        norm(x, global_component, dof_table, mesh, [](double res, double value)
+             { return std::max(res, std::abs(value)); });
 
 #ifdef USE_PETSC
-    double global_result = 0.0;
-    MPI_Allreduce(&res, &global_result, 1, MPI_DOUBLE, MPI_MAX,
-                  PETSC_COMM_WORLD);
-    res = global_result;
+    return BaseLib::MPI::allreduce(res, MPI_MAX, BaseLib::MPI::Mpi{});
 #endif
     return res;
 }
diff --git a/NumLib/DOF/LocalToGlobalIndexMap.cpp b/NumLib/DOF/LocalToGlobalIndexMap.cpp
index 1d37e1a3856f5c26d9ec8a2ea1ae4c4e7948010c..4ac7cea9902d6dc0962073cc9fa8b27a4045f7c3 100644
--- a/NumLib/DOF/LocalToGlobalIndexMap.cpp
+++ b/NumLib/DOF/LocalToGlobalIndexMap.cpp
@@ -16,19 +16,6 @@
 
 namespace NumLib
 {
-namespace
-{
-// Make the cumulative sum of an array, which starts with zero
-template <typename T>
-std::vector<T> to_cumulative(std::vector<T> const& vec)
-{
-    std::vector<T> result(vec.size() + 1, 0);
-    std::partial_sum(vec.begin(), vec.end(), result.begin() + 1);
-
-    return result;
-}
-
-}  // namespace
 
 int LocalToGlobalIndexMap::getGlobalComponent(int const variable_id,
                                               int const component_id) const
@@ -129,7 +116,7 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
     NumLib::ComponentOrder const order)
     : _mesh_subsets(std::move(mesh_subsets)),
       _mesh_component_map(_mesh_subsets, order),
-      _variable_component_offsets(to_cumulative(vec_var_n_components))
+      _variable_component_offsets(BaseLib::sizesToOffsets(vec_var_n_components))
 {
     // For each element of that MeshSubset save a line of global indices.
     for (int variable_id = 0;
@@ -160,7 +147,7 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
     NumLib::ComponentOrder const order)
     : _mesh_subsets(std::move(mesh_subsets)),
       _mesh_component_map(_mesh_subsets, order),
-      _variable_component_offsets(to_cumulative(vec_var_n_components))
+      _variable_component_offsets(BaseLib::sizesToOffsets(vec_var_n_components))
 {
     assert(vec_var_n_components.size() == vec_var_elements.size());
 
diff --git a/ProcessLib/PhaseField/PhaseFieldProcess.cpp b/ProcessLib/PhaseField/PhaseFieldProcess.cpp
index 3020a93570299bd68a0b1945c35cebdef71f1edd..8c150546174e887cfc6594fdf869492efd8271fd 100644
--- a/ProcessLib/PhaseField/PhaseFieldProcess.cpp
+++ b/ProcessLib/PhaseField/PhaseFieldProcess.cpp
@@ -12,6 +12,7 @@
 
 #include <cassert>
 
+#include "BaseLib/MPI.h"
 #include "MeshLib/Utils/getOrCreateMeshProperty.h"
 #include "NumLib/DOF/ComputeSparsityPattern.h"
 #include "PhaseFieldFEM.h"
@@ -307,15 +308,13 @@ void PhaseFieldProcess<DisplacementDim>::postTimestepConcreteProcess(
             _process_data.pressure_work);
 
 #ifdef USE_PETSC
-        double const elastic_energy = _process_data.elastic_energy;
-        MPI_Allreduce(&elastic_energy, &_process_data.elastic_energy, 1,
-                      MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
-        double const surface_energy = _process_data.surface_energy;
-        MPI_Allreduce(&surface_energy, &_process_data.surface_energy, 1,
-                      MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
-        double const pressure_work = _process_data.pressure_work;
-        MPI_Allreduce(&pressure_work, &_process_data.pressure_work, 1,
-                      MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
+        BaseLib::MPI::Mpi mpi{};
+        _process_data.elastic_energy =
+            BaseLib::MPI::allreduce(_process_data.elastic_energy, MPI_SUM, mpi);
+        _process_data.surface_energy =
+            BaseLib::MPI::allreduce(_process_data.surface_energy, MPI_SUM, mpi);
+        _process_data.pressure_work =
+            BaseLib::MPI::allreduce(_process_data.pressure_work, MPI_SUM, mpi);
 #endif
 
         INFO(
@@ -361,9 +360,8 @@ void PhaseFieldProcess<DisplacementDim>::postNonLinearSolverConcreteProcess(
         getActiveElementIDs(), dof_tables, x, t, _process_data.crack_volume);
 
 #ifdef USE_PETSC
-    double const crack_volume = _process_data.crack_volume;
-    MPI_Allreduce(&crack_volume, &_process_data.crack_volume, 1, MPI_DOUBLE,
-                  MPI_SUM, PETSC_COMM_WORLD);
+    _process_data.crack_volume = BaseLib::MPI::allreduce(
+        _process_data.crack_volume, MPI_SUM, BaseLib::MPI::Mpi{});
 #endif
 
     INFO("Integral of crack: {:g}", _process_data.crack_volume);
diff --git a/Tests/BaseLib/MPI.cpp b/Tests/BaseLib/MPI.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..16aadb59b89a3dbde9dafe38a953c5e8d582e1f0
--- /dev/null
+++ b/Tests/BaseLib/MPI.cpp
@@ -0,0 +1,113 @@
+
+#include "BaseLib/MPI.h"
+
+#include <gtest/gtest.h>
+
+using namespace BaseLib::MPI;
+
+#ifdef USE_PETSC
+// Fixture for MPI tests
+struct MPI_BaseLib : public ::testing::Test
+{
+    Mpi mpi{};
+};
+
+TEST_F(MPI_BaseLib, CommunicatorRankAndSize)
+{
+    EXPECT_EQ(mpi.size, 3);
+    EXPECT_GE(mpi.rank, 0);
+    EXPECT_LT(mpi.rank, 3);
+}
+
+TEST_F(MPI_BaseLib, Allgather)
+{
+    std::vector<int> const gathered_values = allgather(mpi.rank, mpi);
+
+    std::vector<int> const expected_values = {0, 1, 2};
+    // Each rank's value should match its rank
+    EXPECT_EQ(expected_values, gathered_values);
+}
+
+TEST_F(MPI_BaseLib, AllgatherVector)
+{
+    std::vector<int> const send_values(2, mpi.rank);
+    std::vector<int> const gathered_values = allgather(send_values, mpi);
+
+    std::vector<int> const expected_values = {0, 0, 1, 1, 2, 2};
+    EXPECT_EQ(expected_values, gathered_values);
+}
+
+TEST_F(MPI_BaseLib, Allgatherv)
+{
+    int const rank_value =
+        mpi.rank + 1;  // Each rank contributes rank+1 elements
+    std::vector<int> const send_data(
+        rank_value, mpi.rank);  // Each element is set to the rank number
+
+    std::vector<int> gathered_values(1 + 2 + 3);
+
+    std::vector<int> const offsets =
+        allgatherv(std::span(send_data), gathered_values, mpi);
+    std::vector<int> const expected_offsets = {0, 1, 3, 6};  // last element is
+                                                             // size.
+    std::vector<int> const expected_gathered_values = {0, 1, 1, 2, 2, 2};
+    EXPECT_EQ(expected_offsets, offsets);
+    EXPECT_EQ(expected_gathered_values, gathered_values);
+}
+
+TEST_F(MPI_BaseLib, Allreduce)
+{
+    int const sum = allreduce(mpi.rank + 1, MPI_SUM, mpi);
+    int const expected_sum = 1 + 2 + 3;
+    EXPECT_EQ(expected_sum, sum);
+
+    int const product = allreduce(mpi.rank + 1, MPI_PROD, mpi);
+    int const expected_product = 1 * 2 * 3;
+    EXPECT_EQ(expected_product, product);
+
+    int const min = allreduce(mpi.rank + 1, MPI_MIN, mpi);
+    int const expected_min = 1;
+    EXPECT_EQ(expected_min, min);
+}
+
+TEST_F(MPI_BaseLib, AllreduceVector)
+{
+    std::vector<int> const send_values = {mpi.rank + 1,
+                                          mpi.size + mpi.rank + 1};
+    std::vector<int> const sum = allreduce(send_values, MPI_SUM, mpi);
+    std::vector<int> const expected_sum = {(1 + 2 + 3), (4 + 5 + 6)};
+    EXPECT_EQ(expected_sum, sum);
+
+    std::vector<int> const product = allreduce(send_values, MPI_PROD, mpi);
+    std::vector<int> const expected_product = {(1 * 2 * 3), (4 * 5 * 6)};
+    EXPECT_EQ(expected_product, product);
+
+    std::vector<int> const min = allreduce(send_values, MPI_MIN, mpi);
+    std::vector<int> const expected_min = {1, 4};
+    EXPECT_EQ(expected_min, min);
+}
+
+TEST_F(MPI_BaseLib, AllreduceVectorInplace)
+{
+    {
+        std::vector<int> values = {mpi.rank + 1, mpi.size + mpi.rank + 1};
+        allreduceInplace(values, MPI_SUM, mpi);
+        std::vector<int> const expected_sum = {(1 + 2 + 3), (4 + 5 + 6)};
+        EXPECT_EQ(expected_sum, values);
+    }
+
+    {
+        std::vector<int> values = {mpi.rank + 1, mpi.size + mpi.rank + 1};
+        allreduceInplace(values, MPI_PROD, mpi);
+        std::vector<int> const expected_product = {(1 * 2 * 3), (4 * 5 * 6)};
+        EXPECT_EQ(expected_product, values);
+    }
+
+    {
+        std::vector<int> values = {mpi.rank + 1, mpi.size + mpi.rank + 1};
+        allreduceInplace(values, MPI_MIN, mpi);
+        std::vector<int> const expected_min = {1, 4};
+        EXPECT_EQ(expected_min, values);
+    }
+}
+#endif
diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt
index 1a450db81a8d57b149b15b08612155b9897a56ea..b92d80da9526040a9ba3575374446b0ca8209dba 100644
--- a/Tests/CMakeLists.txt
+++ b/Tests/CMakeLists.txt
@@ -190,10 +190,10 @@ if(OGS_USE_PETSC)
     if("${HOSTNAME}" MATCHES "frontend.*")
         list(APPEND MPIRUN_ARGS --mca btl_openib_allow_ib 1)
     endif()
-    set(TEST_FILTER_MPI --gtest_filter=-MPITest*)
+    set(TEST_FILTER_MPI --gtest_filter=-MPI*)
     add_custom_target(tests
         mpirun ${MPIRUN_ARGS} -np 1 $<TARGET_FILE:testrunner> ${TESTRUNNER_ADDITIONAL_ARGUMENTS} ${TEST_FILTER_MPI}
-        COMMAND mpirun ${MPIRUN_ARGS} -np 3 $<TARGET_FILE:testrunner> --gtest_filter=MPITest*
+        COMMAND mpirun ${MPIRUN_ARGS} -np 3 $<TARGET_FILE:testrunner> --gtest_filter=MPI*
         DEPENDS testrunner tests-cleanup
     )
 else()
diff --git a/Tests/MathLib/TestGlobalMatrixInterface.cpp b/Tests/MathLib/TestGlobalMatrixInterface.cpp
index 075084ac1efe9a95b80e1f929968d379e0d9f43c..46f3f29cebf124906da12bc23562d499c8f10e05 100644
--- a/Tests/MathLib/TestGlobalMatrixInterface.cpp
+++ b/Tests/MathLib/TestGlobalMatrixInterface.cpp
@@ -17,6 +17,7 @@
 
 #include <Eigen/Core>
 
+#include "BaseLib/MPI.h"
 #include "MathLib/LinAlg/LinAlg.h"
 
 #if defined(USE_PETSC)
@@ -58,24 +59,17 @@ void checkGlobalMatrixInterface(T_MATRIX& m)
 template <class T_MATRIX, class T_VECTOR>
 void checkGlobalMatrixInterfaceMPI(T_MATRIX& m, T_VECTOR& v)
 {
-    int msize;
-    MPI_Comm_size(PETSC_COMM_WORLD, &msize);
-    int mrank;
-    MPI_Comm_rank(PETSC_COMM_WORLD, &mrank);
+    BaseLib::MPI::Mpi mpi{PETSC_COMM_WORLD};
 
-    ASSERT_EQ(3u, msize);
+    ASSERT_EQ(3u, mpi.size);
     ASSERT_EQ(m.getRangeEnd() - m.getRangeBegin(), m.getNumberOfLocalRows());
 
-    int gathered_rows;
-    int local_rows = m.getNumberOfLocalRows();
-    MPI_Allreduce(&local_rows, &gathered_rows, 1, MPI_INT, MPI_SUM,
-                  PETSC_COMM_WORLD);
+    int const gathered_rows =
+        BaseLib::MPI::allreduce(m.getNumberOfLocalRows(), MPI_SUM, mpi);
     ASSERT_EQ(m.getNumberOfRows(), gathered_rows);
 
-    int gathered_cols;
-    int local_cols = m.getNumberOfLocalColumns();
-    MPI_Allreduce(&local_cols, &gathered_cols, 1, MPI_INT, MPI_SUM,
-                  PETSC_COMM_WORLD);
+    int const gathered_cols =
+        BaseLib::MPI::allreduce(m.getNumberOfLocalColumns(), MPI_SUM, mpi);
     ASSERT_EQ(m.getNumberOfColumns(), gathered_cols);
 
     // Add entries
@@ -87,8 +81,8 @@ void checkGlobalMatrixInterfaceMPI(T_MATRIX& m, T_VECTOR& v)
 
     std::vector<GlobalIndexType> row_pos(2);
     std::vector<GlobalIndexType> col_pos(2);
-    row_pos[0] = 2 * mrank;
-    row_pos[1] = 2 * mrank + 1;
+    row_pos[0] = 2 * mpi.rank;
+    row_pos[1] = 2 * mpi.rank + 1;
     col_pos[0] = row_pos[0];
     col_pos[1] = row_pos[1];
 
@@ -112,10 +106,10 @@ void checkGlobalMatrixInterfaceMPI(T_MATRIX& m, T_VECTOR& v)
     ASSERT_EQ(sqrt(3 * (3 * 3 + 7 * 7)), norm2(y));
 
     // set a value
-    m_c.set(2 * mrank, 2 * mrank, 5.0);
+    m_c.set(2 * mpi.rank, 2 * mpi.rank, 5.0);
     MathLib::finalizeMatrixAssembly(m_c);
     // add a value
-    m_c.add(2 * mrank + 1, 2 * mrank + 1, 5.0);
+    m_c.add(2 * mpi.rank + 1, 2 * mpi.rank + 1, 5.0);
     MathLib::finalizeMatrixAssembly(m_c);
 
     matMult(m_c, v, y);
@@ -127,21 +121,16 @@ void checkGlobalMatrixInterfaceMPI(T_MATRIX& m, T_VECTOR& v)
 template <class T_MATRIX, class T_VECTOR>
 void checkGlobalRectangularMatrixInterfaceMPI(T_MATRIX& m, T_VECTOR& v)
 {
-    int mrank;
-    MPI_Comm_rank(PETSC_COMM_WORLD, &mrank);
+    BaseLib::MPI::Mpi mpi{PETSC_COMM_WORLD};
 
     ASSERT_EQ(m.getRangeEnd() - m.getRangeBegin(), m.getNumberOfLocalRows());
 
-    int gathered_rows;
-    int local_rows = m.getNumberOfLocalRows();
-    MPI_Allreduce(&local_rows, &gathered_rows, 1, MPI_INT, MPI_SUM,
-                  PETSC_COMM_WORLD);
+    int const gathered_rows =
+        BaseLib::MPI::allreduce(m.getNumberOfLocalRows(), MPI_SUM, mpi);
     ASSERT_EQ(m.getNumberOfRows(), gathered_rows);
 
-    int gathered_cols;
-    int local_cols = m.getNumberOfLocalColumns();
-    MPI_Allreduce(&local_cols, &gathered_cols, 1, MPI_INT, MPI_SUM,
-                  PETSC_COMM_WORLD);
+    int const gathered_cols =
+        BaseLib::MPI::allreduce(m.getNumberOfLocalColumns(), MPI_SUM, mpi);
     ASSERT_EQ(m.getNumberOfColumns(), gathered_cols);
 
     // Add entries
@@ -155,11 +144,11 @@ void checkGlobalRectangularMatrixInterfaceMPI(T_MATRIX& m, T_VECTOR& v)
 
     std::vector<GlobalIndexType> row_pos(2);
     std::vector<GlobalIndexType> col_pos(3);
-    row_pos[0] = 2 * mrank;
-    row_pos[1] = 2 * mrank + 1;
-    col_pos[0] = 3 * mrank;
-    col_pos[1] = 3 * mrank + 1;
-    col_pos[2] = 3 * mrank + 2;
+    row_pos[0] = 2 * mpi.rank;
+    row_pos[1] = 2 * mpi.rank + 1;
+    col_pos[0] = 3 * mpi.rank;
+    col_pos[1] = 3 * mpi.rank + 1;
+    col_pos[2] = 3 * mpi.rank + 2;
 
     m.add(row_pos, col_pos, loc_m);
 
@@ -178,7 +167,7 @@ void checkGlobalRectangularMatrixInterfaceMPI(T_MATRIX& m, T_VECTOR& v)
 }  // end namespace
 
 #if defined(USE_PETSC)
-TEST(MPITest_Math, CheckInterface_PETScMatrix_Local_Size)
+TEST(MPI_Math, CheckInterface_PETScMatrix_Local_Size)
 {
     MathLib::PETScMatrixOption opt;
     opt.d_nz = 2;
@@ -193,7 +182,7 @@ TEST(MPITest_Math, CheckInterface_PETScMatrix_Local_Size)
     checkGlobalMatrixInterfaceMPI(A, x);
 }
 
-TEST(MPITest_Math, CheckInterface_PETScMatrix_Global_Size)
+TEST(MPI_Math, CheckInterface_PETScMatrix_Global_Size)
 {
     MathLib::PETScMatrixOption opt;
     opt.d_nz = 2;
@@ -205,7 +194,7 @@ TEST(MPITest_Math, CheckInterface_PETScMatrix_Global_Size)
     checkGlobalMatrixInterfaceMPI(A, x);
 }
 
-TEST(MPITest_Math, CheckInterface_PETSc_Rectangular_Matrix_Local_Size)
+TEST(MPI_Math, CheckInterface_PETSc_Rectangular_Matrix_Local_Size)
 {
     MathLib::PETScMatrixOption opt;
     opt.d_nz = 3;
@@ -220,7 +209,7 @@ TEST(MPITest_Math, CheckInterface_PETSc_Rectangular_Matrix_Local_Size)
     checkGlobalRectangularMatrixInterfaceMPI(A, x);
 }
 
-TEST(MPITest_Math, CheckInterface_PETSc_Rectangular_Matrix_Global_Size)
+TEST(MPI_Math, CheckInterface_PETSc_Rectangular_Matrix_Global_Size)
 {
     MathLib::PETScMatrixOption opt;
     opt.d_nz = 3;
diff --git a/Tests/MathLib/TestGlobalVectorInterface.cpp b/Tests/MathLib/TestGlobalVectorInterface.cpp
index dca634e5a221c3aa4ac13fb98e7faad4e4912ab9..4430b29ec106ddf234bff0c4c8c46320f849355f 100644
--- a/Tests/MathLib/TestGlobalVectorInterface.cpp
+++ b/Tests/MathLib/TestGlobalVectorInterface.cpp
@@ -260,7 +260,7 @@ void checkPETScVectorExplictGhostID()
 
 //--------------------------------------------
 #if defined(USE_PETSC)
-TEST(MPITest_Math, PETScVectorPatitionedAutomatically)
+TEST(MPI_Math, PETScVectorPatitionedAutomatically)
 {
     int msize;
     MPI_Comm_size(PETSC_COMM_WORLD, &msize);
@@ -272,7 +272,7 @@ TEST(MPITest_Math, PETScVectorPatitionedAutomatically)
     checkPETScVectorNoExplictGhostID<MathLib::PETScVector>(x, msize);
 }
 
-TEST(MPITest_Math, PETScVectorFixedPartition)
+TEST(MPI_Math, PETScVectorFixedPartition)
 {
     int msize;
     MPI_Comm_size(PETSC_COMM_WORLD, &msize);
@@ -291,7 +291,7 @@ TEST(MPITest_Math, PETScVectorFixedPartition)
     checkPETScVectorNoExplictGhostID<MathLib::PETScVector>(x_fixed_p, msize);
 }
 
-TEST(MPITest_Math, CheckPETScVectorExplictGhostID)
+TEST(MPI_Math, CheckPETScVectorExplictGhostID)
 {
     checkPETScVectorExplictGhostID<MathLib::PETScVector>();
 }
diff --git a/Tests/MathLib/TestLinearSolver.cpp b/Tests/MathLib/TestLinearSolver.cpp
index c3140099cc869cdfd179f9f30fbc99f3c1f79b41..3aada8e8b0db8d268895ef77c75f979e4ce2b644 100644
--- a/Tests/MathLib/TestLinearSolver.cpp
+++ b/Tests/MathLib/TestLinearSolver.cpp
@@ -347,7 +347,7 @@ TEST(Math, CheckInterface_EigenLis)
 #endif
 
 #ifdef USE_PETSC
-TEST(MPITest_Math, CheckInterface_PETSc_Linear_Solver_basic)
+TEST(MPI_Math, CheckInterface_PETSc_Linear_Solver_basic)
 {
     MathLib::PETScMatrixOption opt;
     opt.d_nz = 2;
@@ -377,7 +377,7 @@ TEST(MPITest_Math, CheckInterface_PETSc_Linear_Solver_basic)
         A, b, "" /*prefix, not specified*/, getConfigTree(xml));
 }
 
-TEST(MPITest_Math, CheckInterface_PETSc_Linear_Solver_chebyshev_sor)
+TEST(MPI_Math, CheckInterface_PETSc_Linear_Solver_chebyshev_sor)
 {
     MathLib::PETScMatrixOption opt;
     opt.d_nz = 2;
@@ -407,7 +407,7 @@ TEST(MPITest_Math, CheckInterface_PETSc_Linear_Solver_chebyshev_sor)
         A, b, "" /*prefix, not specified*/, getConfigTree(xml));
 }
 
-TEST(MPITest_Math, CheckInterface_PETSc_Linear_Solver_gmres_amg)
+TEST(MPI_Math, CheckInterface_PETSc_Linear_Solver_gmres_amg)
 {
     MathLib::PETScMatrixOption opt;
     opt.d_nz = 2;
diff --git a/Tests/NumLib/TestComponentNorms.cpp b/Tests/NumLib/TestComponentNorms.cpp
index 22a43bed9beabd2d74593aa0f2a73372974a273c..4373acaca7053a0e1ce61850bab5e250d1388e45 100644
--- a/Tests/NumLib/TestComponentNorms.cpp
+++ b/Tests/NumLib/TestComponentNorms.cpp
@@ -124,7 +124,7 @@ void do_test(unsigned const num_components,
 #ifndef USE_PETSC
 TEST(NumLib, ComponentNormSingleComponent)
 #else
-TEST(MPITest_NumLib, ComponentNormSingleComponent)
+TEST(MPI_NumLib, ComponentNormSingleComponent)
 #endif
 {
     unsigned const num_components = 1;
@@ -145,7 +145,7 @@ TEST(MPITest_NumLib, ComponentNormSingleComponent)
 #ifndef USE_PETSC
 TEST(NumLib, ComponentNormMultiComponent1)
 #else
-TEST(MPITest_NumLib, ComponentNormMultiComponent1)
+TEST(MPI_NumLib, ComponentNormMultiComponent1)
 #endif
 {
     unsigned const num_components = 3;
@@ -162,7 +162,7 @@ TEST(MPITest_NumLib, ComponentNormMultiComponent1)
 #ifndef USE_PETSC
 TEST(NumLib, ComponentNormMultiComponent2)
 #else
-TEST(MPITest_NumLib, ComponentNormMultiComponent2)
+TEST(MPI_NumLib, ComponentNormMultiComponent2)
 #endif
 {
     unsigned const num_components = 3;
@@ -179,7 +179,7 @@ TEST(MPITest_NumLib, ComponentNormMultiComponent2)
 #ifndef USE_PETSC
 TEST(NumLib, ComponentNormMultiComponentInfinity)
 #else
-TEST(MPITest_NumLib, ComponentNormMultiComponentInfinity)
+TEST(MPI_NumLib, ComponentNormMultiComponentInfinity)
 #endif
 {
     unsigned const num_components = 3;