Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17535)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17536)
@@ -3418,8 +3418,10 @@
 
 	/*Initialize vectors*/
-	ElementVector* pe    = element->NewElementVector(FSvelocityEnum);
-	IssmDouble*    Dstar = xNew<IssmDouble>((dim*vnumnodes)*(tausize*tnumnodes));
-	IssmDouble*    tau   = xNew<IssmDouble>(tausize*tnumnodes);
-	IssmDouble*    d     = xNew<IssmDouble>(tausize*tnumnodes);
+	ElementVector* pe      = element->NewElementVector(FSvelocityEnum);
+	IssmDouble*    Dstar   = xNew<IssmDouble>((dim*vnumnodes)*(tausize*tnumnodes));
+	IssmDouble*    tau     = xNew<IssmDouble>(tausize*tnumnodes);
+	IssmDouble*    d       = xNew<IssmDouble>(tausize*tnumnodes);
+	IssmDouble*    vdbasis = xNew<IssmDouble>(dim*vnumnodes);
+	IssmDouble*    tbasis  = xNew<IssmDouble>(tnumnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -3476,4 +3478,7 @@
 
 		/*Create Dstar*/
+		element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
+		element->NodalFunctionsTensor(tbasis,gauss);
+
 		_error_("STOP");
 	}
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 17535)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 17536)
@@ -158,4 +158,5 @@
 		virtual void   NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss)=0;
 		virtual void   NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss)=0;
+		virtual void   NodalFunctionsTensor(IssmDouble* basis, Gauss* gauss)=0;
 		virtual void   NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0;
 		virtual void   NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17535)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17536)
@@ -1725,4 +1725,12 @@
 	_assert_(gauss->Enum()==GaussPentaEnum);
 	this->GetNodalFunctionsPressure(basis,(GaussPenta*)gauss);
+
+}
+/*}}}*/
+/*FUNCTION Penta::NodalFunctionsTensor{{{*/
+void Penta::NodalFunctionsTensor(IssmDouble* basis, Gauss* gauss){
+
+	_assert_(gauss->Enum()==GaussPentaEnum);
+	this->GetNodalFunctionsTensor(basis,(GaussPenta*)gauss);
 
 }
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 17535)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 17536)
@@ -209,4 +209,5 @@
 		void           NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);
 		void           NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss);
+		void           NodalFunctionsTensor(IssmDouble* basis,Gauss* gauss);
 		IssmDouble     MinEdgeLength(IssmDouble* xyz_list);
 		void	         SetClone(int* minranks);
Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 17535)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 17536)
@@ -349,4 +349,15 @@
 	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;
@@ -1192,4 +1203,15 @@
 }
 /*}}}*/
+/*FUNCTION PentaRef::TensorInterpolation{{{*/
+int PentaRef::TensorInterpolation(void){
+
+	switch(this->element_type){
+		case XTaylorHoodEnum:    return P1DGEnum;
+		default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+	}
+
+	return -1;
+}
+/*}}}*/
 /*FUNCTION PentaRef::BasalNodeIndices{{{*/
 void PentaRef::BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){
Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 17535)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 17536)
@@ -26,4 +26,5 @@
 		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);
@@ -55,4 +56,5 @@
 		int  VelocityInterpolation(void);
 		int  PressureInterpolation(void);
+		int  TensorInterpolation(void);
 };
 #endif
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 17535)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 17536)
@@ -84,4 +84,5 @@
 		void        NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
 		void        NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
+		void        NodalFunctionsTensor(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
 		void        NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
 		void        NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp	(revision 17535)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp	(revision 17536)
@@ -512,4 +512,12 @@
 }
 /*}}}*/
+/*FUNCTION Tetra::NodalFunctionsTensor{{{*/
+void Tetra::NodalFunctionsTensor(IssmDouble* basis, Gauss* gauss){
+
+	_assert_(gauss->Enum()==GaussTetraEnum);
+	this->GetNodalFunctionsTensor(basis,(GaussTetra*)gauss);
+
+}
+/*}}}*/
 /*FUNCTION Tetra::NodalFunctionsDerivatives{{{*/
 void Tetra::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){
Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 17535)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 17536)
@@ -88,4 +88,5 @@
 		void        NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);
 		void        NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss);
+		void        NodalFunctionsTensor(IssmDouble* basis,Gauss* gauss);
 		void        NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
 		void        NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp	(revision 17535)
+++ /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp	(revision 17536)
@@ -172,4 +172,19 @@
 }
 /*}}}*/
+/*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");
+	}
+}
+/*}}}*/
 /*FUNCTION TetraRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTetra* gauss){{{*/
 void TetraRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTetra* gauss){
Index: /issm/trunk-jpl/src/c/classes/Elements/TetraRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TetraRef.h	(revision 17535)
+++ /issm/trunk-jpl/src/c/classes/Elements/TetraRef.h	(revision 17536)
@@ -31,4 +31,5 @@
 		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);
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 17535)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 17536)
@@ -1419,4 +1419,12 @@
 }
 /*}}}*/
+/*FUNCTION Tria::NodalFunctionsTensor{{{*/
+void Tria::NodalFunctionsTensor(IssmDouble* basis, Gauss* gauss){
+
+	_assert_(gauss->Enum()==GaussTriaEnum);
+	this->GetNodalFunctionsTensor(basis,(GaussTria*)gauss);
+
+}
+/*}}}*/
 /*FUNCTION Tria::NodalValue {{{*/
 int    Tria::NodalValue(IssmDouble* pvalue, int index, int natureofdataenum){
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 17535)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 17536)
@@ -214,4 +214,5 @@
 		void           NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);
 		void           NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss);
+		void           NodalFunctionsTensor(IssmDouble* basis,Gauss* gauss);
 		void	         SetClone(int* minranks);
 		void           SetTemporaryElementType(int element_type_in){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 17535)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 17536)
@@ -280,4 +280,19 @@
 			this->GetNodalFunctions(basis,gauss);
 			this->element_type = TaylorHoodEnum;
+			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:
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 17535)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 17536)
@@ -31,4 +31,5 @@
 		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);