Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18078)
@@ -26,6 +26,8 @@
 	this->inputs     = NULL;
 	this->parameters = NULL;
+	this->element_type_list=NULL;
 }/*}}}*/
 Element::~Element(){/*{{{*/
+	xDelete<int>(element_type_list);
 	delete inputs;
 }
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 18078)
@@ -46,4 +46,7 @@
 		Matpar      *matpar;
 		Parameters  *parameters;
+
+		int* element_type_list;
+		int  element_type;
 
 	public: 
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18078)
@@ -25,6 +25,5 @@
 /*}}}*/
 Penta::Penta(int penta_id, int penta_sid, int index, IoModel* iomodel,int nummodels)/*{{{*/
-	:PentaRef(nummodels)
-	,ElementHook(nummodels,index+1,NUMVERTICES,iomodel){
+	:ElementHook(nummodels,index+1,NUMVERTICES,iomodel){
 
 	int penta_elements_ids[2];
@@ -35,6 +34,6 @@
 
 	/*id: */
-	this->id=penta_id;
-	this->sid=penta_sid;
+	this->id  = penta_id;
+	this->sid = penta_sid;
 
 	/*Build neighbors list*/
@@ -57,4 +56,7 @@
 	this->matpar            = NULL;
 	this->verticalneighbors = NULL;
+
+	/*Only allocate pointer*/
+	this->element_type_list=xNew<int>(nummodels);
 }
 /*}}}*/
@@ -869,5 +871,5 @@
 
 	_assert_(nodes);
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(this->element_type);
 
 	for(int i=0;i<numnodes;i++){
@@ -879,5 +881,5 @@
 /*}}}*/
 int        Penta::GetNumberOfNodes(void){/*{{{*/
-	return this->NumberofNodes();
+	return this->NumberofNodes(this->element_type);
 }
 /*}}}*/
@@ -888,5 +890,5 @@
 Node* Penta::GetNode(int node_number){/*{{{*/
 	_assert_(node_number>=0); 
-	_assert_(node_number<this->NumberofNodes()); 
+	_assert_(node_number<this->NumberofNodes(this->element_type)); 
 	return this->nodes[node_number];
 }
@@ -1169,5 +1171,5 @@
 		/*Step3: Vertically integrate A COPY of the original*/
 		if(original_input->ObjectEnum()==PentaInputEnum){
-			if(((PentaInput*)original_input)->element_type==P0Enum){
+			if(((PentaInput*)original_input)->interpolation_type==P0Enum){
 				original_input->GetInputValue(&p0top1_list[i]);
 				element_integrated_input= new  PentaInput(original_input->InstanceEnum(),p0top1_list,P1Enum);
@@ -1358,5 +1360,5 @@
 
 	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(this->element_type);
 
 	/*Fetch dof list and allocate solution vector*/
@@ -1645,5 +1647,5 @@
 
 	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctions(basis,(GaussPenta*)gauss);
+	this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->element_type);
 
 }
@@ -1652,5 +1654,5 @@
 
 	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctionsP1(basis,(GaussPenta*)gauss);
+	this->GetNodalFunctions(basis,(GaussPenta*)gauss,P1Enum);
 
 }
@@ -1659,5 +1661,5 @@
 
 	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss);
+	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,this->element_type);
 
 }
@@ -1666,5 +1668,5 @@
 
 	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctionsP1Derivatives(dbasis,xyz_list,(GaussPenta*)gauss);
+	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,P1Enum);
 
 }
@@ -1673,5 +1675,5 @@
 
 	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctionsMINIDerivatives(dbasis,xyz_list,(GaussPenta*)gauss);
+	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,P1bubbleEnum);
 
 }
@@ -2103,5 +2105,5 @@
 
 	/*Recover element type*/
-	this->SetElementType(finiteelement_type,analysis_counter);
+	this->element_type_list[analysis_counter]=finiteelement_type;
 
 	/*Recover vertices ids needed to initialize inputs*/
@@ -2400,5 +2402,5 @@
 
 	GaussPenta* gauss=new GaussPenta();
-	for(int iv=0;iv<this->NumberofNodes();iv++){
+	for(int iv=0;iv<this->NumberofNodes(this->element_type);iv++){
 		gauss->GaussNode(this->element_type,iv);
 		onbase->GetInputValue(&isonbase,gauss);
@@ -2438,5 +2440,5 @@
 /*}}}*/
 void       Penta::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	PentaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss);
+	PentaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss,P1Enum);
 }
 /*}}}*/
@@ -2476,9 +2478,9 @@
 /*}}}*/
 int        Penta::VelocityInterpolation(void){/*{{{*/
-	return PentaRef::VelocityInterpolation();
+	return PentaRef::VelocityInterpolation(this->element_type);
 }
 /*}}}*/
 int        Penta::PressureInterpolation(void){/*{{{*/
-	return PentaRef::PressureInterpolation();
+	return PentaRef::PressureInterpolation(this->element_type);
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 18078)
@@ -28,24 +28,7 @@
 /*Object constructors and destructor*/
 PentaRef::PentaRef(){/*{{{*/
-	this->element_type_list=NULL;
-}
-/*}}}*/
-PentaRef::PentaRef(const int nummodels){/*{{{*/
-
-	/*Only allocate pointer*/
-	element_type_list=xNew<int>(nummodels);
-
 }
 /*}}}*/
 PentaRef::~PentaRef(){/*{{{*/
-	xDelete<int>(element_type_list);
-}
-/*}}}*/
-
-/*Management*/
-void PentaRef::SetElementType(int type,int type_counter){/*{{{*/
-
-	/*initialize element type*/
-	this->element_type_list[type_counter]=type;
 }
 /*}}}*/
@@ -167,12 +150,4 @@
 	/*Invert Jacobian matrix: */
 	Matrix3x3Invert(Jinv,&J[0][0]);
-}
-/*}}}*/
-void PentaRef::GetNodalFunctions(IssmDouble* basis,Gauss* gauss){/*{{{*/
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-	_assert_(basis);
-	GetNodalFunctions(basis,gauss,this->element_type);
-
 }
 /*}}}*/
@@ -321,8 +296,4 @@
 }
 /*}}}*/
-void PentaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss){/*{{{*/
-	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss,this->element_type);
-}
-/*}}}*/
 void PentaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss,int finiteelement){/*{{{*/
 
@@ -356,8 +327,4 @@
 	/*Clean up*/
 	xDelete<IssmDouble>(dbasis_ref);
-}
-/*}}}*/
-void PentaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss){/*{{{*/
-	GetNodalFunctionsDerivativesReference(dbasis,gauss,this->element_type);
 }
 /*}}}*/
@@ -784,155 +751,4 @@
 }
 /*}}}*/
