Index: ../trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp	(revision 15626)
+++ ../trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp	(revision 15627)
@@ -55,8 +55,9 @@
 		else if(isHO){
 			iomodel->Constant(&temp,FlowequationFeHOEnum);
 			switch(temp){
-				case 0 : finiteelement = P1Enum; break;
-				case 1 : finiteelement = P2Enum; break;
+				case 0 : finiteelement = P1Enum;    break;
+				case 1 : finiteelement = P2Enum;    break;
+				case 2 : finiteelement = P2xP1Enum; break;
 				default: _error_("finite element "<<temp<<" not supported");
 			}
 		}
Index: ../trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp	(revision 15626)
+++ ../trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp	(revision 15627)
@@ -84,8 +84,9 @@
 		else if(isHO){
 			iomodel->Constant(&temp,FlowequationFeHOEnum);
 			switch(temp){
-				case 0 : finiteelement = P1Enum; break;
-				case 1 : finiteelement = P2Enum; break;
+				case 0 : finiteelement = P1Enum;    break;
+				case 1 : finiteelement = P2Enum;    break;
+				case 2 : finiteelement = P2xP1Enum; break;
 				default: _error_("finite element "<<temp<<" not supported");
 			}
 		}
Index: ../trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp	(revision 15626)
+++ ../trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp	(revision 15627)
@@ -50,8 +50,9 @@
 			approximation = HOApproximationEnum;
 			iomodel->Constant(&temp,FlowequationFeHOEnum);
 			switch(temp){
-				case 0 : finiteelement = P1Enum; break;
-				case 1 : finiteelement = P2Enum; break;
+				case 0 : finiteelement = P1Enum;    break;
+				case 1 : finiteelement = P2Enum;    break;
+				case 2 : finiteelement = P2xP1Enum; break;
 				default: _error_("finite element "<<temp<<" not supported");
 			}
 		}
Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp	(revision 15626)
+++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp	(revision 15627)
@@ -42,6 +42,25 @@
 			}
 			break;
 
+		case P2xP1Enum:
+			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,approximation));
+				}
+			}
+			int numberofvertices2d;
+			Node* node = NULL;
+			iomodel->Constant(&numberofvertices2d,MeshNumberofvertices2dEnum);
+			/*Create all P2 nodes anyway*/
+			for(i=0;i<iomodel->numberofedges;i++){
+				if(my_edges[i]){
+					node = new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,0,iomodel,analysis,approximation);
+					nodes->AddObject(node);
+				}
+			}
+			break;
+
 		case P2Enum:
 			EdgesPartitioning(&my_edges,iomodel);
 			for(i=0;i<iomodel->numberofvertices;i++){
Index: ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15626)
+++ ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15627)
@@ -1624,9 +1624,10 @@
 int PentaRef::NumberofNodes(void){
 
 	switch(this->element_type){
-		case P1Enum:   return NUMNODESP1;
-		case P2Enum:   return NUMNODESP2;
-		case MINIEnum: return NUMNODESP1b;
+		case P1Enum:    return NUMNODESP1;
+		case P2Enum:    return NUMNODESP2;
+		case P2xP1Enum: return NUMNODESP2xP1;
+		case MINIEnum:  return NUMNODESP1b;
 		default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
 	}
 
Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15626)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15627)
@@ -3177,6 +3177,22 @@
 			penta_node_ids[4]=iomodel->nodecounter+iomodel->elements[6*index+4];
 			penta_node_ids[5]=iomodel->nodecounter+iomodel->elements[6*index+5];
 			break;
+		case P2xP1Enum:
+			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+iomodel->elementtoedgeconnectivity[9*index+3]+1;
+			penta_node_ids[ 7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1;
+			penta_node_ids[ 8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1;
+			penta_node_ids[ 9]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1;
+			penta_node_ids[10]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1;
+			penta_node_ids[11]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1;
+			break;
 		case P2Enum:
 			numnodes         = 15;
 			penta_node_ids   = xNew<int>(numnodes);
@@ -6974,6 +6990,7 @@
 	xDelete<IssmDouble>(D);
 	xDelete<IssmDouble>(Bprime);
 	xDelete<IssmDouble>(B);
+	if(id==1) Ke->Echo();
 	return Ke;
 }
 /*}}}*/
Index: ../trunk-jpl/src/c/classes/gauss/GaussPenta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/gauss/GaussPenta.cpp	(revision 15626)
+++ ../trunk-jpl/src/c/classes/gauss/GaussPenta.cpp	(revision 15627)
@@ -414,6 +414,24 @@
 				default: _error_("node index should be in [0 5]");
 			}
 			break;
+		case 12: //P2xP1 Lagrange element
+			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=0.; coord2=.5; coord3=.5; coord4=-1.;break;
+				case  7: coord1=.5; coord2=0.; coord3=.5; coord4=-1.;break;
+				case  8: coord1=.5; coord2=.5; coord3=0.; coord4=-1.;break;
+				case  9: coord1=0.; coord2=.5; coord3=.5; coord4=+1.;break;
+				case 10: coord1=.5; coord2=0.; coord3=.5; coord4=+1.;break;
+				case 11: coord1=.5; coord2=.5; coord3=0.; coord4=+1.;break;
+				default: _error_("node index should be in [0 5]");
+			}
+			break;
 		case 15: //P2 Lagrange element
 			switch(iv){
 				case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
Index: ../trunk-jpl/src/m/classes/flowequation.m
===================================================================
--- ../trunk-jpl/src/m/classes/flowequation.m	(revision 15626)
+++ ../trunk-jpl/src/m/classes/flowequation.m	(revision 15627)
@@ -79,7 +79,7 @@
 				md = checkfield(md,'flowequation.isHO','numel',[1],'values',[0 1]);
 				md = checkfield(md,'flowequation.isFS','numel',[1],'values',[0 1]);
 				md = checkfield(md,'flowequation.fe_SSA','numel',[1],'values',[0 1]);
-				md = checkfield(md,'flowequation.fe_HO','numel',[1],'values',[0 1]);
+				md = checkfield(md,'flowequation.fe_HO','numel',[1],'values',[0 1 2]);
 				md = checkfield(md,'flowequation.fe_FS','numel',[1],'values',[0]);
 				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]);