Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 15322)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 15323)
@@ -56,7 +56,7 @@
 	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by: 
-	 *       Bi=[ dh/dx ]
-	 *          [ dh/dy ]
-	 * where h is the interpolation function for node i.
+	 *       Bi=[ dN/dx ]
+	 *          [ dN/dy ]
+	 * where N is the interpolation function for node i.
 	 *
 	 * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
@@ -66,5 +66,5 @@
 	int numnodes = this->NumberofNodes();
 
-	/*Get nodal functions*/
+	/*Get nodal functions derivatives*/
 	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
 	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
@@ -96,5 +96,5 @@
 	int numnodes = this->NumberofNodes();
 
-	/*Get nodal functions*/
+	/*Get nodal functions derivatives*/
 	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
 	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
@@ -120,27 +120,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)
-	 */
-
-	/*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=[   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 derivatives*/
+	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
 	/*Build B': */
-	for (int i=0;i<NUMNODESP1;i++){
-		B[NDOF2*NUMNODESP1*0+NDOF2*i+0]=dbasis[0][i]; 
-		B[NDOF2*NUMNODESP1*0+NDOF2*i+1]=0.; 
-		B[NDOF2*NUMNODESP1*1+NDOF2*i+0]=0.; 
-		B[NDOF2*NUMNODESP1*1+NDOF2*i+1]=dbasis[1][i]; 
-		B[NDOF2*NUMNODESP1*2+NDOF2*i+0]=0.5*dbasis[1][i]; 
-		B[NDOF2*NUMNODESP1*2+NDOF2*i+1]=0.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]=0.5*dbasis[1*numnodes+i]; 
+		B[NDOF2*numnodes*2+NDOF2*i+1]=0.5*dbasis[0*numnodes+i]; 
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
 }
 /*}}}*/
@@ -154,11 +158,19 @@
 	 */
 
-	IssmDouble basis[NUMNODESP1];
-	GetNodalFunctions(&basis[0],gauss);
-
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Get nodal functions*/
+	IssmDouble* basis=xNew<IssmDouble>(numnodes);
+	GetNodalFunctions(basis,gauss);
+
+	/*Build B for this segment*/
 	B[0] = +basis[index1];
 	B[1] = +basis[index2];
 	B[2] = -basis[index1];
 	B[3] = -basis[index2];
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(basis);
 }
 /*}}}*/
@@ -172,11 +184,19 @@
 	 */
 
-	IssmDouble basis[NUMNODESP1];
-	GetNodalFunctions(&basis[0],gauss);
-
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Get nodal functions*/
+	IssmDouble* basis=xNew<IssmDouble>(numnodes);
+	GetNodalFunctions(basis,gauss);
+
+	/*Build B'*/
 	Bprime[0] = basis[index1];
 	Bprime[1] = basis[index2];
 	Bprime[2] = basis[index1];
 	Bprime[3] = basis[index2];
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(basis);
 }
 /*}}}*/
@@ -186,21 +206,26 @@
 	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by: 
-	 *       Bi=[ h ]
-	 *          [ h ]
-	 * where h is the interpolation function for node i.
-	 *
-	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*NUMNODESP1)
-	 */
-
-	IssmDouble basis[NUMNODESP1];
-
-	/*Get dh1dh2dh3 in actual coordinate system: */
-	GetNodalFunctions(&basis[0],gauss);
+	 *       Bi=[ N ]
+	 *          [ N ]
+	 * where N is the interpolation function for node i.
+	 *
+	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Get nodal functions*/
+	IssmDouble* basis=xNew<IssmDouble>(numnodes);
+	GetNodalFunctions(basis,gauss);
 
 	/*Build B: */
