Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 25722)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 25723)
@@ -425,6 +425,6 @@
 		IssmDouble        xyz_front[2][3];
 
-		IssmDouble *xyz_list = NULL;
-		this->GetVerticesCoordinates(&xyz_list);
+		IssmDouble  xyz_list[NUMVERTICES][3];
+      ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
 		/*Recover parameters and values*/
@@ -446,10 +446,10 @@
 				pt1 = 1; pt2 = 0;
 			}
-			xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
-			xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
-			xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
-			xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
-			xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
-			xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
+			xyz_front[pt2][0]=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]);
+			xyz_front[pt2][1]=xyz_list[2][1]+s1*(xyz_list[1][1]-xyz_list[2][1]);
+			xyz_front[pt2][2]=xyz_list[2][2]+s1*(xyz_list[1][2]-xyz_list[2][2]);
+			xyz_front[pt1][0]=xyz_list[2][0]+s2*(xyz_list[0][0]-xyz_list[2][0]);
+			xyz_front[pt1][1]=xyz_list[2][1]+s2*(xyz_list[0][1]-xyz_list[2][1]);
+			xyz_front[pt1][2]=xyz_list[2][2]+s2*(xyz_list[0][2]-xyz_list[2][2]);
 		}
 		else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
@@ -462,10 +462,10 @@
 			}
 
-			xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
-			xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
-			xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
-			xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
-			xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
-			xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
+			xyz_front[pt1][0]=xyz_list[0][0]+s1*(xyz_list[1][0]-xyz_list[0][0]);
+			xyz_front[pt1][1]=xyz_list[0][1]+s1*(xyz_list[1][1]-xyz_list[0][1]);
+			xyz_front[pt1][2]=xyz_list[0][2]+s1*(xyz_list[1][2]-xyz_list[0][2]);
+			xyz_front[pt2][0]=xyz_list[0][0]+s2*(xyz_list[2][0]-xyz_list[0][0]);
+			xyz_front[pt2][1]=xyz_list[0][1]+s2*(xyz_list[2][1]-xyz_list[0][1]);
+			xyz_front[pt2][2]=xyz_list[0][2]+s2*(xyz_list[2][2]-xyz_list[0][2]);
 		}
 		else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
@@ -478,10 +478,10 @@
 			}
 
-			xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
-			xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
-			xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
-			xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
-			xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
-			xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
+			xyz_front[pt2][0]=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]);
+			xyz_front[pt2][1]=xyz_list[1][1]+s1*(xyz_list[0][1]-xyz_list[1][1]);
+			xyz_front[pt2][2]=xyz_list[1][2]+s1*(xyz_list[0][2]-xyz_list[1][2]);
+			xyz_front[pt1][0]=xyz_list[1][0]+s2*(xyz_list[2][0]-xyz_list[1][0]);
+			xyz_front[pt1][1]=xyz_list[1][1]+s2*(xyz_list[2][1]-xyz_list[1][1]);
+			xyz_front[pt1][2]=xyz_list[1][2]+s2*(xyz_list[2][2]-xyz_list[1][2]);
 		}
 		else{
@@ -511,5 +511,5 @@
 
 		/*Start looping on Gaussian points*/
-		Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
+		Gauss* gauss=this->NewGaussBase(&xyz_list[0][0],&xyz_front[0][0],3);
 		while(gauss->next()){
 			thickness_input->GetInputValue(&thickness,gauss);
@@ -545,6 +545,6 @@
 		IssmDouble        xyz_front[2][3];
 
-		IssmDouble *xyz_list = NULL;
-		this->GetVerticesCoordinates(&xyz_list);
+		IssmDouble  xyz_list[NUMVERTICES][3];
+      ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
 		/*Recover parameters and values*/
@@ -566,10 +566,10 @@
 				pt1 = 1; pt2 = 0;
 			}
-			xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
-			xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
-			xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
-			xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
-			xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
-			xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
+			xyz_front[pt2][0]=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]);
+			xyz_front[pt2][1]=xyz_list[2][1]+s1*(xyz_list[1][1]-xyz_list[2][1]);
+			xyz_front[pt2][2]=xyz_list[2][2]+s1*(xyz_list[1][2]-xyz_list[2][2]);
+			xyz_front[pt1][0]=xyz_list[2][0]+s2*(xyz_list[0][0]-xyz_list[2][0]);
+			xyz_front[pt1][1]=xyz_list[2][1]+s2*(xyz_list[0][1]-xyz_list[2][1]);
+			xyz_front[pt1][2]=xyz_list[2][2]+s2*(xyz_list[0][2]-xyz_list[2][2]);
 		}
 		else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
