Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 23166)
+++ /issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 23167)
@@ -130,4 +130,6 @@
 
 	this->my_elements=NULL;
+	this->my_faces=NULL;
+	this->my_edges=NULL;
 	this->my_vertices=NULL;
 
@@ -192,4 +194,6 @@
 	/*Initialize permanent data: */
 	this->my_elements = NULL;
+	this->my_faces = NULL;
+	this->my_edges = NULL;
 	this->my_vertices = NULL;
 
@@ -236,4 +240,6 @@
 
 	xDelete<bool>(this->my_elements);
+	xDelete<bool>(this->my_faces);
+	xDelete<bool>(this->my_edges);
 	xDelete<int>(this->my_vertices);
 
Index: /issm/trunk-jpl/src/c/classes/IoModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/IoModel.h	(revision 23166)
+++ /issm/trunk-jpl/src/c/classes/IoModel.h	(revision 23167)
@@ -64,4 +64,6 @@
 		/*Partitioning*/
 		bool *my_elements;
+		bool *my_faces;
+		bool *my_edges;
 		int  *my_vertices;
 
Index: /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp	(revision 23166)
+++ /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp	(revision 23167)
@@ -62,6 +62,4 @@
 	/*Higher-order finite elements*/
 	int   v1,v2;
-	bool *my_edges = NULL;
-	bool *my_faces = NULL;
 	bool *boundaryedge = NULL;
 
@@ -84,20 +82,18 @@
 		case P1xP3Enum:
 		case P1xP4Enum:
-			EdgesPartitioning(&my_edges,iomodel);
+			EdgesPartitioning(iomodel);
 			break;
 		case P2xP1Enum:
-			EdgesPartitioning(&my_edges,iomodel);
+			EdgesPartitioning(iomodel);
 			break;
 		case P2Enum:
-			EdgesPartitioning(&my_edges,iomodel);
-	      if(iomodel->meshelementtype==PentaEnum){
-				FacesPartitioning(&my_faces,iomodel);
-			}
+			EdgesPartitioning(iomodel);
+	      if(iomodel->meshelementtype==PentaEnum) FacesPartitioning(iomodel);
 			EdgeOnBoundaryFlags(&boundaryedge,iomodel);
 			break;
 		case P2bubbleEnum:
-			EdgesPartitioning(&my_edges,iomodel);
+			EdgesPartitioning(iomodel);
 			if(iomodel->meshelementtype==PentaEnum){
-				FacesPartitioning(&my_faces,iomodel);
+				FacesPartitioning(iomodel);
 			}
 			EdgeOnBoundaryFlags(&boundaryedge,iomodel);
@@ -110,6 +106,6 @@
 			break;
 		case P2xP4Enum:
-			EdgesPartitioning(&my_edges,iomodel);
-			FacesPartitioning(&my_faces,iomodel);
+			EdgesPartitioning(iomodel);
+			FacesPartitioning(iomodel);
 			break;
 		default:
