Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16432)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16433)
@@ -10160,14 +10160,13 @@
 void  Penta::InputUpdateFromSolutionStressbalanceSSAFS(IssmDouble* solution){
 
-	const int    numdofm=NDOF2*NUMVERTICES;
-	const int    numdofs=NDOF3*NUMVERTICES;
-	const int    numdof2d=NDOF2*NUMVERTICES2D;
-	const int    numdofpressure=NDOF1*NUMVERTICES;
+	const int    numdofSSA = NDOF2*NUMVERTICES;
+	const int    numdof2d  = NDOF2*NUMVERTICES2D;
+	const int    numdofFSv = NDOF3*NUMVERTICES;
+	const int    numdofFSp = NDOF1*NUMVERTICES;
 
 	int     i;
 	IssmDouble  FSreconditioning;
-	IssmDouble  SSA_values[numdofm];
-	IssmDouble  FS_values[numdofs];
-	IssmDouble  Pressure_values[numdofs];
+	IssmDouble  SSA_values[numdofSSA];
+	IssmDouble  FS_values[numdofFSv+numdofFSp];
 	IssmDouble  vx[NUMVERTICES];
 	IssmDouble  vy[NUMVERTICES];
@@ -10178,8 +10177,13 @@
 	IssmDouble  pressure[NUMVERTICES];
 	IssmDouble  xyz_list[NUMVERTICES][3];
-	int*    doflistm        = NULL;
-	int*    doflists        = NULL;
-	int*    doflistpressure = NULL;
-	Penta   *penta          = NULL;
+	int   *doflistSSA = NULL;
+	int   *doflistFSv = NULL;
+	int   *doflistFSp = NULL;
+	Penta *penta      = NULL;
+
+	/*Prepare coordinate system list*/
+	int* cs_list = xNew<int>(2*NUMVERTICES);
+	for(i=0;i<NUMVERTICES;i++) cs_list[i]             = XYZEnum;
+	for(i=0;i<NUMVERTICES;i++) cs_list[NUMVERTICES+i] = PressureEnum;
 
 	/*OK, we have to add results of this element for SSA 
@@ -10188,7 +10192,7 @@
 
 	/*Get dof listof this element (SSA dofs) and of the penta at base (SSA dofs): */
-	penta->GetDofList(&doflistm,SSAApproximationEnum,GsetEnum);
-	GetDofList(&doflists,FSvelocityEnum,GsetEnum);
-	GetDofListPressure(&doflistpressure,GsetEnum);
+	penta->GetDofList(&doflistSSA,SSAApproximationEnum,GsetEnum);
+	GetDofList(&doflistFSv,FSvelocityEnum,GsetEnum);
+	GetDofListPressure(&doflistFSp,GsetEnum);
 	this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
 
@@ -10198,25 +10202,25 @@
 	/*Use the dof list to index into the solution vector: */
 	for(i=0;i<numdof2d;i++){
-		SSA_values[i]=solution[doflistm[i]];
-		SSA_values[i+numdof2d]=solution[doflistm[i]];
-	}
-	for(i=0;i<numdofs;i++)FS_values[i]=solution[doflists[i]];
-	for(i=0;i<numdofpressure;i++) Pressure_values[i]=solution[doflistpressure[i]];
+		SSA_values[i]=solution[doflistSSA[i]];
+		SSA_values[i+numdof2d]=solution[doflistSSA[i]];
+	}
+	for(i=0;i<numdofFSv;i++) FS_values[i]=solution[doflistFSv[i]];
+	for(i=0;i<numdofFSp;i++) FS_values[numdofFSv+i]=solution[doflistFSp[i]];
 
 	/*Transform solution in Cartesian Space*/
 	TransformSolutionCoord(&SSA_values[0],this->nodes,NUMVERTICES,XYEnum);
