diff --git a/Applications/Utils/ModelPreparation/PartitionMesh/CMakeLists.txt b/Applications/Utils/ModelPreparation/PartitionMesh/CMakeLists.txt
index f1c3a152bd74d64092223744e2f815a7ee35ab7a..83047c47c8a9cbab5bc07ad8fd1644c306a7fb48 100644
--- a/Applications/Utils/ModelPreparation/PartitionMesh/CMakeLists.txt
+++ b/Applications/Utils/ModelPreparation/PartitionMesh/CMakeLists.txt
@@ -1,5 +1,5 @@
 ###
-add_executable(partmesh PartitionMesh.cpp MeshPartitioning.h MeshPartitioning.cpp)
+add_executable(partmesh PartitionMesh.cpp NodeWiseMeshPartitioner.h NodeWiseMeshPartitioner.cpp)
 set_target_properties(partmesh PROPERTIES FOLDER Utilities)
 target_link_libraries(partmesh MeshLib)
 
diff --git a/Applications/Utils/ModelPreparation/PartitionMesh/MeshPartitioning.cpp b/Applications/Utils/ModelPreparation/PartitionMesh/MeshPartitioning.cpp
deleted file mode 100644
index f3719c3b4d2a9d6a45a8a769548caebf8f0e1de6..0000000000000000000000000000000000000000
--- a/Applications/Utils/ModelPreparation/PartitionMesh/MeshPartitioning.cpp
+++ /dev/null
@@ -1,493 +0,0 @@
-/*!
-  \file MeshPartitioning.cpp
-  \date   2016.05
-
-  \brief  Define the members of class MeshPartitioning
-
-  \copyright
-  Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
-             Distributed under a Modified BSD License.
-               See accompanying file LICENSE.txt or
-               http://www.opengeosys.org/project/license
-
-*/
-
-#include "MeshPartitioning.h"
-
-#include <iomanip>
-#include <stdio.h> // for binary output
-
-#include "logog/include/logog.hpp"
-
-#include "MeshLib/Node.h"
-#include "MeshLib/Elements/Element.h"
-
-namespace MeshLib
-{
-void MeshPartitioning :: write2METIS(const std::string& file_name)
-{
-    std::ofstream os(file_name.data(), std::ios::trunc);
-    os << _elements.size() <<" \n";
-    for (const auto elem : _elements)
-    {
-        for (unsigned j=0; j<elem->getNumberOfNodes(); j++)
-        {
-            os << elem->getNodeIndex(j) + 1 <<" ";
-        }
-        os << std::endl;
-    }
-}
-
-void MeshPartitioning :: partitionByNodeMETIS(const std::string& file_name_base,
-        const unsigned npartitions,
-        const bool output_binary)
-{
-    const std::string npartitions_str = std::to_string(npartitions);
-
-    // -- Read partitioned mesh data from METIS
-    std::string fname_parts = file_name_base + ".mesh.npart." + npartitions_str;
-
-    std::ifstream npart_in(fname_parts.data());
-    if (!npart_in.is_open())
-    {
-        ERR("Error: cannot open file %s. It may not exist! ", fname_parts.data());
-        abort();
-    }
-
-    const std::size_t nnodes = _nodes.size();
-
-    std::vector<std::size_t> nodes_partition_ids;
-    nodes_partition_ids.reserve(nnodes);
-
-    for (std::size_t i=0; i<nnodes; i++)
-    {
-        std::size_t part_id;
-        npart_in >> part_id >> std::ws;
-        nodes_partition_ids.push_back(part_id);
-    }
-    npart_in.close();
-    //TEST  std::remove(fname_parts.c_str());
-
-    // -- Open files for output
-    //
-    // for ASCII output
-    std::fstream os_subd;
-    std::fstream os_subd_head;
-    std::fstream os_subd_node;
-    // for binary output
-    FILE *of_bin_cfg = 0;   // File to output headers
-    FILE *of_bin_nod = 0;   // File to output an array of all nodes
-    FILE *of_bin_ele = 0;   // File to output an array of all non-ghost elements
-    FILE *of_bin_ele_g = 0; // File to output an array of all ghost elements
-
-    if (output_binary)
-    {
-        std::string fname = file_name_base + "_partitioned_msh_cfg" + npartitions_str + ".bin";
-        of_bin_cfg = fopen(fname.c_str(), "wb");
-
-        fname = file_name_base + "_partitioned_msh_ele" + npartitions_str + ".bin";
-        of_bin_ele = fopen(fname.c_str(), "wb");
-
-        fname = file_name_base + "_partitioned_msh_ele_g" + npartitions_str + ".bin";
-        of_bin_ele_g = fopen(fname.c_str(), "wb");
-
-        fname = file_name_base + "_partitioned_msh_nod" + npartitions_str + ".bin";
-        of_bin_nod = fopen(fname.c_str(), "wb");
-    }
-    else
-    {
-        std::string fname = file_name_base + "_partitioned_cfg" + npartitions_str + ".msh";
-        os_subd_head.open(fname.c_str(), std::ios::out|std::ios::trunc );
-        const std::string mesh_info = "Subdomain mesh "
-                                      "(Number of non-ghost nodes;  Number of non ghost base nodes; Number of elements; "
-                                      "Number of ghost elements; Number of base nodes; Number of nodes) "
-                                      "Number of the base nodes of the global mesh; Number of the nodes of the global mesh; "
-                                      "Number of integers to define elements; Number of integers to define ghost elements; "
-                                      "Reserved number.";
-        os_subd_head << mesh_info << std::endl;
-        os_subd_head << npartitions << std::endl;
-
-        fname = file_name_base + "_partitioned_elems_" + npartitions_str + ".msh";
-        os_subd.open(fname.c_str(), std::ios::out|std::ios::trunc );
-
-        fname = file_name_base + "_partitioned_nodes_" + npartitions_str + ".msh";
-        os_subd_node.open(fname.c_str(), std::ios::out|std::ios::trunc );
-
-        std::setw(14);
-        os_subd.precision(14);
-        os_subd.setf(std::ios::scientific);
-    }
-
-    PetscInt global_node_id_offset = 0;
-    std::vector<bool> nodes_part_reserved;
-    nodes_part_reserved.reserve(nnodes);
-    std::vector<PetscInt> node_global_ids;
-    node_global_ids.reserve(nnodes);
-    for (std::size_t i=0; i<nnodes; i++)
-    {
-        nodes_part_reserved.push_back(false);
-        node_global_ids.push_back(0);
-    }
-    std::vector<bool> nodes_reseved(nnodes);
-    std::vector<bool> elems_reseved(_elements.size());
-
-    std::vector<MeshLib::Node*> partition_nodes;  // non-ghost nodes of a partition
-    std::vector<std::size_t> partition_start_node_id(npartitions);
-    std::vector<std::size_t> partition_end_node_id(npartitions);
-    std::vector<std::size_t> end_non_ghost_node_id(npartitions);
-
-    const PetscInt num_headers = 14;
-    PetscInt head[num_headers];
-    // node rank offset
-    head[10] = 0;
-    // element rank offset
-    head[11] = 0;
-    // ghost element rank offset
-    head[12] = 0;
-    // Reserved
-    head[13] = 0;
-
-    for (std::size_t ipart=0; ipart<npartitions; ipart++)
-    {
-        INFO("Processing partition: %d", ipart);
-
-        partition_start_node_id[ipart] = partition_nodes.size();
-
-        std::fill(nodes_reseved.begin(), nodes_reseved.end(), false);
-
-        // Find non-ghost nodes in this partition
-        for (std::size_t i=0; i< nnodes; i++)
-        {
-            if ( (nodes_partition_ids[i] == ipart) && (!nodes_part_reserved[i]) )
-            {
-                partition_nodes.push_back(_nodes[i]);
-                nodes_part_reserved[i] = true;
-                nodes_reseved[i] = true;
-            }
-        }
-        end_non_ghost_node_id[ipart] = partition_nodes.size();
-
-        // Find elements in this partition
-        std::vector<const Element*> partition_regular_elements;
-        std::vector<const Element*> partition_ghost_elements;
-        std::fill(elems_reseved.begin(), elems_reseved.end(), false);
-        // For ghost element
-        std::vector< std::vector<unsigned> >
-        elems_non_ghost_nodes_local_ids(_elements.size());
-
-        const PetscInt num_non_ghost_nodes =  end_non_ghost_node_id[ipart]
-                                           - partition_start_node_id[ipart];
-        for (PetscInt i=0; i<num_non_ghost_nodes; i++)
-        {
-            const unsigned node_id = partition_nodes[ i +  partition_start_node_id[ipart] ]->getID();
-
-            // Search the elements connected to this nodes
-            std::vector<Element*> const& connected_elems = _nodes[node_id]->getElements();
-
-            for (const auto elem : connected_elems)
-            {
-                // If checked
-                if (elems_reseved[elem->getID()])
-                    continue;
-
-                std::vector<unsigned> nonghost_nodes_local_ids;
-                unsigned ghost_counter = 0;
-                for (unsigned kk=0; kk<elem->getNumberOfNodes(); kk++)
-                {
-                    if (nodes_reseved[elem->getNodeIndex(kk)])
-                    {
-                        nonghost_nodes_local_ids.push_back(kk);
-                    }
-                    else
-                    {
-                        ghost_counter++;
-                    }
-                }
-
-                // All nodes of this element are inside this partition
-                if (ghost_counter == 0)
-                {
-                    partition_regular_elements.push_back(elem);
-                }
-                else // ghost element
-                {
-                    partition_ghost_elements.push_back(elem);
-                    elems_non_ghost_nodes_local_ids[elem->getID()]
-                        = std::move(nonghost_nodes_local_ids);
-                }
-
-                elems_reseved[elem->getID()] = true;
-            }
-        }
-
-        // Add ghost nodes in ghost elements to the node vector of this partition
-        // nodes_reseved is reused without cleaning
-        // Mark the non-ghost nodes for each ghost element
-        for (const auto ghost_elem : partition_ghost_elements)
-        {
-            for (unsigned k=0; k<ghost_elem->getNumberOfNodes(); k++)
-                nodes_reseved[ghost_elem->getNodeIndex(k)] = false;
-
-            // Mark non-ghost nodes
-            std::vector<unsigned>& local_node_ids = elems_non_ghost_nodes_local_ids[ghost_elem->getID()];
-            for (std::size_t k=0; k<local_node_ids.size(); k++)
-                nodes_reseved[ghost_elem->getNodeIndex(local_node_ids[k])] = true;
-        }
-        // Add the ghost nodes to the node vector of this partition
-        for (const auto ghost_elem : partition_ghost_elements)
-        {
-            for (unsigned k=0; k<ghost_elem->getNumberOfNodes(); k++)
-            {
-                const unsigned node_id = ghost_elem->getNodeIndex(k);
-                if (nodes_reseved[node_id])
-                    continue;
-                nodes_reseved[node_id] = true;
-                partition_nodes.push_back(ghost_elem->getNode(k));
-            }
-        }
-
-        partition_end_node_id[ipart] = partition_nodes.size();
-        // Renumber
-        PetscInt local_id = 0;
-
-        std::vector<unsigned> nodes_local_ids_partition;
-        nodes_local_ids_partition.reserve(nnodes);
-        for (std::size_t i=0; i<nnodes; i++)
-        {
-            nodes_local_ids_partition.push_back(0);
-        }
-
-        for (std::size_t i = partition_start_node_id[ipart];
-                i < end_non_ghost_node_id[ipart]; i++)
-        {
-            const unsigned node_id = partition_nodes[i]->getID();
-            node_global_ids[node_id] = global_node_id_offset;
-            nodes_local_ids_partition[node_id] = local_id;
-            local_id++;
-            global_node_id_offset++;
-        }
-        // Assign local IDs to ghost nodes
-        for (std::size_t i = end_non_ghost_node_id[ipart];
-                i < partition_end_node_id[ipart]; i++)
-        {
-            nodes_local_ids_partition[partition_nodes[i]->getID()] = local_id;
-            local_id++;
-        }
-
-        // Count the number of integer variables used to define non-ghost elements
-        const PetscInt num_non_ghost_elems = partition_regular_elements.size();
-        PetscInt nmb_element_idxs = 3 * num_non_ghost_elems;
-        for (PetscInt j=0; j<num_non_ghost_elems; j++)
-        {
-            nmb_element_idxs += partition_regular_elements[j]->getNumberOfNodes();
-        }
-        const PetscInt num_ghost_elems = partition_ghost_elements.size();
-        PetscInt nmb_element_idxs_g = 3 * num_ghost_elems;
-        for (PetscInt j=0; j<num_ghost_elems; j++)
-        {
-            nmb_element_idxs_g += partition_ghost_elements[j]->getNumberOfNodes();
-        }
-
-        std::string ipart_str = std::to_string(ipart);
-
-        // Number of active elements
-        const PetscInt offset_e = num_non_ghost_elems + nmb_element_idxs;
-        const PetscInt offset_e_g = num_ghost_elems + nmb_element_idxs_g;
-
-        // Reserved for nodes for high order interpolation
-        const PetscInt num_non_ghost_nodes_extra = 0;
-
-        const PetscInt num_nodes_linear_element =   partition_end_node_id[ipart]
-                - partition_start_node_id[ipart];
-        const PetscInt num_nodes_quadratic_element = 0;
-
-        const PetscInt num_nodes_global_mesh = nnodes;
-        // Reserved for nodes for high order interpolation
-        const PetscInt num_nodes_global_mesh_extra = nnodes;
-
-        // Write configuration data and elements of this partition
-        head[0] = num_nodes_linear_element + num_nodes_quadratic_element;
-        head[1] = num_nodes_linear_element;
-        head[2] = num_non_ghost_elems;
-        head[3] = num_ghost_elems;
-        head[4] = num_non_ghost_nodes;
-        head[5] = num_non_ghost_nodes + num_non_ghost_nodes_extra; // active nodes
-        head[6] = num_nodes_global_mesh;
-        head[7] = num_nodes_global_mesh_extra;
-        head[8] = nmb_element_idxs;
-        head[9] = nmb_element_idxs_g;
-        if (output_binary)
-        {
-            fwrite(head, 1, num_headers * sizeof(PetscInt), of_bin_cfg);
-
-            head[10] += head[0] * sizeof(NodeStruct);
-            head[11] += offset_e * sizeof(PetscInt);
-            head[12] += offset_e_g * sizeof(PetscInt);
-
-            // Write non-ghost elements
-            std::vector<PetscInt> ele_info(offset_e);
-
-            PetscInt counter = num_non_ghost_elems;
-
-            for (PetscInt j = 0; j < num_non_ghost_elems; j++)
-            {
-                ele_info[j] = counter;
-                getElementIntegerVariables(*partition_regular_elements[j], nodes_local_ids_partition,
-                                           ele_info, counter);
-            }
-            fwrite(&ele_info[0], 1, (offset_e) * sizeof(PetscInt), of_bin_ele);
-
-            // Write ghost elements
-            ele_info.resize(offset_e_g);
-            counter = num_ghost_elems;
-            for (PetscInt j = 0; j < num_ghost_elems; j++)
-            {
-                ele_info[j] = counter;
-                getElementIntegerVariables(*partition_ghost_elements[j], nodes_local_ids_partition,
-                                           ele_info, counter);
-            }
-            fwrite(&ele_info[0], 1, (offset_e_g) * sizeof(PetscInt), of_bin_ele_g);
-        }
-        else
-        {
-
-            for (int j=0; j<10; j++)
-            {
-                os_subd_head << head[j] << " ";
-            }
-            os_subd_head << " 0 " << std::endl;
-
-            for (const auto elem : partition_regular_elements)
-            {
-                writeLocalElementNodeIndicies(os_subd, *elem, nodes_local_ids_partition);
-            }
-
-            for (const auto elem : partition_ghost_elements)
-            {
-                writeLocalElementNodeIndicies(os_subd, *elem, nodes_local_ids_partition);
-            }
-            os_subd << std::endl;
-        }
-    } // end of for(unsigned ipart=0; ipart<npartitions; ipart++) for partitioning of nodes, elements
-
-    nodes_partition_ids.clear();
-    nodes_part_reserved.clear();
-    nodes_reseved.clear();
-    elems_reseved.clear();
-
-    // Write nodes
-    if (output_binary)
-    {
-        fclose(of_bin_cfg);
-        fclose(of_bin_ele);
-        fclose(of_bin_ele_g);
-
-        for (std::size_t ipart=0; ipart<npartitions; ipart++)
-        {
-            const size_t nnodes_part = partition_end_node_id[ipart] - partition_start_node_id[ipart];
-            std::vector<NodeStruct> nodes_buffer;
-            nodes_buffer.reserve(nnodes_part);
-
-            for (std::size_t i=partition_start_node_id[ipart]; i<partition_end_node_id[ipart]; i++)
-            {
-                double const* coords = partition_nodes[i]->getCoords();
-                NodeStruct node_struct;
-                node_struct.id = node_global_ids[partition_nodes[i]->getID()];
-                node_struct.x = coords[0];
-                node_struct.y = coords[1];
-                node_struct.z = coords[2];
-                nodes_buffer.emplace_back(node_struct);
-            }
-            fwrite(&nodes_buffer[0], sizeof(NodeStruct), nnodes_part, of_bin_nod);
-        }
-        fclose(of_bin_nod);
-    }
-    else
-    {
-        os_subd_head.close();
-        os_subd.close();
-        std::setw(14);
-        os_subd_node.precision(14);
-        os_subd_node.setf(std::ios::scientific);
-
-        for (std::size_t ipart=0; ipart<npartitions; ipart++)
-        {
-            for (std::size_t i=partition_start_node_id[ipart]; i<partition_end_node_id[ipart]; i++)
-            {
-                double const* coords = partition_nodes[i]->getCoords();
-                os_subd_node << node_global_ids[partition_nodes[i]->getID()] << " "
-                             << coords[0] << " " << coords[1] << " " << coords[2]<<"\n";
-            }
-            os_subd_node << std::endl;
-        }
-        os_subd_node.close();
-    }
-
-    // Update the node IDs
-    for (std::size_t i=0; i<nnodes; i++)
-    {
-        _nodes[i]->setID(node_global_ids[i]);
-    }
-    // sort
-    std::sort(_nodes.begin(), _nodes.end(),
-              [](const MeshLib::Node* a, const MeshLib::Node* b)
-                {
-                   return a->getID() < b->getID();
-                }
-             );
-}
-
-ElementType MeshPartitioning::getElementType(const Element& elem)
-{
-    switch ( elem.getCellType() )
-    {
-        case MeshLib::CellType::LINE2:
-            return ElementType::LINE2;
-        case MeshLib::CellType::QUAD4:
-            return ElementType::QUAD4;
-        case MeshLib::CellType::HEX8:
-            return ElementType::HEX8;
-        case MeshLib::CellType::TRI3:
-            return ElementType::TRI3;
-        case MeshLib::CellType::PYRAMID5:
-            return ElementType::PYRAMID5;
-        case MeshLib::CellType::PRISM6:
-            return ElementType::PRISM6;
-        default:
-            ERR("Invalid element type in element %d", elem.getID());
-            abort();
-    }
-}
-
-void MeshPartitioning::getElementIntegerVariables(const Element& elem,
-        const std::vector<unsigned>& local_node_ids,
-        std::vector<PetscInt>& elem_info,
-        PetscInt& counter)
-{
-    unsigned mat_id = 0; // Materical ID to be set from the mesh data
-    const PetscInt nn = elem.getNumberOfNodes();;
-    elem_info[counter++] = mat_id;
-    elem_info[counter++] = static_cast<unsigned>(getElementType(elem)) + 1;
-    elem_info[counter++] = nn;
-
-    for(PetscInt i=0; i<nn; i++)
-    {
-        elem_info[counter++] = local_node_ids[elem.getNodeIndex(i)];
-    }
-}
-
-void MeshPartitioning::writeLocalElementNodeIndicies(std::ostream& os, const Element& elem,
-        const std::vector<unsigned>& local_node_ids)
-{
-    unsigned mat_id = 0; // Materical ID to be set from the mesh data
-    os << mat_id << " " << static_cast<unsigned>(getElementType(elem)) + 1 << " "
-       << elem.getNumberOfNodes() <<" ";
-    for(unsigned i=0; i<elem.getNumberOfNodes(); i++)
-    {
-        os << local_node_ids[elem.getNodeIndex(i)] << " ";
-    }
-    os << "\n";
-}
-
-}   // namespace MeshLib
-
diff --git a/Applications/Utils/ModelPreparation/PartitionMesh/MeshPartitioning.h b/Applications/Utils/ModelPreparation/PartitionMesh/MeshPartitioning.h
deleted file mode 100644
index 9a3a5ed31ddbd8f5d32e9c518b5bb790c50bd98c..0000000000000000000000000000000000000000
--- a/Applications/Utils/ModelPreparation/PartitionMesh/MeshPartitioning.h
+++ /dev/null
@@ -1,86 +0,0 @@
-/*!
-  \file MeshPartitioning.h
-  \date   2016.05
-
-  \brief  Declare a class to perform mesh partitioning
-
-  \copyright
-  Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
-             Distributed under a Modified BSD License.
-               See accompanying file LICENSE.txt or
-               http://www.opengeosys.org/project/license
-
-*/
-
-#ifndef MESH_PARTITIONING_H_
-#define MESH_PARTITIONING_H_
-
-#include <vector>
-#include <string>
-#include <fstream>
-
-#include "MeshLib/Mesh.h"
-
-namespace MeshLib
-{
-
-enum class ElementType : unsigned {LINE2, QUAD4, HEX8, TRI3, TET4, PRISM6, PYRAMID5, INVALID};
-
-/// A subdomain mesh.
-class MeshPartitioning : public Mesh
-{
-        typedef long PetscInt;
-    public:
-        /// Copy constructor
-        MeshPartitioning(const MeshLib::Mesh &mesh) : Mesh(mesh)
-        {}
-
-        /// Write mesh to METIS input file
-        void write2METIS(const std::string& file_name);
-
-        /// Partition by element.
-        void partitionByElementMETIS(const std::string& /*file_name_base*/, const unsigned /*npartitions*/)
-        {
-            /* to be decided */
-        }
-
-        /// Partition by node.
-        void partitionByNodeMETIS(const std::string& file_name_base, const unsigned npartitions,
-                                  const bool output_binary);
-
-    private:
-        ElementType getElementType(const Element& elem);
-
-        /*!
-             \brief get integer variables, which are used to define an element
-             \param elem            Element
-             \param local_node_ids  Local node indicies of a partition
-             \param elem_info       An vector holds all integer variables of element definitions
-             \param counter         Recorder of the number of integer variables.
-        */
-        void getElementIntegerVariables(const Element& elem,
-                                        const std::vector<unsigned>& local_node_ids,
-                                        std::vector<PetscInt>& elem_info,
-                                        PetscInt& counter);
-
-        /*!
-            \brief Write local indicies of element nodes to a ASCII file
-            \param os              Output stream
-            \param elem            Element
-            \param local_node_ids  Local node indicies of a partition
-        */
-        void writeLocalElementNodeIndicies(std::ostream& os, const Element& elem,
-                                           const std::vector<unsigned>& local_node_ids);
-
-        struct NodeStruct
-        {
-            PetscInt id;
-            double x;
-            double y;
-            double z;
-        };
-};
-
-}   // namespace MeshLib
-
-#endif // MESH_PARTITIONING_H_
diff --git a/Applications/Utils/ModelPreparation/PartitionMesh/NodeWiseMeshPartitioner.cpp b/Applications/Utils/ModelPreparation/PartitionMesh/NodeWiseMeshPartitioner.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9e72f7c87505a0e61b582a47fc7900d216307ba8
--- /dev/null
+++ b/Applications/Utils/ModelPreparation/PartitionMesh/NodeWiseMeshPartitioner.cpp
@@ -0,0 +1,537 @@
+/*!
+  \file NodeWiseMeshPartitioner.cpp
+  \date   2016.05
+
+  \brief  Define the members of class NodeWiseMeshPartitioner
+
+  \copyright
+  Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+             Distributed under a Modified BSD License.
+               See accompanying file LICENSE.txt or
+               http://www.opengeosys.org/project/license
+
+*/
+
+#include "NodeWiseMeshPartitioner.h"
+
+#include <limits>
+#include <iomanip>
+#include <cstdio>  // for binary output
+
+#include <logog/include/logog.hpp>
+
+#include "BaseLib/Error.h"
+
+#include "MeshLib/IO/VtkIO/VtuInterface.h"
+
+#include "MeshLib/Node.h"
+#include "MeshLib/Elements/Element.h"
+
+namespace ApplicationUtils
+{
+struct NodeStruct
+{
+    NodeWiseMeshPartitioner::PetscInt id;
+    double x;
+    double y;
+    double z;
+};
+
+void NodeWiseMeshPartitioner::readMetisData(const std::string& file_name_base)
+{
+    const std::string npartitions_str = std::to_string(_npartitions);
+
+    // Read partitioned mesh data from METIS
+    std::string fname_parts = file_name_base + ".mesh.npart." + npartitions_str;
+
+    std::ifstream npart_in(fname_parts.data());
+    if (!npart_in.is_open())
+    {
+        OGS_FATAL(
+            "Error: cannot open file %s. It may not exist! \
+                   \n Run mpmetis beforehand or use option -m",
+            fname_parts.data());
+    }
+
+    const std::size_t nnodes = _mesh->getNumberOfNodes();
+
+    std::size_t counter = 0;
+    while (!npart_in.eof())
+    {
+        npart_in >> _nodes_partition_ids[counter++] >> std::ws;
+        if (counter == nnodes)
+          break;
+    }
+
+    npart_in.close();
+
+    if( counter != nnodes)
+    {
+        OGS_FATAL(
+            "Error: data in %s are less than expected.", fname_parts.data());
+    }
+
+    // TEST  std::remove(fname_parts.c_str());
+}
+
+void NodeWiseMeshPartitioner::partitionByMETIS(const bool is_mixed_hl_elem)
+{
+    std::vector<MeshLib::Node*> const& nodes = _mesh->getNodes();
+    for (std::size_t part_id = 0; part_id < _partitions.size(); part_id++)
+    {
+        auto& partition = _partitions[part_id];
+
+        INFO("Processing partition: %d", part_id);
+
+        // Find non-ghost nodes in this partition
+        // -- Extra nodes for high order elements
+        std::vector<MeshLib::Node*> extra_nodes;
+        for (std::size_t i = 0; i < _mesh->getNumberOfNodes(); i++)
+        {
+            if (_nodes_partition_ids[i] == part_id)
+            {
+                if (is_mixed_hl_elem)
+                {
+                    if (i < _mesh->getNumberOfBaseNodes())
+                        partition.nodes.push_back(nodes[i]);
+                    else
+                        extra_nodes.push_back(nodes[i]);
+                }
+                {
+                    partition.nodes.push_back(nodes[i]);
+                }
+            }
+        }
+        partition.number_of_non_ghost_base_nodes = partition.nodes.size();
+        partition.number_of_non_ghost_nodes =
+            partition.number_of_non_ghost_base_nodes + extra_nodes.size();
+
+        // Find elements that are bellowed to this partition
+        std::vector<MeshLib::Element*> const& elements = _mesh->getElements();
+        for (std::size_t elem_id = 0; elem_id < elements.size(); elem_id++)
+        {
+            const auto* elem = elements[elem_id];
+            if (_elements_status[elem_id])
+                continue;
+
+            std::size_t non_ghost_node_number = 0;
+            for (unsigned i = 0; i < elem->getNumberOfNodes(); i++)
+            {
+                if (_nodes_partition_ids[elem->getNodeIndex(i)] == part_id)
+                {
+                    non_ghost_node_number++;
+                }
+            }
+
+            if (non_ghost_node_number == 0)
+                continue;
+
+            if (non_ghost_node_number == elem->getNumberOfNodes())
+            {
+                partition.regular_elements.push_back(elem);
+                _elements_status[elem_id] = true;
+            }
+            else
+            {
+                partition.ghost_elements.push_back(elem);
+            }
+        }
+
+        // Find the ghost nodes of this partition
+        std::vector<bool> nodes_reserved(_mesh->getNumberOfNodes(), false);
+        for (const auto* ghost_elem : partition.ghost_elements)
+        {
+            for (unsigned i = 0; i < ghost_elem->getNumberOfNodes(); i++)
+            {
+                const unsigned node_id = ghost_elem->getNodeIndex(i);
+                if (nodes_reserved[node_id])
+                    continue;
+
+                if (_nodes_partition_ids[node_id] != part_id)
+                {
+                    if (is_mixed_hl_elem)
+                    {
+                        if (node_id < _mesh->getNumberOfBaseNodes())
+                            partition.nodes.push_back(nodes[node_id]);
+                        else
+                            extra_nodes.push_back(nodes[node_id]);
+                    }
+                    else
+                    {
+                        partition.nodes.push_back(nodes[node_id]);
+                    }
+                    nodes_reserved[node_id] = true;
+                }
+            }
+        }
+        partition.number_of_base_nodes = partition.nodes.size();
+
+        if (is_mixed_hl_elem)
+            partition.nodes.insert(partition.nodes.end(), extra_nodes.begin(),
+                                   extra_nodes.end());
+    }
+
+    renumberNodeIndecies();
+}
+
+void NodeWiseMeshPartitioner::renumberNodeIndecies()
+{
+    std::size_t node_global_id_offset = 0;
+    for (auto& partition : _partitions)
+    {
+        // Renumber the global indecies.
+        // -- Base nodes
+        for (std::size_t i = 0; i < partition.number_of_non_ghost_base_nodes;
+             i++)
+        {
+            _nodes_global_ids[partition.nodes[i]->getID()] =
+                node_global_id_offset;
+            node_global_id_offset++;
+        }
+        // -- Nodes for high order elements.
+        const std::size_t end_id = partition.number_of_base_nodes +
+                                   partition.number_of_non_ghost_nodes -
+                                   partition.number_of_non_ghost_base_nodes;
+        for (std::size_t i = partition.number_of_base_nodes; i < end_id; i++)
+        {
+            _nodes_global_ids[partition.nodes[i]->getID()] =
+                node_global_id_offset;
+            node_global_id_offset++;
+        }
+    }
+}
+
+void NodeWiseMeshPartitioner::writeMETIS(const std::string& file_name)
+{
+    std::ofstream os(file_name.data(), std::ios::trunc);
+    if (!os.is_open())
+    {
+        OGS_FATAL("Error: cannot open file %s. It may not exist! ",
+                  file_name.data());
+    }
+    std::vector<MeshLib::Element*> const& elements = _mesh->getElements();
+    os << elements.size() << " \n";
+    for (const auto* elem : elements)
+    {
+        os << elem->getNodeIndex(0) + 1;
+        for (unsigned j = 1; j < elem->getNumberOfNodes(); j++)
+        {
+            os << " " << elem->getNodeIndex(j) + 1;
+        }
+        os << "\n";
+    }
+    os.flush();
+}
+
+NodeWiseMeshPartitioner::PetscInt
+NodeWiseMeshPartitioner::getNumberOfIntegerVariablesOfElements(
+    const std::vector<const MeshLib::Element*>& elements) const
+{
+    PetscInt nmb_element_idxs = 3 * elements.size();
+    for (const auto* elem : elements)
+    {
+        nmb_element_idxs += elem->getNumberOfNodes();
+    }
+    return nmb_element_idxs;
+}
+
+void NodeWiseMeshPartitioner::writeBinary(const std::string& file_name_base)
+{
+    const std::string npartitions_str = std::to_string(_npartitions);
+
+    // Output configuration data
+    std::string fname =
+        file_name_base + "_partitioned_msh_cfg" + npartitions_str + ".bin";
+    FILE* of_bin_cfg = fopen(fname.c_str(), "wb");
+
+    const PetscInt num_config_data = 14;
+    PetscInt config_data[num_config_data];
+    // node rank offset
+    config_data[10] = 0;
+    // element rank offset
+    config_data[11] = 0;
+    // ghost element rank offset
+    config_data[12] = 0;
+    // Reserved
+    config_data[13] = 0;
+    std::vector<PetscInt> num_elem_integers(_partitions.size());
+    std::vector<PetscInt> num_g_elem_integers(_partitions.size());
+    std::size_t loop_id = 0;
+    for (const auto& partition : _partitions)
+    {
+        config_data[0] = partition.nodes.size();
+        config_data[1] = partition.number_of_base_nodes;
+        config_data[2] = partition.regular_elements.size();
+        config_data[3] = partition.ghost_elements.size();
+        config_data[4] = partition.number_of_non_ghost_base_nodes;
+        config_data[5] = partition.number_of_non_ghost_nodes;
+        config_data[6] = _mesh->getNumberOfBaseNodes();
+        config_data[7] = _mesh->getNumberOfNodes();
+        config_data[8] =
+            getNumberOfIntegerVariablesOfElements(partition.regular_elements);
+        config_data[9] =
+            getNumberOfIntegerVariablesOfElements(partition.ghost_elements);
+
+        fwrite(config_data, 1, num_config_data * sizeof(PetscInt), of_bin_cfg);
+
+        config_data[10] += config_data[0] * sizeof(NodeStruct);
+
+        // Update offsets
+        num_elem_integers[loop_id] =
+            partition.regular_elements.size() + config_data[8];
+        // Offset the ending enrtry of the element integer variales of
+        // the non-ghost elements of this partition in the vector of elem_info.
+        config_data[11] += num_elem_integers[loop_id] * sizeof(PetscInt);
+        // Offset the ending enrtry of the element integer variales of
+        // the ghost elements of this partition in the vector of elem_info.
+        num_g_elem_integers[loop_id] =
+            partition.ghost_elements.size() + config_data[9];
+        config_data[12] += num_g_elem_integers[loop_id] * sizeof(PetscInt);
+
+        loop_id++;
+    }
+    fclose(of_bin_cfg);
+
+    // Output elements
+    fname = file_name_base + "_partitioned_msh_ele" + npartitions_str + ".bin";
+    FILE* of_bin_ele = fopen(fname.c_str(), "wb");
+    fname =
+        file_name_base + "_partitioned_msh_ele_g" + npartitions_str + ".bin";
+    FILE* of_bin_ele_g = fopen(fname.c_str(), "wb");
+    for (std::size_t i = 0; i < _partitions.size(); i++)
+    {
+        const auto& partition = _partitions[i];
+
+        // Set the local node indecies of the current partition.
+        PetscInt node_local_id_offset = 0;
+        std::vector<PetscInt> nodes_local_ids(_mesh->getNumberOfNodes(), -1);
+        for (const auto* node : partition.nodes)
+        {
+            nodes_local_ids[node->getID()] = node_local_id_offset;
+            node_local_id_offset++;
+        }
+
+        // An vector contians all element integer variales of
+        // the non-ghost elements of this partition
+        std::vector<PetscInt> ele_info(num_elem_integers[i]);
+
+        // Non-ghost elements.
+        PetscInt counter = partition.regular_elements.size();
+
+        for (std::size_t j = 0; j < partition.regular_elements.size(); j++)
+        {
+            const auto* elem = partition.regular_elements[j];
+            ele_info[j] = counter;
+            getElementIntegerVariables(*elem, nodes_local_ids, ele_info,
+                                       counter);
+        }
+        // Write vector data of non-ghost elements
+        fwrite(&ele_info[0], 1, (num_elem_integers[i]) * sizeof(PetscInt),
+               of_bin_ele);
+
+        // Ghost elements
+        ele_info.resize(num_g_elem_integers[i]);
+
+        counter = partition.ghost_elements.size();
+
+        for (std::size_t j = 0; j < partition.ghost_elements.size(); j++)
+        {
+            const auto* elem = partition.ghost_elements[j];
+            ele_info[j] = counter;
+            getElementIntegerVariables(*elem, nodes_local_ids, ele_info,
+                                       counter);
+        }
+        // Write vector data of ghost elements
+        fwrite(&ele_info[0], 1, (num_g_elem_integers[i]) * sizeof(PetscInt),
+               of_bin_ele_g);
+    }
+    fclose(of_bin_ele);
+    fclose(of_bin_ele_g);
+
+    // Output an array of all nodes
+    fname = file_name_base + "_partitioned_msh_nod" + npartitions_str + ".bin";
+    FILE* of_bin_nod = fopen(fname.c_str(), "wb");
+    for (const auto& partition : _partitions)
+    {
+        std::vector<NodeStruct> nodes_buffer;
+        nodes_buffer.reserve(partition.nodes.size());
+
+        for (const auto* node : partition.nodes)
+        {
+            double const* coords = node->getCoords();
+            NodeStruct node_struct;
+            node_struct.id = _nodes_global_ids[node->getID()];
+            node_struct.x = coords[0];
+            node_struct.y = coords[1];
+            node_struct.z = coords[2];
+            nodes_buffer.emplace_back(node_struct);
+        }
+        fwrite(&nodes_buffer[0], sizeof(NodeStruct), partition.nodes.size(),
+               of_bin_nod);
+    }
+    fclose(of_bin_nod);
+}
+
+void NodeWiseMeshPartitioner::writeASCII(const std::string& file_name_base)
+{
+    const std::string npartitions_str = std::to_string(_npartitions);
+
+    // Write the configuration data
+    std::string fname =
+        file_name_base + "_partitioned_cfg" + npartitions_str + ".msh";
+    std::fstream os_subd_head(fname.c_str(), std::ios::out | std::ios::trunc);
+    const std::string mesh_info =
+        "Subdomain mesh ("
+        "Number of nodes; Number of base nodes;"
+        " Number of regular elements; Number of ghost elements;"
+        " Number of non-ghost base nodes; Number of non-ghost nodes"
+        " Number of base nodes of the global mesh;"
+        " Number of nodes of the global mesh;"
+        " Number of integer variables to define non-ghost elements;"
+        " Number of integer variables to define ghost elements.)";
+    os_subd_head << mesh_info << "\n";
+    os_subd_head << _npartitions << "\n";
+
+    for (const auto& partition : _partitions)
+    {
+        os_subd_head << partition.nodes.size();
+        os_subd_head << " " << partition.number_of_base_nodes;
+        os_subd_head << " " << partition.regular_elements.size();
+        os_subd_head << " " << partition.ghost_elements.size();
+        os_subd_head << " " << partition.number_of_non_ghost_base_nodes;
+        os_subd_head << " " << partition.number_of_non_ghost_nodes;
+        os_subd_head << " " << _mesh->getNumberOfBaseNodes();
+        os_subd_head << " " << _mesh->getNumberOfNodes();
+        os_subd_head << " " << getNumberOfIntegerVariablesOfElements(
+                                   partition.regular_elements);
+        os_subd_head << " " << getNumberOfIntegerVariablesOfElements(
+                                   partition.ghost_elements)
+                     << " 0 \n";
+    }
+    os_subd_head.close();
+
+    // Write the elements of each partitions
+    fname = file_name_base + "_partitioned_elems_" + npartitions_str + ".msh";
+    std::fstream os_subd(fname.c_str(), std::ios::out | std::ios::trunc);
+    for (const auto& partition : _partitions)
+    {
+        // Set the local node indecies of the current partition.
+        PetscInt node_local_id_offset = 0;
+        std::vector<PetscInt> nodes_local_ids(_mesh->getNumberOfNodes(), -1);
+        for (const auto* node : partition.nodes)
+        {
+            nodes_local_ids[node->getID()] = node_local_id_offset;
+            node_local_id_offset++;
+        }
+
+        for (const auto* elem : partition.regular_elements)
+        {
+            writeLocalElementNodeIndicies(os_subd, *elem, nodes_local_ids);
+        }
+        for (const auto* elem : partition.ghost_elements)
+        {
+            writeLocalElementNodeIndicies(os_subd, *elem, nodes_local_ids);
+        }
+        os_subd << std::endl;
+    }
+    os_subd.close();
+
+    // Write the nodes of each partitions
+    fname = file_name_base + "_partitioned_nodes_" + npartitions_str + ".msh";
+    std::fstream os_subd_node(fname.c_str(), std::ios::out | std::ios::trunc);
+    os_subd_node.precision(std::numeric_limits<double>::digits10);
+    os_subd_node.setf(std::ios::scientific);
+
+    for (const auto& partition : _partitions)
+    {
+        for (const auto* node : partition.nodes)
+        {
+            double const* coords = node->getCoords();
+            os_subd_node << _nodes_global_ids[node->getID()] << " " << coords[0]
+                         << " " << coords[1] << " " << coords[2] << "\n";
+        }
+        os_subd_node << std::endl;
+    }
+    os_subd_node.close();
+}
+
+void NodeWiseMeshPartitioner::resetGlobalNodeIndecis()
+{
+    for (std::size_t i = 0; i < _mesh->_nodes.size(); i++)
+    {
+        _mesh->_nodes[i]->setID(_nodes_global_ids[i]);
+    }
+    // sort
+    std::sort(_mesh->_nodes.begin(), _mesh->_nodes.end(),
+              [](const MeshLib::Node* a, const MeshLib::Node* b) {
+                  return a->getID() < b->getID();
+              });
+}
+
+void NodeWiseMeshPartitioner::writeGlobalMeshVTU(
+    const std::string& file_name_base)
+{
+    resetGlobalNodeIndecis();
+    MeshLib::IO::VtuInterface writer(_mesh.get());
+    const std::string npartitions_str = std::to_string(_npartitions);
+    writer.writeToFile(file_name_base + "_node_id_renumbered_partitions_" +
+                       npartitions_str + ".vtu");
+}
+
+ElementType NodeWiseMeshPartitioner::getElementType(
+    const MeshLib::Element& elem)
+{
+    switch (elem.getCellType())
+    {
+        case MeshLib::CellType::LINE2:
+            return ElementType::LINE2;
+        case MeshLib::CellType::QUAD4:
+            return ElementType::QUAD4;
+        case MeshLib::CellType::HEX8:
+            return ElementType::HEX8;
+        case MeshLib::CellType::TRI3:
+            return ElementType::TRI3;
+        case MeshLib::CellType::PYRAMID5:
+            return ElementType::PYRAMID5;
+        case MeshLib::CellType::PRISM6:
+            return ElementType::PRISM6;
+        default:
+            OGS_FATAL("Invalid element type in element %d", elem.getID());
+    }
+}
+
+void NodeWiseMeshPartitioner::getElementIntegerVariables(
+    const MeshLib::Element& elem,
+    const std::vector<PetscInt>& local_node_ids,
+    std::vector<PetscInt>& elem_info,
+    PetscInt& counter)
+{
+    unsigned mat_id = 0;  // TODO: Materical ID to be set from the mesh data
+    const PetscInt nn = elem.getNumberOfNodes();
+    ;
+    elem_info[counter++] = mat_id;
+    elem_info[counter++] = static_cast<unsigned>(getElementType(elem)) + 1;
+    elem_info[counter++] = nn;
+
+    for (PetscInt i = 0; i < nn; i++)
+    {
+        elem_info[counter++] = local_node_ids[elem.getNodeIndex(i)];
+    }
+}
+
+void NodeWiseMeshPartitioner::writeLocalElementNodeIndicies(
+    std::ostream& os,
+    const MeshLib::Element& elem,
+    const std::vector<PetscInt>& local_node_ids)
+{
+    unsigned mat_id = 0;  // TODO: Materical ID to be set from the mesh data
+    os << mat_id << " " << static_cast<unsigned>(getElementType(elem)) + 1
+       << " " << elem.getNumberOfNodes() << " ";
+    for (unsigned i = 0; i < elem.getNumberOfNodes(); i++)
+    {
+        os << " " << local_node_ids[elem.getNodeIndex(i)];
+    }
+    os << "\n";
+}
+
+}  // namespace MeshLib
diff --git a/Applications/Utils/ModelPreparation/PartitionMesh/NodeWiseMeshPartitioner.h b/Applications/Utils/ModelPreparation/PartitionMesh/NodeWiseMeshPartitioner.h
new file mode 100644
index 0000000000000000000000000000000000000000..1b8bfe853847e677c1e26f3bdc9a0cb5ef556c4e
--- /dev/null
+++ b/Applications/Utils/ModelPreparation/PartitionMesh/NodeWiseMeshPartitioner.h
@@ -0,0 +1,164 @@
+/*!
+  \file NodeWiseMeshPartitioner.h
+  \date   2016.05
+
+  \brief  Declare a class to perform node wise mesh partitioning
+
+  \copyright
+  Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+             Distributed under a Modified BSD License.
+               See accompanying file LICENSE.txt or
+               http://www.opengeosys.org/project/license
+
+*/
+
+#ifndef NODE_WISE_MESH_PARTITIONER_H_
+#define NODE_WISE_MESH_PARTITIONER_H_
+
+#include <memory>
+#include <vector>
+#include <string>
+#include <fstream>
+
+#include "MeshLib/Mesh.h"
+
+namespace ApplicationUtils
+{
+enum class ElementType : unsigned
+{
+    LINE2,
+    QUAD4,
+    HEX8,
+    TRI3,
+    TET4,
+    PRISM6,
+    PYRAMID5,
+    INVALID
+};
+
+///  A subdomain mesh.
+struct Partition
+{
+    std::vector<MeshLib::Node*> nodes;  ///< nodes.
+    std::size_t number_of_non_ghost_base_nodes;
+    std::size_t number_of_non_ghost_nodes;
+    std::size_t number_of_base_nodes;
+    /// Non ghost elements
+    std::vector<const MeshLib::Element*> regular_elements;
+    std::vector<const MeshLib::Element*> ghost_elements;
+};
+
+/// Mesh partitioner.
+class NodeWiseMeshPartitioner
+{
+public:
+    typedef long PetscInt;
+
+public:
+    /*!
+     * \param file_name_base The prefix of the file name.
+     * \param num_partitions Number of partitions,
+     * \param mesh           Pointer to a mesh object.
+     */
+    NodeWiseMeshPartitioner(const PetscInt num_partitions,
+                            std::unique_ptr<MeshLib::Mesh>& mesh)
+        : _npartitions(num_partitions),
+          _partitions(num_partitions),
+          _mesh(std::move(mesh)),
+          _nodes_global_ids(_mesh->getNumberOfNodes()),
+          _nodes_partition_ids(_mesh->getNumberOfNodes()),
+          _elements_status(_mesh->getNumberOfElements(), false)
+    {
+    }
+    ~NodeWiseMeshPartitioner() = default;
+
+    /// Partition by node.
+    /// \param is_mixed_hl_elem Flag to indicate whether the elements of
+    /// a mesh can be used for both linear and high order interpolation
+    void partitionByMETIS(const bool is_mixed_hl_elem);
+
+    void resetGlobalNodeIndecis();
+
+    /// Read metis data
+    /// \param file_name_base The prefix of the file name.
+    void readMetisData(const std::string& file_name_base);
+
+    /// Write mesh to METIS input file
+    /// \param file_name File name with an extension of mesh.
+    void writeMETIS(const std::string& file_name);
+
+    /// Write the partitions into ASCII files
+    /// \param file_name_base The prefix of the file name.
+    void writeASCII(const std::string& file_name_base);
+
+    /// Write the partitions into binary files
+    /// \param file_name_base The prefix of the file name.
+    void writeBinary(const std::string& file_name_base);
+
+    /// Write the global mesh into a VTU file
+    /// \param file_name_base The prefix of the file name.
+    void writeGlobalMeshVTU(const std::string& file_name_base);
+
+private:
+    /// Number of partitions.
+    PetscInt _npartitions;
+
+    /// Data for all  partitions.
+    std::vector<Partition> _partitions;
+
+    /// Pointer to a mesh object.
+    std::unique_ptr<MeshLib::Mesh> _mesh;
+
+    /// Global IDs of all nodes after partitioning.
+    std::vector<std::size_t> _nodes_global_ids;
+
+    /// Partition IDs of each nodes.
+    std::vector<std::size_t> _nodes_partition_ids;
+
+    /// Flags to indicate the status of all elements.
+    std::vector<bool> _elements_status;
+
+    void renumberNodeIndecies();
+
+    /*!
+       Calculate the totoal number of integer variables of an element
+       vector.
+           Each element has three integer variables for material ID,
+       element type, number of nodes of the element. Therefore
+       the total number of the integers in an element vector is
+        3 * vector size + sum (number of nodes of each element)
+     \param is_ghost Flag to indicate ghost elements or not
+   */
+    PetscInt getNumberOfIntegerVariablesOfElements(
+        const std::vector<const MeshLib::Element*>& elements) const;
+
+    ElementType getElementType(const MeshLib::Element& elem);
+
+    /*!
+         \brief get integer variables, which are used to define an element
+         \param elem            Element
+         \param local_node_ids  Local node indicies of a partition
+         \param elem_info       An vector holds all integer variables of element
+       definitions
+         \param counter         Recorder of the number of integer variables.
+    */
+    void getElementIntegerVariables(const MeshLib::Element& elem,
+                                    const std::vector<PetscInt>& local_node_ids,
+                                    std::vector<PetscInt>& elem_info,
+                                    PetscInt& counter);
+
+    /*!
+        \brief Write local indicies of element nodes to a ASCII file
+        \param os              Output stream
+        \param elem            Element
+        \param local_node_ids  Local node indicies of a partition
+    */
+    void writeLocalElementNodeIndicies(
+        std::ostream& os,
+        const MeshLib::Element& elem,
+        const std::vector<PetscInt>& local_node_ids);
+};
+
+}  // namespace MeshLib
+
+#endif  // NODE_WISE_MESH_PARTITIONER_H_
diff --git a/Applications/Utils/ModelPreparation/PartitionMesh/PartitionMesh.cpp b/Applications/Utils/ModelPreparation/PartitionMesh/PartitionMesh.cpp
index 622e51f18d1d089d692614150c4571386d836858..dd6b9f8e2fab61c7baea5754698f7fcfaa35fe8d 100644
--- a/Applications/Utils/ModelPreparation/PartitionMesh/PartitionMesh.cpp
+++ b/Applications/Utils/ModelPreparation/PartitionMesh/PartitionMesh.cpp
@@ -12,49 +12,63 @@
 
 */
 
