Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 17873)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 17874)
@@ -242,4 +242,5 @@
 		virtual int    VelocityInterpolation()=0;
 		virtual int    PressureInterpolation()=0;
+		virtual int    TensorInterpolation()=0;
 		virtual bool   IsZeroLevelset(int levelset_enum)=0;
 		virtual bool   IsIcefront(void)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17873)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17874)
@@ -467,6 +467,7 @@
 /*FUNCTION Penta::Delta18oParameterization{{{*/
 void  Penta::Delta18oParameterization(void){
-        /*Are we on the base? If not, return*/
-        if(!IsOnBase()) return;
+
+	/*Are we on the base? If not, return*/
+	if(!IsOnBase()) return;
 
 	int        i;
@@ -1850,5 +1851,5 @@
 
 	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctionsDerivativesVelocity(dbasis,xyz_list,(GaussPenta*)gauss);
+	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,this->VelocityInterpolation());
 
 }
@@ -1858,5 +1859,5 @@
 
 	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctionsVelocity(basis,(GaussPenta*)gauss);
+	this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->VelocityInterpolation());
 
 }
@@ -1866,5 +1867,5 @@
 
 	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctionsPressure(basis,(GaussPenta*)gauss);
+	this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->PressureInterpolation());
 
 }
@@ -1874,5 +1875,5 @@
 
 	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctionsTensor(basis,(GaussPenta*)gauss);
+	this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->TensorInterpolation());
 
 }
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 17873)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 17874)
@@ -92,4 +92,5 @@
 		int    VelocityInterpolation();
 		int    PressureInterpolation();
+		int    TensorInterpolation(){_error_("not implemented yet");};
 		bool   IsZeroLevelset(int levelset_enum);
 		bool   IsIcefront(void);
Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 17873)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 17874)
@@ -328,43 +328,15 @@
 			return;
 		default:
-			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
 	}
 }
 /*}}}*/
-/*FUNCTION PentaRef::GetNodalFunctionsVelocity{{{*/
-void PentaRef::GetNodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss){
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-	int save = this->element_type;
-	this->element_type = VelocityInterpolation();
-	this->GetNodalFunctions(basis,gauss);
-	this->element_type = save;
-
-}
-/*}}}*/
-/*FUNCTION PentaRef::GetNodalFunctionsPressure{{{*/
-void PentaRef::GetNodalFunctionsPressure(IssmDouble* basis,Gauss* gauss){
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-	int save = this->element_type;
-	this->element_type = PressureInterpolation();
-	this->GetNodalFunctions(basis,gauss);
-	this->element_type = save;
-
-}
-/*}}}*/
-/*FUNCTION PentaRef::GetNodalFunctionsTensor{{{*/
-void PentaRef::GetNodalFunctionsTensor(IssmDouble* basis,Gauss* gauss){
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-	int save = this->element_type;
-	this->element_type = TensorInterpolation();
-	this->GetNodalFunctions(basis,gauss);
-	this->element_type = save;
-
-}
-/*}}}*/
-/*FUNCTION PentaRef::GetNodalFunctionsDerivatives{{{*/
+/*FUNCTION PentaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss){{{*/
 void PentaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss){
+	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss,this->element_type);
+}
+/*}}}*/
+/*FUNCTION PentaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss,int finiteelement){{{*/
+void PentaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss,int finiteelement){
 
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
@@ -373,9 +345,9 @@
 
 	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(finiteelement);
 
 	/*Get nodal functions derivatives in reference triangle*/
 	IssmDouble* dbasis_ref=xNew<IssmDouble>(3*numnodes);
-	GetNodalFunctionsDerivativesReference(dbasis_ref,gauss); 
+	GetNodalFunctionsDerivativesReference(dbasis_ref,gauss,finiteelement);
 
 	/*Get Jacobian invert: */
@@ -397,28 +369,13 @@
 	/*Clean up*/
 	xDelete<IssmDouble>(dbasis_ref);
-
-}
-/*}}}*/
-/*FUNCTION PentaRef::GetNodalFunctionsDerivativesVelocity{{{*/
-void PentaRef::GetNodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss){
-
-	int save = this->element_type;
-	this->element_type = VelocityInterpolation();
-	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-	this->element_type = save;
-
-}
-/*}}}*/
-/*FUNCTION PentaRef::GetNodalFunctionsDerivativesPressure{{{*/
-void PentaRef::GetNodalFunctionsDerivativesPressure(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss){
-
-	int save = this->element_type;
-	this->element_type = PressureInterpolation();
-	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-	this->element_type = save;
-}
-/*}}}*/
-/*FUNCTION PentaRef::GetNodalFunctionsDerivativesReference{{{*/
-void PentaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss_in){
+}
+/*}}}*/
+/*FUNCTION PentaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss){{{*/
+void PentaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss){
+	GetNodalFunctionsDerivativesReference(dbasis,gauss,this->element_type);
+}
+/*}}}*/
+/*FUNCTION PentaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss_in,int finiteelement){{{*/
+void PentaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss_in,int finiteelement){
 
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
@@ -434,5 +391,5 @@
 	IssmDouble zeta=gauss->coord4;
 
-	switch(this->element_type){
+	switch(finiteelement){
 		case P0Enum: 
 			/*Zero derivative*/
@@ -837,5 +794,5 @@
 			return;
 		default:
-			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
 	}
 
Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 17873)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 17874)
@@ -24,11 +24,8 @@
 		void GetNodalFunctions(IssmDouble* basis, Gauss* gauss);
 		void GetNodalFunctions(IssmDouble* basis, Gauss* gauss,int finiteelement);
