Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15655)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15656)
@@ -3282,4 +3282,5 @@
 			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+iomodel->elements[6*index+0];
 			penta_node_ids[ 7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elements[6*index+1];
@@ -3299,4 +3300,5 @@
 			penta_node_ids[ 5]=iomodel->nodecounter+iomodel->elements[6*index+5];
 			penta_node_ids[ 6]=iomodel->nodecounter+iomodel->numberofvertices+index+1;
+
 			penta_node_ids[ 7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+0];
 			penta_node_ids[ 8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+1];
@@ -3305,4 +3307,31 @@
 			penta_node_ids[11]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+4];
 			penta_node_ids[12]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+5];
+			break;
+		case TaylorHoodEnum:
+			numnodes         = 21;
+			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+iomodel->elementtoedgeconnectivity[9*index+0]+1;
+			penta_node_ids[ 7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1;
+			penta_node_ids[ 8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1;
+			penta_node_ids[ 9]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1;
+			penta_node_ids[10]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1;
+			penta_node_ids[11]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1;
+			penta_node_ids[12]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1;
+			penta_node_ids[13]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1;
+			penta_node_ids[14]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1;
+
+			penta_node_ids[15]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[6*index+0];
+			penta_node_ids[16]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[6*index+1];
+			penta_node_ids[17]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[6*index+2];
+			penta_node_ids[18]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[6*index+3];
+			penta_node_ids[19]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[6*index+4];
+			penta_node_ids[20]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[6*index+5];
+			printarray(penta_node_ids,1,21);
 			break;
 		default:
Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15655)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15656)
@@ -1303,4 +1303,9 @@
 			this->element_type = MINIEnum;
 			return;
+		case TaylorHoodEnum:
+			this->element_type = P2Enum;
+			this->GetNodalFunctions(basis,gauss);
+			this->element_type = TaylorHoodEnum;
+			return;
 		default:
 			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
@@ -1332,4 +1337,9 @@
 			this->GetNodalFunctions(basis,gauss);
 			this->element_type = MINIEnum;
+			return;
+		case TaylorHoodEnum:
+			this->element_type = P1Enum;
+			this->GetNodalFunctions(basis,gauss);
+			this->element_type = TaylorHoodEnum;
 			return;
 		default:
@@ -1396,4 +1406,9 @@
 			this->element_type = MINIEnum;
 			return;
+		case TaylorHoodEnum:
+			this->element_type = P2Enum;
+			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+			this->element_type = TaylorHoodEnum;
+			return;
 		default:
 			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
@@ -1423,4 +1438,9 @@
 			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 			this->element_type = MINIEnum;
+			return;
+		case TaylorHoodEnum:
+			this->element_type = P1Enum;
+			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+			this->element_type = TaylorHoodEnum;
 			return;
 		default:
@@ -1911,4 +1931,5 @@
 		case MINIcondensedEnum: return NUMNODESP1*2;
 		case MINIEnum:          return NUMNODESP1b+NUMNODESP1;
+		case TaylorHoodEnum:    return NUMNODESP2+NUMNODESP1;
 		default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
 	}
@@ -1925,4 +1946,5 @@
 		case MINIcondensedEnum: return NUMNODESP1;
 		case MINIEnum:          return NUMNODESP1;
+		case TaylorHoodEnum:    return NUMNODESP1;
 		default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
 	}
@@ -1939,4 +1961,5 @@
 		case MINIcondensedEnum: return NUMNODESP1b;
 		case MINIEnum:          return NUMNODESP1b;
+		case TaylorHoodEnum:    return NUMNODESP2;
 		default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
 	}
@@ -1955,4 +1978,5 @@
 		case MINIcondensedEnum: return NUMNODESP1;
 		case MINIEnum:          return NUMNODESP1b;
+		case TaylorHoodEnum:    return NUMNODESP2;
 		default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
 	}
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp	(revision 15655)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp	(revision 15656)
@@ -156,8 +156,30 @@
 				}
 			}
-			/*P1+ pressure statically condensed*/
+			/*P1 pressure*/
 			for(i=0;i<iomodel->numberofvertices;i++){
 				if(iomodel->my_vertices[i]){
 					nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,iomodel->numberofvertices+iomodel->numberofelements+i,i,iomodel,analysis,FSpressureEnum));