-#include <sstream>
 #include <tclap/CmdLine.h>
 
+#ifdef _MSC_VER
+#include <windows.h>
+#else
+#include <sys/time.h>
+#include <climits>
+#endif
+
 #include "Applications/ApplicationsLib/LogogSetup.h"
 #include "BaseLib/FileTools.h"
-
-#include "MeshLib/IO/readMeshFromFile.h"
-#include "MeshLib/IO/VtkIO/VtuInterface.h"
-
 #include "BaseLib/CPUTime.h"
 #include "BaseLib/RunTime.h"
 
-#include "MeshPartitioning.h"
+#include "MeshLib/IO/readMeshFromFile.h"
+
+#include "NodeWiseMeshPartitioner.h"
 
-int main (int argc, char* argv[])
+int main(int argc, char* argv[])
 {
     ApplicationsLib::LogogSetup logog_setup;
 
-    std::string m_str = "Partition a mesh for parallel computing."
-                        "The tasks of this tool are in twofold:\n"
-                        "1. Convert mesh file to the input file of the partitioning tool,\n"
-                        "2. Partition a mesh using the partitioning tool,\n"
-                        "\tcreate the mesh data of each partition,\n"
-                        "\trenumber the node indicies of each partition,\n"
-                        "\tand output the results for parallel computing.";
+    std::string m_str =
+        "Partition a mesh for parallel computing."
+        "The tasks of this tool are in twofold:\n"
+        "1. Convert mesh file to the input file of the partitioning tool,\n"
+        "2. Partition a mesh using the partitioning tool,\n"
+        "\tcreate the mesh data of each partition,\n"
+        "\trenumber the node indicies of each partition,\n"
+        "\tand output the results for parallel computing.";
 
     TCLAP::CmdLine cmd(m_str, ' ', "0.1");
-    TCLAP::ValueArg<std::string> mesh_input("i", "mesh-input-file",
-                                            "the name of the file containing the input mesh", true,
-                                            "", "file name of input mesh");
+    TCLAP::ValueArg<std::string> mesh_input(
+        "i", "mesh-input-file",
+        "the name of the file containing the input mesh", true, "",
+        "file name of input mesh");
     cmd.add(mesh_input);
-    TCLAP::ValueArg<int> nparts("n", "np", "the number of partitions", false, 2, "integer");
+    TCLAP::ValueArg<int> nparts("n", "np", "the number of partitions", false, 2,
+                                "integer");
     cmd.add(nparts);
 
-    TCLAP::SwitchArg ogs2metis_flag("1", "ogs2metis",
-                                    "Indicator to convert the ogs mesh file to METIS input file", cmd, false);
+    TCLAP::SwitchArg ogs2metis_flag(
+        "s", "ogs2metis",
+        "Indicator to convert the ogs mesh file to METIS input file", cmd,
+        false);
 
-    TCLAP::SwitchArg exe_metis_flag("m","exe_metis","Call mpmetis inside the programme via system().", false);
+    TCLAP::SwitchArg exe_metis_flag(
+        "m", "exe_metis", "Call mpmetis inside the programme via system().",
+        false);
     cmd.add(exe_metis_flag);
-    TCLAP::SwitchArg asci_flag("a","asci","Enable ASCII output.", false);
-    cmd.add(asci_flag);
-    TCLAP::SwitchArg elem_part_flag("e","element_wise","Enable element wise partitioing.", false);
-    cmd.add(elem_part_flag);
+
+    TCLAP::SwitchArg lh_elems_flag(
+        "q", "lh_elements", "Mixed linear and high order elements.", false);
+    cmd.add(lh_elems_flag);
+
+    TCLAP::SwitchArg ascii_flag("a", "ascii", "Enable ASCII output.", false);
+    cmd.add(ascii_flag);
 
     cmd.parse(argc, argv);
 
@@ -65,56 +79,68 @@ int main (int argc, char* argv[])
 
     const std::string ifile_name = mesh_input.getValue();
     const std::string file_name_base = BaseLib::dropFileExtension(ifile_name);
-    MeshLib::MeshPartitioning* mesh = static_cast<MeshLib::MeshPartitioning*>
-                                      (MeshLib::IO::readMeshFromFile(file_name_base + ".vtu"));
-    INFO("Mesh read: %d nodes, %d elements.", mesh->getNumberOfNodes(), mesh->getNumberOfElements());
+    std::unique_ptr<MeshLib::Mesh> mesh_ptr(
+        MeshLib::IO::readMeshFromFile(file_name_base + ".vtu"));
+    INFO("Mesh read: %d nodes, %d elements.",
+         mesh_ptr->getNumberOfNodes(),
+         mesh_ptr->getNumberOfElements());
+
+    ApplicationUtils::NodeWiseMeshPartitioner mesh_partitioner(
+        nparts.getValue(), mesh_ptr);
 
     if (ogs2metis_flag.getValue())
     {
         INFO("Write the mesh into METIS input file.");
-        mesh->write2METIS(BaseLib::dropFileExtension(ifile_name) + ".mesh");
+        mesh_partitioner.writeMETIS(BaseLib::dropFileExtension(ifile_name) +
+                                    ".mesh");
     }
     else
     {
-        if ( elem_part_flag.getValue() )
-        {
-            INFO("Partitioning the mesh in the element wise way ...");
-            mesh->partitionByElementMETIS( file_name_base, nparts.getValue());
-        }
-        else
-        {
-            const int num_partitions = nparts.getValue();
-            std::string str_nparts = std::to_string(num_partitions);
+        const int num_partitions = nparts.getValue();
 
-            // Execute mpmetis via system(...)
-            if (num_partitions > 1 && exe_metis_flag.getValue())
+        // Execute mpmetis via system(...)
+        if (num_partitions > 1 && exe_metis_flag.getValue())
+        {
+            INFO("METIS is running ...");
+            const std::string exe_name = argv[0];
+            std::string exe_path =
+                exe_name.substr(0, exe_name.find_last_of("/\\"));
+            INFO("Path to mtmetis is: \n\t%s", exe_path.c_str());
+
+            std::string mpmetis_com =
+                exe_path + "/mpmetis " + " -gtype=nodal " + file_name_base +
+                ".mesh " + std::to_string(nparts.getValue());
+
+            int status = system(mpmetis_com.c_str());
+            if (status != 0)
             {
-                INFO("METIS is running ...");
-                std::stringstream ss;
-                ss << "mpmetis " << " -gtype=nodal "
-                   << file_name_base + ".mesh "
-                   << nparts.getValue();
-
-                int status = system(ss.str().c_str());
+                INFO("Failed in system calling.");
                 INFO("Return value of system calling %d ", status);
+                return EXIT_FAILURE;
             }
+        }
 
-            INFO("Partitioning the mesh in the node wise way ...");
-            // Binary output is default.
-            mesh->partitionByNodeMETIS(file_name_base, num_partitions,
-                                       !asci_flag.getValue());
+        mesh_partitioner.readMetisData(file_name_base);
 
-            INFO("Write the mesh with renumbered node indicies into VTU");
-            MeshLib::IO::VtuInterface writer(mesh);
-            writer.writeToFile(file_name_base + "_node_id_renumbered_partitions_"
-                               + str_nparts + ".vtu");
+        INFO("Partitioning the mesh in the node wise way ...");
+        mesh_partitioner.partitionByMETIS(lh_elems_flag.getValue());
+        if (ascii_flag.getValue())
+        {
+            INFO("Write the data of partitions into ASCII files ...");
+            mesh_partitioner.writeASCII(file_name_base);
+        }
+        else
+        {
+            INFO("Write the data of partitions into binary files ...");
+            mesh_partitioner.writeBinary(file_name_base);
         }
-    }
 
-    INFO( "Total runtime: %g s.\n", run_timer.elapsed() );
-    INFO( "Total CPU time: %g s.\n", CPU_timer.elapsed() );
+        INFO("Write the mesh with renumbered node indicies into VTU ...");
+        mesh_partitioner.writeGlobalMeshVTU(file_name_base);
+    }
 
-    delete mesh;
+    INFO("Total runtime: %g s.", run_timer.elapsed());
+    INFO("Total CPU time: %g s.\n", CPU_timer.elapsed());
 
     return EXIT_SUCCESS;
 }