-void PentaRef::GetNodalFunctionsMINIDerivatives(IssmDouble* dbasismini,IssmDouble* xyz_list, Gauss* gauss){/*{{{*/
-
-	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
-	 * actual coordinate system): */
-
-	IssmDouble    dbasismini_ref[3][NUMNODESP1b];
-	IssmDouble    Jinv[3][3];
-
-	/*Get derivative values with respect to parametric coordinate system: */
-	GetNodalFunctionsMINIDerivativesReference(&dbasismini_ref[0][0], gauss); 
-
-	/*Get Jacobian invert: */
-	GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
-
-	/*Build dbasis: 
-	 *
-	 * [dhi/dx]= Jinv'*[dhi/dr]
-	 * [dhi/dy]        [dhi/ds]
-	 * [dhi/dz]        [dhi/dzeta]
-	 */
-
-	for(int i=0;i<NUMNODESP1b;i++){
-		*(dbasismini+NUMNODESP1b*0+i)=Jinv[0][0]*dbasismini_ref[0][i]+Jinv[0][1]*dbasismini_ref[1][i]+Jinv[0][2]*dbasismini_ref[2][i];
-		*(dbasismini+NUMNODESP1b*1+i)=Jinv[1][0]*dbasismini_ref[0][i]+Jinv[1][1]*dbasismini_ref[1][i]+Jinv[1][2]*dbasismini_ref[2][i];
-		*(dbasismini+NUMNODESP1b*2+i)=Jinv[2][0]*dbasismini_ref[0][i]+Jinv[2][1]*dbasismini_ref[1][i]+Jinv[2][2]*dbasismini_ref[2][i];
-	}
-
-}
-/*}}}*/
-void PentaRef::GetNodalFunctionsMINIDerivativesReference(IssmDouble* dbasis,Gauss* gauss_in){/*{{{*/
-	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
-	 * natural coordinate system) at the gaussian point. */
-
-	/*Cast gauss to GaussPenta*/
-	_assert_(gauss_in->Enum()==GaussPentaEnum);
-	GaussPenta* gauss = dynamic_cast<GaussPenta*>(gauss_in);
-
-
-	IssmDouble zeta=gauss->coord4;
-
-	/*Nodal function 1*/
-	dbasis[NUMNODESP1b*0+0]=-0.5*(1.0-zeta)/2.0;
-	dbasis[NUMNODESP1b*1+0]=-SQRT3/6.0*(1.0-zeta)/2.0;
-	dbasis[NUMNODESP1b*2+0]=-0.5*gauss->coord1;
-	/*Nodal function 2*/
-	dbasis[NUMNODESP1b*0+1]=0.5*(1.0-zeta)/2.0;
-	dbasis[NUMNODESP1b*1+1]=-SQRT3/6.0*(1.0-zeta)/2.0;
-	dbasis[NUMNODESP1b*2+1]=-0.5*gauss->coord2;
-	/*Nodal function 3*/
-	dbasis[NUMNODESP1b*0+2]=0.;
-	dbasis[NUMNODESP1b*1+2]=SQRT3/3.0*(1.0-zeta)/2.0;
-	dbasis[NUMNODESP1b*2+2]=-0.5*gauss->coord3;
-	/*Nodal function 4*/
-	dbasis[NUMNODESP1b*0+3]=-0.5*(1.0+zeta)/2.0;
-	dbasis[NUMNODESP1b*1+3]=-SQRT3/6.0*(1.0+zeta)/2.0;
-	dbasis[NUMNODESP1b*2+3]=0.5*gauss->coord1;
-	/*Nodal function 5*/
-	dbasis[NUMNODESP1b*0+4]=0.5*(1.0+zeta)/2.0;
-	dbasis[NUMNODESP1b*1+4]=-SQRT3/6.0*(1.0+zeta)/2.0;
-	dbasis[NUMNODESP1b*2+4]=0.5*gauss->coord2;
-	/*Nodal function 6*/
-	dbasis[NUMNODESP1b*0+5]=0.;
-	dbasis[NUMNODESP1b*1+5]=SQRT3/3.0*(1.0+zeta)/2.0;
-	dbasis[NUMNODESP1b*2+5]=0.5*gauss->coord3;
-	/*Nodal function 7*/
-	dbasis[NUMNODESP1b*0+6]=27.*(1.+zeta)*(1.-zeta)*(-.5*gauss->coord2*gauss->coord3 + .5*gauss->coord1*gauss->coord3);
-	dbasis[NUMNODESP1b*1+6]=27.*(1.+zeta)*(1.-zeta)*SQRT3*(-1./6.*gauss->coord2*gauss->coord3 - 1./6.*gauss->coord1*gauss->coord3 +1./3.*gauss->coord1*gauss->coord2);
-	dbasis[NUMNODESP1b*2+6]=27*gauss->coord1*gauss->coord2*gauss->coord3*(-2.0*zeta);
-}
-/*}}}*/
-void PentaRef::GetNodalFunctionsP1(IssmDouble* basis, Gauss* gauss_in){/*{{{*/
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-	/*Cast gauss to GaussPenta*/
-	_assert_(gauss_in->Enum()==GaussPentaEnum);
-	GaussPenta* gauss = dynamic_cast<GaussPenta*>(gauss_in);
-
-	basis[0]=gauss->coord1*(1-gauss->coord4)/2.0;
-	basis[1]=gauss->coord2*(1-gauss->coord4)/2.0;
-	basis[2]=gauss->coord3*(1-gauss->coord4)/2.0;
-	basis[3]=gauss->coord1*(1+gauss->coord4)/2.0;
-	basis[4]=gauss->coord2*(1+gauss->coord4)/2.0;
-	basis[5]=gauss->coord3*(1+gauss->coord4)/2.0;
-
-}
-/*}}}*/
-void PentaRef::GetNodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss){/*{{{*/
-
-	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
-	 * actual coordinate system): */
-	IssmDouble    dbasis_ref[NDOF3][NUMNODESP1];
-	IssmDouble    Jinv[NDOF3][NDOF3];
-
-	/*Get derivative values with respect to parametric coordinate system: */
-	GetNodalFunctionsP1DerivativesReference(&dbasis_ref[0][0], gauss); 
-
-	/*Get Jacobian invert: */
-	GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
-
-	/*Build basis function derivatives: 
-	 *
-	 * [dhi/dx]= Jinv*[dhi/dr]
-	 * [dhi/dy]       [dhi/ds]
-	 * [dhi/dz]       [dhi/dn]
-	 */
-
-	for (int i=0;i<NUMNODESP1;i++){
-		*(dbasis+NUMNODESP1*0+i)=Jinv[0][0]*dbasis_ref[0][i]+Jinv[0][1]*dbasis_ref[1][i]+Jinv[0][2]*dbasis_ref[2][i];
-		*(dbasis+NUMNODESP1*1+i)=Jinv[1][0]*dbasis_ref[0][i]+Jinv[1][1]*dbasis_ref[1][i]+Jinv[1][2]*dbasis_ref[2][i];
-		*(dbasis+NUMNODESP1*2+i)=Jinv[2][0]*dbasis_ref[0][i]+Jinv[2][1]*dbasis_ref[1][i]+Jinv[2][2]*dbasis_ref[2][i];
-	}
-
-}
-/*}}}*/
-void PentaRef::GetNodalFunctionsP1DerivativesReference(IssmDouble* dbasis,Gauss* gauss_in){/*{{{*/
-
-	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
-	 * natural coordinate system) at the gaussian point. Those values vary along xi,eta,z */
-
-	/*Cast gauss to GaussPenta*/
-	_assert_(gauss_in->Enum()==GaussPentaEnum);
-	GaussPenta* gauss = dynamic_cast<GaussPenta*>(gauss_in);
-
-	IssmDouble zeta=gauss->coord4;
-
-	/*Nodal function 1*/
-	dbasis[NUMNODESP1*0+0]=-0.5*(1.0-zeta)/2.0;
-	dbasis[NUMNODESP1*1+0]=-SQRT3/6.0*(1.0-zeta)/2.0;
-	dbasis[NUMNODESP1*2+0]=-0.5*gauss->coord1;
-	/*Nodal function 2*/
-	dbasis[NUMNODESP1*0+1]=0.5*(1.0-zeta)/2.0;
-	dbasis[NUMNODESP1*1+1]=-SQRT3/6.0*(1.0-zeta)/2.0;
-	dbasis[NUMNODESP1*2+1]=-0.5*gauss->coord2;
-	/*Nodal function 3*/
-	dbasis[NUMNODESP1*0+2]=0.;
-	dbasis[NUMNODESP1*1+2]=SQRT3/3.0*(1.0-zeta)/2.0;
-	dbasis[NUMNODESP1*2+2]=-0.5*gauss->coord3;
-	/*Nodal function 4*/
-	dbasis[NUMNODESP1*0+3]=-0.5*(1.0+zeta)/2.0;
-	dbasis[NUMNODESP1*1+3]=-SQRT3/6.0*(1.0+zeta)/2.0;
-	dbasis[NUMNODESP1*2+3]=0.5*gauss->coord1;
-	/*Nodal function 5*/
-	dbasis[NUMNODESP1*0+4]=0.5*(1.0+zeta)/2.0;
-	dbasis[NUMNODESP1*1+4]=-SQRT3/6.0*(1.0+zeta)/2.0;
-	dbasis[NUMNODESP1*2+4]=0.5*gauss->coord2;
-	/*Nodal function 6*/
-	dbasis[NUMNODESP1*0+5]=0.;
-	dbasis[NUMNODESP1*1+5]=SQRT3/3.0*(1.0+zeta)/2.0;
-	dbasis[NUMNODESP1*2+5]=0.5*gauss->coord3;
-}
-/*}}}*/
 void PentaRef::GetQuadJacobianDeterminant(IssmDouble* Jdet,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*This routine returns the values of the nodal functions  at the gaussian point.*/
@@ -960,10 +776,4 @@
 }
 /*}}}*/
-void PentaRef::GetInputValue(IssmDouble* pvalue,IssmDouble* plist,Gauss* gauss){/*{{{*/
-
-	GetInputValue(pvalue,plist,gauss,this->element_type);
-
-}
-/*}}}*/
 void PentaRef::GetInputValue(IssmDouble* pvalue,IssmDouble* plist,Gauss* gauss,int finiteelement){/*{{{*/
 
@@ -987,5 +797,5 @@
 }
 /*}}}*/
-void PentaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss){/*{{{*/
+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
@@ -1004,9 +814,9 @@
 
 	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(finiteelement);
 
 	/*Get nodal functions derivatives*/
 	IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
-	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss,finiteelement);
 
 	/*Calculate parameter for this Gauss point*/
@@ -1021,9 +831,4 @@
 	p[2]=dpz;
 
-}
-/*}}}*/
-int  PentaRef::NumberofNodes(void){/*{{{*/
-
-	return this->NumberofNodes(this->element_type);
 }
 /*}}}*/
@@ -1046,5 +851,5 @@
 		case P2xP4Enum:             return NUMNODESP2xP4;
 		case P1xP3Enum:             return NUMNODESP1xP3;
-		default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+		default:       _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
 	}
 
@@ -1052,7 +857,7 @@
 }
 /*}}}*/
-int  PentaRef::VelocityInterpolation(void){/*{{{*/
-
-	switch(this->element_type){
+int  PentaRef::VelocityInterpolation(int fe_stokes){/*{{{*/
+
+	switch(fe_stokes){
 		case P1P1Enum:          return P1Enum;
 		case P1P1GLSEnum:       return P1Enum;
@@ -1061,5 +866,5 @@
 		case TaylorHoodEnum:    return P2Enum;
 		case OneLayerP4zEnum:   return P2xP4Enum;
-		default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
 	}
 
@@ -1067,7 +872,7 @@
 }
 /*}}}*/
-int  PentaRef::PressureInterpolation(void){/*{{{*/
-
-	switch(this->element_type){
+int  PentaRef::PressureInterpolation(int fe_stokes){/*{{{*/
+
+	switch(fe_stokes){
 		case P1P1Enum:          return P1Enum;
 		case P1P1GLSEnum:       return P1Enum;
@@ -1076,5 +881,5 @@
 		case TaylorHoodEnum:    return P1Enum;
 		case OneLayerP4zEnum:   return P1Enum;
-		default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
 	}
 
@@ -1082,9 +887,9 @@
 }
 /*}}}*/
-int  PentaRef::TensorInterpolation(void){/*{{{*/
-
-	switch(this->element_type){
+int  PentaRef::TensorInterpolation(int fe_stokes){/*{{{*/
+
+	switch(fe_stokes){
 		case XTaylorHoodEnum:    return P1DGEnum;
-		default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+		default: _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
 	}
 
@@ -1158,5 +963,5 @@
 			break;
 		default:
-			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
 	}
 
@@ -1215,5 +1020,5 @@
 			break;
 		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 18077)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 18078)
@@ -11,26 +11,11 @@
 
 	public: 
-		int* element_type_list; //P1CG, P1DG, MINI, P2...
-		int  element_type;
-
 		PentaRef();
-		PentaRef(const int nummodels);
 		~PentaRef();
 
-		/*Management*/
-		void SetElementType(int type,int type_counter);
-
 		/*Numerics*/
-		void GetNodalFunctions(IssmDouble* basis, Gauss* gauss);
 		void GetNodalFunctions(IssmDouble* basis, Gauss* gauss,int finiteelement);
-		void GetNodalFunctionsDerivatives(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* basis, Gauss* gauss);
-		void GetNodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss);
-		void GetNodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss);
-		void GetNodalFunctionsP1DerivativesReference(IssmDouble* dl1dl6,Gauss* gauss);
-		void GetNodalFunctionsMINIDerivativesReference(IssmDouble* dl1dl7,Gauss* gauss);
 		void GetQuadJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
 		void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss);
@@ -40,15 +25,13 @@
 		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);
 		void GetInputValue(IssmDouble* pvalue,IssmDouble* plist, Gauss* gauss,int finiteelement);
-		void GetInputDerivativeValue(IssmDouble* pvalues, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss);
+		void GetInputDerivativeValue(IssmDouble* pvalues, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement);
 
 		void BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement);
 		void SurfaceNodeIndices(int* pnumindices,int** pindices,int finiteelement);
-		int  NumberofNodes(void);
 		int  NumberofNodes(int finiteelement);
-		int  VelocityInterpolation(void);
-		int  PressureInterpolation(void);
-		int  TensorInterpolation(void);
+		int  VelocityInterpolation(int fe_stokes);
+		int  PressureInterpolation(int fe_stokes);
+		int  TensorInterpolation(int fe_stokes);
 };
 #endif
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp	(revision 18078)
@@ -20,5 +20,5 @@
 /*Constructors/destructor/copy*/
 Seg::Seg(int seg_id, int seg_sid, int index, IoModel* iomodel,int nummodels)/*{{{*/
