Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp	(revision 23583)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp	(revision 23584)
@@ -9,9 +9,32 @@
 #include "./ModelProcessorx.h"
 
+#define MAXCONNECTIVITY 5
+
+bool IsNodeInRank(int* nodes_ranks,int* nodes_proc_count,int nid,int rank){/*{{{*/
+
+	/*See if node is already in partition*/
+	for(int k=0;k<nodes_proc_count[nid];k++){
+		if(nodes_ranks[MAXCONNECTIVITY*nid+k] == rank) return true;
+	}
+
+	return false;
+}/*}}}*/
+int  NodeMasterRank(int* nodes_ranks,int nid){/*{{{*/
+	return nodes_ranks[MAXCONNECTIVITY*nid+0];
+}/*}}}*/
+void AddNodeToRank(int* nodes_ranks,int* nodes_proc_count,int nid,int rank){/*{{{*/
+
+	/*See if node is already in partition, return if this is the case*/
+	if(IsNodeInRank(nodes_ranks,nodes_proc_count,nid,rank)) return;
+
+	/*This rank has not been marked for this node just yet so go ahead and add it*/
+	nodes_ranks[MAXCONNECTIVITY*nid+nodes_proc_count[nid]] = rank;
+	nodes_proc_count[nid]++;
+}/*}}}*/
+
 void CreateNodes(Nodes* nodes, IoModel* iomodel,int analysis,int finite_element,int approximation,int* approximations){
 
 	/*Intermediaries*/
 	int        numnodes;
-	const int  MAXCONNECTIVITY = 5;
 	int        element_numnodes;
 	int        element_node_ids[40] = {0};
@@ -22,7 +45,6 @@
 	int*     epart = iomodel->epart;
 
-
 	/*Determine how many nodes we have in total*/
-	/*Define nodes sids for each element*/
+	/*{{{*/
 	if(iomodel->meshelementtype==TriaEnum){
 		switch(finite_element){
@@ -125,4 +147,5 @@
 		_error_("mesh elements "<< EnumToStringx(iomodel->meshelementtype) <<" not supported yet");
 	}
+	/*}}}*/
 
 	/*create matrix that keeps track of all ranks that have node i, and initialize as -1 (Common to all CPUs)*/
@@ -132,8 +155,4 @@
 	/*For all nodes, count how many cpus have node i (initialize with 0)*/
 	int* nodes_proc_count = xNewZeroInit<int>(numnodes);
-
-	/*Create vector of size total numnodes, initialized with -1, that will keep track of local ids*/
-	int* nodes_lids  = xNew<int>(numnodes);
-	for(int i=0;i<numnodes;i++) nodes_lids[i] = -1;
 
 	/*Create vector of approximation per node (used for FS: vel or pressure)*/
@@ -147,5 +166,5 @@
 
 	/*Go through all elements and mark all vertices for all partitions*/
-	int  lid = 0;
+	/*{{{*/
 	for(int i=0;i<iomodel->numberofelements;i++){
 
@@ -478,34 +497,12 @@
 		}
 
+		/*Add rank epart[i] for all nodes belonging to this element*/
 		for(int j=0;j<element_numnodes;j++){
 			int nid = element_node_ids[j]; _assert_(nid<numnodes);
-
-			/*See if it has already been marked*/
-			bool found = false;
-			for(int k=0;k<nodes_proc_count[nid];k++){
-				if(nodes_ranks[MAXCONNECTIVITY*nid+k] == epart[i]){
-					found = true;
-					break; 
-				}
-			}
-
-			/*On go below if this vertex has not been seen yet in this partition*/
-			if(!found){
-				/*This rank has not been marked for this vertex just yet so go ahead and mark it*/
-				nodes_ranks[MAXCONNECTIVITY*nid+nodes_proc_count[nid]] = epart[i];
-				nodes_proc_count[nid]++;
-
-				/*Keep track of local ids!*/
-				if(epart[i]==my_rank){
-					nodes_lids[nid] = lid;
-					lid++;
-				}
-
-				/*Make sure we don't go too far in the table*/
-				if(nodes_proc_count[nid]>MAXCONNECTIVITY) _error_("need to increase MAXCONNECTIVITY (this is hard coded, contact ISSM developer)");
-			}
-		}
-	}
-	if(finite_element==P1DGEnum){/* Special case for DG...{{{*/
+			AddNodeToRank(nodes_ranks,nodes_proc_count,nid,epart[i]);
+		}
+	}
+	/*}}}*/
+	if(finite_element==P1DGEnum){/*Special case for DG...{{{*/
 		int node_list[4];
 		if(!iomodel->domaintype==Domain2DhorizontalEnum) _error_("not implemented yet");
@@ -533,39 +530,6 @@
 					for(int j=0;j<4;j++){
 						int  nid = node_list[j];
-						bool found = false;
-						for(int k=0;k<nodes_proc_count[nid];k++){
-							if(nodes_ranks[MAXCONNECTIVITY*nid+k] == epart[e1]){
-								found = true;
-								break; 
-							}
-						}
-						if(!found){
-							nodes_ranks[MAXCONNECTIVITY*nid+nodes_proc_count[nid]] = epart[e1];
-							nodes_proc_count[nid]++;
-							if(epart[e1]==my_rank){
-								nodes_lids[nid] = lid;
-								lid++;
-							}
-							if(nodes_proc_count[nid]>MAXCONNECTIVITY) _error_("need to increase MAXCONNECTIVITY (this is hard coded, contact ISSM developer)");
-						}
-					}
-					for(int j=0;j<4;j++){
-						int  nid = node_list[j];
-						bool found = false;
-						for(int k=0;k<nodes_proc_count[nid];k++){
-							if(nodes_ranks[MAXCONNECTIVITY*nid+k] == epart[e2]){
-								found = true;
-								break; 
-							}
-						}
-						if(!found){
-							nodes_ranks[MAXCONNECTIVITY*nid+nodes_proc_count[nid]] = epart[e2];
-							nodes_proc_count[nid]++;
-							if(epart[e2]==my_rank){
-								nodes_lids[nid] = lid;
-								lid++;
-							}
-							if(nodes_proc_count[nid]>MAXCONNECTIVITY) _error_("need to increase MAXCONNECTIVITY (this is hard coded, contact ISSM developer)");
-						}
+						AddNodeToRank(nodes_ranks,nodes_proc_count,nid,epart[e1]);
+						AddNodeToRank(nodes_ranks,nodes_proc_count,nid,epart[e2]);
 					}
 				}
