Index: ../trunk-jpl/src/c/classes/Elements/TriaRef.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 18924)
+++ ../trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 18925)
@@ -15,23 +15,22 @@
 		~TriaRef();
 
 		/*Numerics*/
+		void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement);
+		void GetInputValue(IssmDouble* pp, IssmDouble* plist, Gauss* gauss,int finiteelement);
 		void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss);
-		void GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss);
 		void GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss);
 		void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,Gauss* gauss);
 		void GetNodalFunctions(IssmDouble* basis,Gauss* gauss,int finiteelement);
-		void GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss, int index1,int index2,int finiteelement);
+		void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss,int finiteelement);
+		void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss,int finiteelement);
 		void GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2,int finiteelement);
 		void GetSegmentBprimeFlux(IssmDouble* Bprime,Gauss* gauss, int index1,int index2,int finiteelement);
-		void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss,int finiteelement);
-		void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss,int finiteelement);
-		void GetInputValue(IssmDouble* pp, IssmDouble* plist, Gauss* gauss,int finiteelement);
-		void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement);
-
+		void GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss);
+		void GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss, int index1,int index2,int finiteelement);
 		void NodeOnEdgeIndices(int* pnumindices,int** pindices,int index,int finiteelement);
 		int  NumberofNodes(int finiteelement);
-		int  VelocityInterpolation(int fe_stokes);
 		int  PressureInterpolation(int fe_stokes);
 		int  TensorInterpolation(int fe_stokes);
+		int  VelocityInterpolation(int fe_stokes);
 };
 #endif
Index: ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 18924)
+++ ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 18925)
@@ -36,6 +36,150 @@
 /*}}}*/
 
 /*Reference Element numerics*/
+void PentaRef::BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){/*{{{*/
+
+	/*Output*/
+	int  numindices;
+	int* indices = NULL;
+
+	switch(finiteelement){
+		case P1Enum: case P1DGEnum:
+			numindices = 3;
+			indices    = xNew<int>(numindices);
+			indices[0] = 0;
+			indices[1] = 1;
+			indices[2] = 2;
+			break;
+		case P1bubbleEnum: case P1bubblecondensedEnum:
+			numindices = 3;
+			indices    = xNew<int>(numindices);
+			indices[0] = 0;
+			indices[1] = 1;
+			indices[2] = 2;
+			break;
+		case P2xP1Enum:
+			numindices = 6;
+			indices    = xNew<int>(numindices);
+			indices[0] = 0;
+			indices[1] = 1;
+			indices[2] = 2;
+			indices[3] = 6;
+			indices[4] = 7;
+			indices[5] = 8;
+			break;
+		case P1xP2Enum:
+			numindices = 3;
+			indices    = xNew<int>(numindices);
+			indices[0] = 0;
+			indices[1] = 1;
+			indices[2] = 2;
+			break;
+		case P1xP3Enum:
+			numindices = 3;
+			indices    = xNew<int>(numindices);
+			indices[0] = 0;
+			indices[1] = 1;
+			indices[2] = 2;
+			break;
+		case P2Enum:
+			numindices = 6;
+			indices    = xNew<int>(numindices);
+			indices[0] = 0;
+			indices[1] = 1;
+			indices[2] = 2;
+			indices[3] = 9;
+			indices[4] = 10;
+			indices[5] = 11;
+			break;
+		case P2bubbleEnum:
+			numindices = 6;
+			indices    = xNew<int>(numindices);
+			indices[0] = 0;
+			indices[1] = 1;
+			indices[2] = 2;
+			indices[3] = 9;
+			indices[4] = 10;
+			indices[5] = 11;
+			break;
+		case P2xP4Enum:
+			numindices = 6;
+			indices    = xNew<int>(numindices);
+			indices[0] = 0;
+			indices[1] = 1;
+			indices[2] = 2;
+			indices[3] = 9;
+			indices[4] = 10;
+			indices[5] = 11;
+			break;
+		default:
+			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
+	}
+
+	/*Assign output pointer*/
+	*pnumindices = numindices;
+	*pindices    = indices;
+}
+/*}}}*/
+void PentaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement){/*{{{*/
+	/*From node values of parameter p (p_list[0], p_list[1], p_list[2],
+	 * p_list[3], p_list[4] and p_list[4]), return parameter derivative value at
+	 * gaussian point specified by gauss_coord:
+	 *   dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
+	 *   dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy;
+	 *   dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz;
+	 *
+	 *   p is a vector of size 3x1 already allocated.
+	 *
+	 * WARNING: For a significant gain in performance, it is better to use
+	 * static memory allocation instead of dynamic.
+	 */
+
+	/*Allocate derivatives of basis functions*/
+	IssmDouble  dbasis[3*NUMNODESMAX];
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes(finiteelement);
+	_assert_(numnodes<=NUMNODESMAX);
+
+	/*Get basis functions derivatives at this point*/
+	GetNodalFunctionsDerivatives(&dbasis[0],xyz_list,gauss,finiteelement);
+
+	/*Calculate parameter for this Gauss point*/
+	IssmDouble dpx=0.;
+	IssmDouble dpy=0.;
+	IssmDouble dpz=0.;
+	for(int i=0;i<numnodes;i++) dpx += dbasis[0*numnodes+i]*plist[i];
+	for(int i=0;i<numnodes;i++) dpy += dbasis[1*numnodes+i]*plist[i];
+	for(int i=0;i<numnodes;i++) dpz += dbasis[2*numnodes+i]*plist[i];
+
+	/*Assign values*/
+	p[0]=dpx;
+	p[1]=dpy;
+	p[2]=dpz;
+}
+/*}}}*/
+void PentaRef::GetInputValue(IssmDouble* pvalue,IssmDouble* plist,Gauss* gauss,int finiteelement){/*{{{*/
+	/* WARNING: For a significant gain in performance, it is better to use
+	 * static memory allocation instead of dynamic.*/
+
+	/*Allocate basis functions*/
+	IssmDouble  basis[NUMNODESMAX];
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes(finiteelement);
+	_assert_(numnodes<=NUMNODESMAX);
+
+	/*Get basis functions at this point*/
+	GetNodalFunctions(&basis[0],gauss,finiteelement);
+
+	/*Calculate parameter for this Gauss point*/
+	IssmDouble value =0.;
+	for(int i=0;i<numnodes;i++) value += basis[i]*plist[i];
+
+	/*Assign output pointer*/
+	*pvalue = value;
+}
+/*}}}*/
 void PentaRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/
 	/*The Jacobian is constant over the element, discard the gaussian points. 
 	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
@@ -106,41 +250,6 @@
 
 }
 /*}}}*/
-void PentaRef::GetTriaJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
-	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
-
-	IssmDouble x1=xyz_list[3*0+0];
-	IssmDouble y1=xyz_list[3*0+1];
-	IssmDouble z1=xyz_list[3*0+2];
-	IssmDouble x2=xyz_list[3*1+0];
-	IssmDouble y2=xyz_list[3*1+1];
-	IssmDouble z2=xyz_list[3*1+2];
-	IssmDouble x3=xyz_list[3*2+0];
-	IssmDouble y3=xyz_list[3*2+1];
-	IssmDouble z3=xyz_list[3*2+2];
-
-	/*Jdet = norm( AB ^ AC ) / (2 * area of the reference triangle), with areaRef=sqrt(3) */
-	*Jdet=SQRT3/6.*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2),0.5);
-	if(*Jdet<0) _error_("negative jacobian determinant!");
-}
-/*}}}*/
-void PentaRef::GetSegmentJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
-	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
-
-	IssmDouble x1=xyz_list[3*0+0];
-	IssmDouble y1=xyz_list[3*0+1];
-	IssmDouble z1=xyz_list[3*0+2];
-	IssmDouble x2=xyz_list[3*1+0];
-	IssmDouble y2=xyz_list[3*1+1];
-	IssmDouble z2=xyz_list[3*1+2];
-
-	*Jdet=.5*sqrt(pow(x2-x1,2) + pow(y2-y1,2) + pow(z2-z1,2));
-	if(*Jdet<0) _error_("negative jacobian determinant!");
-
-}
-/*}}}*/
 void PentaRef::GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 
 	/*Jacobian*/
