Index: ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15469)
+++ ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15470)
@@ -70,15 +70,15 @@
 	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
 
 	/*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.0;
+	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*1+NDOF2*i+0] = 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]; 
+		B[NDOF2*NUMNODESP1*2+NDOF2*i+0] = .5*dbasis[1][i];
+		B[NDOF2*NUMNODESP1*2+NDOF2*i+1] = .5*dbasis[0][i];
 	}
 }
 /*}}}*/
@@ -96,35 +96,35 @@
 	 * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
 	 */
 
-	int    i;
-	IssmDouble dh1dh7[3][NUMNODESMINI];
-	IssmDouble l1l6[NUMNODESP1];
+	int i;
+	IssmDouble dbasismini[3][NUMNODESMINI];
+	IssmDouble basis[NUMNODESP1];
 
-	/*Get dh1dh6 in actual coordinate system: */
-	GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
-	GetNodalFunctionsP1(l1l6, gauss);
+	/*Get dbasis in actual coordinate system: */
+	GetNodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
+	GetNodalFunctionsP1(basis,gauss);
 
 	/*Build B: */
-	for (i=0;i<NUMNODESMINI;i++){
-		*(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=dh1dh7[0][i]; //B[0][NDOF4*i]=dh1dh6[0][i];
-		*(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=0;
-		*(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0;
-		*(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=0;
-		*(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=dh1dh7[1][i];
-		*(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0;
-		*(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=0.5*dh1dh7[1][i];
-		*(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=0.5*dh1dh7[0][i];
-		*(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=0;
-		*(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=0;
-		*(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=0;
-		*(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0;
+	for(i=0;i<NUMNODESMINI;i++){
+		B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+0] = dbasismini[0][i];
+		B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1] = 0.;
+		B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2] = 0.;
+		B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+0] = 0.;
+		B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1] = dbasismini[1][i];
+		B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2] = 0.;
+		B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+0] = 0.5*dbasismini[1][i];
+		B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1] = 0.5*dbasismini[0][i];
+		B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2] = 0.;
+		B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+0] = 0.;
+		B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1] = 0.;
+		B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2] = 0.;
 	}
 
-	for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
-		*(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0;
-		*(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0;
-		*(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0;
-		*(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=l1l6[i];
+	for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
+		B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3] = 0;
+		B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3] = 0;
+		B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3] = 0;
+		B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3] = basis[i];
 	}
 }
 /*}}}*/
@@ -150,20 +150,20 @@
 
 	/*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.0;
+		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*1+NDOF2*i+0] = 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]; 
+		B[NDOF2*NUMNODESP1*2+NDOF2*i+0] = .5*dbasis[1][i];
+		B[NDOF2*NUMNODESP1*2+NDOF2*i+1] = .5*dbasis[0][i];
 
-		*(B+NDOF2*NUMNODESP1*3+NDOF2*i)=(float).5*dbasis[2][i]; 
-		*(B+NDOF2*NUMNODESP1*3+NDOF2*i+1)=0.0;
+		B[NDOF2*NUMNODESP1*3+NDOF2*i+0] = .5*dbasis[2][i];
+		B[NDOF2*NUMNODESP1*3+NDOF2*i+1] = 0.;
 
-		*(B+NDOF2*NUMNODESP1*4+NDOF2*i)=0.0;
-		*(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=(float).5*dbasis[2][i]; 
+		B[NDOF2*NUMNODESP1*4+NDOF2*i+0] = 0.;
+		B[NDOF2*NUMNODESP1*4+NDOF2*i+1] = .5*dbasis[2][i];
 	}
 
 }
@@ -188,21 +188,21 @@
 	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss_coord);
 
 	/*Build BPrime: */
-	for (int i=0;i<NUMNODESP1;i++){
-		*(B+NDOF2*NUMNODESP1*0+NDOF2*i)=2.0*dbasis[0][i]; 
-		*(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=dbasis[1][i];
+	for(int i=0;i<NUMNODESP1;i++){
+		B[NDOF2*NUMNODESP1*0+NDOF2*i+0]=2.*dbasis[0][i]; 
+		B[NDOF2*NUMNODESP1*0+NDOF2*i+1]=dbasis[1][i];
 
-		*(B+NDOF2*NUMNODESP1*1+NDOF2*i)=dbasis[0][i];
-		*(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=2.0*dbasis[1][i];
+		B[NDOF2*NUMNODESP1*1+NDOF2*i+0]=dbasis[0][i];
+		B[NDOF2*NUMNODESP1*1+NDOF2*i+1]=2.*dbasis[1][i];
 
-		*(B+NDOF2*NUMNODESP1*2+NDOF2*i)=dbasis[1][i]; 
-		*(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=dbasis[0][i]; 
+		B[NDOF2*NUMNODESP1*2+NDOF2*i+0]=dbasis[1][i]; 
+		B[NDOF2*NUMNODESP1*2+NDOF2*i+1]=dbasis[0][i]; 
 
-		*(B+NDOF2*NUMNODESP1*3+NDOF2*i)=dbasis[2][i]; 
-		*(B+NDOF2*NUMNODESP1*3+NDOF2*i+1)=0.0;
+		B[NDOF2*NUMNODESP1*3+NDOF2*i+0]=dbasis[2][i]; 
+		B[NDOF2*NUMNODESP1*3+NDOF2*i+1]=0.;
 
-		*(B+NDOF2*NUMNODESP1*4+NDOF2*i)=0.0;
-		*(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=dbasis[2][i]; 
+		B[NDOF2*NUMNODESP1*4+NDOF2*i+0]=0.;
+		B[NDOF2*NUMNODESP1*4+NDOF2*i+1]=dbasis[2][i]; 
 	}
 }
 /*}}}*/
@@ -220,28 +220,28 @@
 	 */
 
 	int    i;
-	IssmDouble dh1dh7[3][NUMNODESMINI];
+	IssmDouble dbasismini[3][NUMNODESMINI];
 
-	/*Get dh1dh6 in actual coordinate system: */
-	GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
+	/*Get dbasis in actual coordinate system: */
+	GetNodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
 
 	/*Build Bprime: */
-	for (i=0;i<NUMNODESMINI;i++){
-		*(Bprime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=2*dh1dh7[0][i]; //Bprime[0][NDOF4*i]=dh1dh6[0][i];
-		*(Bprime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=dh1dh7[1][i];
-		*(Bprime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0;
-		*(Bprime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=dh1dh7[0][i];
-		*(Bprime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=2*dh1dh7[1][i];
-		*(Bprime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0;
-		*(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=dh1dh7[1][i];
-		*(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=dh1dh7[0][i];
-		*(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=0;
+	for(i=0;i<NUMNODESMINI;i++){
+		Bprime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+0] = 2.*dbasismini[0][i];
+		Bprime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1] = dbasismini[1][i];
+		Bprime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2] = 0.;
+		Bprime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+0] = dbasismini[0][i];
+		Bprime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1] = 2.*dbasismini[1][i];
+		Bprime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2] = 0.;
+		Bprime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+0] = dbasismini[1][i];
+		Bprime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1] = dbasismini[0][i];
+		Bprime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2] = 0.;
 	}
 
-	for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
-		*(Bprime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0;
-		*(Bprime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0;
-		*(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0;
+	for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
+		Bprime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3] = 0.;
+		Bprime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3] = 0.;
+		Bprime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3] = 0.;
 	}
 
 }
@@ -265,50 +265,50 @@
 
 	int i;
 
-	IssmDouble dh1dh7[3][NUMNODESMINI];
-	IssmDouble l1l6[NUMNODESP1];
+	IssmDouble dbasismini[3][NUMNODESMINI];
+	IssmDouble basis[NUMNODESP1];
 
-	/*Get dh1dh7 in actual coordinate system: */
-	GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
-	GetNodalFunctionsP1(l1l6, gauss);
+	/*Get dbasismini in actual coordinate system: */
+	GetNodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
+	GetNodalFunctionsP1(basis, gauss);
 
 	/*Build B: */
-	for (i=0;i<NUMNODESMINI;i++){
-		B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+0] = dh1dh7[0][i+0]; //B[0][NDOF4*i+0] = dh1dh6[0][i+0];
+	for(i=0;i<NUMNODESMINI;i++){
+		B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+0] = dbasismini[0][i+0]; //B[0][NDOF4*i+0] = dbasis[0][i+0];
 		B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1] = 0.;
 		B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2] = 0.;
 		B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+0] = 0.;