-	for (int i=0;i<NUMNODESP1;i++){
-		B[NUMNODESP1*0+i]=basis[i];
-		B[NUMNODESP1*1+i]=basis[i];
-	}
+	for (int i=0;i<numnodes;i++){
+		B[numnodes*0+i]=basis[i];
+		B[numnodes*1+i]=basis[i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(basis);
 }
 /*}}}*/
@@ -227,11 +252,11 @@
 
 	/*Build B': */
-	for (int i=0;i<NUMNODESP1;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];
+	for (int i=0;i<numnodes;i++){
+		Bprime[NDOF2*numnodes*0+NDOF2*i+0] = 2.*dbasis[0*numnodes+i];
+		Bprime[NDOF2*numnodes*0+NDOF2*i+1] =    dbasis[1*numnodes+i];
+		Bprime[NDOF2*numnodes*1+NDOF2*i+0] =    dbasis[0*numnodes+i];
+		Bprime[NDOF2*numnodes*1+NDOF2*i+1] = 2.*dbasis[1*numnodes+i];
+		Bprime[NDOF2*numnodes*2+NDOF2*i+0] =    dbasis[1*numnodes+i];
+		Bprime[NDOF2*numnodes*2+NDOF2*i+1] =    dbasis[0*numnodes+i];
 	}
 
@@ -242,34 +267,37 @@
 /*FUNCTION TriaRef::GetBprimeMacAyealStokes {{{*/
 void TriaRef::GetBprimeMacAyealStokes(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss){
-
 	/*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2. 
 	 * For node i, Bprimei can be expressed in the actual coordinate system
 	 * by: 
-	 *       Bprimei=[  dh/dx    0   ]
-	 *               [    0    dh/dy ]
-	 *               [  dh/dy  dh/dx ]
-	 *               [  dh/dx  dh/dy ]
-	 * where h is the interpolation function for node i.
-	 *
-	 * We assume Bprime 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);
+	 *       Bprimei=[  dN/dx    0   ]
+	 *               [    0    dN/dy ]
+	 *               [  dN/dy  dN/dx ]
+	 N               [  dN/dx  dN/dy ]
+	 * where N is the interpolation function for node i.
+	 *
+	 * We assume Bprime 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 Bprime: */
-	for(int i=0;i<NUMNODESP1;i++){
-		Bprime[NDOF2*NUMNODESP1*0+NDOF2*i+0]=dbasis[0][i]; 
-		Bprime[NDOF2*NUMNODESP1*0+NDOF2*i+1]=0.; 
-		Bprime[NDOF2*NUMNODESP1*1+NDOF2*i+0]=0.; 
-		Bprime[NDOF2*NUMNODESP1*1+NDOF2*i+1]=dbasis[1][i]; 
-		Bprime[NDOF2*NUMNODESP1*2+NDOF2*i+0]=dbasis[1][i]; 
-		Bprime[NDOF2*NUMNODESP1*2+NDOF2*i+1]=dbasis[0][i]; 
-		Bprime[NDOF2*NUMNODESP1*3+NDOF2*i+0]=dbasis[0][i]; 
-		Bprime[NDOF2*NUMNODESP1*3+NDOF2*i+1]=dbasis[1][i]; 
-	}
+	for(int i=0;i<numnodes;i++){
+		Bprime[NDOF2*numnodes*0+NDOF2*i+0]=dbasis[0*numnodes+i]; 
+		Bprime[NDOF2*numnodes*0+NDOF2*i+1]=0.; 
+		Bprime[NDOF2*numnodes*1+NDOF2*i+0]=0.; 
+		Bprime[NDOF2*numnodes*1+NDOF2*i+1]=dbasis[1*numnodes+i]; 
+		Bprime[NDOF2*numnodes*2+NDOF2*i+0]=dbasis[1*numnodes+i]; 
+		Bprime[NDOF2*numnodes*2+NDOF2*i+1]=dbasis[0*numnodes+i]; 
+		Bprime[NDOF2*numnodes*3+NDOF2*i+0]=dbasis[0*numnodes+i]; 
+		Bprime[NDOF2*numnodes*3+NDOF2*i+1]=dbasis[1*numnodes+i]; 
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
 }
 /*}}}*/
@@ -279,22 +307,26 @@
 	 * For node i, Bi' can be expressed in the actual coordinate system
 	 * by: 
-	 *       Bi_prime=[ dh/dx ]
-	 *                [ dh/dy ]
-	 * 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=[ dN/dx ]
+	 *                [ dN/dy ]
+	 * 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 derivatives*/
+	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
 	/*Build B': */
-	for(int i=0;i<NUMNODESP1;i++){
-		Bprime[NUMNODESP1*0+i]=dbasis[0][i]; 
-		Bprime[NUMNODESP1*1+i]=dbasis[1][i]; 
-	}
+	for(int i=0;i<numnodes;i++){
+		Bprime[numnodes*0+i]=dbasis[0*numnodes+i]; 
+		Bprime[numnodes*1+i]=dbasis[1*numnodes+i]; 
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
 }
 /*}}}*/
@@ -308,20 +340,24 @@
 	 * where N is the interpolation function for node i.
 	 *
-	 * We assume B has been allocated already, of size: 2 x (numdof*NUMNODESP1)
-	 */
-
-	int i;
-	IssmDouble basis[3];
-
-	/*Get basis in actual coordinate system: */
+	 * We assume B has been allocated already, of size: 2 x (numdof*numnodes)
+	 */
+
+	/*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 L: */
-	for (i=0;i<NUMNODESP1;i++){
-		B[2*NUMNODESP1*0+2*i+0]=basis[i]; //L[0][NDOF2*i]=dbasis[0][i];
-		B[2*NUMNODESP1*0+2*i+1]=0.;
-		B[2*NUMNODESP1*1+2*i+0]=0.;
-		B[2*NUMNODESP1*1+2*i+1]=basis[i];
-	}
+	for(int i=0;i<numnodes;i++){
+		B[2*numnodes*0+2*i+0]=basis[i];
+		B[2*numnodes*0+2*i+1]=0.;
+		B[2*numnodes*1+2*i+0]=0.;
+		B[2*numnodes*1+2*i+1]=basis[i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(basis);
 }
 /*}}}*/