@@ -885,64 +994,39 @@
 
 }
 /*}}}*/
-void PentaRef::GetInputValue(IssmDouble* pvalue,IssmDouble* plist,Gauss* gauss,int finiteelement){/*{{{*/
-	/* WARNING: For a significant gain in performance, it is better to use
-	 * static memory allocation instead of dynamic.*/
+void PentaRef::GetSegmentJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
+	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
 
-	/*Allocate basis functions*/
-	IssmDouble  basis[NUMNODESMAX];
+	IssmDouble x1=xyz_list[3*0+0];
+	IssmDouble y1=xyz_list[3*0+1];
+	IssmDouble z1=xyz_list[3*0+2];
+	IssmDouble x2=xyz_list[3*1+0];
+	IssmDouble y2=xyz_list[3*1+1];
+	IssmDouble z2=xyz_list[3*1+2];
 
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes(finiteelement);
-	_assert_(numnodes<=NUMNODESMAX);
+	*Jdet=.5*sqrt(pow(x2-x1,2) + pow(y2-y1,2) + pow(z2-z1,2));
+	if(*Jdet<0) _error_("negative jacobian determinant!");
 
-	/*Get basis functions at this point*/
-	GetNodalFunctions(&basis[0],gauss,finiteelement);
-
-	/*Calculate parameter for this Gauss point*/
-	IssmDouble value =0.;
-	for(int i=0;i<numnodes;i++) value += basis[i]*plist[i];
-
-	/*Assign output pointer*/
-	*pvalue = value;
 }
 /*}}}*/
-void PentaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement){/*{{{*/
-	/*From node values of parameter p (p_list[0], p_list[1], p_list[2],
-	 * p_list[3], p_list[4] and p_list[4]), return parameter derivative value at
-	 * gaussian point specified by gauss_coord:
-	 *   dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
-	 *   dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy;
-	 *   dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz;
-	 *
-	 *   p is a vector of size 3x1 already allocated.
-	 *
-	 * WARNING: For a significant gain in performance, it is better to use
-	 * static memory allocation instead of dynamic.
-	 */
+void PentaRef::GetTriaJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
+	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
 
-	/*Allocate derivatives of basis functions*/
-	IssmDouble  dbasis[3*NUMNODESMAX];
+	IssmDouble x1=xyz_list[3*0+0];
+	IssmDouble y1=xyz_list[3*0+1];
+	IssmDouble z1=xyz_list[3*0+2];
+	IssmDouble x2=xyz_list[3*1+0];
+	IssmDouble y2=xyz_list[3*1+1];
+	IssmDouble z2=xyz_list[3*1+2];
+	IssmDouble x3=xyz_list[3*2+0];
+	IssmDouble y3=xyz_list[3*2+1];
+	IssmDouble z3=xyz_list[3*2+2];
 
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes(finiteelement);
-	_assert_(numnodes<=NUMNODESMAX);
-
-	/*Get basis functions derivatives at this point*/
-	GetNodalFunctionsDerivatives(&dbasis[0],xyz_list,gauss,finiteelement);
-
-	/*Calculate parameter for this Gauss point*/
-	IssmDouble dpx=0.;
-	IssmDouble dpy=0.;
-	IssmDouble dpz=0.;
-	for(int i=0;i<numnodes;i++) dpx += dbasis[0*numnodes+i]*plist[i];
-	for(int i=0;i<numnodes;i++) dpy += dbasis[1*numnodes+i]*plist[i];
-	for(int i=0;i<numnodes;i++) dpz += dbasis[2*numnodes+i]*plist[i];
-
-	/*Assign values*/
-	p[0]=dpx;
-	p[1]=dpy;
-	p[2]=dpz;
+	/*Jdet = norm( AB ^ AC ) / (2 * area of the reference triangle), with areaRef=sqrt(3) */
+	*Jdet=SQRT3/6.*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2),0.5);
+	if(*Jdet<0) _error_("negative jacobian determinant!");
 }
 /*}}}*/
 int  PentaRef::NumberofNodes(int finiteelement){/*{{{*/
@@ -976,24 +1060,6 @@
 	return -1;
 }
 /*}}}*/
-int  PentaRef::VelocityInterpolation(int fe_stokes){/*{{{*/
-
-	switch(fe_stokes){
-		case P1P1Enum:              return P1Enum;
-		case P1P1GLSEnum:           return P1Enum;
-		case MINIcondensedEnum:     return P1bubbleEnum;
-		case MINIEnum:              return P1bubbleEnum;
-		case TaylorHoodEnum:        return P2Enum;
-		case LATaylorHoodEnum:      return P2Enum;
-		case OneLayerP4zEnum:       return P2xP4Enum;
-		case CrouzeixRaviartEnum:   return P2bubbleEnum;
-		case LACrouzeixRaviartEnum: return P2bubbleEnum;
-		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
-	}
-
-	return -1;
-}
-/*}}}*/
 int  PentaRef::PressureInterpolation(int fe_stokes){/*{{{*/
 
 	switch(fe_stokes){
@@ -1012,100 +1078,6 @@
 	return -1;
 }
 /*}}}*/
-int  PentaRef::TensorInterpolation(int fe_stokes){/*{{{*/
-
-	switch(fe_stokes){
-		case XTaylorHoodEnum:    return P1DGEnum;
-		default: _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
-	}
-
-	return -1;
-}
-/*}}}*/
-void PentaRef::BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){/*{{{*/
-
-	/*Output*/
-	int  numindices;
-	int* indices = NULL;
-
-	switch(finiteelement){
-		case P1Enum: case P1DGEnum:
-			numindices = 3;
-			indices    = xNew<int>(numindices);
-			indices[0] = 0;
-			indices[1] = 1;
-			indices[2] = 2;
-			break;
-		case P1bubbleEnum: case P1bubblecondensedEnum:
-			numindices = 3;
-			indices    = xNew<int>(numindices);
-			indices[0] = 0;
-			indices[1] = 1;
-			indices[2] = 2;
-			break;
-		case P2xP1Enum:
-			numindices = 6;
-			indices    = xNew<int>(numindices);
-			indices[0] = 0;
-			indices[1] = 1;
-			indices[2] = 2;
-			indices[3] = 6;
-			indices[4] = 7;
-			indices[5] = 8;
-			break;
-		case P1xP2Enum:
-			numindices = 3;
-			indices    = xNew<int>(numindices);
-			indices[0] = 0;
-			indices[1] = 1;
-			indices[2] = 2;
-			break;
-		case P1xP3Enum:
-			numindices = 3;
-			indices    = xNew<int>(numindices);
-			indices[0] = 0;
-			indices[1] = 1;
-			indices[2] = 2;
-			break;
-		case P2Enum:
-			numindices = 6;
-			indices    = xNew<int>(numindices);
-			indices[0] = 0;
-			indices[1] = 1;
-			indices[2] = 2;
-			indices[3] = 9;
-			indices[4] = 10;
-			indices[5] = 11;
-			break;
-		case P2bubbleEnum:
-			numindices = 6;
-			indices    = xNew<int>(numindices);
-			indices[0] = 0;
-			indices[1] = 1;
-			indices[2] = 2;
-			indices[3] = 9;
-			indices[4] = 10;
-			indices[5] = 11;
-			break;
-		case P2xP4Enum:
-			numindices = 6;
-			indices    = xNew<int>(numindices);
-			indices[0] = 0;
-			indices[1] = 1;
-			indices[2] = 2;
-			indices[3] = 9;
-			indices[4] = 10;
-			indices[5] = 11;
-			break;
-		default:
-			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
-	}
-
-	/*Assign output pointer*/
-	*pnumindices = numindices;
-	*pindices    = indices;
-}
-/*}}}*/
 void PentaRef::SurfaceNodeIndices(int* pnumindices,int** pindices,int finiteelement){/*{{{*/
 
 	/*Output*/
@@ -1163,3 +1135,31 @@
 	*pindices    = indices;
 }
 /*}}}*/