@@ -582,10 +582,10 @@
 			}
 
-			xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
-			xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
-			xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
-			xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
-			xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
-			xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
+			xyz_front[pt1][0]=xyz_list[0][0]+s1*(xyz_list[1][0]-xyz_list[0][0]);
+			xyz_front[pt1][1]=xyz_list[0][1]+s1*(xyz_list[1][1]-xyz_list[0][1]);
+			xyz_front[pt1][2]=xyz_list[0][2]+s1*(xyz_list[1][2]-xyz_list[0][2]);
+			xyz_front[pt2][0]=xyz_list[0][0]+s2*(xyz_list[2][0]-xyz_list[0][0]);
+			xyz_front[pt2][1]=xyz_list[0][1]+s2*(xyz_list[2][1]-xyz_list[0][1]);
+			xyz_front[pt2][2]=xyz_list[0][2]+s2*(xyz_list[2][2]-xyz_list[0][2]);
 		}
 		else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
@@ -598,10 +598,10 @@
 			}
 
-			xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
-			xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
-			xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
-			xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
-			xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
-			xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
+			xyz_front[pt2][0]=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]);
+			xyz_front[pt2][1]=xyz_list[1][1]+s1*(xyz_list[0][1]-xyz_list[1][1]);
+			xyz_front[pt2][2]=xyz_list[1][2]+s1*(xyz_list[0][2]-xyz_list[1][2]);
+			xyz_front[pt1][0]=xyz_list[1][0]+s2*(xyz_list[2][0]-xyz_list[1][0]);
+			xyz_front[pt1][1]=xyz_list[1][1]+s2*(xyz_list[2][1]-xyz_list[1][1]);
+			xyz_front[pt1][2]=xyz_list[1][2]+s2*(xyz_list[2][2]-xyz_list[1][2]);
 		}
 		else{
@@ -637,5 +637,5 @@
 
 		/*Start looping on Gaussian points*/
-		Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
+		Gauss* gauss=this->NewGaussBase(&xyz_list[0][0],&xyz_front[0][0],3);
 		while(gauss->next()){
 			thickness_input->GetInputValue(&thickness,gauss);
@@ -1115,5 +1115,4 @@
 
 	/*Intermediaries*/
-	IssmDouble* xyz_list = NULL;
 	IssmDouble  bed_normal[3],base[NUMVERTICES],bed[NUMVERTICES],surface[NUMVERTICES],phi[NUMVERTICES];
 	IssmDouble  water_pressure[NUMVERTICES],pressureice[NUMVERTICES],pressure[NUMVERTICES];
@@ -1131,5 +1130,6 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list);
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
 	/*Retrieve all inputs we will be needing: */
@@ -1144,6 +1144,6 @@
 
 		/*Compute strain rate viscosity and pressure: */
-		this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
-		this->material->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
+		this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
+		this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
 		/*FIXME: this is for Hongju only*/
 		//pressureice[iv]=gravity*rho_ice*(surface[iv]-base[iv]);
@@ -1164,5 +1164,5 @@
 		/*If was grounded*/
 		if (phi[i]>=0.){
-			NormalBase(&bed_normal[0],xyz_list);
+			NormalBase(&bed_normal[0],&xyz_list[0][0]);
 			sigma_nn[i]=-1*(sigmaxx[i]*bed_normal[0]*bed_normal[0] + sigmayy[i]*bed_normal[1]*bed_normal[1] + sigmazz[i]*bed_normal[2]*bed_normal[2]+2*sigmaxy[i]*bed_normal[0]*bed_normal[1]+2*sigmaxz[i]*bed_normal[0]*bed_normal[2]+2*sigmayz[i]*bed_normal[1]*bed_normal[2]);
 			water_pressure[i]=-gravity*rho_water*base[i];
@@ -1183,5 +1183,4 @@
 	/*clean up*/
 	delete gauss;
-	xDelete<IssmDouble>(xyz_list);
 }
 /*}}}*/
@@ -1896,6 +1895,6 @@
 	IssmDouble        xyz_front[2][3];
 
-	IssmDouble *xyz_list = NULL;
-	this->GetVerticesCoordinates(&xyz_list);
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
 	/*Recover parameters and values*/
@@ -1917,10 +1916,10 @@
 			pt1 = 1; pt2 = 0;
 		}
