Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 23569)
+++ /issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 23570)
@@ -136,4 +136,5 @@
 	this->my_elements=NULL;
 	this->my_faces=NULL;
+	this->my_vfaces=NULL;
 	this->my_edges=NULL;
 	this->my_vedges=NULL;
@@ -148,4 +149,5 @@
 	this->numberofelements=-1;
 	this->numberoffaces=-1;
+	this->numberofverticalfaces=-1;
 	this->numberofedges=-1;
 	this->numberofverticaledges=-1;	
@@ -154,8 +156,10 @@
 	this->elements=NULL;
 	this->faces=NULL;
+	this->verticalfaces=NULL;
 	this->edges=NULL;
 	this->verticaledges=NULL;
 	this->horizontaledges=NULL;
 	this->elementtofaceconnectivity           = NULL;
+	this->elementtoverticalfaceconnectivity   = NULL;
 	this->elementtoedgeconnectivity           = NULL;
 	this->elementtoverticaledgeconnectivity   = NULL;
@@ -203,10 +207,11 @@
 	/*Initialize permanent data: */
 	this->my_elements = NULL;
-	this->my_faces = NULL;
-	this->my_edges = NULL;
-	this->my_vedges = NULL;
-	this->my_hedges = NULL;
+	this->my_faces    = NULL;
+	this->my_vfaces   = NULL;
+	this->my_edges    = NULL;
+	this->my_vedges   = NULL;
+	this->my_hedges   = NULL;
 	this->my_vertices = NULL;
-	this->epart = NULL;
+	this->epart       = NULL;
 
 	FindConstant(&this->domaintype,"md.mesh.domain_type");
@@ -219,8 +224,10 @@
 	this->facescols                           = -1;
 	this->faces                               = NULL;
+	this->verticalfaces                       = NULL;
 	this->edges                               = NULL;
 	this->verticaledges                       = NULL;
 	this->horizontaledges                     = NULL;
 	this->elementtofaceconnectivity           = NULL;
+	this->elementtoverticalfaceconnectivity   = NULL;
 	this->elementtoedgeconnectivity           = NULL;
 	this->elementtoverticaledgeconnectivity   = NULL;
@@ -252,4 +259,5 @@
 	xDelete<bool>(this->my_elements);
 	xDelete<bool>(this->my_faces);
+	xDelete<bool>(this->my_vfaces);
 	xDelete<bool>(this->my_edges);
 	xDelete<bool>(this->my_vedges);
@@ -260,8 +268,10 @@
 	xDelete<int>(this->elements);
 	xDelete<int>(this->faces);
+	xDelete<int>(this->verticalfaces);
 	xDelete<int>(this->edges);
 	xDelete<int>(this->verticaledges);
 	xDelete<int>(this->horizontaledges);
 	xDelete<int>(this->elementtofaceconnectivity);
+	xDelete<int>(this->elementtoverticalfaceconnectivity);
 	xDelete<int>(this->elementtoedgeconnectivity);
 	xDelete<int>(this->elementtoverticaledgeconnectivity);
Index: /issm/trunk-jpl/src/c/classes/IoModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/IoModel.h	(revision 23569)
+++ /issm/trunk-jpl/src/c/classes/IoModel.h	(revision 23570)
@@ -65,4 +65,5 @@
 		bool *my_elements;
 		bool *my_faces;
+		bool *my_vfaces;
 		bool *my_edges;
 		bool *my_vedges;
@@ -82,5 +83,7 @@
 		int *elementtohorizontaledgeconnectivity;
 		int *elementtofaceconnectivity;
+		int *elementtoverticalfaceconnectivity;
 		int *faces;
+		int *verticalfaces;
 		int  facescols;
 		int  meshelementtype;
@@ -91,4 +94,5 @@
 		int  numberofelements;
 		int  numberoffaces;
+		int  numberofverticalfaces;
 		int  numberofvertices;
 		int *singlenodetoelementconnectivity;
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp	(revision 23569)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp	(revision 23570)
@@ -132,7 +132,8 @@
 		}
 	}
+
+	/*vertical/horizontal edges*/
 	int nbve = 0;
 	int nbhe = 0;