+int  PentaRef::TensorInterpolation(int fe_stokes){/*{{{*/
+
+	switch(fe_stokes){
+		case XTaylorHoodEnum:    return P1DGEnum;
+		default: _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
+	}
+
+	return -1;
+}
+/*}}}*/
+int  PentaRef::VelocityInterpolation(int fe_stokes){/*{{{*/
+
+	switch(fe_stokes){
+		case P1P1Enum:              return P1Enum;
+		case P1P1GLSEnum:           return P1Enum;
+		case MINIcondensedEnum:     return P1bubbleEnum;
+		case MINIEnum:              return P1bubbleEnum;
+		case TaylorHoodEnum:        return P2Enum;
+		case LATaylorHoodEnum:      return P2Enum;
+		case OneLayerP4zEnum:       return P2xP4Enum;
+		case CrouzeixRaviartEnum:   return P2bubbleEnum;
+		case LACrouzeixRaviartEnum: return P2bubbleEnum;
+		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
+	}
+
+	return -1;
+}
+/*}}}*/
Index: ../trunk-jpl/src/c/classes/Elements/PentaRef.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 18924)
+++ ../trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 18925)
@@ -14,24 +14,23 @@
 		~PentaRef();
 
 		/*Numerics*/
+		void BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement);
+		void GetInputDerivativeValue(IssmDouble* pvalues, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement);
+		void GetInputValue(IssmDouble* pvalue,IssmDouble* plist, Gauss* gauss,int finiteelement);
+		void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss);
+		void GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
+		void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,Gauss* gauss);
+		void GetLprimeFSSSA(IssmDouble* LprimeFSSSA, IssmDouble* xyz_list, Gauss* gauss);
 		void GetNodalFunctions(IssmDouble* basis, Gauss* gauss,int finiteelement);
 		void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss,int finiteelement);
 		void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss,int finiteelement);
 		void GetQuadJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
-		void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss);
-		void GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
+		void GetSegmentJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
 		void GetTriaJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
-		void GetSegmentJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
-		void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,Gauss* gauss);
-		void GetLprimeFSSSA(IssmDouble* LprimeFSSSA, IssmDouble* xyz_list, Gauss* gauss);
-		void GetInputValue(IssmDouble* pvalue,IssmDouble* plist, Gauss* gauss,int finiteelement);
-		void GetInputDerivativeValue(IssmDouble* pvalues, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement);
-
-		void BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement);
-		void SurfaceNodeIndices(int* pnumindices,int** pindices,int finiteelement);
 		int  NumberofNodes(int finiteelement);
-		int  VelocityInterpolation(int fe_stokes);
 		int  PressureInterpolation(int fe_stokes);
+		void SurfaceNodeIndices(int* pnumindices,int** pindices,int finiteelement);
 		int  TensorInterpolation(int fe_stokes);
+		int  VelocityInterpolation(int fe_stokes);
 };
 #endif
Index: ../trunk-jpl/src/c/classes/Elements/ElementHook.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/ElementHook.cpp	(revision 18924)
+++ ../trunk-jpl/src/c/classes/Elements/ElementHook.cpp	(revision 18925)
@@ -74,24 +74,22 @@
 }
 /*}}}*/
 
+void ElementHook::InitHookNeighbors(int* element_ids){/*{{{*/
+	this->hneighbors=new Hook(element_ids,2);
+}
+/*}}}*/
 void ElementHook::SetHookNodes(int* node_ids,int numnodes,int analysis_counter){/*{{{*/
 	if(this->hnodes) this->hnodes[analysis_counter]= new Hook(node_ids,numnodes);
 }
 /*}}}*/
-void ElementHook::InitHookNeighbors(int* element_ids){/*{{{*/
-	this->hneighbors=new Hook(element_ids,2);
-}
-/*}}}*/
-void ElementHook::SpawnTriaHook(ElementHook* triahook,int index1,int index2,int index3){/*{{{*/
+void ElementHook::SpawnSegHook(ElementHook* triahook,int index1,int index2){/*{{{*/
 
-	/*Create arrow of indices depending on location (0=base 1=surface)*/
-	int indices[3];
-	indices[0] = index1;
-	indices[1] = index2;
-	indices[2] = index3;
-
 	triahook->numanalyses=this->numanalyses;
 
+	int indices[2];
+	indices[0]=index1;
+	indices[1]=index2;
+
 	/*Spawn nodes hook*/
 	triahook->hnodes=new Hook*[this->numanalyses];
 	for(int i=0;i<this->numanalyses;i++){
@@ -100,24 +98,26 @@
 			triahook->hnodes[i]=NULL;
 		}
 		else{
-			triahook->hnodes[i]=this->hnodes[i]->Spawn(indices,3);
+			triahook->hnodes[i]=this->hnodes[i]->Spawn(indices,2);
 		}
 	}
 
-	/*do not spawn hmaterial. material will be taken care of by Penta*/
+	/*do not spawn hmaterial. material will be taken care of by Tria*/
 	triahook->hmaterial=NULL;
-	triahook->hvertices=(Hook*)this->hvertices->Spawn(indices,3);
+	triahook->hvertices=(Hook*)this->hvertices->Spawn(indices,2);
 	triahook->hmatpar=(Hook*)this->hmatpar->copy();
 }
 /*}}}*/
-void ElementHook::SpawnSegHook(ElementHook* triahook,int index1,int index2){/*{{{*/
+void ElementHook::SpawnTriaHook(ElementHook* triahook,int index1,int index2,int index3){/*{{{*/
 
+	/*Create arrow of indices depending on location (0=base 1=surface)*/
+	int indices[3];
+	indices[0] = index1;
+	indices[1] = index2;
+	indices[2] = index3;
+
 	triahook->numanalyses=this->numanalyses;
 
-	int indices[2];
-	indices[0]=index1;
-	indices[1]=index2;
-
 	/*Spawn nodes hook*/
 	triahook->hnodes=new Hook*[this->numanalyses];
 	for(int i=0;i<this->numanalyses;i++){
@@ -126,13 +126,13 @@
 			triahook->hnodes[i]=NULL;
 		}
 		else{
-			triahook->hnodes[i]=this->hnodes[i]->Spawn(indices,2);
+			triahook->hnodes[i]=this->hnodes[i]->Spawn(indices,3);
 		}
 	}
 
-	/*do not spawn hmaterial. material will be taken care of by Tria*/
+	/*do not spawn hmaterial. material will be taken care of by Penta*/
 	triahook->hmaterial=NULL;
-	triahook->hvertices=(Hook*)this->hvertices->Spawn(indices,2);
+	triahook->hvertices=(Hook*)this->hvertices->Spawn(indices,3);
 	triahook->hmatpar=(Hook*)this->hmatpar->copy();
 }
 /*}}}*/
Index: ../trunk-jpl/src/c/classes/Elements/SegRef.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/SegRef.cpp	(revision 18924)
+++ ../trunk-jpl/src/c/classes/Elements/SegRef.cpp	(revision 18925)
@@ -28,86 +28,6 @@
 /*}}}*/
 
 /*Reference Element numerics*/