-		xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
-		xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
-		xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
-		xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
-		xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
-		xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
+		xyz_front[pt2][0]=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]);
+		xyz_front[pt2][1]=xyz_list[2][1]+s1*(xyz_list[1][1]-xyz_list[2][1]);
+		xyz_front[pt2][2]=xyz_list[2][2]+s1*(xyz_list[1][2]-xyz_list[2][2]);
+		xyz_front[pt1][0]=xyz_list[2][0]+s2*(xyz_list[0][0]-xyz_list[2][0]);
+		xyz_front[pt1][1]=xyz_list[2][1]+s2*(xyz_list[0][1]-xyz_list[2][1]);
+		xyz_front[pt1][2]=xyz_list[2][2]+s2*(xyz_list[0][2]-xyz_list[2][2]);
 	}
 	else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
@@ -1933,10 +1932,10 @@
 		}
 
-		xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
-		xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
-		xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
-		xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
-		xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
-		xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
+		xyz_front[pt1][0]=xyz_list[0][0]+s1*(xyz_list[1][0]-xyz_list[0][0]);
+		xyz_front[pt1][1]=xyz_list[0][1]+s1*(xyz_list[1][1]-xyz_list[0][1]);
+		xyz_front[pt1][2]=xyz_list[0][2]+s1*(xyz_list[1][2]-xyz_list[0][2]);
+		xyz_front[pt2][0]=xyz_list[0][0]+s2*(xyz_list[2][0]-xyz_list[0][0]);
+		xyz_front[pt2][1]=xyz_list[0][1]+s2*(xyz_list[2][1]-xyz_list[0][1]);
+		xyz_front[pt2][2]=xyz_list[0][2]+s2*(xyz_list[2][2]-xyz_list[0][2]);
 	}
 	else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
@@ -1949,10 +1948,10 @@
 		}
 
-		xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
-		xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
-		xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
-		xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
-		xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
-		xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
+		xyz_front[pt2][0]=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]);
+		xyz_front[pt2][1]=xyz_list[1][1]+s1*(xyz_list[0][1]-xyz_list[1][1]);
+		xyz_front[pt2][2]=xyz_list[1][2]+s1*(xyz_list[0][2]-xyz_list[1][2]);
+		xyz_front[pt1][0]=xyz_list[1][0]+s2*(xyz_list[2][0]-xyz_list[1][0]);
+		xyz_front[pt1][1]=xyz_list[1][1]+s2*(xyz_list[2][1]-xyz_list[1][1]);
+		xyz_front[pt1][2]=xyz_list[1][2]+s2*(xyz_list[2][2]-xyz_list[1][2]);
 	}
 	else{
@@ -1983,5 +1982,5 @@
 
 	/*Start looping on Gaussian points*/
-	Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
+	Gauss* gauss=this->NewGaussBase(&xyz_list[0][0],&xyz_front[0][0],3);
 	while(gauss->next()){
 		thickness_input->GetInputValue(&thickness,gauss);
@@ -1995,7 +1994,5 @@
 	/*Clean up and return*/
 	delete gauss;
-	xDelete<IssmDouble>(xyz_list);
 	return flux;
-
 }
 /*}}}*/
@@ -2014,6 +2011,6 @@
 	IssmDouble        xyz_front[2][3];
 
-	IssmDouble *xyz_list = NULL;
-	this->GetVerticesCoordinates(&xyz_list);
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
 	/*Recover parameters and values*/
@@ -2035,10 +2032,10 @@
 			pt1 = 1; pt2 = 0;
 		}
