Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15636)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15637)
@@ -1096,71 +1096,74 @@
 	_assert_(basis);
 
+	/*Get current coordinates in reference element*/
+	IssmDouble zeta=gauss->coord4;
+
 	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.;
+			basis[0]=gauss->coord1*(1.-zeta)/2.;
+			basis[1]=gauss->coord2*(1.-zeta)/2.;
+			basis[2]=gauss->coord3*(1.-zeta)/2.;
+			basis[3]=gauss->coord1*(1.+zeta)/2.;
+			basis[4]=gauss->coord2*(1.+zeta)/2.;
+			basis[5]=gauss->coord3*(1.+zeta)/2.;
 			return;
 		case P1bubbleEnum:
-			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);
+			basis[0]=gauss->coord1*(1.-zeta)/2.;
+			basis[1]=gauss->coord2*(1.-zeta)/2.;
+			basis[2]=gauss->coord3*(1.-zeta)/2.;
+			basis[3]=gauss->coord1*(1.+zeta)/2.;
+			basis[4]=gauss->coord2*(1.+zeta)/2.;
+			basis[5]=gauss->coord3*(1.+zeta)/2.;
+			basis[6]=27.*gauss->coord1*gauss->coord2*gauss->coord3*(1.+zeta)*(1.-zeta);
 			return;
 		case P2xP1Enum:
 			/*Corner nodes*/
-			basis[ 0]=gauss->coord1*(2.*gauss->coord1-1.)*(1.-gauss->coord4)/2.;
-			basis[ 1]=gauss->coord2*(2.*gauss->coord2-1.)*(1.-gauss->coord4)/2.;
-			basis[ 2]=gauss->coord3*(2.*gauss->coord3-1.)*(1.-gauss->coord4)/2.;
-			basis[ 3]=gauss->coord1*(2.*gauss->coord1-1.)*(1.+gauss->coord4)/2.;
-			basis[ 4]=gauss->coord2*(2.*gauss->coord2-1.)*(1.+gauss->coord4)/2.;
-			basis[ 5]=gauss->coord3*(2.*gauss->coord3-1.)*(1.+gauss->coord4)/2.;
+			basis[ 0]=gauss->coord1*(2.*gauss->coord1-1.)*(1.-zeta)/2.;
+			basis[ 1]=gauss->coord2*(2.*gauss->coord2-1.)*(1.-zeta)/2.;
+			basis[ 2]=gauss->coord3*(2.*gauss->coord3-1.)*(1.-zeta)/2.;
+			basis[ 3]=gauss->coord1*(2.*gauss->coord1-1.)*(1.+zeta)/2.;
+			basis[ 4]=gauss->coord2*(2.*gauss->coord2-1.)*(1.+zeta)/2.;
+			basis[ 5]=gauss->coord3*(2.*gauss->coord3-1.)*(1.+zeta)/2.;
 			/*mid-sides of triangles*/
-			basis[ 6]=4.*gauss->coord3*gauss->coord2*(1.-gauss->coord4)/2.;
-			basis[ 7]=4.*gauss->coord3*gauss->coord1*(1.-gauss->coord4)/2.;
-			basis[ 8]=4.*gauss->coord1*gauss->coord2*(1.-gauss->coord4)/2.;
-			basis[ 9]=4.*gauss->coord3*gauss->coord2*(1.+gauss->coord4)/2.;
-			basis[10]=4.*gauss->coord3*gauss->coord1*(1.+gauss->coord4)/2.;
-			basis[11]=4.*gauss->coord1*gauss->coord2*(1.+gauss->coord4)/2.;
+			basis[ 6]=4.*gauss->coord3*gauss->coord2*(1.-zeta)/2.;
+			basis[ 7]=4.*gauss->coord3*gauss->coord1*(1.-zeta)/2.;
+			basis[ 8]=4.*gauss->coord1*gauss->coord2*(1.-zeta)/2.;
+			basis[ 9]=4.*gauss->coord3*gauss->coord2*(1.+zeta)/2.;
+			basis[10]=4.*gauss->coord3*gauss->coord1*(1.+zeta)/2.;
+			basis[11]=4.*gauss->coord1*gauss->coord2*(1.+zeta)/2.;
 			return;
 		case P1xP2Enum:
 			/*Corner nodes*/