-		B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1] = dh1dh7[1][i+0];
+		B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1] = dbasismini[1][i+0];
 		B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2] = 0.;
 		B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+0] = 0.;
 		B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1] = 0.;
-		B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2] = dh1dh7[2][i+0];
-		B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+0] = (float).5*dh1dh7[1][i+0];
-		B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1] = (float).5*dh1dh7[0][i+0];
+		B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2] = dbasismini[2][i+0];
+		B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+0] = .5*dbasismini[1][i+0];
+		B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1] = .5*dbasismini[0][i+0];
 		B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2] = 0.;
-		B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+0] = (float).5*dh1dh7[2][i+0];
+		B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+0] = .5*dbasismini[2][i+0];
 		B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1] = 0.;
-		B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2] = (float).5*dh1dh7[0][i+0];
+		B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2] = .5*dbasismini[0][i+0];
 		B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+0] = 0.;
-		B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1] = (float).5*dh1dh7[2][i+0];
-		B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2] = (float).5*dh1dh7[1][i+0];
+		B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1] = .5*dbasismini[2][i+0];
+		B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2] = .5*dbasismini[1][i+0];
 		B[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+0] = 0.;
 		B[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1] = 0.;
 		B[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2] = 0.;
-		B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+0] = dh1dh7[0][i+0];
-		B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1] = dh1dh7[1][i+0];
-		B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2] = dh1dh7[2][i+0];
+		B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+0] = dbasismini[0][i+0];
+		B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1] = dbasismini[1][i+0];
+		B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2] = dbasismini[2][i+0];
 	}
 
-	for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
-		*(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0;
-		*(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0;
-		*(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0;
-		*(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=0;
-		*(B+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3)=0;
-		*(B+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3)=0;
-		*(B+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3)=l1l6[i];
-		*(B+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3)=0;
+	for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
+		B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3] = 0.;
+		B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3] = 0.;
+		B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3] = 0.;
+		B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3] = 0.;
+		B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3] = 0.;
+		B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3] = 0.;
+		B[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3] = basis[i];
+		B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3] = 0.;
 	}
 
 }
@@ -331,53 +331,50 @@
 	 */
 
 	int i;
+	IssmDouble dbasis[3][NUMNODESP1];
+	IssmDouble basis[NUMNODESP1];
 
-	IssmDouble dh1dh6[3][NUMNODESP1];
-	IssmDouble l1l6[NUMNODESP1];
+	/*Get dbasismini in actual coordinate system: */
+	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
+	GetNodalFunctionsP1(&basis[0], gauss);
 
-	/*Get dh1dh7 in actual coordinate system: */
-	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
-	GetNodalFunctionsP1(&l1l6[0], gauss);
-
 	/*Build B: */
-	for (i=0;i<NUMNODESP1;i++){
-		*(B+(NDOF4*NUMNODESP1)*0+NDOF4*i)=dh1dh6[0][i]; //B[0][NDOF4*i]=dh1dh6[0][i];
-		*(B+(NDOF4*NUMNODESP1)*0+NDOF4*i+1)=0.;
-		*(B+(NDOF4*NUMNODESP1)*0+NDOF4*i+2)=0.;
-		*(B+(NDOF4*NUMNODESP1)*1+NDOF4*i)=0.;
-		*(B+(NDOF4*NUMNODESP1)*1+NDOF4*i+1)=dh1dh6[1][i];
-		*(B+(NDOF4*NUMNODESP1)*1+NDOF4*i+2)=0.;
-		*(B+(NDOF4*NUMNODESP1)*2+NDOF4*i)=0.;
-		*(B+(NDOF4*NUMNODESP1)*2+NDOF4*i+1)=0.;
-		*(B+(NDOF4*NUMNODESP1)*2+NDOF4*i+2)=dh1dh6[2][i];
-		*(B+(NDOF4*NUMNODESP1)*3+NDOF4*i)=.5*dh1dh6[1][i]; 
-		*(B+(NDOF4*NUMNODESP1)*3+NDOF4*i+1)=.5*dh1dh6[0][i]; 
-		*(B+(NDOF4*NUMNODESP1)*3+NDOF4*i+2)=0.;
-		*(B+(NDOF4*NUMNODESP1)*4+NDOF4*i)=.5*dh1dh6[2][i];
-		*(B+(NDOF4*NUMNODESP1)*4+NDOF4*i+1)=0.;
-		*(B+(NDOF4*NUMNODESP1)*4+NDOF4*i+2)=.5*dh1dh6[0][i];
-		*(B+(NDOF4*NUMNODESP1)*5+NDOF4*i)=0.;
-		*(B+(NDOF4*NUMNODESP1)*5+NDOF4*i+1)=.5*dh1dh6[2][i];
-		*(B+(NDOF4*NUMNODESP1)*5+NDOF4*i+2)=.5*dh1dh6[1][i];
-		*(B+(NDOF4*NUMNODESP1)*6+NDOF4*i)=0.;
-		*(B+(NDOF4*NUMNODESP1)*6+NDOF4*i+1)=0.;
-		*(B+(NDOF4*NUMNODESP1)*6+NDOF4*i+2)=0.;
-		*(B+(NDOF4*NUMNODESP1)*7+NDOF4*i)=dh1dh6[0][i];
-		*(B+(NDOF4*NUMNODESP1)*7+NDOF4*i+1)=dh1dh6[1][i];
-		*(B+(NDOF4*NUMNODESP1)*7+NDOF4*i+2)=dh1dh6[2][i];
-		_assert_((NDOF4*NUMNODESP1)*7+NDOF4*i+2<8*24);
+	for(i=0;i<NUMNODESP1;i++){
+		B[(NDOF4*NUMNODESP1)*0+NDOF4*i+0] = dbasis[0][i];
+		B[(NDOF4*NUMNODESP1)*0+NDOF4*i+1] = 0.;
+		B[(NDOF4*NUMNODESP1)*0+NDOF4*i+2] = 0.;
+		B[(NDOF4*NUMNODESP1)*1+NDOF4*i+0] = 0.;
+		B[(NDOF4*NUMNODESP1)*1+NDOF4*i+1] = dbasis[1][i];
+		B[(NDOF4*NUMNODESP1)*1+NDOF4*i+2] = 0.;
+		B[(NDOF4*NUMNODESP1)*2+NDOF4*i+0] = 0.;
+		B[(NDOF4*NUMNODESP1)*2+NDOF4*i+1] = 0.;
+		B[(NDOF4*NUMNODESP1)*2+NDOF4*i+2] = dbasis[2][i];
+		B[(NDOF4*NUMNODESP1)*3+NDOF4*i+0] = .5*dbasis[1][i];
+		B[(NDOF4*NUMNODESP1)*3+NDOF4*i+1] = .5*dbasis[0][i];
+		B[(NDOF4*NUMNODESP1)*3+NDOF4*i+2] = 0.;
+		B[(NDOF4*NUMNODESP1)*4+NDOF4*i+0] = .5*dbasis[2][i];
+		B[(NDOF4*NUMNODESP1)*4+NDOF4*i+1] = 0.;
+		B[(NDOF4*NUMNODESP1)*4+NDOF4*i+2] = .5*dbasis[0][i];
+		B[(NDOF4*NUMNODESP1)*5+NDOF4*i+0] = 0.;
+		B[(NDOF4*NUMNODESP1)*5+NDOF4*i+1] = .5*dbasis[2][i];
+		B[(NDOF4*NUMNODESP1)*5+NDOF4*i+2] = .5*dbasis[1][i];
+		B[(NDOF4*NUMNODESP1)*6+NDOF4*i+0] = 0.;
+		B[(NDOF4*NUMNODESP1)*6+NDOF4*i+1] = 0.;
+		B[(NDOF4*NUMNODESP1)*6+NDOF4*i+2] = 0.;
+		B[(NDOF4*NUMNODESP1)*7+NDOF4*i+0] = dbasis[0][i];
+		B[(NDOF4*NUMNODESP1)*7+NDOF4*i+1] = dbasis[1][i];
+		B[(NDOF4*NUMNODESP1)*7+NDOF4*i+2] = dbasis[2][i];
 	}
 
-	for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
-		B[(NDOF4*NUMNODESP1)*0+NDOF4*i+3]=0.;
-		B[(NDOF4*NUMNODESP1)*1+NDOF4*i+3]=0.;
-		B[(NDOF4*NUMNODESP1)*2+NDOF4*i+3]=0.;
-		B[(NDOF4*NUMNODESP1)*3+NDOF4*i+3]=0.;
-		B[(NDOF4*NUMNODESP1)*4+NDOF4*i+3]=0.;
-		B[(NDOF4*NUMNODESP1)*5+NDOF4*i+3]=0.;
-		B[(NDOF4*NUMNODESP1)*6+NDOF4*i+3]=l1l6[i];
-		B[(NDOF4*NUMNODESP1)*7+NDOF4*i+3]=0.;
-		_assert_(((NDOF4*NUMNODESP1)*7+NDOF4*i+3)<8*24);
+	for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
+		B[(NDOF4*NUMNODESP1)*0+NDOF4*i+3] = 0.;
+		B[(NDOF4*NUMNODESP1)*1+NDOF4*i+3] = 0.;
+		B[(NDOF4*NUMNODESP1)*2+NDOF4*i+3] = 0.;
+		B[(NDOF4*NUMNODESP1)*3+NDOF4*i+3] = 0.;
+		B[(NDOF4*NUMNODESP1)*4+NDOF4*i+3] = 0.;
+		B[(NDOF4*NUMNODESP1)*5+NDOF4*i+3] = 0.;
+		B[(NDOF4*NUMNODESP1)*6+NDOF4*i+3] = basis[i];
+		B[(NDOF4*NUMNODESP1)*7+NDOF4*i+3] = 0.;
 	}
 
 }
