Index: /issm/trunk-jpl/src/c/classes/Elements/ElementHook.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/ElementHook.cpp	(revision 18924)
+++ /issm/trunk-jpl/src/c/classes/Elements/ElementHook.cpp	(revision 18925)
@@ -75,10 +75,36 @@
 /*}}}*/
 
+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::SpawnSegHook(ElementHook* triahook,int index1,int index2){/*{{{*/
+
+	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++){
+		/*Do not do anything if Hook is empty*/
+		if (!this->hnodes[i] || this->hnodes[i]->GetNum()==0){
+			triahook->hnodes[i]=NULL;
+		}
+		else{
+			triahook->hnodes[i]=this->hnodes[i]->Spawn(indices,2);
+		}
+	}
+
+	/*do not spawn hmaterial. material will be taken care of by Tria*/
+	triahook->hmaterial=NULL;
+	triahook->hvertices=(Hook*)this->hvertices->Spawn(indices,2);
+	triahook->hmatpar=(Hook*)this->hmatpar->copy();
 }
 /*}}}*/
@@ -111,28 +137,2 @@
 }
 /*}}}*/
-void ElementHook::SpawnSegHook(ElementHook* triahook,int index1,int index2){/*{{{*/
-
-	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++){
-		/*Do not do anything if Hook is empty*/
-		if (!this->hnodes[i] || this->hnodes[i]->GetNum()==0){
-			triahook->hnodes[i]=NULL;
-		}
-		else{
-			triahook->hnodes[i]=this->hnodes[i]->Spawn(indices,2);
-		}
-	}
-
-	/*do not spawn hmaterial. material will be taken care of by Tria*/
-	triahook->hmaterial=NULL;
-	triahook->hvertices=(Hook*)this->hvertices->Spawn(indices,2);
-	triahook->hmatpar=(Hook*)this->hmatpar->copy();
-}
-/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/ElementHook.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/ElementHook.h	(revision 18924)
+++ /issm/trunk-jpl/src/c/classes/Elements/ElementHook.h	(revision 18925)
@@ -24,8 +24,8 @@
 		~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
 };
 
Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 18924)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 18925)
@@ -37,4 +37,148 @@
 
 /*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. 
@@ -103,39 +247,4 @@
 	/*Get Determinant*/
 	Matrix3x3Determinant(Jdet,&J[0][0]);
-	if(*Jdet<0) _error_("negative jacobian determinant!");
-
-}
-/*}}}*/
-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!");
 
@@ -886,62 +995,37 @@
 }
 /*}}}*/
-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::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::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::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!");
 }
 /*}}}*/
@@ -977,4 +1061,89 @@
 }
 /*}}}*/
+int  PentaRef::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 OneLayerP4zEnum:       return P1Enum;
+		case CrouzeixRaviartEnum:   return P1DGEnum;
+		case LACrouzeixRaviartEnum: return NoneEnum;
+		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
+	}
+
+	return -1;
+}
+/*}}}*/
+void PentaRef::SurfaceNodeIndices(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] = 3;
+			indices[1] = 4;
+			indices[2] = 5;
+			break;
+		case P1bubbleEnum: case P1bubblecondensedEnum:
+			numindices = 3;
+			indices    = xNew<int>(numindices);
+			indices[0] = 3;
+			indices[1] = 4;
+			indices[2] = 5;
+			break;
+		case P2xP1Enum:
+			numindices = 6;
+			indices    = xNew<int>(numindices);
+			indices[0] = 3;
+			indices[1] = 4;
+			indices[2] = 5;
+			indices[3] = 9;
+			indices[4] = 10;
+			indices[5] = 11;
+			break;
+		case P1xP2Enum:
+			numindices = 3;
+			indices    = xNew<int>(numindices);
+			indices[0] = 3;
+			indices[1] = 4;
+			indices[2] = 5;
+			return;
+		case P2Enum:
+			numindices = 6;
+			indices    = xNew<int>(numindices);
+			indices[0] = 3;
+			indices[1] = 4;
+			indices[2] = 5;
+			indices[3] = 12;
+			indices[4] = 13;
+			indices[5] = 14;
+			break;
+		default:
+			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
+	}
+
+	/*Assign output pointer*/
+	*pnumindices = numindices;
+	*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){/*{{{*/
 
@@ -995,171 +1164,2 @@
 }
 /*}}}*/
