Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17134)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17135)
@@ -2834,34 +2834,34 @@
 			numnodes         = 30;
 			penta_node_ids   = xNew<int>(numnodes);
-			penta_node_ids[ 0]=iomodel->nodecounter+iomodel->elements[6*index+0];
-			penta_node_ids[ 1]=iomodel->nodecounter+iomodel->elements[6*index+1];
-			penta_node_ids[ 2]=iomodel->nodecounter+iomodel->elements[6*index+2];
-			penta_node_ids[ 3]=iomodel->nodecounter+iomodel->elements[6*index+3];
-			penta_node_ids[ 4]=iomodel->nodecounter+iomodel->elements[6*index+4];
-			penta_node_ids[ 5]=iomodel->nodecounter+iomodel->elements[6*index+5];
-			penta_node_ids[ 6]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+0]+1;
-			penta_node_ids[ 7]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+1]+1;
-			penta_node_ids[ 8]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+2]+1;
-			penta_node_ids[ 9]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+3]+1;
-			penta_node_ids[10]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+4]+1;
-			penta_node_ids[11]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+5]+1;
-			penta_node_ids[12]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+6]+1;
-			penta_node_ids[13]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+7]+1;
-			penta_node_ids[14]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+8]+1;
-			penta_node_ids[15]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+0]+2;
-			penta_node_ids[16]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+1]+2;
-			penta_node_ids[17]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+2]+2;
-			penta_node_ids[18]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+0]+3;
-			penta_node_ids[19]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+1]+3;
-			penta_node_ids[20]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+2]+3;
-			penta_node_ids[21]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+3]+2;
-			penta_node_ids[22]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+4]+2;
-			penta_node_ids[23]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+5]+2;
-			penta_node_ids[24]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+3]+3;
-			penta_node_ids[25]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+4]+3;
-			penta_node_ids[26]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+5]+3;
-			penta_node_ids[27]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+3]+4;
-			penta_node_ids[28]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+4]+4;
-			penta_node_ids[29]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+5]+4;
+			penta_node_ids[ 0]=iomodel->nodecounter+iomodel->elements[6*index+0]; /*Vertex 1*/
+			penta_node_ids[ 1]=iomodel->nodecounter+iomodel->elements[6*index+1]; /*Vertex 2*/
+			penta_node_ids[ 2]=iomodel->nodecounter+iomodel->elements[6*index+2]; /*Vertex 3*/
+			penta_node_ids[ 3]=iomodel->nodecounter+iomodel->elements[6*index+3]; /*Vertex 4*/
+			penta_node_ids[ 4]=iomodel->nodecounter+iomodel->elements[6*index+4]; /*Vertex 5*/
+			penta_node_ids[ 5]=iomodel->nodecounter+iomodel->elements[6*index+5]; /*Vertex 6*/
+			penta_node_ids[ 6]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+1; /*mid vertical edge 1*/
+			penta_node_ids[ 7]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+1; /*mid vertical edge 2*/
+			penta_node_ids[ 8]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+1; /*mid vertical edge 3*/
+			penta_node_ids[ 9]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+3]+1; /*mid basal edge 1*/
+			penta_node_ids[10]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+4]+1; /*mid basal edge 2*/
+			penta_node_ids[11]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+5]+1; /*mid basal edge 3*/
+			penta_node_ids[12]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+6]+1; /*mid top edge 1*/
+			penta_node_ids[13]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+7]+1; /*mid top edge 2*/
+			penta_node_ids[14]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+8]+1; /*mid top edge 3*/
+			penta_node_ids[15]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+2; /* 1/4 vertical edge 1*/
+			penta_node_ids[16]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+2; /* 1/4 vertical edge 2*/
+			penta_node_ids[17]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+2; /* 1/4 vertical edge 3*/
+			penta_node_ids[18]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+3; /* 2/4 vertical edge 1*/
+			penta_node_ids[19]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+3; /* 2/4 vertical edge 2*/
+			penta_node_ids[20]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+3; /* 2/4 vertical edge 3*/
+			penta_node_ids[21]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+2]+1; /* 1/4 vertical face 1*/
+			penta_node_ids[22]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+3]+1; /* 1/4 vertical face 2*/
+			penta_node_ids[23]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+4]+1; /* 1/4 vertical face 3*/
+			penta_node_ids[24]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+2]+2; /* 2/4 vertical face 1*/
+			penta_node_ids[25]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+3]+2; /* 2/4 vertical face 2*/
+			penta_node_ids[26]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+4]+2; /* 2/4 vertical face 3*/
+			penta_node_ids[27]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+2]+3; /* 3/4 vertical face 1*/
+			penta_node_ids[28]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+3]+3; /* 3/4 vertical face 2*/
+			penta_node_ids[29]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+4]+3; /* 3/4 vertical face 3*/
 			break;
 		case P2Enum:
Index: /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp	(revision 17134)
+++ /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp	(revision 17135)
@@ -39,7 +39,8 @@
 	bool        spcpresent       = false;
 
-	/*P2 finite elements*/
-	int   v1,v2;
+	/*Higher-order finite elements*/
+	int   v1,v2,v3,v4;
 	bool *my_edges = NULL;