@@ -401,50 +398,50 @@
 	 */
 
 	int i;
-	IssmDouble dh1dh7[3][NUMNODESMINI];
-	IssmDouble l1l6[NUMNODESP1];
+	IssmDouble dbasismini[3][NUMNODESMINI];
+	IssmDouble basis[NUMNODESP1];
 
-	/*Get dh1dh7 in actual coordinate system: */
-	GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
-	GetNodalFunctionsP1(l1l6, gauss);
+	/*Get dbasismini in actual coordinate system: */
+	GetNodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
+	GetNodalFunctionsP1(basis, gauss);
 
 	/*B_primeuild B_prime: */
-	for (i=0;i<NUMNODESMINI;i++){
-		*(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=dh1dh7[0][i]; //B_prime[0][NDOF4*i]=dh1dh6[0][i];
-		*(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=0;
-		*(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0;
-		*(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=0;
-		*(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=dh1dh7[1][i];
-		*(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0;
-		*(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=0;
-		*(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=0;
-		*(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=dh1dh7[2][i];
-		*(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=dh1dh7[1][i]; 
-		*(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=dh1dh7[0][i]; 
-		*(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0;
-		*(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i)=dh1dh7[2][i];
-		*(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1)=0;
-		*(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2)=dh1dh7[0][i];
-		*(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i)=0;
-		*(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1)=dh1dh7[2][i];
-		*(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2)=dh1dh7[1][i];
-		*(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i)=dh1dh7[0][i];
-		*(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1)=dh1dh7[1][i];
-		*(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2)=dh1dh7[2][i];
-		*(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i)=0;
-		*(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1)=0;
-		*(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2)=0;
+	for(i=0;i<NUMNODESMINI;i++){
+		B_prime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+0] = dbasismini[0][i];
+		B_prime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1] = 0.;
+		B_prime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2] = 0.;
+		B_prime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+0] = 0.;
+		B_prime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1] = dbasismini[1][i];
+		B_prime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2] = 0.;
+		B_prime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+0] = 0.;
+		B_prime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1] = 0.;
+		B_prime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2] = dbasismini[2][i];
+		B_prime[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+0] = dbasismini[1][i];
+		B_prime[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1] = dbasismini[0][i];
+		B_prime[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2] = 0.;
+		B_prime[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+0] = dbasismini[2][i];
+		B_prime[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1] = 0.;
+		B_prime[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2] = dbasismini[0][i];
+		B_prime[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+0] = 0.;
+		B_prime[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1] = dbasismini[2][i];
+		B_prime[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2] = dbasismini[1][i];
+		B_prime[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+0] = dbasismini[0][i];
+		B_prime[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1] = dbasismini[1][i];
+		B_prime[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2] = dbasismini[2][i];
+		B_prime[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+0] = 0.;
+		B_prime[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1] = 0.;
+		B_prime[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2] = 0.;
 	}
 
-	for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
-		*(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0;
-		*(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0;
-		*(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0;
-		*(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=0;
-		*(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3)=0;
-		*(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3)=0;
-		*(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3)=0;
-		*(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3)=l1l6[i];
+	for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
+		B_prime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3] = 0.;
+		B_prime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3] = 0.;
+		B_prime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3] = 0.;
+		B_prime[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3] = 0.;
+		B_prime[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3] = 0.;
+		B_prime[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3] = 0.;
+		B_prime[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3] = 0.;
+		B_prime[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3] = basis[i];
 	}
 
 }
@@ -468,52 +465,50 @@
 	 */
 
 	int i;
-	IssmDouble dh1dh6[3][NUMNODESP1];
-	IssmDouble l1l6[NUMNODESP1];
+	IssmDouble dbasis[3][NUMNODESP1];
+	IssmDouble basis[NUMNODESP1];
 
-	/*Get dh1dh7 in actual coordinate system: */
-	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
-	GetNodalFunctionsP1(l1l6, gauss);
+	/*Get dbasismini in actual coordinate system: */
+	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
+	GetNodalFunctionsP1(basis, gauss);
 
 	/*B_primeuild B_prime: */
-	for (i=0;i<NUMNODESP1;i++){
-		*(B_prime+(NDOF4*NUMNODESP1)*0+NDOF4*i)=dh1dh6[0][i]; //B_prime[0][NDOF4*i]=dh1dh6[0][i];
-		*(B_prime+(NDOF4*NUMNODESP1)*0+NDOF4*i+1)=0.;
-		*(B_prime+(NDOF4*NUMNODESP1)*0+NDOF4*i+2)=0.;
-		*(B_prime+(NDOF4*NUMNODESP1)*1+NDOF4*i)=0.;
-		*(B_prime+(NDOF4*NUMNODESP1)*1+NDOF4*i+1)=dh1dh6[1][i];
-		*(B_prime+(NDOF4*NUMNODESP1)*1+NDOF4*i+2)=0.;
-		*(B_prime+(NDOF4*NUMNODESP1)*2+NDOF4*i)=0.;
-		*(B_prime+(NDOF4*NUMNODESP1)*2+NDOF4*i+1)=0.;
-		*(B_prime+(NDOF4*NUMNODESP1)*2+NDOF4*i+2)=dh1dh6[2][i];
-		*(B_prime+(NDOF4*NUMNODESP1)*3+NDOF4*i)=dh1dh6[1][i]; 
-		*(B_prime+(NDOF4*NUMNODESP1)*3+NDOF4*i+1)=dh1dh6[0][i]; 
-		*(B_prime+(NDOF4*NUMNODESP1)*3+NDOF4*i+2)=0.;
-		*(B_prime+(NDOF4*NUMNODESP1)*4+NDOF4*i)=dh1dh6[2][i];
-		*(B_prime+(NDOF4*NUMNODESP1)*4+NDOF4*i+1)=0.;
-		*(B_prime+(NDOF4*NUMNODESP1)*4+NDOF4*i+2)=dh1dh6[0][i];
-		*(B_prime+(NDOF4*NUMNODESP1)*5+NDOF4*i)=0.;
-		*(B_prime+(NDOF4*NUMNODESP1)*5+NDOF4*i+1)=dh1dh6[2][i];
-		*(B_prime+(NDOF4*NUMNODESP1)*5+NDOF4*i+2)=dh1dh6[1][i];
-		*(B_prime+(NDOF4*NUMNODESP1)*6+NDOF4*i)=dh1dh6[0][i];
-		*(B_prime+(NDOF4*NUMNODESP1)*6+NDOF4*i+1)=dh1dh6[1][i];
-		*(B_prime+(NDOF4*NUMNODESP1)*6+NDOF4*i+2)=dh1dh6[2][i];
-		*(B_prime+(NDOF4*NUMNODESP1)*7+NDOF4*i)=0.;
-		*(B_prime+(NDOF4*NUMNODESP1)*7+NDOF4*i+1)=0.;
-		*(B_prime+(NDOF4*NUMNODESP1)*7+NDOF4*i+2)=0.;
-		_assert_((NDOF4*NUMNODESP1)*7+NDOF4*i+2<8*24);
+	for(i=0;i<NUMNODESP1;i++){
+		B_prime[(NDOF4*NUMNODESP1)*0+NDOF4*i+0] = dbasis[0][i];
+		B_prime[(NDOF4*NUMNODESP1)*0+NDOF4*i+1] = 0.;
+		B_prime[(NDOF4*NUMNODESP1)*0+NDOF4*i+2] = 0.;
+		B_prime[(NDOF4*NUMNODESP1)*1+NDOF4*i+0] = 0.;
+		B_prime[(NDOF4*NUMNODESP1)*1+NDOF4*i+1] = dbasis[1][i];
+		B_prime[(NDOF4*NUMNODESP1)*1+NDOF4*i+2] = 0.;
+		B_prime[(NDOF4*NUMNODESP1)*2+NDOF4*i+0] = 0.;
+		B_prime[(NDOF4*NUMNODESP1)*2+NDOF4*i+1] = 0.;
+		B_prime[(NDOF4*NUMNODESP1)*2+NDOF4*i+2] = dbasis[2][i];
+		B_prime[(NDOF4*NUMNODESP1)*3+NDOF4*i+0] = dbasis[1][i];
+		B_prime[(NDOF4*NUMNODESP1)*3+NDOF4*i+1] = dbasis[0][i];
+		B_prime[(NDOF4*NUMNODESP1)*3+NDOF4*i+2] = 0.;
+		B_prime[(NDOF4*NUMNODESP1)*4+NDOF4*i+0] = dbasis[2][i];
+		B_prime[(NDOF4*NUMNODESP1)*4+NDOF4*i+1] = 0.;
+		B_prime[(NDOF4*NUMNODESP1)*4+NDOF4*i+2] = dbasis[0][i];
+		B_prime[(NDOF4*NUMNODESP1)*5+NDOF4*i+0] = 0.;
+		B_prime[(NDOF4*NUMNODESP1)*5+NDOF4*i+1] = dbasis[2][i];
+		B_prime[(NDOF4*NUMNODESP1)*5+NDOF4*i+2] = dbasis[1][i];
+		B_prime[(NDOF4*NUMNODESP1)*6+NDOF4*i+0] = dbasis[0][i];
+		B_prime[(NDOF4*NUMNODESP1)*6+NDOF4*i+1] = dbasis[1][i];
+		B_prime[(NDOF4*NUMNODESP1)*6+NDOF4*i+2] = dbasis[2][i];
+		B_prime[(NDOF4*NUMNODESP1)*7+NDOF4*i+0] = 0.;
+		B_prime[(NDOF4*NUMNODESP1)*7+NDOF4*i+1] = 0.;
+		B_prime[(NDOF4*NUMNODESP1)*7+NDOF4*i+2] = 0.;
 	}
 
-	for (i=0;i<NUMNODESP1;i++){ //last column 
-		*(B_prime+(NDOF4*NUMNODESP1)*0+NDOF4*i+3)=0.;
-		*(B_prime+(NDOF4*NUMNODESP1)*1+NDOF4*i+3)=0.;
-		*(B_prime+(NDOF4*NUMNODESP1)*2+NDOF4*i+3)=0.;
-		*(B_prime+(NDOF4*NUMNODESP1)*3+NDOF4*i+3)=0.;
-		*(B_prime+(NDOF4*NUMNODESP1)*4+NDOF4*i+3)=0.;
-		*(B_prime+(NDOF4*NUMNODESP1)*5+NDOF4*i+3)=0.;
-		*(B_prime+(NDOF4*NUMNODESP1)*6+NDOF4*i+3)=0.;
-		*(B_prime+(NDOF4*NUMNODESP1)*7+NDOF4*i+3)=l1l6[i];
-		_assert_((NDOF4*NUMNODESP1)*7+NDOF4*i+2<8*24);
+	for(i=0;i<NUMNODESP1;i++){ //last column 
+		B_prime[(NDOF4*NUMNODESP1)*0+NDOF4*i+3] = 0.;
+		B_prime[(NDOF4*NUMNODESP1)*1+NDOF4*i+3] = 0.;
+		B_prime[(NDOF4*NUMNODESP1)*2+NDOF4*i+3] = 0.;
+		B_prime[(NDOF4*NUMNODESP1)*3+NDOF4*i+3] = 0.;
+		B_prime[(NDOF4*NUMNODESP1)*4+NDOF4*i+3] = 0.;
+		B_prime[(NDOF4*NUMNODESP1)*5+NDOF4*i+3] = 0.;
+		B_prime[(NDOF4*NUMNODESP1)*6+NDOF4*i+3] = 0.;
+		B_prime[(NDOF4*NUMNODESP1)*7+NDOF4*i+3] = basis[i];
 	}
 
 }
@@ -532,16 +527,16 @@
 	 */
 
 	/*Same thing in the actual coordinate system: */
-	IssmDouble l1l6[6];
+	IssmDouble basis[6];
 
 	/*Get dh1dh2dh3 in actual coordinates system : */
-	GetNodalFunctionsP1(l1l6, gauss);
+	GetNodalFunctionsP1(basis, gauss);
 
 	/*Build B': */
-	for (int i=0;i<NUMNODESP1;i++){
-		*(B_advec+NDOF1*NUMNODESP1*0+NDOF1*i)=l1l6[i]; 
-		*(B_advec+NDOF1*NUMNODESP1*1+NDOF1*i)=l1l6[i]; 
-		*(B_advec+NDOF1*NUMNODESP1*2+NDOF1*i)=l1l6[i]; 
+	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];
 	}
 }
 /*}}}*/
@@ -559,16 +554,16 @@
 	 */
 
 	/*Same thing in the actual coordinate system: */
-	IssmDouble dh1dh6[3][NUMNODESP1];
+	IssmDouble dbasis[3][NUMNODESP1];
 
 	/*Get dh1dh2dh3 in actual coordinates system : */
-	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
+	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
 
 	/*Build B': */
-	for (int i=0;i<NUMNODESP1;i++){
-		*(B_conduct+NDOF1*NUMNODESP1*0+NDOF1*i)=dh1dh6[0][i]; 
-		*(B_conduct+NDOF1*NUMNODESP1*1+NDOF1*i)=dh1dh6[1][i]; 
-		*(B_conduct+NDOF1*NUMNODESP1*2+NDOF1*i)=dh1dh6[2][i]; 
+	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];
 	}
 }
 /*}}}*/