-		xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
-		xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
-		xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
-		xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
-		xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
-		xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
+		xyz_front[pt2][0]=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]);
+		xyz_front[pt2][1]=xyz_list[2][1]+s1*(xyz_list[1][1]-xyz_list[2][1]);
+		xyz_front[pt2][2]=xyz_list[2][2]+s1*(xyz_list[1][2]-xyz_list[2][2]);
+		xyz_front[pt1][0]=xyz_list[2][0]+s2*(xyz_list[0][0]-xyz_list[2][0]);
+		xyz_front[pt1][1]=xyz_list[2][1]+s2*(xyz_list[0][1]-xyz_list[2][1]);
+		xyz_front[pt1][2]=xyz_list[2][2]+s2*(xyz_list[0][2]-xyz_list[2][2]);
 	}
 	else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
@@ -2051,10 +2048,10 @@
 		}
 
-		xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
-		xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
-		xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
-		xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
-		xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
-		xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
+		xyz_front[pt1][0]=xyz_list[0][0]+s1*(xyz_list[1][0]-xyz_list[0][0]);
+		xyz_front[pt1][1]=xyz_list[0][1]+s1*(xyz_list[1][1]-xyz_list[0][1]);
+		xyz_front[pt1][2]=xyz_list[0][2]+s1*(xyz_list[1][2]-xyz_list[0][2]);
+		xyz_front[pt2][0]=xyz_list[0][0]+s2*(xyz_list[2][0]-xyz_list[0][0]);
+		xyz_front[pt2][1]=xyz_list[0][1]+s2*(xyz_list[2][1]-xyz_list[0][1]);
+		xyz_front[pt2][2]=xyz_list[0][2]+s2*(xyz_list[2][2]-xyz_list[0][2]);
 	}
 	else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
@@ -2067,10 +2064,10 @@
 		}
 
-		xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
-		xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
-		xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
-		xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
-		xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
-		xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
+		xyz_front[pt2][0]=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]);
+		xyz_front[pt2][1]=xyz_list[1][1]+s1*(xyz_list[0][1]-xyz_list[1][1]);
+		xyz_front[pt2][2]=xyz_list[1][2]+s1*(xyz_list[0][2]-xyz_list[1][2]);
+		xyz_front[pt1][0]=xyz_list[1][0]+s2*(xyz_list[2][0]-xyz_list[1][0]);
+		xyz_front[pt1][1]=xyz_list[1][1]+s2*(xyz_list[2][1]-xyz_list[1][1]);
+		xyz_front[pt1][2]=xyz_list[1][2]+s2*(xyz_list[2][2]-xyz_list[1][2]);
 	}
 	else{
@@ -2103,5 +2100,5 @@
 
 	/*Start looping on Gaussian points*/
-	Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
+	Gauss* gauss=this->NewGaussBase(&xyz_list[0][0],&xyz_front[0][0],3);
 	while(gauss->next()){
 		thickness_input->GetInputValue(&thickness,gauss);
@@ -2114,8 +2111,6 @@
 
 	/*Clean up and return*/
-	xDelete<IssmDouble>(xyz_list);
 	delete gauss;
 	return flux;
-
 }
 /*}}}*/
@@ -3418,5 +3413,4 @@
 void       Penta::StrainRateparallel(){/*{{{*/
 
-	IssmDouble *xyz_list = NULL;
 	IssmDouble  epsilon[6];
 	GaussPenta* gauss=NULL;
@@ -3428,5 +3422,6 @@
 
 	/* Get node coordinates and dof list: */
-	this->GetVerticesCoordinates(&xyz_list);
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
 	/*Retrieve all inputs we will need*/
@@ -3446,5 +3441,5 @@
 
 		/*Compute strain rate and viscosity: */
-		this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
+		this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
 		strainxx=epsilon[0];
 		strainyy=epsilon[1];
@@ -3460,10 +3455,8 @@
 	/*Clean up and return*/
 	delete gauss;
-	xDelete<IssmDouble>(xyz_list);
 }
 /*}}}*/
 void       Penta::StrainRateperpendicular(){/*{{{*/
 
-	IssmDouble *xyz_list = NULL;
 	IssmDouble  epsilon[6];
 	GaussPenta* gauss=NULL;
@@ -3475,5 +3468,6 @@
 
 	/* Get node coordinates and dof list: */
-	this->GetVerticesCoordinates(&xyz_list);
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
 	/*Retrieve all inputs we will need*/
@@ -3493,5 +3487,5 @@
 
 		/*Compute strain rate and viscosity: */
-		this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
+		this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
 		strainxx=epsilon[0];
 		strainyy=epsilon[1];
@@ -3507,5 +3501,4 @@
 	/*Clean up and return*/
 	delete gauss;
-	xDelete<IssmDouble>(xyz_list);
 }
 /*}}}*/
@@ -3720,6 +3713,6 @@
 	IssmDouble        xyz_front[2][3];
 
-	IssmDouble *xyz_list = NULL;
-	this->GetVerticesCoordinates(&xyz_list);
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
 	/*Recover parameters and values*/
@@ -3741,10 +3734,10 @@
 			pt1 = 1; pt2 = 0;
 		}