-int  PentaRef::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 OneLayerP4zEnum:       return P1Enum;
-		case CrouzeixRaviartEnum:   return P1DGEnum;
-		case LACrouzeixRaviartEnum: return NoneEnum;
-		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
-	}
-
-	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*/
-	int  numindices;
-	int* indices = NULL;
-
-	switch(finiteelement){
-		case P1Enum: case P1DGEnum:
-			numindices = 3;
-			indices    = xNew<int>(numindices);
-			indices[0] = 3;
-			indices[1] = 4;
-			indices[2] = 5;
-			break;
-		case P1bubbleEnum: case P1bubblecondensedEnum:
-			numindices = 3;
-			indices    = xNew<int>(numindices);
-			indices[0] = 3;
-			indices[1] = 4;
-			indices[2] = 5;
-			break;
-		case P2xP1Enum:
-			numindices = 6;
-			indices    = xNew<int>(numindices);
-			indices[0] = 3;
-			indices[1] = 4;
-			indices[2] = 5;
-			indices[3] = 9;
-			indices[4] = 10;
-			indices[5] = 11;
-			break;
-		case P1xP2Enum:
-			numindices = 3;
-			indices    = xNew<int>(numindices);
-			indices[0] = 3;
-			indices[1] = 4;
-			indices[2] = 5;
-			return;
-		case P2Enum:
-			numindices = 6;
-			indices    = xNew<int>(numindices);
-			indices[0] = 3;
-			indices[1] = 4;
-			indices[2] = 5;
-			indices[3] = 12;
-			indices[4] = 13;
-			indices[5] = 14;
-			break;
-		default:
-			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
-	}
-
-	/*Assign output pointer*/
-	*pnumindices = numindices;
-	*pindices    = indices;
-}
-/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 18924)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 18925)
@@ -15,23 +15,22 @@
 
 		/*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);
+		int  NumberofNodes(int finiteelement);
+		int  PressureInterpolation(int fe_stokes);
 		void SurfaceNodeIndices(int* pnumindices,int** pindices,int finiteelement);
-		int  NumberofNodes(int finiteelement);
+		int  TensorInterpolation(int fe_stokes);
 		int  VelocityInterpolation(int fe_stokes);
-		int  PressureInterpolation(int fe_stokes);
-		int  TensorInterpolation(int fe_stokes);
 };
 #endif
Index: /issm/trunk-jpl/src/c/classes/Elements/SegRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/SegRef.cpp	(revision 18924)
+++ /issm/trunk-jpl/src/c/classes/Elements/SegRef.cpp	(revision 18925)
@@ -29,4 +29,88 @@
 
 /*Reference Element numerics*/
+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 
+	 * point specified by gauss_basis:
+	 *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx
+	 *
+	 * p is a vector 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[1*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.;
+	for(int i=0;i<numnodes;i++) dpx += dbasis[0*numnodes+i]*plist[i];
+
+	/*Assign values*/
+	p[0]=dpx;
+}
+/*}}}*/
+void SegRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* 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 SegRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussSeg* 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];
+
+	*J=.5*fabs(x2-x1);
+}
+/*}}}*/
+void SegRef::GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussSeg* gauss){/*{{{*/
+	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
+	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
+
+	/*Call Jacobian routine to get the jacobian:*/
+	GetJacobian(Jdet, xyz_list, gauss);
+	if(*Jdet<0) _error_("negative jacobian determinant!");
+
+}
+/*}}}*/
+void SegRef::GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussSeg* gauss){/*{{{*/
+
+	/*Jacobian*/
+	IssmDouble J;
+
+	/*Call Jacobian routine to get the jacobian:*/
+	GetJacobian(&J, xyz_list, gauss);
+
+	/*Invert Jacobian matrix: */
+	*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.*/
@@ -109,88 +193,4 @@
 }
 /*}}}*/
-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 
-	 * point specified by gauss_basis:
-	 *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx
-	 *
-	 * p is a vector 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[1*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.;
-	for(int i=0;i<numnodes;i++) dpx += dbasis[0*numnodes+i]*plist[i];
-
-	/*Assign values*/
-	p[0]=dpx;
-}
-/*}}}*/
-void SegRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* 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 SegRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussSeg* 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];
-
-	*J=.5*fabs(x2-x1);
-}
-/*}}}*/
-void SegRef::GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussSeg* gauss){/*{{{*/
-	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
-	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
-
-	/*Call Jacobian routine to get the jacobian:*/
-	GetJacobian(Jdet, xyz_list, gauss);
-	if(*Jdet<0) _error_("negative jacobian determinant!");
-
-}
-/*}}}*/
-void SegRef::GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussSeg* gauss){/*{{{*/
-
-	/*Jacobian*/
-	IssmDouble J;
-
-	/*Call Jacobian routine to get the jacobian:*/
-	GetJacobian(&J, xyz_list, gauss);
-
-	/*Invert Jacobian matrix: */
-	*Jinv = 1./J;
-}
-/*}}}*/
 int  SegRef::NumberofNodes(int finiteelement){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/classes/Elements/SegRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/SegRef.h	(revision 18924)
+++ /issm/trunk-jpl/src/c/classes/Elements/SegRef.h	(revision 18925)
@@ -16,4 +16,6 @@
 		~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);