-		void GetNodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss);
-		void GetNodalFunctionsPressure(IssmDouble* basis, Gauss* gauss);
-		void GetNodalFunctionsTensor(IssmDouble* basis, Gauss* gauss);
 		void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
-		void GetNodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
-		void GetNodalFunctionsDerivativesPressure(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
+		void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss,int finiteelement);
 		void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss);
+		void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss,int finiteelement);
 		void GetNodalFunctionsP1(IssmDouble* l1l6, Gauss* gauss);
 		void GetNodalFunctionsMINI(IssmDouble* l1l7, Gauss* gauss);
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 17873)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 17874)
@@ -105,4 +105,5 @@
 		void        ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
 		int         VelocityInterpolation(void){_error_("not implemented yet");};
+		int         TensorInterpolation(void){_error_("not implemented yet");};
 		void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating){_error_("not implemented yet");};
 		IssmDouble  GetGroundedPortion(IssmDouble* xyz_list){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/SegRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/SegRef.cpp	(revision 17873)
+++ /issm/trunk-jpl/src/c/classes/Elements/SegRef.cpp	(revision 17874)
@@ -91,5 +91,5 @@
 	_assert_(basis);
 
-	switch(element_type){
+	switch(finiteelement){
 		case P0Enum:
 			basis[0]=1.;
@@ -105,5 +105,5 @@
 			return;
 		default:
-			_error_("Element type "<<EnumToStringx(element_type)<<" not supported yet");
+			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
 	}
 }
@@ -179,5 +179,5 @@
 			return;
 		default:
-			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
 	}
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp	(revision 17873)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp	(revision 17874)
@@ -590,5 +590,5 @@
 
 	_assert_(gauss->Enum()==GaussTetraEnum);
-	this->GetNodalFunctionsVelocity(basis,(GaussTetra*)gauss);
+	this->GetNodalFunctions(basis,(GaussTetra*)gauss,this->VelocityInterpolation());
 
 }
@@ -598,5 +598,5 @@
 
 	_assert_(gauss->Enum()==GaussTetraEnum);
-	this->GetNodalFunctionsPressure(basis,(GaussTetra*)gauss);
+	this->GetNodalFunctions(basis,(GaussTetra*)gauss,this->PressureInterpolation());
 
 }
@@ -606,5 +606,5 @@
 
 	_assert_(gauss->Enum()==GaussTetraEnum);
-	this->GetNodalFunctionsTensor(basis,(GaussTetra*)gauss);
+	this->GetNodalFunctions(basis,(GaussTetra*)gauss,this->TensorInterpolation());
 
 }
@@ -622,5 +622,5 @@
 
 	_assert_(gauss->Enum()==GaussTetraEnum);
-	this->GetNodalFunctionsDerivativesVelocity(dbasis,xyz_list,(GaussTetra*)gauss);
+	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTetra*)gauss,this->VelocityInterpolation());
 
 }
@@ -992,4 +992,9 @@
 }
 /*}}}*/
+/*FUNCTION Tetra::TensorInterpolation{{{*/
+int Tetra::TensorInterpolation(void){
+	return TetraRef::TensorInterpolation();
+}
+/*}}}*/
 /*FUNCTION Tetra::ZeroLevelsetCoordinates{{{*/
 void Tetra::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){
Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 17873)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 17874)
@@ -110,4 +110,5 @@
 		void        ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
 		int         VelocityInterpolation(void);
