Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15457)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15458)
@@ -1195,6 +1195,4 @@
 
 	/*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;
 
@@ -1204,53 +1202,53 @@
 			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);
+			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*(0.5*r-SQRT3/6.0*s+ONETHIRD);
+			dbasis[NUMNODESP1*2+1]=-0.5*gauss->coord2;
 			/*Nodal function 3*/
-			dbasis[NUMNODESP1*0+2]=0;
+			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);
+			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*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
+			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*(0.5*r-SQRT3/6.0*s+ONETHIRD);
+			dbasis[NUMNODESP1*2+4]=0.5*gauss->coord2;
 			/*Nodal function 6*/
-			dbasis[NUMNODESP1*0+5]=0;
+			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);
+			dbasis[NUMNODESP1*2+5]=0.5*gauss->coord3;
 		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);
+			dbasis[NUMNODESMINI*2+0]=-0.5*gauss->coord1;
 			/*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);
+			dbasis[NUMNODESMINI*2+1]=-0.5*gauss->coord2;
 			/*Nodal function 3*/
-			dbasis[NUMNODESMINI*0+2]=0;
+			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);
+			dbasis[NUMNODESMINI*2+2]=-0.5*gauss->coord3;
 			/*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);
+			dbasis[NUMNODESMINI*2+3]=0.5*gauss->coord1;
 			/*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);
+			dbasis[NUMNODESMINI*2+4]=0.5*gauss->coord2;
 			/*Nodal function 6*/
-			dbasis[NUMNODESMINI*0+5]=0;
+			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);
+			dbasis[NUMNODESMINI*2+5]=0.5*gauss->coord3;
 			/*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*0+6]=27.*(1.+zeta)*(1.-zeta)*(-.5*gauss->coord2*gauss->coord3 + .5*gauss->coord1*gauss->coord3);
+			dbasis[NUMNODESMINI*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[NUMNODESMINI*2+6]=27*gauss->coord1*gauss->coord2*gauss->coord3*(-2.0*zeta);
 			return;
@@ -1282,5 +1280,4 @@
 	 * actual coordinate system): */
 
-	int       i;
 	IssmDouble    dh1dh7_ref[3][NUMNODESMINI];
 	IssmDouble    Jinv[3][3];
@@ -1299,5 +1296,5 @@
 	 */
 
