Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 23567)
+++ /issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 23568)
@@ -137,4 +137,6 @@
 	this->my_faces=NULL;
 	this->my_edges=NULL;
+	this->my_vedges=NULL;
+	this->my_hedges=NULL;
 	this->my_vertices=NULL;
 	this->epart=NULL;
@@ -147,12 +149,18 @@
 	this->numberoffaces=-1;
 	this->numberofedges=-1;
+	this->numberofverticaledges=-1;	
+	this->numberofhorizontaledges=-1;	
 	this->facescols=-1;
 	this->elements=NULL;
 	this->faces=NULL;
 	this->edges=NULL;
-	this->elementtofaceconnectivity      =NULL;
-	this->elementtoedgeconnectivity      =NULL;
-	this->singlenodetoelementconnectivity=NULL;
-	this->numbernodetoelementconnectivity=NULL;
+	this->verticaledges=NULL;
+	this->horizontaledges=NULL;
+	this->elementtofaceconnectivity           = NULL;
+	this->elementtoedgeconnectivity           = NULL;
+	this->elementtoverticaledgeconnectivity   = NULL;
+	this->elementtohorizontaledgeconnectivity = NULL;
+	this->singlenodetoelementconnectivity     = NULL;
+	this->numbernodetoelementconnectivity     = NULL;
 }/*}}}*/
 IoModel::IoModel(FILE* iomodel_handle,int solution_enum_in,bool trace,IssmPDouble* X){/*{{{*/
@@ -197,4 +205,6 @@
 	this->my_faces = NULL;
 	this->my_edges = NULL;
+	this->my_vedges = NULL;
+	this->my_hedges = NULL;
 	this->my_vertices = NULL;
 	this->epart = NULL;
@@ -207,11 +217,15 @@
 	FetchData(&this->numberofelements,"md.mesh.numberofelements");
 	FetchData(&this->elements,NULL,NULL,"md.mesh.elements");
-	this->facescols                       = -1;
-	this->faces                           = NULL;
-	this->edges                           = NULL;
-	this->elementtofaceconnectivity       = NULL;
-	this->elementtoedgeconnectivity       = NULL;
-	this->singlenodetoelementconnectivity = NULL;
-	this->numbernodetoelementconnectivity = NULL;
+	this->facescols                           = -1;
+	this->faces                               = NULL;
+	this->edges                               = NULL;
+	this->verticaledges                       = NULL;
+	this->horizontaledges                     = NULL;
+	this->elementtofaceconnectivity           = NULL;
+	this->elementtoedgeconnectivity           = NULL;
+	this->elementtoverticaledgeconnectivity   = NULL;
+	this->elementtohorizontaledgeconnectivity = NULL;
+	this->singlenodetoelementconnectivity     = NULL;
+	this->numbernodetoelementconnectivity     = NULL;
 }/*}}}*/
 IoModel::~IoModel(){/*{{{*/
@@ -239,4 +253,6 @@
 	xDelete<bool>(this->my_faces);
 	xDelete<bool>(this->my_edges);
+	xDelete<bool>(this->my_vedges);
+	xDelete<bool>(this->my_hedges);
 	xDelete<bool>(this->my_vertices);
 	xDelete<int>(this->epart);
@@ -245,6 +261,10 @@
 	xDelete<int>(this->faces);
 	xDelete<int>(this->edges);
+	xDelete<int>(this->verticaledges);
+	xDelete<int>(this->horizontaledges);
 	xDelete<int>(this->elementtofaceconnectivity);
 	xDelete<int>(this->elementtoedgeconnectivity);
+	xDelete<int>(this->elementtoverticaledgeconnectivity);
+	xDelete<int>(this->elementtohorizontaledgeconnectivity);
 	xDelete<int>(this->singlenodetoelementconnectivity);
 	xDelete<int>(this->numbernodetoelementconnectivity);
Index: /issm/trunk-jpl/src/c/classes/IoModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/IoModel.h	(revision 23567)
+++ /issm/trunk-jpl/src/c/classes/IoModel.h	(revision 23568)
@@ -66,4 +66,6 @@
 		bool *my_faces;
 		bool *my_edges;
+		bool *my_vedges;
+		bool *my_hedges;
 		bool *my_vertices;
 		int  *epart;
@@ -74,5 +76,9 @@
 		int *elements;
 		int *edges;
+		int *verticaledges;
+		int *horizontaledges;
 		int *elementtoedgeconnectivity;
+		int *elementtoverticaledgeconnectivity;
+		int *elementtohorizontaledgeconnectivity;
 		int *elementtofaceconnectivity;
 		int *faces;
@@ -81,4 +87,6 @@
 		int *numbernodetoelementconnectivity;
 		int  numberofedges;
+		int  numberofverticaledges;
+		int  numberofhorizontaledges;
 		int  numberofelements;
 		int  numberoffaces;
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp	(revision 23567)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp	(revision 23568)
@@ -17,7 +17,6 @@
 
 	/*Intermediaries*/
-	bool exist;
 	int  i,j,v1,v2,v3;
-	int  maxnbe,nbe,elementnbe,elementnbv;
+	int  elementnbe,elementnbv;
 	int *elementedges         = NULL;
 	int *elementedges_markers = NULL;
@@ -30,7 +29,7 @@
 			elementedges         = xNew<int>(elementnbe*2);
 			elementedges_markers = xNew<int>(elementnbe);
-			elementedges[2*0+0] = 1; elementedges[2*0+1] = 2; elementedges_markers[0] = 1;
-			elementedges[2*1+0] = 2; elementedges[2*1+1] = 0; elementedges_markers[1] = 1;
-			elementedges[2*2+0] = 0; elementedges[2*2+1] = 1; elementedges_markers[2] = 1;
+			elementedges[2*0+0] = 1; elementedges[2*0+1] = 2; elementedges_markers[0] = -1;
+			elementedges[2*1+0] = 2; elementedges[2*1+1] = 0; elementedges_markers[1] = -1;
+			elementedges[2*2+0] = 0; elementedges[2*2+1] = 1; elementedges_markers[2] = -1;
 			break;
 		case TetraEnum:
@@ -39,10 +38,10 @@
 			elementedges         = xNew<int>(elementnbe*2);
 			elementedges_markers = xNew<int>(elementnbe);
-			elementedges[2*0+0] = 1; elementedges[2*0+1] = 2; elementedges_markers[0] = 1;
-			elementedges[2*1+0] = 0; elementedges[2*1+1] = 2; elementedges_markers[1] = 1;
-			elementedges[2*2+0] = 0; elementedges[2*2+1] = 1; elementedges_markers[2] = 1;
-			elementedges[2*3+0] = 1; elementedges[2*3+1] = 3; elementedges_markers[3] = 1;
-			elementedges[2*4+0] = 2; elementedges[2*4+1] = 3; elementedges_markers[4] = 1;
-			elementedges[2*5+0] = 0; elementedges[2*5+1] = 3; elementedges_markers[5] = 1;
+			elementedges[2*0+0] = 1; elementedges[2*0+1] = 2; elementedges_markers[0] = -1;
+			elementedges[2*1+0] = 0; elementedges[2*1+1] = 2; elementedges_markers[1] = -1;
+			elementedges[2*2+0] = 0; elementedges[2*2+1] = 1; elementedges_markers[2] = -1;
+			elementedges[2*3+0] = 1; elementedges[2*3+1] = 3; elementedges_markers[3] = -1;
+			elementedges[2*4+0] = 2; elementedges[2*4+1] = 3; elementedges_markers[4] = -1;
+			elementedges[2*5+0] = 0; elementedges[2*5+1] = 3; elementedges_markers[5] = -1;
 			break;
 		case PentaEnum:
@@ -66,9 +65,17 @@
 
 	/*Maximum number of edges*/
-	maxnbe = elementnbe*iomodel->numberofelements;
+	int maxnbe = elementnbe*iomodel->numberofelements;
 
 	/*Initialize intermediaries*/
-	int *edgestemp                 = xNew<int>(maxnbe*3);                             /*format: [vertex1 vertex2 marker]       */
-	int *element_edge_connectivity = xNew<int>(iomodel->numberofelements*elementnbe); /*format: [edge1 edge2 ... edgen] */
+	int *edgestemp                  = xNew<int>(maxnbe*3);                             /*format: [vertex1 vertex2 marker]*/
+	int *vedgestemp                 = xNew<int>(maxnbe*2);                             /*format: [vertex1 vertex2]       */
+	int *hedgestemp                 = xNew<int>(maxnbe*2);                             /*format: [vertex1 vertex2]       */
+	int *element_edge_connectivity  = xNew<int>(iomodel->numberofelements*elementnbe); /*format: [edge1 edge2 ... edgen] */
+	int *element_vedge_connectivity = NULL;
+	int *element_hedge_connectivity = NULL;
+	if(iomodel->meshelementtype==PentaEnum){
+		element_vedge_connectivity  = xNew<int>(iomodel->numberofelements*3); /*format: [edge1 edge2 ... edgen] */
+		element_hedge_connectivity  = xNew<int>(iomodel->numberofelements*6); /*format: [edge1 edge2 ... edgen] */
+	}
 
 	/*Initialize chain*/
@@ -78,5 +85,5 @@
 
 	/*Initialize number of edges*/
-	nbe = 0;
+	int nbe  = 0;
 
 	for(i=0;i<iomodel->numberofelements;i++){
@@ -93,5 +100,5 @@
 
 			/*This edge a priori has not been processed yet*/
-			exist = false;
+			bool exist = false;
 
 			/*Go through all processed edges connected to v1 and check whether we have seen this edge yet*/
@@ -125,22 +132,83 @@
 		}
 	}