+		int         TensorInterpolation(void);
 		void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating){_error_("not implemented yet");};
 		IssmDouble  GetGroundedPortion(IssmDouble* xyz_list){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp	(revision 17873)
+++ /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp	(revision 17874)
@@ -98,90 +98,5 @@
 			return;
 		default:
-			_error_("Element type "<<EnumToStringx(element_type)<<" not supported yet");
-	}
-}
-/*}}}*/
-/*FUNCTION TetraRef::GetNodalFunctionsVelocity{{{*/
-void TetraRef::GetNodalFunctionsVelocity(IssmDouble* basis,GaussTetra* gauss){
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-	switch(this->element_type){
-		case P1P1Enum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctions(basis,gauss);
-			this->element_type = P1P1Enum;
-			return;
-		case P1P1GLSEnum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctions(basis,gauss);
-			this->element_type = P1P1GLSEnum;
-			return;
-		case MINIcondensedEnum:
-			this->element_type = P1bubbleEnum;
-			this->GetNodalFunctions(basis,gauss);
-			this->element_type = MINIcondensedEnum;
-			return;
-		case MINIEnum:
-			this->element_type = P1bubbleEnum;
-			this->GetNodalFunctions(basis,gauss);
-			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");
-	}
-}
-/*}}}*/
-/*FUNCTION TetraRef::GetNodalFunctionsPressure{{{*/
-void TetraRef::GetNodalFunctionsPressure(IssmDouble* basis,GaussTetra* gauss){
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-	switch(this->element_type){
-		case P1P1Enum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctions(basis,gauss);
-			this->element_type = P1P1Enum;
-			return;
-		case P1P1GLSEnum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctions(basis,gauss);
-			this->element_type = P1P1GLSEnum;
-			return;
-		case MINIcondensedEnum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctions(basis,gauss);
-			this->element_type = MINIcondensedEnum;
-			return;
-		case MINIEnum:
-			this->element_type = P1Enum;
-			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:
-			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
-	}
-}
-/*}}}*/
-/*FUNCTION TetraRef::GetNodalFunctionsTensor{{{*/
-void TetraRef::GetNodalFunctionsTensor(IssmDouble* basis,GaussTetra* gauss){
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-	switch(this->element_type){
-		case XTaylorHoodEnum:
-			this->element_type = P1DGEnum;
-			this->GetNodalFunctions(basis,gauss);
-			this->element_type = XTaylorHoodEnum;
-			return;
-		default:
-			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
 	}
 }
@@ -226,70 +141,4 @@
 }
 /*}}}*/
-/*FUNCTION TetraRef::GetNodalFunctionsDerivativesPressure         (THIS ONE and following ones){{{*/
-void TetraRef::GetNodalFunctionsDerivativesPressure(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTetra* gauss){
-	switch(this->element_type){
-		case P1P1Enum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-			this->element_type = P1P1Enum;
-			return;
-		case P1P1GLSEnum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-			this->element_type = P1P1GLSEnum;
-			return;
-		case MINIcondensedEnum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-			this->element_type = MINIcondensedEnum;
-			return;
-		case MINIEnum:
-			this->element_type = P1Enum;
-			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:
-			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
-	}
-}
-/*}}}*/
-/*FUNCTION TetraRef::GetNodalFunctionsDerivativesVelocity{{{*/
-void TetraRef::GetNodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTetra* gauss){
-	switch(this->element_type){
-		case P1P1Enum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-			this->element_type = P1P1Enum;
-			return;
-		case P1P1GLSEnum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-			this->element_type = P1P1GLSEnum;
-			return;
-		case MINIcondensedEnum:
-			this->element_type = P1bubbleEnum;
-			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-			this->element_type = MINIcondensedEnum;
-			return;
-		case MINIEnum:
-			this->element_type = P1bubbleEnum;
-			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-			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");
-	}
-}
-/*}}}*/
 /*FUNCTION TetraRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTetra* gauss){{{*/
 void TetraRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTetra* gauss){
@@ -395,5 +244,5 @@
 			return;
 		default:
-			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
 	}
 
@@ -629,4 +478,13 @@
 
 	return -1;
-}
-/*}}}*/
+}/*}}}*/
+/*FUNCTION TetraRef::TensorInterpolation{{{*/
+int TetraRef::TensorInterpolation(void){
+	/*This routine returns the values of the nodal functions  at the gaussian point.*/
+
+	switch(this->element_type){
+		case XTaylorHoodEnum: return P1DGEnum;
+		default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+	}
+}
+/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/TetraRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TetraRef.h	(revision 17873)
+++ /issm/trunk-jpl/src/c/classes/Elements/TetraRef.h	(revision 17874)
@@ -29,11 +29,6 @@
 		void GetNodalFunctions(IssmDouble* basis,GaussTetra* gauss);
 		void GetNodalFunctions(IssmDouble* basis,GaussTetra* gauss,int finiteelement);