-		:SegRef(nummodels),ElementHook(nummodels,index+1,NUMVERTICES,iomodel){
+		:ElementHook(nummodels,index+1,NUMVERTICES,iomodel){
 
 			/*id: */
@@ -37,4 +37,7 @@
 			this->material = NULL;
 			this->matpar   = NULL;
+
+			/*Only allocate pointer*/
+			this->element_type_list=xNew<int>(nummodels);
 		}
 /*}}}*/
@@ -100,5 +103,5 @@
 }/*}}}*/
 int        Seg::GetNumberOfNodes(void){/*{{{*/
-	return this->NumberofNodes();
+	return this->NumberofNodes(this->element_type);
 }
 /*}}}*/
@@ -177,5 +180,5 @@
 
 	_assert_(gauss->Enum()==GaussSegEnum);
-	this->GetNodalFunctions(basis,(GaussSeg*)gauss);
+	this->GetNodalFunctions(basis,(GaussSeg*)gauss,this->element_type);
 
 }
@@ -184,5 +187,5 @@
 
 	_assert_(gauss->Enum()==GaussSegEnum);
-	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussSeg*)gauss);
+	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussSeg*)gauss,this->element_type);
 
 }
Index: /issm/trunk-jpl/src/c/classes/Elements/SegRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/SegRef.cpp	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Elements/SegRef.cpp	(revision 18078)
@@ -21,36 +21,11 @@
 /*Object constructors and destructor*/
 SegRef::SegRef(){/*{{{*/
-	this->element_type_list=NULL;
-}
-/*}}}*/
-SegRef::SegRef(const int nummodels){/*{{{*/
-
-	/*Only allocate pointer*/
-	element_type_list=xNew<int>(nummodels);
-
 }
 /*}}}*/
 SegRef::~SegRef(){/*{{{*/
-	xDelete<int>(element_type_list);
-}
-/*}}}*/
-
-/*Management*/
-void SegRef::SetElementType(int type,int type_counter){/*{{{*/
-
-	/*initialize element type*/
-	this->element_type_list[type_counter]=type;
 }
 /*}}}*/
 
 /*Reference Element numerics*/
-void SegRef::GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss){/*{{{*/
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-	_assert_(basis);
-
-	GetNodalFunctions(basis,gauss,this->element_type);
-}
-/*}}}*/
 void SegRef::GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss,int finiteelement){/*{{{*/
 	/*This routine returns the values of the nodal functions  at the gaussian point.*/
@@ -74,10 +49,4 @@
 			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
 	}
-}
-/*}}}*/
-void SegRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussSeg* gauss){/*{{{*/
-
-	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss,this->element_type);
-
 }
 /*}}}*/
@@ -107,12 +76,4 @@
 	/*Clean up*/
 	xDelete<IssmDouble>(dbasis_ref);
-
-}
-/*}}}*/
-void SegRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussSeg* gauss){/*{{{*/
-	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
-	 * natural coordinate system) at the gaussian point. */
-
-	GetNodalFunctionsDerivativesReference(dbasis,gauss,this->element_type);
 
 }
@@ -147,5 +108,5 @@
 }
 /*}}}*/
-void SegRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussSeg* gauss){/*{{{*/
+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 
@@ -160,9 +121,9 @@
 
 	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(finiteelement);
 
 	/*Get nodal functions derivatives*/
 	IssmDouble* dbasis=xNew<IssmDouble>(1*numnodes);
-	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss,finiteelement);
 
 	/*Calculate parameter for this Gauss point*/
@@ -175,9 +136,4 @@
 }
 /*}}}*/
-void SegRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* gauss){/*{{{*/
-
-	GetInputValue(p,plist,gauss,this->element_type);
-}
-/*}}}*/
 void SegRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* gauss,int finiteelement){/*{{{*/
 
@@ -230,9 +186,4 @@
 	/*Invert Jacobian matrix: */
 	*Jinv = 1./J;
-}
-/*}}}*/
-int  SegRef::NumberofNodes(void){/*{{{*/
-
-	return this->NumberofNodes(this->element_type);
 }
 /*}}}*/
@@ -243,5 +194,5 @@
 		case P1Enum:                return NUMNODESP1;
 		case P1DGEnum:              return NUMNODESP1;
-		default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+		default: _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
 	}
 
Index: /issm/trunk-jpl/src/c/classes/Elements/SegRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/SegRef.h	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Elements/SegRef.h	(revision 18078)
@@ -13,27 +13,15 @@
 
 	public: 
-		int* element_type_list;
-		int  element_type;
-
 		SegRef();
-		SegRef(const int nummodels);
 		~SegRef();
 
-		/*Management*/
-		void SetElementType(int type,int type_counter);
 		void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussSeg* gauss);
 		void GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussSeg* gauss);
 		void GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussSeg* gauss);
-		void GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss);
 		void GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss,int finiteelement);
-		void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussSeg* gauss);
 		void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussSeg* gauss,int finiteelement);
-		void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussSeg* gauss);
 		void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussSeg* gauss,int finiteelement);
-		void GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussSeg* gauss);
-		void GetInputValue(IssmDouble* p, IssmDouble* plist, GaussSeg* gauss);
+		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(void);
 		int  NumberofNodes(int finiteelement);
 };
Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp	(revision 18078)
@@ -21,5 +21,5 @@
 /*Constructors/destructor/copy*/
 Tetra::Tetra(int seg_id, int seg_sid, int index, IoModel* iomodel,int nummodels)/*{{{*/
-		:TetraRef(nummodels),ElementHook(nummodels,index+1,NUMVERTICES,iomodel){
+		:ElementHook(nummodels,index+1,NUMVERTICES,iomodel){
 
 			/*id: */
@@ -38,4 +38,7 @@
 			this->material = NULL;
 			this->matpar   = NULL;
+
+			/*Only allocate pointer*/
+			this->element_type_list=xNew<int>(nummodels);
 		}
 /*}}}*/
@@ -202,5 +205,5 @@
 
 	_assert_(nodes);
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(this->element_type);
 
 	for(int i=0;i<numnodes;i++){
@@ -212,5 +215,5 @@
 /*}}}*/
 int      Tetra::GetNumberOfNodes(void){/*{{{*/
-	return this->NumberofNodes();
+	return this->NumberofNodes(this->element_type);
 }
 /*}}}*/
@@ -400,5 +403,5 @@
 
 	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(this->element_type);
 
 	/*Fetch dof list and allocate solution vector*/
@@ -502,5 +505,5 @@
 
 	_assert_(gauss->Enum()==GaussTetraEnum);
-	this->GetNodalFunctions(basis,(GaussTetra*)gauss);
+	this->GetNodalFunctions(basis,(GaussTetra*)gauss,this->element_type);
 
 }
@@ -530,5 +533,5 @@
 
 	_assert_(gauss->Enum()==GaussTetraEnum);
-	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTetra*)gauss);
+	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTetra*)gauss,this->element_type);
 
 }
@@ -788,5 +791,5 @@
 
 	/*Recover element type*/
-	this->SetElementType(finiteelement_type,analysis_counter);
+	this->element_type_list[analysis_counter]=finiteelement_type;
 
 	/*Recover vertices ids needed to initialize inputs*/
@@ -872,5 +875,5 @@
 /*}}}*/
 int      Tetra::VelocityInterpolation(void){/*{{{*/
-	return TetraRef::VelocityInterpolation();
+	return TetraRef::VelocityInterpolation(this->element_type);
 }
 /*}}}*/
@@ -892,9 +895,9 @@
 /*}}}*/
 int      Tetra::PressureInterpolation(void){/*{{{*/
-	return TetraRef::PressureInterpolation();
+	return TetraRef::PressureInterpolation(this->element_type);
 }
 /*}}}*/
 int      Tetra::TensorInterpolation(void){/*{{{*/
-	return TetraRef::TensorInterpolation();
+	return TetraRef::TensorInterpolation(this->element_type);
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp	(revision 18078)
@@ -4,5 +4,5 @@
 
 /*Headers:*/
-/*{{{*//*{{{*/
+/*{{{*/
 #ifdef HAVE_CONFIG_H
 #include <config.h>
@@ -23,34 +23,11 @@
 /*Object constructors and destructor*/
 TetraRef::TetraRef(){/*{{{*/
-	this->element_type_list=NULL;
-}
-/*}}}*/
-TetraRef::TetraRef(const int nummodels){/*{{{*/
-
-	/*Only allocate pointer*/
-	element_type_list=xNew<int>(nummodels);
-
 }
 /*}}}*/
 TetraRef::~TetraRef(){/*{{{*/
-	xDelete<int>(element_type_list);
-}
-/*}}}*/
-
-/*Management*/
-void TetraRef::SetElementType(int type,int type_counter){/*{{{*/
-
-	/*initialize element type*/
-	this->element_type_list[type_counter]=type;
 }
 /*}}}*/
 
 /*Reference Element numerics*/
-void TetraRef::GetNodalFunctions(IssmDouble* basis,GaussTetra* gauss){/*{{{*/
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-	_assert_(basis);
-	GetNodalFunctions(basis,gauss,this->element_type);
-}
-/*}}}*/
 void TetraRef::GetNodalFunctions(IssmDouble* basis,GaussTetra* gauss,int finiteelement){/*{{{*/
 	/*This routine returns the values of the nodal functions  at the gaussian point.*/
@@ -96,8 +73,4 @@
 }
 /*}}}*/
-void TetraRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTetra* gauss){/*{{{*/
-	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss,this->element_type);
-}
-/*}}}*/
 void TetraRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTetra* gauss,int finiteelement){/*{{{*/
 
@@ -133,12 +106,4 @@
 }
 /*}}}*/
-void TetraRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTetra* gauss){/*{{{*/
-	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
-	 * natural coordinate system) at the gaussian point. */
-
-	GetNodalFunctionsDerivativesReference(dbasis,gauss,this->element_type);
-
-}
-/*}}}*/
 void TetraRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTetra* gauss,int finiteelement){/*{{{*/
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
@@ -239,5 +204,5 @@
 }
 /*}}}*/
-void TetraRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussTetra* gauss){/*{{{*/
+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
@@ -256,9 +221,9 @@
 
 	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(finiteelement);
 
 	/*Get nodal functions derivatives*/
 	IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
-	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss,finiteelement);
 
 	/*Calculate parameter for this Gauss point*/
@@ -272,9 +237,4 @@
 	p[1]=dpy;
 	p[2]=dpz;
-}
-/*}}}*/
-void TetraRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTetra* gauss){/*{{{*/
-
-	GetInputValue(p,plist,gauss,this->element_type);
 }
 /*}}}*/
@@ -374,9 +334,4 @@
 	/*Invert Jacobian matrix: */
 	Matrix3x3Invert(Jinv,&J[0][0]);
