Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 15313)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 15314)
@@ -63,15 +63,19 @@
 	 */
 
-	int i;
-	IssmDouble dbasis[NDOF2][NUMNODESP1];
-
-	/*Get dh1dh2dh3 in actual coordinate system: */
-	GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);
+	/*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[NDOF1*NUMNODESP1*0+NDOF1*i]=dbasis[0][i]; 
-		B[NDOF1*NUMNODESP1*1+NDOF1*i]=dbasis[1][i]; 
-	}
+	for(int i=0;i<numnodes;i++){
+		B[numnodes*0+i]=dbasis[0*numnodes+i]; 
+		B[numnodes*1+i]=dbasis[1*numnodes+i]; 
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
 }
 /*}}}*/
@@ -132,10 +136,10 @@
 	/*Build B': */
 	for (int i=0;i<NUMNODESP1;i++){
-		*(B+NDOF2*NUMNODESP1*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)=0.5*dbasis[1][i]; 
-		*(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=0.5*dbasis[0][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]; 
 	}
 }
@@ -150,12 +154,11 @@
 	 */
 
-	IssmDouble l1l3[NUMNODESP1];
-
-	GetNodalFunctions(&l1l3[0],gauss);
-
-	B[0] = +l1l3[index1];
-	B[1] = +l1l3[index2];
-	B[2] = -l1l3[index1];
-	B[3] = -l1l3[index2];
+	IssmDouble basis[NUMNODESP1];
+	GetNodalFunctions(&basis[0],gauss);
+
+	B[0] = +basis[index1];
+	B[1] = +basis[index2];
+	B[2] = -basis[index1];
+	B[3] = -basis[index2];
 }
 /*}}}*/
@@ -169,16 +172,15 @@
 	 */
 
-	IssmDouble l1l3[NUMNODESP1];
-
-	GetNodalFunctions(&l1l3[0],gauss);
-
-	Bprime[0] = l1l3[index1];
-	Bprime[1] = l1l3[index2];
-	Bprime[2] = l1l3[index1];
-	Bprime[3] = l1l3[index2];
+	IssmDouble basis[NUMNODESP1];
+	GetNodalFunctions(&basis[0],gauss);
+
+	Bprime[0] = basis[index1];
+	Bprime[1] = basis[index2];
+	Bprime[2] = basis[index1];
+	Bprime[3] = basis[index2];
 }
 /*}}}*/
 /*FUNCTION TriaRef::GetBPrognostic{{{*/
-void TriaRef::GetBPrognostic(IssmDouble* B_prog, IssmDouble* xyz_list, GaussTria* gauss){
+void TriaRef::GetBPrognostic(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){
 	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
 	 * For node i, Bi can be expressed in the actual coordinate system
@@ -196,8 +198,8 @@
 	GetNodalFunctions(&basis[0],gauss);
 
-	/*Build B_prog: */
+	/*Build B: */
 	for (int i=0;i<NUMNODESP1;i++){
-		*(B_prog+NDOF1*NUMNODESP1*0+NDOF1*i)=basis[i];
-		*(B_prog+NDOF1*NUMNODESP1*1+NDOF1*i)=basis[i];
+		B[NUMNODESP1*0+i]=basis[i];
+		B[NUMNODESP1*1+i]=basis[i];
 	}
 }
@@ -220,5 +222,5 @@
 	int numnodes = this->NumberofNodes();
 
-	/*Get nodal functions*/
+	/*Get nodal functions derivatives*/
 	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
 	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
@@ -260,18 +262,18 @@
 
 	/*Build Bprime: */
-	for (int i=0;i<NUMNODESP1;i++){
-		*(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i]; 
-		*(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0; 
-		*(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i)=0; 
-		*(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i+1)=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*3+NDOF2*i)=dbasis[0][i]; 
-		*(Bprime+NDOF2*NUMNODESP1*3+NDOF2*i+1)=dbasis[1][i]; 
+	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]; 
 	}
 }
 /*}}}*/
 /*FUNCTION TriaRef::GetBprimePrognostic{{{*/
