Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18530)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18531)
@@ -167,6 +167,6 @@
 
 	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;
 	}
@@ -174,5 +174,6 @@
 		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];*/
@@ -192,25 +193,27 @@
 
 		/* Start looping on the number of vertices: */
-		Gauss* gauss=NewGaussBase(1);
-		gauss->GaussPoint(0);
-
-		/*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];
-
-		/*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];
+		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 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);
+
+			/*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*/
@@ -1745,4 +1748,5 @@
 			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