-}
-/*}}}*/
-int  TetraRef::NumberofNodes(void){/*{{{*/
-
-	return this->NumberofNodes(this->element_type);
 }
 /*}}}*/
@@ -395,5 +350,5 @@
 		case MINIEnum:              return NUMNODESP1b+NUMNODESP1;
 		case TaylorHoodEnum:        return NUMNODESP2+NUMNODESP1;
-		default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+		default: _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
 	}
 
@@ -401,7 +356,7 @@
 }
 /*}}}*/
-int  TetraRef::VelocityInterpolation(void){/*{{{*/
-
-	switch(this->element_type){
+int  TetraRef::VelocityInterpolation(int fe_stokes){/*{{{*/
+
+	switch(fe_stokes){
 		case P1P1Enum:          return P1Enum;
 		case P1P1GLSEnum:       return P1Enum;
@@ -409,5 +364,5 @@
 		case MINIEnum:          return P1bubbleEnum;
 		case TaylorHoodEnum:    return P2Enum;
-		default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
 	}
 
@@ -415,7 +370,7 @@
 }
 /*}}}*/
-int TetraRef::PressureInterpolation(void){/*{{{*/
-
-	switch(this->element_type){
+int TetraRef::PressureInterpolation(int fe_stokes){/*{{{*/
+
+	switch(fe_stokes){
 		case P1P1Enum:          return P1Enum;
 		case P1P1GLSEnum:       return P1Enum;
@@ -423,16 +378,16 @@
 		case MINIEnum:          return P1Enum;
 		case TaylorHoodEnum:    return P1Enum;
-		default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
 	}
 
 	return -1;
 }/*}}}*/
-int  TetraRef::TensorInterpolation(void){/*{{{*/
+int  TetraRef::TensorInterpolation(int fe_stokes){/*{{{*/
 	/*This routine returns the values of the nodal functions  at the gaussian point.*/
 
-	switch(this->element_type){
+	switch(fe_stokes){
 		case XTaylorHoodEnum: return P1DGEnum;
-		default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
-	}
-}
-/*}}}*/
+		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 18077)
+++ /issm/trunk-jpl/src/c/classes/Elements/TetraRef.h	(revision 18078)
@@ -13,32 +13,21 @@
 
 	public: 
-		int* element_type_list;
-		int  element_type;
-
 		TetraRef();
-		TetraRef(const int nummodels);
 		~TetraRef();
 
-		/*Management*/
-		void SetElementType(int type,int type_counter);
 		void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussTetra* gauss);
 		void GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTetra* gauss);
 		void GetJacobianDeterminantFace(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTetra* gauss);
 		void GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussTetra* gauss);
-		void GetNodalFunctions(IssmDouble* basis,GaussTetra* gauss);
 		void GetNodalFunctions(IssmDouble* basis,GaussTetra* gauss,int finiteelement);
-		void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTetra* gauss);
 		void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTetra* gauss,int finiteelement);
-		void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTetra* gauss);
 		void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTetra* gauss,int finiteelement);
-		void GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussTetra* gauss);
-		void GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTetra* gauss);
+		void GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussTetra* gauss,int finiteelement);
 		void GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTetra* gauss,int finiteelement);
 
-		int  NumberofNodes(void);
 		int  NumberofNodes(int finiteelement);
-		int  VelocityInterpolation(void);
-		int  PressureInterpolation(void);
-		int  TensorInterpolation(void);
+		int  VelocityInterpolation(int fe_stokes);
+		int  PressureInterpolation(int fe_stokes);
+		int  TensorInterpolation(int fe_stokes);
 };
 #endif
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18078)
@@ -25,5 +25,5 @@
 /*Constructors/destructor/copy*/
 Tria::Tria(int tria_id, int tria_sid, int index, IoModel* iomodel,int nummodels)/*{{{*/
-	:TriaRef(nummodels),ElementHook(nummodels,index+1,NUMVERTICES,iomodel){
+	:ElementHook(nummodels,index+1,NUMVERTICES,iomodel){
 
 		/*id: */
@@ -42,4 +42,5 @@
 		this->material = NULL;
 		this->matpar   = NULL;
+		this->element_type_list=xNew<int>(nummodels);
 }
 /*}}}*/
@@ -900,5 +901,5 @@
 /*}}}*/
 int        Tria::GetNumberOfNodes(void){/*{{{*/
-	return this->NumberofNodes();
+	return this->NumberofNodes(this->element_type);
 }
 /*}}}*/
@@ -921,5 +922,5 @@
 Node*      Tria::GetNode(int node_number){/*{{{*/
 	_assert_(node_number>=0); 
-	_assert_(node_number<this->NumberofNodes()); 
+	_assert_(node_number<this->NumberofNodes(this->element_type)); 
 	return this->nodes[node_number];
 
@@ -1059,5 +1060,5 @@
 
 	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(this->element_type);
 
 	/*Fetch dof list and allocate solution vector*/
@@ -1110,5 +1111,5 @@
 
 		/*Get number of nodes and dof list: */
-		numnodes = this->NumberofNodes();
+		numnodes = this->NumberofNodes(this->element_type);
 		values   = xNew<IssmDouble>(numnodes);
 		GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
@@ -1124,5 +1125,5 @@
 
 		/*Get number of nodes and dof list: */
-		numnodes = this->NumberofNodes();
+		numnodes = this->NumberofNodes(this->element_type);
 		values   = xNew<IssmDouble>(numnodes);
 
@@ -1453,5 +1454,5 @@
 
 	_assert_(gauss->Enum()==GaussTriaEnum);
-	this->GetNodalFunctions(basis,(GaussTria*)gauss);
+	this->GetNodalFunctions(basis,(GaussTria*)gauss,this->element_type);
 
 }
@@ -1467,5 +1468,5 @@
 
 	_assert_(gauss->Enum()==GaussTriaEnum);
-	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss);
+	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->element_type);
 
 }
@@ -1573,13 +1574,13 @@
 /*}}}*/
 int        Tria::VelocityInterpolation(void){/*{{{*/
-	return TriaRef::VelocityInterpolation();
+	return TriaRef::VelocityInterpolation(this->element_type);
 }
 /*}}}*/
 int        Tria::PressureInterpolation(void){/*{{{*/
-	return TriaRef::PressureInterpolation();
+	return TriaRef::PressureInterpolation(this->element_type);
 }
 /*}}}*/
 int        Tria::TensorInterpolation(void){/*{{{*/
-	return TriaRef::TensorInterpolation();
+	return TriaRef::TensorInterpolation(this->element_type);
 }
 /*}}}*/
@@ -1890,5 +1891,5 @@
 
 	/*Recover element type*/
-	this->SetElementType(finiteelement_type,analysis_counter);
+	this->element_type_list[analysis_counter]=finiteelement_type;
 
 	/*Recover nodes ids needed to initialize the node hook.*/
@@ -1988,5 +1989,5 @@
 
 	GaussTria* gauss=new GaussTria();
-	for(int iv=0;iv<this->NumberofNodes();iv++){
+	for(int iv=0;iv<this->NumberofNodes(this->element_type);iv++){
 		gauss->GaussNode(this->element_type,iv);
 		onbase->GetInputValue(&isonbase,gauss);
@@ -2027,5 +2028,5 @@
 /*}}}*/
 void       Tria::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	TriaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss);
+	TriaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss,P1Enum);
 }
 /*}}}*/
@@ -2828,5 +2829,5 @@
 
 	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(this->element_type);
 
 	/*Fetch dof list and allocate solution vector*/
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 18078)
@@ -23,29 +23,12 @@
 /*Object constructors and destructor*/
 TriaRef::TriaRef(){/*{{{*/
-	this->element_type_list=NULL;
-}
-/*}}}*/
-TriaRef::TriaRef(const int nummodels){/*{{{*/
-
-	/*Only allocate pointer*/
-	element_type_list=xNew<int>(nummodels);
-
 }
 /*}}}*/
 TriaRef::~TriaRef(){/*{{{*/
-	xDelete<int>(element_type_list);
-}
-/*}}}*/
-
-/*Management*/
-void TriaRef::SetElementType(int type,int type_counter){/*{{{*/
-
-	/*initialize element type*/
-	this->element_type_list[type_counter]=type;
 }
 /*}}}*/
 
 /*Reference Element numerics*/
-void TriaRef::GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2){/*{{{*/
+void TriaRef::GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/
 	/*Compute B  matrix. B=[phi1 phi2 -phi3 -phi4]
 	 *
@@ -56,9 +39,9 @@
 
 	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(finiteelement);
 
 	/*Get nodal functions*/
 	IssmDouble* basis=xNew<IssmDouble>(numnodes);
-	GetNodalFunctions(basis,gauss);
+	GetNodalFunctions(basis,gauss,finiteelement);
 
 	/*Build B for this segment*/
@@ -72,5 +55,5 @@
 }
 /*}}}*/
-void TriaRef::GetSegmentBprimeFlux(IssmDouble* Bprime,Gauss* gauss, int index1,int index2){/*{{{*/
+void TriaRef::GetSegmentBprimeFlux(IssmDouble* Bprime,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/
 	/*Compute Bprime  matrix. Bprime=[phi1 phi2 phi3 phi4]
 	 *
@@ -81,9 +64,9 @@
 
 	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(finiteelement);
 
 	/*Get nodal functions*/
 	IssmDouble* basis=xNew<IssmDouble>(numnodes);
-	GetNodalFunctions(basis,gauss);
+	GetNodalFunctions(basis,gauss,finiteelement);
 
 	/*Build B'*/
@@ -152,12 +135,4 @@
 	/*Invert Jacobian matrix: */
 	Matrix2x2Invert(Jinv,&J[0][0]);
-
-}
-/*}}}*/
-void TriaRef::GetNodalFunctions(IssmDouble* basis,Gauss* gauss){/*{{{*/
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-	_assert_(basis);
-	GetNodalFunctions(basis,gauss,this->element_type);
 
 }
@@ -204,5 +179,5 @@
 }
 /*}}}*/
-void TriaRef::GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss,int index1,int index2){/*{{{*/
+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.*/
 
@@ -211,11 +186,11 @@
 
 	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(finiteelement);
 
 	/*Get nodal functions*/
 	IssmDouble* triabasis=xNew<IssmDouble>(numnodes);
-	GetNodalFunctions(triabasis,gauss);
-
-	switch(this->element_type){
+	GetNodalFunctions(triabasis,gauss,finiteelement);
+
+	switch(finiteelement){
 		case P1Enum: case P1DGEnum:
 			basis[0]=triabasis[index1];
@@ -236,15 +211,9 @@
 			return;
 		default:
-			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+			_error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
 	}
 
 	/*Clean up*/
 	xDelete<IssmDouble>(triabasis);
-}
-/*}}}*/
-void TriaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss){/*{{{*/
-
-	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss,this->element_type);
-
 }
 /*}}}*/
