Index: /issm/trunk-jpl/src/c/classes/IoModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/IoModel.h	(revision 17121)
+++ /issm/trunk-jpl/src/c/classes/IoModel.h	(revision 17122)
@@ -40,4 +40,5 @@
 		int *edges;
 		int *elementtoedgeconnectivity;
+		int *elementtofaceconnectivity;
 		int *singlenodetoelementconnectivity;
 		int *numbernodetoelementconnectivity;
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateFaces.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateFaces.cpp	(revision 17121)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateFaces.cpp	(revision 17122)
@@ -5,14 +5,26 @@
 #include "../../classes/classes.h"
 #include "../../shared/shared.h"
-
-void CreateFaces(IoModel* iomodel){
+#include "./ModelProcessorx.h"
+
+void CreateFaces(IoModel* iomodel){/*{{{*/
 
 	/*If faces are already present, exit*/
 	if(iomodel->faces) return;
 
-	/*Check Iomodel properties*/
-	if(iomodel->meshtype!=Mesh2DhorizontalEnum)             _error_("only 2d model are supported");
+	/*Some checks*/
 	if(iomodel->numberofvertices<3) _error_("not enough elements in mesh");
 	_assert_(iomodel->elements);
+
+	/*Check Iomodel properties*/
+	if(iomodel->meshtype==Mesh2DhorizontalEnum){
+		/*Keep going*/
+	}
+	else if(iomodel->meshtype==Mesh3DEnum){
+		CreateFaces3d(iomodel);
+		return;
+	}
+	else{
+		_error_("mesh dimension not supported yet");
+	}
 
 	/*Intermediaries*/
@@ -31,5 +43,5 @@
 	/*Initialize chain*/
 	int* head_minv = xNew<int>(iomodel->numberofvertices);
-	int* next_edge = xNew<int>(maxnbf);
+	int* next_face = xNew<int>(maxnbf);
 	for(i=0;i<iomodel->numberofvertices;i++) head_minv[i]=-1;
 
@@ -57,5 +69,5 @@
 			/*Go through all processed faces connected to v1 and check whether we have seen this edge yet*/
 			_assert_(v1>=0 && v1<iomodel->numberofvertices);
-			for(int e=head_minv[v1]; e!=-1; e=next_edge[e]){
+			for(int e=head_minv[v1]; e!=-1; e=next_face[e]){
 				if(facestemp[e*4+1]==v2){
 					exist = true;
@@ -76,5 +88,5 @@
 
 				/*Update chain*/
-				next_edge[nbf] = head_minv[v1];
+				next_face[nbf] = head_minv[v1];
 				head_minv[v1]  = nbf;
 
@@ -87,5 +99,5 @@
 	/*Clean up*/
 	xDelete<int>(head_minv);
-	xDelete<int>(next_edge);
+	xDelete<int>(next_face);
 
 	/*Create final faces*/
@@ -109,3 +121,118 @@
 	iomodel->faces         = faces;
 	iomodel->numberoffaces = nbf;
-}
+}/*}}}*/
+void CreateFaces3d(IoModel* iomodel){/*{{{*/
+
+	/*Intermediaries*/
+	bool exist;
+	int  i,j,k,v0,cols;
+	int  maxnbf,nbf,elementnbf,elementnbv,facemaxnbv;
+	int *elementfaces         = NULL;
+	int *elementfaces_markers = NULL;
+
+	/*Maximum number of faces*/
+	maxnbf = 6*iomodel->numberofelements;
+
+	/*Initialize intermediaries*/
+	int*  facestemp = xNew<int>(maxnbf*6);         /*format: [vertex1 vertex2 vertex3 element1 element2 marker]    */
+	for(i=0;i<maxnbf;i++) facestemp[i*6+4]=-1;     /*Initialize last column of faces as -1 (boundary face)         */
+
+	/*Initialize chain*/
+	int* head_minv = xNew<int>(iomodel->numberofvertices);
+	int* next_face = xNew<int>(maxnbf);
+	for(i=0;i<iomodel->numberofvertices;i++) head_minv[i]=-1;
+
+	/*Mesh specific face indexing per element*/
+	if(iomodel->meshtype==Mesh3DEnum){
+		elementnbv = 6; /*Number of vertices per element*/
+		elementnbf = 5; /*Number of faces per element*/
+		facemaxnbv = 4; /*Maximum number of vertices per face*/
+		cols       = facemaxnbv + 1;
+		elementfaces         = xNew<int>(elementnbf*cols);
+		elementfaces_markers = xNew<int>(elementnbf);
+		/*2 triangles*/
+		elementfaces_markers[0] = 1;
+		elementfaces_markers[1] = 1;
+		elementfaces[cols*0+0] = 3; elementfaces[cols*0+1] = 0;  elementfaces[cols*0+2] = 1; elementfaces[cols*0+3] = 2;
+		elementfaces[cols*1+0] = 3; elementfaces[cols*1+1] = 3;  elementfaces[cols*1+2] = 4; elementfaces[cols*1+3] = 5;
+		/*3 quads*/
+		elementfaces_markers[2] = 2;
+		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;
+	}
+	else{
+		_error_("mesh dimension not supported yet");
+	}
+	int *element_face_connectivity = xNewZeroInit<int>(iomodel->numberofelements*elementnbf); /*format: [face1 face2 ...] */
+
+	/*Initialize number of faces and list of vertex indices*/
+	nbf = 0;
+	int* v = xNew<int>(facemaxnbv);
+
+	for(i=0;i<iomodel->numberofelements;i++){
+		for(j=0;j<elementnbf;j++){
+
+			/*Get indices of current face*/
+			for(k=0;k<elementfaces[cols*j+0];k++){
+				v[k] = iomodel->elements[i*elementnbv + elementfaces[cols*j+k+1]] - 1;
+			}
+
+			/*Sort list of vertices*/
+			HeapSort(v,elementfaces[cols*j+0]);
+			v0 = v[0];
+
+			/*This face a priori has not been processed yet*/
+			exist = false;
+
+			/*Go through all processed faces connected to v1 and check whether we have seen this edge yet*/
+			_assert_(v0>=0 && v0<iomodel->numberofvertices);
+			for(int e=head_minv[v0]; e!=-1; e=next_face[e]){
+				if(facestemp[e*6+1]==v[1]+1 && facestemp[e*6+2]==v[2]+1){
+					exist = true;
+					facestemp[e*6+4]=i+1;
+					element_face_connectivity[i*elementnbf+j]=e;
+					break;
+				}
+			}
+
+			/*If this face is new, add it to the lists*/
+			if(!exist){
+				_assert_(nbf<maxnbf);
+
+				/*Update faces*/
+				facestemp[nbf*6+0] = v[0]+1;
+				facestemp[nbf*6+1] = v[1]+1;
+				facestemp[nbf*6+2] = v[2]+1;
+				facestemp[nbf*6+3] = i+1;
+				facestemp[nbf*6+5] = elementfaces_markers[j];
+
+				/*Update Connectivity*/
+				element_face_connectivity[i*elementnbf+j]=nbf;
+
+				/*Update chain*/
+				next_face[nbf] = head_minv[v0];
+				head_minv[v0]  = nbf;
+
+				/*Increase number of faces*/
+				nbf++;
+			}
+		}
+	}
+
+	/*Clean up*/
+	xDelete<int>(head_minv);
+	xDelete<int>(next_face);
+
+	/*Create final faces (now that we have the correct size)*/
+	int* faces = xNew<int>(nbf*6); /*vertex1 vertex2 vertex3 element1 element2 marker*/
+	xMemCpy<int>(faces,facestemp,nbf*6);
+	xDelete<int>(facestemp);
+
+	/*Assign output pointers*/
+	iomodel->faces                     = faces;
+	iomodel->numberoffaces             = nbf;
+	iomodel->elementtofaceconnectivity = element_face_connectivity;
+}/*}}}*/
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/FacesPartitioning.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/FacesPartitioning.cpp	(revision 17121)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/FacesPartitioning.cpp	(revision 17122)
@@ -11,29 +11,34 @@
 
 	/*Intermediaries*/
-	int  el1,el2;
-	bool my_face;
-
-	/*Check Iomodel properties*/
-	if(iomodel->meshtype!=Mesh2DhorizontalEnum) _error_("only 2d model are supported");
+	int elementnbf;
 
 	/*Get faces and elements*/
 	CreateFaces(iomodel);
+	_assert_(iomodel->elementtofaceconnectivity);
 
+	/*Mesh dependent variables*/
+	if(iomodel->meshtype==Mesh2DhorizontalEnum){
+		elementnbf = 3;
+	}
+	else if(iomodel->meshtype==Mesh2DverticalEnum){
+		elementnbf = 3;
+	}
+	else if(iomodel->meshtype==Mesh3DEnum){
+		elementnbf = 5;
+	}
+	else{
+		_error_("mesh dimension not supported yet");
+	}
 	/*output: */
-	bool* my_faces=xNew<bool>(iomodel->numberoffaces);
+	bool* my_faces=xNewZeroInit<bool>(iomodel->numberoffaces);
 
-	for(int i=0;i<iomodel->numberoffaces;i++){
-
-		/*Get left and right elements*/
-		el1=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
-		el2=iomodel->faces[4*i+3]-1; //faces are [node1 node2 elem1 elem2]
-
-		/*Check whether we should include this face (el2 is -2 for boundary faces)*/
-		my_face = iomodel->my_elements[el1];
-		if(!my_face && el2>=0){
-			my_face = iomodel->my_elements[el2];
+	for(int i=0;i<iomodel->numberofelements;i++){
+		if(iomodel->my_elements[i]){
+			for(int j=0;j<elementnbf;j++){
+				_assert_(iomodel->elementtofaceconnectivity[i*elementnbf+j] >= 0);
+				_assert_(iomodel->elementtofaceconnectivity[i*elementnbf+j] <  iomodel->numberoffaces);
+				my_faces[iomodel->elementtofaceconnectivity[i*elementnbf+j]] = true;
+			}
 		}
-
-		my_faces[i] = my_face;
 	}
 
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h	(revision 17121)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h	(revision 17122)
@@ -30,6 +30,7 @@
 
 /*Mesh properties*/
+void CreateEdges(IoModel* iomodel);
 void CreateFaces(IoModel* iomodel);
-void CreateEdges(IoModel* iomodel);
+void CreateFaces3d(IoModel* iomodel);
 
 /*Connectivity*/