-void SegRef::GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss,int finiteelement){/*{{{*/
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-	_assert_(basis);
-
-	switch(finiteelement){
-		case P0Enum:
-			basis[0]=1.;
-			return;
-		case P1Enum: case P1DGEnum:
-			basis[0]=(1.-gauss->coord1)/2.;
-			basis[1]=(1.+gauss->coord1)/2.;
-			return;
-		case P2Enum:
-			basis[0]=(gauss->coord1-1.)*gauss->coord1/2.;
-			basis[1]=gauss->coord1*(1.+gauss->coord1)/2.;
-			basis[2]=(1.-gauss->coord1)*(1.+gauss->coord1);
-			return;
-		default:
-			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
-	}
-}
-/*}}}*/
-void SegRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussSeg* gauss,int finiteelement){/*{{{*/
-
-	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
-	 * actual coordinate system): */
-	IssmDouble    Jinv;
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes(finiteelement);
-
-	/*Get nodal functions derivatives in reference triangle*/
-	IssmDouble* dbasis_ref=xNew<IssmDouble>(numnodes);
-	GetNodalFunctionsDerivativesReference(dbasis_ref,gauss,finiteelement); 
-
-	/*Get Jacobian invert: */
-	GetJacobianInvert(&Jinv, xyz_list, gauss);
-
-	/*Build dbasis: 
-	 * [dhi/dx]= Jinv*[dhi/dr]
-	 */
-	for(int i=0;i<numnodes;i++){
-		dbasis[i] = Jinv*dbasis_ref[i];
-	}
-
-	/*Clean up*/
-	xDelete<IssmDouble>(dbasis_ref);
-
-}
-/*}}}*/
-void SegRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussSeg* gauss,int finiteelement){/*{{{*/
-	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
-	 * natural coordinate system) at the gaussian point. */
-
-	_assert_(dbasis && gauss);
-
-	switch(finiteelement){
-		case P0Enum:
-			/*Nodal function 1*/
-			dbasis[0] = 0.;
-			break;
-		case P1Enum: case P1DGEnum:
-			/*Nodal function 1*/
-			dbasis[0] = -0.5;
-			/*Nodal function 2*/
-			dbasis[1] = 0.5;
-			return;
-		case P2Enum:
-			/*Nodal function 1*/
-			dbasis[0] = (gauss->coord1-1.)/2. + gauss->coord1/2.;
-			dbasis[1] = (1.+gauss->coord1)/2. + gauss->coord1/2.;
-			dbasis[2] = -2.*gauss->coord1;
-			return;
-		default:
-			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
-	}
-
-}
-/*}}}*/
 void SegRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussSeg* gauss,int finiteelement){/*{{{*/
 
 	/*From node values of parameter p (plist[0],plist[1]), return parameter derivative value at gaussian 
@@ -192,6 +112,86 @@
 	*Jinv = 1./J;
 }
 /*}}}*/
+void SegRef::GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss,int finiteelement){/*{{{*/
+	/*This routine returns the values of the nodal functions  at the gaussian point.*/
+
+	_assert_(basis);
+
+	switch(finiteelement){
+		case P0Enum:
+			basis[0]=1.;
+			return;
+		case P1Enum: case P1DGEnum:
+			basis[0]=(1.-gauss->coord1)/2.;
+			basis[1]=(1.+gauss->coord1)/2.;
+			return;
+		case P2Enum:
+			basis[0]=(gauss->coord1-1.)*gauss->coord1/2.;
+			basis[1]=gauss->coord1*(1.+gauss->coord1)/2.;
+			basis[2]=(1.-gauss->coord1)*(1.+gauss->coord1);
+			return;
+		default:
+			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
+	}
+}
+/*}}}*/
+void SegRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussSeg* gauss,int finiteelement){/*{{{*/
+
+	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
+	 * actual coordinate system): */
+	IssmDouble    Jinv;
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes(finiteelement);
+
+	/*Get nodal functions derivatives in reference triangle*/
+	IssmDouble* dbasis_ref=xNew<IssmDouble>(numnodes);
+	GetNodalFunctionsDerivativesReference(dbasis_ref,gauss,finiteelement); 
+
+	/*Get Jacobian invert: */
+	GetJacobianInvert(&Jinv, xyz_list, gauss);
+
+	/*Build dbasis: 
+	 * [dhi/dx]= Jinv*[dhi/dr]
+	 */
+	for(int i=0;i<numnodes;i++){
+		dbasis[i] = Jinv*dbasis_ref[i];
+	}
+
+	/*Clean up*/
+	xDelete<IssmDouble>(dbasis_ref);
+
+}
+/*}}}*/
+void SegRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussSeg* gauss,int finiteelement){/*{{{*/
+	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
+	 * natural coordinate system) at the gaussian point. */
+
+	_assert_(dbasis && gauss);
+
+	switch(finiteelement){
+		case P0Enum:
+			/*Nodal function 1*/
+			dbasis[0] = 0.;
+			break;
+		case P1Enum: case P1DGEnum:
+			/*Nodal function 1*/
+			dbasis[0] = -0.5;
+			/*Nodal function 2*/
+			dbasis[1] = 0.5;
+			return;
+		case P2Enum:
+			/*Nodal function 1*/
+			dbasis[0] = (gauss->coord1-1.)/2. + gauss->coord1/2.;
+			dbasis[1] = (1.+gauss->coord1)/2. + gauss->coord1/2.;
+			dbasis[2] = -2.*gauss->coord1;
+			return;
+		default:
+			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
+	}
+
+}
+/*}}}*/
 int  SegRef::NumberofNodes(int finiteelement){/*{{{*/
 
 	switch(finiteelement){
Index: ../trunk-jpl/src/c/classes/Elements/TetraRef.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/TetraRef.cpp	(revision 18924)
+++ ../trunk-jpl/src/c/classes/Elements/TetraRef.cpp	(revision 18925)
@@ -30,6 +30,143 @@
 /*}}}*/
 
 /*Reference Element numerics*/
+void TetraRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussTetra* gauss,int finiteelement){/*{{{*/
+	/*From node values of parameter p (p_list[0], p_list[1], p_list[2],
+	 * p_list[3], p_list[4] and p_list[4]), return parameter derivative value at
+	 * gaussian point specified by gauss_coord:
+	 *   dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
+	 *   dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy;
+	 *   dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz;
+	 *
+	 *   p is a vector of size 3x1 already allocated.
+	 *
+	 * WARNING: For a significant gain in performance, it is better to use
+	 * static memory allocation instead of dynamic.
+	 */
+
+	/*Allocate derivatives of basis functions*/
+	IssmDouble  dbasis[3*NUMNODESMAX];
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes(finiteelement);
+	_assert_(numnodes<=NUMNODESMAX);
+
+	/*Get basis functions derivatives at this point*/
+	GetNodalFunctionsDerivatives(&dbasis[0],xyz_list,gauss,finiteelement);
+
+	/*Calculate parameter for this Gauss point*/
+	IssmDouble dpx=0.;
+	IssmDouble dpy=0.;
+	IssmDouble dpz=0.;
+	for(int i=0;i<numnodes;i++) dpx += dbasis[0*numnodes+i]*plist[i];
+	for(int i=0;i<numnodes;i++) dpy += dbasis[1*numnodes+i]*plist[i];
+	for(int i=0;i<numnodes;i++) dpz += dbasis[2*numnodes+i]*plist[i];
+
+	/*Assign values*/
+	p[0]=dpx;
+	p[1]=dpy;
+	p[2]=dpz;
+}
+/*}}}*/
+void TetraRef::GetInputValue(IssmDouble* p, IssmDouble* plist, Gauss* gauss,int finiteelement){/*{{{*/
+	/* WARNING: For a significant gain in performance, it is better to use
+	 * static memory allocation instead of dynamic.*/
+
+	/*Allocate basis functions*/
+	IssmDouble  basis[NUMNODESMAX];
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes(finiteelement);
+	_assert_(numnodes<=NUMNODESMAX);
+
+	/*Get basis functions at this point*/
+	GetNodalFunctions(&basis[0],gauss,finiteelement);
+
+	/*Calculate parameter for this Gauss point*/
+	IssmDouble value =0.;
+	for(int i=0;i<numnodes;i++) value += basis[i]*plist[i];
+
+	/*Assign output pointer*/
+	*p = value;
+}
+/*}}}*/
+void TetraRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussTetra* gauss){/*{{{*/
+	/*The Jacobian is constant over the element, discard the gaussian points. 
+	 * J is assumed to have been allocated of size 1*/
+
+	IssmDouble x1=xyz_list[3*0+0];
+	IssmDouble x2=xyz_list[3*1+0];
+	IssmDouble x3=xyz_list[3*2+0];
+	IssmDouble x4=xyz_list[3*3+0];
+
+	IssmDouble y1=xyz_list[3*0+1];
+	IssmDouble y2=xyz_list[3*1+1];
+	IssmDouble y3=xyz_list[3*2+1];
+	IssmDouble y4=xyz_list[3*3+1];
+
+	IssmDouble z1=xyz_list[3*0+2];
+	IssmDouble z2=xyz_list[3*1+2];
+	IssmDouble z3=xyz_list[3*2+2];
+	IssmDouble z4=xyz_list[3*3+2];
+
+	J[NDOF3*0+0] = x2-x1;
+	J[NDOF3*0+1] = y2-y1;
+	J[NDOF3*0+2] = z2-z1;
+
+	J[NDOF3*1+0] = x3-x1;
+	J[NDOF3*1+1] = y3-y1;
+	J[NDOF3*1+2] = z3-z1;
+
+	J[NDOF3*2+0] = x4-x1;
+	J[NDOF3*2+1] = y4-y1;
+	J[NDOF3*2+2] = z4-z1;
+}
+/*}}}*/
+void TetraRef::GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTetra* gauss){/*{{{*/
+	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
+	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
+	IssmDouble J[3][3];
+
+	/*Call Jacobian routine to get the jacobian:*/
+	GetJacobian(&J[0][0],xyz_list, gauss);
+
+	/*Get Determinant*/
+	Matrix3x3Determinant(Jdet,&J[0][0]);
+	if(*Jdet<0) _error_("negative jacobian determinant!");
+
+}
+/*}}}*/
+void TetraRef::GetJacobianDeterminantFace(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTetra* gauss){/*{{{*/
+	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
+	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
+
+	IssmDouble x1=xyz_list[3*0+0];
+	IssmDouble y1=xyz_list[3*0+1];
+	IssmDouble z1=xyz_list[3*0+2];
+	IssmDouble x2=xyz_list[3*1+0];
+	IssmDouble y2=xyz_list[3*1+1];
+	IssmDouble z2=xyz_list[3*1+2];
+	IssmDouble x3=xyz_list[3*2+0];
+	IssmDouble y3=xyz_list[3*2+1];
+	IssmDouble z3=xyz_list[3*2+2];
+
+	/*Jdet = norm( AB ^ AC ) / (2 * area of the reference triangle), with areaRef=sqrt(3) */
+	*Jdet=SQRT3/6.*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2),0.5);
+	if(*Jdet<0) _error_("negative jacobian determinant!");
+}
+/*}}}*/
+void TetraRef::GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussTetra* gauss){/*{{{*/
+
+	/*Jacobian*/
+	IssmDouble J[3][3];
+
+	/*Call Jacobian routine to get the jacobian:*/
+	GetJacobian(&J[0][0], xyz_list, gauss);
+
+	/*Invert Jacobian matrix: */
+	Matrix3x3Invert(Jinv,&J[0][0]);
+}
+/*}}}*/
 void TetraRef::GetNodalFunctions(IssmDouble* basis,Gauss* gauss_in,int finiteelement){/*{{{*/
 	/*This routine returns the values of the nodal functions  at the gaussian point.*/
 
@@ -208,143 +345,6 @@
 
 }
 /*}}}*/
