Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15325)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15326)
@@ -986,4 +986,142 @@
 }
 /*}}}*/
+/*FUNCTION PentaRef::GetNodalFunctions{{{*/
+void PentaRef::GetNodalFunctions(IssmDouble* basis,GaussPenta* gauss){
+	/*This routine returns the values of the nodal functions  at the gaussian point.*/
+
+	_assert_(basis);
+
+	switch(this->element_type){
+		case P1Enum:
+		case P1DGEnum:
+			basis[0]=gauss->coord1*(1.-gauss->coord4)/2.;
+			basis[1]=gauss->coord2*(1.-gauss->coord4)/2.;
+			basis[2]=gauss->coord3*(1.-gauss->coord4)/2.;
+			basis[3]=gauss->coord1*(1.+gauss->coord4)/2.;
+			basis[4]=gauss->coord2*(1.+gauss->coord4)/2.;
+			basis[5]=gauss->coord3*(1.+gauss->coord4)/2.;
+			return;
+		case MINIEnum:
+			basis[0]=gauss->coord1*(1.-gauss->coord4)/2.;
+			basis[1]=gauss->coord2*(1.-gauss->coord4)/2.;
+			basis[2]=gauss->coord3*(1.-gauss->coord4)/2.;
+			basis[3]=gauss->coord1*(1.+gauss->coord4)/2.;
+			basis[4]=gauss->coord2*(1.+gauss->coord4)/2.;
+			basis[5]=gauss->coord3*(1.+gauss->coord4)/2.;
+			basis[6]=27.*gauss->coord1*gauss->coord2*gauss->coord3*(1.+gauss->coord4)*(1.-gauss->coord4);
+			return;
+		default:
+			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+	}
+}
+/*}}}*/
+/*FUNCTION PentaRef::GetNodalFunctionsDerivatives{{{*/
+void PentaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussPenta* gauss){
+
+	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
+	 * actual coordinate system): */
+	IssmDouble    Jinv[3][3];
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Get nodal functions derivatives in reference triangle*/
+	IssmDouble* dbasis_ref=xNew<IssmDouble>(3*numnodes);
+	GetNodalFunctionsDerivativesReference(dbasis_ref,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<numnodes;i++){
+		dbasis[numnodes*0+i]=Jinv[0][0]*dbasis_ref[0*numnodes+i]+Jinv[0][1]*dbasis_ref[1*numnodes+i]+Jinv[0][2]*dbasis_ref[2*numnodes+i];
+		dbasis[numnodes*1+i]=Jinv[1][0]*dbasis_ref[0*numnodes+i]+Jinv[1][1]*dbasis_ref[1*numnodes+i]+Jinv[1][2]*dbasis_ref[2*numnodes+i];
+		dbasis[numnodes*2+i]=Jinv[2][0]*dbasis_ref[0*numnodes+i]+Jinv[2][1]*dbasis_ref[1*numnodes+i]+Jinv[2][2]*dbasis_ref[2*numnodes+i];
+	}
+
+}
+/*}}}*/
+/*FUNCTION PentaRef::GetNodalFunctionsDerivativesReference{{{*/
+void PentaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussPenta* gauss){
+
+	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
+	 * natural coordinate system) at the gaussian point. */
+
+	_assert_(dbasis && gauss);
+
+	/*Get current coordinates in reference element*/
+	IssmDouble r=gauss->coord2-gauss->coord1;
+	IssmDouble s=-3.0/SQRT3*(gauss->coord1+gauss->coord2-2.0/3.0);
+	IssmDouble zeta=gauss->coord4;
+
+	switch(this->element_type){
+		case P1Enum: case P1DGEnum:
+			/*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*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
+			/*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*(0.5*r-SQRT3/6.0*s+ONETHIRD);
+			/*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*(SQRT3/3.0*s+ONETHIRD);
+			/*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*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
+			/*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*(0.5*r-SQRT3/6.0*s+ONETHIRD);
+			/*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*(SQRT3/3.0*s+ONETHIRD);
+		case MINIEnum:
+			/*Nodal function 1*/
+			dbasis[NUMNODESMINI*0+0]=-0.5*(1.0-zeta)/2.0;
+			dbasis[NUMNODESMINI*1+0]=-SQRT3/6.0*(1.0-zeta)/2.0;
+			dbasis[NUMNODESMINI*2+0]=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
+			/*Nodal function 2*/
+			dbasis[NUMNODESMINI*0+1]=0.5*(1.0-zeta)/2.0;
+			dbasis[NUMNODESMINI*1+1]=-SQRT3/6.0*(1.0-zeta)/2.0;
+			dbasis[NUMNODESMINI*2+1]=-0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
+			/*Nodal function 3*/
+			dbasis[NUMNODESMINI*0+2]=0;
+			dbasis[NUMNODESMINI*1+2]=SQRT3/3.0*(1.0-zeta)/2.0;
+			dbasis[NUMNODESMINI*2+2]=-0.5*(SQRT3/3.0*s+ONETHIRD);
+			/*Nodal function 4*/
+			dbasis[NUMNODESMINI*0+3]=-0.5*(1.0+zeta)/2.0;
+			dbasis[NUMNODESMINI*1+3]=-SQRT3/6.0*(1.0+zeta)/2.0;
+			dbasis[NUMNODESMINI*2+3]=0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
+			/*Nodal function 5*/
+			dbasis[NUMNODESMINI*0+4]=0.5*(1.0+zeta)/2.0;
+			dbasis[NUMNODESMINI*1+4]=-SQRT3/6.0*(1.0+zeta)/2.0;
+			dbasis[NUMNODESMINI*2+4]=0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
+			/*Nodal function 6*/
+			dbasis[NUMNODESMINI*0+5]=0;
+			dbasis[NUMNODESMINI*1+5]=SQRT3/3.0*(1.0+zeta)/2.0;
+			dbasis[NUMNODESMINI*2+5]=0.5*(SQRT3/3.0*s+ONETHIRD);
+			/*Nodal function 7*/
+			dbasis[NUMNODESMINI*0+6]=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(SQRT3*s+1.0);
+			dbasis[NUMNODESMINI*1+6]=9.0/4.0*(1+zeta)*(1-zeta)*(SQRT3*pow(s,2.0)-2.0*s-SQRT3*pow(r,2.0));
+			dbasis[NUMNODESMINI*2+6]=27*gauss->coord1*gauss->coord2*gauss->coord3*(-2.0*zeta);
+			return;
+		default:
+			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+	}
+
+
+}
+/*}}}*/
 /*FUNCTION PentaRef::GetNodalFunctionsMINI{{{*/
 void PentaRef::GetNodalFunctionsMINI(IssmDouble* l1l7, GaussPenta* gauss){
@@ -1136,4 +1274,7 @@
 	*(dl1dl6+NUMNODESP1*1+0)=-0.5/SQRT3*(1.0-z)/2.0;
 	*(dl1dl6+NUMNODESP1*2+0)=-0.5*A1;
+	//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*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
 
 	/*Second nodal function: The corresponding nodal function is N=A2*(1-z)/2. Its derivatives follow*/
@@ -1250,4 +1391,5 @@
 	switch(this->element_type){
 		case P1Enum:   return NUMNODESP1;
+		case MINIEnum: return NUMNODESMINI;
 		default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
 	}
Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 15325)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 15326)
@@ -22,4 +22,7 @@
 
 		/*Numerics*/