-	TransformSolutionCoord(&FS_values[0],this->nodes,NUMVERTICES,XYZEnum);
+	TransformSolutionCoord(&FS_values[0],this->nodes,2*NUMVERTICES,cs_list);
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
 	for(i=0;i<NUMVERTICES;i++){
-		vx[i]=FS_values[i*NDOF3+0]+SSA_values[i*NDOF2+0];
-		vy[i]=FS_values[i*NDOF3+1]+SSA_values[i*NDOF2+1];
-		vzFS[i]=FS_values[i*NDOF3+2];
-		pressure[i]=Pressure_values[i*NDOF1]*FSreconditioning;
+		vx[i]       = FS_values[i*NDOF3+0]+SSA_values[i*NDOF2+0];
+		vy[i]       = FS_values[i*NDOF3+1]+SSA_values[i*NDOF2+1];
+		vzFS[i]     = FS_values[i*NDOF3+2];
+		pressure[i] = FS_values[NUMVERTICES*NDOF3+i]*FSreconditioning;
 
 		/*Check solution*/
 		if(xIsNan<IssmDouble>(vx[i]))       _error_("NaN found in solution vector");
 		if(xIsNan<IssmDouble>(vy[i]))       _error_("NaN found in solution vector");
-		if(xIsNan<IssmDouble>(vzFS[i])) _error_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vzFS[i]))     _error_("NaN found in solution vector");
 		if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
 	}
@@ -10236,6 +10240,6 @@
 	/*Now Compute vel*/
 	for(i=0;i<NUMVERTICES;i++) {
-		vz[i]=vzSSA[i]+vzFS[i];
-		vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+		vz[i]  = vzSSA[i]+vzFS[i];
+		vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
 	}
 
@@ -10256,7 +10260,8 @@
 
 	/*Free ressources:*/
-	xDelete<int>(doflistm);
-	xDelete<int>(doflists);
-	xDelete<int>(doflistpressure);
+	xDelete<int>(doflistSSA);
+	xDelete<int>(doflistFSv);
+	xDelete<int>(doflistFSp);
+	xDelete<int>(cs_list);
 }
 /*}}}*/
@@ -10418,12 +10423,11 @@
 void  Penta::InputUpdateFromSolutionStressbalanceHOFS(IssmDouble* solution){
 
-	const int    numdofp=NDOF2*NUMVERTICES;
-	const int    numdofs=NDOF3*NUMVERTICES;
-	const int    numdofpressure=NDOF1*NUMVERTICES;
+	const int    numdofHO  = NDOF2*NUMVERTICES;
+	const int    numdofFSv = NDOF3*NUMVERTICES;
+	const int    numdofFSp = NDOF1*NUMVERTICES;
 
 	int        i;
-	IssmDouble HO_values[numdofp];
-	IssmDouble FS_values[numdofs];
-	IssmDouble Pressure_values[numdofpressure];
+	IssmDouble HO_values[numdofHO];
+	IssmDouble FS_values[numdofFSv+numdofFSp];
 	IssmDouble vx[NUMVERTICES];
 	IssmDouble vy[NUMVERTICES];
@@ -10435,12 +10439,17 @@
 	IssmDouble xyz_list[NUMVERTICES][3];
 	IssmDouble FSreconditioning;
-	int*       doflistp        = NULL;
-	int*       doflists        = NULL;
-	int*       doflistpressure = NULL;
+	int*       doflistHO        = NULL;
+	int*       doflistFSv        = NULL;
+	int*       doflistFSp = NULL;
+
+	/*Prepare coordinate system list*/
+	int* cs_list = xNew<int>(2*NUMVERTICES);
+	for(i=0;i<NUMVERTICES;i++) cs_list[i]             = XYZEnum;
+	for(i=0;i<NUMVERTICES;i++) cs_list[NUMVERTICES+i] = PressureEnum;
 
 	/*Get dof listof this element (HO dofs) and of the penta at base (SSA dofs): */
-	GetDofList(&doflistp,HOApproximationEnum,GsetEnum);
-	GetDofList(&doflists,FSvelocityEnum,GsetEnum);
-	GetDofListPressure(&doflistpressure,GsetEnum);
+	GetDofList(&doflistHO,HOApproximationEnum,GsetEnum);
+	GetDofList(&doflistFSv,FSvelocityEnum,GsetEnum);
+	GetDofListPressure(&doflistFSp,GsetEnum);
 	this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
 
@@ -10449,23 +10458,23 @@
 
 	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<numdofp;i++) HO_values[i]=solution[doflistp[i]];