-void TetraRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussTetra* gauss,int finiteelement){/*{{{*/
-	/*From node values of parameter p (p_list[0], p_list[1], p_list[2],
-	 * p_list[3], p_list[4] and p_list[4]), return parameter derivative value at
-	 * gaussian point specified by gauss_coord:
-	 *   dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
-	 *   dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy;
-	 *   dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz;
-	 *
-	 *   p is a vector of size 3x1 already allocated.
-	 *
-	 * WARNING: For a significant gain in performance, it is better to use
-	 * static memory allocation instead of dynamic.
-	 */
-
-	/*Allocate derivatives of basis functions*/
-	IssmDouble  dbasis[3*NUMNODESMAX];
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes(finiteelement);
-	_assert_(numnodes<=NUMNODESMAX);
-
-	/*Get basis functions derivatives at this point*/
-	GetNodalFunctionsDerivatives(&dbasis[0],xyz_list,gauss,finiteelement);
-
-	/*Calculate parameter for this Gauss point*/
-	IssmDouble dpx=0.;
-	IssmDouble dpy=0.;
-	IssmDouble dpz=0.;
-	for(int i=0;i<numnodes;i++) dpx += dbasis[0*numnodes+i]*plist[i];
-	for(int i=0;i<numnodes;i++) dpy += dbasis[1*numnodes+i]*plist[i];
-	for(int i=0;i<numnodes;i++) dpz += dbasis[2*numnodes+i]*plist[i];
-
-	/*Assign values*/
-	p[0]=dpx;
-	p[1]=dpy;
-	p[2]=dpz;
-}
-/*}}}*/
-void TetraRef::GetInputValue(IssmDouble* p, IssmDouble* plist, Gauss* gauss,int finiteelement){/*{{{*/
-	/* WARNING: For a significant gain in performance, it is better to use
-	 * static memory allocation instead of dynamic.*/
-
-	/*Allocate basis functions*/
-	IssmDouble  basis[NUMNODESMAX];
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes(finiteelement);
-	_assert_(numnodes<=NUMNODESMAX);
-
-	/*Get basis functions at this point*/
-	GetNodalFunctions(&basis[0],gauss,finiteelement);
-
-	/*Calculate parameter for this Gauss point*/
-	IssmDouble value =0.;
-	for(int i=0;i<numnodes;i++) value += basis[i]*plist[i];
-
-	/*Assign output pointer*/
-	*p = value;
-}
-/*}}}*/
-void TetraRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussTetra* gauss){/*{{{*/
-	/*The Jacobian is constant over the element, discard the gaussian points. 
-	 * J is assumed to have been allocated of size 1*/
-
-	IssmDouble x1=xyz_list[3*0+0];
-	IssmDouble x2=xyz_list[3*1+0];
-	IssmDouble x3=xyz_list[3*2+0];
-	IssmDouble x4=xyz_list[3*3+0];
-
-	IssmDouble y1=xyz_list[3*0+1];
-	IssmDouble y2=xyz_list[3*1+1];
-	IssmDouble y3=xyz_list[3*2+1];
-	IssmDouble y4=xyz_list[3*3+1];
-
-	IssmDouble z1=xyz_list[3*0+2];
-	IssmDouble z2=xyz_list[3*1+2];
-	IssmDouble z3=xyz_list[3*2+2];
-	IssmDouble z4=xyz_list[3*3+2];
-
-	J[NDOF3*0+0] = x2-x1;
-	J[NDOF3*0+1] = y2-y1;
-	J[NDOF3*0+2] = z2-z1;
-
-	J[NDOF3*1+0] = x3-x1;
-	J[NDOF3*1+1] = y3-y1;
-	J[NDOF3*1+2] = z3-z1;
-
-	J[NDOF3*2+0] = x4-x1;
-	J[NDOF3*2+1] = y4-y1;
-	J[NDOF3*2+2] = z4-z1;
-}
-/*}}}*/
-void TetraRef::GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTetra* gauss){/*{{{*/
-	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
-	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
-	IssmDouble J[3][3];
-
-	/*Call Jacobian routine to get the jacobian:*/
-	GetJacobian(&J[0][0],xyz_list, gauss);
-
-	/*Get Determinant*/
-	Matrix3x3Determinant(Jdet,&J[0][0]);
-	if(*Jdet<0) _error_("negative jacobian determinant!");
-
-}
-/*}}}*/
-void TetraRef::GetJacobianDeterminantFace(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTetra* gauss){/*{{{*/
-	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
-	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
-
-	IssmDouble x1=xyz_list[3*0+0];
-	IssmDouble y1=xyz_list[3*0+1];
-	IssmDouble z1=xyz_list[3*0+2];
-	IssmDouble x2=xyz_list[3*1+0];
-	IssmDouble y2=xyz_list[3*1+1];
-	IssmDouble z2=xyz_list[3*1+2];
-	IssmDouble x3=xyz_list[3*2+0];
-	IssmDouble y3=xyz_list[3*2+1];
-	IssmDouble z3=xyz_list[3*2+2];
-
-	/*Jdet = norm( AB ^ AC ) / (2 * area of the reference triangle), with areaRef=sqrt(3) */
-	*Jdet=SQRT3/6.*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2),0.5);
-	if(*Jdet<0) _error_("negative jacobian determinant!");
-}
-/*}}}*/
-void TetraRef::GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussTetra* gauss){/*{{{*/
-
-	/*Jacobian*/
-	IssmDouble J[3][3];
-
-	/*Call Jacobian routine to get the jacobian:*/
-	GetJacobian(&J[0][0], xyz_list, gauss);
-
-	/*Invert Jacobian matrix: */
-	Matrix3x3Invert(Jinv,&J[0][0]);
-}
-/*}}}*/
 int  TetraRef::NumberofNodes(int finiteelement){/*{{{*/
 
 	switch(finiteelement){
@@ -367,22 +367,6 @@
 	return -1;
 }
 /*}}}*/
-int  TetraRef::VelocityInterpolation(int fe_stokes){/*{{{*/
-
-	switch(fe_stokes){
-		case P1P1Enum:          return P1Enum;
-		case P1P1GLSEnum:       return P1Enum;
-		case MINIcondensedEnum: return P1bubbleEnum;
-		case MINIEnum:          return P1bubbleEnum;
-		case TaylorHoodEnum:    return P2Enum;
-		case LATaylorHoodEnum:  return P2Enum;
-		case XTaylorHoodEnum:   return P2Enum;
-		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
-	}
-
-	return -1;
-}
-/*}}}*/
 int  TetraRef::PressureInterpolation(int fe_stokes){/*{{{*/
 
 	switch(fe_stokes){
@@ -407,3 +391,19 @@
 	}
 }
 /*}}}*/
+int  TetraRef::VelocityInterpolation(int fe_stokes){/*{{{*/
+
+	switch(fe_stokes){
+		case P1P1Enum:          return P1Enum;
+		case P1P1GLSEnum:       return P1Enum;
+		case MINIcondensedEnum: return P1bubbleEnum;
+		case MINIEnum:          return P1bubbleEnum;
+		case TaylorHoodEnum:    return P2Enum;
+		case LATaylorHoodEnum:  return P2Enum;
+		case XTaylorHoodEnum:   return P2Enum;
+		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
+	}
+
+	return -1;
+}
+/*}}}*/
Index: ../trunk-jpl/src/c/classes/Elements/ElementHook.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/ElementHook.h	(revision 18924)
+++ ../trunk-jpl/src/c/classes/Elements/ElementHook.h	(revision 18925)
@@ -23,10 +23,10 @@
 		ElementHook(int in_numanalyses,int material_id,int numvertices,IoModel* iomodel);
 		~ElementHook();
 
+		void InitHookNeighbors(int* element_ids);               //3d only
 		void SetHookNodes(int* node_ids,int numnodes,int analysis_counter);
+		void SpawnSegHook(ElementHook* triahook,int ndex1,int index2); //2d only
 		void SpawnTriaHook(ElementHook* triahook,int index1,int index2,int index3); //3d only
-		void SpawnSegHook(ElementHook* triahook,int ndex1,int index2); //2d only
-		void InitHookNeighbors(int* element_ids);               //3d only
 };
 
 #endif //ifndef _ELEMENTHOOK_H_