@@ -573,5 +537,31 @@
 		}
 	}/*}}}*/
-	//if(my_rank==0) printarray(nodes_ranks,numnodes,MAXCONNECTIVITY);
+	/*Vertex pairing for stressbalance{{{*/
+	if(analysis==StressbalanceAnalysisEnum || analysis==StressbalanceVerticalAnalysisEnum){
+		int *vertex_pairing = NULL;
+		int  numvertex_pairing;
+		iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,"md.stressbalance.vertex_pairing");
+		_assert_(numvertex_pairing==0 || finite_element==P1Enum);
+		for(int i=0;i<numvertex_pairing;i++){
+			int nid1 = vertex_pairing[2*i+0]-1;
+			int nid2 = vertex_pairing[2*i+1]-1;
+			for(int j=0;j<nodes_proc_count[nid1];j++) AddNodeToRank(nodes_ranks,nodes_proc_count,nid2,nodes_ranks[MAXCONNECTIVITY*nid1+j]);
+			for(int j=0;j<nodes_proc_count[nid2];j++) AddNodeToRank(nodes_ranks,nodes_proc_count,nid1,nodes_ranks[MAXCONNECTIVITY*nid2+j]);
+		}
+		xDelete<int>(vertex_pairing);
+	}
+	/*}}}*/
+
+	/*Create vector of size total numnodes, initialized with -1, that will keep track of local ids*/
+	int  lid = 0;
+	int* nodes_lids  = xNew<int>(numnodes);
+	for(int i=0;i<numnodes;i++){
+		if(IsNodeInRank(nodes_ranks,nodes_proc_count,i,my_rank)){
+			nodes_lids[i] = lid++;
+		}
+		else{
+			nodes_lids[i] = -1;
+		}
+	}
 
 	/*Now, Count how many clones we have with other partitions*/
@@ -687,2 +677,3 @@
 	nodes->common_recv_ids=common_recv_ids;
 }
+