Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 25721)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 25722)
@@ -366,6 +366,4 @@
 void       Penta::CalvingRateLevermann(){/*{{{*/
 
-	IssmDouble  xyz_list[NUMVERTICES][3];
-	GaussPenta* gauss=NULL;
 	IssmDouble  vx,vy,vel;
 	IssmDouble  strainparallel;
@@ -376,7 +374,4 @@
 	IssmDouble  calvingrate[NUMVERTICES];
 
-	/* Get node coordinates and dof list: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-
 	/*Retrieve all inputs and parameters we will need*/
 	Input* vx_input=this->GetInput(VxEnum);																		_assert_(vx_input);
@@ -387,5 +382,5 @@
 
 	/* Start looping on the number of vertices: */
-	gauss=new GaussPenta();
+	GaussPenta* gauss=new GaussPenta();
 	for (int iv=0;iv<NUMVERTICES;iv++){
 		gauss->GaussVertex(iv);
@@ -415,7 +410,5 @@
 	/*Clean up and return*/
 	delete gauss;
-
-}
-/*}}}*/
+}/*}}}*/
 void       Penta::CalvingFluxLevelset(){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 25721)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 25722)
@@ -406,5 +406,4 @@
 void       Tria::CalvingCrevasseDepth(){/*{{{*/
 
-	IssmDouble  xyz_list[NUMVERTICES][3];
 	IssmDouble  calvingrate[NUMVERTICES];
 	IssmDouble  vx,vy,vel;
@@ -414,7 +413,4 @@
 	IssmDouble  s_xx,s_xy,s_yy,s1,s2,stmp;
 	int crevasse_opening_stress;
-
-	/* Get node coordinates and dof list: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
 	/*retrieve the type of crevasse_opening_stress*/
@@ -504,5 +500,4 @@
 void       Tria::CalvingRateLevermann(){/*{{{*/
 
-	IssmDouble  xyz_list[NUMVERTICES][3];
 	IssmDouble  vx,vy,vel;
 	IssmDouble  strainparallel;
@@ -513,14 +508,11 @@
 	IssmDouble  calvingrate[NUMVERTICES];
 
-	/* Get node coordinates and dof list: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-
 	/*Retrieve all inputs and parameters we will need*/
-	Input* vx_input=this->GetInput(VxEnum);													_assert_(vx_input);
-	Input* vy_input=this->GetInput(VyEnum);													_assert_(vy_input);
-	Input* bs_input = this->GetInput(BaseEnum);                                 _assert_(bs_input);
-	Input* strainparallel_input=this->GetInput(StrainRateparallelEnum);			_assert_(strainparallel_input);
-	Input* strainperpendicular_input=this->GetInput(StrainRateperpendicularEnum);_assert_(strainperpendicular_input);
-	Input* levermanncoeff_input=this->GetInput(CalvinglevermannCoeffEnum);      _assert_(levermanncoeff_input);
+	Input *vx_input                  = this->GetInput(VxEnum);                      _assert_(vx_input);
+	Input *vy_input                  = this->GetInput(VyEnum);                      _assert_(vy_input);
+	Input *bs_input                  = this->GetInput(BaseEnum);                    _assert_(bs_input);
+	Input *strainparallel_input      = this->GetInput(StrainRateparallelEnum);      _assert_(strainparallel_input);
+	Input *strainperpendicular_input = this->GetInput(StrainRateperpendicularEnum); _assert_(strainperpendicular_input);
+	Input *levermanncoeff_input      = this->GetInput(CalvinglevermannCoeffEnum);   _assert_(levermanncoeff_input);
 
 	/* Start looping on the number of vertices: */
@@ -556,7 +548,5 @@
 	/*Clean up and return*/
 	delete gauss;
-
-}
-/*}}}*/
+}/*}}}*/
 void       Tria::CalvingFluxLevelset(){/*{{{*/
 
@@ -573,6 +563,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*/
@@ -599,10 +589,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
@@ -615,10 +605,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
@@ -631,10 +621,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{
@@ -673,5 +663,5 @@
 
 		/*Start looping on Gaussian points*/
-		Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
+		Gauss* gauss=this->NewGauss(&xyz_list[0][0],&xyz_front[0][0],3);
 		while(gauss->next()){
 			thickness_input->GetInputValue(&thickness,gauss);
@@ -707,7 +697,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*/
@@ -734,10 +723,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
@@ -750,10 +739,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
@@ -766,10 +755,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{
@@ -814,5 +803,5 @@
 
 		/*Start looping on Gaussian points*/
-		Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
+		Gauss* gauss=this->NewGauss(&xyz_list[0][0],&xyz_front[0][0],3);
 		while(gauss->next()){
 			thickness_input->GetInputValue(&thickness,gauss);
@@ -1372,5 +1361,4 @@
 	int         domaintype;
 	IssmDouble  phi,scalefactor,floatingarea;
-	IssmDouble *xyz_list  = NULL;
 
 	if(!IsIceInElement())return 0.;
@@ -1380,6 +1368,7 @@
 	if(domaintype!=Domain2DhorizontalEnum && domaintype!=Domain3DEnum) _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
 
-	this->GetVerticesCoordinates(&xyz_list);
-	phi=this->GetGroundedPortion(xyz_list);
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	phi=this->GetGroundedPortion(&xyz_list[0][0]);
 	floatingarea=(1-phi)*this->GetArea();
 	if(scaled==true){
@@ -1390,5 +1379,4 @@
 
 	/*Clean up and return*/
-	xDelete<IssmDouble>(xyz_list);
 	return floatingarea;
 }
@@ -1405,5 +1393,4 @@
 	}
 	/*Intermediaries*/
-	IssmDouble* xyz_list = NULL;
 	IssmDouble  bed_normal[2],base[NUMVERTICES],bed[NUMVERTICES],surface[NUMVERTICES],phi[NUMVERTICES];
 	IssmDouble  water_pressure[NUMVERTICES],pressureice[NUMVERTICES],pressure[NUMVERTICES];
@@ -1419,6 +1406,7 @@
 	IssmDouble gravity   = FindParam(ConstantsGEnum);
 
-	/* 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: */
 	Input* vx_input       = this->GetInput(VxEnum);       _assert_(vx_input);
@@ -1431,6 +1419,6 @@
 
 		/*Compute strain rate viscosity and pressure: */
-		this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
-		this->material->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL);
+		this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+		this->material->ViscosityFS(&viscosity,2,&xyz_list[0][0],gauss,vx_input,vy_input,NULL);
 		/*FIXME: this is for Hongju only*/
 	//	pressureice[iv]=gravity*rho_ice*(surface[iv]-base[iv]);
@@ -1447,5 +1435,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]+2*sigmaxy[i]*bed_normal[0]*bed_normal[1]);
 			water_pressure[i]=-gravity*rho_water*base[i];
@@ -1465,5 +1453,4 @@
 	/*clean up*/
 	delete gauss;
-	xDelete<IssmDouble>(xyz_list);
 }
 /*}}}*/