@@ -577,15 +572,13 @@
 	/*	Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
 		where hi is the interpolation function for node i.*/
 
-	int i;
-	IssmDouble dh1dh6[3][NUMNODESP1];
+	/*Get dbasis in actual coordinate system: */
+	IssmDouble dbasis[3][NUMNODESP1];
+	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
 
-	/*Get dh1dh6 in actual coordinate system: */
-	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
-
 	/*Build B: */
-	for (i=0;i<NUMNODESP1;i++){
-		B[i]=dh1dh6[2][i];  
+	for(int i=0;i<NUMNODESP1;i++){
+		B[i] = dbasis[2][i];  
 	}
 
 }
@@ -603,23 +596,20 @@
 	 * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
 	 */
 
-	/*Same thing in the actual coordinate system: */
-	IssmDouble dh1dh6[3][NUMNODESP1];
+	/*Get nodal function derivatives in actual coordinates system : */
+	IssmDouble dbasis[3][NUMNODESP1];
+	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
 
-	/*Get dh1dh2dh3 in actual coordinates system : */
-	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
-
 	/*Build B': */
-	for (int i=0;i<NUMNODESP1;i++){
-		*(Bprime_advec+NDOF1*NUMNODESP1*0+NDOF1*i)=dh1dh6[0][i]; 
-		*(Bprime_advec+NDOF1*NUMNODESP1*1+NDOF1*i)=dh1dh6[1][i]; 
-		*(Bprime_advec+NDOF1*NUMNODESP1*2+NDOF1*i)=dh1dh6[2][i]; 
+	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];
 	}
 }
 /*}}}*/
 /*FUNCTION PentaRef::GetBprimeVert{{{*/
 void PentaRef::GetBprimeVert(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
-	/* Compute Bprime  matrix. Bprime=[L1 L2 L3 L4 L5 L6] where Li is the nodal function for node i*/
 
 	GetNodalFunctionsP1(B, gauss);
 
@@ -637,23 +627,21 @@
 	 ** We assume B has been allocated already, of size: 2 (2 x numnodes)
 	 **/
 
+	/*Get basis in actual coordinate system: */
 	IssmDouble basis[6];
-
-	/*Get l1l6 in actual coordinate system: */
 	GetNodalFunctionsP1(&basis[0],gauss);
 
 	for(int i=0;i<NUMNODESP1;i++){
-		B[2*NUMNODESP1*0+2*i+0]=basis[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];
+		B[2*NUMNODESP1*0+2*i+0] = basis[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];
 	}
 } 
 /*}}}*/
 /*FUNCTION PentaRef::GetLStokes{{{*/
 void PentaRef::GetLStokes(IssmDouble* LStokes, GaussPenta* gauss){
-	/*
-	 * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
+	/* Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
 	 * For node i, Li can be expressed in the actual coordinate system
 	 * by: 
 	 *       Li=[ h 0 ]
@@ -672,24 +660,22 @@
 	l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
 
 	/*Build LStokes: */
-	for (int i=0;i<3;i++){
-		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+0)=l1l2l3[i];
-		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0.;
-		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0.;
-		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0.;
+	for(int i=0;i<3;i++){
+		LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
+		LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
+		LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
+		LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
 
-		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+0)=0.;
-		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
-		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0.;
-		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0.;
+		LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
+		LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
+		LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
+		LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
 	}
 }
 /*}}}*/
 /*FUNCTION PentaRef::GetLprimeStokes {{{*/
 void PentaRef::GetLprimeStokes(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){
-
-	/*
-	 * Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 
+	/* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 
 	 * For node i, Lpi can be expressed in the actual coordinate system
 	 * by: 
 	 *       Lpi=[ h    0    0   0]1
@@ -727,73 +713,73 @@
 	int num_dof=4;
 
 	IssmDouble l1l2l3[NUMNODESP1_2d];
-	IssmDouble dh1dh6[3][NUMNODESP1];
+	IssmDouble dbasis[3][NUMNODESP1];
 
 	/*Get l1l2l3 in actual coordinate system: */
 	l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
 	l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
 	l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
 
-	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
+	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
 
 	/*Build LprimeStokes: */
-	for (i=0;i<3;i++){
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=l1l2l3[i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=l1l2l3[i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+2)=l1l2l3[i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+2)=l1l2l3[i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+2)=dh1dh6[2][i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+2)=dh1dh6[2][i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+2)=dh1dh6[2][i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i)=dh1dh6[2][i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+2)=dh1dh6[0][i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+1)=dh1dh6[2][i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+2)=dh1dh6[1][i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+2)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+3)=l1l2l3[i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+2)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+3)=l1l2l3[i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+2)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+3)=l1l2l3[i];
+	for(int i=0;i<3;i++){
+		LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0]  = l1l2l3[i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+2]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+3]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1]  = l1l2l3[i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+2]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+3]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0]  = l1l2l3[i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+2]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+3]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1]  = l1l2l3[i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+2]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+3]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+0]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+1]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+2]  = l1l2l3[i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+3]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+0]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+1]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+2]  = l1l2l3[i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+3]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+0]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+1]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+2]  = dbasis[2][i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+3]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+0]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+1]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+2]  = dbasis[2][i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+3]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*8+num_dof*i+0]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*8+num_dof*i+1]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*8+num_dof*i+2]  = dbasis[2][i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*8+num_dof*i+3]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*9+num_dof*i+0]  = dbasis[2][i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*9+num_dof*i+1]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*9+num_dof*i+2]  = dbasis[0][i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*9+num_dof*i+3]  = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*10+num_dof*i+0] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*10+num_dof*i+1] = dbasis[2][i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*10+num_dof*i+2] = dbasis[1][i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*10+num_dof*i+3] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*11+num_dof*i+0] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*11+num_dof*i+1] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*11+num_dof*i+2] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*11+num_dof*i+3] = l1l2l3[i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*12+num_dof*i+0] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*12+num_dof*i+1] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*12+num_dof*i+2] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*12+num_dof*i+3] = l1l2l3[i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*13+num_dof*i+0] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*13+num_dof*i+1] = 0;
+		LprimeStokes[num_dof*NUMNODESP1_2d*13+num_dof*i+2] = 0;
+		LprimeStokes[num_dof*NUMNODESP1_2d*13+num_dof*i+3] = l1l2l3[i];
 	}
 }
 /*}}}*/
@@ -814,9 +800,7 @@
 	 * where h is the interpolation function for node i.
 	 */
 
-	int i;
 	int num_dof=2;
-
 	IssmDouble l1l2l3[NUMNODESP1_2d];
 
 	/*Get l1l2l3 in actual coordinate system: */
@@ -825,32 +809,29 @@
 	l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
 
 	/*Build LStokes: */
-	for (i=0;i<3;i++){
-		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
-		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
-		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
-		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
-		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=l1l2l3[i];
-		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
-		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
-		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=l1l2l3[i];
-		*(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=l1l2l3[i];
-		*(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0;
-		*(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0;
-		*(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=l1l2l3[i];
-		*(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=l1l2l3[i];
-		*(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0;
-		*(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0;
-		*(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=l1l2l3[i];
-
+	for(int i=0;i<3;i++){
+		LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
+		LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0;
+		LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0;
+		LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
+		LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = l1l2l3[i];
+		LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0;
+		LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0;
+		LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = l1l2l3[i];
+		LStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+0] = l1l2l3[i];
+		LStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+1] = 0;
+		LStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+0] = 0;
+		LStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+1] = l1l2l3[i];
+		LStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+0] = l1l2l3[i];
+		LStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+1] = 0;
+		LStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+0] = 0;
+		LStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+1] = l1l2l3[i];
 	}
 }
 /*}}}*/
 /*FUNCTION PentaRef::GetLprimeMacAyealStokes {{{*/
 void PentaRef::GetLprimeMacAyealStokes(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){
-
-	/*
-	 * Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 
+	/* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 
 	 * For node i, Lpi can be expressed in the actual coordinate system
 	 * by: 
 	 *       Lpi=[ h    0    0   0]
@@ -863,60 +844,57 @@
 	 *           [ 0    0    0   h]
 	 * where h is the interpolation function for node i.
 	 */
-	int i;
 	int num_dof=4;
-
 	IssmDouble l1l2l3[NUMNODESP1_2d];
-	IssmDouble dh1dh6[3][NUMNODESP1];
+	IssmDouble dbasis[3][NUMNODESP1];
 
 	/*Get l1l2l3 in actual coordinate system: */
 	l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
 	l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
 	l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
 
-	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
+	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
 
 	/*Build LprimeStokes: */
-	for (i=0;i<3;i++){
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=l1l2l3[i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=l1l2l3[i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+2)=dh1dh6[2][i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+2)=dh1dh6[2][i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+2)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+3)=l1l2l3[i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+2)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+3)=l1l2l3[i];
+	for(int i=0;i<3;i++){
+		LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = l1l2l3[i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+3] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = l1l2l3[i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+3] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+0] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+1] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+2] = dbasis[2][i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+3] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+0] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+1] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+2] = dbasis[2][i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+3] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+0] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+1] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+2] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+3] = l1l2l3[i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+0] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+1] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+2] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+3] = l1l2l3[i];
 	}
 }
 /*}}}*/
 /*FUNCTION PentaRef::GetLStokesMacAyeal {{{*/
 void PentaRef::GetLStokesMacAyeal(IssmDouble* LStokes, GaussPenta* gauss){
-	/*
-	 * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
+	/* Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
 	 * For node i, Li can be expressed in the actual coordinate system
 	 * by: 
 	 *       Li=[ h    0    0   0]
@@ -926,9 +904,7 @@
 	 * where h is the interpolation function for node i.
 	 */
 
-	int i;
 	int num_dof=4;
-
 	IssmDouble l1l2l3[NUMNODESP1_2d];
 
 	/*Get l1l2l3 in actual coordinate system: */
@@ -937,32 +913,29 @@
 	l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
 
 	/*Build LStokes: */
-	for (i=0;i<3;i++){
-		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
-		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
-		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0;
-		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0;
-		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
-		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
-		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0;
-		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0;
-		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=0;
-		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
-		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=l1l2l3[i];
-		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0;
-		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
-		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0;
-		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=l1l2l3[i];
-		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0;
-
+	for(int i=0;i<3;i++){
+		LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
+		LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
+		LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
+		LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
+		LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
+		LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
+		LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
+		LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
+		LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = 0.;
+		LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
+		LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = l1l2l3[i];
+		LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+3] = 0.;
+		LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
+		LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = 0.;
+		LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = l1l2l3[i];
+		LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+3] = 0.;
 	}
 }
 /*}}}*/
 /*FUNCTION PentaRef::GetLprimeStokesMacAyeal {{{*/
 void PentaRef::GetLprimeStokesMacAyeal(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){
-
-	/*
-	 * Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 
+	/* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 
 	 * For node i, Lpi can be expressed in the actual coordinate system
 	 * by: 
 	 *       Lpi=[ h    0 ]
@@ -971,87 +944,80 @@
 	 *		       [ 0    h ]
 	 * where h is the interpolation function for node i.
 	 */
-	int i;
 	int num_dof=2;
-
 	IssmDouble l1l2l3[NUMNODESP1_2d];
-	IssmDouble dh1dh6[3][NUMNODESP1];
+	IssmDouble dbasis[3][NUMNODESP1];
 
 	/*Get l1l2l3 in actual coordinate system: */
 	l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
 	l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
 	l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
+	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
 
-	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
-
 	/*Build LprimeStokes: */
-	for (i=0;i<3;i++){
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=l1l2l3[i];
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
-		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=l1l2l3[i];
+	for(int i=0;i<3;i++){
+		LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = l1l2l3[i];
+		LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
+		LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = l1l2l3[i];
 	}
 }
 /*}}}*/
 /*FUNCTION PentaRef::GetJacobian {{{*/
 void PentaRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussPenta* gauss){
-
 	/*The Jacobian is constant over the element, discard the gaussian points. 
 	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
 
-	IssmDouble A1,A2,A3; //area coordinates
-	IssmDouble xi,eta,zi; //parametric coordinates
-
+	IssmDouble A1,A2,A3;  // area coordinates
+	IssmDouble xi,eta,zi; // parametric coordinates
 	IssmDouble x1,x2,x3,x4,x5,x6;
 	IssmDouble y1,y2,y3,y4,y5,y6;
 	IssmDouble z1,z2,z3,z4,z5,z6;
 
 	/*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */
-	A1=gauss->coord1;
-	A2=gauss->coord2;
-	A3=gauss->coord3;
-
+	A1  = gauss->coord1;
+	A2  = gauss->coord2;
+	A3  = gauss->coord3;
 	xi  = A2-A1;
 	eta = SQRT3*A3;
 	zi  = gauss->coord4;
 
-	x1=*(xyz_list+3*0+0);
-	x2=*(xyz_list+3*1+0);
-	x3=*(xyz_list+3*2+0);
-	x4=*(xyz_list+3*3+0);
-	x5=*(xyz_list+3*4+0);
-	x6=*(xyz_list+3*5+0);
+	x1=xyz_list[3*0+0];
+	x2=xyz_list[3*1+0];
+	x3=xyz_list[3*2+0];
+	x4=xyz_list[3*3+0];
+	x5=xyz_list[3*4+0];
+	x6=xyz_list[3*5+0];
 
-	y1=*(xyz_list+3*0+1);
-	y2=*(xyz_list+3*1+1);
-	y3=*(xyz_list+3*2+1);
-	y4=*(xyz_list+3*3+1);
-	y5=*(xyz_list+3*4+1);
-	y6=*(xyz_list+3*5+1);
+	y1=xyz_list[3*0+1];
+	y2=xyz_list[3*1+1];
+	y3=xyz_list[3*2+1];
+	y4=xyz_list[3*3+1];
+	y5=xyz_list[3*4+1];
+	y6=xyz_list[3*5+1];
 
-	z1=*(xyz_list+3*0+2);
-	z2=*(xyz_list+3*1+2);
-	z3=*(xyz_list+3*2+2);
-	z4=*(xyz_list+3*3+2);
-	z5=*(xyz_list+3*4+2);
-	z6=*(xyz_list+3*5+2);
+	z1=xyz_list[3*0+2];
+	z2=xyz_list[3*1+2];
+	z3=xyz_list[3*2+2];
+	z4=xyz_list[3*3+2];
+	z5=xyz_list[3*4+2];
+	z6=xyz_list[3*5+2];
 
-	*(J+NDOF3*0+0)=0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5);
-	*(J+NDOF3*1+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*zi+SQRT3/12.0*(-x1-x2+2*x3-x4-x5+2*x6);
-	*(J+NDOF3*2+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*eta+1/4*(x1-x2-x4+x5)*xi +0.25*(-x1+x5-x2+x4);
+	J[NDOF3*0+0] = 0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5);
+	J[NDOF3*1+0] = SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*zi+SQRT3/12.0*(-x1-x2+2*x3-x4-x5+2*x6);
+	J[NDOF3*2+0] = SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*eta+1/4*(x1-x2-x4+x5)*xi +0.25*(-x1+x5-x2+x4);
 
-	*(J+NDOF3*0+1)=0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5);
-	*(J+NDOF3*1+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*zi+SQRT3/12.0*(-y1-y2+2*y3-y4-y5+2*y6);
-	*(J+NDOF3*2+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*eta+0.25*(y1-y2-y4+y5)*xi+0.25*(y4-y1+y5-y2);
+	J[NDOF3*0+1] = 0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5);
+	J[NDOF3*1+1] = SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*zi+SQRT3/12.0*(-y1-y2+2*y3-y4-y5+2*y6);
+	J[NDOF3*2+1] = SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*eta+0.25*(y1-y2-y4+y5)*xi+0.25*(y4-y1+y5-y2);
 
-	*(J+NDOF3*0+2)=0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5);
-	*(J+NDOF3*1+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*zi+SQRT3/12.0*(-z1-z2+2*z3-z4-z5+2*z6);
-	*(J+NDOF3*2+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4);
-
+	J[NDOF3*0+2] = 0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5);
+	J[NDOF3*1+2] = SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*zi+SQRT3/12.0*(-z1-z2+2*z3-z4-z5+2*z6);
+	J[NDOF3*2+2] = SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4);
 }
 /*}}}*/
 /*FUNCTION PentaRef::GetJacobianDeterminant {{{*/
@@ -1074,20 +1040,18 @@
 	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
 	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
 
-	IssmDouble x1,x2,x3,y1,y2,y3,z1,z2,z3;
+	IssmDouble x1=xyz_list[3*0+0];
+	IssmDouble y1=xyz_list[3*0+1];
+	IssmDouble z1=xyz_list[3*0+2];
+	IssmDouble x2=xyz_list[3*1+0];
+	IssmDouble y2=xyz_list[3*1+1];
+	IssmDouble z2=xyz_list[3*1+2];
+	IssmDouble x3=xyz_list[3*2+0];
+	IssmDouble y3=xyz_list[3*2+1];
+	IssmDouble z3=xyz_list[3*2+2];
 
-	x1=*(xyz_list+3*0+0);
-	y1=*(xyz_list+3*0+1);
-	z1=*(xyz_list+3*0+2);
-	x2=*(xyz_list+3*1+0);
-	y2=*(xyz_list+3*1+1);
-	z2=*(xyz_list+3*1+2);
-	x3=*(xyz_list+3*2+0);
-	y3=*(xyz_list+3*2+1);
-	z3=*(xyz_list+3*2+2);
-
 	/*Jdet = norm( AB ^ AC ) / (2 * area of the reference triangle), with areaRef=sqrt(3) */
-	*Jdet=SQRT3/6.0*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2.0)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2.0)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2.0),0.5);
+	*Jdet=SQRT3/6.*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2),0.5);
 	if(*Jdet<0) _error_("negative jacobian determinant!");
 }
 /*}}}*
@@ -1096,16 +1060,14 @@
 	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
 	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
 
-	IssmDouble x1,x2,y1,y2,z1,z2;
+	IssmDouble x1=xyz_list[3*0+0];
+	IssmDouble y1=xyz_list[3*0+1];
+	IssmDouble z1=xyz_list[3*0+2];
+	IssmDouble x2=xyz_list[3*1+0];
+	IssmDouble y2=xyz_list[3*1+1];
+	IssmDouble z2=xyz_list[3*1+2];
 
-	x1=*(xyz_list+3*0+0);
-	y1=*(xyz_list+3*0+1);
-	z1=*(xyz_list+3*0+2);
-	x2=*(xyz_list+3*1+0);
-	y2=*(xyz_list+3*1+1);
-	z2=*(xyz_list+3*1+2);
-
-	*Jdet=1.0/2.0*sqrt(pow(x2-x1,2.) + pow(y2-y1,2.) + pow(z2-z1,2.));
+	*Jdet=.5*sqrt(pow(x2-x1,2) + pow(y2-y1,2) + pow(z2-z1,2));
 	if(*Jdet<0) _error_("negative jacobian determinant!");
 
 }
@@ -1274,21 +1236,21 @@
 }
 /*}}}*/
 /*FUNCTION PentaRef::GetNodalFunctionsMINIDerivatives{{{*/
-void PentaRef::GetNodalFunctionsMINIDerivatives(IssmDouble* dh1dh7,IssmDouble* xyz_list, GaussPenta* gauss){
+void PentaRef::GetNodalFunctionsMINIDerivatives(IssmDouble* dbasismini,IssmDouble* xyz_list, GaussPenta* gauss){
 
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
 	 * actual coordinate system): */
 
-	IssmDouble    dh1dh7_ref[3][NUMNODESMINI];
+	IssmDouble    dbasismini_ref[3][NUMNODESMINI];
 	IssmDouble    Jinv[3][3];
 
 	/*Get derivative values with respect to parametric coordinate system: */
-	GetNodalFunctionsMINIDerivativesReference(&dh1dh7_ref[0][0], gauss); 
+	GetNodalFunctionsMINIDerivativesReference(&dbasismini_ref[0][0], gauss); 
 
 	/*Get Jacobian invert: */
 	GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
 
-	/*Build dh1dh6: 
+	/*Build dbasis: 
 	 *
 	 * [dhi/dx]= Jinv'*[dhi/dr]
 	 * [dhi/dy]        [dhi/ds]
@@ -1296,9 +1258,9 @@
 	 */
 
 	for(int i=0;i<NUMNODESMINI;i++){
-		*(dh1dh7+NUMNODESMINI*0+i)=Jinv[0][0]*dh1dh7_ref[0][i]+Jinv[0][1]*dh1dh7_ref[1][i]+Jinv[0][2]*dh1dh7_ref[2][i];
-		*(dh1dh7+NUMNODESMINI*1+i)=Jinv[1][0]*dh1dh7_ref[0][i]+Jinv[1][1]*dh1dh7_ref[1][i]+Jinv[1][2]*dh1dh7_ref[2][i];
-		*(dh1dh7+NUMNODESMINI*2+i)=Jinv[2][0]*dh1dh7_ref[0][i]+Jinv[2][1]*dh1dh7_ref[1][i]+Jinv[2][2]*dh1dh7_ref[2][i];
+		*(dbasismini+NUMNODESMINI*0+i)=Jinv[0][0]*dbasismini_ref[0][i]+Jinv[0][1]*dbasismini_ref[1][i]+Jinv[0][2]*dbasismini_ref[2][i];
+		*(dbasismini+NUMNODESMINI*1+i)=Jinv[1][0]*dbasismini_ref[0][i]+Jinv[1][1]*dbasismini_ref[1][i]+Jinv[1][2]*dbasismini_ref[2][i];
+		*(dbasismini+NUMNODESMINI*2+i)=Jinv[2][0]*dbasismini_ref[0][i]+Jinv[2][1]*dbasismini_ref[1][i]+Jinv[2][2]*dbasismini_ref[2][i];
 	}
 
 }
@@ -1341,28 +1303,28 @@
 }
 /*}}}*/
 /*FUNCTION PentaRef::GetNodalFunctionsP1 {{{*/
-void PentaRef::GetNodalFunctionsP1(IssmDouble* l1l6, GaussPenta* gauss){
+void PentaRef::GetNodalFunctionsP1(IssmDouble* basis, GaussPenta* gauss){
 	/*This routine returns the values of the nodal functions  at the gaussian point.*/
 
-	l1l6[0]=gauss->coord1*(1-gauss->coord4)/2.0;
-	l1l6[1]=gauss->coord2*(1-gauss->coord4)/2.0;
-	l1l6[2]=gauss->coord3*(1-gauss->coord4)/2.0;
-	l1l6[3]=gauss->coord1*(1+gauss->coord4)/2.0;
-	l1l6[4]=gauss->coord2*(1+gauss->coord4)/2.0;
-	l1l6[5]=gauss->coord3*(1+gauss->coord4)/2.0;
+	basis[0]=gauss->coord1*(1-gauss->coord4)/2.0;
+	basis[1]=gauss->coord2*(1-gauss->coord4)/2.0;
+	basis[2]=gauss->coord3*(1-gauss->coord4)/2.0;
+	basis[3]=gauss->coord1*(1+gauss->coord4)/2.0;
+	basis[4]=gauss->coord2*(1+gauss->coord4)/2.0;
+	basis[5]=gauss->coord3*(1+gauss->coord4)/2.0;
 
 }
 /*}}}*/
 /*FUNCTION PentaRef::GetNodalFunctionsP1Derivatives {{{*/
-void PentaRef::GetNodalFunctionsP1Derivatives(IssmDouble* dh1dh6,IssmDouble* xyz_list, GaussPenta* gauss){
+void PentaRef::GetNodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussPenta* gauss){
 
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
 	 * actual coordinate system): */
-	IssmDouble    dh1dh6_ref[NDOF3][NUMNODESP1];
+	IssmDouble    dbasis_ref[NDOF3][NUMNODESP1];
 	IssmDouble    Jinv[NDOF3][NDOF3];
 
 	/*Get derivative values with respect to parametric coordinate system: */
-	GetNodalFunctionsP1DerivativesReference(&dh1dh6_ref[0][0], gauss); 
+	GetNodalFunctionsP1DerivativesReference(&dbasis_ref[0][0], gauss); 
 
 	/*Get Jacobian invert: */
 	GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
@@ -1375,9 +1337,9 @@
 	 */
 
 	for (int i=0;i<NUMNODESP1;i++){
-		*(dh1dh6+NUMNODESP1*0+i)=Jinv[0][0]*dh1dh6_ref[0][i]+Jinv[0][1]*dh1dh6_ref[1][i]+Jinv[0][2]*dh1dh6_ref[2][i];
-		*(dh1dh6+NUMNODESP1*1+i)=Jinv[1][0]*dh1dh6_ref[0][i]+Jinv[1][1]*dh1dh6_ref[1][i]+Jinv[1][2]*dh1dh6_ref[2][i];
-		*(dh1dh6+NUMNODESP1*2+i)=Jinv[2][0]*dh1dh6_ref[0][i]+Jinv[2][1]*dh1dh6_ref[1][i]+Jinv[2][2]*dh1dh6_ref[2][i];
+		*(dbasis+NUMNODESP1*0+i)=Jinv[0][0]*dbasis_ref[0][i]+Jinv[0][1]*dbasis_ref[1][i]+Jinv[0][2]*dbasis_ref[2][i];
+		*(dbasis+NUMNODESP1*1+i)=Jinv[1][0]*dbasis_ref[0][i]+Jinv[1][1]*dbasis_ref[1][i]+Jinv[1][2]*dbasis_ref[2][i];
+		*(dbasis+NUMNODESP1*2+i)=Jinv[2][0]*dbasis_ref[0][i]+Jinv[2][1]*dbasis_ref[1][i]+Jinv[2][2]*dbasis_ref[2][i];
 	}
 
 }
@@ -1467,34 +1429,36 @@
 	/*P1 interpolation on Gauss point*/
 
 	/*intermediary*/
-	IssmDouble l1l6[6];
+	IssmDouble basis[6];
 
 	/*nodal functions: */
-	GetNodalFunctionsP1(&l1l6[0],gauss);
+	GetNodalFunctionsP1(&basis[0],gauss);
 
 	/*Assign output pointers:*/
