Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15469)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15470)
@@ -71,13 +71,13 @@
 
 	/*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*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]; 
+	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] = .5*dbasis[1][i];
+		B[NDOF2*NUMNODESP1*2+NDOF2*i+1] = .5*dbasis[0][i];
 	}
 }
@@ -97,33 +97,33 @@
 	 */
 
-	int    i;
-	IssmDouble dh1dh7[3][NUMNODESMINI];
-	IssmDouble l1l6[NUMNODESP1];
-
-	/*Get dh1dh6 in actual coordinate system: */
-	GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
-	GetNodalFunctionsP1(l1l6, gauss);
+	int i;
+	IssmDouble dbasismini[3][NUMNODESMINI];
+	IssmDouble basis[NUMNODESP1];
+
+	/*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<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<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] = basis[i];
 	}
 }
@@ -151,18 +151,18 @@
 	/*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*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*3+NDOF2*i)=(float).5*dbasis[2][i]; 
-		*(B+NDOF2*NUMNODESP1*3+NDOF2*i+1)=0.0;
-
-		*(B+NDOF2*NUMNODESP1*4+NDOF2*i)=0.0;
-		*(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=(float).5*dbasis[2][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] = .5*dbasis[1][i];
+		B[NDOF2*NUMNODESP1*2+NDOF2*i+1] = .5*dbasis[0][i];
+
+		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] = .5*dbasis[2][i];
 	}
 
@@ -189,19 +189,19 @@
 
 	/*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];
-
-		*(B+NDOF2*NUMNODESP1*1+NDOF2*i)=dbasis[0][i];
-		*(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=2.0*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*3+NDOF2*i)=dbasis[2][i]; 
-		*(B+NDOF2*NUMNODESP1*3+NDOF2*i+1)=0.0;
-
-		*(B+NDOF2*NUMNODESP1*4+NDOF2*i)=0.0;
-		*(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=dbasis[2][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+0]=dbasis[0][i];
+		B[NDOF2*NUMNODESP1*1+NDOF2*i+1]=2.*dbasis[1][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+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]; 
 	}
 }
@@ -221,26 +221,26 @@
 
 	int    i;
-	IssmDouble dh1dh7[3][NUMNODESMINI];
-
-	/*Get dh1dh6 in actual coordinate system: */
-	GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
+	IssmDouble dbasismini[3][NUMNODESMINI];
+
+	/*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<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<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.;
 	}
 
@@ -266,48 +266,48 @@
 	int i;
 
-	IssmDouble dh1dh7[3][NUMNODESMINI];
-	IssmDouble l1l6[NUMNODESP1];
-
-	/*Get dh1dh7 in actual coordinate system: */
-	GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
-	GetNodalFunctionsP1(l1l6, gauss);
+	IssmDouble dbasismini[3][NUMNODESMINI];
+	IssmDouble basis[NUMNODESP1];
+
+	/*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];
-	}
-
-	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;
+		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] = basis[i];
+		B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3] = 0.;
 	}
 
@@ -332,51 +332,48 @@
 
 	int i;
-
-	IssmDouble dh1dh6[3][NUMNODESP1];
-	IssmDouble l1l6[NUMNODESP1];
-
-	/*Get dh1dh7 in actual coordinate system: */
-	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
-	GetNodalFunctionsP1(&l1l6[0], gauss);
+	IssmDouble dbasis[3][NUMNODESP1];
+	IssmDouble basis[NUMNODESP1];
+
+	/*Get dbasismini in actual coordinate system: */
+	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
+	GetNodalFunctionsP1(&basis[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++){ //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++){
+		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] = basis[i];
+		B[(NDOF4*NUMNODESP1)*7+NDOF4*i+3] = 0.;
 	}
 
@@ -402,48 +399,48 @@
 
 	int i;
-	IssmDouble dh1dh7[3][NUMNODESMINI];
-	IssmDouble l1l6[NUMNODESP1];
-
-	/*Get dh1dh7 in actual coordinate system: */
-	GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
-	GetNodalFunctionsP1(l1l6, gauss);
+	IssmDouble dbasismini[3][NUMNODESMINI];
+	IssmDouble basis[NUMNODESP1];
+
+	/*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<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<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] = basis[i];
 	}
 
@@ -469,50 +466,48 @@
 
 	int i;
-	IssmDouble dh1dh6[3][NUMNODESP1];
-	IssmDouble l1l6[NUMNODESP1];
-
-	/*Get dh1dh7 in actual coordinate system: */
-	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
-	GetNodalFunctionsP1(l1l6, gauss);
+	IssmDouble dbasis[3][NUMNODESP1];
+	IssmDouble basis[NUMNODESP1];
+
+	/*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++){ //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++){
+		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] = basis[i];
 	}
 
@@ -533,14 +528,14 @@
 
 	/*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];
 	}
 }
@@ -560,14 +555,14 @@
 
 	/*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];
 	}
 }
@@ -578,13 +573,11 @@
 		where hi is the interpolation function for node i.*/
 
