diff --git a/Applications/Utils/PartitionMesh/MeshPartitioning.cpp b/Applications/Utils/PartitionMesh/MeshPartitioning.cpp
index ead10df8c646d945b769d8e58368f512bbb90ba3..199aedee75b75c8fe45245743708db42d2488623 100644
--- a/Applications/Utils/PartitionMesh/MeshPartitioning.cpp
+++ b/Applications/Utils/PartitionMesh/MeshPartitioning.cpp
@@ -14,8 +14,12 @@
 
 #include "MeshPartitioning.h"
 
-#include <fstream>
+#include <iomanip>
+#include <stdio.h> // for binary output
 
+#include "logog/include/logog.hpp"
+
+#include "MeshLib/Node.h"
 #include "MeshLib/Elements/Element.h"
 
 namespace MeshLib
@@ -34,17 +38,464 @@ void MeshPartitioning :: write2METIS(const std::string& file_name)
     }
 }
 
-void MeshPartitioning :: partitionByNode()
+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();
+    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 = fname + "_partitioned_msh_ele" + npartitions_str + ".bin";
+        of_bin_ele = fopen(fname.c_str(), "wb");
+
+        fname = fname + "_partitioned_msh_ele_g" + npartitions_str + ".bin";
+        of_bin_ele_g = fopen(fname.c_str(), "wb");
+
+        fname = fname + "_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);
+    }
+
+    findConnectedElements();
+
+    MyInt global_node_id_offset = 0;
+    std::vector<bool> nodes_part_reserved;
+    nodes_part_reserved.reserve(nnodes);
+    std::vector<unsigned> nodes_local_ids_partition;
+    nodes_local_ids_partition.reserve(nnodes);
+    for (std::size_t i=0; i<nnodes; i++)
+    {
+        nodes_part_reserved.push_back(false);
+        nodes_local_ids_partition.push_back(0);
+    }
+    std::vector<bool> nodes_reseved(nnodes);
+    std::vector<bool> elems_reseved(_elements.size());
+    // For ghost element
+    std::vector< std::vector<unsigned> > elems_non_ghost_nodes_local_ids(_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);
+
+    /// Assume that the maximum node number element amoung all element types are 100.
+    std::vector<unsigned> nonghost_nodes_local_ids(100);
+    std::vector<unsigned> ghost_nodes_local_ids(100);
 
