From a282e79c8b0e8f6a3e736a4c39f364e65be0edb6 Mon Sep 17 00:00:00 2001
From: Christian Godenschwager <christian.godenschwager@fau.de>
Date: Fri, 19 Feb 2016 13:35:39 +0100
Subject: [PATCH] Added support for ParMetis library

Also added wrapper for Metis & ParMetis to prevent it from polluting the global name space
---
 CMakeLists.txt                              |  23 +++
 src/blockforest/GlobalLoadBalancing.h       | 118 ++++++------
 src/blockforest/MetisWrapper.h              |  64 -------
 src/core/load_balancing/MetisWrapper.cpp    | 188 ++++++++++++++++++++
 src/core/load_balancing/MetisWrapper.h      |  72 ++++++++
 src/core/load_balancing/ParMetisWrapper.cpp | 180 +++++++++++++++++++
 src/core/load_balancing/ParMetisWrapper.h   |  55 ++++++
 src/waLBerlaDefinitions.in.h                |   1 +
 tests/core/CMakeLists.txt                   |  16 ++
 tests/core/load_balancing/MetisTest.cpp     | 149 ++++++++++++++++
 tests/core/load_balancing/ParMetisTest.cpp  | 170 ++++++++++++++++++
 11 files changed, 909 insertions(+), 127 deletions(-)
 delete mode 100644 src/blockforest/MetisWrapper.h
 create mode 100644 src/core/load_balancing/MetisWrapper.cpp
 create mode 100644 src/core/load_balancing/MetisWrapper.h
 create mode 100644 src/core/load_balancing/ParMetisWrapper.cpp
 create mode 100644 src/core/load_balancing/ParMetisWrapper.h
 create mode 100644 tests/core/load_balancing/MetisTest.cpp
 create mode 100644 tests/core/load_balancing/ParMetisTest.cpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 961780993..09791376e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -70,6 +70,7 @@ option ( WALBERLA_BUILD_TUTORIALS           "Build Tutorials"
 
 option ( WALBERLA_BUILD_WITH_MPI            "Build with MPI"                                  ON )
 option ( WALBERLA_BUILD_WITH_METIS          "Build with metis graph partitioner"             OFF )
+option ( WALBERLA_BUILD_WITH_PARMETIS       "Build with ParMetis graph partitioner"          OFF )                                            
 
 option ( WALBERLA_BUILD_WITH_GPROF          "Enables gprof"                                      )
 option ( WALBERLA_BUILD_WITH_GCOV           "Enables gcov"                                       )
@@ -853,6 +854,28 @@ else()
     set ( METIS_FOUND OFF CACHE BOOL "Metis found" FORCE )
 endif()
 
+
+if ( WALBERLA_BUILD_WITH_PARMETIS )
+   find_path(PARMETIS_INCLUDE_DIR parmetis.h
+      /usr/local/include
+      /usr/include
+      ${PARMETIS_ROOT}/include
+      $ENV{PARMETIS_ROOT}/include
+   )
+  
+  find_library(PARMETIS_LIBRARY parmetis
+    /usr/local/lib
+    /usr/lib
+    ${PARMETIS_ROOT}/lib
+    $ENV{PARMETIS_ROOT}/lib
+  )
+  
+  if( PARMETIS_INCLUDE_DIR AND PARMETIS_LIBRARY AND METIS_LIBRARY )
+    include_directories( ${PARMETIS_INCLUDE_DIR} )
+    list ( APPEND SERVICE_LIBS ${PARMETIS_LIBRARY} ${METIS_LIBRARY} )
+  endif()
+endif()
+
 ############################################################################################################################
 
 
diff --git a/src/blockforest/GlobalLoadBalancing.h b/src/blockforest/GlobalLoadBalancing.h
index 8ad791cd0..23f79d7b5 100644
--- a/src/blockforest/GlobalLoadBalancing.h
+++ b/src/blockforest/GlobalLoadBalancing.h
@@ -24,12 +24,12 @@
 #include "waLBerlaDefinitions.h" // for macro WALBERLA_BUILD_WITH_METIS
 
 #include "BlockID.h"
-#include "MetisWrapper.h"
 #include "Types.h"
 
 #include "core/Abort.h"
 #include "core/debug/Debug.h"
 #include "core/logging/Logging.h"
+#include "core/load_balancing/MetisWrapper.h"
 #include "core/math/KahanSummation.h"
 
 #include <boost/function.hpp>
@@ -144,10 +144,10 @@ public:
    static void metis2( const std::vector< BLOCK* >& blocks, const uint_t numberOfProcesses, const CommFunction & commFunction );
 private:
 
-   static inline uint_t metisAdaptPartVector( std::vector< idx_t >& part, const uint_t numberOfProcesses );
+   static inline uint_t metisAdaptPartVector( std::vector< int64_t >& part, const uint_t numberOfProcesses );
 
    template< typename BLOCK >
-   static memory_t metisMaxMemory( const std::vector< BLOCK* >& blocks, const uint_t numberOfProcesses, const std::vector< idx_t >& part );
+   static memory_t metisMaxMemory( const std::vector< BLOCK* >& blocks, const uint_t numberOfProcesses, const std::vector< int64_t >& part );
 
    static inline std::ostream &  metisErrorCodeToStream( int errorCode, std::ostream & oss );
    static inline std::string     metisErrorCodeToString( int errorCode );
@@ -781,51 +781,51 @@ uint_t GlobalLoadBalancing::metis( const std::vector< BLOCK* >& blocks, const me
 
    // translate block forest to metis graph structure
 
-   idx_t nvtxs = numeric_cast<idx_t>( blocks.size() ); // number of vertices in the graph
-   idx_t ncon  = 2;                                    // number of balancing constraints
+   int64_t nvtxs = numeric_cast<int64_t>( blocks.size() ); // number of vertices in the graph
+   int64_t ncon  = 2;                                    // number of balancing constraints
 
    for( uint_t i = 0; i != blocks.size(); ++i )
       blocks[i]->setIndex(i);
 
-   std::vector< idx_t > xadj;                         // adjacency structure ...
-   std::vector< idx_t > adjncy;                       // ... of the graph
-   std::vector< idx_t > vwgt( uint_c(ncon * nvtxs) ); // weights of the vertices
-   std::vector< idx_t > adjwgt;                       // weights of the edges
+   std::vector< int64_t > xadj;                         // adjacency structure ...
+   std::vector< int64_t > adjncy;                       // ... of the graph
+   std::vector< int64_t > vwgt( uint_c(ncon * nvtxs) ); // weights of the vertices
+   std::vector< int64_t > adjwgt;                       // weights of the edges
 
    xadj.push_back(0);
    for( uint_t i = 0; i != blocks.size(); ++i ) {
       const BLOCK* block = blocks[i];
 
-      idx_t next = xadj.back() + numeric_cast< idx_t >( block->getNeighborhoodSize() );
+      int64_t next = xadj.back() + numeric_cast< int64_t >( block->getNeighborhoodSize() );
       xadj.push_back( next );
 
       for( uint_t j = 0; j != block->getNeighborhoodSize(); ++j ) {
-         adjncy.push_back( numeric_cast<idx_t>( block->getNeighbor(j)->getIndex() ) );
+         adjncy.push_back( numeric_cast<int64_t>( block->getNeighbor(j)->getIndex() ) );
          adjwgt.push_back( metisConfig.communicationFunction().empty() ? 1 :
-                           ( numeric_cast<idx_t>( static_cast< memory_t >(0.5) + metisConfig.communicationFunction()( block, block->getNeighbor(j) ) ) ) );
+                           ( numeric_cast<int64_t>( static_cast< memory_t >(0.5) + metisConfig.communicationFunction()( block, block->getNeighbor(j) ) ) ) );
       }
 
-      vwgt[ i * 2 ]     = numeric_cast< idx_t >( static_cast< workload_t >(0.5) + blocks[i]->getWorkload() );
-      vwgt[ i * 2 + 1 ] = numeric_cast< idx_t >( static_cast<  memory_t  >(0.5) + blocks[i]->getMemory() );
+      vwgt[ i * 2 ]     = numeric_cast< int64_t >( static_cast< workload_t >(0.5) + blocks[i]->getWorkload() );
+      vwgt[ i * 2 + 1 ] = numeric_cast< int64_t >( static_cast<  memory_t  >(0.5) + blocks[i]->getMemory() );
    }
 
-   idx_t nparts =  numeric_cast< idx_t >( numberOfProcesses ); // number of parts to partition the graph
-   idx_t objval;
+   int64_t nparts =  numeric_cast< int64_t >( numberOfProcesses ); // number of parts to partition the graph
+   int64_t objval;
 
-   std::vector< idx_t > part( uint_c(nvtxs) );
+   std::vector< int64_t > part( uint_c(nvtxs) );
 
    real_t maxUbvec = metisConfig.maxUbvec();
    real_t ubvec[]  = { real_c(1.01), maxUbvec };
 
    // first call to METIS: always try to balance the workload as good as possible, but allow large imbalances concerning the memory (-> ubvec[1])
 
-   int ret = METIS_PartGraphRecursive( &nvtxs, &ncon, &(xadj[0]), &(adjncy[0]), &(vwgt[0]), NULL, &(adjwgt[0]), &nparts, NULL,
-                                       &(ubvec[0]), NULL /*idx t *options*/, &objval, &(part[0]) );
+   int ret = core::METIS_PartGraphRecursive( &nvtxs, &ncon, &(xadj[0]), &(adjncy[0]), &(vwgt[0]), NULL, &(adjwgt[0]), &nparts, NULL,
+                                             &(ubvec[0]), NULL /*idx t *options*/, &objval, &(part[0]) );
 
    // if METIS was successful AND the memory limit of each process is not violated (which is highly unlikely due to a large value for ubvec[1])
    // then the algorithm is finished
 
-   if( ret == METIS_OK && metisMaxMemory( blocks, numberOfProcesses, part ) <= memoryLimit ) {
+   if( ret == core::METIS_OK && metisMaxMemory( blocks, numberOfProcesses, part ) <= memoryLimit ) {
 
       nProcesses = metisAdaptPartVector( part, numberOfProcesses );
 
@@ -844,16 +844,16 @@ uint_t GlobalLoadBalancing::metis( const std::vector< BLOCK* >& blocks, const me
    real_t minUbvec = real_t(1);
    ubvec[1] = minUbvec;
 
-   ret = METIS_PartGraphRecursive( &nvtxs, &ncon, &(xadj[0]), &(adjncy[0]), &(vwgt[0]), NULL, &(adjwgt[0]), &nparts, NULL,
-                                   &(ubvec[0]), NULL /*idx t *options*/, &objval, &(part[0]) );
+   ret = core::METIS_PartGraphRecursive( &nvtxs, &ncon, &(xadj[0]), &(adjncy[0]), &(vwgt[0]), NULL, &(adjwgt[0]), &nparts, NULL,
+                                         &(ubvec[0]), NULL /*idx t *options*/, &objval, &(part[0]) );
 
    // ... if this doesn't work OR if the memory limit is still violated then METIS is unable to find a valid partitioning
 
-   if( ret != METIS_OK ) {
+   if( ret != core::METIS_OK ) {
 
       std::string error( "METIS_ERROR" );
-      if( ret == METIS_ERROR_INPUT ) error.assign( "METIS_ERROR_INPUT" );
-      else if( ret == METIS_ERROR_MEMORY ) error.assign( "METIS_ERROR_MEMORY" );
+      if( ret == core::METIS_ERROR_INPUT ) error.assign( "METIS_ERROR_INPUT" );
+      else if( ret == core::METIS_ERROR_MEMORY ) error.assign( "METIS_ERROR_MEMORY" );
 
       // DEBUG_LOGGING_SECTION { std::cout << "ERROR: static load balancing with METIS failed (" << error << ")" << std::endl; }
       return 0;
@@ -881,10 +881,10 @@ uint_t GlobalLoadBalancing::metis( const std::vector< BLOCK* >& blocks, const me
 
       ubvec[1] = ( maxUbvec + minUbvec ) / real_c(2);
 
-      ret = METIS_PartGraphRecursive( &nvtxs, &ncon, &(xadj[0]), &(adjncy[0]), &(vwgt[0]), NULL, &(adjwgt[0]), &nparts, NULL,
-                                      &(ubvec[0]), NULL /*idx t *options*/, &objval, &(part[0]) );
+      ret = core::METIS_PartGraphRecursive( &nvtxs, &ncon, &(xadj[0]), &(adjncy[0]), &(vwgt[0]), NULL, &(adjwgt[0]), &nparts, NULL,
+                                            &(ubvec[0]), NULL /*idx t *options*/, &objval, &(part[0]) );
 
-      if( ret == METIS_OK && metisMaxMemory( blocks, numberOfProcesses, part ) <= memoryLimit ) {
+      if( ret == core::METIS_OK && metisMaxMemory( blocks, numberOfProcesses, part ) <= memoryLimit ) {
 
          nProcesses = metisAdaptPartVector( part, numberOfProcesses );
 
@@ -913,8 +913,8 @@ void GlobalLoadBalancing::metis2( const std::vector< BLOCK* >& blocks, const uin
    if( blocks.empty() )
       return;
   
-   idx_t nvtxs = 0; // number of vertices in the graph
-   idx_t ncon  = 1; // number of balancing constraints
+   int64_t nvtxs = 0; // number of vertices in the graph
+   int64_t ncon  = 1; // number of balancing constraints
    
    uint_t j = 0;
    for( uint_t i = 0; i != blocks.size(); ++i )
@@ -947,10 +947,10 @@ void GlobalLoadBalancing::metis2( const std::vector< BLOCK* >& blocks, const uin
 
    commFunction( blockPairs, communicationWeights );
    
-   std::vector< idx_t > xadj;   // adjacency structure ...
-   std::vector< idx_t > adjncy; // ... of the graph
-   std::vector< idx_t > vwgt;   // weights of the vertices
-   std::vector< idx_t > adjwgt; // weights of the edges
+   std::vector< int64_t > xadj;   // adjacency structure ...
+   std::vector< int64_t > adjncy; // ... of the graph
+   std::vector< int64_t > vwgt;   // weights of the vertices
+   std::vector< int64_t > adjwgt; // weights of the edges
    
    xadj.push_back( 0 );
    uint_t commIdx = 0;
@@ -963,7 +963,7 @@ void GlobalLoadBalancing::metis2( const std::vector< BLOCK* >& blocks, const uin
    
       ++nvtxs;
       xadj.push_back( xadj.back() );
-      vwgt.push_back( numeric_cast<idx_t>( blocks[ i ]->getWorkload() ) );
+      vwgt.push_back( numeric_cast<int64_t>( blocks[ i ]->getWorkload() ) );
    
       for( uint_t k = 0; k != block->getNeighborhoodSize(); ++k ) 
       {
@@ -971,10 +971,10 @@ void GlobalLoadBalancing::metis2( const std::vector< BLOCK* >& blocks, const uin
          {
             if( communicationWeights[ commIdx ] > real_t(0) )
             {
-               adjncy.push_back( numeric_cast<idx_t>( block->getNeighbor( k )->getIndex() ) );
+               adjncy.push_back( numeric_cast<int64_t>( block->getNeighbor( k )->getIndex() ) );
                WALBERLA_ASSERT_LESS( commIdx, communicationWeights.size() );
-               WALBERLA_ASSERT_GREATER( numeric_cast<idx_t>( communicationWeights[ commIdx ] ), idx_t(0) );
-               adjwgt.push_back( numeric_cast<idx_t>( communicationWeights[ commIdx ] ) );
+               WALBERLA_ASSERT_GREATER( numeric_cast<int64_t>( communicationWeights[ commIdx ] ), int64_t(0) );
+               adjwgt.push_back( numeric_cast<int64_t>( communicationWeights[ commIdx ] ) );
                xadj.back() += 1;
             }
             ++commIdx;
@@ -987,21 +987,21 @@ void GlobalLoadBalancing::metis2( const std::vector< BLOCK* >& blocks, const uin
    WALBERLA_ASSERT_EQUAL( adjncy.size(), adjwgt.size() );
    WALBERLA_ASSERT_EQUAL( adjncy.size(), communicationWeights.size() );
    
-   idx_t nparts = numeric_cast<idx_t>( numberOfProcesses ); // number of parts to partition the graph
-   idx_t objval;
+   int64_t nparts = numeric_cast<int64_t>( numberOfProcesses ); // number of parts to partition the graph
+   int64_t objval;
    
-   std::vector< idx_t > part( uint_c(nvtxs) );
+   std::vector< int64_t > part( uint_c(nvtxs) );
    
-   idx_t options[ METIS_NOPTIONS ];
-   METIS_SetDefaultOptions( options );
-   options[ METIS_OPTION_NITER ] = 1000;
-   options[ METIS_OPTION_NSEPS ] = 100;
-   options[ METIS_OPTION_NCUTS ] = 100;
+   int64_t options[ METIS_NOPTIONS ];
+   core::METIS_SetDefaultOptions( options );
+   options[ core::METIS_OPTION_NITER ] = 1000;
+   options[ core::METIS_OPTION_NSEPS ] = 100;
+   options[ core::METIS_OPTION_NCUTS ] = 100;
    
-   int ret = METIS_PartGraphKway( &nvtxs, &ncon, &( xadj[ 0 ] ), &( adjncy[ 0 ] ), &( vwgt[ 0 ] ), NULL, &( adjwgt[0] ), 
+   int ret = core::METIS_PartGraphKway( &nvtxs, &ncon, &( xadj[ 0 ] ), &( adjncy[ 0 ] ), &( vwgt[ 0 ] ), NULL, &( adjwgt[0] ), 
                                   &nparts, NULL, NULL, options, &objval, &( part[ 0 ] ) );
 
-   if( ret != METIS_OK ) 
+   if( ret != core::METIS_OK ) 
    {
       WALBERLA_ABORT( "METIS partitioning failed! Error: " << metisErrorCodeToString( ret ) );
    }
@@ -1028,24 +1028,16 @@ void GlobalLoadBalancing::metis2( const std::vector< BLOCK* >& blocks, const uin
 
 std::ostream & GlobalLoadBalancing::metisErrorCodeToStream( int errorCode, std::ostream & oss )
 {
-   switch( errorCode )
-   {
-   case METIS_OK:
+   if( errorCode == core::METIS_OK)
       oss << "OK, no METIS error";
-      break;
-   case METIS_ERROR_INPUT:
+   else if( errorCode == core::METIS_ERROR_INPUT )
       oss << "Error in METIS input";
-      break;
-   case METIS_ERROR_MEMORY:
+   else if (errorCode == core::METIS_ERROR_MEMORY )
       oss << "METIS could not allocate enough memory";
-      break;
-   case METIS_ERROR:
+   else if (errorCode == core::METIS_ERROR )
       oss << "Unknown type of error";
-      break;
-   default:
+   else
       oss << "Unknown error code";
-      break;
-   }
 
    return oss;
 }
@@ -1060,7 +1052,7 @@ std::string GlobalLoadBalancing::metisErrorCodeToString( int errorCode )
 
 
 
-uint_t GlobalLoadBalancing::metisAdaptPartVector( std::vector< idx_t >& part, const uint_t numberOfProcesses )
+uint_t GlobalLoadBalancing::metisAdaptPartVector( std::vector< int64_t >& part, const uint_t numberOfProcesses )
 {
    std::vector<bool> hasBlock( numberOfProcesses, false );
    for( uint_t i = 0; i != part.size(); ++i )
@@ -1090,7 +1082,7 @@ uint_t GlobalLoadBalancing::metisAdaptPartVector( std::vector< idx_t >& part, co
 
 template< typename BLOCK >
 memory_t GlobalLoadBalancing::metisMaxMemory( const std::vector< BLOCK* >& blocks, const uint_t numberOfProcesses,
-                                              const std::vector< idx_t >& part ) {
+                                              const std::vector< int64_t >& part ) {
 
    WALBERLA_ASSERT_EQUAL( blocks.size(), part.size() );
 
diff --git a/src/blockforest/MetisWrapper.h b/src/blockforest/MetisWrapper.h
deleted file mode 100644
index 03cb3dfaa..000000000
--- a/src/blockforest/MetisWrapper.h
+++ /dev/null
@@ -1,64 +0,0 @@
-//======================================================================================================================
-//
-//  This file is part of waLBerla. waLBerla is free software: you can 
-//  redistribute it and/or modify it under the terms of the GNU General Public
-//  License as published by the Free Software Foundation, either version 3 of 
-//  the License, or (at your option) any later version.
-//  
-//  waLBerla is distributed in the hope that it will be useful, but WITHOUT 
-//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
-//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
-//  for more details.
-//  
-//  You should have received a copy of the GNU General Public License along
-//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
-//
-//! \file MetisWrapper.h
-//! \ingroup blockforest
-//! \author Christian Godenschwager <christian.godenschwager@fau.de>
-//
-//======================================================================================================================
-
-#pragma once
-
-#include "waLBerlaDefinitions.h" // for macro WALBERLA_BUILD_WITH_METIS
-
-#ifdef _MSC_VER
-   #pragma push_macro( "INT32_MIN" )
-   #pragma push_macro( "INT32_MAX" )
-   #pragma push_macro( "INT64_MIN" )
-   #pragma push_macro( "INT64_MAX" )
-
-   #ifdef INT32_MIN
-      #undef INT32_MIN
-   #endif
-
-   #ifdef INT32_MAX
-      #undef INT32_MAX
-   #endif
-
-   #ifdef INT64_MIN
-      #undef INT64_MIN
-   #endif
-
-   #ifdef INT64_MAX
-      #undef INT64_MAX
-   #endif
-#endif
-
-#ifdef WALBERLA_BUILD_WITH_METIS
-// external software includes
-#include "metis.h"
-
-#include <boost/type_traits/is_same.hpp>
-
-static_assert( boost::is_same< ::real_t, ::walberla::real_t >::value , "The width of waLBerla's real_t data type does not match the width of METIS' real_t data type. " \
-                                                                       "Please adapt either waLBerla's CMake option WALBERLA_DOUBLE_ACCURACY or change REALTYPEWIDTH in metis.h." );													  
-#endif
-
-#ifdef _MSC_VER
-#pragma pop_macro( "INT64_MAX" )
-#pragma pop_macro( "INT64_MIN" )
-#pragma pop_macro( "INT32_MAX" )
-#pragma pop_macro( "INT32_MIN" )
-#endif
\ No newline at end of file
diff --git a/src/core/load_balancing/MetisWrapper.cpp b/src/core/load_balancing/MetisWrapper.cpp
new file mode 100644
index 000000000..4d8d08a11
--- /dev/null
+++ b/src/core/load_balancing/MetisWrapper.cpp
@@ -0,0 +1,188 @@
+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file MetisWrapper.cpp
+//! \ingroup blockforest
+//! \author Christian Godenschwager <christian.godenschwager@fau.de>
+//
+//======================================================================================================================
+
+#include "MetisWrapper.h"
+#define WALBERLA_METIS_NOPTIONS METIS_NOPTIONS
+#undef METIS_NOPTIONS
+
+#include "waLBerlaDefinitions.h"
+
+#include "core/Abort.h"
+
+#include <boost/type_traits/is_same.hpp>
+
+#ifdef WALBERLA_BUILD_WITH_METIS
+
+#ifdef _MSC_VER
+#  pragma push_macro( "INT32_MIN" )
+#  pragma push_macro( "INT32_MAX" )
+#  pragma push_macro( "INT64_MIN" )
+#  pragma push_macro( "INT64_MAX" )
+
+#  ifdef INT32_MIN
+#     undef INT32_MIN
+#  endif
+
+#  ifdef INT32_MAX
+#     undef INT32_MAX
+#  endif
+
+#  ifdef INT64_MIN
+#     undef INT64_MIN
+#  endif
+
+#  ifdef INT64_MAX
+#     undef INT64_MAX
+#  endif
+#endif
+
+#ifdef WALBERLA_BUILD_WITH_METIS
+#  include "metis.h"
+#endif
+
+#ifdef _MSC_VER
+#  pragma pop_macro( "INT64_MAX" )
+#  pragma pop_macro( "INT64_MIN" )
+#  pragma pop_macro( "INT32_MAX" )
+#  pragma pop_macro( "INT32_MIN" )
+#endif
+
+static_assert(WALBERLA_METIS_NOPTIONS == METIS_NOPTIONS, "The macro METIS_NOPTIONS defined in blockforest/loadbalancing/MetisWRapper.h does not match the value in metis.h!");
+
+#endif // WALBERLA_BUILD_WITH_METIS
+
+namespace walberla {
+namespace core {
+
+#ifdef WALBERLA_BUILD_WITH_METIS
+
+int METIS_PartGraphKway( ::walberla::int64_t *nvtxs, ::walberla::int64_t *ncon, ::walberla::int64_t *xadj, ::walberla::int64_t *adjncy,
+                         ::walberla::int64_t *vwgt, ::walberla::int64_t *vsize, ::walberla::int64_t *adjwgt, ::walberla::int64_t *nparts,
+                         double *tpwgts, double *ubvec, ::walberla::int64_t *options, ::walberla::int64_t *edgecut, ::walberla::int64_t *part)
+{
+   static_assert(boost::is_same< ::walberla::int64_t, ::idx_t >::value, "You have to compile the metis library with 64-bit wide integer type support!");
+   static_assert(boost::is_same< double, ::real_t >::value, "You have to compile the metis library with 64-bit wide floating-point type support!");
+
+   return ::METIS_PartGraphKway( nvtxs, ncon, xadj, adjncy, vwgt, vsize, adjwgt, nparts, tpwgts, ubvec, options, edgecut, part );
+}
+
+
+int METIS_PartGraphRecursive( ::walberla::int64_t *nvtxs, ::walberla::int64_t *ncon, ::walberla::int64_t *xadj, ::walberla::int64_t *adjncy,
+                              ::walberla::int64_t *vwgt, ::walberla::int64_t *vsize, ::walberla::int64_t *adjwgt, ::walberla::int64_t *nparts,
+                              double *tpwgts, double *ubvec, ::walberla::int64_t *options, ::walberla::int64_t *edgecut, ::walberla::int64_t *part)
+{
+   static_assert(boost::is_same< ::walberla::int64_t, ::idx_t >::value, "You have to compile the metis library with 64-bit wide integer type support!");
+   static_assert(boost::is_same< double, ::real_t >::value, "You have to compile the metis library with 64-bit wide floating-point type support!");
+
+   return ::METIS_PartGraphRecursive( nvtxs, ncon, xadj, adjncy, vwgt, vsize, adjwgt, nparts, tpwgts, ubvec, options, edgecut, part );
+}
+
+int METIS_SetDefaultOptions( ::walberla::int64_t *options )
+{
+   static_assert(boost::is_same< ::walberla::int64_t, ::idx_t >::value, "You have to compile the metis library with 64-bit wide integer type support!");
+   
+   return ::METIS_SetDefaultOptions( options );
+}
+
+
+
+const int METIS_OK           = ::METIS_OK;
+const int METIS_ERROR        = ::METIS_ERROR;
+const int METIS_ERROR_INPUT  = ::METIS_ERROR_INPUT;
+const int METIS_ERROR_MEMORY = ::METIS_ERROR_MEMORY;
+
+const int METIS_OPTION_PTYPE      = ::METIS_OPTION_PTYPE;   
+const int METIS_OPTION_OBJTYPE    = ::METIS_OPTION_OBJTYPE;
+const int METIS_OPTION_CTYPE      = ::METIS_OPTION_CTYPE; 
+const int METIS_OPTION_IPTYPE     = ::METIS_OPTION_IPTYPE;   
+const int METIS_OPTION_RTYPE      = ::METIS_OPTION_RTYPE;  
+const int METIS_OPTION_DBGLVL     = ::METIS_OPTION_DBGLVL;   
+const int METIS_OPTION_NITER      = ::METIS_OPTION_NITER;  
+const int METIS_OPTION_NCUTS      = ::METIS_OPTION_NCUTS;   
+const int METIS_OPTION_SEED       = ::METIS_OPTION_SEED;   
+const int METIS_OPTION_NO2HOP     = ::METIS_OPTION_NO2HOP;   
+const int METIS_OPTION_MINCONN    = ::METIS_OPTION_MINCONN;  
+const int METIS_OPTION_CONTIG     = ::METIS_OPTION_CONTIG; 
+const int METIS_OPTION_COMPRESS   = ::METIS_OPTION_COMPRESS; 
+const int METIS_OPTION_CCORDER    = ::METIS_OPTION_CCORDER;
+const int METIS_OPTION_PFACTOR    = ::METIS_OPTION_PFACTOR; 
+const int METIS_OPTION_NSEPS      = ::METIS_OPTION_NSEPS; 
+const int METIS_OPTION_UFACTOR    = ::METIS_OPTION_UFACTOR;  
+const int METIS_OPTION_NUMBERING  = ::METIS_OPTION_NUMBERING;
+const int METIS_OPTION_HELP       = ::METIS_OPTION_HELP;
+const int METIS_OPTION_TPWGTS     = ::METIS_OPTION_TPWGTS;   
+const int METIS_OPTION_NCOMMON    = ::METIS_OPTION_NCOMMON;  
+const int METIS_OPTION_NOOUTPUT   = ::METIS_OPTION_NOOUTPUT; 
+const int METIS_OPTION_BALANCE    = ::METIS_OPTION_BALANCE;
+const int METIS_OPTION_GTYPE      = ::METIS_OPTION_GTYPE;    
+const int METIS_OPTION_UBVEC      = ::METIS_OPTION_UBVEC;   
+
+
+#else // build without Metis
+
+int METIS_PartGraphKway( ::walberla::int64_t * /*nvtxs*/, ::walberla::int64_t * /*ncon*/, ::walberla::int64_t * /*xadj*/, ::walberla::int64_t * /*adjncy*/,
+                         ::walberla::int64_t * /*vwgt*/, ::walberla::int64_t * /*vsize*/, ::walberla::int64_t * /*adjwgt*/, ::walberla::int64_t * /*nparts*/,
+                         double * /*tpwgts*/, double * /*ubvec*/, ::walberla::int64_t * /*options*/, ::walberla::int64_t * /*edgecut*/, ::walberla::int64_t * /*part*/)
+{
+    WALBERLA_ABORT( "You are trying to use Metis functionality but waLBerla is not configured to use it. Set 'WALBERLA_BUILD_WITH_METIS' to 'ON' in your CMake cache to build against an installed version of Metis!" );
+}
+
+int METIS_SetDefaultOptions( ::walberla::int64_t * /*options*/ )
+{
+   WALBERLA_ABORT("You are trying to use Metis functionality but waLBerla is not configured to use it. Set 'WALBERLA_BUILD_WITH_METIS' to 'ON' in your CMake cache to build against an installed version of Metis!");
+}
+
+const int METIS_OK           = 0;
+const int METIS_ERROR        = 0;
+const int METIS_ERROR_INPUT  = 0;
+const int METIS_ERROR_MEMORY = 0;
+
+const int METIS_OPTION_PTYPE      = 0;
+const int METIS_OPTION_OBJTYPE    = 0;
+const int METIS_OPTION_CTYPE      = 0;
+const int METIS_OPTION_IPTYPE     = 0; 
+const int METIS_OPTION_RTYPE      = 0;
+const int METIS_OPTION_DBGLVL     = 0; 
+const int METIS_OPTION_NITER      = 0;
+const int METIS_OPTION_NCUTS      = 0;
+const int METIS_OPTION_SEED       = 0;
+const int METIS_OPTION_NO2HOP     = 0; 
+const int METIS_OPTION_MINCONN    = 0; 
+const int METIS_OPTION_CONTIG     = 0;
+const int METIS_OPTION_COMPRESS   = 0; 
+const int METIS_OPTION_CCORDER    = 0;
+const int METIS_OPTION_PFACTOR    = 0;
+const int METIS_OPTION_NSEPS      = 0;
+const int METIS_OPTION_UFACTOR    = 0; 
+const int METIS_OPTION_NUMBERING  = 0;
+const int METIS_OPTION_HELP       = 0;
+const int METIS_OPTION_TPWGTS     = 0; 
+const int METIS_OPTION_NCOMMON    = 0; 
+const int METIS_OPTION_NOOUTPUT   = 0; 
+const int METIS_OPTION_BALANCE    = 0;
+const int METIS_OPTION_GTYPE      = 0; 
+const int METIS_OPTION_UBVEC      = 0;
+
+#endif
+
+
+} // namespace core
+} // namespace walberla
\ No newline at end of file
diff --git a/src/core/load_balancing/MetisWrapper.h b/src/core/load_balancing/MetisWrapper.h
new file mode 100644
index 000000000..3b4ecbfcd
--- /dev/null
+++ b/src/core/load_balancing/MetisWrapper.h
@@ -0,0 +1,72 @@
+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file MetisWrapper.h
+//! \ingroup blockforest
+//! \author Christian Godenschwager <christian.godenschwager@fau.de>
+//
+//======================================================================================================================
+
+#include "core/DataTypes.h"
+#include "core/mpi/MPIWrapper.h"
+
+namespace walberla {
+namespace core {
+
+int METIS_PartGraphKway( int64_t *nvtxs, int64_t *ncon, int64_t *xadj, int64_t *adjncy, int64_t *vwgt, int64_t *vsize, int64_t *adjwgt,
+                         int64_t *nparts, double *tpwgts, double *ubvec, int64_t *options, int64_t *edgecut, int64_t *part );
+
+int METIS_PartGraphRecursive( int64_t *nvtxs, int64_t *ncon, int64_t *xadj, int64_t *adjncy, int64_t *vwgt, int64_t *vsize, int64_t *adjwgt,
+                              int64_t *nparts, double *tpwgts, double *ubvec, int64_t *options, int64_t *edgecut, int64_t *part );
+
+int METIS_SetDefaultOptions(int64_t *options );
+
+#ifndef METIS_NOPTIONS
+#  define METIS_NOPTIONS 40
+#endif
+
+extern const int METIS_OK;
+extern const int METIS_ERROR;
+extern const int METIS_ERROR_INPUT;
+extern const int METIS_ERROR_MEMORY;
+
+extern const int METIS_OPTION_PTYPE;
+extern const int METIS_OPTION_OBJTYPE;
+extern const int METIS_OPTION_CTYPE;
+extern const int METIS_OPTION_IPTYPE;
+extern const int METIS_OPTION_RTYPE;
+extern const int METIS_OPTION_DBGLVL;
+extern const int METIS_OPTION_NITER;
+extern const int METIS_OPTION_NCUTS;
+extern const int METIS_OPTION_SEED;
+extern const int METIS_OPTION_NO2HOP;
+extern const int METIS_OPTION_MINCONN;
+extern const int METIS_OPTION_CONTIG;
+extern const int METIS_OPTION_COMPRESS;
+extern const int METIS_OPTION_CCORDER;
+extern const int METIS_OPTION_PFACTOR;
+extern const int METIS_OPTION_NSEPS;
+extern const int METIS_OPTION_UFACTOR;
+extern const int METIS_OPTION_NUMBERING;
+extern const int METIS_OPTION_HELP;
+extern const int METIS_OPTION_TPWGTS;
+extern const int METIS_OPTION_NCOMMON;
+extern const int METIS_OPTION_NOOUTPUT;
+extern const int METIS_OPTION_BALANCE;
+extern const int METIS_OPTION_GTYPE;
+extern const int METIS_OPTION_UBVEC;
+
+} // namespace core
+} // namespace walberla
\ No newline at end of file
diff --git a/src/core/load_balancing/ParMetisWrapper.cpp b/src/core/load_balancing/ParMetisWrapper.cpp
new file mode 100644
index 000000000..3c8b525a0
--- /dev/null
+++ b/src/core/load_balancing/ParMetisWrapper.cpp
@@ -0,0 +1,180 @@
+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file ParMetisWrapper.cpp
+//! \ingroup blockforest
+//! \author Christian Godenschwager <christian.godenschwager@fau.de>
+//
+//======================================================================================================================
+
+#include "waLBerlaDefinitions.h"
+
+#ifdef WALBERLA_BUILD_WITH_PARMETIS
+
+#ifdef _MSC_VER
+#  pragma push_macro( "INT32_MIN" )
+#  pragma push_macro( "INT32_MAX" )
+#  pragma push_macro( "INT64_MIN" )
+#  pragma push_macro( "INT64_MAX" )
+
+#  ifdef INT32_MIN
+#     undef INT32_MIN
+#  endif
+
+#  ifdef INT32_MAX
+#     undef INT32_MAX
+#  endif
+
+#  ifdef INT64_MIN
+#     undef INT64_MIN
+#  endif
+
+#  ifdef INT64_MAX
+#     undef INT64_MAX
+#  endif
+#endif
+
+#ifdef WALBERLA_BUILD_WITH_METIS
+#  include "metis.h"
+#endif
+
+#ifdef _MSC_VER
+#  pragma pop_macro( "INT64_MAX" )
+#  pragma pop_macro( "INT64_MIN" )
+#  pragma pop_macro( "INT32_MAX" )
+#  pragma pop_macro( "INT32_MIN" )
+#endif
+
+#include <parmetis.h>
+#endif
+
+#include "ParMetisWrapper.h"
+
+#include "core/Abort.h"
+
+#include <boost/type_traits/is_same.hpp>
+
+
+namespace walberla {
+namespace core {
+
+#ifdef WALBERLA_BUILD_WITH_PARMETIS
+
+int ParMETIS_V3_AdaptiveRepart(
+   ::walberla::int64_t *vtxdist, ::walberla::int64_t *xadj, ::walberla::int64_t *adjncy, ::walberla::int64_t *vwgt,
+   ::walberla::int64_t *vsize, ::walberla::int64_t *adjwgt, ::walberla::int64_t *wgtflag, ::walberla::int64_t *numflag, int64_t *ncon,
+   ::walberla::int64_t *nparts, double *tpwgts, double *ubvec, double *ipc2redist,
+   ::walberla::int64_t *options, ::walberla::int64_t *edgecut, ::walberla::int64_t *part, MPI_Comm *comm )
+{
+
+   static_assert( boost::is_same< ::walberla::int64_t, ::idx_t >::value, "You have to compile the metis library with 64-bit wide integer type support!"        );
+   static_assert( boost::is_same< double, ::real_t >::value,             "You have to compile the metis library with 64-bit wide floating-point type support!" );
+
+   return ::ParMETIS_V3_AdaptiveRepart( vtxdist, xadj, adjncy, vwgt,
+                                        vsize, adjwgt, wgtflag, numflag, ncon,
+                                        nparts, tpwgts, ubvec, ipc2redist,
+                                        options, edgecut, part, comm );
+}
+
+int ParMETIS_V3_PartKway(
+   ::walberla::int64_t *vtxdist, ::walberla::int64_t *xadj, ::walberla::int64_t *adjncy, ::walberla::int64_t *vwgt,
+   ::walberla::int64_t *adjwgt, ::walberla::int64_t *wgtflag, ::walberla::int64_t *numflag, ::walberla::int64_t *ncon, ::walberla::int64_t *nparts,
+   double *tpwgts, double *ubvec, ::walberla::int64_t *options, ::walberla::int64_t *edgecut, ::walberla::int64_t *part,
+   MPI_Comm *comm )
+{
+   static_assert( boost::is_same< ::walberla::int64_t, ::idx_t >::value, "You have to compile the metis library with 64-bit wide integer type support!" );
+   static_assert( boost::is_same< double, ::real_t >::value, "You have to compile the metis library with 64-bit wide floating-point type support!" );
+
+   return ::ParMETIS_V3_PartKway( vtxdist, xadj, adjncy, vwgt,
+                                  adjwgt, wgtflag, numflag, ncon, nparts,
+                                  tpwgts, ubvec, options, edgecut, part, comm );
+
+}
+
+int ParMETIS_V3_PartGeomKway(
+   ::walberla::int64_t *vtxdist, ::walberla::int64_t *xadj, ::walberla::int64_t *adjncy, ::walberla::int64_t *vwgt,
+   ::walberla::int64_t *adjwgt, ::walberla::int64_t *wgtflag, ::walberla::int64_t *numflag, ::walberla::int64_t *ndims, double *xyz,
+   ::walberla::int64_t *ncon, ::walberla::int64_t *nparts, double *tpwgts, double *ubvec, ::walberla::int64_t *options,
+   ::walberla::int64_t *edgecut, ::walberla::int64_t *part, MPI_Comm *comm )
+{
+   static_assert( boost::is_same< ::walberla::int64_t, ::idx_t >::value, "You have to compile the metis library with 64-bit wide integer type support!" );
+   static_assert( boost::is_same< double, ::real_t >::value, "You have to compile the metis library with 64-bit wide floating-point type support!" );
+
+   return ::ParMETIS_V3_PartGeomKway( vtxdist, xadj, adjncy, vwgt,
+                                      adjwgt, wgtflag, numflag, ndims, xyz,
+                                      ncon, nparts, tpwgts, ubvec, options,
+                                      edgecut, part, comm );
+}
+
+int ParMETIS_V3_RefineKway(
+   ::walberla::int64_t *vtxdist, ::walberla::int64_t *xadj, ::walberla::int64_t *adjncy, ::walberla::int64_t *vwgt,
+   ::walberla::int64_t *adjwgt, ::walberla::int64_t *wgtflag, ::walberla::int64_t *numflag, ::walberla::int64_t *ncon, ::walberla::int64_t *nparts,
+   double *tpwgts, double *ubvec, ::walberla::int64_t *options, ::walberla::int64_t *edgecut,
+   ::walberla::int64_t *part, MPI_Comm *comm )
+{
+   static_assert( boost::is_same< ::walberla::int64_t, ::idx_t >::value, "You have to compile the metis library with 64-bit wide integer type support!" );
+   static_assert( boost::is_same< double, ::real_t >::value, "You have to compile the metis library with 64-bit wide floating-point type support!" );
+
+   return ::ParMETIS_V3_RefineKway( vtxdist, xadj, adjncy, vwgt,
+                                    adjwgt, wgtflag, numflag, ncon, nparts,
+                                    tpwgts, ubvec, options, edgecut,
+                                    part, comm );
+}
+
+#else // build without ParMetis
+
+int ParMETIS_V3_AdaptiveRepart(
+   ::walberla::int64_t *, ::walberla::int64_t *, ::walberla::int64_t *, ::walberla::int64_t *,
+   ::walberla::int64_t *, ::walberla::int64_t *, ::walberla::int64_t *, ::walberla::int64_t *, int64_t *,
+   ::walberla::int64_t *, double *, double *, double *,
+   ::walberla::int64_t *, ::walberla::int64_t *, ::walberla::int64_t *, MPI_Comm * )
+{
+
+   WALBERLA_ABORT( "You are trying to use ParMetis functionality but waLBerla is not configured to use it. Set 'WALBERLA_BUILD_WITH_PARMETIS' to 'ON' in your CMake cache to build against an installed version of ParMetis!" );
+}
+
+int ParMETIS_V3_PartKway(
+   ::walberla::int64_t *, ::walberla::int64_t *, ::walberla::int64_t *, ::walberla::int64_t *,
+   ::walberla::int64_t *, ::walberla::int64_t *, ::walberla::int64_t *, ::walberla::int64_t *, ::walberla::int64_t *,
+   double *, double *, ::walberla::int64_t *, ::walberla::int64_t *, ::walberla::int64_t *,
+   MPI_Comm * )
+{
+   WALBERLA_ABORT( "You are trying to use ParMetis functionality but waLBerla is not configured to use it. Set 'WALBERLA_BUILD_WITH_PARMETIS' to 'ON' in your CMake cache to build against an installed version of ParMetis!" );
+
+}
+
+int ParMETIS_V3_PartGeomKway(
+   ::walberla::int64_t *, ::walberla::int64_t *, ::walberla::int64_t *, ::walberla::int64_t *,
+   ::walberla::int64_t *, ::walberla::int64_t *, ::walberla::int64_t *, ::walberla::int64_t *, double *,
+   ::walberla::int64_t *, ::walberla::int64_t *, double *, double *, ::walberla::int64_t *,
+   ::walberla::int64_t *, ::walberla::int64_t *, MPI_Comm * )
+{
+   WALBERLA_ABORT( "You are trying to use ParMetis functionality but waLBerla is not configured to use it. Set 'WALBERLA_BUILD_WITH_PARMETIS' to 'ON' in your CMake cache to build against an installed version of ParMetis!" );
+}
+
+int ParMETIS_V3_RefineKway(
+   ::walberla::int64_t *, ::walberla::int64_t *, ::walberla::int64_t *, ::walberla::int64_t *,
+   ::walberla::int64_t *, ::walberla::int64_t *, ::walberla::int64_t *, ::walberla::int64_t *, ::walberla::int64_t *,
+   double *, double *, ::walberla::int64_t *, ::walberla::int64_t *,
+   ::walberla::int64_t *, MPI_Comm * )
+{
+   WALBERLA_ABORT( "You are trying to use ParMetis functionality but waLBerla is not configured to use it. Set 'WALBERLA_BUILD_WITH_PARMETIS' to 'ON' in your CMake cache to build against an installed version of ParMetis!" );
+}
+
+#endif
+
+
+} // namespace core
+} // namespace walberla
\ No newline at end of file
diff --git a/src/core/load_balancing/ParMetisWrapper.h b/src/core/load_balancing/ParMetisWrapper.h
new file mode 100644
index 000000000..09a7601e8
--- /dev/null
+++ b/src/core/load_balancing/ParMetisWrapper.h
@@ -0,0 +1,55 @@
+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file ParMetisWrapper.h
+//! \ingroup blockforest
+//! \author Christian Godenschwager <christian.godenschwager@fau.de>
+//
+//======================================================================================================================
+
+#include "MetisWrapper.h"
+
+#include "core/DataTypes.h"
+#include "core/mpi/MPIWrapper.h"
+
+namespace walberla {
+namespace core {
+
+int ParMETIS_V3_AdaptiveRepart(
+   int64_t *vtxdist, int64_t *xadj, int64_t *adjncy, int64_t *vwgt,
+   int64_t *vsize, int64_t *adjwgt, int64_t *wgtflag, int64_t *numflag, int64_t *ncon,
+   int64_t *nparts, double *tpwgts, double *ubvec, double *ipc2redist,
+   int64_t *options, int64_t *edgecut, int64_t *part, MPI_Comm *comm );
+
+int ParMETIS_V3_PartKway(
+   int64_t *vtxdist, int64_t *xadj, int64_t *adjncy, int64_t *vwgt,
+   int64_t *adjwgt, int64_t *wgtflag, int64_t *numflag, int64_t *ncon, int64_t *nparts,
+   double *tpwgts, double *ubvec, int64_t *options, int64_t *edgecut, int64_t *part,
+   MPI_Comm *comm );
+
+int ParMETIS_V3_PartGeomKway(
+   int64_t *vtxdist, int64_t *xadj, int64_t *adjncy, int64_t *vwgt,
+   int64_t *adjwgt, int64_t *wgtflag, int64_t *numflag, int64_t *ndims, double *xyz,
+   int64_t *ncon, int64_t *nparts, double *tpwgts, double *ubvec, int64_t *options,
+   int64_t *edgecut, int64_t *part, MPI_Comm *comm );
+
+int ParMETIS_V3_RefineKway(
+   int64_t *vtxdist, int64_t *xadj, int64_t *adjncy, int64_t *vwgt,
+   int64_t *adjwgt, int64_t *wgtflag, int64_t *numflag, int64_t *ncon, int64_t *nparts,
+   double *tpwgts, double *ubvec, int64_t *options, int64_t *edgecut,
+   int64_t *part, MPI_Comm *comm );
+
+} // namespace core
+} // namespace walberla
\ No newline at end of file
diff --git a/src/waLBerlaDefinitions.in.h b/src/waLBerlaDefinitions.in.h
index ce7d276a1..15a6df259 100644
--- a/src/waLBerlaDefinitions.in.h
+++ b/src/waLBerlaDefinitions.in.h
@@ -21,6 +21,7 @@
 // External libraries
 #cmakedefine WALBERLA_BUILD_WITH_MPI
 #cmakedefine WALBERLA_BUILD_WITH_METIS
+#cmakedefine WALBERLA_BUILD_WITH_PARMETIS
 
 #cmakedefine WALBERLA_BUILD_WITH_BOOST_THREAD
 #cmakedefine WALBERLA_BUILD_WITH_PYTHON
diff --git a/tests/core/CMakeLists.txt b/tests/core/CMakeLists.txt
index 1a10d3d01..214ec7d5f 100644
--- a/tests/core/CMakeLists.txt
+++ b/tests/core/CMakeLists.txt
@@ -200,3 +200,19 @@ waLBerla_execute_test( NAME UNIQUEID PROCESSES 4)
 
 waLBerla_compile_test( FILES VersionTest.cpp )
 waLBerla_execute_test( NAME VersionTest )
+
+##################
+# load_balancing #
+##################
+
+if( WALBERLA_BUILD_WITH_METIS )
+   waLBerla_compile_test( NAME MetisTest FILES load_balancing/MetisTest.cpp DEPENDS field )
+   waLBerla_execute_test( NAME MetisTest COMMAND $<TARGET_FILE:MetisTest> 64 64  4 --no-vtk )
+endif()
+
+if( WALBERLA_BUILD_WITH_PARMETIS )
+   waLBerla_compile_test( NAME ParMetisTest FILES load_balancing/ParMetisTest.cpp DEPENDS blockforest field stencil vtk  )
+   waLBerla_execute_test( NAME ParMetisTest1 COMMAND $<TARGET_FILE:ParMetisTest> 64 64  4 --no-vtk )
+   waLBerla_execute_test( NAME ParMetisTest2 COMMAND $<TARGET_FILE:ParMetisTest> 64 64  8 --no-vtk PROCESSES 2 )
+   waLBerla_execute_test( NAME ParMetisTest4 COMMAND $<TARGET_FILE:ParMetisTest> 64 64 16 --no-vtk PROCESSES 4 )
+endif()
\ No newline at end of file
diff --git a/tests/core/load_balancing/MetisTest.cpp b/tests/core/load_balancing/MetisTest.cpp
new file mode 100644
index 000000000..a69cffb99
--- /dev/null
+++ b/tests/core/load_balancing/MetisTest.cpp
@@ -0,0 +1,149 @@
+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can 
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT 
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
+//  for more details.
+//  
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file MetisTest.cpp
+//! \ingroup core
+//! \author Christian Godenschwager <christian.godenschwager@fau.de>
+//! \brief Test for MetisWrapper
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+#include "blockforest/communication/UniformBufferedScheme.h"
+
+#include "core/Environment.h"
+#include "core/debug/TestSubsystem.h"
+#include "core/load_balancing/MetisWrapper.h"
+#include "core/math/Vector2.h"
+#include "core/logging/Logging.h"
+
+#include "field/AddToStorage.h"
+#include "field/GhostLayerField.h"
+#include "field/vtk/VTKWriter.h"
+#include "field/communication/PackInfo.h"
+
+#include "stencil/D2Q9.h"
+
+#include "vtk/VTKOutput.h"
+#include "vtk/BlockCellDataWriter.h"
+
+#include <boost/lexical_cast.hpp>
+#include <boost/random.hpp>
+
+int main( int argc, char * argv[] )
+{
+   using namespace walberla;
+
+   mpi::Environment env( argc, argv );
+
+   debug::enterTestMode();
+
+   std::vector< std::string > args( argv, argv + argc );
+   Vector2<uint_t> fieldSize;
+   uint_t          partitions;
+   bool            vtk = true;
+
+   auto argIt = std::find( args.begin(), args.end(), "--no-vtk" );
+   if(argIt != args.end())
+   {
+      vtk = false;
+      args.erase( argIt );
+   }
+
+   try {
+      fieldSize.set( boost::lexical_cast< uint_t >( args.at(1) ), boost::lexical_cast< uint_t >( args.at(2) ) );
+      partitions = boost::lexical_cast< uint_t >( args.at(3) );
+   }
+   catch( std::exception & e )
+   {
+      WALBERLA_ABORT( "USAGE: " << args[0] << " SizeX SizeY #Partitions\n\n" << e.what() );
+   }
+
+   auto blocks = blockforest::createUniformBlockGrid( uint_t(1), uint_t(1), uint_t(1),
+                                                      fieldSize[0], fieldSize[1], uint_t(1),
+                                                      real_t(1),
+                                                      uint_t(1), uint_t(1), uint_t(1),
+                                                      true, true, false );
+
+   typedef field::GhostLayerField< int64_t, 1 > FieldType;
+
+   auto domainId    = field::addToStorage< FieldType >( blocks, "domain", int64_t(-1), field::zyxf, uint_t(1) );
+   auto partFieldId = field::addToStorage< FieldType >( blocks, "partitions", int64_t(-1), field::zyxf, uint_t(1) );
+
+   auto & domain    = *( blocks->begin()->getData< FieldType >( domainId    ) );
+   auto & partField = *( blocks->begin()->getData< FieldType >( partFieldId ) );
+
+   int64_t ctr = 0;
+   for(auto it = domain.begin(); it != domain.end(); ++it)
+   {
+      *it = ctr++;
+   }
+   WALBERLA_CHECK_GREATER( ctr, 0 );
+
+   blockforest::communication::UniformBufferedScheme< stencil::D2Q9 > scheme( blocks );
+   scheme.addPackInfo( make_shared< field::communication::PackInfo< FieldType > >( domainId ) );
+   scheme();
+
+   int64_t nvtxs = ctr;
+   int64_t ncon  = int64_t(1);
+   std::vector< int64_t > xadj, adjncy;
+   xadj.push_back( int64_t(0) );
+
+   Cell c( cell_idx_t(0), cell_idx_t(0), cell_idx_t(0) );
+   for( c[0] = 0; c[0] < cell_idx_c( domain.xSize() ); ++c[0] )
+      for( c[1] = 0; c[1] < cell_idx_c( domain.ySize() ); ++c[1] )
+      {
+         for( auto dirIt = stencil::D2Q9::beginNoCenter(); dirIt != stencil::D2Q9::end(); ++dirIt )
+            adjncy.push_back( domain.get( c + *dirIt ) );
+       
+         xadj.push_back( int64_c( adjncy.size() ) );         
+      }
+
+   WALBERLA_CHECK_EQUAL( xadj.size(), numeric_cast<size_t>( nvtxs ) + size_t(1) );
+
+   int64_t nparts = int64_c( partitions );
+   int64_t edgecut;
+   std::vector< int64_t > part( numeric_cast<size_t>( nvtxs ) );
+
+   int result =  core::METIS_PartGraphRecursive( &nvtxs, &ncon, &(xadj.front()), &(adjncy.front()), NULL, NULL, NULL,
+                                                 &nparts, NULL, NULL, NULL, &edgecut, &(part.front()) );
+
+   WALBERLA_CHECK_EQUAL( result, core::METIS_OK );  
+
+   c = Cell( cell_idx_t(0), cell_idx_t(0), cell_idx_t(0) );
+   auto it = part.begin();
+   for( c[0] = 0; c[0] < cell_idx_c( domain.xSize() ); ++c[0] )
+      for( c[1] = 0; c[1] < cell_idx_c( domain.ySize() ); ++c[1] )
+      {
+         if( domain.get( c ) == int64_t( -1 ) )
+            continue;
+
+         WALBERLA_CHECK_UNEQUAL( it, part.end() );
+
+         partField.get( c ) = *it++;
+
+      }
+   WALBERLA_CHECK_EQUAL( it, part.end() );
+
+   if( vtk )
+   {
+      auto vtkOutput = vtk::createVTKOutput_BlockData( blocks, "Metis", 1, 0 );
+      vtkOutput->addCellDataWriter( make_shared< field::VTKWriter< FieldType, int64_t > >( domainId, "domain" ) );
+      vtkOutput->addCellDataWriter( make_shared< field::VTKWriter< FieldType, int64_t > >( partFieldId, "partitions" ) );
+      vtk::writeFiles( vtkOutput )();
+   }
+
+   return EXIT_SUCCESS;
+}
diff --git a/tests/core/load_balancing/ParMetisTest.cpp b/tests/core/load_balancing/ParMetisTest.cpp
new file mode 100644
index 000000000..ce43ce97b
--- /dev/null
+++ b/tests/core/load_balancing/ParMetisTest.cpp
@@ -0,0 +1,170 @@
+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can 
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT 
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
+//  for more details.
+//  
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file ParMetisTest.cpp
+//! \ingroup core
+//! \author Christian Godenschwager <christian.godenschwager@fau.de>
+//! \brief Test for MetisWrapper
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+#include "blockforest/communication/UniformBufferedScheme.h"
+
+#include "core/Environment.h"
+#include "core/debug/TestSubsystem.h"
+#include "core/load_balancing/ParMetisWrapper.h"
+#include "core/math/Vector2.h"
+#include "core/math/IntegerFactorization.h"
+#include "core/logging/Logging.h"
+
+#include "field/AddToStorage.h"
+#include "field/GhostLayerField.h"
+#include "field/vtk/VTKWriter.h"
+#include "field/communication/PackInfo.h"
+
+#include "stencil/D2Q9.h"
+#include "stencil/D3Q27.h"
+
+#include "vtk/VTKOutput.h"
+#include "vtk/BlockCellDataWriter.h"
+
+#include <boost/lexical_cast.hpp>
+#include <boost/random.hpp>
+
+int main( int argc, char * argv[] )
+{
+   using namespace walberla;
+
+   mpi::Environment env( argc, argv );
+
+   debug::enterTestMode();
+
+   std::vector< std::string > args( argv, argv + argc );
+   Vector2<uint_t> fieldSize;
+   uint_t          partitions;
+   bool            vtk = true;
+
+   try {
+      fieldSize.set( boost::lexical_cast< uint_t >( args.at(1) ), boost::lexical_cast< uint_t >( args.at(2) ) );
+      partitions = boost::lexical_cast< uint_t >( args.at(3) );
+
+      auto it = std::find( args.begin(), args.end(), "--no-vtk" );
+      if(it != args.end())
+      {
+         vtk = false;
+         args.erase( it );
+      }
+   }
+   catch( std::exception & e )
+   {
+      WALBERLA_ABORT( "USAGE: " << args[0] << " SizeX SizeY #Partitions\n\n" << e.what() );
+   }
+
+   int numProcesses = MPIManager::instance()->numProcesses();
+   auto gridSize = math::getFactors( uint_c( numProcesses ), uint_t(2) );
+
+   auto blocks = blockforest::createUniformBlockGrid( gridSize[0],  gridSize[1], uint_t(1),
+                                                      fieldSize[0], fieldSize[1], uint_t(1),
+                                                      real_t(1),
+                                                      true,
+                                                      true, true, false );
+
+   typedef field::GhostLayerField< int64_t, 1 > FieldType;
+
+   auto domainId    = field::addToStorage< FieldType >( blocks, "domain", int64_t(-1), field::zyxf, uint_t(1) );
+   auto partFieldId = field::addToStorage< FieldType >( blocks, "partitions", int64_t(-1), field::zyxf, uint_t(1) );
+
+   auto & domain    = *( blocks->begin()->getData< FieldType >( domainId    ) );
+   auto & partField = *( blocks->begin()->getData< FieldType >( partFieldId ) );
+
+   WALBERLA_CHECK_EQUAL( std::distance( blocks->begin(), blocks->end() ), 1 );
+
+   int64_t ctr = int64_c( fieldSize[0] * fieldSize[1] ) * int64_c( MPIManager::instance()->rank() );
+   Cell c( cell_idx_t(0), cell_idx_t(0), cell_idx_t(0) );
+   for( c[0] = 0; c[0] < cell_idx_c( domain.xSize() ); ++c[0] )
+      for( c[1] = 0; c[1] < cell_idx_c( domain.ySize() ); ++c[1] )
+      {
+         domain.get( c ) = ctr++;
+      }
+
+   blockforest::communication::UniformBufferedScheme< stencil::D2Q9 > scheme( blocks );
+   scheme.addPackInfo( make_shared< field::communication::PackInfo< FieldType > >( domainId ) );
+   scheme();
+
+   std::vector< int64_t > vtxdist;
+   for( int i = 0; i < MPIManager::instance()->numProcesses() + 1; ++i )
+   {
+      vtxdist.push_back( int64_c(i) * int64_c( fieldSize[0] * fieldSize[1] ) );
+   }
+
+   int64_t ncon  = int64_t(1);
+   std::vector< int64_t > xadj, adjncy;
+   xadj.push_back( int64_t(0) );
+
+   c = Cell( cell_idx_t(0), cell_idx_t(0), cell_idx_t(0) );
+   for( c[0] = 0; c[0] < cell_idx_c( domain.xSize() ); ++c[0] )
+      for( c[1] = 0; c[1] < cell_idx_c( domain.ySize() ); ++c[1] )
+      {
+         for( auto dirIt = stencil::D2Q9::beginNoCenter(); dirIt != stencil::D2Q9::end(); ++dirIt )
+            adjncy.push_back( domain.get( c + *dirIt ) );
+       
+         xadj.push_back( int64_c( adjncy.size() ) );         
+      }
+
+   WALBERLA_CHECK_EQUAL( xadj.size(), fieldSize[0] * fieldSize[1] + uint_t(1) );
+   WALBERLA_CHECK_EQUAL( adjncy.size(), fieldSize[0] * fieldSize[1] * uint_t(8) );
+
+   int64_t wgtflag = 0;
+   int64_t numflag = 0;
+   int64_t nparts = int64_c( partitions );
+   std::vector< double > tpwgts( partitions, 1.0 / numeric_cast<double>( partitions ) );
+   std::vector< double > ubvec( ncon, 1.05 );
+   int64_t options[] = {0,0,0};
+   int64_t edgecut;
+   std::vector< int64_t > part( fieldSize[0] * fieldSize[1] );
+   MPI_Comm comm = MPIManager::instance()->comm();
+
+   int result = core::ParMETIS_V3_PartKway( &(vtxdist.front()), &(xadj.front()), &(adjncy.front()), NULL, NULL, &wgtflag, &numflag, &ncon, &nparts,
+                                            &(tpwgts.front()), &(ubvec.front()), options, &edgecut, &(part.front()), &comm );
+
+
+   WALBERLA_CHECK_EQUAL( result, core::METIS_OK );  
+
+   c = Cell( cell_idx_t(0), cell_idx_t(0), cell_idx_t(0) );
+   auto it = part.begin();
+   for( c[0] = 0; c[0] < cell_idx_c( domain.xSize() ); ++c[0] )
+      for( c[1] = 0; c[1] < cell_idx_c( domain.ySize() ); ++c[1] )
+      {
+         if( domain.get( c ) == int64_t( -1 ) )
+            continue;
+
+         WALBERLA_CHECK_UNEQUAL( it, part.end() );
+
+         partField.get( c ) = *it++;
+
+      }
+   WALBERLA_CHECK_EQUAL( it, part.end() );
+
+   if( vtk )
+   {
+      auto vtkOutput = vtk::createVTKOutput_BlockData( blocks, "Metis", 1, 0 );
+      vtkOutput->addCellDataWriter( make_shared< field::VTKWriter< FieldType, int64_t > >( domainId, "domain" ) );
+      vtkOutput->addCellDataWriter( make_shared< field::VTKWriter< FieldType, int64_t > >( partFieldId, "partitions" ) );
+      vtk::writeFiles( vtkOutput )();
+   }
+
+   return EXIT_SUCCESS;
+}
-- 
GitLab