-	*pvalue=l1l6[0]*plist[0]+l1l6[1]*plist[1]+l1l6[2]*plist[2]+l1l6[3]*plist[3]+l1l6[4]*plist[4]+l1l6[5]*plist[5];
+	*pvalue=basis[0]*plist[0]+basis[1]*plist[1]+basis[2]*plist[2]+basis[3]*plist[3]+basis[4]*plist[4]+basis[5]*plist[5];
 
 }
 /*}}}*/
 /*FUNCTION PentaRef::GetInputDerivativeValue{{{*/
 void PentaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussPenta* gauss){
-	/*From node values of parameter p (p_list[0], p_list[1], p_list[2], p_list[3], p_list[4] and p_list[4]), return parameter derivative value at gaussian point specified by gauss_coord:
+	/*From node values of parameter p (p_list[0], p_list[1], p_list[2],
+	 * p_list[3], p_list[4] and p_list[4]), return parameter derivative value at
+	 * gaussian point specified by gauss_coord:
 	 *   dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
 	 *   dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy;
 	 *   dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz;
 	 *
 	 *   p is a vector of size 3x1 already allocated.
 	 */
-	IssmDouble dh1dh6[3][NUMNODESP1];
+	IssmDouble dbasis[3][NUMNODESP1];
 
 	/*Get nodal funnctions derivatives in actual coordinate system: */
-	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
+	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
 
 	/*Assign output*/
-	p[0]=plist[0]*dh1dh6[0][0]+plist[1]*dh1dh6[0][1]+plist[2]*dh1dh6[0][2]+plist[3]*dh1dh6[0][3]+plist[4]*dh1dh6[0][4]+plist[5]*dh1dh6[0][5];
-	p[1]=plist[0]*dh1dh6[1][0]+plist[1]*dh1dh6[1][1]+plist[2]*dh1dh6[1][2]+plist[3]*dh1dh6[1][3]+plist[4]*dh1dh6[1][4]+plist[5]*dh1dh6[1][5];
-	p[2]=plist[0]*dh1dh6[2][0]+plist[1]*dh1dh6[2][1]+plist[2]*dh1dh6[2][2]+plist[3]*dh1dh6[2][3]+plist[4]*dh1dh6[2][4]+plist[5]*dh1dh6[2][5];
+	p[0]=plist[0]*dbasis[0][0]+plist[1]*dbasis[0][1]+plist[2]*dbasis[0][2]+plist[3]*dbasis[0][3]+plist[4]*dbasis[0][4]+plist[5]*dbasis[0][5];
+	p[1]=plist[0]*dbasis[1][0]+plist[1]*dbasis[1][1]+plist[2]*dbasis[1][2]+plist[3]*dbasis[1][3]+plist[4]*dbasis[1][4]+plist[5]*dbasis[1][5];
+	p[2]=plist[0]*dbasis[2][0]+plist[1]*dbasis[2][1]+plist[2]*dbasis[2][2]+plist[3]*dbasis[2][3]+plist[4]*dbasis[2][4]+plist[5]*dbasis[2][5];
 
 }
 /*}}}*/
Index: ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 15469)
+++ ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 15470)
@@ -71,8 +71,8 @@
 
 	/*Build B: */
 	for(int i=0;i<numnodes;i++){
-		B[numnodes*0+i]=dbasis[0*numnodes+i]; 
-		B[numnodes*1+i]=dbasis[1*numnodes+i]; 
+		B[numnodes*0+i] = dbasis[0*numnodes+i];
+		B[numnodes*1+i] = dbasis[1*numnodes+i];
 	}
 
 	/*Clean-up*/
@@ -101,12 +101,12 @@
 
 	/*Build B: */
 	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]; 
+		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*/
@@ -135,13 +135,13 @@
 	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
 	/*Build B': */
-	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]; 
+	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*/
@@ -220,9 +220,9 @@
 	GetNodalFunctions(basis,gauss);
 
 	/*Build B: */
-	for (int i=0;i<numnodes;i++){
-		B[numnodes*0+i]=basis[i];
-		B[numnodes*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*/
@@ -251,7 +251,7 @@
 	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
 	/*Build B': */
-	for (int i=0;i<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];
@@ -287,14 +287,14 @@
 
 	/*Build Bprime: */
 	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]; 
+		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*/
@@ -322,8 +322,8 @@
 
 	/*Build B': */
 	for(int i=0;i<numnodes;i++){
-		Bprime[numnodes*0+i]=dbasis[0*numnodes+i]; 
-		Bprime[numnodes*1+i]=dbasis[1*numnodes+i]; 
+		Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
+		Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
 	}
 
 	/*Clean-up*/
@@ -351,10 +351,10 @@
 
 	/*Build L: */
 	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];
+		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*/
@@ -507,8 +507,8 @@
 	 * [dhi/dy]       [dhi/ds]
 	 */
 	for(int i=0;i<numnodes;i++){
-		dbasis[numnodes*0+i]=Jinv[0][0]*dbasis_ref[0*numnodes+i]+Jinv[0][1]*dbasis_ref[1*numnodes+i];
-		dbasis[numnodes*1+i]=Jinv[1][0]*dbasis_ref[0*numnodes+i]+Jinv[1][1]*dbasis_ref[1*numnodes+i];
+		dbasis[numnodes*0+i] = Jinv[0][0]*dbasis_ref[0*numnodes+i]+Jinv[0][1]*dbasis_ref[1*numnodes+i];
+		dbasis[numnodes*1+i] = Jinv[1][0]*dbasis_ref[0*numnodes+i]+Jinv[1][1]*dbasis_ref[1*numnodes+i];
 	}
 
 	/*Clean up*/
@@ -526,34 +526,34 @@
 	switch(this->element_type){
 		case P1Enum: case P1DGEnum:
 			/*Nodal function 1*/
-			dbasis[NUMNODESP1*0+0]=-0.5; 
-			dbasis[NUMNODESP1*1+0]=-1.0/(2.0*SQRT3);
+			dbasis[NUMNODESP1*0+0] = -0.5;
+			dbasis[NUMNODESP1*1+0] = -1.0/(2.0*SQRT3);
 			/*Nodal function 2*/
-			dbasis[NUMNODESP1*0+1]=0.5;
-			dbasis[NUMNODESP1*1+1]=-1.0/(2.0*SQRT3);
+			dbasis[NUMNODESP1*0+1] = 0.5;
+			dbasis[NUMNODESP1*1+1] = -1.0/(2.0*SQRT3);
 			/*Nodal function 3*/
-			dbasis[NUMNODESP1*0+2]=0;
-			dbasis[NUMNODESP1*1+2]=1.0/SQRT3;
+			dbasis[NUMNODESP1*0+2] = 0;
+			dbasis[NUMNODESP1*1+2] = 1.0/SQRT3;
 			return;
 		case P2Enum:
 			/*Nodal function 1*/
-			dbasis[NUMNODESP2*0+0]=-2.*gauss->coord1 + 0.5;
-			dbasis[NUMNODESP2*1+0]=-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.;
+			dbasis[NUMNODESP2*0+0] = -2.*gauss->coord1 + 0.5;
+			dbasis[NUMNODESP2*1+0] = -2.*SQRT3/3.*gauss->coord1 + SQRT3/6.;
 			/*Nodal function 2*/
-			dbasis[NUMNODESP2*0+1]=+2.*gauss->coord2 - 0.5;
-			dbasis[NUMNODESP2*1+1]=-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.;
+			dbasis[NUMNODESP2*0+1] = +2.*gauss->coord2 - 0.5;
+			dbasis[NUMNODESP2*1+1] = -2.*SQRT3/3.*gauss->coord2 + SQRT3/6.;
 			/*Nodal function 3*/
-			dbasis[NUMNODESP2*0+2]=0.;
-			dbasis[NUMNODESP2*1+2]=+4.*SQRT3/3.*gauss->coord3 - SQRT3/3.;
+			dbasis[NUMNODESP2*0+2] = 0.;
+			dbasis[NUMNODESP2*1+2] = +4.*SQRT3/3.*gauss->coord3 - SQRT3/3.;
 			/*Nodal function 4*/
-			dbasis[NUMNODESP2*0+3]=+2.*gauss->coord3;
-			dbasis[NUMNODESP2*1+3]=+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3;
+			dbasis[NUMNODESP2*0+3] = +2.*gauss->coord3;
+			dbasis[NUMNODESP2*1+3] = +4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3;
 			/*Nodal function 5*/
-			dbasis[NUMNODESP2*0+4]=-2.*gauss->coord3;
-			dbasis[NUMNODESP2*1+4]=+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3;
+			dbasis[NUMNODESP2*0+4] = -2.*gauss->coord3;
+			dbasis[NUMNODESP2*1+4] = +4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3;
 			/*Nodal function 6*/
-			dbasis[NUMNODESP2*0+5]=2.*(gauss->coord1-gauss->coord2);
-			dbasis[NUMNODESP2*1+5]=-2.*SQRT3/3.*(gauss->coord1+gauss->coord2);
+			dbasis[NUMNODESP2*0+5] = 2.*(gauss->coord1-gauss->coord2);
+			dbasis[NUMNODESP2*1+5] = -2.*SQRT3/3.*(gauss->coord1+gauss->coord2);
 			return;
 		default:
 			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
@@ -589,8 +589,8 @@
 
 	/*Assign values*/
 	xDelete<IssmDouble>(dbasis);
-	*(p+0)=dpx;
-	*(p+1)=dpy;
+	p[0]=dpx;
+	p[1]=dpy;
 
 }
 /*}}}*/