+    for (std::size_t ipart=0; ipart<npartitions; ipart++)
+    {
+        INFO("Processing partition: %d", ipart);
+
+        partition_start_node_id[ipart] = partition_nodes.size();
+
+        for (auto node_reserved_flag : nodes_reseved)
+        {
+            node_reserved_flag = 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();
+
+        for (auto elem_reserved_flag : elems_reseved)
+        {
+            elem_reserved_flag = false;
+        }
+
+        // Find elements in this partition
+        std::vector<const Element*> partition_regular_elements;
+        std::vector<const Element*> partition_ghost_elements;
+
+        const MyInt num_non_ghost_nodes =  end_non_ghost_node_id[ipart]
+                                           - partition_start_node_id[ipart];
+        for (MyInt 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
+            for (const auto elem : _node_connected_elements[node_id])
+            {
+                // If checked
+                if (elems_reseved[elem->getID()])
+                    continue;
+
+                unsigned non_ghost_counter = 0;
+                unsigned ghost_counter = 0;
+                for (unsigned kk=0; kk<elem->getNNodes(); kk++)
+                {
+                    if (nodes_reseved[elem->getNodeIndex(kk)])
+                    {
+                        nonghost_nodes_local_ids[non_ghost_counter] = kk;
+                        non_ghost_counter++;
+                    }
+                    else
+                    {
+                        ghost_nodes_local_ids[ghost_counter] = kk;
+                        ghost_counter++;
+                    }
+                }
+
+                // All nodes of this element are inside this partition
+                if (ghost_nodes_local_ids.size() == 0)
+                {
+                    partition_regular_elements.push_back(elem);
+                }
+                else if (ghost_nodes_local_ids.size() != elem->getNNodes()) // ghost element
+                {
+                    partition_ghost_elements.push_back(elem);
+
+                    elems_non_ghost_nodes_local_ids[elem->getID()].resize(non_ghost_counter);
+                    std::vector<unsigned>& local_node_ids = elems_non_ghost_nodes_local_ids[elem->getID()];
+                    for (unsigned kk=0; kk<non_ghost_counter; kk++)
+                        local_node_ids[kk] = nonghost_nodes_local_ids[kk];
+                }
+
+                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->getNNodes(); 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[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->getNNodes(); k++)
+            {
+                if (nodes_reseved[ghost_elem->getNodeIndex(k)])
+                    continue;
+                nodes_reseved[ghost_elem->getNodeIndex(k)] = true;
+                partition_nodes.push_back(ghost_elem->getNode(k));
+            }
+        }
+
+        partition_end_node_id[ipart] = partition_nodes.size();
+        // Renumber
+        MyInt local_id = 0;
+        for (std::size_t i = partition_start_node_id[ipart];
+                         i < end_non_ghost_node_id[ipart]; i++)
+        {
+            MeshLib::Node* a_node = partition_nodes[i];
+            a_node->setID(global_node_id_offset);
+            nodes_local_ids_partition[i] = 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[i] = local_id;
+            local_id++;
+        }
+
+        // Count the number of integer variables used to define non-ghost elements
+        const MyInt num_non_ghost_elems = partition_regular_elements.size();
+        MyInt nmb_element_idxs = 3 * num_non_ghost_elems;
+        for (MyInt j=0; j<num_non_ghost_elems; j++)
+        {
+            nmb_element_idxs += partition_regular_elements[j]->getNNodes();
+        }
+        const MyInt num_ghost_elems = partition_ghost_elements.size();
+        MyInt nmb_element_idxs_g = 3 * num_ghost_elems;
+        for (MyInt j=0; j<num_ghost_elems; j++)
+        {
+            nmb_element_idxs_g += partition_ghost_elements[j]->getNNodes();
+        }
+
+        std::string ipart_str = std::to_string(ipart);
+
+        // Number of active elements
+        const MyInt offset_e = num_non_ghost_elems + nmb_element_idxs;
+        const MyInt offset_e_g = num_ghost_elems + nmb_element_idxs_g;
+
+        // Reserved for nodes for high order interpolation
+        const MyInt num_non_ghost_nodes_extra = 0;
+
+        const MyInt num_nodes_linear_element =   partition_end_node_id[ipart]
+                                               - partition_start_node_id[ipart];
+        const MyInt num_nodes_quadratic_element = 0;
+
+        const MyInt num_nodes_global_mesh = nnodes;
+        // Reserved for nodes for high order interpolation
+        const MyInt num_nodes_global_mesh_extra = nnodes;
+
+        // Write configuration data and elements of this partition
+        if (output_binary)
+        {
+            const MyInt num_headers = 14;
+            MyInt head[num_headers];
+            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;
+            head[10] = 0;
+            head[11] = 0;
+            head[12] = 0;
+            head[13] = 0;
+
+            fwrite(head, 1, num_headers * sizeof(MyInt), of_bin_cfg);
+
+            head[10] += head[0] * sizeof(NodeStruct);
+            head[11] += offset_e * sizeof(MyInt);
+            head[12] += offset_e_g * sizeof(MyInt);
+
+            // Write non-ghost elements
+            std::vector<MyInt> ele_info(offset_e);
+
+            MyInt counter = num_non_ghost_elems;
+
+            for (MyInt 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(MyInt), of_bin_ele);
+            ele_info.clear();
+
+            // Write ghost elements
+            ele_info.resize(offset_e_g);
+            counter = num_ghost_elems;
+            for (MyInt 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(MyInt), of_bin_ele_g);
+        }
+        else
+        {
+            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_local_ids_partition.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 = partition_end_node_id[ipart] - partition_start_node_id[ipart];
+            std::vector<NodeStruct> nodes_buffer;
+            nodes_buffer.reserve(nnodes);
+
+            for (std::size_t i=partition_end_node_id[ipart]; i<partition_start_node_id[ipart]; i++)
+            {
+                double const* coords = partition_nodes[i]->getCoords();
+                NodeStruct node_struct;
+                node_struct.id = 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, 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_end_node_id[ipart]; i<partition_start_node_id[ipart]; i++)
+            {
+                double const* coords = partition_nodes[i]->getCoords();
+                os_subd_node << partition_nodes[i]->getID() << " "
+                             << coords[0] << " " << coords[1] << " " << coords[2]<<"\n";
+            }
+            os_subd_node << std::endl;
+        }
+        os_subd_node.close();
+    }
+}
+
+void MeshPartitioning::findConnectedElements()
+{
+    _node_connected_elements.resize(_nodes.size());
+
+    for (auto elem : _elements)
+    {
+        for (unsigned i=0; i<elem->getNNodes(); i++)
+        {
+            bool done = false;
+            unsigned node_id = elem->getNodeIndex(i);
+            for (std::size_t j=0; j<_node_connected_elements[node_id].size(); j++)
+            {
+                if ( elem->getID() == _node_connected_elements[node_id][j]->getID())
+                {
+                    done = true;
+                    break;
+                }
+            }
+            if(!done)
+                _node_connected_elements[node_id].push_back(elem);
+        }
+    }
+}
+
+ElementType MeshPartitioning::getElementType(const Element& elem)
+{
+    switch ( elem.getCellType() )
+    {
+        case MeshLib::CellType::LINE2:
+            return LINE2;
+        case MeshLib::CellType::QUAD4:
+            return QUAD4;
+        case MeshLib::CellType::HEX8:
+            return HEX8;
+        case MeshLib::CellType::TRI3:
+            return TRI3;
+        case MeshLib::CellType::PYRAMID5:
+            return PYRAMID5;
+        case MeshLib::CellType::PRISM6:
+            return PRISM6;
+        default:
+            ERR("Invalid element type in element %d", elem.getID());
+            abort();
+    }
 }
 
-void MeshPartitioning :: writeNodePartitionedMeshASCII()
+void MeshPartitioning::getElementIntegerVariables(const Element& elem,
+                                                   const std::vector<unsigned>& local_node_ids,
+                                                   std::vector<MyInt>& elem_info,
+                                                   MyInt& counter)
 {
+    unsigned mat_id = 0; // Materical ID to be set from the mesh data
+    const MyInt nn = elem.getNNodes();;
+    elem_info[counter++] = mat_id;
+    elem_info[counter++] = getElementType(elem) + 1;
+    elem_info[counter++] = nn;
+
+    for(MyInt i=0; i<nn; i++)
+    {
+        elem_info[counter++] = local_node_ids[elem.getNodeIndex(i)];
+    }
 }
 
-void MeshPartitioning :: writeNodePartitionedMeshBinary()
+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 << " " << getElementType(elem) + 1 << " "
+    << elem.getNNodes() <<" ";
+    for(unsigned i=0; i<elem.getNNodes(); i++)
+    {
+        os << local_node_ids[elem.getNodeIndex(i)] << " ";
+    }
+    os << "\n";
 }
 
 }   // namespace MeshLib