@@ -276,12 +245,4 @@
 	/*Clean up*/
 	xDelete<IssmDouble>(dbasis_ref);
-
-}
-/*}}}*/
-void TriaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss){/*{{{*/
-	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
-	 * natural coordinate system) at the gaussian point. */
-
-	GetNodalFunctionsDerivativesReference(dbasis,gauss,this->element_type);
 
 }
@@ -354,5 +315,5 @@
 }
 /*}}}*/
-void TriaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss){/*{{{*/
+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 
@@ -369,9 +330,9 @@
 
 	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(finiteelement);
 
 	/*Get nodal functions derivatives*/
 	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
-	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss,finiteelement);
 
 	/*Calculate parameter for this Gauss point*/
@@ -386,9 +347,4 @@
 }
 /*}}}*/
-void TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, Gauss* gauss){/*{{{*/
-
-	GetInputValue(p,plist,gauss,this->element_type);
-}
-/*}}}*/
 void TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, Gauss* gauss,int finiteelement){/*{{{*/
 
@@ -409,9 +365,4 @@
 	xDelete<IssmDouble>(basis);
 	*p = value;
-}
-/*}}}*/
-int  TriaRef::NumberofNodes(void){/*{{{*/
-
-	return this->NumberofNodes(this->element_type);
 }
 /*}}}*/
@@ -431,5 +382,5 @@
 		case TaylorHoodEnum:        return NUMNODESP2+NUMNODESP1;
 		case XTaylorHoodEnum:       return NUMNODESP2+NUMNODESP1;
-		default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+		default: _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
 	}
 
@@ -437,7 +388,7 @@
 }
 /*}}}*/
-int  TriaRef::VelocityInterpolation(void){/*{{{*/
-
-	switch(this->element_type){
+int  TriaRef::VelocityInterpolation(int fe_stokes){/*{{{*/
+
+	switch(fe_stokes){
 		case P1P1Enum:          return P1Enum;
 		case P1P1GLSEnum:       return P1Enum;
@@ -446,5 +397,5 @@
 		case TaylorHoodEnum:    return P2Enum;
 		case XTaylorHoodEnum:   return P2Enum;
-		default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
 	}
 
@@ -452,7 +403,7 @@
 }
 /*}}}*/
-int  TriaRef::PressureInterpolation(void){/*{{{*/
-
-	switch(this->element_type){
+int  TriaRef::PressureInterpolation(int fe_stokes){/*{{{*/
+
+	switch(fe_stokes){
 		case P1P1Enum:          return P1Enum;
 		case P1P1GLSEnum:       return P1Enum;
@@ -461,5 +412,5 @@
 		case TaylorHoodEnum:    return P1Enum;
 		case XTaylorHoodEnum:   return P1Enum;
-		default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
 	}
 
@@ -467,10 +418,10 @@
 }
 /*}}}*/
-int  TriaRef::TensorInterpolation(void){/*{{{*/
+int  TriaRef::TensorInterpolation(int fe_stokes){/*{{{*/
 	/*This routine returns the values of the nodal functions  at the gaussian point.*/
 
-	switch(this->element_type){
+	switch(fe_stokes){
 		case XTaylorHoodEnum: return P1DGEnum;
-		default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+		default: _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
 	}
 }
@@ -527,5 +478,5 @@
 			break;
 		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/TriaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 18078)
@@ -12,13 +12,6 @@
 
 	public: 
-		int* element_type_list; //P1CG, P1DG, MINI, P2...
-		int  element_type;
-
 		TriaRef();
-		TriaRef(const int nummodels);
 		~TriaRef();
-
-		/*Management*/
-		void SetElementType(int type,int type_counter);
 
 		/*Numerics*/
@@ -27,23 +20,18 @@
 		void GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss);
 		void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,Gauss* gauss);
-		void GetNodalFunctions(IssmDouble* basis,Gauss* gauss);
 		void GetNodalFunctions(IssmDouble* basis,Gauss* gauss,int finiteelement);
-		void GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss, int index1,int index2);
-		void GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2);
-		void GetSegmentBprimeFlux(IssmDouble* Bprime,Gauss* gauss, int index1,int index2);
-		void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss);
+		void GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss, int index1,int index2,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);
 		void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss,int finiteelement);
-		void GetInputValue(IssmDouble* pp, IssmDouble* plist, Gauss* gauss);
 		void GetInputValue(IssmDouble* pp, IssmDouble* plist, Gauss* gauss,int finiteelement);
-		void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss);
+		void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement);
 
 		void NodeOnEdgeIndices(int* pnumindices,int** pindices,int index,int finiteelement);
-		int  NumberofNodes(void);
 		int  NumberofNodes(int finiteelement);
-		int  VelocityInterpolation(void);
-		int  PressureInterpolation(void);
-		int  TensorInterpolation(void);
+		int  VelocityInterpolation(int fe_stokes);
+		int  PressureInterpolation(int fe_stokes);
+		int  TensorInterpolation(int fe_stokes);
 };
 #endif
Index: /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp	(revision 18078)
@@ -17,18 +17,13 @@
 }
 /*}}}*/
-PentaInput::PentaInput(int in_enum_type,IssmDouble* in_values,int element_type_in)/*{{{*/
-		:PentaRef(1)
-{
-
-	/*Set PentaRef*/
-	this->SetElementType(element_type_in,0);
-	this->element_type=element_type_in;
+PentaInput::PentaInput(int in_enum_type,IssmDouble* in_values,int interpolation_type_in){/*{{{*/
 
 	/*Set Enum*/
 	enum_type=in_enum_type;
+	this->interpolation_type=interpolation_type_in;
 
 	/*Set values*/
-	this->values=xNew<IssmDouble>(this->NumberofNodes());
-	for(int i=0;i<this->NumberofNodes();i++) values[i]=in_values[i];
+	this->values=xNew<IssmDouble>(this->NumberofNodes(this->interpolation_type));
+	for(int i=0;i<this->NumberofNodes(this->interpolation_type);i++) values[i]=in_values[i];
 }
 /*}}}*/
@@ -46,5 +41,5 @@
 
 	_printf_(setw(15)<<"   PentaInput "<<setw(25)<<left<<EnumToStringx(this->enum_type)<<" [");
-	for(int i=0;i<this->NumberofNodes();i++) _printf_(" "<<this->values[i]);
+	for(int i=0;i<this->NumberofNodes(this->interpolation_type);i++) _printf_(" "<<this->values[i]);
 	_printf_("]\n");
 }
@@ -60,5 +55,5 @@
 Object* PentaInput::copy() {/*{{{*/
 
-	return new PentaInput(this->enum_type,this->values,this->element_type);
+	return new PentaInput(this->enum_type,this->values,this->interpolation_type);
 
 }
@@ -77,5 +72,5 @@
 	TriaInput* outinput=NULL;
 