Index: ../trunk-jpl/src/c/classes/Elements/SegRef.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/SegRef.h	(revision 18924)
+++ ../trunk-jpl/src/c/classes/Elements/SegRef.h	(revision 18925)
@@ -15,14 +15,14 @@
 		SegRef();
 		~SegRef();
 
+		void GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussSeg* gauss,int finiteelement);
+		void GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* gauss,int finiteelement);
 		void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussSeg* gauss);
 		void GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussSeg* gauss);
 		void GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussSeg* gauss);
 		void GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss,int finiteelement);
 		void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussSeg* gauss,int finiteelement);
 		void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussSeg* gauss,int finiteelement);
-		void GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussSeg* gauss,int finiteelement);
-		void GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* gauss,int finiteelement);
 		int  NumberofNodes(int finiteelement);
 };
 #endif
Index: ../trunk-jpl/src/c/classes/Elements/TetraRef.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/TetraRef.h	(revision 18924)
+++ ../trunk-jpl/src/c/classes/Elements/TetraRef.h	(revision 18925)
@@ -15,6 +15,8 @@
 		TetraRef();
 		~TetraRef();
 
+		void GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussTetra* gauss,int finiteelement);
+		void GetInputValue(IssmDouble* p, IssmDouble* plist, Gauss* gauss,int finiteelement);
 		void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussTetra* gauss);
 		void GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTetra* gauss);
 		void GetJacobianDeterminantFace(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTetra* gauss);
@@ -22,12 +24,9 @@
 		void GetNodalFunctions(IssmDouble* basis,Gauss* gauss_in,int finiteelement);
 		void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTetra* gauss,int finiteelement);
 		void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTetra* gauss,int finiteelement);
-		void GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussTetra* gauss,int finiteelement);
-		void GetInputValue(IssmDouble* p, IssmDouble* plist, Gauss* gauss,int finiteelement);
-
 		int  NumberofNodes(int finiteelement);
-		int  VelocityInterpolation(int fe_stokes);
 		int  PressureInterpolation(int fe_stokes);
 		int  TensorInterpolation(int fe_stokes);
+		int  VelocityInterpolation(int fe_stokes);
 };
 #endif
Index: ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 18924)
+++ ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 18925)
@@ -31,54 +31,60 @@
 /*}}}*/
 
 /*Reference Element numerics*/
-void TriaRef::GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/
-	/*Compute B  matrix. B=[phi1 phi2 -phi3 -phi4]
+void TriaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement){/*{{{*/
+	/*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter derivative value at gaussian 
+	 * point specified by gauss_basis:
+	 *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
+	 *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
 	 *
-	 * and phi1=phi3 phi2=phi4
+	 * p is a vector already allocated.
 	 *
-	 * We assume B has been allocated already, of size: 1x4
+	 * WARNING: For a significant gain in performance, it is better to use
+	 * static memory allocation instead of dynamic.
 	 */
 
+	/*Allocate derivatives of basis functions*/
+	IssmDouble  dbasis[2*NUMNODESMAX];
+
 	/*Fetch number of nodes for this finite element*/
 	int numnodes = this->NumberofNodes(finiteelement);
+	_assert_(numnodes<=NUMNODESMAX);
 
-	/*Get nodal functions*/
-	IssmDouble* basis=xNew<IssmDouble>(numnodes);
-	GetNodalFunctions(basis,gauss,finiteelement);
+	/*Get basis functions derivatives at this point*/
+	GetNodalFunctionsDerivatives(&dbasis[0],xyz_list,gauss,finiteelement);
 
-	/*Build B for this segment*/
-	B[0] = +basis[index1];
-	B[1] = +basis[index2];
-	B[2] = -basis[index1];
-	B[3] = -basis[index2];
+	/*Calculate parameter for this Gauss point*/
+	IssmDouble dpx=0.;
+	IssmDouble dpy=0.;
+	for(int i=0;i<numnodes;i++) dpx += dbasis[0*numnodes+i]*plist[i];
+	for(int i=0;i<numnodes;i++) dpy += dbasis[1*numnodes+i]*plist[i];
 
-	/*Clean-up*/
-	xDelete<IssmDouble>(basis);
+	/*Assign values*/
+	p[0]=dpx;
+	p[1]=dpy;
+
 }
 /*}}}*/