@@ -140,5 +136,5 @@
 				}
 				for(i=0;i<iomodel->numberofedges;i++){
-					if(my_edges[i] && boundaryedge[i]){
+					if(iomodel->my_edges[i] && boundaryedge[i]){
 						v1 = iomodel->edges[3*i+0]-1;
 						v2 = iomodel->edges[3*i+1]-1;
@@ -153,5 +149,5 @@
 					for(i=0;i<iomodel->numberoffaces;i++){
 						if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
-							if(my_faces[i]){
+							if(iomodel->my_faces[i]){
 								numfacevertices = iomodel->faces[i*iomodel->facescols+3];
 								value=0.;
@@ -181,5 +177,5 @@
 				}
 				for(i=0;i<iomodel->numberofedges;i++){
-					if(my_edges[i] && boundaryedge[i]){
+					if(iomodel->my_edges[i] && boundaryedge[i]){
 						v1 = iomodel->edges[3*i+0]-1;
 						v2 = iomodel->edges[3*i+1]-1;
@@ -194,5 +190,5 @@
 					for(i=0;i<iomodel->numberoffaces;i++){
 						if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
-							if(my_faces[i]){
+							if(iomodel->my_faces[i]){
 								numfacevertices = iomodel->faces[i*iomodel->facescols+3];
 								value=0.;
@@ -238,5 +234,5 @@
 				for(i=0;i<iomodel->numberofedges;i++){
 					if(iomodel->edges[i*3+2]==2){/*Vertical edges*/
-						if(my_edges[i]){
+						if(iomodel->my_edges[i]){
 							v1 = iomodel->edges[3*i+0]-1;
 							v2 = iomodel->edges[3*i+1]-1;
@@ -253,5 +249,5 @@
 					}
 					if(iomodel->edges[i*3+2]==1){/*Horizontal edges*/
-						if(my_edges[i]){
+						if(iomodel->my_edges[i]){
 							v1 = iomodel->edges[3*i+0]-1;
 							v2 = iomodel->edges[3*i+1]-1;
@@ -266,5 +262,5 @@
 				for(i=0;i<iomodel->numberoffaces;i++){
 					if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
-						if(my_faces[i]){
+						if(iomodel->my_faces[i]){
 							numfacevertices = iomodel->faces[i*iomodel->facescols+3];
 							value=0.;
@@ -332,5 +328,5 @@
 				for(i=0;i<iomodel->numberofedges;i++){
 					if(iomodel->edges[i*3+2]==2){
-						if(my_edges[i]){
+						if(iomodel->my_edges[i]){
 							v1 = iomodel->edges[3*i+0]-1;
 							v2 = iomodel->edges[3*i+1]-1;
@@ -355,5 +351,5 @@
 				for(i=0;i<iomodel->numberofedges;i++){
 					if(iomodel->edges[i*3+2]==2){
-						if(my_edges[i]){
+						if(iomodel->my_edges[i]){
 							v1 = iomodel->edges[3*i+0]-1;
 							v2 = iomodel->edges[3*i+1]-1;
@@ -380,5 +376,5 @@
 				for(i=0;i<iomodel->numberofedges;i++){
 					if(iomodel->edges[i*3+2]==2){/*Vertical edges*/
-						if(my_edges[i]){
+						if(iomodel->my_edges[i]){
 							v1 = iomodel->edges[3*i+0]-1;
 							v2 = iomodel->edges[3*i+1]-1;
@@ -407,5 +403,5 @@
 				for(i=0;i<iomodel->numberofedges;i++){
 					if(iomodel->edges[i*3+2]!=2){
-						if(my_edges[i]){
+						if(iomodel->my_edges[i]){
 							v1 = iomodel->edges[3*i+0]-1;
 							v2 = iomodel->edges[3*i+1]-1;
@@ -472,5 +468,5 @@
 				}
 				for(i=0;i<iomodel->numberofedges;i++){
-					if(my_edges[i]){
+					if(iomodel->my_edges[i]){
 						v1 = iomodel->edges[3*i+0]-1;
 						v2 = iomodel->edges[3*i+1]-1;
@@ -511,5 +507,5 @@
 				for(i=0;i<iomodel->numberofedges;i++){
 					if(iomodel->edges[i*3+2]==2){
-						if(my_edges[i]){
+						if(iomodel->my_edges[i]){
 							v1 = iomodel->edges[3*i+0]-1;
 							v2 = iomodel->edges[3*i+1]-1;
@@ -551,5 +547,5 @@
 				for(i=0;i<iomodel->numberofedges;i++){
 					if(iomodel->edges[i*3+2]==2){
-						if(my_edges[i]){
+						if(iomodel->my_edges[i]){
 							v1 = iomodel->edges[3*i+0]-1;
 							v2 = iomodel->edges[3*i+1]-1;
@@ -601,5 +597,5 @@
 				for(i=0;i<iomodel->numberofedges;i++){
 					if(iomodel->edges[i*3+2]!=2){
-						if(my_edges[i]){
+						if(iomodel->my_edges[i]){
 							v1 = iomodel->edges[3*i+0]-1;
 							v2 = iomodel->edges[3*i+1]-1;
@@ -634,6 +630,4 @@
 	xDelete<IssmDouble>(times);
 	xDelete<IssmDouble>(values);
-	xDelete<bool>(my_edges);
-	xDelete<bool>(my_faces);
 	xDelete<bool>(boundaryedge);
 }/*}}}*/
@@ -649,6 +643,4 @@
 	/*Higher-order finite elements*/
 	int   v1,v2;
-	bool *my_edges = NULL;
-	bool *my_faces = NULL;
 	bool *boundaryedge = NULL;
 
@@ -669,23 +661,23 @@
 			break;
 		case P1xP2Enum:
-			EdgesPartitioning(&my_edges,iomodel);
+			EdgesPartitioning(iomodel);
 			break;
 		case P1xP3Enum:
-			EdgesPartitioning(&my_edges,iomodel);
+			EdgesPartitioning(iomodel);
 			break;
 		case P2xP1Enum:
-			EdgesPartitioning(&my_edges,iomodel);
+			EdgesPartitioning(iomodel);
 			break;
 		case P2Enum:
-			EdgesPartitioning(&my_edges,iomodel);
+			EdgesPartitioning(iomodel);
 			if(iomodel->meshelementtype==PentaEnum){
-				FacesPartitioning(&my_faces,iomodel);
+				FacesPartitioning(iomodel);
 			}
 			EdgeOnBoundaryFlags(&boundaryedge,iomodel);
 			break;
 		case P2bubbleEnum:
-			EdgesPartitioning(&my_edges,iomodel);
+			EdgesPartitioning(iomodel);
 			if(iomodel->meshelementtype==PentaEnum){
-				FacesPartitioning(&my_faces,iomodel);
+				FacesPartitioning(iomodel);
 			}
 			EdgeOnBoundaryFlags(&boundaryedge,iomodel);
@@ -698,6 +690,6 @@
 			break;
 		case P2xP4Enum:
-			EdgesPartitioning(&my_edges,iomodel);
-			FacesPartitioning(&my_faces,iomodel);
+			EdgesPartitioning(iomodel);
+			FacesPartitioning(iomodel);
 			break;
 		default:
@@ -729,5 +721,5 @@
 				for(i=0;i<iomodel->numberofedges;i++){
 					if(iomodel->edges[i*3+2]==2){
-						if(my_edges[i]){
+						if(iomodel->my_edges[i]){
 							v1 = iomodel->edges[3*i+0]-1;
 							v2 = iomodel->edges[3*i+1]-1;
@@ -752,5 +744,5 @@
 				for(i=0;i<iomodel->numberofedges;i++){
 					if(iomodel->edges[i*3+2]==2){
-						if(my_edges[i]){
+						if(iomodel->my_edges[i]){
 							v1 = iomodel->edges[3*i+0]-1;
 							v2 = iomodel->edges[3*i+1]-1;
@@ -780,6 +772,4 @@
 	xDelete<IssmDouble>(times);
 	xDelete<IssmDouble>(values);
-	xDelete<bool>(my_edges);
-	xDelete<bool>(my_faces);
 	xDelete<bool>(boundaryedge);
 }/*}}}*/
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp	(revision 23166)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp	(revision 23167)
@@ -15,6 +15,4 @@
 	int   numberoffaces,elementnbv;
 	int   id0 = iomodel->nodecounter;
-	bool *my_faces = NULL;
-	bool *my_edges = NULL;
 	bool *my_nodes = NULL;
 	Node *node     = NULL;
@@ -70,5 +68,5 @@
 
 		case P1xP2Enum:
-			EdgesPartitioning(&my_edges,iomodel);
+			EdgesPartitioning(iomodel);
 			for(i=0;i<iomodel->numberofvertices;i++){
 				if(iomodel->my_vertices[i]){
@@ -80,5 +78,5 @@
 			for(i=0;i<iomodel->numberofedges;i++){
 				if(iomodel->edges[i*3+2]==2){
-					if(my_edges[i]){
+					if(iomodel->my_edges[i]){
 						node = new Node(id0+iomodel->numberofvertices+i+1,counter+1,lid++,0,iomodel,analysis,approximation);
 						nodes->AddObject(node);
@@ -90,5 +88,5 @@
 
 		case P1xP3Enum:
-			EdgesPartitioning(&my_edges,iomodel);
+			EdgesPartitioning(iomodel);
 			for(i=0;i<iomodel->numberofvertices;i++){
 				if(iomodel->my_vertices[i]){
@@ -100,5 +98,5 @@
 			for(i=0;i<iomodel->numberofedges;i++){
 				if(iomodel->edges[i*3+2]==2){
-					if(my_edges[i]){
+					if(iomodel->my_edges[i]){
 						node = new Node(id0+iomodel->numberofvertices+2*i+1,counter+1,lid++,0,iomodel,analysis,approximation);
 						nodes->AddObject(node);
@@ -111,5 +109,5 @@
 			break;
 		case P1xP4Enum:
-			EdgesPartitioning(&my_edges,iomodel);
+			EdgesPartitioning(iomodel);
 			for(i=0;i<iomodel->numberofvertices;i++){
 				if(iomodel->my_vertices[i]){
@@ -120,5 +118,5 @@
 			for(i=0;i<iomodel->numberofedges;i++){
 				if(iomodel->edges[i*3+2]==2){/*Vertical edges*/
-					if(my_edges[i]){
+					if(iomodel->my_edges[i]){
 						node = new Node(id0+iomodel->numberofvertices+3*i+1,counter+1,lid++,0,iomodel,analysis,approximation);
 						nodes->AddObject(node);
@@ -134,5 +132,5 @@
 
 		case P2xP1Enum:
-			EdgesPartitioning(&my_edges,iomodel);
+			EdgesPartitioning(iomodel);
 			for(i=0;i<iomodel->numberofvertices;i++){
 				if(iomodel->my_vertices[i]){
@@ -144,5 +142,5 @@
 			for(i=0;i<iomodel->numberofedges;i++){
 				if(iomodel->edges[i*3+2]!=2){
-					if(my_edges[i]){
+					if(iomodel->my_edges[i]){
 						node = new Node(id0+iomodel->numberofvertices+i+1,counter+1,lid++,0,iomodel,analysis,approximation);
 						nodes->AddObject(node);
@@ -154,12 +152,12 @@
 
 		case P2Enum:
-			EdgesPartitioning(&my_edges,iomodel);
-			for(i=0;i<iomodel->numberofvertices;i++){
-				if(iomodel->my_vertices[i]){
-					nodes->AddObject(new Node(id0+i+1,i,lid++,i,iomodel,analysis,approximation));
-				}
-			}
-			for(i=0;i<iomodel->numberofedges;i++){
-				if(my_edges[i]){
+			EdgesPartitioning(iomodel);
+			for(i=0;i<iomodel->numberofvertices;i++){
+				if(iomodel->my_vertices[i]){
+					nodes->AddObject(new Node(id0+i+1,i,lid++,i,iomodel,analysis,approximation));
+				}
+			}
+			for(i=0;i<iomodel->numberofedges;i++){
+				if(iomodel->my_edges[i]){
 					nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,analysis,approximation));
 				}
@@ -167,8 +165,8 @@
 			id0 = id0+iomodel->numberofvertices+iomodel->numberofedges;
 	      if(iomodel->meshelementtype==PentaEnum){
-				FacesPartitioning(&my_faces,iomodel);
+				FacesPartitioning(iomodel);
 				for(i=0;i<iomodel->numberoffaces;i++){
 					if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
-						if(my_faces[i]){
+						if(iomodel->my_faces[i]){
 							node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,iomodel,analysis,approximation);
 							nodes->AddObject(node);
@@ -185,12 +183,12 @@
 			break;
 		case P2bubbleEnum:
-			EdgesPartitioning(&my_edges,iomodel);
-			for(i=0;i<iomodel->numberofvertices;i++){
-				if(iomodel->my_vertices[i]){
-					nodes->AddObject(new Node(id0+i+1,i,lid++,i,iomodel,analysis,approximation));
-				}
-			}
-			for(i=0;i<iomodel->numberofedges;i++){
-				if(my_edges[i]){
+			EdgesPartitioning(iomodel);
+			for(i=0;i<iomodel->numberofvertices;i++){
+				if(iomodel->my_vertices[i]){
+					nodes->AddObject(new Node(id0+i+1,i,lid++,i,iomodel,analysis,approximation));
+				}
+			}
+			for(i=0;i<iomodel->numberofedges;i++){
+				if(iomodel->my_edges[i]){
 					nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,analysis,approximation));
 				}
@@ -198,8 +196,8 @@
 			id0 = id0+iomodel->numberofvertices+iomodel->numberofedges;
 			if(iomodel->meshelementtype==PentaEnum){
-				FacesPartitioning(&my_faces,iomodel);
+				FacesPartitioning(iomodel);
 				for(i=0;i<iomodel->numberoffaces;i++){
 					if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
-						if(my_faces[i]){
+						if(iomodel->my_faces[i]){
 							node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,iomodel,analysis,approximation);
 							nodes->AddObject(node);
@@ -222,6 +220,6 @@
 			break;
 		case P2xP4Enum:
-			EdgesPartitioning(&my_edges,iomodel);
-			FacesPartitioning(&my_faces,iomodel);
+			EdgesPartitioning(iomodel);
+			FacesPartitioning(iomodel);
 			for(i=0;i<iomodel->numberofvertices;i++){
 				if(iomodel->my_vertices[i]){
@@ -232,5 +230,5 @@
 			for(i=0;i<iomodel->numberofedges;i++){
 				if(iomodel->edges[i*3+2]==2){/*Vertical edges*/
-					if(my_edges[i]){
+					if(iomodel->my_edges[i]){
 						node = new Node(id0+iomodel->numberofvertices+3*i+1,counter+1,lid++,0,iomodel,analysis,approximation);
 						nodes->AddObject(node);
@@ -243,5 +241,5 @@
 				}
 				else if(iomodel->edges[i*3+2]==1){/*Horizontal edges*/
-					if(my_edges[i]){
+					if(iomodel->my_edges[i]){
 						node = new Node(id0+iomodel->numberofvertices+3*i+1,counter+1,lid++,0,iomodel,analysis,approximation);
 						nodes->AddObject(node);
@@ -256,5 +254,5 @@
 			for(i=0;i<iomodel->numberoffaces;i++){
 				if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
-					if(my_faces[i]){
+					if(iomodel->my_faces[i]){
 						node = new Node(id0+3*i+1,counter+1,lid++,0,iomodel,analysis,approximation);
 						nodes->AddObject(node);
@@ -356,12 +354,12 @@
 			_assert_(approximation==FSApproximationEnum);
 			/*P2 velocity*/
-			EdgesPartitioning(&my_edges,iomodel);
-			for(i=0;i<iomodel->numberofvertices;i++){
-				if(iomodel->my_vertices[i]){
-					nodes->AddObject(new Node(id0+i+1,i,lid++,i,iomodel,analysis,FSvelocityEnum));
-				}
-			}
-			for(i=0;i<iomodel->numberofedges;i++){
-				if(my_edges[i]){
+			EdgesPartitioning(iomodel);
+			for(i=0;i<iomodel->numberofvertices;i++){
+				if(iomodel->my_vertices[i]){
+					nodes->AddObject(new Node(id0+i+1,i,lid++,i,iomodel,analysis,FSvelocityEnum));
+				}
+			}
+			for(i=0;i<iomodel->numberofedges;i++){
+				if(iomodel->my_edges[i]){
 					nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,analysis,FSvelocityEnum));
 				}
@@ -369,8 +367,8 @@
 			id0 = id0+iomodel->numberofvertices+iomodel->numberofedges;
 	      if(iomodel->meshelementtype==PentaEnum){
-				FacesPartitioning(&my_faces,iomodel);
+				FacesPartitioning(iomodel);
 				for(i=0;i<iomodel->numberoffaces;i++){
 					if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
-						if(my_faces[i]){
+						if(iomodel->my_faces[i]){
 							node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,iomodel,analysis,FSvelocityEnum);
 							nodes->AddObject(node);
@@ -403,12 +401,12 @@
 			_assert_(approximation==FSApproximationEnum);
 			/*P2 velocity*/
-			EdgesPartitioning(&my_edges,iomodel);
-			for(i=0;i<iomodel->numberofvertices;i++){
-				if(iomodel->my_vertices[i]){
-					nodes->AddObject(new Node(id0+i+1,i,lid++,i,iomodel,analysis,FSvelocityEnum));
-				}
-			}
-			for(i=0;i<iomodel->numberofedges;i++){
-				if(my_edges[i]){
+			EdgesPartitioning(iomodel);
+			for(i=0;i<iomodel->numberofvertices;i++){
+				if(iomodel->my_vertices[i]){
+					nodes->AddObject(new Node(id0+i+1,i,lid++,i,iomodel,analysis,FSvelocityEnum));
+				}
+			}
+			for(i=0;i<iomodel->numberofedges;i++){
+				if(iomodel->my_edges[i]){
 					nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,analysis,FSvelocityEnum));
 				}
@@ -416,8 +414,8 @@
 			id0 = id0+iomodel->numberofvertices+iomodel->numberofedges;
 	      if(iomodel->meshelementtype==PentaEnum){
-				FacesPartitioning(&my_faces,iomodel);
+				FacesPartitioning(iomodel);
 				for(i=0;i<iomodel->numberoffaces;i++){
 					if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
-						if(my_faces[i]){
+						if(iomodel->my_faces[i]){
 							node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,iomodel,analysis,FSvelocityEnum);
 							nodes->AddObject(node);
@@ -437,6 +435,6 @@
 		case OneLayerP4zEnum:
 			_assert_(approximation==FSApproximationEnum);
-			EdgesPartitioning(&my_edges,iomodel);
-			FacesPartitioning(&my_faces,iomodel);
+			EdgesPartitioning(iomodel);
+			FacesPartitioning(iomodel);
 			/*P2xP4 velocity*/
 			for(i=0;i<iomodel->numberofvertices;i++){
@@ -448,5 +446,5 @@
 			for(i=0;i<iomodel->numberofedges;i++){
 				if(iomodel->edges[i*3+2]==2){/*Vertical edges*/
-					if(my_edges[i]){
+					if(iomodel->my_edges[i]){
 						node = new Node(id0+iomodel->numberofvertices+3*i+1,counter+1,lid++,0,iomodel,analysis,FSvelocityEnum);
 						nodes->AddObject(node);
@@ -459,5 +457,5 @@
 				}
 				else if(iomodel->edges[i*3+2]==1){/*Horizontal edges*/
-					if(my_edges[i]){
+					if(iomodel->my_edges[i]){
 						node = new Node(id0+iomodel->numberofvertices+3*i+1,counter+1,lid++,0,iomodel,analysis,FSvelocityEnum);
 						nodes->AddObject(node);
@@ -472,5 +470,5 @@
 			for(i=0;i<iomodel->numberoffaces;i++){
 				if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
-					if(my_faces[i]){
+					if(iomodel->my_faces[i]){
 						node = new Node(id0+3*i+1,counter+1,lid++,0,iomodel,analysis,FSvelocityEnum);
 						nodes->AddObject(node);
@@ -502,12 +500,12 @@
 			_assert_(approximation==FSApproximationEnum);
 			/*P2b velocity*/
-			EdgesPartitioning(&my_edges,iomodel);
-			for(i=0;i<iomodel->numberofvertices;i++){
-				if(iomodel->my_vertices[i]){
-					nodes->AddObject(new Node(id0+i+1,i,lid++,i,iomodel,analysis,FSvelocityEnum));
-				}
-			}
-			for(i=0;i<iomodel->numberofedges;i++){
-				if(my_edges[i]){
+			EdgesPartitioning(iomodel);
+			for(i=0;i<iomodel->numberofvertices;i++){
+				if(iomodel->my_vertices[i]){
+					nodes->AddObject(new Node(id0+i+1,i,lid++,i,iomodel,analysis,FSvelocityEnum));
+				}
+			}
+			for(i=0;i<iomodel->numberofedges;i++){
+				if(iomodel->my_edges[i]){
 					nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,analysis,FSvelocityEnum));
 				}
@@ -515,8 +513,8 @@
 			id0 = id0+iomodel->numberofvertices+iomodel->numberofedges;
 			if(iomodel->meshelementtype==PentaEnum){
-				FacesPartitioning(&my_faces,iomodel);
+				FacesPartitioning(iomodel);
 				for(i=0;i<iomodel->numberoffaces;i++){
 					if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
-						if(my_faces[i]){
+						if(iomodel->my_faces[i]){
 							node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,iomodel,analysis,FSvelocityEnum);
 							nodes->AddObject(node);
@@ -558,12 +556,12 @@
 			_assert_(approximation==FSApproximationEnum);
 			/*P2b velocity*/
-			EdgesPartitioning(&my_edges,iomodel);
-			for(i=0;i<iomodel->numberofvertices;i++){
-				if(iomodel->my_vertices[i]){
-					nodes->AddObject(new Node(id0+i+1,i,lid++,i,iomodel,analysis,FSvelocityEnum));
-				}
-			}
-			for(i=0;i<iomodel->numberofedges;i++){
-				if(my_edges[i]){
+			EdgesPartitioning(iomodel);
+			for(i=0;i<iomodel->numberofvertices;i++){
+				if(iomodel->my_vertices[i]){
+					nodes->AddObject(new Node(id0+i+1,i,lid++,i,iomodel,analysis,FSvelocityEnum));
+				}
+			}
+			for(i=0;i<iomodel->numberofedges;i++){
+				if(iomodel->my_edges[i]){
 					nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,analysis,FSvelocityEnum));
 				}
@@ -571,8 +569,8 @@
 			id0 = id0+iomodel->numberofvertices+iomodel->numberofedges;
 			if(iomodel->meshelementtype==PentaEnum){
-				FacesPartitioning(&my_faces,iomodel);
+				FacesPartitioning(iomodel);
 				for(i=0;i<iomodel->numberoffaces;i++){
 					if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
-						if(my_faces[i]){
+						if(iomodel->my_faces[i]){
 							node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,iomodel,analysis,FSvelocityEnum);
 							nodes->AddObject(node);
@@ -602,6 +600,4 @@
 
 	/*Clean up*/
-	xDelete<bool>(my_faces);
-	xDelete<bool>(my_edges);
 	xDelete<bool>(my_nodes);
 }
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/EdgesPartitioning.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/EdgesPartitioning.cpp	(revision 23166)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/EdgesPartitioning.cpp	(revision 23167)
@@ -8,8 +8,8 @@
 #include "./ModelProcessorx.h"
 
-void EdgesPartitioning(bool** pmy_edges,IoModel* iomodel){
+void EdgesPartitioning(IoModel* iomodel){
 
-	/*Intermediaries*/
-	int elementnbe;
+	/*If faces are already present, exit*/
+	if(iomodel->my_edges) return;
 
 	/*Get edges and elements*/
@@ -18,4 +18,5 @@
 
 	/*Mesh dependent variables*/
+	int elementnbe;
 	switch(iomodel->meshelementtype){
 		case TriaEnum:  elementnbe = 3; break;
@@ -26,15 +27,12 @@
 
 	/*output: */
-	bool* my_edges=xNewZeroInit<bool>(iomodel->numberofedges);
+	iomodel->my_edges=xNewZeroInit<bool>(iomodel->numberofedges);
 
 	for(int i=0;i<iomodel->numberofelements;i++){
 		if(iomodel->my_elements[i]){
 			for(int j=0;j<elementnbe;j++){
-				my_edges[iomodel->elementtoedgeconnectivity[i*elementnbe+j]] = true;
+				iomodel->my_edges[iomodel->elementtoedgeconnectivity[i*elementnbe+j]] = true;
 			}
 		}
 	}
-
-	/*Free data and assign output pointers */
-	*pmy_edges=my_edges;
 }
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/FacesPartitioning.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/FacesPartitioning.cpp	(revision 23166)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/FacesPartitioning.cpp	(revision 23167)
@@ -8,8 +8,8 @@
 #include "./ModelProcessorx.h"
 
-void FacesPartitioning(bool** pmy_faces,IoModel* iomodel){
+void FacesPartitioning(IoModel* iomodel){
 
-	/*Intermediaries*/
-	int elementnbf;
+	/*If faces are already present, exit*/
+	if(iomodel->my_faces) return;
 
 	/*Get faces and elements*/
@@ -18,4 +18,5 @@
 
 	/*Mesh dependent variables*/
+	int elementnbf;
 	if(iomodel->domaintype==Domain2DhorizontalEnum){
 		elementnbf = 3;
@@ -31,5 +32,5 @@
 	}
 	/*output: */
-	bool* my_faces=xNewZeroInit<bool>(iomodel->numberoffaces);
+	iomodel->my_faces=xNewZeroInit<bool>(iomodel->numberoffaces);
 
 	for(int i=0;i<iomodel->numberofelements;i++){
@@ -38,10 +39,7 @@
 				_assert_(iomodel->elementtofaceconnectivity[i*elementnbf+j] >= 0);
 				_assert_(iomodel->elementtofaceconnectivity[i*elementnbf+j] <  iomodel->numberoffaces);
-				my_faces[iomodel->elementtofaceconnectivity[i*elementnbf+j]] = true;
+				iomodel->my_faces[iomodel->elementtofaceconnectivity[i*elementnbf+j]] = true;
 			}
 		}
 	}
-
-	/*Free data and assign output pointers */
-	*pmy_faces=my_faces;
 }
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h	(revision 23166)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h	(revision 23167)
@@ -26,6 +26,6 @@
 void ElementsAndVerticesPartitioning(bool** pmy_elements, int** pmy_vertices, IoModel* iomodel);
 void NodesPartitioning(bool** pmy_nodes,bool* my_elements, int* my_vertices,  IoModel* iomodel, bool continuous);
-void FacesPartitioning(bool** pmy_faces,IoModel* iomodel);
-void EdgesPartitioning(bool** pmy_edges,IoModel* iomodel);
+void FacesPartitioning(IoModel* iomodel);
+void EdgesPartitioning(IoModel* iomodel);
 
 /*Mesh properties*/
Index: /issm/trunk-jpl/src/c/solutionsequences/convergence.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/convergence.cpp	(revision 23166)
+++ /issm/trunk-jpl/src/c/solutionsequences/convergence.cpp	(revision 23167)
@@ -37,15 +37,18 @@
 	if (VerboseConvergence()){
 
-		//compute KUF = KU - F = K*U - F
+		/*compute KUF = KU - F = K*U - F*/
 		KU=uf->Duplicate(); Kff->MatMult(uf,KU);
 		KUF=KU->Duplicate(); KU->Copy(KUF); KUF->AYPX(pf,-1.0);
 
-		//compute norm(KUF), norm(F) and residue
+		/*compute norm(KUF), norm(F) and residue*/
 		nKUF=KUF->Norm(NORM_TWO);
 		nF=pf->Norm(NORM_TWO);
 		solver_residue=nKUF/nF;
 		_printf0_("\n" << "   solver residue: norm(KU-F)/norm(F)=" << solver_residue << "\n");
+		if(xIsNan<IssmDouble>(solver_residue)){
+			//Kff->Echo();
+		}
 
-		//clean up
+		/*clean up*/
 		delete KU;
 		delete KUF;
@@ -54,5 +57,5 @@
 	/*Force equilibrium (Mandatory)*/
 
-	//compute K[n]U[n-1]F = K[n]U[n-1] - F
+	/*compute K[n]U[n-1]F = K[n]U[n-1] - F*/
 	_assert_(uf); _assert_(Kff);
 	KUold=uf->Duplicate(); Kff->MatMult(old_uf,KUold);