-			basis[ 0]=gauss->coord1*gauss->coord4*(gauss->coord4-1.)/2.;
-			basis[ 1]=gauss->coord2*gauss->coord4*(gauss->coord4-1.)/2.;
-			basis[ 2]=gauss->coord3*gauss->coord4*(gauss->coord4-1.)/2.;
-			basis[ 3]=gauss->coord1*gauss->coord4*(gauss->coord4+1.)/2.;
-			basis[ 4]=gauss->coord2*gauss->coord4*(gauss->coord4+1.)/2.;
-			basis[ 5]=gauss->coord3*gauss->coord4*(gauss->coord4+1.)/2.;
+			basis[ 0]=gauss->coord1*zeta*(zeta-1.)/2.;
+			basis[ 1]=gauss->coord2*zeta*(zeta-1.)/2.;
+			basis[ 2]=gauss->coord3*zeta*(zeta-1.)/2.;
+			basis[ 3]=gauss->coord1*zeta*(zeta+1.)/2.;
+			basis[ 4]=gauss->coord2*zeta*(zeta+1.)/2.;
+			basis[ 5]=gauss->coord3*zeta*(zeta+1.)/2.;
 			/*mid-sides of quads*/
-			basis[ 6]=gauss->coord1*(1.-gauss->coord4*gauss->coord4);
-			basis[ 7]=gauss->coord2*(1.-gauss->coord4*gauss->coord4);
-			basis[ 8]=gauss->coord3*(1.-gauss->coord4*gauss->coord4);
+			basis[ 6]=gauss->coord1*(1.-zeta*zeta);
+			basis[ 7]=gauss->coord2*(1.-zeta*zeta);
+			basis[ 8]=gauss->coord3*(1.-zeta*zeta);
 			return;
 		case P2Enum:
 			/*Corner nodes*/
-			basis[ 0]=gauss->coord1*(2.*gauss->coord1-1.)*gauss->coord4*(gauss->coord4-1.)/2.;
-			basis[ 1]=gauss->coord2*(2.*gauss->coord2-1.)*gauss->coord4*(gauss->coord4-1.)/2.;
-			basis[ 2]=gauss->coord3*(2.*gauss->coord3-1.)*gauss->coord4*(gauss->coord4-1.)/2.;
-			basis[ 3]=gauss->coord1*(2.*gauss->coord1-1.)*gauss->coord4*(gauss->coord4+1.)/2.;
-			basis[ 4]=gauss->coord2*(2.*gauss->coord2-1.)*gauss->coord4*(gauss->coord4+1.)/2.;
-			basis[ 5]=gauss->coord3*(2.*gauss->coord3-1.)*gauss->coord4*(gauss->coord4+1.)/2.;
+			basis[ 0]=gauss->coord1*(2.*gauss->coord1-1.)*zeta*(zeta-1.)/2.;
+			basis[ 1]=gauss->coord2*(2.*gauss->coord2-1.)*zeta*(zeta-1.)/2.;
+			basis[ 2]=gauss->coord3*(2.*gauss->coord3-1.)*zeta*(zeta-1.)/2.;
+			basis[ 3]=gauss->coord1*(2.*gauss->coord1-1.)*zeta*(zeta+1.)/2.;
+			basis[ 4]=gauss->coord2*(2.*gauss->coord2-1.)*zeta*(zeta+1.)/2.;
+			basis[ 5]=gauss->coord3*(2.*gauss->coord3-1.)*zeta*(zeta+1.)/2.;
 			/*mid-sides of quads*/
-			basis[ 6]=gauss->coord1*(1.-gauss->coord4*gauss->coord4);
-			basis[ 7]=gauss->coord2*(1.-gauss->coord4*gauss->coord4);
-			basis[ 8]=gauss->coord3*(1.-gauss->coord4*gauss->coord4);
+			basis[ 6]=gauss->coord1*(1.-zeta*zeta);
+			basis[ 7]=gauss->coord2*(1.-zeta*zeta);
+			basis[ 8]=gauss->coord3*(1.-zeta*zeta);
 			/*mid-sides of triangles*/
-			basis[ 9]=2.*gauss->coord3*gauss->coord2*gauss->coord4*(gauss->coord4-1.);
-			basis[10]=2.*gauss->coord3*gauss->coord1*gauss->coord4*(gauss->coord4-1.);
-			basis[11]=2.*gauss->coord1*gauss->coord2*gauss->coord4*(gauss->coord4-1.);
-			basis[12]=2.*gauss->coord3*gauss->coord2*gauss->coord4*(gauss->coord4+1.);
-			basis[13]=2.*gauss->coord3*gauss->coord1*gauss->coord4*(gauss->coord4+1.);
-			basis[14]=2.*gauss->coord1*gauss->coord2*gauss->coord4*(gauss->coord4+1.);
+			basis[ 9]=2.*gauss->coord3*gauss->coord2*zeta*(zeta-1.);
+			basis[10]=2.*gauss->coord3*gauss->coord1*zeta*(zeta-1.);
+			basis[11]=2.*gauss->coord1*gauss->coord2*zeta*(zeta-1.);
+			basis[12]=2.*gauss->coord3*gauss->coord2*zeta*(zeta+1.);
+			basis[13]=2.*gauss->coord3*gauss->coord1*zeta*(zeta+1.);
+			basis[14]=2.*gauss->coord1*gauss->coord2*zeta*(zeta+1.);
 			return;
 		default:
@@ -1325,27 +1328,27 @@
 		case P1xP2Enum:
 			/*Nodal function 1*/
-			dbasis[NUMNODESP1xP2*0+0]   = zeta*(zeta-1.)/4.;
-			dbasis[NUMNODESP1xP2*1+0]   = SQRT3/12.*zeta*(zeta-1.);
-			dbasis[NUMNODESP1xP2*2+0]   = .5*(1.-2.*zeta)*gauss->coord1;
+			dbasis[NUMNODESP1xP2*0+0]   = -zeta*(zeta-1.)/4.;
+			dbasis[NUMNODESP1xP2*1+0]   = -SQRT3/12.*zeta*(zeta-1.);
+			dbasis[NUMNODESP1xP2*2+0]   = .5*(2.*zeta-1.)*gauss->coord1;
 			/*Nodal function 2*/
-			dbasis[NUMNODESP1xP2*0+1]   = zeta*(1.-zeta)/4.;
-			dbasis[NUMNODESP1xP2*1+1]   = SQRT3/12.*zeta*(zeta-1);
-			dbasis[NUMNODESP1xP2*2+1]   = .5*(1.-2.*zeta)*gauss->coord2;
+			dbasis[NUMNODESP1xP2*0+1]   = zeta*(zeta-1.)/4.;
+			dbasis[NUMNODESP1xP2*1+1]   = -SQRT3/12.*zeta*(zeta-1);
+			dbasis[NUMNODESP1xP2*2+1]   = .5*(2.*zeta-1.)*gauss->coord2;
 			/*Nodal function 3*/
 			dbasis[NUMNODESP1xP2*0+2]   = 0.;
-			dbasis[NUMNODESP1xP2*1+2]   = SQRT3/6.*zeta*(1.-zeta);
-			dbasis[NUMNODESP1xP2*2+2]   = .5*(1.-2.*zeta)*gauss->coord3;
+			dbasis[NUMNODESP1xP2*1+2]   = SQRT3/6.*zeta*(zeta-1.);
+			dbasis[NUMNODESP1xP2*2+2]   = .5*(2.*zeta-1.)*gauss->coord3;
 			/*Nodal function 4*/
-			dbasis[NUMNODESP1xP2*0+3]   = -zeta*(1.+zeta)/4.;
-			dbasis[NUMNODESP1xP2*1+3]   = -SQRT3/12.*(1.+zeta);
-			dbasis[NUMNODESP1xP2*2+3]   = .5*(1.+2.*zeta)*gauss->coord1;
+			dbasis[NUMNODESP1xP2*0+3]   = -zeta*(zeta+1)/4.;
+			dbasis[NUMNODESP1xP2*1+3]   = -SQRT3/12.*zeta*(zeta+1.);
+			dbasis[NUMNODESP1xP2*2+3]   = .5*(2.*zeta+1.)*gauss->coord1;
 			/*Nodal function 5*/
-			dbasis[NUMNODESP1xP2*0+4]   = zeta*(1.+zeta)/4.;
-			dbasis[NUMNODESP1xP2*1+4]   = -SQRT3/12.*zeta*(1.+zeta);
-			dbasis[NUMNODESP1xP2*2+4]   = .5*(1.+2.*zeta)*gauss->coord2;
+			dbasis[NUMNODESP1xP2*0+4]   = zeta*(zeta+1.)/4.;
+			dbasis[NUMNODESP1xP2*1+4]   = -SQRT3/12.*zeta*(zeta+1.);
+			dbasis[NUMNODESP1xP2*2+4]   = .5*(2.*zeta+1.)*gauss->coord2;
 			/*Nodal function 6*/
 			dbasis[NUMNODESP1xP2*0+5]   = 0.;
-			dbasis[NUMNODESP1xP2*1+5]   = SQRT3/6.*zeta*(1.+zeta);
-			dbasis[NUMNODESP1xP2*2+5]   = .5*(1.+2.*zeta)*gauss->coord3;
+			dbasis[NUMNODESP1xP2*1+5]   = SQRT3/6.*zeta*(zeta+1.);
+			dbasis[NUMNODESP1xP2*2+5]   = .5*(2.*zeta+1.)*gauss->coord3;
 
 			/*Nodal function 7*/