+	int nbve = 0;
+	int nbhe = 0;
+	/*vertical/horizontal edges*/
+	if(iomodel->meshelementtype==PentaEnum){
+		for(i=0;i<iomodel->numberofvertices;i++) head_minv[i]=-1;
+		for(i=0;i<iomodel->numberofelements;i++){
+			for(j=0;j<3;j++){
+				v1 = iomodel->elements[i*elementnbv+elementedges[2*j+0]]-1; _assert_(v1>=0 & v1<iomodel->numberofvertices);
+				v2 = iomodel->elements[i*elementnbv+elementedges[2*j+1]]-1; _assert_(v2>=0 & v2<iomodel->numberofvertices);
+				if(v2<v1){ v3=v2; v2=v1; v1=v3;}
+				bool exist = false;
+				for(int e=head_minv[v1]; e!=-1; e=next_edge[e]){
+					if(vedgestemp[e*2+1]==v2+1){
+						exist = true;
+						element_vedge_connectivity[i*3+j]=e;
+						break;
+					}
+				}
+				if(!exist){
+					vedgestemp[nbve*2+0] = v1+1;
+					vedgestemp[nbve*2+1] = v2+1;
+					element_vedge_connectivity[i*3+j]=nbve;
+					head_minv[v1] = nbve;
+					nbve++;
+				}
+			}
+		}
+		for(i=0;i<iomodel->numberofvertices;i++) head_minv[i]=-1;
+		for(i=0;i<iomodel->numberofelements;i++){
+			for(j=3;j<9;j++){
+				v1 = iomodel->elements[i*elementnbv+elementedges[2*j+0]]-1; _assert_(v1>=0 & v1<iomodel->numberofvertices);
+				v2 = iomodel->elements[i*elementnbv+elementedges[2*j+1]]-1; _assert_(v2>=0 & v2<iomodel->numberofvertices);
+				if(v2<v1){ v3=v2; v2=v1; v1=v3;}
+				bool exist = false;
+				for(int e=head_minv[v1]; e!=-1; e=next_edge[e]){
+					if(hedgestemp[e*2+1]==v2+1){
+						exist = true;
+						element_hedge_connectivity[i*6+(j-3)]=e;
+						break;
+					}
+				}
+				if(!exist){
+					hedgestemp[nbhe*2+0] = v1+1;
+					hedgestemp[nbhe*2+1] = v2+1;
+					element_hedge_connectivity[i*6+(j-3)]=nbhe;
+					head_minv[v1] = nbhe;
+					nbhe++;
+				}
+			}
+		}
+	}
 
 	/*Clean up*/
 	xDelete<int>(head_minv);
 	xDelete<int>(next_edge);