-		void GetNodalFunctionsVelocity(IssmDouble* basis, GaussTetra* gauss);
-		void GetNodalFunctionsPressure(IssmDouble* basis, GaussTetra* gauss);
-		void GetNodalFunctionsTensor(IssmDouble* basis, GaussTetra* gauss);
 		void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTetra* gauss);
 		void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTetra* gauss,int finiteelement);
-		void GetNodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,GaussTetra* gauss);
-		void GetNodalFunctionsDerivativesPressure(IssmDouble* dbasis,IssmDouble* xyz_list,GaussTetra* gauss);
 		void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTetra* gauss);
 		void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTetra* gauss,int finiteelement);
@@ -48,4 +43,5 @@
 		int  VelocityInterpolation(void);
 		int  PressureInterpolation(void);
+		int  TensorInterpolation(void);
 };
 #endif
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 17873)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 17874)
@@ -1597,5 +1597,5 @@
 
 	_assert_(gauss->Enum()==GaussTriaEnum);
-	this->GetNodalFunctionsDerivativesVelocity(dbasis,xyz_list,(GaussTria*)gauss);
+	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->VelocityInterpolation());
 
 }
@@ -1605,5 +1605,5 @@
 
 	_assert_(gauss->Enum()==GaussTriaEnum);
-	this->GetNodalFunctionsVelocity(basis,(GaussTria*)gauss);
+	this->GetNodalFunctions(basis,(GaussTria*)gauss,this->VelocityInterpolation());
 
 }
@@ -1613,5 +1613,5 @@
 
 	_assert_(gauss->Enum()==GaussTriaEnum);
-	this->GetNodalFunctionsPressure(basis,(GaussTria*)gauss);
+	this->GetNodalFunctions(basis,(GaussTria*)gauss,this->PressureInterpolation());
 
 }
@@ -1621,5 +1621,5 @@
 
 	_assert_(gauss->Enum()==GaussTriaEnum);
-	this->GetNodalFunctionsTensor(basis,(GaussTria*)gauss);
+	this->GetNodalFunctions(basis,(GaussTria*)gauss,this->TensorInterpolation());
 
 }
@@ -1702,4 +1702,9 @@
 int Tria::PressureInterpolation(void){
 	return TriaRef::PressureInterpolation();
+}
+/*}}}*/
+/*FUNCTION Tria::TensorInterpolation{{{*/
+int Tria::TensorInterpolation(void){
+	return TriaRef::TensorInterpolation();
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 17873)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 17874)
@@ -100,4 +100,5 @@
 		int         VelocityInterpolation();
 		int         PressureInterpolation();
+		int         TensorInterpolation();
 		IssmDouble  SurfaceArea(void);
 		void        Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 17873)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 17874)
@@ -212,100 +212,5 @@
 			return;
 		default:
-			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
-	}
-}
-/*}}}*/
-/*FUNCTION TriaRef::GetNodalFunctionsVelocity{{{*/
-void TriaRef::GetNodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss){
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-	switch(this->element_type){
-		case P1P1Enum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctions(basis,gauss);
-			this->element_type = P1P1Enum;
-			return;
-		case P1P1GLSEnum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctions(basis,gauss);
-			this->element_type = P1P1GLSEnum;
-			return;
-		case MINIcondensedEnum:
-			this->element_type = P1bubbleEnum;
-			this->GetNodalFunctions(basis,gauss);
-			this->element_type = MINIcondensedEnum;
-			return;
-		case MINIEnum:
-			this->element_type = P1bubbleEnum;
-			this->GetNodalFunctions(basis,gauss);
-			this->element_type = MINIEnum;
-			return;
-		case TaylorHoodEnum:
-			this->element_type = P2Enum;
-			this->GetNodalFunctions(basis,gauss);
-			this->element_type = TaylorHoodEnum;
-			return;
-		case XTaylorHoodEnum:
-			this->element_type = P2Enum;
-			this->GetNodalFunctions(basis,gauss);
-			this->element_type = XTaylorHoodEnum;
-			return;
-		default:
-			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
-	}
-}
-/*}}}*/
-/*FUNCTION TriaRef::GetNodalFunctionsPressure{{{*/
-void TriaRef::GetNodalFunctionsPressure(IssmDouble* basis,Gauss* gauss){
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-	switch(this->element_type){
-		case P1P1Enum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctions(basis,gauss);
-			this->element_type = P1P1Enum;
-			return;
-		case P1P1GLSEnum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctions(basis,gauss);
-			this->element_type = P1P1GLSEnum;
-			return;
-		case MINIcondensedEnum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctions(basis,gauss);
-			this->element_type = MINIcondensedEnum;
-			return;
-		case MINIEnum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctions(basis,gauss);
-			this->element_type = MINIEnum;
-			return;
-		case TaylorHoodEnum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctions(basis,gauss);
-			this->element_type = TaylorHoodEnum;
-			return;
-		case XTaylorHoodEnum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctions(basis,gauss);
-			this->element_type = XTaylorHoodEnum;
-			return;
-		default:
-			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
-	}
-}
-/*}}}*/
-/*FUNCTION TriaRef::GetNodalFunctionsTensor{{{*/
-void TriaRef::GetNodalFunctionsTensor(IssmDouble* basis,Gauss* gauss){
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-	switch(this->element_type){
-		case XTaylorHoodEnum:
-			this->element_type = P1DGEnum;
-			this->GetNodalFunctions(basis,gauss);
-			this->element_type = XTaylorHoodEnum;
-			return;
-		default:
-			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
 	}
 }
@@ -387,71 +292,4 @@
 	xDelete<IssmDouble>(dbasis_ref);
 
-}
-/*}}}*/
-/*FUNCTION TriaRef::GetNodalFunctionsDerivativesPressure{{{*/
-void TriaRef::GetNodalFunctionsDerivativesPressure(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss){
-	switch(this->element_type){
-		case P1P1Enum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-			this->element_type = P1P1Enum;
-			return;
-		case P1P1GLSEnum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-			this->element_type = P1P1GLSEnum;
-			return;
-		case MINIcondensedEnum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-			this->element_type = MINIcondensedEnum;
-			return;
-		case MINIEnum:
-			this->element_type = P1Enum;
-			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:
-			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
-	}
-}
-/*}}}*/
-/*FUNCTION TriaRef::GetNodalFunctionsDerivativesVelocity{{{*/
-void TriaRef::GetNodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss){
-	switch(this->element_type){
-		case P1P1Enum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-			this->element_type = P1P1Enum;
-			return;
-		case P1P1GLSEnum:
-			this->element_type = P1Enum;
-			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-			this->element_type = P1P1GLSEnum;
-			return;
-		case MINIcondensedEnum:
-			this->element_type = P1bubbleEnum;
-			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-			this->element_type = MINIcondensedEnum;
-			return;
-		case MINIEnum:
-			this->element_type = P1bubbleEnum;
-			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-			this->element_type = MINIEnum;
-			return;
-		case TaylorHoodEnum:
-		case XTaylorHoodEnum:
-			this->element_type = P2Enum;
-			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-			this->element_type = XTaylorHoodEnum;
-			return;
-		default:
-			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
-	}
 }
 /*}}}*/
@@ -528,5 +366,5 @@
 			return;
 		default:
-			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
 	}
 
@@ -683,4 +521,14 @@
 
 	return -1;
+}
+/*}}}*/
+/*FUNCTION TriaRef::TensorInterpolation{{{*/
+int TriaRef::TensorInterpolation(void){
+	/*This routine returns the values of the nodal functions  at the gaussian point.*/
+
+	switch(this->element_type){
+		case XTaylorHoodEnum: return P1DGEnum;
+		default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+	}
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 17873)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 17874)
@@ -29,7 +29,4 @@
 		void GetNodalFunctions(IssmDouble* basis,Gauss* gauss);
 		void GetNodalFunctions(IssmDouble* basis,Gauss* gauss,int finiteelement);
-		void GetNodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss);
-		void GetNodalFunctionsPressure(IssmDouble* basis, Gauss* gauss);
-		void GetNodalFunctionsTensor(IssmDouble* basis, Gauss* gauss);
 		void GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss, int index1,int index2);
 		void GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2);
@@ -37,6 +34,4 @@
 		void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss);
 		void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss,int finiteelement);
-		void GetNodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
-		void GetNodalFunctionsDerivativesPressure(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
 		void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss);
 		void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss,int finiteelement);
@@ -52,4 +47,5 @@
 		int  VelocityInterpolation(void);
 		int  PressureInterpolation(void);
+		int  TensorInterpolation(void);
 };
 #endif