Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15641)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15642)
@@ -144,5 +144,4 @@
 	 * We assume B has been allocated already, of size: 5x(NDOF2*numnodes)
 	 */
-
 
 	/*Fetch number of nodes for this finite element*/
@@ -531,16 +530,20 @@
 	 */
 
-	/*Same thing in the actual coordinate system: */
-	IssmDouble basis[6];
-
-	/*Get dh1dh2dh3 in actual coordinates system : */
-	GetNodalFunctionsP1(basis, gauss);
-
-	/*Build B': */
-	for(int i=0;i<NUMNODESP1;i++){
-		B_advec[NDOF1*NUMNODESP1*0+NDOF1*i] = basis[i];
-		B_advec[NDOF1*NUMNODESP1*1+NDOF1*i] = basis[i];
-		B_advec[NDOF1*NUMNODESP1*2+NDOF1*i] = basis[i];
-	}
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* basis=xNew<IssmDouble>(numnodes);
+	GetNodalFunctions(basis,gauss);
+
+	/*Build B: */
+	for(int i=0;i<numnodes;i++){
+		B_advec[numnodes*0+i] = basis[i];
+		B_advec[numnodes*1+i] = basis[i];
+		B_advec[numnodes*2+i] = basis[i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(basis);
 }
 /*}}}*/
@@ -555,19 +558,23 @@
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
-	 */
-
-	/*Same thing in the actual coordinate system: */
-	IssmDouble dbasis[3][NUMNODESP1];
-
-	/*Get dh1dh2dh3 in actual coordinates system : */
-	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
-
-	/*Build B': */
-	for(int i=0;i<NUMNODESP1;i++){
-		B_conduct[NDOF1*NUMNODESP1*0+NDOF1*i] = dbasis[0][i];
-		B_conduct[NDOF1*NUMNODESP1*1+NDOF1*i] = dbasis[1][i];
-		B_conduct[NDOF1*NUMNODESP1*2+NDOF1*i] = dbasis[2][i];
-	}
+	 * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
+	 */
+
+	/*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);
+
+	/*Build B: */
+	for(int i=0;i<numnodes;i++){
+		B_conduct[numnodes*0+i] = dbasis[0*numnodes+i];
+		B_conduct[numnodes*1+i] = dbasis[1*numnodes+i];
+		B_conduct[numnodes*2+i] = dbasis[2*numnodes+i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
 }
 /*}}}*/
@@ -577,13 +584,18 @@
 		where hi is the interpolation function for node i.*/
 
-	/*Get dbasis in actual coordinate system: */
-	IssmDouble dbasis[3][NUMNODESP1];
-	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
+	/*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);
 
 	/*Build B: */
-	for(int i=0;i<NUMNODESP1;i++){
-		B[i] = dbasis[2][i];  
-	}
-
+	for(int i=0;i<numnodes;i++){
+		B[i] = dbasis[2*numnodes+i];  
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
 }
 /*}}}*/
@@ -598,17 +610,23 @@
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
-	 */
-
-	/*Get nodal function derivatives in actual coordinates system : */
-	IssmDouble dbasis[3][NUMNODESP1];
-	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
+	 * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
+	 */
+
+	/*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);
 
 	/*Build B': */
-	for(int i=0;i<NUMNODESP1;i++){
-		Bprime_advec[NDOF1*NUMNODESP1*0+NDOF1*i] = dbasis[0][i];
-		Bprime_advec[NDOF1*NUMNODESP1*1+NDOF1*i] = dbasis[1][i];
-		Bprime_advec[NDOF1*NUMNODESP1*2+NDOF1*i] = dbasis[2][i];
-	}
+	for(int i=0;i<numnodes;i++){
+		Bprime_advec[numnodes*0+i] = dbasis[0*numnodes+i];
+		Bprime_advec[numnodes*1+i] = dbasis[1*numnodes+i];
+		Bprime_advec[numnodes*2+i] = dbasis[2*numnodes+i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
 }
 /*}}}*/
@@ -616,5 +634,5 @@
 void PentaRef::GetBprimeVert(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
 
-	GetNodalFunctionsP1(B, gauss);
+	GetNodalFunctions(B,gauss);
 
 }
@@ -631,5 +649,4 @@
 	 ** We assume B has been allocated already, of size: 2 (2 x numnodes)
 	 **/
-
 
 	/*Fetch number of nodes for this finite element*/