-	/*vertical/horizontal edges*/
 	if(iomodel->meshelementtype==PentaEnum){
 		for(i=0;i<iomodel->numberofvertices;i++) head_minv[i]=-1;
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateFaces.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateFaces.cpp	(revision 23569)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateFaces.cpp	(revision 23570)
@@ -125,7 +125,6 @@
 
 	/*Intermediaries*/
-	bool exist;
-	int  i,j,k,v0,cols,facescols;
-	int  maxnbf,nbf,elementnbf,elementnbv,facenbv,facemaxnbv;
+	int  i,j,k,v0,cols,facemaxnbv;
+	int  elementnbf,elementnbv,facenbv;
 	int *elementfaces         = NULL;
 	int *elementfaces_markers = NULL;
@@ -175,12 +174,17 @@
 
 	/*Allocate connectivity*/
-	int *element_face_connectivity = xNew<int>(iomodel->numberofelements*elementnbf); /*format: [face1 face2 ...] */
+	int *element_face_connectivity  = xNew<int>(iomodel->numberofelements*elementnbf); /*format: [face1 face2 ...] */
+	int *element_vface_connectivity = NULL;
+	if(iomodel->meshelementtype==PentaEnum){
+		element_vface_connectivity = xNew<int>(iomodel->numberofelements*3); /*format: [face1 face2 face3] */
+	}
 
 	/*Maximum number of faces for initial allocation*/
-	maxnbf     = elementnbf*iomodel->numberofelements;
-	facescols  = 4+facemaxnbv; _assert_(facescols>6);
+	int maxnbf     = elementnbf*iomodel->numberofelements;
+	int facescols  = 4+facemaxnbv; _assert_(facescols>6);
 
 	/*Initialize intermediaries*/
-	int* facestemp = xNew<int>(maxnbf*facescols);        /*format: [element1 element2 marker nbv vertex1 vertex2 vertex3 ...]    */
+	int* facestemp  = xNew<int>(maxnbf*facescols);        /*format: [element1 element2 marker nbv vertex1 vertex2 vertex3 ...]    */
+	int* vfacestemp = xNew<int>(maxnbf*4);
 	for(i=0;i<maxnbf;i++) facestemp[i*facescols+1]=-1;   /*Initialize second column of faces as -1 (boundary face)               */
 
@@ -191,7 +195,6 @@
 
 	/*Initialize number of faces and list of vertex indices*/
-	nbf = 0;
+	int nbf = 0;
 	int* v = xNew<int>(facemaxnbv);
-
 	for(i=0;i<iomodel->numberofelements;i++){
 		for(j=0;j<elementnbf;j++){
@@ -208,5 +211,5 @@
 
 			/*This face a priori has not been processed yet*/
-			exist = false;
+			bool exist = false;
 
 			/*Go through all processed faces connected to v0 and check whether we have seen this face yet*/
@@ -243,4 +246,32 @@
 	}
 
+	/*Vertical faces*/
+	int nbvf = 0;
+	if(iomodel->meshelementtype==PentaEnum){
+		for(i=0;i<iomodel->numberofvertices;i++) head_minv[i]=-1;
+		for(i=0;i<iomodel->numberofelements;i++){
+			for(j=2;j<5;j++){
+				for(k=0;k<4;k++) v[k] = iomodel->elements[i*elementnbv + elementfaces[cols*j+k+1]] - 1;
+				HeapSort(v,4);
+				v0 = v[0]; _assert_(v0>=0 && v0<iomodel->numberofvertices);
+				bool exist = false;
+				for(int f=head_minv[v0]; f!=-1; f=next_face[f]){
+					if(vfacestemp[f*4+1]==v[1]+1 && vfacestemp[f*4+2]==v[2]+1){
+						exist = true;
+						element_vface_connectivity[i*3+(j-2)]=f;
+						break;
+					}
+				}
+				if(!exist){ _assert_(nbvf<maxnbf);
+					for(k=0;k<4;k++) vfacestemp[nbvf*4+k] = v[k]+1;
+					element_vface_connectivity[i*3+(j-2)]=nbvf;
+					next_face[nbvf] = head_minv[v0];
+					head_minv[v0]   = nbvf;
+					nbvf++;
+				}
+			}
+		}
+	}
+
 	/*Clean up*/
 	xDelete<int>(head_minv);
@@ -254,9 +285,15 @@
 	xMemCpy<int>(faces,facestemp,nbf*facescols);
 	xDelete<int>(facestemp);
+	int* vfaces = xNew<int>(nbvf*4);
+	xMemCpy<int>(vfaces,vfacestemp,nbvf*4);
+	xDelete<int>(vfacestemp);
 
 	/*Assign output pointers*/
-	iomodel->faces                     = faces;
-	iomodel->numberoffaces             = nbf;
-	iomodel->facescols                 = facescols;
-	iomodel->elementtofaceconnectivity = element_face_connectivity;
+	iomodel->faces                             = faces;
+	iomodel->verticalfaces                     = vfaces;
+	iomodel->numberoffaces                     = nbf;
+	iomodel->numberofverticalfaces             = nbvf;
+	iomodel->facescols                         = facescols;
+	iomodel->elementtofaceconnectivity         = element_face_connectivity;
+	iomodel->elementtoverticalfaceconnectivity = element_vface_connectivity;
 }/*}}}*/