diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp
index 50ad0995f..5336c7685 100644
--- a/opm/simulators/flow/BlackoilModelEbos.hpp
+++ b/opm/simulators/flow/BlackoilModelEbos.hpp
@@ -266,7 +266,11 @@ namespace Opm {
         {
             // Create partitions.
             const auto& [partition_vector, num_domains] =
-                this->partitionCells();
+                partitionCells(this->grid_,
+                               this->ebosSimulator_.vanguard().schedule().getWellsatEnd(),
+                               this->param_.local_domain_partition_method_,
+                               this->param_.num_local_domains_,
+                               this->param_.local_domain_partition_imbalance_);
 
             // Scan through partitioning to get correct size for each.
             std::vector<int> sizes(num_domains, 0);
@@ -1875,19 +1879,6 @@ namespace Opm {
         }
 
     private:
-        template <typename, class = void>
-        struct HasZoltanPartitioning : public std::false_type {};
-
-        template <typename GridType>
-        struct HasZoltanPartitioning<
-            GridType,
-            std::void_t<decltype(std::declval<GridType>().zoltanPartitionWithoutScatter
-                                 (std::declval<const std::vector<Well>*>(),
-                                  std::declval<const double*>(),
-                                  std::declval<const int>(),
-                                  std::declval<const double>()))>
-            > : public std::true_type {};
-
         template<class T>
         bool isNumericalAquiferCell(const Dune::CpGrid& grid, const T& elem)
         {
@@ -1914,51 +1905,6 @@ namespace Opm {
         double maxResidualAllowed() const { return param_.max_residual_allowed_; }
         double linear_solve_setup_time_;
 
-        std::pair<std::vector<int>, int> partitionCells() const
-        {
-            const std::string& method = this->param_.local_domain_partition_method_;
-            if (method == "zoltan") {
-                if constexpr (HasZoltanPartitioning<Grid>::value) {
-                    return this->partitionCellsZoltan();
-                } else {
-                    OPM_THROW(std::runtime_error, "Zoltan requested for local domain partitioning, "
-                              "but is not available for the current grid type.");
-                }
-            } else if (method == "simple") {
-                const int num_cells = detail::countLocalInteriorCells(this->grid_);
-                return partitionCellsSimple(num_cells, this->param_.num_local_domains_);
-            } else if (method.size() > 10 && method.substr(method.size() - 10, 10) == ".partition") {
-                // Method string ends with ".partition", interpret as filename for partitioning.
-                const int num_cells = detail::countLocalInteriorCells(this->grid_);
-                return partitionCellsFromFile(method, num_cells);
-            } else {
-                OPM_THROW(std::runtime_error, "Unknown local domain partitioning method requested: " + method);
-            }
-        }
-
-        std::pair<std::vector<int>, int> partitionCellsZoltan() const
-        {
-            const auto wells = this->ebosSimulator_.vanguard().schedule().getWellsatEnd();
-
-            auto partition_vector = this->grid_.zoltanPartitionWithoutScatter
-                (&wells, nullptr, this->param_.num_local_domains_,
-                 this->param_.local_domain_partition_imbalance_);
-
-            return this->countDomains(std::move(partition_vector));
-        }
-
-        std::pair<std::vector<int>, int>
-        countDomains(std::vector<int> partition_vector) const
-        {
-            auto maxPos = std::max_element(partition_vector.begin(),
-                                           partition_vector.end());
-
-            const auto num_domains = (maxPos == partition_vector.end())
-                ? 0 : *maxPos + 1;
-
-            return { std::move(partition_vector), num_domains };
-        }
-
     public:
         std::vector<bool> wasSwitched_;
     };
diff --git a/opm/simulators/flow/partitionCells.cpp b/opm/simulators/flow/partitionCells.cpp
index 7a39b5d7e..f3e2d71aa 100644
--- a/opm/simulators/flow/partitionCells.cpp
+++ b/opm/simulators/flow/partitionCells.cpp
@@ -17,8 +17,20 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
+#include <config.h>
 #include <opm/simulators/flow/partitionCells.hpp>
 
+#include <opm/grid/CpGrid.hpp>
+#include <opm/grid/polyhedralgrid.hh>
+
+#include <opm/input/eclipse/Schedule/Well/Well.hpp>
+
+#include <opm/simulators/flow/countGlobalCells.hpp>
+
+#if HAVE_DUNE_ALUGRID
+#include <dune/alugrid/grid.hh>
+#endif // HAVE_DUNE_ALUGRID
+
 #include <fmt/format.h>
 
 #include <algorithm>
@@ -28,47 +40,143 @@
 #include <numeric>
 #include <stdexcept>
 #include <string>
+#include <type_traits>
 
-namespace Opm
+namespace {
+
+std::pair<std::vector<int>, int>
+countDomains(std::vector<int> partition_vector)
 {
+    auto maxPos = std::max_element(partition_vector.begin(),
+                                   partition_vector.end());
 
+    const auto num_domains = (maxPos == partition_vector.end())
+        ? 0 : *maxPos + 1;
 
-    // Read from file, containing one number per cell.
-    std::pair<std::vector<int>, int> partitionCellsFromFile(const std::string& filename, const int num_cells)
-    {
-        // Read file into single vector.
-        std::ifstream is(filename);
-        const std::vector<int> cellpart{std::istream_iterator<int>(is), std::istream_iterator<int>()};
-        if (cellpart.size() != size_t(num_cells)) {
-            auto msg = fmt::format("Partition file contains {} entries, but there are {} cells.",
-                                   cellpart.size(), num_cells);
-            throw std::runtime_error(msg);
+    return { std::move(partition_vector), num_domains };
+}
+
+template <typename, class = void>
+struct HasZoltanPartitioning : public std::false_type {};
+
+template <typename GridType>
+struct HasZoltanPartitioning<
+    GridType,
+    std::void_t<decltype(std::declval<const GridType&>().zoltanPartitionWithoutScatter
+                         (std::declval<const std::vector<Opm::Well>*>(),
+                          std::declval<const double*>(),
+                          std::declval<const int>(),
+                          std::declval<const double>()))>
+    > : public std::true_type {};
+
+} // anonymous namespace
+
+namespace Opm {
+
+template<class Grid>
+std::pair<std::vector<int>, int> partitionCells(const Grid& grid,
+                                                const std::vector<Well>& wells,
+                                                const std::string& method,
+                                                const int num_local_domains,
+                                                const double partition_imbalance)
+{
+    if (method == "zoltan") {
+        if constexpr (HasZoltanPartitioning<Grid>::value) {
+            return partitionCellsZoltan(grid, wells, num_local_domains, partition_imbalance);
+        } else {
+            OPM_THROW(std::runtime_error, "Zoltan requested for local domain partitioning, "
+                      "but is not available for the current grid type.");
         }
+    } else if (method == "simple") {
+        const int num_cells = detail::countLocalInteriorCells(grid);
+        return partitionCellsSimple(num_cells, num_local_domains);
+    } else if (method.size() > 10 && method.substr(method.size() - 10, 10) == ".partition") {
+        // Method string ends with ".partition", interpret as filename for partitioning.
+        const int num_cells = detail::countLocalInteriorCells(grid);
+        return partitionCellsFromFile(method, num_cells);
+    } else {
+        OPM_THROW(std::runtime_error, "Unknown local domain partitioning method requested: " + method);
+    }
+}
 
-        // Create and return the output domain vector.
-        const int num_domains = (*std::max_element(cellpart.begin(), cellpart.end())) + 1;
-        return { cellpart, num_domains };
+
+// Read from file, containing one number per cell.
+std::pair<std::vector<int>, int> partitionCellsFromFile(const std::string& filename, const int num_cells)
+{
+    // Read file into single vector.
+    std::ifstream is(filename);
+    const std::vector<int> cellpart{std::istream_iterator<int>(is), std::istream_iterator<int>()};
+    if (cellpart.size() != size_t(num_cells)) {
+        auto msg = fmt::format("Partition file contains {} entries, but there are {} cells.",
+                               cellpart.size(), num_cells);
+        throw std::runtime_error(msg);
     }
 
+    // Create and return the output domain vector.
+    const int num_domains = (*std::max_element(cellpart.begin(), cellpart.end())) + 1;
+    return { cellpart, num_domains };
+}
 
-    // Trivially simple partitioner
-    std::pair<std::vector<int>, int> partitionCellsSimple(const int num_cells, const int num_domains)
-    {
-        // Build the partitions.
-        const int dom_sz = num_cells / num_domains;
-        std::vector<int> bounds(num_domains + 1, dom_sz);
-        bounds[0] = 0;
-        for (int i = 0; i < num_cells % num_domains; ++i) {
-            ++bounds[i + 1];
-        }
-        std::partial_sum(bounds.begin(), bounds.end(), bounds.begin());
-        assert(bounds[num_domains] == num_cells);
-        std::vector<int> part(num_cells);
-        for (int i = 0; i < num_domains; ++i) {
-            std::fill(part.begin() + bounds[i], part.begin() + bounds[i + 1], i);
-        }
-        return { part, num_domains };
+
+// Trivially simple partitioner
+std::pair<std::vector<int>, int> partitionCellsSimple(const int num_cells, const int num_domains)
+{
+    // Build the partitions.
+    const int dom_sz = num_cells / num_domains;
+    std::vector<int> bounds(num_domains + 1, dom_sz);
+    bounds[0] = 0;
+    for (int i = 0; i < num_cells % num_domains; ++i) {
+        ++bounds[i + 1];
     }
+    std::partial_sum(bounds.begin(), bounds.end(), bounds.begin());
+    assert(bounds[num_domains] == num_cells);
+    std::vector<int> part(num_cells);
+    for (int i = 0; i < num_domains; ++i) {
+        std::fill(part.begin() + bounds[i], part.begin() + bounds[i + 1], i);
+    }
+    return { part, num_domains };
+}
 
 
+template<class Grid>
+std::pair<std::vector<int>, int> partitionCellsZoltan(const Grid& grid,
+                                                      const std::vector<Well>& wells,
+                                                      const int num_domains,
+                                                      const double domain_imbalance)
+{
+    auto partition_vector = grid.zoltanPartitionWithoutScatter
+        (&wells, nullptr, num_domains, domain_imbalance);
+
+    return countDomains(std::move(partition_vector));
+}
+
+template std::pair<std::vector<int>,int>
+partitionCells<Dune::CpGrid>(const Dune::CpGrid&,
+                             const std::vector<Well>&,
+                             const std::string&,
+                             const int,
+                             const double);
+
+template std::pair<std::vector<int>,int>
+partitionCells<Dune::PolyhedralGrid<3,3,double>>(const Dune::PolyhedralGrid<3,3,double>&,
+                                                 const std::vector<Well>&,
+                                                 const std::string&,
+                                                 const int,
+                                                 const double);
+
+#if HAVE_DUNE_ALUGRID
+#if HAVE_MPI
+using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
+#else
+using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
+#endif //HAVE_MPI
+
+template std::pair<std::vector<int>,int>
+partitionCells<ALUGrid3CN>(const ALUGrid3CN&,
+                           const std::vector<Well>&,
+                           const std::string&,
+                           const int,
+                           const double);
+#endif
+
 } // namespace Opm
diff --git a/opm/simulators/flow/partitionCells.hpp b/opm/simulators/flow/partitionCells.hpp
index b39077efc..5dfd9b84a 100644
--- a/opm/simulators/flow/partitionCells.hpp
+++ b/opm/simulators/flow/partitionCells.hpp
@@ -24,18 +24,34 @@
 #include <utility>
 #include <vector>
 
-namespace Opm
-{
+namespace Opm {
 
-    /// Read a partitioning from file, assumed to contain one number per cell, its partition number.
-    /// \return pair containing a partition vector (partition number for each cell), and the number of partitions.
-    std::pair<std::vector<int>, int> partitionCellsFromFile(const std::string& filename, const int num_cells);
+class Well;
 
-    /// Trivially simple partitioner assigning partitions en bloc, consecutively by cell index.
-    /// \return pair containing a partition vector (partition number for each cell), and the number of partitions.
-    std::pair<std::vector<int>, int> partitionCellsSimple(const int num_cells, const int num_domains);
+/// Partitions the grid using the specified method.
+/// \return pair containing a partition vector (partition number for each cell), and the number of partitions.
+template<class Grid>
+std::pair<std::vector<int>, int> partitionCells(const Grid& grid,
+                                                const std::vector<Well>& wells,
+                                                const std::string& method,
+                                                const int num_local_domains,
+                                                const double partition_imbalance);
 
+/// Read a partitioning from file, assumed to contain one number per cell, its partition number.
+/// \return pair containing a partition vector (partition number for each cell), and the number of partitions.
+std::pair<std::vector<int>, int> partitionCellsFromFile(const std::string& filename, const int num_cells);
+
+/// Trivially simple partitioner assigning partitions en bloc, consecutively by cell index.
+/// \return pair containing a partition vector (partition number for each cell), and the number of partitions.
+std::pair<std::vector<int>, int> partitionCellsSimple(const int num_cells, const int num_domains);
+
+/// Partitions the grid using the Zoltan graph partitioner.
+/// \return pair containing a partition vector (partition number for each cell), and the number of partitions.
+template<class Grid>
+std::pair<std::vector<int>, int> partitionCellsZoltan(const Grid& grid,
+                                                      const std::vector<Well>& wells,
+                                                      const int num_domains,
+                                                      const double domain_imbalance);
 } // namespace Opm
 
-
 #endif // OPM_ASPINPARTITION_HEADER_INCLUDED