@@ -2479,5 +2466,4 @@
 	int         domaintype;
 	IssmDouble  phi,scalefactor,groundedarea;
-	IssmDouble *xyz_list  = NULL;
 
 	if(!IsIceInElement())return 0.;
@@ -2487,6 +2473,7 @@
 	if(domaintype!=Domain2DhorizontalEnum && domaintype!=Domain3DEnum) _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
 
-	this->GetVerticesCoordinates(&xyz_list);
-	phi=this->GetGroundedPortion(xyz_list);
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	phi=this->GetGroundedPortion(&xyz_list[0][0]);
 	groundedarea=phi*this->GetArea();
 	if(scaled==true){
@@ -2497,5 +2484,4 @@
 
 	/*Clean up and return*/
-	xDelete<IssmDouble>(xyz_list);
 	return groundedarea;
 }
@@ -2556,8 +2542,8 @@
 
 	/*Get ice front coordinates*/
-	IssmDouble *xyz_list = NULL;
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 	IssmDouble* xyz_front = NULL;
-	this->GetVerticesCoordinates(&xyz_list);
-	this->GetIcefrontCoordinates(&xyz_front,xyz_list,MaskIceLevelsetEnum);
+	this->GetIcefrontCoordinates(&xyz_front,&xyz_list[0][0],MaskIceLevelsetEnum);
 
 	/*Get normal vector*/
@@ -2584,5 +2570,5 @@
 
 	/*Start looping on Gaussian points*/