+				}
+			}
+			break;
+		case TaylorHoodEnum:
+			_assert_(approximation==FSApproximationEnum);
+			/*P2 velocity*/
+			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,i,iomodel,analysis,FSvelocityEnum));
+				}
+			}
+			for(i=0;i<iomodel->numberofedges;i++){
+				if(my_edges[i]){
+					nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,0,iomodel,analysis,FSvelocityEnum));
+				}
+			}
+
+			/*P1 pressure*/
+			for(i=0;i<iomodel->numberofvertices;i++){
+				if(iomodel->my_vertices[i]){
+					nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,i,iomodel,analysis,FSpressureEnum));
 				}
 			}
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp	(revision 15655)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp	(revision 15656)
@@ -96,8 +96,9 @@
 			iomodel->Constant(&temp,FlowequationFeFSEnum);
 			switch(temp){
-				case 0 : finiteelement = P1Enum; break;//P1P1
-				case 1 : finiteelement = P1Enum; break;//P1P1GSL
-				case 2 : finiteelement = P1Enum; break;//MINIcondensedEnum
-				case 3 : finiteelement = P1bubbleEnum; break;//MINIEnum
+				case 0 : finiteelement = P1Enum;       break;//P1P1
+				case 1 : finiteelement = P1Enum;       break;//P1P1GSL
+				case 2 : finiteelement = P1Enum;       break;//MINIcondensed
+				case 3 : finiteelement = P1bubbleEnum; break;//MINI
+				case 4 : finiteelement = P2Enum;       break;//TaylorHood (P2P1)
 				default: _error_("finite element "<<temp<<" not supported");
 			}
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp	(revision 15655)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp	(revision 15656)
@@ -66,4 +66,5 @@
 				case 2 : finiteelement = MINIcondensedEnum; break;
 				case 3 : finiteelement = MINIEnum;          break;
+				case 4 : finiteelement = TaylorHoodEnum;    break;
 				default: _error_("finite element "<<temp<<" not supported");
 			}
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp	(revision 15655)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp	(revision 15656)
@@ -70,4 +70,5 @@
 				case 2 : finiteelement = MINIcondensedEnum; break;
 				case 3 : finiteelement = MINIEnum;          break;
+				case 4 : finiteelement = TaylorHoodEnum;    break;
 				default: _error_("finite element "<<temp<<" not supported");
 			}
Index: /issm/trunk-jpl/src/m/classes/flowequation.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/flowequation.m	(revision 15655)
+++ /issm/trunk-jpl/src/m/classes/flowequation.m	(revision 15656)
@@ -83,5 +83,5 @@
 				md = checkfield(md,'flowequation.fe_SSA','numel',[1],'values',[0 1]);
 				md = checkfield(md,'flowequation.fe_HO','numel',[1],'values',[0:3]);
-				md = checkfield(md,'flowequation.fe_FS','numel',[1],'values',[0:3]);
+				md = checkfield(md,'flowequation.fe_FS','numel',[1],'values',[0:4]);
 				md = checkfield(md,'flowequation.borderSSA','size',[md.mesh.numberofvertices 1],'values',[0 1]);
 				md = checkfield(md,'flowequation.borderHO','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 15655)
+++ /issm/trunk-jpl/src/m/classes/flowequation.py	(revision 15656)
@@ -68,5 +68,5 @@
 			md = checkfield(md,'flowequation.fe_SSA','numel',[1],'values',[0,1])
 			md = checkfield(md,'flowequation.fe_HO','numel',[1],'values',numpy.arange(0,3+1))
-			md = checkfield(md,'flowequation.fe_FS','numel',[1],'values',numpy.arange(0,3+1))
+			md = checkfield(md,'flowequation.fe_FS','numel',[1],'values',numpy.arange(0,4+1))
 			md = checkfield(md,'flowequation.borderSSA','size',[md.mesh.numberofvertices],'values',[0,1])
 			md = checkfield(md,'flowequation.borderHO','size',[md.mesh.numberofvertices],'values',[0,1])