+		void GetNodalFunctions(IssmDouble* basis, GaussPenta* gauss);
+		void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,GaussPenta* gauss);
+		void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussPenta* gauss);
 		void GetNodalFunctionsP1(IssmDouble* l1l6, GaussPenta* gauss);
 		void GetNodalFunctionsMINI(IssmDouble* l1l7, GaussPenta* gauss);
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 15325)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 15326)
@@ -491,5 +491,5 @@
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
 	 * actual coordinate system): */
-	IssmDouble    Jinv[NDOF2][NDOF2];
+	IssmDouble    Jinv[2][2];
 
 	/*Fetch number of nodes for this finite element*/
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 15325)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 15326)
@@ -27,6 +27,6 @@
 		void GetBprimeMacAyeal(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
 		void GetBprimeMacAyealStokes(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
-		void GetBprimePrognostic(IssmDouble* Bprime_prog, IssmDouble* xyz_list, GaussTria* gauss);
-		void GetBPrognostic(IssmDouble* B_prog, IssmDouble* xyz_list, GaussTria* gauss);
+		void GetBprimePrognostic(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
+		void GetBPrognostic(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss);
 		void GetBHydro(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss);
 		void GetBMacAyealFriction(IssmDouble* L, IssmDouble* xyz_list,GaussTria* gauss);
@@ -39,6 +39,6 @@
 		void GetSegmentBFlux(IssmDouble* B,GaussTria* gauss, int index1,int index2);
 		void GetSegmentBprimeFlux(IssmDouble* Bprime,GaussTria* gauss, int index1,int index2);
-		void GetNodalFunctionsDerivatives(IssmDouble* basis,IssmDouble* xyz_list, GaussTria* gauss);
-		void GetNodalFunctionsDerivativesReference(IssmDouble* dl1dl3,GaussTria* gauss);
+		void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss);
+		void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTria* gauss);
 		void GetInputValue(IssmDouble* pp, IssmDouble* plist, GaussTria* gauss);
 		void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss);