diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index 8ab9b1ec872a1b7413d33610144b33ae6b499bd0..212d0a6764c02c7c700ab077fd144da831572ebf 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -32,6 +32,7 @@
 #include "Options.h"
 #include "polynomialBasis.h"
 #include "Gauss.h"
+#include "meshPartitionOptions.h"
 
 #if defined(HAVE_OPENGL)
 #include "drawContext.h"
@@ -418,6 +419,7 @@ binding::binding()
   Msg::registerBindings(this);
   polynomialBasis::registerBindings(this);
   gaussIntegration::registerBindings(this);
+  meshPartitionOptions::registerBindings(this);
 #if defined(HAVE_SOLVER)
   function::registerBindings(this);
   linearSystem<double>::registerBindings(this);
diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index e3c98321b90c179876b7e1d082713c8eb606d881..7d2e4a57b258df6c49c38a578d8d4d078ddcc2f1 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -23,7 +23,8 @@ set(SRC
     BDS.cpp 
     HighOrder.cpp 
       highOrderSmoother.cpp 
-    meshPartition.cpp 
+    meshPartition.cpp
+    meshPartitionOptions.cpp
     meshRefine.cpp
     multiscalePartition.cpp
 )
diff --git a/Mesh/meshPartition.cpp b/Mesh/meshPartition.cpp
index 8e88e8cf462103c4f5c7c7b014b77afb892a66d5..7b1283f9132c638e930b751563c861c4ddc003ee 100644
--- a/Mesh/meshPartition.cpp
+++ b/Mesh/meshPartition.cpp
@@ -55,6 +55,12 @@ extern "C" void METIS_PartGraphKway
 (int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *,
  int *, idxtype *);
 extern "C" void METIS_NodeND(int *n,int *xadj,int *adjncy,int *numflag,int *options,int *perm,int *iperm);
+extern "C" void METIS_mCPartGraphRecursive
+(int *, int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *,
+ int *, idxtype *);
+extern "C" void METIS_mCPartGraphKway
+(int *, int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *,
+ int *, idxtype *);
 
 
 /*==============================================================================
@@ -429,6 +435,7 @@ int PartitionGraph(Graph &graph, meshPartitionOptions &options)
         int n = graph.getNumVertex();
         int wgtflag = 0;
         int numflag = 0;
+        // if metisOptions[0]=0 then default options
         int metisOptions[5];
         int edgeCut;
         const int iSec = 0;
@@ -460,7 +467,7 @@ int PartitionGraph(Graph &graph, meshPartitionOptions &options)
           }
           break;
         case 3:  // Nodal weight
-          //printf("METIS with weights\n");
+          printf("METIS with weights\n");
           metisOptions[0] = 1;
           metisOptions[1] = options.edge_matching;
           metisOptions[2] = 1;
@@ -477,7 +484,39 @@ int PartitionGraph(Graph &graph, meshPartitionOptions &options)
                &graph.partition[graph.section[iSec]]);
           }
           break;
+        case 4:  // Vertices Multi-Contrained Recursive
+          Msg::Info("Vertices Multi-Constrained Recursive Algorithm Used");
+          wgtflag = 2;
+          metisOptions[0] = 1;
+          metisOptions[1] = options.edge_matching;
+          metisOptions[2] = 1;
+          metisOptions[3] = 1;
+          metisOptions[4] = 0;
+          graph.fillWithMultipleWeights(options.ncon,options.getWeightMap());
+          METIS_mCPartGraphRecursive
+            (&n,&options.ncon,&graph.xadj[graph.section[iSec]],
+             &graph.adjncy[graph.section[iSec]], &graph.vwgts[graph.section[iSec]], NULL, &wgtflag, &numflag,
+             &options.num_partitions, metisOptions, &edgeCut,
+             &graph.partition[graph.section[iSec]]);
+          break;
+        case 5:  // Vertices Multi-Constrained K-way
+          Msg::Info("Vertices Multi-Constrained K-way Algorithm Used");
+          wgtflag = 2;
+          metisOptions[0] = 1;
+          metisOptions[1] = options.edge_matching;
+          metisOptions[2] = 1;
+          metisOptions[3] = options.refine_algorithm;
+          metisOptions[4] = 0;
+          if (options.num_partitions > 1) {
+            METIS_mCPartGraphKway
+              (&n,&options.ncon,&graph.xadj[graph.section[iSec]],
+               &graph.adjncy[graph.section[iSec]], &graph.vwgts[graph.section[iSec]], NULL, &wgtflag, &numflag,
+               &options.num_partitions, metisOptions, &edgeCut,
+               &graph.partition[graph.section[iSec]]);
+          }
+          break;
         }
+        Msg::Info("Number of Edges Cut : %d", edgeCut);
       }
       catch(...) {
         Msg::Error("METIS threw an exception");
@@ -574,9 +613,7 @@ int MakeGraph(GModel *const model, Graph &graph, meshPartitionOptions &options,
 /*--------------------------------------------------------------------*
  * Make a graph for the entire domain
  *--------------------------------------------------------------------*/
