Index: /issm/trunk/src/c/objects/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Penta.cpp	(revision 2955)
+++ /issm/trunk/src/c/objects/Penta.cpp	(revision 2956)
@@ -416,5 +416,5 @@
 	int     ig1,ig2;
 	double  gauss_weight1,gauss_weight2;
-	double  gauss_l1l2l3l4[4];
+	double  gauss_coord[4];
 	int     order_area_gauss;
 	int     num_vert_gauss;
@@ -571,13 +571,13 @@
 
 
-				gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1); 
-				gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
-				gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
-				gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
+				gauss_coord[0]=*(first_gauss_area_coord+ig1); 
+				gauss_coord[1]=*(second_gauss_area_coord+ig1);
+				gauss_coord[2]=*(third_gauss_area_coord+ig1);
+				gauss_coord[3]=*(fourth_gauss_vert_coord+ig2);
 
 
 				/*Get strain rate from velocity: */
-				GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3l4);
-				GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3l4);
+				GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_coord);
+				GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_coord);
 
 				/*Get viscosity: */
@@ -586,9 +586,9 @@
 
 				/*Get B and Bprime matrices: */
-				GetB(&B[0][0], &xyz_list[0][0], gauss_l1l2l3l4);
-				GetBPrime(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3l4);
+				GetB(&B[0][0], &xyz_list[0][0], gauss_coord);
+				GetBPrime(&Bprime[0][0], &xyz_list[0][0], gauss_coord);
 
 				/* Get Jacobian determinant: */
-				GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
+				GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
 
 				/*Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant 
@@ -960,5 +960,5 @@
 	int     ig1,ig2;
 	double  gauss_weight1,gauss_weight2;
-	double  gauss_l1l2l3l4[4];
+	double  gauss_coord[4];
 	int     order_area_gauss;
 	int     num_vert_gauss;
@@ -1035,15 +1035,15 @@
 			gauss_weight=gauss_weight1*gauss_weight2;
 
-			gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1); 
-			gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
-			gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
-			gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
+			gauss_coord[0]=*(first_gauss_area_coord+ig1); 
+			gauss_coord[1]=*(second_gauss_area_coord+ig1);
+			gauss_coord[2]=*(third_gauss_area_coord+ig1);
+			gauss_coord[3]=*(fourth_gauss_vert_coord+ig2);
 
 			/*Get B and Bprime matrices: */
-			GetB_vert(&B[0][0], &xyz_list[0][0], gauss_l1l2l3l4);
-			GetBPrime_vert(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3l4);
+			GetB_vert(&B[0][0], &xyz_list[0][0], gauss_coord);
+			GetBPrime_vert(&Bprime[0][0], &xyz_list[0][0], gauss_coord);
 
 			/* Get Jacobian determinant: */
-			GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
+			GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
 			DL_scalar=gauss_weight*Jdet;
 
@@ -1527,5 +1527,5 @@
 	double* area_gauss_weights      =  NULL;
 	double* vert_gauss_weights      =  NULL;
-	double  gauss_l1l2l3l4[4];
+	double  gauss_coord[4];
 	int     order_area_gauss;
 	int     num_vert_gauss;
@@ -1539,5 +1539,5 @@
 
 	/*nodal functions: */
-	double l1l2l3l4l5l6[6];
+	double l1l6[6];
 
 	/*element vector at the gaussian points: */
@@ -1609,20 +1609,20 @@
 				gauss_weight=gauss_weight1*gauss_weight2;
 
-				gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1); 
-				gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
-				gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
-				gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
+				gauss_coord[0]=*(first_gauss_area_coord+ig1); 
+				gauss_coord[1]=*(second_gauss_area_coord+ig1);
+				gauss_coord[2]=*(third_gauss_area_coord+ig1);
+				gauss_coord[3]=*(fourth_gauss_vert_coord+ig2);
 
 				/*Compute thickness at gaussian point: */
-				GetParameterValue(&thickness, &h[0],gauss_l1l2l3l4);
+				GetParameterValue(&thickness, &h[0],gauss_coord);
 
 				/*Compute slope at gaussian point: */
-				GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3l4);
+				GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_coord);
 
 				/* Get Jacobian determinant: */
-				GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
+				GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
 
 				/*Get nodal functions: */
-				GetNodalFunctions(l1l2l3l4l5l6, gauss_l1l2l3l4);
+				GetNodalFunctions(l1l6, gauss_coord);
 
 				/*Compute driving stress: */
@@ -1632,5 +1632,5 @@
 				for (i=0;i<numgrids;i++){
 					for (j=0;j<NDOF2;j++){
-						pe_g_gaussian[i*NDOF2+j]=-driving_stress_baseline*slope[j]*Jdet*gauss_weight*l1l2l3l4l5l6[i];
+						pe_g_gaussian[i*NDOF2+j]=-driving_stress_baseline*slope[j]*Jdet*gauss_weight*l1l6[i];
 					}
 				}
@@ -1921,5 +1921,5 @@
 	double* area_gauss_weights      =  NULL;
 	double* vert_gauss_weights      =  NULL;
-	double  gauss_l1l2l3l4[4];
+	double  gauss_coord[4];
 	int     order_area_gauss;
 	int     num_vert_gauss;
@@ -1935,5 +1935,5 @@
 	double  pe_g[numdof];
 	double  pe_g_gaussian[numdof];
-	double l1l2l3l4l5l6[6];
+	double l1l6[6];
 
 	/*Spawning: */
@@ -2000,12 +2000,12 @@
 			gauss_weight=gauss_weight1*gauss_weight2;
 
-			gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1); 
-			gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
-			gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
-			gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
+			gauss_coord[0]=*(first_gauss_area_coord+ig1); 
+			gauss_coord[1]=*(second_gauss_area_coord+ig1);
+			gauss_coord[2]=*(third_gauss_area_coord+ig1);
+			gauss_coord[3]=*(fourth_gauss_vert_coord+ig2);
 
 			/*Get velocity derivative, with respect to x and y: */