-void TriaRef::GetSegmentBprimeFlux(IssmDouble* Bprime,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/
-	/*Compute Bprime  matrix. Bprime=[phi1 phi2 phi3 phi4]
-	 *
-	 * and phi1=phi3 phi2=phi4
-	 *
-	 * We assume Bprime has been allocated already, of size: 1x4
-	 */
+void TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, Gauss* gauss,int finiteelement){/*{{{*/
+	/* WARNING: For a significant gain in performance, it is better to use
+	 * static memory allocation instead of dynamic.*/
 
+	/*Allocate basis functions*/
+	IssmDouble  basis[NUMNODESMAX];
+
 	/*Fetch number of nodes for this finite element*/
 	int numnodes = this->NumberofNodes(finiteelement);
+	_assert_(numnodes<=NUMNODESMAX);
 
-	/*Get nodal functions*/
-	IssmDouble* basis=xNew<IssmDouble>(numnodes);
-	GetNodalFunctions(basis,gauss,finiteelement);
+	/*Get basis functions at this point*/
+	GetNodalFunctions(&basis[0],gauss,finiteelement);
 
-	/*Build B'*/
-	Bprime[0] = basis[index1];
-	Bprime[1] = basis[index2];
-	Bprime[2] = basis[index1];
-	Bprime[3] = basis[index2];
+	/*Calculate parameter for this Gauss point*/
+	IssmDouble value =0.;
+	for(int i=0;i<numnodes;i++) value += basis[i]*plist[i];
 
-	/*Clean-up*/
-	xDelete<IssmDouble>(basis);
+	/*Assign output pointer*/
+	*p = value;
 }
 /*}}}*/
 void TriaRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
@@ -98,20 +104,6 @@
 	J[2*1+1] = SQRT3/6.0*(2*y3-y1-y2);
 }
 /*}}}*/
-void TriaRef::GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
-	 * J is assumed to have been allocated*/
-
-	IssmDouble x1 = xyz_list[3*0+0];
-	IssmDouble y1 = xyz_list[3*0+1];
-	IssmDouble x2 = xyz_list[3*1+0];
-	IssmDouble y2 = xyz_list[3*1+1];
-
-	*Jdet = .5*sqrt(pow(x2-x1,2) + pow(y2-y1,2));
-	if(*Jdet<0) _error_("negative jacobian determinant!");
-
-}
-/*}}}*/
 void TriaRef::GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
 	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
@@ -194,45 +186,6 @@
 	}
 }
 /*}}}*/
-void TriaRef::GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss,int index1,int index2,int finiteelement){/*{{{*/
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-	_assert_(index1>=0 && index1<3);
-	_assert_(index2>=0 && index2<3);
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes(finiteelement);
-
-	/*Get nodal functions*/
-	IssmDouble* triabasis=xNew<IssmDouble>(numnodes);
-	GetNodalFunctions(triabasis,gauss,finiteelement);
-
-	switch(finiteelement){
-		case P1Enum: case P1DGEnum:
-			basis[0]=triabasis[index1];
-			basis[1]=triabasis[index2];
-			xDelete<IssmDouble>(triabasis);
-			return;
-		case P1bubbleEnum: case P1bubblecondensedEnum:
-			basis[0]=triabasis[index1];
-			basis[1]=triabasis[index2];
-			xDelete<IssmDouble>(triabasis);
-			return;
-		case P2Enum:
-			_assert_(index2<index1);
-			basis[0]=triabasis[index1];
-			basis[1]=triabasis[index2];
-			basis[2]=triabasis[3+index2-1];
-			xDelete<IssmDouble>(triabasis);
-			return;
-		default:
-			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
-	}
-
-	/*Clean up*/
-	xDelete<IssmDouble>(triabasis);
-}
-/*}}}*/
 void TriaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss,int finiteelement){/*{{{*/
 
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
@@ -353,134 +306,109 @@
 
 }
 /*}}}*/
-void TriaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement){/*{{{*/
-	/*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter derivative value at gaussian 
-	 * point specified by gauss_basis:
-	 *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
-	 *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
+void TriaRef::GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/
+	/*Compute B  matrix. B=[phi1 phi2 -phi3 -phi4]
 	 *
-	 * p is a vector already allocated.
+	 * and phi1=phi3 phi2=phi4
 	 *
-	 * WARNING: For a significant gain in performance, it is better to use
-	 * static memory allocation instead of dynamic.
+	 * We assume B has been allocated already, of size: 1x4
 	 */
 
-	/*Allocate derivatives of basis functions*/
-	IssmDouble  dbasis[2*NUMNODESMAX];
-
 	/*Fetch number of nodes for this finite element*/
 	int numnodes = this->NumberofNodes(finiteelement);
-	_assert_(numnodes<=NUMNODESMAX);
 
-	/*Get basis functions derivatives at this point*/
-	GetNodalFunctionsDerivatives(&dbasis[0],xyz_list,gauss,finiteelement);
+	/*Get nodal functions*/
+	IssmDouble* basis=xNew<IssmDouble>(numnodes);
+	GetNodalFunctions(basis,gauss,finiteelement);
 
-	/*Calculate parameter for this Gauss point*/
-	IssmDouble dpx=0.;
-	IssmDouble dpy=0.;
-	for(int i=0;i<numnodes;i++) dpx += dbasis[0*numnodes+i]*plist[i];
-	for(int i=0;i<numnodes;i++) dpy += dbasis[1*numnodes+i]*plist[i];
+	/*Build B for this segment*/
+	B[0] = +basis[index1];
+	B[1] = +basis[index2];
+	B[2] = -basis[index1];
+	B[3] = -basis[index2];
 
-	/*Assign values*/
-	p[0]=dpx;
-	p[1]=dpy;
-
+	/*Clean-up*/
+	xDelete<IssmDouble>(basis);
 }
 /*}}}*/