+	xDelete<int>(elementedges);
 	xDelete<int>(elementedges_markers);
 
 	/*Create final edges*/
-	int* edges = xNew<int>(nbe*3); /*format: [vertex1 vertex2 marker]*/
+	int* edges = xNew<int>(nbe*3);
 	for(int i=0;i<3*nbe;i++) edges[i] = edgestemp[i];
-
-	/*Clean up*/
 	xDelete<int>(edgestemp);
-	xDelete<int>(elementedges);
+	int* vedges = xNew<int>(nbve*2);
+	for(int i=0;i<2*nbve;i++) vedges[i] = vedgestemp[i];
+	xDelete<int>(vedgestemp);
+	int* hedges = xNew<int>(nbhe*2);
+	for(int i=0;i<2*nbhe;i++) hedges[i] = hedgestemp[i];
+	xDelete<int>(hedgestemp);
 
 	/*Assign output pointers*/
-	iomodel->edges                     = edges;
+	iomodel->edges           = edges;
+	iomodel->verticaledges   = vedges;
+	iomodel->horizontaledges = hedges;
 	iomodel->elementtoedgeconnectivity = element_edge_connectivity;
+	iomodel->elementtoverticaledgeconnectivity = element_vedge_connectivity;
+	iomodel->elementtohorizontaledgeconnectivity = element_hedge_connectivity;
 	iomodel->numberofedges             = nbe;
+	iomodel->numberofverticaledges     = nbve;
+	iomodel->numberofhorizontaledges   = nbhe;
 }/*}}}*/
 void EdgeOnBoundaryFlags(bool** pflags,IoModel* iomodel){/*{{{*/
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/EdgesPartitioning.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/EdgesPartitioning.cpp	(revision 23567)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/EdgesPartitioning.cpp	(revision 23568)
@@ -27,5 +27,5 @@
 
 	/*output: */
-	iomodel->my_edges=xNewZeroInit<bool>(iomodel->numberofedges);
+	iomodel->my_edges  = xNewZeroInit<bool>(iomodel->numberofedges);
 
 	for(int i=0;i<iomodel->numberofelements;i++){
@@ -36,3 +36,14 @@
 		}
 	}
+
+	if(iomodel->meshelementtype==PentaEnum){
+		iomodel->my_vedges = xNewZeroInit<bool>(iomodel->numberofverticaledges);
+		iomodel->my_hedges = xNewZeroInit<bool>(iomodel->numberofhorizontaledges);
+		for(int i=0;i<iomodel->numberofelements;i++){
+			if(iomodel->my_elements[i]){
+				for(int j=0;j<3;j++) iomodel->my_vedges[iomodel->elementtoverticaledgeconnectivity[i*3+j]]   = true;
+				for(int j=0;j<6;j++) iomodel->my_hedges[iomodel->elementtohorizontaledgeconnectivity[i*6+j]] = true;
+			}
+		}
+	}
 }