Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18530)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18531)
@@ -166,14 +166,15 @@
 void       Tria::ComputeSigmaNN(){/*{{{*/
 
 	if(!IsOnBase()){
-		IssmDouble sigma_nn=0;
-		this->inputs->AddInput(new TriaInput(SigmaNNEnum,&sigma_nn,P0Enum));
+		IssmDouble sigma_nn[3]={0.};
+		this->inputs->AddInput(new TriaInput(SigmaNNEnum,&sigma_nn[0],P1Enum));
 		return;
 	}
 	else{
 		IssmDouble* xyz_list=NULL;
 		IssmDouble *xyz_list_base=NULL;
-		IssmDouble  pressure,viscosity,sigma_nn;
+		IssmDouble  pressure,viscosity;
+		IssmDouble  sigma_nn[3];
 		IssmDouble  sigma_xx,sigma_xy,sigma_yy;
 		IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
 		IssmDouble  base_normal[2]; 
@@ -191,27 +192,29 @@
 		Input* vy_input=inputs->GetInput(VyEnum);             _assert_(vy_input);
 
 		/* Start looping on the number of vertices: */
-		Gauss* gauss=NewGaussBase(1);
-		gauss->GaussPoint(0);
+		Gauss* gauss = this->NewGauss();
+		for(int i=0;i<NUMVERTICES;i++){
+			gauss->GaussNode(P1Enum,i);
 
-		/*Compute strain rate viscosity and pressure: */
-		this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
-		this->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,NULL);
-		pressure_input->GetInputValue(&pressure,gauss);
+			/*Compute strain rate viscosity and pressure: */
+			this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
+			this->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,NULL);
+			pressure_input->GetInputValue(&pressure,gauss);
 
-		/*Compute Stress*/
-		sigma_xx=2*viscosity*epsilon[0]-pressure; // sigma = nu eps - pressure
-		sigma_yy=2*viscosity*epsilon[1]-pressure;
-		sigma_xy=2*viscosity*epsilon[2];
+			/*Compute Stress*/
+			sigma_xx=2*viscosity*epsilon[0]-pressure; // sigma = nu eps - pressure
+			sigma_yy=2*viscosity*epsilon[1]-pressure;
+			sigma_xy=2*viscosity*epsilon[2];
 
-		/*Get normal vector to the bed */
-		NormalBase(&base_normal[0],xyz_list_base);
+			/*Get normal vector to the bed */
+			NormalBase(&base_normal[0],xyz_list_base);
 
-		/*Compute sigma_nn*/
-		sigma_nn=sigma_xx*base_normal[0]*base_normal[0] + 2*sigma_xy*base_normal[0]*base_normal[1] + sigma_yy*base_normal[1]*base_normal[1];
+			/*Compute sigma_nn*/
+			sigma_nn[i]=sigma_xx*base_normal[0]*base_normal[0] + 2*sigma_xy*base_normal[0]*base_normal[1] + sigma_yy*base_normal[1]*base_normal[1];
+		}
 
 		/*Add Stress tensor components into inputs*/
-		this->inputs->AddInput(new TriaInput(SigmaNNEnum,&sigma_nn,P0Enum));
+		this->inputs->AddInput(new TriaInput(SigmaNNEnum,&sigma_nn[0],P1Enum));
 
 		/*Clean up and return*/
 		xDelete<IssmDouble>(xyz_list);
@@ -1744,6 +1747,7 @@
 			xz_plane[0]=cos(theta);       xz_plane[3]=0.;  
 			xz_plane[1]=sin(theta);       xz_plane[4]=0.;  
 			xz_plane[2]=0.;               xz_plane[5]=1.;          
+
 			if(groundedice>=0){
 				this->nodes[i]->DofInSSet(1); //vy
 			}
