Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp	(revision 740)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp	(revision 741)
@@ -51,5 +51,5 @@
 	#ifdef _PARALLEL_
 	/*keep only this partition's nodes:*/
-	if((model->npart[i]==1)){
+	if((model->my_grids[i])){
 	#endif
 
@@ -76,5 +76,5 @@
 
 	#ifdef _PARALLEL_
-	} //if((npart[i]==1))
+	} //if((my_grids[i]))
 	#endif
 	}
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp	(revision 740)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp	(revision 741)
@@ -97,4 +97,14 @@
 	int node_numdofs;
 
+	#ifdef _PARALLEL_
+	/*Metis partitioning: */
+	int  range;
+	Vec  gridborder=NULL;
+	int  my_numgrids;
+	int* all_numgrids=NULL;
+	int  gridcount;
+	int  count;
+	#endif
+
 	/*First create the elements, nodes and material properties: */
 	elements  = new DataSet(ElementsEnum());
@@ -134,9 +144,65 @@
 	xfree((void**)&model->elements2d);
 
-			
-	#endif
-
-	/*Hutter elements can be partitioned using npart, instead of epart for more classic tria and penta elements. The reason is 
-	 * that each hutter elements either lies on a node (in 2d), or a pair of vertically juxtaposed nodes (in 3d): */
+	/*Used later on: */
+	my_grids=(int*)xcalloc(model->numberofnodes,sizeof(int));
+	#endif
+
+	#ifdef _PARALLEL_
+	if(strcmp(model->meshtype,"2d")==0){
+		ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
+		for (i=0;i<model->numberofelements;i++){
+			if(my_rank==epart[i]){ 
+				/*Now that we are here, we can also start building the list of grids belonging to this node partition: we use 
+				 *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 
+				 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids 
+				 will hold which grids belong to this partition*/
+
+				my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
+				my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
+				my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
+			}
+		}
+	}
+	else{
+		ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
+		for (i=0;i<model->numberofelements;i++){
+			if(my_rank==epart[i]){ 
+				my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
+				my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
+				my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
+				my_grids[(int)*(model->elements+elements_width*i+3)-1]=1;
+				my_grids[(int)*(model->elements+elements_width*i+4)-1]=1;
+				my_grids[(int)*(model->elements+elements_width*i+5)-1]=1;
+			}
+		}
+	}
+
+	/*From the element partitioning, we can determine which grids are on the inside of this cpu's 
+	 *element partition, and which are on its border with other nodes:*/
+	gridborder=NewVec(model->numberofnodes);
+
+	for (i=0;i<model->numberofnodes;i++){
+		if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);
+	}
+	VecAssemblyBegin(gridborder);
+	VecAssemblyEnd(gridborder);
+
+	#ifdef _DEBUG_
+	VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD);
+	#endif
+	
+	VecToMPISerial(&my_bordergrids,gridborder);
+
+	#ifdef _DEBUG_
+	if(my_rank==0){
+		for (i=0;i<model->numberofnodes;i++){
+			printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);
+		}
+	}
+	#endif
+	#endif
+
+	/*Hutter elements can be partitioned using epart, even if
+	 * each hutter elements either lies on a node (in 2d), or a pair of vertically juxtaposed nodes (in 3d): */
 
 	/*Fetch data temporarily needed: */
@@ -158,5 +224,5 @@
 		#ifdef _PARALLEL_
 		/*keep only this partition's nodes:*/
-		if(my_rank==npart[i]){
+		if(my_grids[i]){
 		#endif
 
@@ -200,5 +266,5 @@
 		#ifdef _PARALLEL_
 		/*keep only this partition's nodes:*/
-		if(my_rank==npart[i]){
+		if(my_grids[i]){
 		#endif
 
@@ -253,4 +319,5 @@
 
 	/*Free data: */
+	xfree((void**)&model->elements);
 	xfree((void**)&model->gridonhutter);
 	xfree((void**)&model->thickness);
@@ -305,10 +372,21 @@
 	#ifdef _PARALLEL_
 	/*keep only this partition's nodes:*/
-	if(npart[i]==my_rank){
+	if(my_grids[i]){
 	#endif
 
 		node_id=i+1; //matlab indexing
 			
-		node_partitionborder=0;
+
+		#ifdef _PARALLEL_
+		if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
+			node_partitionborder=1;
+		}
+		else{
+			node_partitionborder=0;
+		}
+		#else
+			node_partitionborder=0;
+		#endif
+
 
 		node_x[0]=model->x[i];
@@ -368,5 +446,14 @@
 
 	/*Keep partitioning information into model*/
-	xfree((void**)&epart);
+	model->epart=epart;
+	model->my_grids=my_grids;
+	model->my_bordergrids=my_bordergrids;
+
+	/*Keep partitioning information into model*/
+	#ifdef _PARALLEL_
+	xfree((void**)&all_numgrids);
+	xfree((void**)&npart);
+	VecFree(&gridborder);
+	#endif
 	model->npart=npart;