-void TriaRef::GetBprimePrognostic(IssmDouble* Bprime_prog, IssmDouble* xyz_list, GaussTria* gauss){
+void TriaRef::GetBprimePrognostic(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss){
 	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
 	 * For node i, Bi' can be expressed in the actual coordinate system
@@ -291,11 +293,11 @@
 
 	/*Build B': */
-	for (int i=0;i<NUMNODESP1;i++){
-		*(Bprime_prog+NDOF1*NUMNODESP1*0+NDOF1*i)=dbasis[0][i]; 
-		*(Bprime_prog+NDOF1*NUMNODESP1*1+NDOF1*i)=dbasis[1][i]; 
-	}
-}
-/*}}}*/
-/*FUNCTION TriaRef::GetL{{{*/
+	for(int i=0;i<NUMNODESP1;i++){
+		Bprime[NUMNODESP1*0+i]=dbasis[0][i]; 
+		Bprime[NUMNODESP1*1+i]=dbasis[1][i]; 
+	}
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetBMacAyealFriction{{{*/
 void TriaRef::GetBMacAyealFriction(IssmDouble* B, IssmDouble* xyz_list,GaussTria* gauss){
 	/*Compute B  matrix. B=[B1 B2 B3] where Bi is square and of size 2. 
@@ -317,8 +319,8 @@
 	/*Build L: */
 	for (i=0;i<NUMNODESP1;i++){
-		*(B+2*NUMNODESP1*0+2*i)=basis[i]; //L[0][NDOF2*i]=dbasis[0][i];
-		*(B+2*NUMNODESP1*0+2*i+1)=0;
-		*(B+2*NUMNODESP1*1+2*i)=0;
-		*(B+2*NUMNODESP1*1+2*i+1)=basis[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];
 	}
 }
@@ -328,17 +330,16 @@
 	/*The Jacobian is constant over the element, discard the gaussian points. 
 	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
-	IssmDouble x1,y1,x2,y2,x3,y3;
-
-	x1=*(xyz_list+NUMNODESP1*0+0);
-	y1=*(xyz_list+NUMNODESP1*0+1);
-	x2=*(xyz_list+NUMNODESP1*1+0);
-	y2=*(xyz_list+NUMNODESP1*1+1);
-	x3=*(xyz_list+NUMNODESP1*2+0);
-	y3=*(xyz_list+NUMNODESP1*2+1);
-
-	*(J+NDOF2*0+0)=0.5*(x2-x1);
-	*(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2);
-	*(J+NDOF2*0+1)=0.5*(y2-y1);
-	*(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2);
+
+	IssmDouble x1 = xyz_list[3*0+0];
+	IssmDouble y1 = xyz_list[3*0+1];
+	IssmDouble x2 = xyz_list[3*1+0];
+	IssmDouble y2 = xyz_list[3*1+1];
+	IssmDouble x3 = xyz_list[3*2+0];
+	IssmDouble y3 = xyz_list[3*2+1];
+
+	J[2*0+0] = 0.5*(x2-x1);
+	J[2*1+0] = SQRT3/6.0*(2*x3-x1-x2);
+	J[2*0+1] = 0.5*(y2-y1);
+	J[2*1+1] = SQRT3/6.0*(2*y3-y1-y2);
 }
 /*}}}*/
@@ -347,12 +348,11 @@
 	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
 	 * J is assumed to have been allocated*/
-	IssmDouble x1,y1,x2,y2;
-
-	x1=*(xyz_list+3*0+0);
-	y1=*(xyz_list+3*0+1);
-	x2=*(xyz_list+3*1+0);
-	y2=*(xyz_list+3*1+1);
-
-	*Jdet=1.0/2.0*sqrt(pow(x2-x1,2.) + pow(y2-y1,2.));
+
+	IssmDouble x1 = xyz_list[3*0+0];
+	IssmDouble y1 = xyz_list[3*0+1];
+	IssmDouble x2 = xyz_list[3*1+0];
+	IssmDouble y2 = xyz_list[3*1+1];
+
+	*Jdet = .5*sqrt(pow(x2-x1,2) + pow(y2-y1,2));
 	if(*Jdet<0) _error_("negative jacobian determinant!");
 
@@ -391,4 +391,6 @@
 void TriaRef::GetNodalFunctions(IssmDouble* basis,GaussTria* gauss){
 	/*This routine returns the values of the nodal functions  at the gaussian point.*/
+
+	_assert_(basis);
 
 	switch(this->element_type){
@@ -483,4 +485,6 @@
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
 	 * natural coordinate system) at the gaussian point. */
+
+	_assert_(dbasis && gauss);
 
 	switch(this->element_type){