@@ -22,6 +24,4 @@
 		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);
 };
Index: /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp	(revision 18924)
+++ /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp	(revision 18925)
@@ -31,4 +31,141 @@
 
 /*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.*/
@@ -207,141 +344,4 @@
 	}
 
-}
-/*}}}*/
-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]);
 }
 /*}}}*/
@@ -368,4 +368,28 @@
 }
 /*}}}*/
+int  TetraRef::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;
+		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
+	}
+
+	return -1;
+}/*}}}*/
+int  TetraRef::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  TetraRef::VelocityInterpolation(int fe_stokes){/*{{{*/
 
@@ -384,26 +408,2 @@
 }
 /*}}}*/
-int  TetraRef::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;
-		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
-	}
-
-	return -1;
-}/*}}}*/
-int  TetraRef::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");
-	}
-}
-/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/TetraRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TetraRef.h	(revision 18924)
+++ /issm/trunk-jpl/src/c/classes/Elements/TetraRef.h	(revision 18925)
@@ -16,4 +16,6 @@
 		~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);
@@ -23,11 +25,8 @@
 		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: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 18924)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 18925)
@@ -32,52 +32,58 @@
 
 /*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);
-
-	/*Get nodal functions*/
-	IssmDouble* basis=xNew<IssmDouble>(numnodes);
-	GetNodalFunctions(basis,gauss,finiteelement);
-
-	/*Build B for this segment*/
-	B[0] = +basis[index1];
-	B[1] = +basis[index2];
-	B[2] = -basis[index1];
-	B[3] = -basis[index2];
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(basis);
-}
-/*}}}*/
-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
-	 */
+	_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.;
+	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];
+
+	/*Assign values*/
+	p[0]=dpx;
+	p[1]=dpy;
+
+}
+/*}}}*/
+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);
-
-	/*Get nodal functions*/
-	IssmDouble* basis=xNew<IssmDouble>(numnodes);
-	GetNodalFunctions(basis,gauss,finiteelement);
-
-	/*Build B'*/
-	Bprime[0] = basis[index1];
-	Bprime[1] = basis[index2];
-	Bprime[2] = basis[index1];
-	Bprime[3] = basis[index2];
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(basis);
+	_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;
 }
 /*}}}*/
@@ -97,18 +103,4 @@
 	J[2*0+1] = 0.5*(y2-y1);
 	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!");
-
 }
 /*}}}*/
@@ -193,43 +185,4 @@
 			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
 	}
-}
-/*}}}*/
-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);
 }
 /*}}}*/
@@ -354,130 +307,105 @@
 }
 /*}}}*/
-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);
-
-	/*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];
-
-	/*Assign values*/
-	p[0]=dpx;
-	p[1]=dpy;
-
-}
-/*}}}*/
-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];
+
+	/*Get nodal functions*/
+	IssmDouble* basis=xNew<IssmDouble>(numnodes);
+	GetNodalFunctions(basis,gauss,finiteelement);
+
+	/*Build B for this segment*/
+	B[0] = +basis[index1];
+	B[1] = +basis[index2];
+	B[2] = -basis[index1];
+	B[3] = -basis[index2];
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(basis);
+}
+/*}}}*/
+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
+	 */
 
 	/*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;
-}
-/*}}}*/
-int  TriaRef::NumberofNodes(int finiteelement){/*{{{*/
+
+	/*Get nodal functions*/
+	IssmDouble* basis=xNew<IssmDouble>(numnodes);
+	GetNodalFunctions(basis,gauss,finiteelement);
+
+	/*Build B'*/
+	Bprime[0] = basis[index1];
+	Bprime[1] = basis[index2];
+	Bprime[2] = basis[index1];
+	Bprime[3] = basis[index2];
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(basis);
+}
+/*}}}*/
+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::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 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::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;
-}
-/*}}}*/
-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");
-	}
+		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);
 }
 /*}}}*/
@@ -541,2 +469,74 @@
 }
 /*}}}*/
+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;
+}
+/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 18924)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 18925)
@@ -16,22 +16,21 @@
 
 		/*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