-		xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
-		xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
-		xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
-		xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
-		xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
-		xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
+		xyz_front[pt2][0]=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]);
+		xyz_front[pt2][1]=xyz_list[2][1]+s1*(xyz_list[1][1]-xyz_list[2][1]);
+		xyz_front[pt2][2]=xyz_list[2][2]+s1*(xyz_list[1][2]-xyz_list[2][2]);
+		xyz_front[pt1][0]=xyz_list[2][0]+s2*(xyz_list[0][0]-xyz_list[2][0]);
+		xyz_front[pt1][1]=xyz_list[2][1]+s2*(xyz_list[0][1]-xyz_list[2][1]);
+		xyz_front[pt1][2]=xyz_list[2][2]+s2*(xyz_list[0][2]-xyz_list[2][2]);
 	}
 	else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
@@ -3757,10 +3750,10 @@
 		}
 
-		xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
-		xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
-		xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
-		xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
-		xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
-		xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
+		xyz_front[pt1][0]=xyz_list[0][0]+s1*(xyz_list[1][0]-xyz_list[0][0]);
+		xyz_front[pt1][1]=xyz_list[0][1]+s1*(xyz_list[1][1]-xyz_list[0][1]);
+		xyz_front[pt1][2]=xyz_list[0][2]+s1*(xyz_list[1][2]-xyz_list[0][2]);
+		xyz_front[pt2][0]=xyz_list[0][0]+s2*(xyz_list[2][0]-xyz_list[0][0]);
+		xyz_front[pt2][1]=xyz_list[0][1]+s2*(xyz_list[2][1]-xyz_list[0][1]);
+		xyz_front[pt2][2]=xyz_list[0][2]+s2*(xyz_list[2][2]-xyz_list[0][2]);
 	}
 	else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
@@ -3773,10 +3766,10 @@
 		}
 