diff --git a/MathLib/Point3dWithID.h b/MathLib/Point3dWithID.h
index fdf7b9360b11012f6a40cf84bbe57bfd36f896d8..2bcf36ec2a3ba000b2963f719711de67f62ca534 100644
--- a/MathLib/Point3dWithID.h
+++ b/MathLib/Point3dWithID.h
@@ -63,8 +63,6 @@ public:
 
     std::size_t getID() const { return _id; }
 
-    void setID(const std::size_t id) { _id = id; }
-
 protected:
     std::size_t _id;
 };
diff --git a/MeshLib/Elements/Element.cpp b/MeshLib/Elements/Element.cpp
index adaa43d8996569a768313b64d0b275bd85648f26..f6d0871ab7550d2660892e2563adc280927bea17 100644
--- a/MeshLib/Elements/Element.cpp
+++ b/MeshLib/Elements/Element.cpp
@@ -140,7 +140,7 @@ unsigned Element::getNodeIDinElement(const MeshLib::Node* node) const
     return std::numeric_limits<unsigned>::max();
 }
 
-Node* Element::getNode(unsigned i) const
+const Node* Element::getNode(unsigned i) const
 {
 #ifndef NDEBUG
     if (i < getNumberOfNodes())
diff --git a/MeshLib/Elements/Element.h b/MeshLib/Elements/Element.h
index 467da896c7639a927927a1a104f91a08bd09be4e..84fe08ae654366cc4e5ee8a885c36f02be558805 100644
--- a/MeshLib/Elements/Element.h
+++ b/MeshLib/Elements/Element.h
@@ -66,7 +66,7 @@ public:
      * modifiable by the user) instance of class Node or a NULL pointer
      * @sa Element::getNodeIndex()
      */
