Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 17224)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 17225)
@@ -178,6 +178,15 @@
 }
 /*}}}*/
-/*FUNCTION PentaRef::GetNodalFunctions{{{*/
-void PentaRef::GetNodalFunctions(IssmDouble* basis,Gauss* gauss_in){
+/*FUNCTION PentaRef::GetNodalFunctions(IssmDouble* basis,Gauss* gauss){{{*/
+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);
+
+}
+/*}}}*/
+/*FUNCTION PentaRef::GetNodalFunctions(IssmDouble* basis,Gauss* gauss_in,int finiteelement){{{*/
+void PentaRef::GetNodalFunctions(IssmDouble* basis,Gauss* gauss_in,int finiteelement){
 	/*This routine returns the values of the nodal functions  at the gaussian point.*/
 
@@ -191,5 +200,8 @@
 	IssmDouble zeta=gauss->coord4;
 
-	switch(this->element_type){
+	switch(finiteelement){
+		case P0Enum: 
+			basis[0]=1.;
+			return;
 		case P1Enum: case P1DGEnum:
 			basis[0]=gauss->coord1*(1.-zeta)/2.;
@@ -408,4 +420,10 @@
 
 	switch(this->element_type){
+		case P0Enum: 
+			/*Zero derivative*/
+			dbasis[NUMNODESP0*0+0]   = 0.;
+			dbasis[NUMNODESP0*1+0]   = 0.;
+			dbasis[NUMNODESP0*2+0]   = 0.;
+			return;
 		case P1Enum: case P1DGEnum:
 			/*Nodal function 1*/
@@ -996,16 +1014,30 @@
 }
 /*}}}*/
-/*FUNCTION PentaRef::GetInputValue{{{*/
+/*FUNCTION PentaRef::GetInputValue(IssmDouble* pvalue,IssmDouble* plist,Gauss* gauss){{{*/
 void PentaRef::GetInputValue(IssmDouble* pvalue,IssmDouble* plist,Gauss* gauss){
-	/*P1 interpolation on Gauss point*/
-
-	/*intermediary*/
-	IssmDouble basis[6];
-
-	/*nodal functions: */
-	GetNodalFunctionsP1(&basis[0],gauss);
-
-	/*Assign output pointers:*/
-	*pvalue=basis[0]*plist[0]+basis[1]*plist[1]+basis[2]*plist[2]+basis[3]*plist[3]+basis[4]*plist[4]+basis[5]*plist[5];
+
+	GetInputValue(pvalue,plist,gauss,this->element_type);
+
+}
+/*}}}*/
+/*FUNCTION PentaRef::GetInputValue(IssmDouble* pvalue,IssmDouble* plist,Gauss* gauss,int finiteelement){{{*/
+void PentaRef::GetInputValue(IssmDouble* pvalue,IssmDouble* plist,Gauss* gauss,int finiteelement){
+
+	/*Output*/
+	IssmDouble value =0.;
+
+	/*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);
+
+	/*Calculate parameter for this Gauss point*/
+	for(int i=0;i<numnodes;i++) value += basis[i]*plist[i];
+
+	/*Assign output pointer*/
+	xDelete<IssmDouble>(basis);
+	*pvalue = value;
 
 }
@@ -1022,20 +1054,40 @@
 	 *   p is a vector of size 3x1 already allocated.
 	 */
-	IssmDouble dbasis[3][NUMNODESP1];
-
-	/*Get nodal funnctions derivatives in actual coordinate system: */
-	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
-
-	/*Assign output*/
-	p[0]=plist[0]*dbasis[0][0]+plist[1]*dbasis[0][1]+plist[2]*dbasis[0][2]+plist[3]*dbasis[0][3]+plist[4]*dbasis[0][4]+plist[5]*dbasis[0][5];
-	p[1]=plist[0]*dbasis[1][0]+plist[1]*dbasis[1][1]+plist[2]*dbasis[1][2]+plist[3]*dbasis[1][3]+plist[4]*dbasis[1][4]+plist[5]*dbasis[1][5];
-	p[2]=plist[0]*dbasis[2][0]+plist[1]*dbasis[2][1]+plist[2]*dbasis[2][2]+plist[3]*dbasis[2][3]+plist[4]*dbasis[2][4]+plist[5]*dbasis[2][5];
-
-}
-/*}}}*/
-/*FUNCTION PentaRef::NumberofNodes{{{*/
+
+	/*Output*/
+	IssmDouble dpx=0.;
+	IssmDouble dpy=0.;
+	IssmDouble dpz=0.;
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
+	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+	/*Calculate parameter for this Gauss point*/
+	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*/
+	xDelete<IssmDouble>(dbasis);
+	p[0]=dpx;
+	p[1]=dpy;
+	p[2]=dpz;
+
+}
+/*}}}*/
+/*FUNCTION PentaRef::NumberofNodes(){{{*/
 int PentaRef::NumberofNodes(void){
 
-	switch(this->element_type){
+	return this->NumberofNodes(this->element_type);
+}
+/*}}}*/
+/*FUNCTION PentaRef::NumberofNodes(int finiteelement){{{*/
+int PentaRef::NumberofNodes(int finiteelement){
+
+	switch(finiteelement){
 		case P0Enum:                return NUMNODESP0;
 		case P1Enum:                return NUMNODESP1;
Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 17224)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 17225)
@@ -23,4 +23,5 @@
 		/*Numerics*/
 		void GetNodalFunctions(IssmDouble* basis, Gauss* gauss);
+		void GetNodalFunctions(IssmDouble* basis, Gauss* gauss,int finiteelement);
 		void GetNodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss);
 		void GetNodalFunctionsPressure(IssmDouble* basis, Gauss* gauss);
@@ -43,4 +44,5 @@
 		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);
 
@@ -48,4 +50,5 @@
 		void SurfaceNodeIndices(int* pnumindices,int** pindices,int finiteelement);
 		int  NumberofNodes(void);
+		int  NumberofNodes(int finiteelement);
 		int  NumberofNodesVelocity(void);
 		int  NumberofNodesPressure(void);
Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 17224)
+++ /issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 17225)
@@ -1482,5 +1482,5 @@
 	for(int i=0;i<elements->Size();i++){
 		Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
-		if(!doublearray) element->InputUpdateFromConstant(default_value,vector_enum); 
+		if(!doublearray) element->AddInput(vector_enum,&default_value,P0Enum); 
 		else             element->InputCreate(doublearray,this,M,N,vector_layout,vector_enum,code);//we need i to index into elements.
 	}