-	int i;
-	IssmDouble dh1dh6[3][NUMNODESP1];
-
-	/*Get dh1dh6 in actual coordinate system: */
-	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
+	/*Get dbasis in actual coordinate system: */
+	IssmDouble dbasis[3][NUMNODESP1];
+	GetNodalFunctionsP1Derivatives(&dbasis[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];  
 	}
 
@@ -604,15 +597,13 @@
 	 */
 
-	/*Same thing in the actual coordinate system: */
-	IssmDouble dh1dh6[3][NUMNODESP1];
-
-	/*Get dh1dh2dh3 in actual coordinates system : */
-	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
+	/*Get nodal function derivatives in actual coordinates system : */
+	IssmDouble dbasis[3][NUMNODESP1];
+	GetNodalFunctionsP1Derivatives(&dbasis[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];
 	}
 }
@@ -620,5 +611,4 @@
 /*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);
@@ -638,14 +628,13 @@
 	 **/
 
+	/*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];
 	}
 } 
@@ -653,6 +642,5 @@
 /*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: 
@@ -673,14 +661,14 @@
 
 	/*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.;
-
-		*(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.;
+	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.;
 	}
 }
@@ -688,7 +676,5 @@
 /*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: 
@@ -728,5 +714,5 @@
 
 	IssmDouble l1l2l3[NUMNODESP1_2d];
-	IssmDouble dh1dh6[3][NUMNODESP1];
+	IssmDouble dbasis[3][NUMNODESP1];
 
 	/*Get l1l2l3 in actual coordinate system: */
@@ -735,64 +721,64 @@
 	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];
 	}
 }
@@ -815,7 +801,5 @@
 	 */
 
-	int i;
 	int num_dof=2;
-
 	IssmDouble l1l2l3[NUMNODESP1_2d];
 
@@ -826,22 +810,21 @@
 
 	/*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];
 	}
 }
@@ -849,7 +832,5 @@
 /*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: 
@@ -864,9 +845,7 @@
 	 * 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: */
@@ -875,40 +854,40 @@
 	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];
 	}
 }
@@ -916,6 +895,5 @@
 /*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: 
@@ -927,7 +905,5 @@
 	 */
 
-	int i;
 	int num_dof=4;
-
 	IssmDouble l1l2l3[NUMNODESP1_2d];
 
@@ -938,22 +914,21 @@
 
 	/*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.;
 	}
 }
@@ -961,7 +936,5 @@
 /*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: 
@@ -972,9 +945,7 @@
 	 * 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: */
@@ -982,17 +953,16 @@
 	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*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];
 	}
 }
@@ -1000,11 +970,9 @@
 /*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;
@@ -1012,45 +980,43 @@
 
 	/*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);
-
-	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);
-
-	*(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+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);
-
+	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];
+
+	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+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);
 }
 /*}}}*/
@@ -1075,18 +1041,16 @@
 	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
 
-	IssmDouble x1,x2,x3,y1,y2,y3,z1,z2,z3;
-
-	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);
+	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];
 
 	/*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!");
 }
@@ -1097,14 +1061,12 @@
 	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
 
-	IssmDouble x1,x2,y1,y2,z1,z2;
-
-	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.));
+	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];
+
+	*Jdet=.5*sqrt(pow(x2-x1,2) + pow(y2-y1,2) + pow(z2-z1,2));
 	if(*Jdet<0) _error_("negative jacobian determinant!");
 
@@ -1275,19 +1237,19 @@
 /*}}}*/
 /*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]
@@ -1297,7 +1259,7 @@
 
 	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];
 	}
 
@@ -1342,26 +1304,26 @@
 /*}}}*/
 /*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: */
@@ -1376,7 +1338,7 @@
 
 	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];
 	}
 
@@ -1468,11 +1430,11 @@
 
 	/*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];
 
 }
@@ -1480,5 +1442,7 @@
 /*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;
@@ -1487,13 +1451,13 @@
 	 *   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: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 15469)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 15470)
@@ -72,6 +72,6 @@
 	/*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];
 	}
 
@@ -102,10 +102,10 @@
 	/*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];
 	}
 
@@ -136,11 +136,11 @@
 
 	/*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];
 	}
 
@@ -221,7 +221,7 @@
 
 	/*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];
 	}
 
@@ -252,5 +252,5 @@
 
 	/*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];
@@ -288,12 +288,12 @@
 	/*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];
 	}
 
@@ -323,6 +323,6 @@
 	/*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];
 	}
 
@@ -352,8 +352,8 @@
 	/*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];
 	}
 
@@ -508,6 +508,6 @@
 	 */
 	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];
 	}
 
@@ -527,32 +527,32 @@
 		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:
@@ -590,6 +590,6 @@
 	/*Assign values*/
 	xDelete<IssmDouble>(dbasis);
-	*(p+0)=dpx;
-	*(p+1)=dpy;
+	p[0]=dpx;
+	p[1]=dpy;
 
 }