-	for(i=0;i<numdofs;i++) FS_values[i]=solution[doflists[i]];
-	for(i=0;i<numdofpressure;i++) Pressure_values[i]=solution[doflistpressure[i]];
+	for(i=0;i<numdofHO;i++)  HO_values[i]=solution[doflistHO[i]];
+	for(i=0;i<numdofFSv;i++) FS_values[i]=solution[doflistFSv[i]];
+	for(i=0;i<numdofFSp;i++) FS_values[numdofFSv+i]=solution[doflistFSp[i]];
 
 	/*Transform solution in Cartesian Space*/
 	TransformSolutionCoord(&HO_values[0],this->nodes,NUMVERTICES,XYEnum);
-	TransformSolutionCoord(&FS_values[0],this->nodes,NUMVERTICES,XYZEnum);
+	TransformSolutionCoord(&FS_values[0],this->nodes,2*NUMVERTICES,cs_list);
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
 	for(i=0;i<NUMVERTICES;i++){
-		vx[i]=FS_values[i*NDOF3+0]+HO_values[i*NDOF2+0];
-		vy[i]=FS_values[i*NDOF3+1]+HO_values[i*NDOF2+1];
-		vzFS[i]=FS_values[i*NDOF3+2];
-		pressure[i]=Pressure_values[i*NDOF1]*FSreconditioning;
+		vx[i]       = FS_values[i*NDOF3+0]+HO_values[i*NDOF2+0];
+		vy[i]       = FS_values[i*NDOF3+1]+HO_values[i*NDOF2+1];
+		vzFS[i]     = FS_values[i*NDOF3+2];
+		pressure[i] = FS_values[NUMVERTICES*NDOF3+i]*FSreconditioning;
 
 		/*Check solution*/
 		if(xIsNan<IssmDouble>(vx[i]))       _error_("NaN found in solution vector");
 		if(xIsNan<IssmDouble>(vy[i]))       _error_("NaN found in solution vector");
-		if(xIsNan<IssmDouble>(vzFS[i])) _error_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vzFS[i]))     _error_("NaN found in solution vector");
 		if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
 	}
@@ -10485,6 +10494,6 @@
 	/*Now Compute vel*/
 	for(i=0;i<NUMVERTICES;i++) {
-		vz[i]=vzHO[i]+vzFS[i];
-		vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+		vz[i]  = vzHO[i]+vzFS[i];
+		vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
 	}
 
@@ -10505,7 +10514,8 @@
 
 	/*Free ressources:*/
-	xDelete<int>(doflistp);
-	xDelete<int>(doflists);
-	xDelete<int>(doflistpressure);
+	xDelete<int>(doflistHO);
+	xDelete<int>(doflistFSv);
+	xDelete<int>(doflistFSp);
+	xDelete<int>(cs_list);
 }
 /*}}}*/
@@ -10702,5 +10712,5 @@
 	/*Prepare coordinate system list*/
 	int* cs_list = xNew<int>(vnumnodes+pnumnodes);
-	for(i=0;i<vnumnodes;i++) cs_list[i]           = XYEnum;
+	for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
 	for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16432)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16433)
@@ -2455,5 +2455,4 @@
 	GaussTria* gauss = new GaussTria();
 	for(int i=0;i<numindices;i++){//FIXME
-
 
 		gauss->GaussNode(this->VelocityInterpolation(),indices[i]);