-void TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, Gauss* gauss,int finiteelement){/*{{{*/
-	/* WARNING: For a significant gain in performance, it is better to use
-	 * static memory allocation instead of dynamic.*/
+void TriaRef::GetSegmentBprimeFlux(IssmDouble* Bprime,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/
+	/*Compute Bprime  matrix. Bprime=[phi1 phi2 phi3 phi4]
+	 *
+	 * and phi1=phi3 phi2=phi4
+	 *
+	 * We assume Bprime has been allocated already, of size: 1x4
+	 */
 
-	/*Allocate basis functions*/
-	IssmDouble  basis[NUMNODESMAX];
-
 	/*Fetch number of nodes for this finite element*/
 	int numnodes = this->NumberofNodes(finiteelement);
-	_assert_(numnodes<=NUMNODESMAX);
 
-	/*Get basis functions at this point*/
-	GetNodalFunctions(&basis[0],gauss,finiteelement);
+	/*Get nodal functions*/
+	IssmDouble* basis=xNew<IssmDouble>(numnodes);
+	GetNodalFunctions(basis,gauss,finiteelement);
 
-	/*Calculate parameter for this Gauss point*/
-	IssmDouble value =0.;
-	for(int i=0;i<numnodes;i++) value += basis[i]*plist[i];
+	/*Build B'*/
+	Bprime[0] = basis[index1];
+	Bprime[1] = basis[index2];
+	Bprime[2] = basis[index1];
+	Bprime[3] = basis[index2];
 
-	/*Assign output pointer*/
-	*p = value;
+	/*Clean-up*/
+	xDelete<IssmDouble>(basis);
 }
 /*}}}*/
-int  TriaRef::NumberofNodes(int finiteelement){/*{{{*/
+void TriaRef::GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
+	 * J is assumed to have been allocated*/
 
-	switch(finiteelement){
-		case NoneEnum:                return 0;
-		case P0Enum:                  return NUMNODESP0;
-		case P1Enum:                  return NUMNODESP1;
-		case P1DGEnum:                return NUMNODESP1;
-		case P1bubbleEnum:            return NUMNODESP1b;
-		case P1bubblecondensedEnum:   return NUMNODESP1b;
-		case P2Enum:                  return NUMNODESP2;
-		case P2bubbleEnum:            return NUMNODESP2b;
-		case P2bubblecondensedEnum:   return NUMNODESP2b;
-		case P1P1Enum:                return NUMNODESP1*2;
-		case P1P1GLSEnum:             return NUMNODESP1*2;
-		case MINIcondensedEnum:       return NUMNODESP1b+NUMNODESP1;
-		case MINIEnum:                return NUMNODESP1b+NUMNODESP1;
-		case TaylorHoodEnum:          return NUMNODESP2+NUMNODESP1;
-		case LATaylorHoodEnum:        return NUMNODESP2;
-		case XTaylorHoodEnum:         return NUMNODESP2+NUMNODESP1;
-		case CrouzeixRaviartEnum:     return NUMNODESP2b+NUMNODESP1;
-		case LACrouzeixRaviartEnum:   return NUMNODESP2b;
-		default: _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
-	}
+	IssmDouble x1 = xyz_list[3*0+0];
+	IssmDouble y1 = xyz_list[3*0+1];
+	IssmDouble x2 = xyz_list[3*1+0];
+	IssmDouble y2 = xyz_list[3*1+1];
 
-	return -1;
+	*Jdet = .5*sqrt(pow(x2-x1,2) + pow(y2-y1,2));
+	if(*Jdet<0) _error_("negative jacobian determinant!");
+
 }
 /*}}}*/
-int  TriaRef::VelocityInterpolation(int fe_stokes){/*{{{*/
+void TriaRef::GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss,int index1,int index2,int finiteelement){/*{{{*/
+	/*This routine returns the values of the nodal functions  at the gaussian point.*/
 
-	switch(fe_stokes){
-		case P1P1Enum:              return P1Enum;
-		case P1P1GLSEnum:           return P1Enum;
-		case MINIcondensedEnum:     return P1bubbleEnum;
-		case MINIEnum:              return P1bubbleEnum;
-		case TaylorHoodEnum:        return P2Enum;
-		case LATaylorHoodEnum:      return P2Enum;
-		case XTaylorHoodEnum:       return P2Enum;
-		case CrouzeixRaviartEnum:   return P2bubbleEnum;
-		case LACrouzeixRaviartEnum: return P2bubbleEnum;
-		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
-	}
+	_assert_(index1>=0 && index1<3);
+	_assert_(index2>=0 && index2<3);
 
-	return -1;
-}
-/*}}}*/
-int  TriaRef::PressureInterpolation(int fe_stokes){/*{{{*/
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes(finiteelement);
 
-	switch(fe_stokes){
-		case P1P1Enum:              return P1Enum;
-		case P1P1GLSEnum:           return P1Enum;
-		case MINIcondensedEnum:     return P1Enum;
-		case MINIEnum:              return P1Enum;
-		case TaylorHoodEnum:        return P1Enum;
-		case LATaylorHoodEnum:      return NoneEnum;
-		case XTaylorHoodEnum:       return P1Enum;
-		case CrouzeixRaviartEnum:   return P1DGEnum;
-		case LACrouzeixRaviartEnum: return NoneEnum;
-		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
+	/*Get nodal functions*/
+	IssmDouble* triabasis=xNew<IssmDouble>(numnodes);
+	GetNodalFunctions(triabasis,gauss,finiteelement);
+
+	switch(finiteelement){
+		case P1Enum: case P1DGEnum:
+			basis[0]=triabasis[index1];
+			basis[1]=triabasis[index2];
+			xDelete<IssmDouble>(triabasis);
+			return;
+		case P1bubbleEnum: case P1bubblecondensedEnum:
+			basis[0]=triabasis[index1];
+			basis[1]=triabasis[index2];
+			xDelete<IssmDouble>(triabasis);
+			return;
+		case P2Enum:
+			_assert_(index2<index1);
+			basis[0]=triabasis[index1];
+			basis[1]=triabasis[index2];
+			basis[2]=triabasis[3+index2-1];
+			xDelete<IssmDouble>(triabasis);
+			return;
+		default:
+			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
 	}
 
-	return -1;
+	/*Clean up*/
+	xDelete<IssmDouble>(triabasis);
 }
 /*}}}*/
-int  TriaRef::TensorInterpolation(int fe_stokes){/*{{{*/
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-	switch(fe_stokes){
-		case XTaylorHoodEnum: return P1DGEnum;
-		default: _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
-	}
-}
-/*}}}*/
 void TriaRef::NodeOnEdgeIndices(int* pnumindices,int** pindices,int index,int finiteelement){/*{{{*/
 
 	/*Output*/
@@ -540,3 +468,75 @@
 	*pindices    = indices;
 }
 /*}}}*/
+int  TriaRef::NumberofNodes(int finiteelement){/*{{{*/
+
+	switch(finiteelement){
+		case NoneEnum:                return 0;
+		case P0Enum:                  return NUMNODESP0;
+		case P1Enum:                  return NUMNODESP1;
+		case P1DGEnum:                return NUMNODESP1;
+		case P1bubbleEnum:            return NUMNODESP1b;
+		case P1bubblecondensedEnum:   return NUMNODESP1b;
+		case P2Enum:                  return NUMNODESP2;
+		case P2bubbleEnum:            return NUMNODESP2b;
+		case P2bubblecondensedEnum:   return NUMNODESP2b;
+		case P1P1Enum:                return NUMNODESP1*2;
+		case P1P1GLSEnum:             return NUMNODESP1*2;
+		case MINIcondensedEnum:       return NUMNODESP1b+NUMNODESP1;
+		case MINIEnum:                return NUMNODESP1b+NUMNODESP1;
+		case TaylorHoodEnum:          return NUMNODESP2+NUMNODESP1;
+		case LATaylorHoodEnum:        return NUMNODESP2;
+		case XTaylorHoodEnum:         return NUMNODESP2+NUMNODESP1;
+		case CrouzeixRaviartEnum:     return NUMNODESP2b+NUMNODESP1;
+		case LACrouzeixRaviartEnum:   return NUMNODESP2b;
+		default: _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
+	}
+
+	return -1;
+}
+/*}}}*/
+int  TriaRef::PressureInterpolation(int fe_stokes){/*{{{*/
+
+	switch(fe_stokes){
+		case P1P1Enum:              return P1Enum;
+		case P1P1GLSEnum:           return P1Enum;
+		case MINIcondensedEnum:     return P1Enum;
+		case MINIEnum:              return P1Enum;
+		case TaylorHoodEnum:        return P1Enum;
+		case LATaylorHoodEnum:      return NoneEnum;
+		case XTaylorHoodEnum:       return P1Enum;
+		case CrouzeixRaviartEnum:   return P1DGEnum;
+		case LACrouzeixRaviartEnum: return NoneEnum;
+		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
+	}
+
+	return -1;
+}
+/*}}}*/
+int  TriaRef::TensorInterpolation(int fe_stokes){/*{{{*/
+	/*This routine returns the values of the nodal functions  at the gaussian point.*/
+
+	switch(fe_stokes){
+		case XTaylorHoodEnum: return P1DGEnum;
+		default: _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
+	}
+}
+/*}}}*/
+int  TriaRef::VelocityInterpolation(int fe_stokes){/*{{{*/
+
+	switch(fe_stokes){
+		case P1P1Enum:              return P1Enum;
+		case P1P1GLSEnum:           return P1Enum;
+		case MINIcondensedEnum:     return P1bubbleEnum;
+		case MINIEnum:              return P1bubbleEnum;
+		case TaylorHoodEnum:        return P2Enum;
+		case LATaylorHoodEnum:      return P2Enum;
+		case XTaylorHoodEnum:       return P2Enum;
+		case CrouzeixRaviartEnum:   return P2bubbleEnum;
+		case LACrouzeixRaviartEnum: return P2bubbleEnum;
+		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
+	}
+
+	return -1;
+}
+/*}}}*/