-	if(this->element_type==P0Enum){ 
+	if(this->interpolation_type==P0Enum){ 
 		outinput=new TriaInput(this->enum_type,&this->values[0],P0Enum);
 	}
@@ -109,5 +104,5 @@
 int  PentaInput::GetResultInterpolation(void){/*{{{*/
 
-	if(this->element_type==P0Enum){
+	if(this->interpolation_type==P0Enum){
 		return P0Enum;
 	}
@@ -118,5 +113,5 @@
 int  PentaInput::GetResultNumberOfNodes(void){/*{{{*/
 
-	return this->NumberofNodes();;
+	return this->NumberofNodes(this->interpolation_type);;
 
 }
@@ -124,5 +119,5 @@
 void PentaInput::ResultToPatch(IssmDouble* values,int nodesperelement,int sid){/*{{{*/
 
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(this->interpolation_type);
 
 	/*Some checks*/
@@ -138,5 +133,5 @@
 void PentaInput::GetInputValue(IssmDouble* pvalue){/*{{{*/
 
-	if(this->element_type==P0Enum){
+	if(this->interpolation_type==P0Enum){
 		pvalue=&values[0];
 	}
@@ -148,5 +143,5 @@
 	/*Call PentaRef function*/
 	_assert_(gauss->Enum()==GaussPentaEnum);
-	PentaRef::GetInputValue(pvalue,&values[0],(GaussPenta*)gauss);
+	PentaRef::GetInputValue(pvalue,&values[0],(GaussPenta*)gauss,this->interpolation_type);
 
 }
@@ -156,5 +151,5 @@
 	/*Call PentaRef function*/
 	_assert_(gauss->Enum()==GaussPentaEnum);
-	PentaRef::GetInputDerivativeValue(p,&values[0],xyz_list,(GaussPenta*)gauss);
+	PentaRef::GetInputDerivativeValue(p,&values[0],xyz_list,(GaussPenta*)gauss,this->interpolation_type);
 }
 /*}}}*/
@@ -165,5 +160,5 @@
 void PentaInput::GetInputAverage(IssmDouble* pvalue){/*{{{*/
 
-	int        numnodes  = this->NumberofNodes();
+	int        numnodes  = this->NumberofNodes(this->interpolation_type);
 	IssmDouble numnodesd = reCast<int,IssmDouble>(numnodes);
 	IssmDouble value     = 0.;
@@ -179,5 +174,5 @@
 void PentaInput::SquareMin(IssmDouble* psquaremin,Parameters* parameters){/*{{{*/
 
-	int        numnodes=this->NumberofNodes();
+	int        numnodes=this->NumberofNodes(this->interpolation_type);
 	IssmDouble squaremin;
 
@@ -193,5 +188,5 @@
 void PentaInput::ConstrainMin(IssmDouble minimum){/*{{{*/
 
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(this->interpolation_type);
 	for(int i=0;i<numnodes;i++) if (values[i]<minimum) values[i]=minimum;
 }
@@ -201,5 +196,5 @@
 	/*Output*/
 	IssmDouble norm=0.;
-	int numnodes=this->NumberofNodes();
+	int numnodes=this->NumberofNodes(this->interpolation_type);
 
 	for(int i=0;i<numnodes;i++) if(fabs(values[i])>norm) norm=fabs(values[i]);
@@ -209,5 +204,5 @@
 IssmDouble PentaInput::Max(void){/*{{{*/
 
-	int  numnodes=this->NumberofNodes();
+	int  numnodes=this->NumberofNodes(this->interpolation_type);
 	IssmDouble max=values[0];
 
@@ -220,5 +215,5 @@
 IssmDouble PentaInput::MaxAbs(void){/*{{{*/
 
-	int  numnodes=this->NumberofNodes();
+	int  numnodes=this->NumberofNodes(this->interpolation_type);
 	IssmDouble max=fabs(values[0]);
 
@@ -231,5 +226,5 @@
 IssmDouble PentaInput::Min(void){/*{{{*/
 
-	const int  numnodes=this->NumberofNodes();
+	const int  numnodes=this->NumberofNodes(this->interpolation_type);
 	IssmDouble min=values[0];
 
@@ -242,5 +237,5 @@
 IssmDouble PentaInput::MinAbs(void){/*{{{*/
 
-	const int  numnodes=this->NumberofNodes();
+	const int  numnodes=this->NumberofNodes(this->interpolation_type);
 	IssmDouble min=fabs(values[0]);
 
@@ -253,5 +248,5 @@
 void PentaInput::Scale(IssmDouble scale_factor){/*{{{*/
 
-	const int numnodes=this->NumberofNodes();
+	const int numnodes=this->NumberofNodes(this->interpolation_type);
 	for(int i=0;i<numnodes;i++)values[i]=values[i]*scale_factor;
 }
@@ -259,5 +254,5 @@
 void PentaInput::AXPY(Input* xinput,IssmDouble scalar){/*{{{*/
 
-	const int numnodes=this->NumberofNodes();
+	const int numnodes=this->NumberofNodes(this->interpolation_type);
 	PentaInput* xpentainput=NULL;
 
@@ -271,5 +266,5 @@
 	  _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
 	xpentainput=(PentaInput*)xinput;
-	if(xpentainput->element_type!=this->element_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xpentainput->element_type));
+	if(xpentainput->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xpentainput->interpolation_type));
 
 	/*Carry out the AXPY operation depending on type:*/
@@ -281,5 +276,5 @@
 
 	int i;
-	const int numnodes=this->NumberofNodes();
+	const int numnodes=this->NumberofNodes(this->interpolation_type);
 
 	if(!xIsNan<IssmDouble>(cm_min)) for(i=0;i<numnodes;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
@@ -290,10 +285,10 @@
 void PentaInput::Extrude(void){/*{{{*/
 
-	switch(this->element_type){
+	switch(this->interpolation_type){
 		case P1Enum:
 			for(int i=0;i<3;i++) this->values[3+i]=this->values[i];
 			break;
 		default:
-			_error_("not supported yet for type "<<EnumToStringx(this->element_type));
+			_error_("not supported yet for type "<<EnumToStringx(this->interpolation_type));
 	}
 }
@@ -308,10 +303,10 @@
 
 	/*vertically integrate depending on type (and use P1 interpolation from now on)*/
-	switch(this->element_type){
+	switch(this->interpolation_type){
 		case P1Enum:
 		case P1bubbleEnum:
 		case P2Enum:
 			  {
-				this->element_type=P1Enum;
+				this->interpolation_type=P1Enum;
 				GaussPenta *gauss=new GaussPenta();
 				for(int iv=0;iv<3;iv++){
@@ -325,5 +320,5 @@
 			  }
 		default:
-			_error_("not supported yet for type "<<EnumToStringx(this->element_type));
+			_error_("not supported yet for type "<<EnumToStringx(this->interpolation_type));
 	}
 }
@@ -336,10 +331,10 @@
 	/*Intermediaries*/
 	PentaInput *xinputB  = NULL;
-	const int   numnodes = this->NumberofNodes();
+	const int   numnodes = this->NumberofNodes(this->interpolation_type);
 
 	/*Check that inputB is of the same type*/
 	if(inputB->ObjectEnum()!=PentaInputEnum)     _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
 	xinputB=(PentaInput*)inputB;
-	if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
+	if(xinputB->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->interpolation_type));
 
 	/*Allocate intermediary*/
@@ -353,5 +348,5 @@
 
 	/*Create new Penta vertex input (copy of current input)*/
-	outinput=new PentaInput(this->enum_type,AdotBvalues,this->element_type);
+	outinput=new PentaInput(this->enum_type,AdotBvalues,this->interpolation_type);
 
 	/*Return output pointer*/
@@ -369,5 +364,5 @@
 	int         i;
 	PentaInput  *xinputB   = NULL;
-	const int   numnodes  = this->NumberofNodes();
+	const int   numnodes  = this->NumberofNodes(this->interpolation_type);
 	IssmDouble *minvalues = xNew<IssmDouble>(numnodes);
 
@@ -375,5 +370,5 @@
 	if(inputB->ObjectEnum()!=PentaInputEnum)       _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
 	xinputB=(PentaInput*)inputB;
-	if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
+	if(xinputB->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->interpolation_type));
 
 	/*Create point wise min*/
@@ -384,5 +379,5 @@
 
 	/*Create new Penta vertex input (copy of current input)*/
-	outinput=new PentaInput(this->enum_type,&minvalues[0],this->element_type);
+	outinput=new PentaInput(this->enum_type,&minvalues[0],this->interpolation_type);
 
 	/*Return output pointer*/
@@ -399,5 +394,5 @@
 	int         i;
 	PentaInput  *xinputB   = NULL;
-	const int   numnodes  = this->NumberofNodes();
+	const int   numnodes  = this->NumberofNodes(this->interpolation_type);
 	IssmDouble *maxvalues = xNew<IssmDouble>(numnodes);
 
@@ -405,5 +400,5 @@
 	if(inputB->ObjectEnum()!=PentaInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
 	xinputB=(PentaInput*)inputB;
-	if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
+	if(xinputB->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->interpolation_type));
 
 	/*Create point wise max*/
@@ -414,5 +409,5 @@
 
 	/*Create new Penta vertex input (copy of current input)*/
-	outinput=new PentaInput(this->enum_type,&maxvalues[0],this->element_type);
+	outinput=new PentaInput(this->enum_type,&maxvalues[0],this->interpolation_type);
 
 	/*Return output pointer*/
Index: /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h	(revision 18078)
@@ -17,5 +17,6 @@
 
 	public:
-		int        enum_type;
+		int         enum_type;
+		int         interpolation_type;
 		IssmDouble* values;
 
Index: /issm/trunk-jpl/src/c/classes/Inputs/SegInput.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/SegInput.cpp	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Inputs/SegInput.cpp	(revision 18078)
@@ -17,18 +17,13 @@
 }
 /*}}}*/
-SegInput::SegInput(int in_enum_type,IssmDouble* in_values,int element_type_in)/*{{{*/
-	:SegRef(1)
-{
-
-	/*Set SegRef*/
-	this->SetElementType(element_type_in,0);
-	this->element_type=element_type_in;
+SegInput::SegInput(int in_enum_type,IssmDouble* in_values,int interpolation_type_in){/*{{{*/
 
 	/*Set Enum*/
 	enum_type=in_enum_type;
+	this->interpolation_type=interpolation_type_in;
 
 	/*Set values*/
-	this->values=xNew<IssmDouble>(this->NumberofNodes());
-	for(int i=0;i<this->NumberofNodes();i++) values[i]=in_values[i];
+	this->values=xNew<IssmDouble>(this->NumberofNodes(this->interpolation_type));
+	for(int i=0;i<this->NumberofNodes(this->interpolation_type);i++) values[i]=in_values[i];
 }
 /*}}}*/
@@ -46,5 +41,5 @@
 
 	_printf_(setw(15)<<"   SegInput "<<setw(25)<<left<<EnumToStringx(this->enum_type)<<" [");
-	for(int i=0;i<this->NumberofNodes();i++) _printf_(" "<<this->values[i]);
+	for(int i=0;i<this->NumberofNodes(this->interpolation_type);i++) _printf_(" "<<this->values[i]);
 	_printf_("]\n");
 }
@@ -60,5 +55,5 @@
 Object* SegInput::copy() {/*{{{*/
 
-	return new SegInput(this->enum_type,this->values,this->element_type);
+	return new SegInput(this->enum_type,this->values,this->interpolation_type);
 
 }
@@ -76,5 +71,5 @@
 void SegInput::GetInputAverage(IssmDouble* pvalue){/*{{{*/
 
-	int        numnodes  = this->NumberofNodes();
+	int        numnodes  = this->NumberofNodes(this->interpolation_type);
 	IssmDouble numnodesd = reCast<int,IssmDouble>(numnodes);
 	IssmDouble value     = 0.;
@@ -90,5 +85,5 @@
 	/*Call SegRef function*/
 	_assert_(gauss->Enum()==GaussSegEnum);
-	SegRef::GetInputValue(pvalue,&values[0],(GaussSeg*)gauss);
+	SegRef::GetInputValue(pvalue,&values[0],(GaussSeg*)gauss,this->interpolation_type);
 
 }
@@ -98,5 +93,5 @@
 	/*Call SegRef function*/
 	_assert_(gauss->Enum()==GaussSegEnum);
-	SegRef::GetInputDerivativeValue(p,&values[0],xyz_list,(GaussSeg*)gauss);
+	SegRef::GetInputDerivativeValue(p,&values[0],xyz_list,(GaussSeg*)gauss,this->interpolation_type);
 }
 /*}}}*/
@@ -107,5 +102,5 @@
 IssmDouble SegInput::Min(void){/*{{{*/
 
-	const int  numnodes=this->NumberofNodes();
+	const int  numnodes=this->NumberofNodes(this->interpolation_type);
 	IssmDouble min=values[0];
 
Index: /issm/trunk-jpl/src/c/classes/Inputs/SegInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/SegInput.h	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Inputs/SegInput.h	(revision 18078)
@@ -18,4 +18,5 @@
 	public:
 		int         enum_type;
+		int         interpolation_type;
 		IssmDouble* values;
 
Index: /issm/trunk-jpl/src/c/classes/Inputs/TetraInput.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/TetraInput.cpp	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Inputs/TetraInput.cpp	(revision 18078)
@@ -17,18 +17,13 @@
 }
 /*}}}*/
-TetraInput::TetraInput(int in_enum_type,IssmDouble* in_values,int element_type_in)/*{{{*/
-	:TetraRef(1)
-{
-
-	/*Set TetraRef*/
-	this->SetElementType(element_type_in,0);
-	this->element_type=element_type_in;
+TetraInput::TetraInput(int in_enum_type,IssmDouble* in_values,int interpolation_type_in){/*{{{*/
 
 	/*Set Enum*/
 	enum_type=in_enum_type;
+	this->interpolation_type=interpolation_type_in;
 
 	/*Set values*/
-	this->values=xNew<IssmDouble>(this->NumberofNodes());
-	for(int i=0;i<this->NumberofNodes();i++) values[i]=in_values[i];
+	this->values=xNew<IssmDouble>(this->NumberofNodes(this->interpolation_type));
+	for(int i=0;i<this->NumberofNodes(this->interpolation_type);i++) values[i]=in_values[i];
 }
 /*}}}*/
@@ -46,6 +41,6 @@
 
 	_printf_(setw(15)<<"   TetraInput "<<setw(25)<<left<<EnumToStringx(this->enum_type)<<" [");
-	for(int i=0;i<this->NumberofNodes();i++) _printf_(" "<<this->values[i]);
-	_printf_("] ("<<EnumToStringx(this->element_type)<<")\n");
+	for(int i=0;i<this->NumberofNodes(this->interpolation_type);i++) _printf_(" "<<this->values[i]);
+	_printf_("] ("<<EnumToStringx(this->interpolation_type)<<")\n");
 }
 /*}}}*/
@@ -60,5 +55,5 @@
 Object* TetraInput::copy() {/*{{{*/
 
-	return new TetraInput(this->enum_type,this->values,this->element_type);
+	return new TetraInput(this->enum_type,this->values,this->interpolation_type);
 
 }
@@ -74,5 +69,5 @@
 int  TetraInput::GetResultInterpolation(void){/*{{{*/
 
-	if(this->element_type==P0Enum){
+	if(this->interpolation_type==P0Enum){
 		return P0Enum;
 	}
@@ -83,5 +78,5 @@
 int  TetraInput::GetResultNumberOfNodes(void){/*{{{*/
 
-	return this->NumberofNodes();
+	return this->NumberofNodes(this->interpolation_type);
 
 }
@@ -89,5 +84,5 @@
 void TetraInput::ResultToPatch(IssmDouble* values,int nodesperelement,int sid){/*{{{*/
 
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(this->interpolation_type);
 
 	/*Some checks*/
@@ -105,5 +100,5 @@
 	/*Call TetraRef function*/
 	_assert_(gauss->Enum()==GaussTetraEnum);
-	TetraRef::GetInputValue(pvalue,&values[0],(GaussTetra*)gauss);
+	TetraRef::GetInputValue(pvalue,&values[0],(GaussTetra*)gauss,this->interpolation_type);
 
 }
@@ -113,5 +108,5 @@
 	/*Call TetraRef function*/
 	_assert_(gauss->Enum()==GaussTetraEnum);
-	TetraRef::GetInputDerivativeValue(p,&values[0],xyz_list,(GaussTetra*)gauss);
+	TetraRef::GetInputDerivativeValue(p,&values[0],xyz_list,(GaussTetra*)gauss,this->interpolation_type);
 }
 /*}}}*/
@@ -122,5 +117,5 @@
 void TetraInput::GetInputAverage(IssmDouble* pvalue){/*{{{*/
 
-	int        numnodes  = this->NumberofNodes();
+	int        numnodes  = this->NumberofNodes(this->interpolation_type);
 	IssmDouble numnodesd = reCast<int,IssmDouble>(numnodes);
 	IssmDouble value     = 0.;
@@ -175,5 +170,5 @@
 	TriaInput* outinput=NULL;
 
-	if(this->element_type==P0Enum){ 
+	if(this->interpolation_type==P0Enum){ 
 		outinput=new TriaInput(this->enum_type,&this->values[0],P0Enum);
 	}
@@ -204,5 +199,5 @@
 void TetraInput::SquareMin(IssmDouble* psquaremin,Parameters* parameters){/*{{{*/
 
-	int        numnodes=this->NumberofNodes();
+	int        numnodes=this->NumberofNodes(this->interpolation_type);
 	IssmDouble squaremin;
 
@@ -218,5 +213,5 @@
 void TetraInput::ConstrainMin(IssmDouble minimum){/*{{{*/
 
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(this->interpolation_type);
 	for(int i=0;i<numnodes;i++) if (values[i]<minimum) values[i]=minimum;
 }
@@ -226,5 +221,5 @@
 	/*Output*/
 	IssmDouble norm=0.;
-	int numnodes=this->NumberofNodes();
+	int numnodes=this->NumberofNodes(this->interpolation_type);
 
 	for(int i=0;i<numnodes;i++) if(fabs(values[i])>norm) norm=fabs(values[i]);
@@ -234,5 +229,5 @@
 IssmDouble TetraInput::Max(void){/*{{{*/
 
-	int  numnodes=this->NumberofNodes();
+	int  numnodes=this->NumberofNodes(this->interpolation_type);
 	IssmDouble max=values[0];
 
@@ -245,5 +240,5 @@
 IssmDouble TetraInput::MaxAbs(void){/*{{{*/
 
-	int  numnodes=this->NumberofNodes();
+	int  numnodes=this->NumberofNodes(this->interpolation_type);
 	IssmDouble max=fabs(values[0]);
 
@@ -256,5 +251,5 @@
 IssmDouble TetraInput::Min(void){/*{{{*/
 
-	const int  numnodes=this->NumberofNodes();
+	const int  numnodes=this->NumberofNodes(this->interpolation_type);
 	IssmDouble min=values[0];
 
@@ -267,5 +262,5 @@
 IssmDouble TetraInput::MinAbs(void){/*{{{*/
 
-	const int  numnodes=this->NumberofNodes();
+	const int  numnodes=this->NumberofNodes(this->interpolation_type);
 	IssmDouble min=fabs(values[0]);
 
@@ -278,5 +273,5 @@
 void TetraInput::Scale(IssmDouble scale_factor){/*{{{*/
 
-	const int numnodes=this->NumberofNodes();
+	const int numnodes=this->NumberofNodes(this->interpolation_type);
 	for(int i=0;i<numnodes;i++)values[i]=values[i]*scale_factor;
 }
@@ -284,5 +279,5 @@
 void TetraInput::Set(IssmDouble setvalue){/*{{{*/
 
-	const int numnodes=this->NumberofNodes();
+	const int numnodes=this->NumberofNodes(this->interpolation_type);
 	for(int i=0;i<numnodes;i++)values[i]=setvalue;
 }
@@ -290,5 +285,5 @@
 void TetraInput::AXPY(Input* xinput,IssmDouble scalar){/*{{{*/
 
-	const int numnodes=this->NumberofNodes();
+	const int numnodes=this->NumberofNodes(this->interpolation_type);
 	TetraInput*  xtriainput=NULL;
 
@@ -296,5 +291,5 @@
 	if(xinput->ObjectEnum()!=TetraInputEnum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
 	xtriainput=(TetraInput*)xinput;
-	if(xtriainput->element_type!=this->element_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
+	if(xtriainput->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
 
 	/*Carry out the AXPY operation depending on type:*/
@@ -306,5 +301,5 @@
 
 	int i;
-	const int numnodes=this->NumberofNodes();
+	const int numnodes=this->NumberofNodes(this->interpolation_type);
 
 	if(!xIsNan<IssmDouble>(cm_min)) for(i=0;i<numnodes;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
@@ -325,5 +320,5 @@
 	int         i;
 	TetraInput  *xinputB   = NULL;
-	const int   numnodes  = this->NumberofNodes();
+	const int   numnodes  = this->NumberofNodes(this->interpolation_type);
 	IssmDouble *minvalues = xNew<IssmDouble>(numnodes);
 
@@ -331,5 +326,5 @@
 	if(inputB->ObjectEnum()!=TetraInputEnum)       _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
 	xinputB=(TetraInput*)inputB;
-	if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
+	if(xinputB->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->interpolation_type));
 
 	/*Create point wise min*/
@@ -340,5 +335,5 @@
 
 	/*Create new Tetra vertex input (copy of current input)*/
-	outinput=new TetraInput(this->enum_type,&minvalues[0],this->element_type);
+	outinput=new TetraInput(this->enum_type,&minvalues[0],this->interpolation_type);
 
 	/*Return output pointer*/
@@ -356,5 +351,5 @@
 	int         i;
 	TetraInput  *xinputB   = NULL;
-	const int   numnodes  = this->NumberofNodes();
+	const int   numnodes  = this->NumberofNodes(this->interpolation_type);
 	IssmDouble *maxvalues = xNew<IssmDouble>(numnodes);
 
@@ -362,5 +357,5 @@
 	if(inputB->ObjectEnum()!=TetraInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
 	xinputB=(TetraInput*)inputB;
-	if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
+	if(xinputB->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->interpolation_type));
 
 	/*Create point wise max*/
@@ -371,5 +366,5 @@
 
 	/*Create new Tetra vertex input (copy of current input)*/
-	outinput=new TetraInput(this->enum_type,&maxvalues[0],this->element_type);
+	outinput=new TetraInput(this->enum_type,&maxvalues[0],this->interpolation_type);
 
 	/*Return output pointer*/
@@ -386,10 +381,10 @@
 	/*Intermediaries*/
 	TetraInput *xinputB  = NULL;
-	const int   numnodes = this->NumberofNodes();
+	const int   numnodes = this->NumberofNodes(this->interpolation_type);
 
 	/*Check that inputB is of the same type*/
 	if(inputB->ObjectEnum()!=TetraInputEnum)     _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
 	xinputB=(TetraInput*)inputB;
-	if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
+	if(xinputB->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->interpolation_type));
 
 	/*Allocate intermediary*/
@@ -403,5 +398,5 @@
 
 	/*Create new Tetra vertex input (copy of current input)*/
-	outinput=new TetraInput(this->enum_type,AdotBvalues,this->element_type);
+	outinput=new TetraInput(this->enum_type,AdotBvalues,this->interpolation_type);
 
 	/*Return output pointer*/
Index: /issm/trunk-jpl/src/c/classes/Inputs/TetraInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/TetraInput.h	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Inputs/TetraInput.h	(revision 18078)
@@ -18,4 +18,5 @@
 	public:
 		int         enum_type;
+		int         interpolation_type;
 		IssmDouble* values;
 
Index: /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp	(revision 18078)
@@ -17,18 +17,13 @@
 }
 /*}}}*/
-TriaInput::TriaInput(int in_enum_type,IssmDouble* in_values,int element_type_in)/*{{{*/
-	:TriaRef(1)
-{
-
-	/*Set TriaRef*/
-	this->SetElementType(element_type_in,0);
-	this->element_type=element_type_in;
+TriaInput::TriaInput(int in_enum_type,IssmDouble* in_values,int interpolation_type_in){/*{{{*/
 
 	/*Set Enum*/
 	enum_type=in_enum_type;
+	this->interpolation_type=interpolation_type_in;
 
 	/*Set values*/
-	this->values=xNew<IssmDouble>(this->NumberofNodes());
-	for(int i=0;i<this->NumberofNodes();i++) values[i]=in_values[i];
+	this->values=xNew<IssmDouble>(this->NumberofNodes(this->interpolation_type));
+	for(int i=0;i<this->NumberofNodes(this->interpolation_type);i++) values[i]=in_values[i];
 }
 /*}}}*/
@@ -46,6 +41,6 @@
 
 	_printf_(setw(15)<<"   TriaInput "<<setw(25)<<left<<EnumToStringx(this->enum_type)<<" [");
-	for(int i=0;i<this->NumberofNodes();i++) _printf_(" "<<this->values[i]);
-	_printf_("] ("<<EnumToStringx(this->element_type)<<")\n");
+	for(int i=0;i<this->NumberofNodes(this->interpolation_type);i++) _printf_(" "<<this->values[i]);
+	_printf_("] ("<<EnumToStringx(this->interpolation_type)<<")\n");
 }
 /*}}}*/
@@ -60,5 +55,5 @@
 Object* TriaInput::copy() {/*{{{*/
 
-	return new TriaInput(this->enum_type,this->values,this->element_type);
+	return new TriaInput(this->enum_type,this->values,this->interpolation_type);
 
 }
@@ -78,5 +73,5 @@
 
 	/*Create new Tria input (copy of current input)*/
-	outinput=new TriaInput(this->enum_type,&this->values[0],this->element_type);
+	outinput=new TriaInput(this->enum_type,&this->values[0],this->interpolation_type);
 
 	/*Assign output*/
@@ -90,5 +85,5 @@
 	SegInput* outinput=NULL;
 
-	if(this->element_type==P0Enum){ 
+	if(this->interpolation_type==P0Enum){ 
 		outinput=new SegInput(this->enum_type,&this->values[0],P0Enum);
 	}
@@ -112,5 +107,5 @@
 int  TriaInput::GetResultInterpolation(void){/*{{{*/
 
-	if(this->element_type==P0Enum){
+	if(this->interpolation_type==P0Enum){
 		return P0Enum;
 	}
@@ -121,5 +116,5 @@
 int  TriaInput::GetResultNumberOfNodes(void){/*{{{*/
 
-	return this->NumberofNodes();
+	return this->NumberofNodes(this->interpolation_type);
 
 }
@@ -127,5 +122,5 @@
 void TriaInput::ResultToPatch(IssmDouble* values,int nodesperelement,int sid){/*{{{*/
 
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(this->interpolation_type);
 
 	/*Some checks*/
@@ -143,5 +138,5 @@
 	/*Call TriaRef function*/
 	_assert_(gauss->Enum()==GaussTriaEnum);
-	TriaRef::GetInputValue(pvalue,&values[0],(GaussTria*)gauss);
+	TriaRef::GetInputValue(pvalue,&values[0],(GaussTria*)gauss,this->interpolation_type);
 
 }
@@ -151,5 +146,5 @@
 	/*Call TriaRef function*/
 	_assert_(gauss->Enum()==GaussTriaEnum);
-	TriaRef::GetInputDerivativeValue(p,&values[0],xyz_list,(GaussTria*)gauss);
+	TriaRef::GetInputDerivativeValue(p,&values[0],xyz_list,(GaussTria*)gauss,this->interpolation_type);
 }
 /*}}}*/
@@ -160,5 +155,5 @@
 void TriaInput::GetInputAverage(IssmDouble* pvalue){/*{{{*/
 
-	int        numnodes  = this->NumberofNodes();
+	int        numnodes  = this->NumberofNodes(this->interpolation_type);
 	IssmDouble numnodesd = reCast<int,IssmDouble>(numnodes);
 	IssmDouble value     = 0.;
@@ -212,5 +207,5 @@
 void TriaInput::SquareMin(IssmDouble* psquaremin,Parameters* parameters){/*{{{*/
 
-	int        numnodes=this->NumberofNodes();
+	int        numnodes=this->NumberofNodes(this->interpolation_type);
 	IssmDouble squaremin;
 
@@ -226,5 +221,5 @@
 void TriaInput::ConstrainMin(IssmDouble minimum){/*{{{*/
 
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(this->interpolation_type);
 	for(int i=0;i<numnodes;i++) if (values[i]<minimum) values[i]=minimum;
 }
@@ -234,5 +229,5 @@
 	/*Output*/
 	IssmDouble norm=0.;
-	int numnodes=this->NumberofNodes();
+	int numnodes=this->NumberofNodes(this->interpolation_type);
 
 	for(int i=0;i<numnodes;i++) if(fabs(values[i])>norm) norm=fabs(values[i]);
@@ -242,5 +237,5 @@
 IssmDouble TriaInput::Max(void){/*{{{*/
 
-	int  numnodes=this->NumberofNodes();
+	int  numnodes=this->NumberofNodes(this->interpolation_type);
 	IssmDouble max=values[0];
 
@@ -253,5 +248,5 @@
 IssmDouble TriaInput::MaxAbs(void){/*{{{*/
 
-	int  numnodes=this->NumberofNodes();
+	int  numnodes=this->NumberofNodes(this->interpolation_type);
 	IssmDouble max=fabs(values[0]);
 
@@ -264,5 +259,5 @@
 IssmDouble TriaInput::Min(void){/*{{{*/
 
-	const int  numnodes=this->NumberofNodes();
+	const int  numnodes=this->NumberofNodes(this->interpolation_type);
 	IssmDouble min=values[0];
 
@@ -275,5 +270,5 @@
 IssmDouble TriaInput::MinAbs(void){/*{{{*/
 
-	const int  numnodes=this->NumberofNodes();
+	const int  numnodes=this->NumberofNodes(this->interpolation_type);
 	IssmDouble min=fabs(values[0]);
 
@@ -286,5 +281,5 @@
 void TriaInput::Scale(IssmDouble scale_factor){/*{{{*/
 
-	const int numnodes=this->NumberofNodes();
+	const int numnodes=this->NumberofNodes(this->interpolation_type);
 	for(int i=0;i<numnodes;i++)values[i]=values[i]*scale_factor;
 }
@@ -292,5 +287,5 @@
 void TriaInput::Set(IssmDouble setvalue){/*{{{*/
 
-	const int numnodes=this->NumberofNodes();
+	const int numnodes=this->NumberofNodes(this->interpolation_type);
 	for(int i=0;i<numnodes;i++)values[i]=setvalue;
 }
@@ -298,5 +293,5 @@
 void TriaInput::AXPY(Input* xinput,IssmDouble scalar){/*{{{*/
 
-	const int numnodes=this->NumberofNodes();
+	const int numnodes=this->NumberofNodes(this->interpolation_type);
 	TriaInput*  xtriainput=NULL;
 
@@ -304,5 +299,5 @@
 	if(xinput->ObjectEnum()!=TriaInputEnum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
 	xtriainput=(TriaInput*)xinput;
-	if(xtriainput->element_type!=this->element_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
+	if(xtriainput->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
 
 	/*Carry out the AXPY operation depending on type:*/
@@ -314,5 +309,5 @@
 
 	int i;
-	const int numnodes=this->NumberofNodes();
+	const int numnodes=this->NumberofNodes(this->interpolation_type);
 
 	if(!xIsNan<IssmDouble>(cm_min)) for(i=0;i<numnodes;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
@@ -333,5 +328,5 @@
 	int         i;
 	TriaInput  *xinputB   = NULL;
-	const int   numnodes  = this->NumberofNodes();
+	const int   numnodes  = this->NumberofNodes(this->interpolation_type);
 	IssmDouble *minvalues = xNew<IssmDouble>(numnodes);
 
@@ -339,5 +334,5 @@
 	if(inputB->ObjectEnum()!=TriaInputEnum)       _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
 	xinputB=(TriaInput*)inputB;
-	if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
+	if(xinputB->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->interpolation_type));
 
 	/*Create point wise min*/
@@ -348,5 +343,5 @@
 
 	/*Create new Tria vertex input (copy of current input)*/
-	outinput=new TriaInput(this->enum_type,&minvalues[0],this->element_type);
+	outinput=new TriaInput(this->enum_type,&minvalues[0],this->interpolation_type);
 
 	/*Return output pointer*/
@@ -364,5 +359,5 @@
 	int         i;
 	TriaInput  *xinputB   = NULL;
-	const int   numnodes  = this->NumberofNodes();
+	const int   numnodes  = this->NumberofNodes(this->interpolation_type);
 	IssmDouble *maxvalues = xNew<IssmDouble>(numnodes);
 
@@ -370,5 +365,5 @@
 	if(inputB->ObjectEnum()!=TriaInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
 	xinputB=(TriaInput*)inputB;
-	if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
+	if(xinputB->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->interpolation_type));
 
 	/*Create point wise max*/
@@ -379,5 +374,5 @@
 
 	/*Create new Tria vertex input (copy of current input)*/
-	outinput=new TriaInput(this->enum_type,&maxvalues[0],this->element_type);
+	outinput=new TriaInput(this->enum_type,&maxvalues[0],this->interpolation_type);
 
 	/*Return output pointer*/
@@ -394,10 +389,10 @@
 	/*Intermediaries*/
 	TriaInput *xinputB  = NULL;
-	const int   numnodes = this->NumberofNodes();
+	const int   numnodes = this->NumberofNodes(this->interpolation_type);
 
 	/*Check that inputB is of the same type*/
 	if(inputB->ObjectEnum()!=TriaInputEnum)     _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
 	xinputB=(TriaInput*)inputB;
-	if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
+	if(xinputB->interpolation_type!=this->interpolation_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->interpolation_type));
 
 	/*Allocate intermediary*/
@@ -411,5 +406,5 @@
 
 	/*Create new Tria vertex input (copy of current input)*/
-	outinput=new TriaInput(this->enum_type,AdotBvalues,this->element_type);
+	outinput=new TriaInput(this->enum_type,AdotBvalues,this->interpolation_type);
 
 	/*Return output pointer*/
Index: /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h	(revision 18078)
@@ -18,4 +18,5 @@
 	public:
 		int         enum_type;
+		int         interpolation_type;
 		IssmDouble* values;
 
Index: /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp	(revision 18077)
+++ /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp	(revision 18078)
@@ -451,6 +451,6 @@
 		gauss->GaussPoint(ig);
 
-		tria->GetSegmentBFlux(&B[0],gauss,index1,index2);
-		tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2);
+		tria->GetSegmentBFlux(&B[0],gauss,index1,index2,tria->FiniteElement());
+		tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2,tria->FiniteElement());
 
 		vxaverage_input->GetInputValue(&vx,gauss);
@@ -529,5 +529,5 @@
 		gauss->GaussPoint(ig);
 
-		tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);
+		tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement());
 
 		vxaverage_input->GetInputValue(&vx,gauss);
@@ -597,6 +597,6 @@
 		gauss->GaussPoint(ig);
 
-		tria->GetSegmentBFlux(&B[0],gauss,index1,index2);
-		tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2);
+		tria->GetSegmentBFlux(&B[0],gauss,index1,index2,tria->FiniteElement());
+		tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2,tria->FiniteElement());
 
 		vxaverage_input->GetInputValue(&vx,gauss);
@@ -674,5 +674,5 @@
 		gauss->GaussPoint(ig);
 
-		tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);
+		tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement());
 
 		vxaverage_input->GetInputValue(&vx,gauss);
@@ -790,5 +790,5 @@
 		gauss->GaussPoint(ig);
 
-		tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);
+		tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement());
 
 		vxaverage_input->GetInputValue(&vx,gauss);
@@ -876,5 +876,5 @@
 		gauss->GaussPoint(ig);
 
-		tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);
+		tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement());
 
 		vxaverage_input->GetInputValue(&vx,gauss);