-			GetParameterDerivativeValue(&du[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3l4);
-			GetParameterDerivativeValue(&dv[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3l4);
+			GetParameterDerivativeValue(&du[0], &vx_list[0],&xyz_list[0][0], gauss_coord);
+			GetParameterDerivativeValue(&dv[0], &vy_list[0],&xyz_list[0][0], gauss_coord);
 			dudx=du[0];
 			dvdy=dv[1];
@@ -2013,5 +2013,5 @@
 
 			/* Get Jacobian determinant: */
-			GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
+			GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
 #ifdef _ISSM_DEBUG_ 
 			printf("Element id %i Jacobian determinant: %lf\n",GetId(),Jdet);
@@ -2019,9 +2019,9 @@
 
 			/*Get nodal functions: */
-			GetNodalFunctions(l1l2l3l4l5l6, gauss_l1l2l3l4);
+			GetNodalFunctions(l1l6, gauss_coord);
 
 			/*Build pe_g_gaussian vector: */
 			for (i=0;i<numgrids;i++){
-				pe_g_gaussian[i]=(dudx+dvdy)*Jdet*gauss_weight*l1l2l3l4l5l6[i];
+				pe_g_gaussian[i]=(dudx+dvdy)*Jdet*gauss_weight*l1l6[i];
 			}
 
@@ -2547,10 +2547,10 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::GetB"
-void Penta::GetB(double* B, double* xyz_list, double* gauss_l1l2l3l4){
+void Penta::GetB(double* B, double* xyz_list, double* gauss_coord){
 
 	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. 
-	 * For grid i, Bi can be expressed in the basic coordinate system
+	 * For grid i, Bi can be expressed in the actual coordinate system
 	 * by: 
-	 *       Bi_basic=[ dh/dx          0      ]
+	 *       Bi=[ dh/dx          0      ]
 	 *                [   0           dh/dy   ]
 	 *                [ 1/2*dh/dy  1/2*dh/dx  ]
@@ -2567,12 +2567,12 @@
 	const int NDOF2=2;
 
-	double dh1dh2dh3dh4dh5dh6_basic[NDOF3][numgrids];
-
-	/*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */
-	GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
+	double dh1dh6[NDOF3][numgrids];
+
+	/*Get dh1dh6 in actual coordinate system: */
+	GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list, gauss_coord);
 
 #ifdef _ISSM_DEBUG_ 
 	for (i=0;i<numgrids;i++){
-		printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh2dh3dh4dh5dh6_basic[0][i],dh1dh2dh3dh4dh5dh6_basic[1][i],dh1dh2dh3dh4dh5dh6_basic[2][i]);
+		printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh6[0][i],dh1dh6[1][i],dh1dh6[2][i]);
 	}
 #endif
@@ -2580,18 +2580,18 @@
 	/*Build B: */
 	for (i=0;i<numgrids;i++){
-		*(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh2dh3dh4dh5dh6_basic[0][i]; 
+		*(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh6[0][i]; 
 		*(B+NDOF2*numgrids*0+NDOF2*i+1)=0.0;
 
 		*(B+NDOF2*numgrids*1+NDOF2*i)=0.0;
-		*(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh2dh3dh4dh5dh6_basic[1][i];
-
-		*(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh2dh3dh4dh5dh6_basic[1][i]; 
-		*(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh2dh3dh4dh5dh6_basic[0][i]; 
-
-		*(B+NDOF2*numgrids*3+NDOF2*i)=(float).5*dh1dh2dh3dh4dh5dh6_basic[2][i]; 
+		*(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh6[1][i];
+
+		*(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh6[1][i]; 
+		*(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh6[0][i]; 
+
+		*(B+NDOF2*numgrids*3+NDOF2*i)=(float).5*dh1dh6[2][i]; 
 		*(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
 
 		*(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
-		*(B+NDOF2*numgrids*4+NDOF2*i+1)=(float).5*dh1dh2dh3dh4dh5dh6_basic[2][i]; 
+		*(B+NDOF2*numgrids*4+NDOF2*i+1)=(float).5*dh1dh6[2][i]; 
 	}
 
@@ -2604,7 +2604,7 @@
 
 	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 
-	 * For grid i, Bi' can be expressed in the basic coordinate system
+	 * For grid i, Bi' can be expressed in the actual coordinate system
 	 * by: 
-	 *       Bi_artdiff_basic=[ dh/dx ]
+	 *       Bi_artdiff=[ dh/dx ]
 	 *                       [ dh/dy ]
 	 * where h is the interpolation function for grid i.
@@ -2618,14 +2618,14 @@
 	int DOFPERGRID=1;
 
-	/*Same thing in the basic coordinate system: */
-	double dh1dh6_basic[calculationdof][numgrids];
-
-	/*Get dh1dh2dh3 in basic coordinates system : */
-	GetNodalFunctionsDerivativesBasic(&dh1dh6_basic[0][0],xyz_list,gauss_coord);
+	/*Same thing in the actual coordinate system: */
+	double dh1dh6[calculationdof][numgrids];
+
+	/*Get dh1dh2dh3 in actual coordinates system : */
+	GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list,gauss_coord);
 
 	/*Build B': */
 	for (i=0;i<numgrids;i++){
-		*(B_artdiff+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6_basic[0][i]; 
-		*(B_artdiff+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6_basic[1][i]; 
+		*(B_artdiff+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i]; 
+		*(B_artdiff+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i]; 
 	}
 }
@@ -2637,7 +2637,7 @@
 
 	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 
-	 * For grid i, Bi' can be expressed in the basic coordinate system
+	 * For grid i, Bi' can be expressed in the actual coordinate system
 	 * by: 
-	 *       Bi_advec_basic =[ h ]
+	 *       Bi_advec =[ h ]
 	 *                       [ h ]
 	 *                       [ h ]
@@ -2652,8 +2652,8 @@
 	int DOFPERGRID=1;
 
-	/*Same thing in the basic coordinate system: */
+	/*Same thing in the actual coordinate system: */
 	double l1l6[6];
 
-	/*Get dh1dh2dh3 in basic coordinates system : */
+	/*Get dh1dh2dh3 in actual coordinates system : */
 	GetNodalFunctions(l1l6, gauss_coord);
 
@@ -2672,7 +2672,7 @@
 
 	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 
-	 * For grid i, Bi' can be expressed in the basic coordinate system
+	 * For grid i, Bi' can be expressed in the actual coordinate system
 	 * by: 
-	 *       Bi_conduct_basic=[ dh/dx ]
+	 *       Bi_conduct=[ dh/dx ]
 	 *                       [ dh/dy ]
 	 *                       [ dh/dz ]
@@ -2687,15 +2687,15 @@
 	int DOFPERGRID=1;
 
-	/*Same thing in the basic coordinate system: */
-	double dh1dh6_basic[calculationdof][numgrids];
-
-	/*Get dh1dh2dh3 in basic coordinates system : */
-	GetNodalFunctionsDerivativesBasic(&dh1dh6_basic[0][0],xyz_list,gauss_coord);
+	/*Same thing in the actual coordinate system: */
+	double dh1dh6[calculationdof][numgrids];
+
+	/*Get dh1dh2dh3 in actual coordinates system : */
+	GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list,gauss_coord);
 
 	/*Build B': */
 	for (i=0;i<numgrids;i++){
-		*(B_conduct+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6_basic[0][i]; 
-		*(B_conduct+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6_basic[1][i]; 
-		*(B_conduct+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6_basic[2][i]; 
+		*(B_conduct+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i]; 
+		*(B_conduct+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i]; 
+		*(B_conduct+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6[2][i]; 
 	}
 }
@@ -2704,5 +2704,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta:GetB_vert"
-void Penta::GetB_vert(double* B, double* xyz_list, double* gauss_l1l2l3l4){
+void Penta::GetB_vert(double* B, double* xyz_list, double* gauss_coord){
 
 
@@ -2713,12 +2713,12 @@
 	const int NDOF3=3;
 	const int numgrids=6;
-	double dh1dh2dh3dh4dh5dh6_basic[NDOF3][numgrids];
-
-	/*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */
-	GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
+	double dh1dh6[NDOF3][numgrids];
+
+	/*Get dh1dh6 in actual coordinate system: */
+	GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list, gauss_coord);
 
 #ifdef _ISSM_DEBUG_ 
 	for (i=0;i<numgrids;i++){
-		printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh2dh3dh4dh5dh6_basic[0][i],dh1dh2dh3dh4dh5dh6_basic[1][i],dh1dh2dh3dh4dh5dh6_basic[2][i]);
+		printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh6[0][i],dh1dh6[1][i],dh1dh6[2][i]);
 	}
 #endif
@@ -2726,5 +2726,5 @@
 	/*Build B: */
 	for (i=0;i<numgrids;i++){
-		B[i]=dh1dh2dh3dh4dh5dh6_basic[2][i];  
+		B[i]=dh1dh6[2][i];  
 	}
 
@@ -2742,10 +2742,10 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::GetBPrime"
-void Penta::GetBPrime(double* B, double* xyz_list, double* gauss_l1l2l3l4){
+void Penta::GetBPrime(double* B, double* xyz_list, double* gauss_coord){
 
 	/*Compute B  prime matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. 
-	 * For grid i, Bi can be expressed in the basic coordinate system
+	 * For grid i, Bi can be expressed in the actual coordinate system
 	 * by: 
-	 *       Bi_basic=[ 2*dh/dx     dh/dy   ]
+	 *       Bi=[ 2*dh/dx     dh/dy   ]
 	 *                [   dh/dx    2*dh/dy  ]
 	 *                [ dh/dy      dh/dx    ]
@@ -2762,12 +2762,12 @@
 	const int numgrids=6;
 
-	double dh1dh2dh3dh4dh5dh6_basic[NDOF3][numgrids];
-
-	/*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */
-	GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
+	double dh1dh6[NDOF3][numgrids];
+
+	/*Get dh1dh6 in actual coordinate system: */
+	GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list, gauss_coord);
 
 #ifdef _ISSM_DEBUG_ 
 	for (i=0;i<numgrids;i++){
-		printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh2dh3dh4dh5dh6_basic[0][i],dh1dh2dh3dh4dh5dh6_basic[1][i],dh1dh2dh3dh4dh5dh6_basic[2][i]);
+		printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh6[0][i],dh1dh6[1][i],dh1dh6[2][i]);
 	}
 #endif
@@ -2775,18 +2775,18 @@
 	/*Build BPrime: */
 	for (i=0;i<numgrids;i++){
-		*(B+NDOF2*numgrids*0+NDOF2*i)=2.0*dh1dh2dh3dh4dh5dh6_basic[0][i]; 
-		*(B+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh2dh3dh4dh5dh6_basic[1][i];
-
-		*(B+NDOF2*numgrids*1+NDOF2*i)=dh1dh2dh3dh4dh5dh6_basic[0][i];
-		*(B+NDOF2*numgrids*1+NDOF2*i+1)=2.0*dh1dh2dh3dh4dh5dh6_basic[1][i];
-
-		*(B+NDOF2*numgrids*2+NDOF2*i)=dh1dh2dh3dh4dh5dh6_basic[1][i]; 
-		*(B+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh2dh3dh4dh5dh6_basic[0][i]; 
-
-		*(B+NDOF2*numgrids*3+NDOF2*i)=dh1dh2dh3dh4dh5dh6_basic[2][i]; 
+		*(B+NDOF2*numgrids*0+NDOF2*i)=2.0*dh1dh6[0][i]; 
+		*(B+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh6[1][i];
+
+		*(B+NDOF2*numgrids*1+NDOF2*i)=dh1dh6[0][i];
+		*(B+NDOF2*numgrids*1+NDOF2*i+1)=2.0*dh1dh6[1][i];
+
+		*(B+NDOF2*numgrids*2+NDOF2*i)=dh1dh6[1][i]; 
+		*(B+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh6[0][i]; 
+
+		*(B+NDOF2*numgrids*3+NDOF2*i)=dh1dh6[2][i]; 
 		*(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
 
 		*(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
-		*(B+NDOF2*numgrids*4+NDOF2*i+1)=dh1dh2dh3dh4dh5dh6_basic[2][i]; 
+		*(B+NDOF2*numgrids*4+NDOF2*i+1)=dh1dh6[2][i]; 
 	}
 }
@@ -2799,5 +2799,5 @@
 
 	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 
-	 * For grid i, Bi' can be expressed in the basic coordinate system
+	 * For grid i, Bi' can be expressed in the actual coordinate system
 	 * by: 
 	 *       Biprime_advec=[ dh/dx ]
@@ -2814,15 +2814,15 @@
 	int DOFPERGRID=1;
 
-	/*Same thing in the basic coordinate system: */
-	double dh1dh6_basic[calculationdof][numgrids];
-
-	/*Get dh1dh2dh3 in basic coordinates system : */
-	GetNodalFunctionsDerivativesBasic(&dh1dh6_basic[0][0],xyz_list,gauss_coord);
+	/*Same thing in the actual coordinate system: */
+	double dh1dh6[calculationdof][numgrids];
+
+	/*Get dh1dh2dh3 in actual coordinates system : */
+	GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list,gauss_coord);
 
 	/*Build B': */
 	for (i=0;i<numgrids;i++){
-		*(Bprime_advec+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6_basic[0][i]; 
-		*(Bprime_advec+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6_basic[1][i]; 
-		*(Bprime_advec+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6_basic[2][i]; 
+		*(Bprime_advec+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i]; 
+		*(Bprime_advec+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i]; 
+		*(Bprime_advec+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6[2][i]; 
 	}
 }
@@ -2831,5 +2831,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta:GetBPrime_vert"
-void Penta::GetBPrime_vert(double* B, double* xyz_list, double* gauss_l1l2l3l4){
+void Penta::GetBPrime_vert(double* B, double* xyz_list, double* gauss_coord){
 
 	// Compute Bprime  matrix. Bprime=[L1 L2 L3 L4 L5 L6] where Li is the nodal function for grid i
@@ -2837,5 +2837,5 @@
 	int i;
 
-	GetNodalFunctions(B, gauss_l1l2l3l4);
+	GetNodalFunctions(B, gauss_coord);
 
 }
@@ -2847,7 +2847,7 @@
 
 	/*	Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. 
-	 *	For grid i, Bi' can be expressed in the basic coordinate system
+	 *	For grid i, Bi' can be expressed in the actual coordinate system
 	 *	by: 
-	 *				Bi_basic'=[ dh/dx   0          0       0]
+	 *				Bi'=[ dh/dx   0          0       0]
 	 *					 [   0      dh/dy      0       0]
 	 *					 [   0      0         dh/dz    0]
@@ -2867,9 +2867,9 @@
 	int DOFPERGRID=4;
 
-	double dh1dh7_basic[calculationdof][numgrids+1];
+	double dh1dh7[calculationdof][numgrids+1];
 	double l1l6[numgrids];
 
-	/*Get dh1dh7 in basic coordinate system: */
-	GetNodalFunctionsDerivativesBasicStokes(&dh1dh7_basic[0][0],xyz_list, gauss_coord);
+	/*Get dh1dh7 in actual coordinate system: */
+	GetNodalFunctionsDerivativesStokes(&dh1dh7[0][0],xyz_list, gauss_coord);
 
 	GetNodalFunctions(l1l6, gauss_coord);
@@ -2877,5 +2877,5 @@
 #ifdef _ISSM_DEBUG_ 
 	for (i=0;i<6;i++){
-		printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7_basic[0][i],dh1dh7_basic[1][i]);
+		printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7[0][i],dh1dh7[1][i]);
 	}
 
@@ -2884,25 +2884,25 @@
 	/*B_primeuild B_prime: */
 	for (i=0;i<numgrids+1;i++){
-		*(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7_basic[0][i]; //B_prime[0][DOFPERGRID*i]=dh1dh6_basic[0][i];
+		*(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B_prime[0][DOFPERGRID*i]=dh1dh6[0][i];
 		*(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;
 		*(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;
 		*(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7_basic[1][i];
+		*(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i];
 		*(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;
 		*(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;
 		*(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7_basic[2][i];
-		*(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=dh1dh7_basic[1][i]; 
-		*(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=dh1dh7_basic[0][i]; 
+		*(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i];
+		*(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=dh1dh7[1][i]; 
+		*(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=dh1dh7[0][i]; 
 		*(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=dh1dh7_basic[2][i];
+		*(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=dh1dh7[2][i];
 		*(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=dh1dh7_basic[0][i];
+		*(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=dh1dh7[0][i];
 		*(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;
-		*(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=dh1dh7_basic[2][i];
-		*(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=dh1dh7_basic[1][i];
-		*(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=dh1dh7_basic[0][i];
-		*(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=dh1dh7_basic[1][i];
-		*(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=dh1dh7_basic[2][i];
+		*(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=dh1dh7[2][i];
+		*(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=dh1dh7[1][i];
+		*(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=dh1dh7[0][i];
+		*(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=dh1dh7[1][i];
+		*(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=dh1dh7[2][i];
 		*(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=0;
 		*(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=0;
@@ -2929,6 +2929,6 @@
 
 	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*DOFPERGRID. 
-	 * For grid i, Bi can be expressed in the basic coordinate system
-	 * by: 				Bi_basic=[ dh/dx          0             0       0  ]
+	 * For grid i, Bi can be expressed in the actual coordinate system
+	 * by: 				Bi=[ dh/dx          0             0       0  ]
 	 *					[   0           dh/dy           0       0  ]
 	 *					[   0             0           dh/dy     0  ]
@@ -2947,10 +2947,10 @@
 	int DOFPERGRID=4;
 
-	double dh1dh7_basic[calculationdof][numgrids+1];
+	double dh1dh7[calculationdof][numgrids+1];
 	double l1l6[numgrids];
 
 
-	/*Get dh1dh7 in basic coordinate system: */
-	GetNodalFunctionsDerivativesBasicStokes(&dh1dh7_basic[0][0],xyz_list, gauss_coord);
+	/*Get dh1dh7 in actual coordinate system: */
+	GetNodalFunctionsDerivativesStokes(&dh1dh7[0][0],xyz_list, gauss_coord);
 
 	GetNodalFunctions(l1l6, gauss_coord);
@@ -2958,5 +2958,5 @@
 #ifdef _ISSM_DEBUG_ 
 	for (i=0;i<7;i++){
-		printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7_basic[0][i],dh1dh7_basic[1][i],dh1dh7_basic[2][i]);
+		printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7[0][i],dh1dh7[1][i],dh1dh7[2][i]);
 	}
 
@@ -2965,28 +2965,28 @@
 	/*Build B: */
 	for (i=0;i<numgrids+1;i++){
-		*(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7_basic[0][i]; //B[0][DOFPERGRID*i]=dh1dh6_basic[0][i];
+		*(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B[0][DOFPERGRID*i]=dh1dh6[0][i];
 		*(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;
 		*(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;
 		*(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;
-		*(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7_basic[1][i];
+		*(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i];
 		*(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;
 		*(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;
 		*(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;
-		*(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7_basic[2][i];
-		*(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=(float).5*dh1dh7_basic[1][i]; 
-		*(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=(float).5*dh1dh7_basic[0][i]; 
+		*(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i];
+		*(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=(float).5*dh1dh7[1][i]; 
+		*(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=(float).5*dh1dh7[0][i]; 
 		*(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;
-		*(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=(float).5*dh1dh7_basic[2][i];
+		*(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=(float).5*dh1dh7[2][i];
 		*(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;
-		*(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=(float).5*dh1dh7_basic[0][i];
+		*(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=(float).5*dh1dh7[0][i];
 		*(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;
-		*(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=(float).5*dh1dh7_basic[2][i];
-		*(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=(float).5*dh1dh7_basic[1][i];
+		*(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=(float).5*dh1dh7[2][i];
+		*(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=(float).5*dh1dh7[1][i];
 		*(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=0;
 		*(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=0;
 		*(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=0;
-		*(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=dh1dh7_basic[0][i];
-		*(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=dh1dh7_basic[1][i];
-		*(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=dh1dh7_basic[2][i];
+		*(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=dh1dh7[0][i];
+		*(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=dh1dh7[1][i];
+		*(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=dh1dh7[2][i];
 	}
 
@@ -3041,5 +3041,5 @@
 #undef __FUNCT__
 #define __FUNCT__ "Penta::GetJacobian" 
-void Penta::GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3l4){
+void Penta::GetJacobian(double* J, double* xyz_list,double* gauss_coord){
 
 	const int NDOF3=3;
@@ -3059,11 +3059,11 @@
 
 	/*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */
-	A1=gauss_l1l2l3l4[0];
-	A2=gauss_l1l2l3l4[1];
-	A3=gauss_l1l2l3l4[2];
+	A1=gauss_coord[0];
+	A2=gauss_coord[1];
+	A3=gauss_coord[2];
 
 	xi=A2-A1;
 	eta=sqrt3*A3;
-	zi=gauss_l1l2l3l4[3];
+	zi=gauss_coord[3];
 
 	x1=*(xyz_list+3*0+0);
@@ -3114,5 +3114,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::GetJacobianDeterminant" 
-void Penta::GetJacobianDeterminant(double*  Jdet, double* xyz_list,double* gauss_l1l2l3l4){
+void Penta::GetJacobianDeterminant(double*  Jdet, double* xyz_list,double* gauss_coord){
 
 	/*On a penta, Jacobian varies according to coordinates. We need to get the Jacobian, and take 
@@ -3122,5 +3122,5 @@
 	double J[NDOF3][NDOF3];
 
-	GetJacobian(&J[0][0],xyz_list,gauss_l1l2l3l4);
+	GetJacobian(&J[0][0],xyz_list,gauss_coord);
 
 	*Jdet= J[0][0]*J[1][1]*J[2][2]-J[0][0]*J[1][2]*J[2][1]-J[1][0]*J[0][1]*J[2][2]+J[1][0]*J[0][2]*J[2][1]+J[2][0]*J[0][1]*J[1][2]-J[2][0]*J[0][2]*J[1][1];
@@ -3137,7 +3137,7 @@
 	/*
 	 * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
-	 * For grid i, Li can be expressed in the basic coordinate system
+	 * For grid i, Li can be expressed in the actual coordinate system
 	 * by: 
-	 *       Li_basic=[ h    0    0   0]
+	 *       Li=[ h    0    0   0]
 	 *	 	 [ 0    h    0   0]
 	 *		 [ 0    0    h   0]
@@ -3163,5 +3163,5 @@
 
 
-	/*Get l1l2l3 in basic coordinate system: */
+	/*Get l1l2l3 in actual coordinate system: */
 	l1l2l3[0]=gauss_coord_tria[0];
 	l1l2l3[1]=gauss_coord_tria[1];
@@ -3177,5 +3177,5 @@
 	/*Build LStokes: */
 	for (i=0;i<3;i++){
-		*(LStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh2dh3_basic[0][i];
+		*(LStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
 		*(LStokes+num_dof*numgrids2d*0+num_dof*i+1)=0;
 		*(LStokes+num_dof*numgrids2d*0+num_dof*i+2)=0;
@@ -3244,7 +3244,7 @@
 	/*
 	 * Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 
-	 * For grid i, Lpi can be expressed in the basic coordinate system
+	 * For grid i, Lpi can be expressed in the actual coordinate system
 	 * by: 
-	 *       Lpi_basic=[ h    0    0   0]
+	 *       Lpi=[ h    0    0   0]
 	 *		 [ 0    h    0   0]
 	 *		 [ h    0    0   0]
@@ -3268,13 +3268,13 @@
 
 	double l1l2l3[numgrids2d];
-	double dh1dh6_basic[3][6];
-
-
-	/*Get l1l2l3 in basic coordinate system: */
+	double dh1dh6[3][6];
+
+
+	/*Get l1l2l3 in actual coordinate system: */
 	l1l2l3[0]=gauss_coord_tria[0];
 	l1l2l3[1]=gauss_coord_tria[1];
 	l1l2l3[2]=gauss_coord_tria[2];
 
-	GetNodalFunctionsDerivativesBasic(&dh1dh6_basic[0][0],xyz_list,gauss_coord);
+	GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list,gauss_coord);
 
 #ifdef _DELUG_ 
@@ -3286,5 +3286,5 @@
 	/*Build LprimeStokes: */
 	for (i=0;i<3;i++){
-		*(LprimeStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh2dh3_basic[0][i];
+		*(LprimeStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
 		*(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+1)=0;
 		*(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+2)=0;
@@ -3312,21 +3312,21 @@
 		*(LprimeStokes+num_dof*numgrids2d*6+num_dof*i)=0;
 		*(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+2)=dh1dh6_basic[2][i];
+		*(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+2)=dh1dh6[2][i];
 		*(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+3)=0;
 		*(LprimeStokes+num_dof*numgrids2d*7+num_dof*i)=0;
 		*(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+2)=dh1dh6_basic[2][i];
+		*(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+2)=dh1dh6[2][i];
 		*(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+3)=0;
 		*(LprimeStokes+num_dof*numgrids2d*8+num_dof*i)=0;
 		*(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+2)=dh1dh6_basic[2][i];
+		*(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+2)=dh1dh6[2][i];
 		*(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+3)=0;
-		*(LprimeStokes+num_dof*numgrids2d*9+num_dof*i)=dh1dh6_basic[2][i];
+		*(LprimeStokes+num_dof*numgrids2d*9+num_dof*i)=dh1dh6[2][i];
 		*(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+1)=0;
-		*(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+2)=dh1dh6_basic[0][i];
+		*(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+2)=dh1dh6[0][i];
 		*(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+3)=0;
 		*(LprimeStokes+num_dof*numgrids2d*10+num_dof*i)=0;
-		*(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+1)=dh1dh6_basic[2][i];
-		*(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+2)=dh1dh6_basic[1][i];
+		*(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+1)=dh1dh6[2][i];
+		*(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+2)=dh1dh6[1][i];
 		*(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+3)=0;
 		*(LprimeStokes+num_dof*numgrids2d*11+num_dof*i)=0;
@@ -3386,10 +3386,10 @@
 }
 /*}}}*/
-/*FUNCTION GetNodalFunctionsDerivativesBasic {{{1*/
+/*FUNCTION GetNodalFunctionsDerivatives {{{1*/
 #undef __FUNCT__
-#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesBasic" 
-void Penta::GetNodalFunctionsDerivativesBasic(double* dh1dh2dh3dh4dh5dh6_basic,double* xyz_list, double* gauss_l1l2l3l4){
-
-	/*This routine returns the values of the nodal functions derivatives  (with respect to the basic coordinate system: */
+#define __FUNCT__ "Penta::GetNodalFunctionsDerivatives" 
+void Penta::GetNodalFunctionsDerivatives(double* dh1dh6,double* xyz_list, double* gauss_coord){
+
+	/*This routine returns the values of the nodal functions derivatives  (with respect to the actual coordinate system: */
 
 
@@ -3398,14 +3398,14 @@
 	const int numgrids=6;
 
-	double dh1dh2dh3dh4dh5dh6_param[NDOF3][numgrids];
+	double dh1dh6_ref[NDOF3][numgrids];
 	double Jinv[NDOF3][NDOF3];
 
 	/*Get derivative values with respect to parametric coordinate system: */
-	GetNodalFunctionsDerivativesParams(&dh1dh2dh3dh4dh5dh6_param[0][0], gauss_l1l2l3l4); 
+	GetNodalFunctionsDerivativesReference(&dh1dh6_ref[0][0], gauss_coord); 
 
 	/*Get Jacobian invert: */
-	GetJacobianInvert(&Jinv[0][0], xyz_list, gauss_l1l2l3l4);
-
-	/*Build dh1dh2dh3_basic: 
+	GetJacobianInvert(&Jinv[0][0], xyz_list, gauss_coord);
+
+	/*Build dh1dh3: 
 	 *
 	 * [dhi/dx]= Jinv*[dhi/dr]
@@ -3415,33 +3415,33 @@
 
 	for (i=0;i<numgrids;i++){
-		*(dh1dh2dh3dh4dh5dh6_basic+numgrids*0+i)=Jinv[0][0]*dh1dh2dh3dh4dh5dh6_param[0][i]+Jinv[0][1]*dh1dh2dh3dh4dh5dh6_param[1][i]+Jinv[0][2]*dh1dh2dh3dh4dh5dh6_param[2][i];
-		*(dh1dh2dh3dh4dh5dh6_basic+numgrids*1+i)=Jinv[1][0]*dh1dh2dh3dh4dh5dh6_param[0][i]+Jinv[1][1]*dh1dh2dh3dh4dh5dh6_param[1][i]+Jinv[1][2]*dh1dh2dh3dh4dh5dh6_param[2][i];
-		*(dh1dh2dh3dh4dh5dh6_basic+numgrids*2+i)=Jinv[2][0]*dh1dh2dh3dh4dh5dh6_param[0][i]+Jinv[2][1]*dh1dh2dh3dh4dh5dh6_param[1][i]+Jinv[2][2]*dh1dh2dh3dh4dh5dh6_param[2][i];
-	}
-
-}
-/*}}}*/
-/*FUNCTION GetNodalFunctionsDerivativesBasicStokes {{{1*/
-#undef __FUNCT__ 
-#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesBasicStokes"
-void Penta::GetNodalFunctionsDerivativesBasicStokes(double* dh1dh7_basic,double* xyz_list, double* gauss_coord){
+		*(dh1dh6+numgrids*0+i)=Jinv[0][0]*dh1dh6_ref[0][i]+Jinv[0][1]*dh1dh6_ref[1][i]+Jinv[0][2]*dh1dh6_ref[2][i];
+		*(dh1dh6+numgrids*1+i)=Jinv[1][0]*dh1dh6_ref[0][i]+Jinv[1][1]*dh1dh6_ref[1][i]+Jinv[1][2]*dh1dh6_ref[2][i];
+		*(dh1dh6+numgrids*2+i)=Jinv[2][0]*dh1dh6_ref[0][i]+Jinv[2][1]*dh1dh6_ref[1][i]+Jinv[2][2]*dh1dh6_ref[2][i];
+	}
+
+}
+/*}}}*/
+/*FUNCTION GetNodalFunctionsDerivativesStokes {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesStokes"
+void Penta::GetNodalFunctionsDerivativesStokes(double* dh1dh7,double* xyz_list, double* gauss_coord){
 
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
-	 * basic coordinate system: */
+	 * actual coordinate system: */
 
 	int i;
 
 	const  int numgrids=7;
-	double dh1dh7_param[3][numgrids];
+	double dh1dh7_ref[3][numgrids];
 	double Jinv[3][3];
 
 
 	/*Get derivative values with respect to parametric coordinate system: */
-	GetNodalFunctionsDerivativesParamsStokes(&dh1dh7_param[0][0], gauss_coord); 
+	GetNodalFunctionsDerivativesReferenceStokes(&dh1dh7_ref[0][0], gauss_coord); 
 
 	/*Get Jacobian invert: */
 	GetJacobianInvert(&Jinv[0][0], xyz_list, gauss_coord);
 
-	/*Build dh1dh6_basic: 
+	/*Build dh1dh6: 
 	 *
 	 * [dhi/dx]= Jinv'*[dhi/dr]
@@ -3451,15 +3451,15 @@
 
 	for (i=0;i<numgrids;i++){
-		*(dh1dh7_basic+numgrids*0+i)=Jinv[0][0]*dh1dh7_param[0][i]+Jinv[0][1]*dh1dh7_param[1][i]+Jinv[0][2]*dh1dh7_param[2][i];
-		*(dh1dh7_basic+numgrids*1+i)=Jinv[1][0]*dh1dh7_param[0][i]+Jinv[1][1]*dh1dh7_param[1][i]+Jinv[1][2]*dh1dh7_param[2][i];
-		*(dh1dh7_basic+numgrids*2+i)=Jinv[2][0]*dh1dh7_param[0][i]+Jinv[2][1]*dh1dh7_param[1][i]+Jinv[2][2]*dh1dh7_param[2][i];
-	}
-
-}
-/*}}}*/
-/*FUNCTION GetNodalFunctionsDerivativesParams {{{1*/
+		*(dh1dh7+numgrids*0+i)=Jinv[0][0]*dh1dh7_ref[0][i]+Jinv[0][1]*dh1dh7_ref[1][i]+Jinv[0][2]*dh1dh7_ref[2][i];
+		*(dh1dh7+numgrids*1+i)=Jinv[1][0]*dh1dh7_ref[0][i]+Jinv[1][1]*dh1dh7_ref[1][i]+Jinv[1][2]*dh1dh7_ref[2][i];
+		*(dh1dh7+numgrids*2+i)=Jinv[2][0]*dh1dh7_ref[0][i]+Jinv[2][1]*dh1dh7_ref[1][i]+Jinv[2][2]*dh1dh7_ref[2][i];
+	}
+
+}
+/*}}}*/
+/*FUNCTION GetNodalFunctionsDerivativesReference {{{1*/
 #undef __FUNCT__
-#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesParams" 
-void Penta::GetNodalFunctionsDerivativesParams(double* dl1dl2dl3dl4dl5dl6,double* gauss_l1l2l3l4){
+#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesReference" 
+void Penta::GetNodalFunctionsDerivativesReference(double* dl1dl6,double* gauss_coord){
 
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
@@ -3470,45 +3470,45 @@
 	double sqrt3=sqrt(3.0);
 
-	A1=gauss_l1l2l3l4[0]; //first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*sqrt(3));
-	A2=gauss_l1l2l3l4[1]; //second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*sqrt(3));
-	A3=gauss_l1l2l3l4[2]; //third area coordinate value  In term of xi and eta: A3=y/sqrt(3);
-	z=gauss_l1l2l3l4[3]; //fourth vertical coordinate value. Corresponding nodal function: (1-z)/2 and (1+z)/2
+	A1=gauss_coord[0]; //first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*sqrt(3));
+	A2=gauss_coord[1]; //second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*sqrt(3));
+	A3=gauss_coord[2]; //third area coordinate value  In term of xi and eta: A3=y/sqrt(3);
+	z=gauss_coord[3]; //fourth vertical coordinate value. Corresponding nodal function: (1-z)/2 and (1+z)/2
 
 
 	/*First nodal function derivatives. The corresponding nodal function is N=A1*(1-z)/2. Its derivatives follow*/
-	*(dl1dl2dl3dl4dl5dl6+numgrids*0+0)=-1.0/2.0*(1.0-z)/2.0;
-	*(dl1dl2dl3dl4dl5dl6+numgrids*1+0)=-1.0/2.0/sqrt3*(1.0-z)/2.0;
-	*(dl1dl2dl3dl4dl5dl6+numgrids*2+0)=-1.0/2.0*A1;
+	*(dl1dl6+numgrids*0+0)=-1.0/2.0*(1.0-z)/2.0;
+	*(dl1dl6+numgrids*1+0)=-1.0/2.0/sqrt3*(1.0-z)/2.0;
+	*(dl1dl6+numgrids*2+0)=-1.0/2.0*A1;
 
 	/*Second nodal function: The corresponding nodal function is N=A2*(1-z)/2. Its derivatives follow*/
-	*(dl1dl2dl3dl4dl5dl6+numgrids*0+1)=1.0/2.0*(1.0-z)/2.0;
-	*(dl1dl2dl3dl4dl5dl6+numgrids*1+1)=-1.0/2.0/sqrt3*(1.0-z)/2.0;
-	*(dl1dl2dl3dl4dl5dl6+numgrids*2+1)=-1.0/2.0*A2;
+	*(dl1dl6+numgrids*0+1)=1.0/2.0*(1.0-z)/2.0;
+	*(dl1dl6+numgrids*1+1)=-1.0/2.0/sqrt3*(1.0-z)/2.0;
+	*(dl1dl6+numgrids*2+1)=-1.0/2.0*A2;
 
 	/*Third nodal function: The corresponding nodal function is N=A3*(1-z)/2. Its derivatives follow*/
-	*(dl1dl2dl3dl4dl5dl6+numgrids*0+2)=0.0;
-	*(dl1dl2dl3dl4dl5dl6+numgrids*1+2)=1.0/sqrt3*(1.0-z)/2.0;
-	*(dl1dl2dl3dl4dl5dl6+numgrids*2+2)=-1.0/2.0*A3;
+	*(dl1dl6+numgrids*0+2)=0.0;
+	*(dl1dl6+numgrids*1+2)=1.0/sqrt3*(1.0-z)/2.0;
+	*(dl1dl6+numgrids*2+2)=-1.0/2.0*A3;
 
 	/*Fourth nodal function: The corresponding nodal function is N=A1*(1+z)/2. Its derivatives follow*/
-	*(dl1dl2dl3dl4dl5dl6+numgrids*0+3)=-1.0/2.0*(1.0+z)/2.0;
-	*(dl1dl2dl3dl4dl5dl6+numgrids*1+3)=-1.0/2.0/sqrt3*(1.0+z)/2.0;
-	*(dl1dl2dl3dl4dl5dl6+numgrids*2+3)=1.0/2.0*A1;
+	*(dl1dl6+numgrids*0+3)=-1.0/2.0*(1.0+z)/2.0;
+	*(dl1dl6+numgrids*1+3)=-1.0/2.0/sqrt3*(1.0+z)/2.0;
+	*(dl1dl6+numgrids*2+3)=1.0/2.0*A1;
 
 	/*Fifth nodal function: The corresponding nodal function is N=A2*(1+z)/2. Its derivatives follow*/
-	*(dl1dl2dl3dl4dl5dl6+numgrids*0+4)=1.0/2.0*(1.0+z)/2.0;
-	*(dl1dl2dl3dl4dl5dl6+numgrids*1+4)=-1.0/2.0/sqrt3*(1.0+z)/2.0;
-	*(dl1dl2dl3dl4dl5dl6+numgrids*2+4)=1.0/2.0*A2;
+	*(dl1dl6+numgrids*0+4)=1.0/2.0*(1.0+z)/2.0;
+	*(dl1dl6+numgrids*1+4)=-1.0/2.0/sqrt3*(1.0+z)/2.0;
+	*(dl1dl6+numgrids*2+4)=1.0/2.0*A2;
 
 	/*Sixth nodal function: The corresponding nodal function is N=A3*(1+z)/2. Its derivatives follow*/
-	*(dl1dl2dl3dl4dl5dl6+numgrids*0+5)=0.0;
-	*(dl1dl2dl3dl4dl5dl6+numgrids*1+5)=1.0/sqrt3*(1.0+z)/2.0;
-	*(dl1dl2dl3dl4dl5dl6+numgrids*2+5)=1.0/2.0*A3;
-}
-/*}}}*/
-/*FUNCTION GetNodalFunctionsDerivativesParamsStokes {{{1*/
-#undef __FUNCT__ 
-#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesParamsStokes"
-void Penta::GetNodalFunctionsDerivativesParamsStokes(double* dl1dl7,double* gauss_coord){
+	*(dl1dl6+numgrids*0+5)=0.0;
+	*(dl1dl6+numgrids*1+5)=1.0/sqrt3*(1.0+z)/2.0;
+	*(dl1dl6+numgrids*2+5)=1.0/2.0*A3;
+}
+/*}}}*/
+/*FUNCTION GetNodalFunctionsDerivativesReferenceStokes {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesReferenceStokes"
+void Penta::GetNodalFunctionsDerivativesReferenceStokes(double* dl1dl7,double* gauss_coord){
 
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
@@ -3592,7 +3592,7 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::GetParameterDerivativeValue"
-void Penta::GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4){
-
-	/*From grid 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_l1l2l3l4:
+void Penta::GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_coord){
+
+	/*From grid 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;
@@ -3604,14 +3604,14 @@
 	const int NDOF3=3;
 	const int numgrids=6;
-	double dh1dh2dh3dh4dh5dh6_basic[NDOF3][numgrids];
-
-	/*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */
-	GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
-
-	*(p+0)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[0][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[0][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[0][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[0][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[0][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[0][5];
+	double dh1dh6[NDOF3][numgrids];
+
+	/*Get dh1dh6 in actual coordinate system: */
+	GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list, gauss_coord);
+
+	*(p+0)=p_list[0]*dh1dh6[0][0]+p_list[1]*dh1dh6[0][1]+p_list[2]*dh1dh6[0][2]+p_list[3]*dh1dh6[0][3]+p_list[4]*dh1dh6[0][4]+p_list[5]*dh1dh6[0][5];
 	;
-	*(p+1)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[1][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[1][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[1][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[1][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[1][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[1][5];
-
-	*(p+2)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[2][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[2][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[2][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[2][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[2][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[2][5];
+	*(p+1)=p_list[0]*dh1dh6[1][0]+p_list[1]*dh1dh6[1][1]+p_list[2]*dh1dh6[1][2]+p_list[3]*dh1dh6[1][3]+p_list[4]*dh1dh6[1][4]+p_list[5]*dh1dh6[1][5];
+
+	*(p+2)=p_list[0]*dh1dh6[2][0]+p_list[1]*dh1dh6[2][1]+p_list[2]*dh1dh6[2][2]+p_list[3]*dh1dh6[2][3]+p_list[4]*dh1dh6[2][4]+p_list[5]*dh1dh6[2][5];
 
 }
@@ -3620,12 +3620,12 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::GetParameterValue"
-void Penta::GetParameterValue(double* pvalue, double* v_list,double* gauss_l1l2l3l4){
+void Penta::GetParameterValue(double* pvalue, double* v_list,double* gauss_coord){
 
 	const int numgrids=6;
-	double l1l2l3l4l5l6[numgrids];
-
-	GetNodalFunctions(&l1l2l3l4l5l6[0], gauss_l1l2l3l4);
-
-	*pvalue=l1l2l3l4l5l6[0]*v_list[0]+l1l2l3l4l5l6[1]*v_list[1]+l1l2l3l4l5l6[2]*v_list[2]+l1l2l3l4l5l6[3]*v_list[3]+l1l2l3l4l5l6[4]*v_list[4]+l1l2l3l4l5l6[5]*v_list[5];
+	double l1l6[numgrids];
+
+	GetNodalFunctions(&l1l6[0], gauss_coord);
+
+	*pvalue=l1l6[0]*v_list[0]+l1l6[1]*v_list[1]+l1l6[2]*v_list[2]+l1l6[3]*v_list[3]+l1l6[4]*v_list[4]+l1l6[5]*v_list[5];
 }
 /*}}}*/
@@ -3633,5 +3633,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::GetJacobianInvert"
-void Penta::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_l1l2l3l4){
+void Penta::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_coord){
 
 	double Jdet;
@@ -3639,5 +3639,5 @@
 
 	/*Call Jacobian routine to get the jacobian:*/
-	GetJacobian(Jinv, xyz_list, gauss_l1l2l3l4);
+	GetJacobian(Jinv, xyz_list, gauss_coord);
 
 	/*Invert Jacobian matrix: */
@@ -3653,19 +3653,19 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::GetNodalFunctions"
-void Penta::GetNodalFunctions(double* l1l2l3l4l5l6, double* gauss_l1l2l3l4){
+void Penta::GetNodalFunctions(double* l1l6, double* gauss_coord){
 
 	/*This routine returns the values of the nodal functions  at the gaussian point.*/
 
-	l1l2l3l4l5l6[0]=gauss_l1l2l3l4[0]*(1-gauss_l1l2l3l4[3])/2.0;
-
-	l1l2l3l4l5l6[1]=gauss_l1l2l3l4[1]*(1-gauss_l1l2l3l4[3])/2.0;
-
-	l1l2l3l4l5l6[2]=gauss_l1l2l3l4[2]*(1-gauss_l1l2l3l4[3])/2.0;
-
-	l1l2l3l4l5l6[3]=gauss_l1l2l3l4[0]*(1+gauss_l1l2l3l4[3])/2.0;
-
-	l1l2l3l4l5l6[4]=gauss_l1l2l3l4[1]*(1+gauss_l1l2l3l4[3])/2.0;
-
-	l1l2l3l4l5l6[5]=gauss_l1l2l3l4[2]*(1+gauss_l1l2l3l4[3])/2.0;
+	l1l6[0]=gauss_coord[0]*(1-gauss_coord[3])/2.0;
+
+	l1l6[1]=gauss_coord[1]*(1-gauss_coord[3])/2.0;
+
+	l1l6[2]=gauss_coord[2]*(1-gauss_coord[3])/2.0;
+
+	l1l6[3]=gauss_coord[0]*(1+gauss_coord[3])/2.0;
+
+	l1l6[4]=gauss_coord[1]*(1+gauss_coord[3])/2.0;
+
+	l1l6[5]=gauss_coord[2]*(1+gauss_coord[3])/2.0;
 
 }
@@ -3730,5 +3730,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::GetStrainRate"
-void Penta::GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_l1l2l3l4){
+void Penta::GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord){
 
 	int i;
@@ -3739,5 +3739,5 @@
 
 	/*Get B matrix: */
-	GetB(&B[0][0], xyz_list, gauss_l1l2l3l4);
+	GetB(&B[0][0], xyz_list, gauss_coord);
 
 #ifdef _ISSM_DEBUG_
Index: /issm/trunk/src/c/objects/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Penta.h	(revision 2955)
+++ /issm/trunk/src/c/objects/Penta.h	(revision 2956)
@@ -95,19 +95,19 @@
 		void*  SpawnTria(int g0, int g1, int g2);
 
-		void  GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_l1l2l3l4);
-		void  GetB(double* pB, double* xyz_list, double* gauss_l1l2l3l4);
-		void  GetBPrime(double* B, double* xyz_list, double* gauss_l1l2l3l4);
-		void  GetB_vert(double* B, double* xyz_list, double* gauss_l1l2l3l4);
-		void  GetBPrime_vert(double* B, double* xyz_list, double* gauss_l1l2l3l4);
-		void  GetJacobianDeterminant(double*  Jdet, double* xyz_list,double* gauss_l1l2l3l4);
-		void  GetNodalFunctionsDerivativesBasic(double* dh1dh2dh3dh4dh5dh6_basic,double* xyz_list, double* gauss_l1l2l3l4);
-		void  GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3l4);
-		void  GetNodalFunctionsDerivativesParams(double* dl1dl2dl3dl4dl5dl6,double* gauss_l1l2l3l4);
-		void  GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_l1l2l3l4);
+		void  GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord);
+		void  GetB(double* pB, double* xyz_list, double* gauss_coord);
+		void  GetBPrime(double* B, double* xyz_list, double* gauss_coord);
+		void  GetB_vert(double* B, double* xyz_list, double* gauss_coord);
+		void  GetBPrime_vert(double* B, double* xyz_list, double* gauss_coord);
+		void  GetJacobianDeterminant(double*  Jdet, double* xyz_list,double* gauss_coord);
+		void  GetNodalFunctionsDerivatives(double* dh1dh6,double* xyz_list, double* gauss_coord);
+		void  GetJacobian(double* J, double* xyz_list,double* gauss_coord);
+		void  GetNodalFunctionsDerivativesReference(double* dl1dl6,double* gauss_coord);
+		void  GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_coord);
 		void  CreatePVectorDiagnosticHoriz( Vec pg, void* inputs,int analysis_type,int sub_analysis_type);
 		void  CreatePVectorDiagnosticVert( Vec pg, void* inputs,int analysis_type,int sub_analysis_type);
-		void  GetParameterValue(double* pvalue, double* v_list,double* gauss_l1l2l3l4);
-		void  GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4);
-		void  GetNodalFunctions(double* l1l2l3l4l5l6, double* gauss_l1l2l3l4);
+		void  GetParameterValue(double* pvalue, double* v_list,double* gauss_coord);
+		void  GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_coord);
+		void  GetNodalFunctions(double* l1l6, double* gauss_coord);
 		void  FieldExtrude(Vec field,double* field_serial,char* field_name, int iscollapsed);
 		void  ComputePressure(Vec p_g);
@@ -131,6 +131,6 @@
 		void  GetLStokes(double* LStokes, double* gauss_coord_tria);
 		void  GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss_coord_tria, double* gauss_coord);
-		void  GetNodalFunctionsDerivativesBasicStokes(double* dh1dh7_basic,double* xyz_list, double* gauss_coord);
-		void  GetNodalFunctionsDerivativesParamsStokes(double* dl1dl7,double* gauss_coord);
+		void  GetNodalFunctionsDerivativesStokes(double* dh1dh7,double* xyz_list, double* gauss_coord);
+		void  GetNodalFunctionsDerivativesReferenceStokes(double* dl1dl7,double* gauss_coord);
 		void  ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp);
 		void  GetNodalFunctionsStokes(double* l1l7, double* gauss_coord);
Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 2955)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 2956)
@@ -3054,7 +3054,7 @@
 
 	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
-	 * For grid i, Bi can be expressed in the basic coordinate system
+	 * For grid i, Bi can be expressed in the actual coordinate system
 	 * by: 
-	 *       Bi_basic=[ dh/dx    0    ]
+	 *       Bi=[ dh/dx    0    ]
 	 *                [   0    dh/dy  ]
 	 *                [ 1/2*dh/dy  1/2*dh/dx  ]
@@ -3068,13 +3068,13 @@
 	const int numgrids=3;
 
-	double dh1dh2dh3_basic[NDOF2][numgrids];
-
-
-	/*Get dh1dh2dh3 in basic coordinate system: */
-	GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],xyz_list, gauss_l1l2l3);
+	double dh1dh3[NDOF2][numgrids];
+
+
+	/*Get dh1dh2dh3 in actual coordinate system: */
+	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list, gauss_l1l2l3);
 
 	#ifdef _ISSM_DEBUG_ 
 	for (i=0;i<3;i++){
-		printf("Node %i  dh/dx=%lf dh/dy=%lf \n",i,dh1dh2dh3_basic[0][i],dh1dh2dh3_basic[1][i]);
+		printf("Node %i  dh/dx=%lf dh/dy=%lf \n",i,dh1dh3[0][i],dh1dh3[1][i]);
 	}
 	#endif
@@ -3082,10 +3082,10 @@
 	/*Build B: */
 	for (i=0;i<numgrids;i++){
-		*(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh2dh3_basic[0][i]; //B[0][NDOF2*i]=dh1dh2dh3_basic[0][i];
+		*(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh3[0][i]; //B[0][NDOF2*i]=dh1dh3[0][i];
 		*(B+NDOF2*numgrids*0+NDOF2*i+1)=0;
 		*(B+NDOF2*numgrids*1+NDOF2*i)=0;
-		*(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh2dh3_basic[1][i];
-		*(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh2dh3_basic[1][i]; 
-		*(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh2dh3_basic[0][i]; 
+		*(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh3[1][i];
+		*(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh3[1][i]; 
+		*(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh3[0][i]; 
 	}
 }
@@ -3098,7 +3098,7 @@
 
 	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
-	 * For grid i, Bi can be expressed in the basic coordinate system
+	 * For grid i, Bi can be expressed in the actual coordinate system
 	 * by: 
-	 *       Bi_basic=[ h ]
+	 *       Bi=[ h ]
 	 *                [ h ]
 	 * where h is the interpolation function for grid i.
@@ -3114,5 +3114,5 @@
 
 
-	/*Get dh1dh2dh3 in basic coordinate system: */
+	/*Get dh1dh2dh3 in actual coordinate system: */
 	GetNodalFunctions(&l1l2l3[0],gauss_l1l2l3);
 
@@ -3145,7 +3145,7 @@
 
 	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
-	 * For grid i, Bi' can be expressed in the basic coordinate system
+	 * For grid i, Bi' can be expressed in the actual coordinate system
 	 * by: 
-	 *       Bi_prime__basic=[ 2*dh/dx dh/dy ]
+	 *       Bi_prime=[ 2*dh/dx dh/dy ]
 	 *                       [ dh/dx  2*dh/dy]
 	 *                       [dh/dy dh/dx]
@@ -3159,19 +3159,19 @@
 	const int numgrids=3;
 
-	/*Same thing in the basic coordinate system: */
-	double dh1dh2dh3_basic[NDOF2][numgrids];
-
-
-	/*Get dh1dh2dh3 in basic coordinates system : */
-	GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],xyz_list,gauss_l1l2l3);
+	/*Same thing in the actual coordinate system: */
+	double dh1dh3[NDOF2][numgrids];
+
+
+	/*Get dh1dh2dh3 in actual coordinates system : */
+	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss_l1l2l3);
 
 	/*Build B': */
 	for (i=0;i<numgrids;i++){
-		*(Bprime+NDOF2*numgrids*0+NDOF2*i)=2*dh1dh2dh3_basic[0][i]; 
-		*(Bprime+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh2dh3_basic[1][i]; 
-		*(Bprime+NDOF2*numgrids*1+NDOF2*i)=dh1dh2dh3_basic[0][i]; 
-		*(Bprime+NDOF2*numgrids*1+NDOF2*i+1)=2*dh1dh2dh3_basic[1][i]; 
-		*(Bprime+NDOF2*numgrids*2+NDOF2*i)=dh1dh2dh3_basic[1][i]; 
-		*(Bprime+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh2dh3_basic[0][i]; 
+		*(Bprime+NDOF2*numgrids*0+NDOF2*i)=2*dh1dh3[0][i]; 
+		*(Bprime+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh3[1][i]; 
+		*(Bprime+NDOF2*numgrids*1+NDOF2*i)=dh1dh3[0][i]; 
+		*(Bprime+NDOF2*numgrids*1+NDOF2*i+1)=2*dh1dh3[1][i]; 
+		*(Bprime+NDOF2*numgrids*2+NDOF2*i)=dh1dh3[1][i]; 
+		*(Bprime+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh3[0][i]; 
 	}
 }
@@ -3184,7 +3184,7 @@
 
 	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
-	 * For grid i, Bi' can be expressed in the basic coordinate system
+	 * For grid i, Bi' can be expressed in the actual coordinate system
 	 * by: 
-	 *       Bi_prime__basic=[ dh/dx ]
+	 *       Bi_prime=[ dh/dx ]
 	 *                       [ dh/dy ]
 	 * where h is the interpolation function for grid i.
@@ -3198,14 +3198,14 @@
 	const int numgrids=3;
 
-	/*Same thing in the basic coordinate system: */
-	double dh1dh2dh3_basic[NDOF2][numgrids];
-
-	/*Get dh1dh2dh3 in basic coordinates system : */
-	GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],xyz_list,gauss_l1l2l3);
+	/*Same thing in the actual coordinate system: */
+	double dh1dh3[NDOF2][numgrids];
+
+	/*Get dh1dh2dh3 in actual coordinates system : */
+	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss_l1l2l3);
 
 	/*Build B': */
 	for (i=0;i<numgrids;i++){
-		*(Bprime_prog+NDOF1*numgrids*0+NDOF1*i)=dh1dh2dh3_basic[0][i]; 
-		*(Bprime_prog+NDOF1*numgrids*1+NDOF1*i)=dh1dh2dh3_basic[1][i]; 
+		*(Bprime_prog+NDOF1*numgrids*0+NDOF1*i)=dh1dh3[0][i]; 
+		*(Bprime_prog+NDOF1*numgrids*1+NDOF1*i)=dh1dh3[1][i]; 
 	}
 }
@@ -3351,10 +3351,10 @@
 
 	/*Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
-	 * For grid i, Li can be expressed in the basic coordinate system
+	 * For grid i, Li can be expressed in the actual coordinate system
 	 * by: 
 	 *       numdof=1: 
-	 *       Li_basic=h;
+	 *       Li=h;
 	 *       numdof=2:
-	 *       Li_basic=[ h    0    ]
+	 *       Li=[ h    0    ]
 	 *                [   0   h  ]
 	 * where h is the interpolation function for grid i.
@@ -3370,5 +3370,5 @@
 
 
-	/*Get l1l2l3 in basic coordinate system: */
+	/*Get l1l2l3 in actual coordinate system: */
 	GetNodalFunctions(l1l2l3, gauss_l1l2l3);
 
@@ -3387,5 +3387,5 @@
 	else{
 		for (i=0;i<numgrids;i++){
-			*(L+numdof*numgrids*0+numdof*i)=l1l2l3[i]; //L[0][NDOF2*i]=dh1dh2dh3_basic[0][i];
+			*(L+numdof*numgrids*0+numdof*i)=l1l2l3[i]; //L[0][NDOF2*i]=dh1dh3[0][i];
 			*(L+numdof*numgrids*0+numdof*i+1)=0;
 			*(L+numdof*numgrids*1+numdof*i)=0;
@@ -3423,11 +3423,11 @@
 }
 /*}}}*/
-/*FUNCTION GetNodalFunctionsDerivativesBasic {{{1*/
-#undef __FUNCT__ 
-#define __FUNCT__ "Tria::GetNodalFunctionsDerivativesBasic"
-void Tria::GetNodalFunctionsDerivativesBasic(double* dh1dh2dh3_basic,double* xyz_list, double* gauss_l1l2l3){
+/*FUNCTION GetNodalFunctionsDerivatives {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::GetNodalFunctionsDerivatives"
+void Tria::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss_l1l2l3){
 	
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
-	 * basic coordinate system: */
+	 * actual coordinate system: */
 
 	int i;
@@ -3435,15 +3435,15 @@
 	const int numgrids=3;
 
-	double dh1dh2dh3_param[NDOF2][numgrids];
+	double dh1dh3_ref[NDOF2][numgrids];
 	double Jinv[NDOF2][NDOF2];
 
 
 	/*Get derivative values with respect to parametric coordinate system: */
-	GetNodalFunctionsDerivativesParams(&dh1dh2dh3_param[0][0], gauss_l1l2l3); 
+	GetNodalFunctionsDerivativesReference(&dh1dh3_ref[0][0], gauss_l1l2l3); 
 
 	/*Get Jacobian invert: */
 	GetJacobianInvert(&Jinv[0][0], xyz_list, gauss_l1l2l3);
 
-	/*Build dh1dh2dh3_basic: 
+	/*Build dh1dh3: 
 	 *
 	 * [dhi/dx]= Jinv*[dhi/dr]
@@ -3452,14 +3452,14 @@
 
 	for (i=0;i<numgrids;i++){
-		*(dh1dh2dh3_basic+numgrids*0+i)=Jinv[0][0]*dh1dh2dh3_param[0][i]+Jinv[0][1]*dh1dh2dh3_param[1][i];
-		*(dh1dh2dh3_basic+numgrids*1+i)=Jinv[1][0]*dh1dh2dh3_param[0][i]+Jinv[1][1]*dh1dh2dh3_param[1][i];
-	}
-
-}
-/*}}}*/
-/*FUNCTION GetNodalFunctionsDerivativesParams {{{1*/
-#undef __FUNCT__ 
-#define __FUNCT__ "Tria::GetNodalFunctionsDerivativesParams"
-void Tria::GetNodalFunctionsDerivativesParams(double* dl1dl2dl3,double* gauss_l1l2l3){
+		*(dh1dh3+numgrids*0+i)=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i];
+		*(dh1dh3+numgrids*1+i)=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i];
+	}
+
+}
+/*}}}*/
+/*FUNCTION GetNodalFunctionsDerivativesReference {{{1*/
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::GetNodalFunctionsDerivativesReference"
+void Tria::GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss_l1l2l3){
 	
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
@@ -3472,14 +3472,14 @@
 
 	/*First nodal function: */
-	*(dl1dl2dl3+numgrids*0+0)=-1.0/2.0; 
-	*(dl1dl2dl3+numgrids*1+0)=-1.0/(2.0*sqrt3);
+	*(dl1dl3+numgrids*0+0)=-1.0/2.0; 
+	*(dl1dl3+numgrids*1+0)=-1.0/(2.0*sqrt3);
 
 	/*Second nodal function: */
-	*(dl1dl2dl3+numgrids*0+1)=1.0/2.0;
-	*(dl1dl2dl3+numgrids*1+1)=-1.0/(2.0*sqrt3);
+	*(dl1dl3+numgrids*0+1)=1.0/2.0;
+	*(dl1dl3+numgrids*1+1)=-1.0/(2.0*sqrt3);
 
 	/*Third nodal function: */
-	*(dl1dl2dl3+numgrids*0+2)=0;
-	*(dl1dl2dl3+numgrids*1+2)=1.0/sqrt3;
+	*(dl1dl3+numgrids*0+2)=0;
+	*(dl1dl3+numgrids*1+2)=1.0/sqrt3;
 
 }
@@ -3515,11 +3515,11 @@
 	 */
 	
-	double dh1dh2dh3_basic[NDOF2][numgrids]; //nodal derivative functions in basic coordinate system.
-
-	/*Get dh1dh2dh3 in basic coordinate system: */
-	GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],xyz_list, gauss_l1l2l3);
-
-	*(p+0)=plist[0]*dh1dh2dh3_basic[0][0]+plist[1]*dh1dh2dh3_basic[0][1]+plist[2]*dh1dh2dh3_basic[0][2];
-	*(p+1)=plist[0]*dh1dh2dh3_basic[1][0]+plist[1]*dh1dh2dh3_basic[1][1]+plist[2]*dh1dh2dh3_basic[1][2];
+	double dh1dh3[NDOF2][numgrids]; //nodal derivative functions in actual coordinate system.
+
+	/*Get dh1dh2dh3 in actual coordinate system: */
+	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list, gauss_l1l2l3);
+
+	*(p+0)=plist[0]*dh1dh3[0][0]+plist[1]*dh1dh3[0][1]+plist[2]*dh1dh3[0][2];
+	*(p+1)=plist[0]*dh1dh3[1][0]+plist[1]*dh1dh3[1][1]+plist[2]*dh1dh3[1][2];
 
 }
@@ -3611,5 +3611,5 @@
 	double       xyz_list[numgrids][3];
 	int          doflist1[numgrids];
-	double       dh1dh2dh3_basic[NDOF2][numgrids];
+	double       dh1dh3[NDOF2][numgrids];
 
 	/* grid data: */
@@ -3733,5 +3733,5 @@
 
 		/*Get nodal functions derivatives*/
-		GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],&xyz_list[0][0],gauss_l1l2l3);
+		GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3);
 
 		/*Get B derivative: dB/dx */
@@ -3745,5 +3745,5 @@
 
 			//Add regularization term
-			grade_g_gaussian[i]-=numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dB[0]+dh1dh2dh3_basic[1][i]*dB[1]);
+			grade_g_gaussian[i]-=numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dB[0]+dh1dh3[1][i]*dB[1]);
 
 			//min dampening
@@ -3787,5 +3787,5 @@
 	double       xyz_list[numgrids][3];
 	int          doflist1[numgrids];
-	double       dh1dh2dh3_basic[NDOF2][numgrids];
+	double       dh1dh3[NDOF2][numgrids];
 
 	/* grid data: */
@@ -3953,5 +3953,5 @@
 
 		/*Get nodal functions derivatives*/
-		GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],&xyz_list[0][0],gauss_l1l2l3);
+		GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3);
 
 		/*Get k derivative: dk/dx */
@@ -3965,5 +3965,5 @@
 
 			//noise dampening d/dki(1/2*(dk/dx)^2)
-			grade_g_gaussian[i]+=-numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dk[0]+dh1dh2dh3_basic[1][i]*dk[1]);
+			grade_g_gaussian[i]+=-numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
 			
 			//min dampening
@@ -4005,5 +4005,5 @@
 	double       xyz_list[numgrids][3];
 	int          doflist1[numgrids];
-	double       dh1dh2dh3_basic[NDOF2][numgrids];
+	double       dh1dh3[NDOF2][numgrids];
 
 	/* grid data: */
@@ -4186,5 +4186,5 @@
 
 		/*Get nodal functions derivatives*/
-		GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],&xyz_list[0][0],gauss_l1l2l3);
+		GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3);
 
 		/*Get k derivative: dk/dx */
@@ -4201,5 +4201,5 @@
 
 			//Add regularization term
-			grade_g_gaussian[i]+= - numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dk[0]+dh1dh2dh3_basic[1][i]*dk[1]);
+			grade_g_gaussian[i]+= - numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
 
 			//min dampening
Index: /issm/trunk/src/c/objects/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Tria.h	(revision 2955)
+++ /issm/trunk/src/c/objects/Tria.h	(revision 2956)
@@ -89,6 +89,6 @@
 		void  GetBPrime_prog(double* Bprime_prog, double* xyz_list, double* gauss_l1l2l3);
 		void  GetNodalFunctions(double* l1l2l3, double* gauss_l1l2l3);
-		void  GetNodalFunctionsDerivativesBasic(double* dh1dh2dh3_basic,double* xyz_list, double* gauss_l1l2l3);
-		void  GetNodalFunctionsDerivativesParams(double* dl1dl2dl3,double* gauss_l1l2l3);
+		void  GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss_l1l2l3);
+		void  GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss_l1l2l3);
 		void  GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_l1l2l3);
 		void  GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3);