-		xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
-		xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
-		xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
-		xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
-		xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
-		xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
+		xyz_front[pt2][0]=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]);
+		xyz_front[pt2][1]=xyz_list[1][1]+s1*(xyz_list[0][1]-xyz_list[1][1]);
+		xyz_front[pt2][2]=xyz_list[1][2]+s1*(xyz_list[0][2]-xyz_list[1][2]);
+		xyz_front[pt1][0]=xyz_list[1][0]+s2*(xyz_list[2][0]-xyz_list[1][0]);
+		xyz_front[pt1][1]=xyz_list[1][1]+s2*(xyz_list[2][1]-xyz_list[1][1]);
+		xyz_front[pt1][2]=xyz_list[1][2]+s2*(xyz_list[2][2]-xyz_list[1][2]);
 	}
 	else{
@@ -3806,5 +3799,5 @@
 
 	/*Start looping on Gaussian points*/
-	Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
+	Gauss* gauss=this->NewGaussBase(&xyz_list[0][0],&xyz_front[0][0],3);
 	while(gauss->next()){
 		thickness_input->GetInputValue(&thickness,gauss);
@@ -3819,5 +3812,4 @@
 	delete gauss;
 	return flux;
-
 }
 /*}}}*/
@@ -3836,6 +3828,6 @@
 	IssmDouble        xyz_front[2][3];
 
-	IssmDouble *xyz_list = NULL;
-	this->GetVerticesCoordinates(&xyz_list);
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
 	/*Recover parameters and values*/
@@ -3857,10 +3849,10 @@
 			pt1 = 1; pt2 = 0;
 		}
-		xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
-		xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
-		xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
-		xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
-		xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
-		xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
+		xyz_front[pt2][0]=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]);
+		xyz_front[pt2][1]=xyz_list[2][1]+s1*(xyz_list[1][1]-xyz_list[2][1]);
+		xyz_front[pt2][2]=xyz_list[2][2]+s1*(xyz_list[1][2]-xyz_list[2][2]);
+		xyz_front[pt1][0]=xyz_list[2][0]+s2*(xyz_list[0][0]-xyz_list[2][0]);
+		xyz_front[pt1][1]=xyz_list[2][1]+s2*(xyz_list[0][1]-xyz_list[2][1]);
+		xyz_front[pt1][2]=xyz_list[2][2]+s2*(xyz_list[0][2]-xyz_list[2][2]);
 	}
 	else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
@@ -3873,10 +3865,10 @@
 		}
 
-		xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
-		xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
-		xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
-		xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
-		xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
-		xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
+		xyz_front[pt1][0]=xyz_list[0][0]+s1*(xyz_list[1][0]-xyz_list[0][0]);
+		xyz_front[pt1][1]=xyz_list[0][1]+s1*(xyz_list[1][1]-xyz_list[0][1]);
+		xyz_front[pt1][2]=xyz_list[0][2]+s1*(xyz_list[1][2]-xyz_list[0][2]);
+		xyz_front[pt2][0]=xyz_list[0][0]+s2*(xyz_list[2][0]-xyz_list[0][0]);
+		xyz_front[pt2][1]=xyz_list[0][1]+s2*(xyz_list[2][1]-xyz_list[0][1]);
+		xyz_front[pt2][2]=xyz_list[0][2]+s2*(xyz_list[2][2]-xyz_list[0][2]);
 	}
 	else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
@@ -3889,10 +3881,10 @@
 		}
 
-		xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
-		xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
-		xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
-		xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
-		xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
-		xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
+		xyz_front[pt2][0]=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]);
+		xyz_front[pt2][1]=xyz_list[1][1]+s1*(xyz_list[0][1]-xyz_list[1][1]);
+		xyz_front[pt2][2]=xyz_list[1][2]+s1*(xyz_list[0][2]-xyz_list[1][2]);
+		xyz_front[pt1][0]=xyz_list[1][0]+s2*(xyz_list[2][0]-xyz_list[1][0]);
+		xyz_front[pt1][1]=xyz_list[1][1]+s2*(xyz_list[2][1]-xyz_list[1][1]);
+		xyz_front[pt1][2]=xyz_list[1][2]+s2*(xyz_list[2][2]-xyz_list[1][2]);
 	}
 	else{
@@ -3931,5 +3923,5 @@
 
 	/*Start looping on Gaussian points*/
-	Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
+	Gauss* gauss=this->NewGaussBase(&xyz_list[0][0],&xyz_front[0][0],3);
 	while(gauss->next()){
 		thickness_input->GetInputValue(&thickness,gauss);
@@ -3950,5 +3942,4 @@
 	delete gauss;
 	return flux;
-
 }
 /*}}}*/