Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 15311)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 15312)
@@ -81,27 +81,31 @@
 	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by: 
-	 *       Bi=[ dh/dx           0    ]
-	 *          [   0           dh/dy  ]
-	 *          [ 1/2*dh/dy  1/2*dh/dx ]
-	 * where h is the interpolation function for node i.
-	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
-	 */
-
-	int i;
-	IssmDouble dbasis[NDOF2][NUMNODESP1];
-
-	/*Get dh1dh2dh3 in actual coordinate system: */
-	GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);
+	 *       Bi=[ dN/dx           0    ]
+	 *          [   0           dN/dy  ]
+	 *          [ 1/2*dN/dy  1/2*dN/dx ]
+	 * where N is the interpolation function for node i.
+	 *
+	 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Get nodal functions*/
+	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
 	/*Build B: */
-	for (i=0;i<NUMNODESP1;i++){
-		*(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i]; //B[0][NDOF2*i]=dbasis[0][i];
-		*(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0;
-		*(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0;
-		*(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i];
-		*(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dbasis[1][i]; 
-		*(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dbasis[0][i]; 
-	}
+	for(int i=0;i<numnodes;i++){
+		B[NDOF2*numnodes*0+NDOF2*i+0]=dbasis[0*numnodes+i];
+		B[NDOF2*numnodes*0+NDOF2*i+1]=0.;
+		B[NDOF2*numnodes*1+NDOF2*i+0]=0.;
+		B[NDOF2*numnodes*1+NDOF2*i+1]=dbasis[1*numnodes+i];
+		B[NDOF2*numnodes*2+NDOF2*i+0]=.5*dbasis[1*numnodes+i]; 
+		B[NDOF2*numnodes*2+NDOF2*i+1]=.5*dbasis[0*numnodes+i]; 
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
 }
 /*}}}*/
@@ -205,27 +209,31 @@
 	 * For node i, Bi' can be expressed in the actual coordinate system
 	 * by: 
-	 *       Bi_prime=[ 2*dh/dx    dh/dy ]
-	 *                [   dh/dx  2*dh/dy ]
-	 *                [   dh/dy    dh/dx ]
-	 * where h is the interpolation function for node i.
-	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
-	 */
-
-	/*Same thing in the actual coordinate system: */
-	IssmDouble dbasis[NDOF2][NUMNODESP1];
-
-	/*Get dh1dh2dh3 in actual coordinates system : */
-	GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);
+	 *       Bi_prime=[ 2*dN/dx    dN/dy ]
+	 *                [   dN/dx  2*dN/dy ]
+	 *                [   dN/dy    dN/dx ]
+	 * where hNis the interpolation function for node i.
+	 *
+	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Get nodal functions*/
+	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
 	/*Build B': */
 	for (int i=0;i<NUMNODESP1;i++){
-		*(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i)   = 2.*dbasis[0][i];
-		*(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i+1) = dbasis[1][i];
-		*(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i)   = dbasis[0][i];
-		*(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i+1) = 2.*dbasis[1][i];
-		*(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i)   = dbasis[1][i];
-		*(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i+1) = dbasis[0][i];
-	}
+		Bprime[NDOF2*NUMNODESP1*0+NDOF2*i+0] = 2.*dbasis[0*numnodes+i];
+		Bprime[NDOF2*NUMNODESP1*0+NDOF2*i+1] =    dbasis[1*numnodes+i];
+		Bprime[NDOF2*NUMNODESP1*1+NDOF2*i+0] =    dbasis[0*numnodes+i];
+		Bprime[NDOF2*NUMNODESP1*1+NDOF2*i+1] = 2.*dbasis[1*numnodes+i];
+		Bprime[NDOF2*NUMNODESP1*2+NDOF2*i+0] =    dbasis[1*numnodes+i];
+		Bprime[NDOF2*NUMNODESP1*2+NDOF2*i+1] =    dbasis[0*numnodes+i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
 }
 /*}}}*/
@@ -425,20 +433,20 @@
 
 	/*Get nodal functions*/
-	IssmDouble* BasisFunctions=xNew<IssmDouble>(numnodes);
-	GetNodalFunctions(BasisFunctions,gauss);
+	IssmDouble* triabasis=xNew<IssmDouble>(numnodes);
+	GetNodalFunctions(triabasis,gauss);
 
 	switch(this->element_type){
 		case P1Enum:
 		case P1DGEnum:
-			basis[0]=BasisFunctions[index1];
-			basis[1]=BasisFunctions[index2];
-			xDelete<IssmDouble>(BasisFunctions);
+			basis[0]=triabasis[index1];
+			basis[1]=triabasis[index2];
+			xDelete<IssmDouble>(triabasis);
 			return;
 		case P2Enum:
 			_assert_(index2<index1);
-			basis[0]=BasisFunctions[index1];
-			basis[1]=BasisFunctions[index2];
-			basis[2]=BasisFunctions[3+index2-1];
-			xDelete<IssmDouble>(BasisFunctions);
+			basis[0]=triabasis[index1];
+			basis[1]=triabasis[index2];
+			basis[2]=triabasis[3+index2-1];
+			xDelete<IssmDouble>(triabasis);
 			return;
 		default:
@@ -446,5 +454,6 @@
 	}
 
-	xDelete<IssmDouble>(BasisFunctions);
+	/*Clean up*/
+	xDelete<IssmDouble>(triabasis);
 }
 /*}}}*/