Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17040)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17041)
@@ -2769,4 +2769,20 @@
 			penta_node_ids[ 8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1;
 			break;
+		case P1xP3Enum:
+			numnodes         = 12;
+			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+2*iomodel->elementtoedgeconnectivity[9*index+0]+1;
+			penta_node_ids[ 7]=iomodel->nodecounter+iomodel->numberofvertices+2*iomodel->elementtoedgeconnectivity[9*index+1]+1;
+			penta_node_ids[ 8]=iomodel->nodecounter+iomodel->numberofvertices+2*iomodel->elementtoedgeconnectivity[9*index+2]+1;
+			penta_node_ids[ 9]=iomodel->nodecounter+iomodel->numberofvertices+2*iomodel->elementtoedgeconnectivity[9*index+0]+2;
+			penta_node_ids[10]=iomodel->nodecounter+iomodel->numberofvertices+2*iomodel->elementtoedgeconnectivity[9*index+1]+2;
+			penta_node_ids[11]=iomodel->nodecounter+iomodel->numberofvertices+2*iomodel->elementtoedgeconnectivity[9*index+2]+2;
+			break;
 		case P2xP1Enum:
 			numnodes         = 12;
Index: /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp	(revision 17040)
+++ /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp	(revision 17041)
@@ -587,4 +587,22 @@
 			}
 			break;
+		case P1xP3Enum: 
+			switch(iv){
+				case 0 : coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
+				case 1 : coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
+				case 2 : coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
+				case 3 : coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
+				case 4 : coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
+				case 5 : coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
+
+				case 6 : coord1=1.; coord2=0.; coord3=0.; coord4=-1./3.; break;
+				case 7 : coord1=0.; coord2=1.; coord3=0.; coord4=-1./3.; break;
+				case 8 : coord1=0.; coord2=0.; coord3=1.; coord4=-1./3.; break;
+				case 9 : coord1=1.; coord2=0.; coord3=0.; coord4=+1./3.; break;
+				case 10: coord1=0.; coord2=1.; coord3=0.; coord4=+1./3.; break;
+				case 11: coord1=0.; coord2=0.; coord3=1.; coord4=+1./3.; break;
+				default: _error_("node index should be in [0 11]");
+			}
+			break;
 		case P2xP1Enum: 
 			switch(iv){
Index: /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp	(revision 17040)
+++ /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp	(revision 17041)
@@ -61,8 +61,14 @@
 			EdgesPartitioning(&my_edges,iomodel);
 			break;
+		case P1xP3Enum:
+			EdgesPartitioning(&my_edges,iomodel);
+			break;
 		case P2xP1Enum:
 			EdgesPartitioning(&my_edges,iomodel);
 			break;
 		case P2Enum:
+			EdgesPartitioning(&my_edges,iomodel);
+			break;
+		case P2xP4Enum:
 			EdgesPartitioning(&my_edges,iomodel);
 			break;
@@ -160,4 +166,29 @@
 				}
 				break;
+			case P1xP3Enum:
+				for(i=0;i<iomodel->numberofvertices;i++){
+					if((iomodel->my_vertices[i])){
+						if (!xIsNan<IssmDouble>(spcdata[i])){
+							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,dof,spcdata[i],analysis_type));
+							count++;
+						}
+					}
+				}
+				for(i=0;i<iomodel->numberofedges;i++){
+					if(iomodel->edges[i*3+2]==2){
+						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+2*count+1,iomodel->nodecounter+iomodel->numberofvertices+2*i+1,
+												dof,2./3.*spcdata[v1]+1./3.*spcdata[v2],analysis_type));
+								constraints->AddObject(new SpcStatic(iomodel->constraintcounter+2*count+2,iomodel->nodecounter+iomodel->numberofvertices+2*i+2,
+												dof,1./3.*spcdata[v1]+2./3.*spcdata[v2],analysis_type));
+								count=count+2;
+							}
+						}
+					}
+				}
+				break;
 			case P2xP1Enum:
 				for(i=0;i<iomodel->numberofvertices;i++){
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp	(revision 17040)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp	(revision 17041)
@@ -85,4 +85,27 @@
 			}
 			break;
+
+		case P1xP3Enum:
+			EdgesPartitioning(&my_edges,iomodel);
+			for(i=0;i<iomodel->numberofvertices;i++){
+				if(iomodel->my_vertices[i]){
+					nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,analysis,approximation));
+				}
+			}
+
+			counter = iomodel->numberofvertices;
+			for(i=0;i<iomodel->numberofedges;i++){
+				if(iomodel->edges[i*3+2]==2){
+					if(my_edges[i]){
+						node = new Node(iomodel->nodecounter+iomodel->numberofvertices+2*i+1,counter+1,lid++,0,iomodel,analysis,approximation);
+						nodes->AddObject(node);
+						node = new Node(iomodel->nodecounter+iomodel->numberofvertices+2*i+2,counter+2,lid++,0,iomodel,analysis,approximation);
+						nodes->AddObject(node);
+					}
+					counter=counter+2;
+				}
+			}
+			break;
+
 		case P2xP1Enum:
 			EdgesPartitioning(&my_edges,iomodel);
@@ -213,4 +236,11 @@
 				}
 			}
+			break;
+		case OneLayerP4zEnum:
+			_assert_(approximation==FSApproximationEnum);
+			/*P2xP4 velocity*/
+			EdgesPartitioning(&my_edges,iomodel);
+			/*P1xP3 pressure*/
+			_error_("Not supported");
 			break;
 
Index: /issm/trunk-jpl/src/m/classes/flowequation.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/flowequation.m	(revision 17040)
+++ /issm/trunk-jpl/src/m/classes/flowequation.m	(revision 17041)
@@ -86,5 +86,5 @@
 				md = checkfield(md,'fieldname','flowequation.isFS','numel',[1],'values',[0 1]);
 				md = checkfield(md,'fieldname','flowequation.fe_SSA','values',{'P1','P1bubble','P1bubblecondensed','P2'});
-				md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',{'P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2'});
+				md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',{'P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P1xP3','P2xP4'});
 				md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',{'P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','OneLayerP4z'});
 				md = checkfield(md,'fieldname','flowequation.borderSSA','size',[md.mesh.numberofvertices 1],'values',[0 1]);
Index: /issm/trunk-jpl/src/m/classes/flowequation.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/flowequation.py	(revision 17040)
+++ /issm/trunk-jpl/src/m/classes/flowequation.py	(revision 17041)
@@ -75,5 +75,5 @@
 			md = checkfield(md,'fieldname','flowequation.fe_SSA','values',['P1','P1bubble','P1bubblecondensed','P2'])
 			md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',['P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2'])
-			md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',['P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','OneLayerP4z'])
+			md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',['P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','OneLayerP4z','P1xP3','P2xP4'])
 			md = checkfield(md,'fieldname','flowequation.borderSSA','size',[md.mesh.numberofvertices],'values',[0,1])
 			md = checkfield(md,'fieldname','flowequation.borderHO','size',[md.mesh.numberofvertices],'values',[0,1])