diff --git a/Applications/Utils/PartitionMesh/MeshPartitioning.h b/Applications/Utils/PartitionMesh/MeshPartitioning.h
index ae1270479bed591e6592d70f7f1c125d4334fd62..922d31b7e2363ccb4a23472ff6b70765416395d9 100644
--- a/Applications/Utils/PartitionMesh/MeshPartitioning.h
+++ b/Applications/Utils/PartitionMesh/MeshPartitioning.h
@@ -17,15 +17,19 @@
 
 #include <vector>
 #include <string>
+#include <fstream>
 
 #include "MeshLib/Mesh.h"
 
 namespace MeshLib
 {
 
+enum ElementType {LINE2, QUAD4, HEX8, TRI3, TET4, PRISM6, PYRAMID5, INVALID};
+
 /// A subdomain mesh.
 class MeshPartitioning : public Mesh
 {
+        typedef int MyInt; // for PetscInt
     public:
         /// Copy constructor
         MeshPartitioning(const MeshLib::Mesh &mesh) : Mesh(mesh)
@@ -35,22 +39,52 @@ class MeshPartitioning : public Mesh
         void write2METIS(const std::string& file_name);
 
         /// Partition by element.
-        void partitionByElement()
+        void partitionByElementMETIS(const std::string& /*file_name_base*/, const unsigned /*npartitions*/)
         {
             /* to be decided */
         }
 
         /// Partition by node.
-        void partitionByNode();
-
-        /// Write the partitioned mesh into ASCII files.
-        void writeNodePartitionedMeshASCII();
-
-        /// Write the partitioned mesh into binary files.
-        void writeNodePartitionedMeshBinary();
+        void partitionByNodeMETIS(const std::string& file_name_base, const unsigned npartitions,
+                                  const bool output_binary);
 
     private:
-
+        /// Elements connected to each nodes.
+        std::vector< std::vector<Element*> > _node_connected_elements;
+
+        /// Find connected elements for each nodes.
+        void findConnectedElements();
+
+        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<MyInt>& elem_info,
+                                    MyInt& 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
+        {
+            MyInt id;
+            double x;
+            double y;
+            double z;
+        };
 };
 
 }   // namespace MeshLib
diff --git a/Applications/Utils/PartitionMesh/PartitionMesh.cpp b/Applications/Utils/PartitionMesh/PartitionMesh.cpp
index b277abddb8e869ecce300165eb70ad67e7a54933..7ba122e4bf2a94d229413abe971971188ceecc58 100644
--- a/Applications/Utils/PartitionMesh/PartitionMesh.cpp
+++ b/Applications/Utils/PartitionMesh/PartitionMesh.cpp
@@ -19,8 +19,11 @@
 #include "LogogSimpleFormatter.h"
 #include "BaseLib/FileTools.h"
 