-	Gauss* gauss=this->NewGauss(xyz_list,xyz_front,3);
+	Gauss* gauss=this->NewGauss(&xyz_list[0][0],xyz_front,3);
 	while(gauss->next()){
 		thickness_input->GetInputValue(&thickness,gauss);
@@ -2595,5 +2581,4 @@
 
 	/*Cleanup and return*/
-	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(xyz_front);
 	delete gauss;
@@ -2615,7 +2600,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*/
@@ -2642,10 +2626,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
@@ -2658,10 +2642,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
@@ -2674,10 +2658,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{
@@ -2715,5 +2699,5 @@
 
 	/*Start looping on Gaussian points*/
-	Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
+	Gauss* gauss=this->NewGauss(&xyz_list[0][0],&xyz_front[0][0],3);
 	while(gauss->next()){
 		thickness_input->GetInputValue(&thickness,gauss);
@@ -2727,5 +2711,4 @@
 	/*Cleanup and return*/
 	delete gauss;
-	xDelete<IssmDouble>(xyz_list);
 	return flux;
 }
@@ -2746,6 +2729,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*/
@@ -2772,10 +2755,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
@@ -2788,10 +2771,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
@@ -2804,10 +2787,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{
@@ -2843,5 +2826,5 @@
 
 	/*Start looping on Gaussian points*/
-	Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
+	Gauss* gauss=this->NewGauss(&xyz_list[0][0],&xyz_front[0][0],3);
 	while(gauss->next()){
 		thickness_input->GetInputValue(&thickness,gauss);
@@ -2855,5 +2838,4 @@
 	/*Cleanup and return*/
 	delete gauss;
-	xDelete<IssmDouble>(xyz_list);
 	return flux;
 }
@@ -3282,5 +3264,4 @@
 	IssmDouble  volume;
 	IssmDouble  rho_ice;
-	IssmDouble* xyz_list=NULL;
 	int         point1;
 	IssmDouble  fraction1,fraction2;
@@ -3291,5 +3272,6 @@
 
 	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list);
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
 	/*Retrieve inputs required:*/
@@ -3313,5 +3295,5 @@
 	while(gauss->next()){
 
-		this->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		this->JacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
 		thickness_input->GetInputValue(&thickness, gauss);
 
@@ -3320,5 +3302,4 @@
 
 	/* clean up and Return: */
-	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(values);
 	delete gauss;
@@ -4046,5 +4027,4 @@
 void       Tria::StrainRateparallel(){/*{{{*/
 
-	IssmDouble *xyz_list = NULL;
 	IssmDouble  epsilon[3];
 	GaussTria* gauss=NULL;
@@ -4056,9 +4036,10 @@
 
 	/* 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*/
-	Input* vx_input=this->GetInput(VxEnum);                                  _assert_(vx_input);
-	Input* vy_input=this->GetInput(VyEnum);                                  _assert_(vy_input);
+	Input *vx_input = this->GetInput(VxEnum); _assert_(vx_input);
+	Input *vy_input = this->GetInput(VyEnum); _assert_(vy_input);
 
 	/* Start looping on the number of vertices: */
@@ -4073,5 +4054,5 @@
 
 		/*Compute strain rate viscosity and pressure: */
-		this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
+		this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
 		strainxx=epsilon[0];
 		strainyy=epsilon[1];
@@ -4087,10 +4068,8 @@
 	/*Clean up and return*/
 	delete gauss;
-	xDelete<IssmDouble>(xyz_list);
 }
 /*}}}*/
 void       Tria::StrainRateperpendicular(){/*{{{*/
 
-	IssmDouble *xyz_list = NULL;
 	GaussTria* gauss=NULL;
 	IssmDouble  epsilon[3];
@@ -4102,9 +4081,10 @@
 
 	/* 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*/
-	Input* vx_input=this->GetInput(VxEnum);                                  _assert_(vx_input);
-	Input* vy_input=this->GetInput(VyEnum);                                  _assert_(vy_input);
+	Input *vx_input = this->GetInput(VxEnum); _assert_(vx_input);
+	Input *vy_input = this->GetInput(VyEnum); _assert_(vy_input);
 
 	/* Start looping on the number of vertices: */
@@ -4119,5 +4099,5 @@
 
 		/*Compute strain rate viscosity and pressure: */
-		this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
+		this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
 		strainxx=epsilon[0];
 		strainyy=epsilon[1];
@@ -4133,5 +4113,4 @@
 	/*Clean up and return*/
 	delete gauss;
-	xDelete<IssmDouble>(xyz_list);
 }
 /*}}}*/
@@ -4219,7 +4198,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*/
@@ -4246,10 +4224,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
@@ -4262,10 +4240,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
@@ -4278,10 +4256,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{
@@ -4319,5 +4297,5 @@
 
 	/*Start looping on Gaussian points*/
-	Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
+	Gauss* gauss=this->NewGauss(&xyz_list[0][0],&xyz_front[0][0],3);
 	while(gauss->next()){
 		thickness_input->GetInputValue(&thickness,gauss);
@@ -4331,5 +4309,4 @@
 	/*Clean up and return*/
 	delete gauss;
-	xDelete<IssmDouble>(xyz_list);
 	return flux;
 }
@@ -4349,7 +4326,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*/
@@ -4376,10 +4352,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
@@ -4392,10 +4368,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
@@ -4408,10 +4384,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{
@@ -4455,5 +4431,5 @@
 
 	/*Start looping on Gaussian points*/
-	Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
+	Gauss* gauss=this->NewGauss(&xyz_list[0][0],&xyz_front[0][0],3);
 	while(gauss->next()){
 		thickness_input->GetInputValue(&thickness,gauss);