-
     
-
 //--Get the dimension of the mesh and count the numbers of elements
 
       unsigned numElem[5];
@@ -595,8 +632,10 @@ int MakeGraph(GModel *const model, Graph &graph, meshPartitionOptions &options,
             // Allocate memory for the graph
             const int numGrVert = numElem[ElemTypeTri] + numElem[ElemTypeQuad]
                                   + numElem[ElemTypePolyg];
+            // maximum possible number of corresponding edges for the mesh
             const int maxGrEdge = (numElem[ElemTypeTri]*3 + numElem[ElemTypeQuad]*4
                                    + numElem[ElemTypePolyg]*4)/2;
+  
             graph.allocate(numGrVert, maxGrEdge);
             // Make the graph
             MakeGraphDIM<2>(model->firstFace(), model->lastFace(),
@@ -839,9 +878,9 @@ void MakeGraphDIM(const EntIter begin, const EntIter end,
 
   typedef typename DimTr<DIM>::FaceT FaceT;
   typedef typename LFaceTr<FaceT>::FaceMap FaceMap;
-
+  
   graph.markSection();
-
+  
   FaceMap faceMap;
   GrVertexMap grVertMap;
 
@@ -1425,6 +1464,7 @@ void createPartitionFaces(GModel *model, GFaceCompound *gf, int N,
 //--Explicit instantiations of the routine for adding elements from a
 //--container of entities
 
+// fiter= iterator on faces,  eiter= iterator on edges
 template void MakeGraphDIM<2, GModel::fiter, GModel::eiter>
 (const GModel::fiter begin, const GModel::fiter end,
  const GModel::eiter beginBE, const GModel::eiter endBE,
diff --git a/Mesh/meshPartition.h b/Mesh/meshPartition.h
index 928b5ea020ca19c278fe1e04c95dc3d0e91238c6..928fa8d939f6d20e870590d95162621d2e4d337f 100644
--- a/Mesh/meshPartition.h
+++ b/Mesh/meshPartition.h
@@ -14,7 +14,7 @@
 class GModel;
 class GFace;
 class Graph;
-struct meshPartitionOptions;
+class meshPartitionOptions;
 
 /*******************************************************************************
  *
diff --git a/Mesh/meshPartitionObjects.h b/Mesh/meshPartitionObjects.h
index 134cb58658d730b0fc29359b80560920963e06fe..97830e1f10bffc972a131653e24c20be4847d78e 100644
--- a/Mesh/meshPartitionObjects.h
+++ b/Mesh/meshPartitionObjects.h
@@ -11,6 +11,7 @@
 #include "MElement.h"
 #include "GmshMessage.h"
 #include "Context.h"
+#include "fullMatrix.h"
 
 
 /*******************************************************************************
@@ -96,6 +97,8 @@ class Graph
                                         // partitioner
   std::vector<MElement*> element;       // The element corresponding to each
                                         // graph vertex in 'xadj'
+  fullMatrix<int> *loads;                // Matrix of loads on each partition
+                                        
  private:
   unsigned cIndex;                      // An index for created graph vertices
                                         // (used externally)
@@ -148,6 +151,29 @@ class Graph
     // Translated vertex numbers start from 1
     c2w[grVertMapIt->second.index] = i + 1;
   }
+  void computeLoads(int numParts, int numCon){
+    loads=new fullMatrix<int> (numParts,numCon);
+    for(int i=0; i<numParts;i++){
+      for(int j=0; j<partition.size(); j++){
+        if(partition[j]==i){
+          for(int k=0; k<numCon;k++){
+            (*loads)(i, k)+=vwgts[j*numCon+k];
+          }
+        }
+      }
+    }
+  }
+  void checkLoads(int numParts,  int numCon){
+    printf("******************************************************* \n");
+    for(int i=0; i<numParts;i++){
+      printf("------- PARTITION %d: [", i+1);
+      for(int j=0; j<numCon; j++){
+        printf(" %d", (*loads)(i, j));
+      }
+      printf("] -------\n");
+    }
+    printf("******************************************************* \n");
+  }
   void fillWeights(std::vector<int> wgts)
   {
     int num = 0;
@@ -156,6 +182,19 @@ class Graph
        num++;
     }
   }
+  // Add multiple weights on vertices of the graph given in a map between original element Numbers and their corresponding vector of weights 
+  void fillWithMultipleWeights(int ncon, std::map<int, std::vector<int> > weightMap)
+  {
+    std::vector<MElement*>::iterator eIt;
+    vwgts.resize(element.size()*ncon);
+    int localElNum=0;
+    for(eIt=element.begin();eIt !=element.end();eIt++){
+      for(int i=0; i<ncon; i++){
+        vwgts[localElNum*ncon+i]=weightMap[(*eIt)->getNum()][i];
+      }
+      localElNum+=1;
+    }
+  }
 
   // Add weights per element, as defined in options 
   void fillDefaultWeights() 
diff --git a/Mesh/meshPartitionOptions.cpp b/Mesh/meshPartitionOptions.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..101a7ea3e6ff42ad81928c373f5c18c038c50185
--- /dev/null
+++ b/Mesh/meshPartitionOptions.cpp
@@ -0,0 +1,82 @@
+#include "meshPartitionOptions.h"
+#include "meshPartition.h"
+#include "GModel.h"
+#include <map>
+
+
+meshPartitionOptions::meshPartitionOptions()
+{
+  setDefaults();
+}
+
+  //--Default values
+
+void meshPartitionOptions::setDefaults()
+{
+  partitioner = 2;
+  num_partitions=1;
+  ncon = 0;
+  renumber = 0;
+  global_method = 1;
+  architecture = 1;
+  ndims_tot = 2;
+  mesh_dims[0] = 4;
+  mesh_dims[1] = 1;
+  mesh_dims[2] = 1;
+  goal = 0;
+  local_method = 1;
+  rqi_flag = 1;
+  vmax = 250;
+  ndims = 1;
+  eigtol = 1.E-3;
+  seed = 7654321L;
+  refine_partition = false;
+  internal_vertices = false;
+  refine_map = true;
+  terminal_propogation = false;
+  algorithm = 1;
+  edge_matching = 3;
+  refine_algorithm = 3;
+  createPartitionBoundaries = true;
+  createGhostCells = true;
+  partitionByExtrusion =false;
+  triWeight = 1;
+  quaWeight = 1;
+  tetWeight = 1;
+  priWeight = 1;
+  pyrWeight = 1;
+  hexWeight = 1; 
+}
+
+void meshPartitionOptions::partition(GModel *model)
+{
+  PartitionMesh(model, *this);
+}
+   
+void meshPartitionOptions::setAlgorithm(int algo)
+{
+  algorithm=algo;
+}
+
+
+#include "Bindings.h"
+void meshPartitionOptions::registerBindings(binding *b){
+  classBinding *cb = b->addClass<meshPartitionOptions>("meshPartitionOptions");
+  cb->setDescription("Defines the options for mesh partitioning (CHACO/METIS)");
+  methodBinding *cm;
+  cm = cb->setConstructor<meshPartitionOptions>();
+  cm->setDescription("Build an object that contains all the options needed to perform a partitioning, when creating object default options are applied");
+  cm->setArgNames(NULL);
+  cm = cb->addMethod("partition", &meshPartitionOptions::partition);
+  cm->setDescription("Partition with the specified option the GModel given as argument");
+  cm->setArgNames("GModel", NULL);
+  cm = cb->addMethod("setNumOfPartitions", &meshPartitionOptions::setNumOfPartitions);
+  cm->setDescription("Define the number of partitions desired");
+  cm->setArgNames("numPart", NULL);
+  cm = cb->addMethod("getNumConstraints", &meshPartitionOptions::getNumConstraints);
+  cm->setDescription("Get the number of different weights on the elements");
+  cm->setArgNames(NULL);
+  cm = cb->addMethod("setAlgorithm", &meshPartitionOptions::setAlgorithm);
+  cm->setDescription("Set the partitionning algorithm");
+  cm->setArgNames("algo", NULL);
+}
diff --git a/Mesh/meshPartitionOptions.h b/Mesh/meshPartitionOptions.h
index 56aa905f546e2241cf14b7efb5a326e7b4ea5a05..b90bdd9ecc2221b8093b6c06ac9ab8aba56ef0ae 100644
--- a/Mesh/meshPartitionOptions.h
+++ b/Mesh/meshPartitionOptions.h
@@ -2,134 +2,119 @@
 //
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to <gmsh@geuz.org>.
-
+#include <map>
+#include "GModel.h"
+#include "meshPartition.h"
+#include "fullMatrix.h" 
 #ifndef _MESH_PARTITION_OPTIONS_H_
 #define _MESH_PARTITION_OPTIONS_H_
 
-struct meshPartitionOptions
+class meshPartitionOptions
 {
+  private:
 
-//--General
-  int partitioner;                      // 1 - Chaco
-                                        // 2 - METIS
-  int num_partitions;
-
-  int renumber;
+  public:
 
-  bool createPartitionBoundaries;
-  bool createGhostCells;
-
-//--Chaco
+//--General
 
-  int global_method;                    // 1 - Multilevel-KL
-                                        // 2 - Spectral
-                                        // 3 - Inertial (Disabled)
-                                        // 4 - Linear
-                                        // 5 - Random
-                                        // 6 - Scattered
-  int architecture;                     // 0   - hypercube
-                                        // 1-3 - dimensional mesh parallel
-  int ndims_tot;                        // Number of hypercube dimensions */
-  int mesh_dims[3];                     // Number of partitions in each
-                                        // dimension for a dimensional
-                                        // architecture
-  double *goal;                         // Not-implemented (Weights vertices in
-                                        // a set)
-  int local_method;                     // 1 - Kernighan-Lin
-                                        // 2 - None
-  int rqi_flag;                         // Eigensolver for spectral partitioner
-                                        // 0 - Lanczos
-                                        // 1 - Multilevel RQI/Symmlq
-  int vmax;                             // Number of vertices in the coarsest
-                                        // graph.
-  int ndims;                            // Number of divisions at each stage
-                                        // 1 - Bisection
-                                        // 2 - Quadrisection
-                                        // 3 - Octasection
-  double eigtol;                        // Tolerance to the eigensolver
-  long seed;                            // RNG seed
+    int partitioner;                      // 1 - Chaco
+                                          // 2 - METIS
+    int num_partitions;
+     
+    int ncon;                             // Number of constraints/different weights
 
-  // parameters
-  int refine_partition;                 // Refine partitions at each level
-  int internal_vertices;                // Increase internal vertices
-  int refine_map;                       // Refine processor mapping
-  int terminal_propogation;             // Run with terminal propogation
+    int renumber;
 
-//--METIS
+    bool createPartitionBoundaries;
+    bool createGhostCells;
 
-  int algorithm;
-                                        // 1 - Recursive
-                                        // 2 - K-way
-                                        // 3 - Nodal weight
-  int edge_matching;                    // 1 - Random matching
-                                        // 2 - Heavy-edge matching
-                                        // 3 - Sorted heavy-edge matching
-  int refine_algorithm;                 // 1 - Random boundary refinement
-                                        // 2 - Greedy boundary refinement
-                                        // 3 - Random boundary refinement (with
-                                        //     minimization of connectivity
-                                        //     along sub-domains)
-  int partitionByExtrusion;            // if true, all extruded elements belong
-                                       // to the same partition as the source element
+//--Chaco
 
-  // element weights for load-balancing (currently used in METIS algorithm 3) 
+    int global_method;                    // 1 - Multilevel-KL
+                                          // 2 - Spectral
+                                          // 3 - Inertial (Disabled)
+                                          // 4 - Linear
+                                          // 5 - Random
+                                          // 6 - Scattered
+    int architecture;                     // 0   - hypercube
+                                          // 1-3 - dimensional mesh parallel
+    int ndims_tot;                        // Number of hypercube dimensions */
+    int mesh_dims[3];                     // Number of partitions in each
+                                          // dimension for a dimensional
+                                          // architecture
+    double *goal;                         // Not-implemented (Weights vertices in
+                                          // a set)
+    int local_method;                     // 1 - Kernighan-Lin
+                                          // 2 - None
+    int rqi_flag;                         // Eigensolver for spectral partitioner
+                                          // 0 - Lanczos
+                                          // 1 - Multilevel RQI/Symmlq
+    int vmax;                             // Number of vertices in the coarsest
+                                          // graph.
+    int ndims;                            // Number of divisions at each stage
+                                          // 1 - Bisection
+                                          // 2 - Quadrisection
+                                          // 3 - Octasection
+    double eigtol;                        // Tolerance to the eigensolver
+    long seed;                            // RNG seed
   
-  int triWeight;  
-  int quaWeight;
-  int tetWeight;
-  int priWeight;
-  int pyrWeight;
-  int hexWeight;
-
-//--NODAL WEIGHT
-  std::vector<int> nodalWeights;
-
-//--Constructor
+    // parameters
+    int refine_partition;                 // Refine partitions at each level
+    int internal_vertices;                // Increase internal vertices
+    int refine_map;                       // Refine processor mapping
+    int terminal_propogation;             // Run with terminal propogation
+  
+  //--METIS
+  
+    int algorithm;
+                                          // 1 - Recursive
+                                          // 2 - K-way
+                                          // 3 - Nodal weight
+                                          // 4 - Multi-Constrained Recursive
+                                          // 5 - Multi-Constrained K-way
+    int edge_matching;                    // 1 - Random matching
+                                          // 2 - Heavy-edge matching
+                                          // 3 - Sorted heavy-edge matching
+    int refine_algorithm;                 // 1 - Random boundary refinement
+                                          // 2 - Greedy boundary refinement
+                                          // 3 - Random boundary refinement (with
+                                          //     minimization of connectivity
+                                          //     along sub-domains)
+    int partitionByExtrusion;            // if true, all extruded elements belong
+                                         // to the same partition as the source element
+  
+    // element weights for load-balancing (currently used in METIS algorithm 3) 
+    
+    int triWeight;  
+    int quaWeight;
+    int tetWeight;
+    int priWeight;
+    int pyrWeight;
+    int hexWeight;
+  
+  //--NODAL WEIGHT
+    std::vector<int> nodalWeights;
 
-  meshPartitionOptions()
-  {
-    setDefaults();
-  }
+    std::map<int, std::vector<int> > weightMap;
+  
+  
+  public:
+  
+  //--Constructor
 
-//--Default values
+    meshPartitionOptions();
 
-  void setDefaults()
-  {
-    partitioner = 2;
-    num_partitions = 1;
-    renumber = 0;
-    global_method = 1;
-    architecture = 1;
-    ndims_tot = 2;
-    mesh_dims[0] = 4;
-    mesh_dims[1] = 1;
-    mesh_dims[2] = 1;
-    goal = 0;
-    local_method = 1;
-    rqi_flag = 1;
-    vmax = 250;
-    ndims = 1;
-    eigtol = 1.E-3;
-    seed = 7654321L;
-    refine_partition = false;
-    internal_vertices = false;
-    refine_map = true;
-    terminal_propogation = false;
-    algorithm = 1;
-    edge_matching = 3;
-    refine_algorithm = 3;
-    createPartitionBoundaries = true;
-    createGhostCells = true;
-    partitionByExtrusion =false;
-    triWeight = 1;
-    quaWeight = 1;
-    tetWeight = 1;
-    priWeight = 1;
-    pyrWeight = 1;
-    hexWeight = 1;
-    
-  }
+  //--Default values
 
+    void setDefaults();
+    void setAlgorithm(int algo);
+    void setNumOfPartitions(int numPart){num_partitions=numPart;};
+    int getNumOfPartitions(){return num_partitions;};
+    void setNumConstraints(int numConst){ncon=numConst;};
+    int getNumConstraints(){return ncon;};
+    void partition(GModel *model);
+    std::map<int,  std::vector<int> >& getWeightMap(){return weightMap;};
+    static void registerBindings(binding *b);
 };
 
 #endif