-    Node* getNode(unsigned i) const;
+    const Node* getNode(unsigned i) const;
 
     /**
      * (Re)Sets the node of the element.
diff --git a/MeshLib/Mesh.h b/MeshLib/Mesh.h
index d138bf3553a2c6a33e0b93c40db07bab7b5545d3..bb43ae634ae6816f318ebf81765afd9b0da7ea40 100644
--- a/MeshLib/Mesh.h
+++ b/MeshLib/Mesh.h
@@ -24,6 +24,11 @@
 #include "MeshEnums.h"
 #include "Properties.h"
 
+namespace ApplicationUtils
+{
+    class NodeWiseMeshPartitioner;
+}
+
 namespace MeshLib
 {
     class Node;
@@ -37,6 +42,8 @@ class Mesh : BaseLib::Counter<Mesh>
     /* friend functions: */
     friend void removeMeshNodes(MeshLib::Mesh &mesh, const std::vector<std::size_t> &nodes);
 
+    friend class ApplicationUtils::NodeWiseMeshPartitioner;
+
 public:
     /// Constructor using a mesh name and an array of nodes and elements
     /// @param name          Mesh name.
diff --git a/MeshLib/Node.h b/MeshLib/Node.h
index 6bccb298703ffb58c1d532f25809c315af89a8f0..de493e4f2d462d3d1c92ac334638dc90d6e0b48f 100644
--- a/MeshLib/Node.h
+++ b/MeshLib/Node.h
@@ -68,7 +68,7 @@ public:
     }
 
     /// Sets the ID of a node to the given value.
-    void setID(std::size_t id) { this->_id = id; }
+    void setID(std::size_t id) { _id = id; }
 
 protected:
     /// Update coordinates of a node.