+#include "FileIO/readMeshFromFile.h"
 #include "FileIO/VtkIO/VtuInterface.h"
-#include "Legacy/MeshIO.h"
+
+#include "BaseLib/CPUTime.h"
+#include "BaseLib/RunTime.h"
 
 #include "MeshPartitioning.h"
 
@@ -40,30 +43,64 @@ int main (int argc, char* argv[])
                         "\tand output the results for parallel computing.";
 
     TCLAP::CmdLine cmd(m_str, ' ', "0.1");
-    TCLAP::ValueArg<std::string> mesh_in("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");
+    cmd.add(nparts);
 
     TCLAP::SwitchArg ogs2metis_flag("1", "ogs2metis",
                                     "Indicator to convert the ogs mesh file to METIS input file", cmd, false);
 
+    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 output_part_flag("o","output_parititon","Enable output partitions.", false);
+    cmd.add(output_part_flag);
+
     cmd.parse(argc, argv);
 
-    const std::string ifile_name = mesh_in.getValue();
+    BaseLib::RunTime run_timer;
+    run_timer.start();
+    BaseLib::CPUTime CPU_timer;
+    CPU_timer.start();
 
-    MeshLib::MeshPartitioning* mesh = dynamic_cast<MeshLib::MeshPartitioning*>
-                                      (FileIO::VtuInterface::readVTUFile(ifile_name));
+    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*>
+                                      (FileIO::readMeshFromFile(file_name_base + ".vtu"));
     INFO("Mesh read: %d nodes, %d elements.", mesh->getNNodes(), mesh->getNElements());
 
-    if ( ogs2metis_flag.getValue() )
+    if (ogs2metis_flag.getValue())
     {
+        INFO("Write the mesh into METIS input file.");
         mesh->write2METIS(BaseLib::dropFileExtension(ifile_name) + ".mesh");
     }
     else
     {
-
+        if ( elem_part_flag.getValue() )
+        {
+            INFO("Partition the mesh in the element wise way.");
+            mesh->partitionByElementMETIS( file_name_base, nparts.getValue());
+        }
+        else
+        {
+            const int num_partitions = nparts.getValue();
+
+            INFO("Partition the mesh in the node wise way.");
+            mesh->partitionByNodeMETIS(file_name_base, num_partitions, asci_flag.getValue());
+
+            INFO("Write the mesh with renumbered node indicies into VTU");
+            FileIO::VtuInterface writer(mesh);
+            writer.writeToFile(file_name_base + "_" + std::to_string(num_partitions) + ".vtu");
+        }
     }
 
+    INFO( "Total runtime: %g s.\n", run_timer.elapsed() );
+    INFO( "Total CPU time: %g s.\n", CPU_timer.elapsed() );
+
     delete custom_format;
     delete logog_cout;
     LOGOG_SHUTDOWN();
diff --git a/MathLib/Point3dWithID.h b/MathLib/Point3dWithID.h
index 2bcf36ec2a3ba000b2963f719711de67f62ca534..fdf7b9360b11012f6a40cf84bbe57bfd36f896d8 100644
--- a/MathLib/Point3dWithID.h
+++ b/MathLib/Point3dWithID.h
@@ -63,6 +63,8 @@ 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 f6d0871ab7550d2660892e2563adc280927bea17..adaa43d8996569a768313b64d0b275bd85648f26 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();
 }
 
-const Node* Element::getNode(unsigned i) 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 84fe08ae654366cc4e5ee8a885c36f02be558805..467da896c7639a927927a1a104f91a08bd09be4e 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()
      */
-    const Node* getNode(unsigned i) const;
+    Node* getNode(unsigned i) const;
 
     /**
      * (Re)Sets the node of the element.
diff --git a/MeshLib/Node.h b/MeshLib/Node.h
index ee7f8db349702cd25ef3cd95934d4bcbb0651da3..6bccb298703ffb58c1d532f25809c315af89a8f0 100644
--- a/MeshLib/Node.h
+++ b/MeshLib/Node.h
@@ -67,6 +67,9 @@ public:
         return Node(_x[0]-v[0], _x[1]-v[1], _x[2]-v[2]);
     }
 
+    /// Sets the ID of a node to the given value.
+    void setID(std::size_t id) { this->_id = id; }
+
 protected:
     /// Update coordinates of a node.
     /// This method automatically also updates the areas/volumes of all connected elements.
@@ -89,9 +92,6 @@ protected:
         _connected_nodes = connected_nodes;
     }
 
-    /// Sets the ID of a node to the given value.
-    void setID(std::size_t id) { this->_id = id; }
-
     std::vector<Node*> _connected_nodes;
     std::vector<Element*> _elements;
 }; /* class */