+	bool *my_faces = NULL;
 
 	switch(finite_element){
@@ -72,4 +73,5 @@
 		case P2xP4Enum:
 			EdgesPartitioning(&my_edges,iomodel);
+			FacesPartitioning(&my_faces,iomodel);
 			break;
 		default:
@@ -126,9 +128,9 @@
 							v2 = iomodel->edges[3*i+1]-1;
 							if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
-								constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+4*i+1,
-												dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
-								constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+2,iomodel->nodecounter+iomodel->numberofvertices+4*i+2,
-												dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
-								constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+3,iomodel->nodecounter+iomodel->numberofvertices+4*i+3,
+								constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+3*i+1,
+												dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
+								constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+2,iomodel->nodecounter+iomodel->numberofvertices+3*i+2,
+												dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
+								constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+3,iomodel->nodecounter+iomodel->numberofvertices+3*i+3,
 												dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
 								count=count+3;
@@ -136,29 +138,29 @@
 						}
 					}
-					if(iomodel->edges[i*3+2]==1){/*Basal edges*/
-						if(my_edges[i]){
-							v1 = iomodel->edges[3*i+0]-1;
-							v2 = iomodel->edges[3*i+1]-1;
-							if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
-								constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+4*i+1,
-												dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
-								constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+2,iomodel->nodecounter+iomodel->numberofvertices+4*i+2,
-												dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
-								constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+3,iomodel->nodecounter+iomodel->numberofvertices+4*i+3,
-												dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
-								constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+4,iomodel->nodecounter+iomodel->numberofvertices+4*i+4,
-												dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
-								count=count+4;
-							}
-						}
-					}
-					if(iomodel->edges[i*3+2]==3){/*Top edges*/
-						if(my_edges[i]){
-							v1 = iomodel->edges[3*i+0]-1;
-							v2 = iomodel->edges[3*i+1]-1;
-							if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
-								constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+4*i+1,
+					if(iomodel->edges[i*3+2]==1){/*Horizontal edges*/
+						if(my_edges[i]){
+							v1 = iomodel->edges[3*i+0]-1;
+							v2 = iomodel->edges[3*i+1]-1;
+							if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
+								constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+3*i+1,
 												dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
 								count=count+1;
+							}
+						}
+					}
+				}
+				for(i=0;i<iomodel->numberoffaces;i++){
+					if(iomodel->faces[i*6+5]==2){/*Vertical quads*/
+						if(my_faces[i]){
+							v1 = iomodel->faces[6*i+0]-1;
+							v2 = iomodel->faces[6*i+1]-1;
+							if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
+								constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*i+1,
+												dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
+								constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+2,iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*i+2,
+												dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
+								constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+3,iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*i+3,
+												dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
+								count=count+3;
 							}
 						}
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp	(revision 17134)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp	(revision 17135)
@@ -52,7 +52,7 @@
 		elementedges[2*4+0] = 2; elementedges[2*4+1] = 0; elementedges_markers[4] = 1;
 		elementedges[2*5+0] = 0; elementedges[2*5+1] = 1; elementedges_markers[5] = 1;
-		elementedges[2*6+0] = 4; elementedges[2*6+1] = 5; elementedges_markers[6] = 3;
-		elementedges[2*7+0] = 5; elementedges[2*7+1] = 3; elementedges_markers[7] = 3;
-		elementedges[2*8+0] = 3; elementedges[2*8+1] = 4; elementedges_markers[8] = 3;
+		elementedges[2*6+0] = 4; elementedges[2*6+1] = 5; elementedges_markers[6] = 1;
+		elementedges[2*7+0] = 5; elementedges[2*7+1] = 3; elementedges_markers[7] = 1;
+		elementedges[2*8+0] = 3; elementedges[2*8+1] = 4; elementedges_markers[8] = 1;
 	}
 	else{
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateFaces.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateFaces.cpp	(revision 17134)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateFaces.cpp	(revision 17135)
@@ -160,7 +160,7 @@
 		elementfaces_markers[3] = 2;
 		elementfaces_markers[4] = 2;
-		elementfaces[cols*2+0] = 4; elementfaces[cols*2+1] = 0;  elementfaces[cols*2+2] = 1; elementfaces[cols*2+3] = 4;  elementfaces[cols*2+4] = 3;
-		elementfaces[cols*3+0] = 4; elementfaces[cols*3+1] = 1;  elementfaces[cols*3+2] = 2; elementfaces[cols*3+3] = 5;  elementfaces[cols*3+4] = 4;
-		elementfaces[cols*4+0] = 4; elementfaces[cols*4+1] = 2;  elementfaces[cols*4+2] = 0; elementfaces[cols*4+3] = 3;  elementfaces[cols*4+4] = 5;
+		elementfaces[cols*2+0] = 4; elementfaces[cols*2+1] = 1;  elementfaces[cols*2+2] = 2; elementfaces[cols*2+3] = 5;  elementfaces[cols*2+4] = 4;
+		elementfaces[cols*3+0] = 4; elementfaces[cols*3+1] = 2;  elementfaces[cols*3+2] = 0; elementfaces[cols*3+3] = 3;  elementfaces[cols*3+4] = 5;
+		elementfaces[cols*4+0] = 4; elementfaces[cols*4+1] = 0;  elementfaces[cols*4+2] = 1; elementfaces[cols*4+3] = 4;  elementfaces[cols*4+4] = 3;
 	}
 	else{
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp	(revision 17134)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp	(revision 17135)
@@ -14,4 +14,5 @@
 	int   i,j,counter,vnodes,lid=0;
 	int   id0 = iomodel->nodecounter;
+	bool *my_faces = NULL;
 	bool *my_edges = NULL;
 	bool *my_nodes = NULL;
@@ -144,4 +145,5 @@
 		case P2xP4Enum:
 			EdgesPartitioning(&my_edges,iomodel);
+			FacesPartitioning(&my_faces,iomodel);
 			for(i=0;i<iomodel->numberofvertices;i++){
 				if(iomodel->my_vertices[i]){
@@ -153,32 +155,39 @@
 				if(iomodel->edges[i*3+2]==2){/*Vertical edges*/
 					if(my_edges[i]){
-						node = new Node(id0+iomodel->numberofvertices+4*i+1,counter+1,lid++,0,iomodel,analysis,approximation);
-						nodes->AddObject(node);
-						node = new Node(id0+iomodel->numberofvertices+4*i+2,counter+2,lid++,0,iomodel,analysis,approximation);
-						nodes->AddObject(node);
-						node = new Node(id0+iomodel->numberofvertices+4*i+3,counter+3,lid++,0,iomodel,analysis,approximation);
+						node = new Node(id0+iomodel->numberofvertices+3*i+1,counter+1,lid++,0,iomodel,analysis,approximation);
+						nodes->AddObject(node);
+						node = new Node(id0+iomodel->numberofvertices+3*i+2,counter+2,lid++,0,iomodel,analysis,approximation);
+						nodes->AddObject(node);
+						node = new Node(id0+iomodel->numberofvertices+3*i+3,counter+3,lid++,0,iomodel,analysis,approximation);
 						nodes->AddObject(node);
 					}
 					counter=counter+3;
 				}
-				else if(iomodel->edges[i*3+2]==1){/*Basal edges*/
-					if(my_edges[i]){
-						node = new Node(id0+iomodel->numberofvertices+4*i+1,counter+1,lid++,0,iomodel,analysis,approximation);
-						nodes->AddObject(node);
-						node = new Node(id0+iomodel->numberofvertices+4*i+2,counter+2,lid++,0,iomodel,analysis,approximation);
-						nodes->AddObject(node);
-						node = new Node(id0+iomodel->numberofvertices+4*i+3,counter+3,lid++,0,iomodel,analysis,approximation);
-						nodes->AddObject(node);
-						node = new Node(id0+iomodel->numberofvertices+4*i+4,counter+4,lid++,0,iomodel,analysis,approximation);
-						nodes->AddObject(node);
-					}
-					counter=counter+4;
-				}
-				else if(iomodel->edges[i*3+2]==3){/*Top edges*/
-					if(my_edges[i]){
-						node = new Node(id0+iomodel->numberofvertices+4*i+1,counter+1,lid++,0,iomodel,analysis,approximation);
+				else if(iomodel->edges[i*3+2]==1){/*Horizontal edges*/
+					if(my_edges[i]){
+						node = new Node(id0+iomodel->numberofvertices+3*i+1,counter+1,lid++,0,iomodel,analysis,approximation);
 						nodes->AddObject(node);
 					}
 					counter=counter+1;
+				}
+				else{
+					_error_("not supported");
+				}
+			}
+			id0 = id0+iomodel->numberofvertices+3*iomodel->numberofedges;
+			for(i=0;i<iomodel->numberoffaces;i++){
+				if(iomodel->faces[i*6+5]==2){/*Vertical quads*/
+					if(my_faces[i]){
+						node = new Node(id0+3*i+1,counter+1,lid++,0,iomodel,analysis,approximation);
+						nodes->AddObject(node);
+						node = new Node(id0+3*i+2,counter+2,lid++,0,iomodel,analysis,approximation);
+						nodes->AddObject(node);
+						node = new Node(id0+3*i+3,counter+3,lid++,0,iomodel,analysis,approximation);
+						nodes->AddObject(node);
+					}
+					counter=counter+3;
+				}
+				else if(iomodel->faces[i*6+5]==1){/*Triangular base/top*/
+					/*Nothing*/
 				}
 				else{
@@ -291,4 +300,5 @@
 			_assert_(approximation==FSApproximationEnum);
 			EdgesPartitioning(&my_edges,iomodel);
+			FacesPartitioning(&my_faces,iomodel);
 			/*P2xP4 velocity*/
 			for(i=0;i<iomodel->numberofvertices;i++){
@@ -360,4 +370,5 @@
 
 	/*Clean up*/
+	xDelete<bool>(my_faces);
 	xDelete<bool>(my_edges);
 	xDelete<bool>(my_nodes);