-	for (i=0;i<NUMNODESMINI;i++){
+	for(int i=0;i<NUMNODESMINI;i++){
 		*(dh1dh7+NUMNODESMINI*0+i)=Jinv[0][0]*dh1dh7_ref[0][i]+Jinv[0][1]*dh1dh7_ref[1][i]+Jinv[0][2]*dh1dh7_ref[2][i];
 		*(dh1dh7+NUMNODESMINI*1+i)=Jinv[1][0]*dh1dh7_ref[0][i]+Jinv[1][1]*dh1dh7_ref[1][i]+Jinv[1][2]*dh1dh7_ref[2][i];
@@ -1308,47 +1305,38 @@
 /*}}}*/
 /*FUNCTION PentaRef::GetNodalFunctionsMINIDerivativesReference{{{*/
-void PentaRef::GetNodalFunctionsMINIDerivativesReference(IssmDouble* dl1dl7,GaussPenta* gauss){
+void PentaRef::GetNodalFunctionsMINIDerivativesReference(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. */
-	IssmDouble r=gauss->coord2-gauss->coord1;
-	IssmDouble s=-3.0/SQRT3*(gauss->coord1+gauss->coord2-2.0/3.0);
 	IssmDouble zeta=gauss->coord4;
 
-	/*First nodal function: */
-	*(dl1dl7+NUMNODESMINI*0+0)=-0.5*(1.0-zeta)/2.0;
-	*(dl1dl7+NUMNODESMINI*1+0)=-SQRT3/6.0*(1.0-zeta)/2.0;
-	*(dl1dl7+NUMNODESMINI*2+0)=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
-
-	/*Second nodal function: */
-	*(dl1dl7+NUMNODESMINI*0+1)=0.5*(1.0-zeta)/2.0;
-	*(dl1dl7+NUMNODESMINI*1+1)=-SQRT3/6.0*(1.0-zeta)/2.0;
-	*(dl1dl7+NUMNODESMINI*2+1)=-0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
-
-	/*Third nodal function: */
-	*(dl1dl7+NUMNODESMINI*0+2)=0;
-	*(dl1dl7+NUMNODESMINI*1+2)=SQRT3/3.0*(1.0-zeta)/2.0;
-	*(dl1dl7+NUMNODESMINI*2+2)=-0.5*(SQRT3/3.0*s+ONETHIRD);
-
-	/*Fourth nodal function: */
-	*(dl1dl7+NUMNODESMINI*0+3)=-0.5*(1.0+zeta)/2.0;
-	*(dl1dl7+NUMNODESMINI*1+3)=-SQRT3/6.0*(1.0+zeta)/2.0;
-	*(dl1dl7+NUMNODESMINI*2+3)=0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
-
-	/*Fith nodal function: */
-	*(dl1dl7+NUMNODESMINI*0+4)=0.5*(1.0+zeta)/2.0;
-	*(dl1dl7+NUMNODESMINI*1+4)=-SQRT3/6.0*(1.0+zeta)/2.0;
-	*(dl1dl7+NUMNODESMINI*2+4)=0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
-
-	/*Sixth nodal function: */
-	*(dl1dl7+NUMNODESMINI*0+5)=0;
-	*(dl1dl7+NUMNODESMINI*1+5)=SQRT3/3.0*(1.0+zeta)/2.0;
-	*(dl1dl7+NUMNODESMINI*2+5)=0.5*(SQRT3/3.0*s+ONETHIRD);
-
-	/*Seventh nodal function: */
-	*(dl1dl7+NUMNODESMINI*0+6)=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(SQRT3*s+1.0);
-	*(dl1dl7+NUMNODESMINI*1+6)=9.0/4.0*(1+zeta)*(1-zeta)*(SQRT3*pow(s,2.0)-2.0*s-SQRT3*pow(r,2.0));
-	*(dl1dl7+NUMNODESMINI*2+6)=27*gauss->coord1*gauss->coord2*gauss->coord3*(-2.0*zeta);
-
+	/*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*gauss->coord1;
+	/*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*gauss->coord2;
+	/*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*gauss->coord3;
+	/*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*gauss->coord1;
+	/*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*gauss->coord2;
+	/*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*gauss->coord3;
+	/*Nodal function 7*/
+	dbasis[NUMNODESMINI*0+6]=27.*(1.+zeta)*(1.-zeta)*(-.5*gauss->coord2*gauss->coord3 + .5*gauss->coord1*gauss->coord3);
+	dbasis[NUMNODESMINI*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[NUMNODESMINI*2+6]=27*gauss->coord1*gauss->coord2*gauss->coord3*(-2.0*zeta);
 }
 /*}}}*/
@@ -1396,48 +1384,35 @@
 /*}}}*/
 /*FUNCTION PentaRef::GetNodalFunctionsP1DerivativesReference {{{*/
-void PentaRef::GetNodalFunctionsP1DerivativesReference(IssmDouble* dl1dl6,GaussPenta* gauss){
+void PentaRef::GetNodalFunctionsP1DerivativesReference(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. Those values vary along xi,eta,z */
 
-	IssmDouble A1,A2,A3,z;
-
-	A1=gauss->coord1; _assert_(A1>=0 && A1<=1);//first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*SQRT3);
-	A2=gauss->coord2; _assert_(A2>=0 && A2<=1);//second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*SQRT3);
-	A3=gauss->coord3; _assert_(A3>=0 && A3<=1);//third area coordinate value  In term of xi and eta: A3=y/SQRT3;
-	z =gauss->coord4; _assert_(z>=-1 &&  z<=1);//fourth vertical coordinate value. Corresponding nodal function: (1-z)/2 and (1+z)/2
-
-	/*First nodal function derivatives. The corresponding nodal function is N=A1*(1-z)/2. Its derivatives follow*/
-	*(dl1dl6+NUMNODESP1*0+0)=-0.5*(1.0-z)/2.0;
-	*(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*/
-	*(dl1dl6+NUMNODESP1*0+1)=0.5*(1.0-z)/2.0;
-	*(dl1dl6+NUMNODESP1*1+1)=-0.5/SQRT3*(1.0-z)/2.0;
-	*(dl1dl6+NUMNODESP1*2+1)=-0.5*A2;
-
-	/*Third nodal function: The corresponding nodal function is N=A3*(1-z)/2. Its derivatives follow*/
-	*(dl1dl6+NUMNODESP1*0+2)=0.0;
-	*(dl1dl6+NUMNODESP1*1+2)=1.0/SQRT3*(1.0-z)/2.0;
-	*(dl1dl6+NUMNODESP1*2+2)=-0.5*A3;
-
-	/*Fourth nodal function: The corresponding nodal function is N=A1*(1+z)/2. Its derivatives follow*/
-	*(dl1dl6+NUMNODESP1*0+3)=-0.5*(1.0+z)/2.0;
-	*(dl1dl6+NUMNODESP1*1+3)=-0.5/SQRT3*(1.0+z)/2.0;
-	*(dl1dl6+NUMNODESP1*2+3)=0.5*A1;
-
-	/*Fifth nodal function: The corresponding nodal function is N=A2*(1+z)/2. Its derivatives follow*/
-	*(dl1dl6+NUMNODESP1*0+4)=0.5*(1.0+z)/2.0;
-	*(dl1dl6+NUMNODESP1*1+4)=-0.5/SQRT3*(1.0+z)/2.0;
-	*(dl1dl6+NUMNODESP1*2+4)=0.5*A2;
-
-	/*Sixth nodal function: The corresponding nodal function is N=A3*(1+z)/2. Its derivatives follow*/
-	*(dl1dl6+NUMNODESP1*0+5)=0.0;
-	*(dl1dl6+NUMNODESP1*1+5)=1.0/SQRT3*(1.0+z)/2.0;
-	*(dl1dl6+NUMNODESP1*2+5)=0.5*A3;
+	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;
 }
 /*}}}*/