Index: /issm/trunk/src/c/ComputeBasalStressx/ComputeBasalStressx.cpp
===================================================================
--- /issm/trunk/src/c/ComputeBasalStressx/ComputeBasalStressx.cpp	(revision 3529)
+++ /issm/trunk/src/c/ComputeBasalStressx/ComputeBasalStressx.cpp	(revision 3529)
@@ -0,0 +1,45 @@
+/*!\file ComputeBasalStressx
+ * \brief: compute pressure according to each element
+ */
+
+#include "./ComputeBasalStressx.h"
+
+#include "../shared/shared.h"
+#include "../include/macros.h"
+#include "../toolkits/toolkits.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+void	ComputeBasalStressx( Vec* psigma,DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials,DataSet* parameters,
+			ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
+
+	int i;
+
+	int  found=0;
+	double numberofelements;
+
+	/*output: */
+	Vec sigma=NULL;
+
+	/*Recover numberofelements: */
+	found= parameters->FindParam(&numberofelements,"numberofelements");
+	if (!found) ISSMERROR("numberofelements not provided in parameters");
+
+	/*Allocate sigma on numberofelements: */
+	sigma=NewVec((int)numberofelements);
+
+	/*Get elements configured: */
+	elements->Configure(elements,loads, nodes,vertices, materials,parameters);
+	nodes->Configure(elements,loads, nodes,vertices, materials,parameters);
+	parameters->Configure(elements,loads, nodes,vertices, materials,parameters);
+
+	/*Call on dataset driver: */
+	elements->ComputeBasalStress(sigma,inputs,analysis_type,sub_analysis_type);
+
+	/*Assemble vector: */
+	VecAssemblyBegin(sigma);
+	VecAssemblyEnd(sigma);
+
+	/*Assign output pointers: */
+	*psigma=sigma;
+	
+}
Index: /issm/trunk/src/c/ComputeBasalStressx/ComputeBasalStressx.h
===================================================================
--- /issm/trunk/src/c/ComputeBasalStressx/ComputeBasalStressx.h	(revision 3529)
+++ /issm/trunk/src/c/ComputeBasalStressx/ComputeBasalStressx.h	(revision 3529)
@@ -0,0 +1,16 @@
+/*!\file:  ComputeBasalStressx.h
+ * \brief header file pressure computation
+ */ 
+
+#ifndef _COMPUTEBASALSTRESSX_H
+#define _COMPUTEBASALSTRESSX_H
+
+#include "../objects/objects.h"
+#include "../DataSet/DataSet.h"
+
+/* local prototypes: */
+void	ComputeBasalStressx( Vec* pp_g,DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials,  DataSet* parameters,
+			ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
+
+#endif  /* _COMPUTEBASALSTRESSX_H */
+
Index: /issm/trunk/src/c/ComputePressurex/ComputePressurex.cpp
===================================================================
--- /issm/trunk/src/c/ComputePressurex/ComputePressurex.cpp	(revision 3528)
+++ /issm/trunk/src/c/ComputePressurex/ComputePressurex.cpp	(revision 3529)
@@ -10,5 +10,6 @@
 #include "../EnumDefinitions/EnumDefinitions.h"
 
-void	ComputePressurex( Vec* pp_g,DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials,DataSet* parameters){
+void	ComputePressurex( Vec* pp_g,DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials,DataSet* parameters,
+			ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
 
 	int i;
@@ -31,5 +32,5 @@
 
 	/*Call on dataset driver: */
-	elements->ComputePressure(p_g);
+	elements->ComputePressure(p_g,inputs,analysis_type,sub_analysis_type);
 
 	/*Assemble vector: */
Index: /issm/trunk/src/c/ComputePressurex/ComputePressurex.h
===================================================================
--- /issm/trunk/src/c/ComputePressurex/ComputePressurex.h	(revision 3528)
+++ /issm/trunk/src/c/ComputePressurex/ComputePressurex.h	(revision 3529)
@@ -7,7 +7,9 @@
 
 #include "../DataSet/DataSet.h"
+#include "../objects/objects.h"
 
 /* local prototypes: */
-void	ComputePressurex( Vec* pp_g,DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials,  DataSet* parameters);
+void	ComputePressurex( Vec* pp_g,DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials,  DataSet* parameters,
+			ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
 
 #endif  /* _COMPUTEPRESSUREX_H */
Index: /issm/trunk/src/c/ComputeStrainRatex/ComputeStrainRatex.cpp
===================================================================
--- /issm/trunk/src/c/ComputeStrainRatex/ComputeStrainRatex.cpp	(revision 3529)
+++ /issm/trunk/src/c/ComputeStrainRatex/ComputeStrainRatex.cpp	(revision 3529)
@@ -0,0 +1,44 @@
+/*!\file ComputeStrainRatex
+ * \brief: compute pressure according to each element
+ */
+
+#include "./ComputeStrainRatex.h"
+
+#include "../shared/shared.h"
+#include "../include/macros.h"
+#include "../toolkits/toolkits.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+void	ComputeStrainRatex( Vec* peps,DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials,DataSet* parameters,
+			ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
+
+	int i;
+
+	int  found=0;
+	double numberofelements;
+
+	/*output: */
+	Vec eps=NULL;
+
+	/*Recover numberofelements: */
+	parameters->FindParam(&numberofelements,"numberofelements");
+
+	/*Allocate eps on numberofelements (only 1 dof): */
+	eps=NewVec(numberofelements);
+
+	/*Get elements configured: */
+	elements->Configure(elements,loads, nodes,vertices, materials,parameters);
+	nodes->Configure(elements,loads, nodes,vertices, materials,parameters);
+	parameters->Configure(elements,loads, nodes,vertices, materials,parameters);
+
+	/*Call on dataset driver: */
+	elements->ComputeStrainRate(eps,inputs,analysis_type,sub_analysis_type);
+
+	/*Assemble vector: */
+	VecAssemblyBegin(eps);
+	VecAssemblyEnd(eps);
+
+	/*Assign output pointers: */
+	*peps=eps;
+	
+}
Index: /issm/trunk/src/c/ComputeStrainRatex/ComputeStrainRatex.h
===================================================================
--- /issm/trunk/src/c/ComputeStrainRatex/ComputeStrainRatex.h	(revision 3529)
+++ /issm/trunk/src/c/ComputeStrainRatex/ComputeStrainRatex.h	(revision 3529)
@@ -0,0 +1,16 @@
+/*!\file:  ComputeStrainRate.h
+ * \brief header file pressure computation
+ */ 
+
+#ifndef _COMPUTESTRAINRATEX_H
+#define _COMPUTESTRAINRATEX_H
+
+#include "../objects/objects.h"
+#include "../DataSet/DataSet.h"
+
+/* local prototypes: */
+void	ComputeStrainRatex(Vec* eps_g,DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials,  DataSet* parameters,
+			ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
+
+#endif  /* _COMPUTESTRAINRATEX_H */
+
Index: /issm/trunk/src/c/DataSet/DataSet.cpp
===================================================================
--- /issm/trunk/src/c/DataSet/DataSet.cpp	(revision 3528)
+++ /issm/trunk/src/c/DataSet/DataSet.cpp	(revision 3529)
@@ -789,6 +789,23 @@
 
 /*Objects methods*/
+/*FUNCTION DataSet::ComputeBasalStress{{{1*/
+void DataSet::ComputeBasalStress(Vec sigma_b,void* inputs,int analysis_type,int sub_analysis_type){
+
+	vector<Object*>::iterator object;
+	Element* element=NULL;
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		if(EnumIsElement((*object)->Enum())){
+
+			element=(Element*)(*object);
+			element->ComputeBasalStress(sigma_b,inputs,analysis_type,sub_analysis_type);
+		}
+	}
+
+}
+/*}}}*/
 /*FUNCTION DataSet::ComputePressure{{{1*/
-void DataSet::ComputePressure(Vec p_g){
+void DataSet::ComputePressure(Vec p_g,void* inputs,int analysis_type,int sub_analysis_type){
 
 	vector<Object*>::iterator object;
@@ -800,5 +817,22 @@
 
 			element=(Element*)(*object);
-			element->ComputePressure(p_g);
+			element->ComputePressure(p_g,inputs,analysis_type,sub_analysis_type);
+		}
+	}
+
+}
+/*}}}*/
+/*FUNCTION DataSet::ComputeStrainRate{{{1*/
+void DataSet::ComputeStrainRate(Vec eps,void* inputs,int analysis_type,int sub_analysis_type){
+
+	vector<Object*>::iterator object;
+	Element* element=NULL;
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		if(EnumIsElement((*object)->Enum())){
+
+			element=(Element*)(*object);
+			element->ComputeStrainRate(eps,inputs,analysis_type,sub_analysis_type);
 		}
 	}
Index: /issm/trunk/src/c/DataSet/DataSet.h
===================================================================
--- /issm/trunk/src/c/DataSet/DataSet.h	(revision 3528)
+++ /issm/trunk/src/c/DataSet/DataSet.h	(revision 3529)
@@ -91,5 +91,7 @@
 		void  FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname);
 		int   DeleteObject(Object* object);
-		void  ComputePressure(Vec p_g);
+		void  ComputeBasalStress(Vec sigma_b,void* inputs,int analysis_type,int sub_analysis_type);
+		void  ComputePressure(Vec p_g,void* inputs,int analysis_type,int sub_analysis_type);
+		void  ComputeStrainRate(Vec eps,void* inputs,int analysis_type,int sub_analysis_type);
 		int   FindResult(void* pvalue, char* name);
 		void  FieldExtrude(Vec field,double* field_serial,char* field_name, int collapse);
Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 3528)
+++ /issm/trunk/src/c/Makefile.am	(revision 3529)
@@ -274,6 +274,10 @@
 					./ConfigureObjectsx/ConfigureObjectsx.h\
 					./ConfigureObjectsx/ConfigureObjectsx.cpp\
+					./ComputeBasalStressx/ComputeBasalStressx.h\
+					./ComputeBasalStressx/ComputeBasalStressx.cpp\
 					./ComputePressurex/ComputePressurex.h\
 					./ComputePressurex/ComputePressurex.cpp\
+					./ComputeStrainRatex/ComputeStrainRatex.h\
+					./ComputeStrainRatex/ComputeStrainRatex.cpp\
 					./BuildNodeSetsx/BuildNodeSetsx.h\
 					./BuildNodeSetsx/BuildNodeSetsx.cpp\
@@ -670,6 +674,10 @@
 					./ConfigureObjectsx/ConfigureObjectsx.h\
 					./ConfigureObjectsx/ConfigureObjectsx.cpp\
+					./ComputeBasalStressx/ComputeBasalStressx.h\
+					./ComputeBasalStressx/ComputeBasalStressx.cpp\
 					./ComputePressurex/ComputePressurex.h\
 					./ComputePressurex/ComputePressurex.cpp\
+					./ComputeStrainRatex/ComputeStrainRatex.h\
+					./ComputeStrainRatex/ComputeStrainRatex.cpp\
 					./BuildNodeSetsx/BuildNodeSetsx.h\
 					./BuildNodeSetsx/BuildNodeSetsx.cpp\
Index: /issm/trunk/src/c/issm.h
===================================================================
--- /issm/trunk/src/c/issm.h	(revision 3528)
+++ /issm/trunk/src/c/issm.h	(revision 3529)
@@ -58,5 +58,7 @@
 #include "./FieldAverageOntoVerticesx/FieldAverageOntoVerticesx.h"
 #include "./FieldDepthAveragex/FieldDepthAveragex.h"
+#include "./ComputeBasalStressx/ComputeBasalStressx.h"
 #include "./ComputePressurex/ComputePressurex.h"
+#include "./ComputeStrainRatex/ComputeStrainRatex.h"
 #include "./FieldExtrudex/FieldExtrudex.h"
 #include "./Qmux/Qmux.h"
@@ -68,4 +70,3 @@
 #include "./BamgConvertMeshx/BamgConvertMeshx.h"
 
-
 #endif
Index: /issm/trunk/src/c/objects/Beam.cpp
===================================================================
--- /issm/trunk/src/c/objects/Beam.cpp	(revision 3528)
+++ /issm/trunk/src/c/objects/Beam.cpp	(revision 3529)
@@ -306,6 +306,13 @@
 
 /*Object functions*/
+/*FUNCTION Beam::ComputeBasalStress{{{1*/
+void  Beam::ComputeBasalStress(Vec eps,void* inputs,int analysis_type,int sub_analysis_type){
+
+	ISSMERROR("Not implemented yet");
+
+}
+/*}}}*/
 /*FUNCTION Beam::ComputePressure{{{1*/
-void  Beam::ComputePressure(Vec p_g){
+void  Beam::ComputePressure(Vec p_g,void* inputs,int analysis_type,int sub_analysis_type){
 
 	int i;
@@ -346,4 +353,11 @@
 	/*plug local pressure values into global pressure vector: */
 	VecSetValues(p_g,numgrids,doflist,(const double*)pressure,INSERT_VALUES);
+
+}
+/*}}}*/
+/*FUNCTION Beam::ComputeStrainRate{{{1*/
+void  Beam::ComputeStrainRate(Vec eps,void* inputs,int analysis_type,int sub_analysis_type){
+
+	ISSMERROR("Not implemented yet");
 
 }
Index: /issm/trunk/src/c/objects/Beam.h
===================================================================
--- /issm/trunk/src/c/objects/Beam.h	(revision 3528)
+++ /issm/trunk/src/c/objects/Beam.h	(revision 3529)
@@ -74,5 +74,7 @@
 		void  MaticeConfiguration(Matice* matice,int matice_offset);
 		void  MatparConfiguration(Matpar* matpar,int matpar_offset);
-		void  ComputePressure(Vec p_g);
+		void  ComputeBasalStress(Vec sigma_b,void* inputs,int analysis_type,int sub_analysis_type);
+		void  ComputePressure(Vec p_g,void* inputs,int analysis_type,int sub_analysis_type);
+		void  ComputeStrainRate(Vec eps,void* inputs,int analysis_type,int sub_analysis_type);
 		void  GetNodes(void** vpnodes);
 		/*}}}*/
Index: /issm/trunk/src/c/objects/Element.h
===================================================================
--- /issm/trunk/src/c/objects/Element.h	(revision 3528)
+++ /issm/trunk/src/c/objects/Element.h	(revision 3529)
@@ -34,5 +34,7 @@
 		virtual double CostFunction(void* inputs,int analysis_type,int sub_analysis_type)=0;
 		virtual double SurfaceArea(void* inputs,int analysis_type,int sub_analysis_type)=0;
-		virtual void   ComputePressure(Vec p_g)=0;
+		virtual void   ComputeBasalStress(Vec sigma_b,void* inputs,int analysis_type,int sub_analysis_type)=0;
+		virtual void   ComputePressure(Vec p_g,       void* inputs,int analysis_type,int sub_analysis_type)=0;
+		virtual void   ComputeStrainRate(Vec eps,     void* inputs,int analysis_type,int sub_analysis_type)=0;
 		virtual double MassFlux(double* segment,double* ug)=0;
 
Index: /issm/trunk/src/c/objects/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Penta.cpp	(revision 3528)
+++ /issm/trunk/src/c/objects/Penta.cpp	(revision 3529)
@@ -336,6 +336,138 @@
 
 /*Object functions*/
+/*FUNCTION ComputeBasalStress {{{1*/
+void  Penta::ComputeBasalStress(Vec sigma_b,void* vinputs,int analysis_type,int sub_analysis_type){
+
+	int i,j;
+	const int numgrids=6;
+	int    doflist[numgrids];
+	double xyz_list[numgrids][3];
+	double xyz_list_tria[3][3];
+
+	/*Parameters*/
+	double rho_ice,gravity;
+	double surface_normal[3];
+	double bed_normal[3];
+	double bed;
+	double basalforce[3];
+	double vxvyvz_list[numgrids][3];
+	double pressure_list[numgrids];
+	double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
+	double devstresstensor[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
+	double stresstensor[6]={0.0};
+	double viscosity;
+
+	int  dofv[3]={0,1,2};
+	int  dofp[1]={3};
+	ParameterInputs* inputs=NULL;
+	double Jdet2d;
+	Tria* tria=NULL;
+
+	/*Gauss*/
+	int     num_gauss,ig;
+	double* first_gauss_area_coord  =  NULL;
+	double* second_gauss_area_coord =  NULL;
+	double* third_gauss_area_coord  =  NULL;
+	double* gauss_weights           =  NULL;
+	double  gauss_weight;
+	double  gauss_coord[4];
+
+	/*Output*/
+	double pressure;
+	double sigma_xx,sigma_yy,sigma_zz;
+	double sigma_xy,sigma_xz,sigma_yz;
+	double surface=0;
+	double value=0;
+
+	/*dynamic objects pointed to by hooks: */
+	Node**  nodes=NULL;
+	Matpar* matpar=NULL;
+	Matice* matice=NULL;
+	Numpar* numpar=NULL;
+
+	/*Check analysis_types*/
+	if (analysis_type!=DiagnosticAnalysisEnum() || sub_analysis_type!=StokesAnalysisEnum()) ISSMERROR("Not supported yet!");
+
+	/*recover pointers: */
+	inputs=(ParameterInputs*)vinputs;
+
+	/*recover objects from hooks: */
+	nodes=(Node**)hnodes.deliverp();
+	matpar=(Matpar*)hmatpar.delivers();
+	matice=(Matice*)hmatice.delivers();
+	numpar=(Numpar*)hnumpar.delivers();
+
+	if(!this->properties.onbed){
+		//put zero
+		VecSetValue(sigma_b,id-1,0.0,INSERT_VALUES);
+		return;
+	}
+
+	/*recovre material parameters: */
+	rho_ice=matpar->GetRhoIce();
+	gravity=matpar->GetG();
+
+	/*recover extra inputs from users, at current convergence iteration: */
+	inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofv,numgrids,(void**)nodes);
+	inputs->Recover("velocity",&pressure_list[0] ,1,dofp,numgrids,(void**)nodes);
+
+	/* Get node coordinates and dof list: */
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+	for(i=0;i<3;i++){
+		for(j=0;j<3;j++){
+			xyz_list_tria[i][j]=xyz_list[i][j];
+		}
+	}
+
+	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+	GaussTria(&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights,2);
+
+	/* Start  looping on the number of gaussian points: */
+	for (ig=0; ig<num_gauss; ig++){
+
+			/*Pick up the gaussian point: */
+			gauss_weight=*(gauss_weights+ig);
+			gauss_coord[0]=*(first_gauss_area_coord+ig); 
+			gauss_coord[1]=*(second_gauss_area_coord+ig);
+			gauss_coord[2]=*(third_gauss_area_coord+ig);
+			gauss_coord[3]=-1.0; //take the base
+
+			/*Compute strain rate viscosity and pressure: */
+			GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
+			matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+			GetParameterValue(&pressure,&pressure_list[0],&gauss_coord[0]);
+
+			/*Compute Stress*/
+			sigma_xx=viscosity*epsilon[0]-pressure*numpar->stokesreconditioning; // sigma = nu eps - pressure
+			sigma_yy=viscosity*epsilon[1]-pressure*numpar->stokesreconditioning;
+			sigma_zz=viscosity*epsilon[2]-pressure*numpar->stokesreconditioning;
+			sigma_xy=viscosity*epsilon[3];
+			sigma_xz=viscosity*epsilon[4];
+			sigma_yz=viscosity*epsilon[5];
+
+			/*Get normal vector to the bed */
+			SurfaceNormal(&surface_normal[0],xyz_list_tria);
+			bed_normal[0] = - surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result
+			bed_normal[1] = - surface_normal[1];
+			bed_normal[2] = - surface_normal[2];
+
+			/*basalforce*/
+			basalforce[0] += sigma_xx*bed_normal[0] + sigma_xy*bed_normal[1] + sigma_xz*bed_normal[2];
+			basalforce[1] += sigma_xy*bed_normal[0] + sigma_yy*bed_normal[1] + sigma_yz*bed_normal[2];
+			basalforce[2] += sigma_xz*bed_normal[0] + sigma_yz*bed_normal[1] + sigma_zz*bed_normal[2];
+
+			/*Get the Jacobian determinant */
+			tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0],gauss_coord);
+			value+=sigma_zz*Jdet2d*gauss_weight;
+			surface+=Jdet2d*gauss_weight;
+	}
+	value=value/surface;
+
+	/*Add value to output*/
+	VecSetValue(sigma_b,id-1,(const double)value,INSERT_VALUES);
+}
+/*}}}*/
 /*FUNCTION ComputePressure {{{1*/
-void  Penta::ComputePressure(Vec pg){
+void  Penta::ComputePressure(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type){
 
 	int i;
@@ -375,4 +507,11 @@
 	/*plug local pressure values into global pressure vector: */
 	VecSetValues(pg,numgrids,doflist,(const double*)pressure,INSERT_VALUES);
+
+}
+/*}}}*/
+/*FUNCTION ComputeStrainRate {{{1*/
+void  Penta::ComputeStrainRate(Vec eps,void* vinputs,int analysis_type,int sub_analysis_type){
+
+	ISSMERROR("Not implemented yet");
 
 }
@@ -4057,5 +4196,5 @@
 }
 /*}}}*/
-/**FUNCTION ReduceVectorStokes {{{1*/
+/*FUNCTION ReduceVectorStokes {{{1*/
 void Penta::ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp){
 
@@ -4099,5 +4238,5 @@
 }
 /*}}}*/
-/*FUNCTION Penta::SetClone {{{1*/
+/*FUNCTION SetClone {{{1*/
 void  Penta::SetClone(int* minranks){
 
Index: /issm/trunk/src/c/objects/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Penta.h	(revision 3528)
+++ /issm/trunk/src/c/objects/Penta.h	(revision 3529)
@@ -103,5 +103,7 @@
 		void  GetNodalFunctions(double* l1l6, double* gauss_coord);
 		void  FieldExtrude(Vec field,double* field_serial,char* field_name, int iscollapsed);
-		void  ComputePressure(Vec p_g);
+		void  ComputeBasalStress(Vec sigma_b,void* inputs,int analysis_type,int sub_analysis_type);
+		void  ComputePressure(Vec p_g,void* inputs,int analysis_type,int sub_analysis_type);
+		void  ComputeStrainRate(Vec eps,void* inputs,int analysis_type,int sub_analysis_type);
 		void  CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
 		void  CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
Index: /issm/trunk/src/c/objects/Sing.cpp
===================================================================
--- /issm/trunk/src/c/objects/Sing.cpp	(revision 3528)
+++ /issm/trunk/src/c/objects/Sing.cpp	(revision 3529)
@@ -290,6 +290,13 @@
 		
 /*Object functions*/
+/*FUNCTION Sing::ComputeBasalStress {{{1*/
+void  Sing::ComputeBasalStress(Vec p_g,void* inputs,int analysis_type,int sub_analysis_type){
+
+	ISSMERROR("Not implemented yet");
+
+}
+/*}}}*/
 /*FUNCTION Sing::ComputePressure {{{1*/
-void  Sing::ComputePressure(Vec p_g){
+void  Sing::ComputePressure(Vec p_g,void* inputs,int analysis_type,int sub_analysis_type){
 
 	int i;
@@ -321,4 +328,11 @@
 	/*plug local pressure values into global pressure vector: */
 	VecSetValues(p_g,numgrids,doflist,(const double*)pressure,INSERT_VALUES);
+
+}
+/*}}}*/
+/*FUNCTION Sing::ComputeStrainRate {{{1*/
+void  Sing::ComputeStrainRate(Vec p_g,void* inputs,int analysis_type,int sub_analysis_type){
+
+	ISSMERROR("Not implemented yet");
 
 }
Index: /issm/trunk/src/c/objects/Sing.h
===================================================================
--- /issm/trunk/src/c/objects/Sing.h	(revision 3528)
+++ /issm/trunk/src/c/objects/Sing.h	(revision 3529)
@@ -68,5 +68,7 @@
 		void  MaticeConfiguration(Matice* matice,int matice_offset);
 		void  MatparConfiguration(Matpar* matpar,int matpar_offset);
-		void  ComputePressure(Vec p_g);
+		void  ComputeBasalStress(Vec sigma_b,void* inputs,int analysis_type,int sub_analysis_type);
+		void  ComputePressure(Vec p_g,void* inputs,int analysis_type,int sub_analysis_type);
+		void  ComputeStrainRate(Vec eps,void* inputs,int analysis_type,int sub_analysis_type);
 		void  GetNodes(void** vpnodes);
 		/*}}}*/
Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 3528)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 3529)
@@ -349,6 +349,32 @@
 
 /*Object functions*/
+/*FUNCTION Tria::ComputeBasalStress {{{1*/
+void  Tria::ComputeBasalStress(Vec eps,void* vinputs,int analysis_type,int sub_analysis_type){
+
+	int i;
+	const int numgrids=3;
+	int doflist[numgrids];
+	double value;
+
+	/*dynamic objects pointed to by hooks: */
+	Node**  nodes=NULL;
+	Matpar* matpar=NULL;
+	Matice* matice=NULL;
+	Numpar* numpar=NULL;
+
+	/*recover objects from hooks: */
+	nodes=(Node**)hnodes.deliverp();
+	matpar=(Matpar*)hmatpar.delivers();
+	matice=(Matice*)hmatice.delivers();
+	numpar=(Numpar*)hnumpar.delivers();
+
+	/*plug local pressure values into global pressure vector: */
+	ISSMERROR("Not Implemented yet");
+	//VecSetValues(eps,1,2,(const double*)&value,INSERT_VALUES);
+
+}
+/*}}}*/
 /*FUNCTION Tria::ComputePressure {{{1*/
-void  Tria::ComputePressure(Vec pg){
+void  Tria::ComputePressure(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type){
 
 	int i;
@@ -383,4 +409,30 @@
 	/*plug local pressure values into global pressure vector: */
 	VecSetValues(pg,numgrids,doflist,(const double*)pressure,INSERT_VALUES);
+
+}
+/*}}}*/
+/*FUNCTION Tria::ComputeStrainRate {{{1*/
+void  Tria::ComputeStrainRate(Vec eps,void* vinputs,int analysis_type,int sub_analysis_type){
+
+	int i;
+	const int numgrids=3;
+	int doflist[numgrids];
+	double value;
+
+	/*dynamic objects pointed to by hooks: */
+	Node**  nodes=NULL;
+	Matpar* matpar=NULL;
+	Matice* matice=NULL;
+	Numpar* numpar=NULL;
+
+	/*recover objects from hooks: */
+	nodes=(Node**)hnodes.deliverp();
+	matpar=(Matpar*)hmatpar.delivers();
+	matice=(Matice*)hmatice.delivers();
+	numpar=(Numpar*)hnumpar.delivers();
+
+	/*plug local pressure values into global pressure vector: */
+	ISSMERROR("Not Implemented yet");
+	//VecSetValues(eps,1,2,(const double*)&value,INSERT_VALUES);
 
 }
Index: /issm/trunk/src/c/objects/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Tria.h	(revision 3528)
+++ /issm/trunk/src/c/objects/Tria.h	(revision 3529)
@@ -104,5 +104,7 @@
 		void  GetThicknessList(double* thickness_list);
 		void  GetBedList(double* bed_list);
-		void  ComputePressure(Vec p_g);
+		void  ComputeBasalStress(Vec sigma_b,void* inputs,int analysis_type,int sub_analysis_type);
+		void  ComputePressure(Vec p_g,void* inputs,int analysis_type,int sub_analysis_type);
+		void  ComputeStrainRate(Vec eps,void* inputs,int analysis_type,int sub_analysis_type);
 		void  CreateKMatrixThermal(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
 		void  CreateKMatrixMelting(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
Index: /issm/trunk/src/c/parallel/ControlInitialization.cpp
===================================================================
--- /issm/trunk/src/c/parallel/ControlInitialization.cpp	(revision 3528)
+++ /issm/trunk/src/c/parallel/ControlInitialization.cpp	(revision 3529)
@@ -113,5 +113,5 @@
 	//Create 4d u_g
 	if(verbose)_printf_("%s\n"," computing pressure according to Pattyn...");
-	ComputePressurex( &pg,fem_dh->elements, fem_dh->nodes, fem_dh->vertices,fem_dh->loads,  fem_dh->materials, fem_dh->parameters);
+	ComputePressurex( &pg,fem_dh->elements, fem_dh->nodes, fem_dh->vertices,fem_dh->loads,  fem_dh->materials, fem_dh->parameters,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
 	VecScale(pg,1.0/stokesreconditioning);
 	ug_stokes=NewVec(fem_ds->nodesets->GetGSize());
Index: /issm/trunk/src/c/parallel/ProcessResults.cpp
===================================================================
--- /issm/trunk/src/c/parallel/ProcessResults.cpp	(revision 3528)
+++ /issm/trunk/src/c/parallel/ProcessResults.cpp	(revision 3529)
@@ -100,8 +100,11 @@
 	double* gradient=NULL;
 
+	Vec     sigma_zz=NULL;
+	double* sigma_zz_serial=NULL;
+
 	Vec     v_g=NULL;
 	double* v_g_serial=NULL;
 
-	int numberofnodes;
+	int numberofnodes,numberofelements;
 
 	/*Initialize new results: */
@@ -516,4 +519,18 @@
 
 		}
+		else if(strcmp(result->GetFieldName(),"sigma_zz")==0){
+			/*easy, param_g is of size numberofelements, on 1 dof, just repartition: */
+			fem_ds->parameters->FindParam(&numberofelements,"numberofelements");
+			result->GetField(&sigma_zz);
+			VecToMPISerial(&sigma_zz_serial,sigma_zz);
+
+			/*Ok, add parameter to newresults: */
+			newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"sigma_zz",sigma_zz_serial,numberofelements);
+			newresults->AddObject(newresult);
+
+			/*do some cleanup: */
+			xfree((void**)&sigma_zz_serial);
+
+		}
 		else{
 			/*Just copy the result into the new results dataset: */
Index: /issm/trunk/src/c/parallel/diagnostic_core.cpp
===================================================================
--- /issm/trunk/src/c/parallel/diagnostic_core.cpp	(revision 3528)
+++ /issm/trunk/src/c/parallel/diagnostic_core.cpp	(revision 3529)
@@ -109,5 +109,5 @@
 
 		if(verbose)_printf_("%s\n"," computing pressure according to MacAyeal...");
-		ComputePressurex( &pg,fem_dhu->elements,fem_dhu->nodes,fem_dhu->vertices,fem_dhu->loads,fem_dhu->materials, fem_dhu->parameters);
+		ComputePressurex(&pg,fem_dhu->elements,fem_dhu->nodes, fem_dhu->vertices,fem_dhu->loads,fem_dhu->materials,fem_dhu->parameters,inputs,DiagnosticAnalysisEnum(),HutterAnalysisEnum());
 
 		if(verbose)_printf_("%s\n"," update boundary conditions for macyeal pattyn using hutter results...");
@@ -127,5 +127,5 @@
 		if(dim==2){
 			if(verbose)_printf_("%s\n"," computing pressure according to MacAyeal...");
-			ComputePressurex( &pg,fem_dh->elements, fem_dh->nodes, fem_dh->vertices,fem_dh->loads,  fem_dh->materials, fem_dh->parameters);
+			ComputePressurex(&pg,fem_dh->elements,fem_dh->nodes, fem_dh->vertices,fem_dh->loads,fem_dh->materials,fem_dh->parameters,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
 		}
 
@@ -150,5 +150,5 @@
 
 		if(verbose)_printf_("%s\n"," computing pressure according to Pattyn...");
-		ComputePressurex( &pg,fem_dh->elements, fem_dh->nodes, fem_dh->vertices,fem_dh->loads,  fem_dh->materials,fem_dh->parameters);
+		ComputePressurex(&pg,fem_dh->elements, fem_dh->nodes, fem_dh->vertices,fem_dh->loads,  fem_dh->materials,fem_dh->parameters,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
 		
 		if (isstokes){
Index: /issm/trunk/src/m/solutions/jpl/ControlInitialization.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/ControlInitialization.m	(revision 3528)
+++ /issm/trunk/src/m/solutions/jpl/ControlInitialization.m	(revision 3529)
@@ -55,5 +55,5 @@
 % get pressure (reconditionned) and create 4d u_g
 displaystring(verbose,'\n%s',['computing pressure according to Pattyn...']);
-p_g=ComputePressure(m_dh.elements,m_dh.nodes,mdh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,inputs);
+p_g=ComputePressure(m_dh.elements,m_dh.nodes,mdh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
 p_g=p_g/m_ds.parameters.stokesreconditioning;
 u_g_stokes=zeros(m_ds.nodesets.gsize,1);
Index: /issm/trunk/src/m/solutions/jpl/diagnostic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/diagnostic_core.m	(revision 3528)
+++ /issm/trunk/src/m/solutions/jpl/diagnostic_core.m	(revision 3529)
@@ -43,5 +43,5 @@
 
 	displaystring(verbose,'\n%s',['computing pressure according to MacAyeal...']);
-	p_g=ComputePressure(m_dhu.elements,m_dhu.nodes,m_dhu.vertices,m_dhu.loads,m_dhu.materials,m_dhu.parameters,inputs);
+	p_g=ComputePressure(m_dhu.elements,m_dhu.nodes,m_dhu.vertices,m_dhu.loads,m_dhu.materials,m_dhu.parameters,inputs,DiagnosticAnalysisEnum(),HutterAnalysisEnum());
 
 	displaystring(verbose,'\n%s',['update boundary conditions for macyeal pattyn using hutter results...']);
@@ -59,5 +59,5 @@
 
 	displaystring(verbose,'\n%s',['computing pressure according to MacAyeal...']);
-	p_g=ComputePressure(m_dh.elements,m_dh.nodes,m_dh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,inputs);
+	p_g=ComputePressure(m_dh.elements,m_dh.nodes,m_dh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
 end
 	
@@ -77,5 +77,5 @@
 
 	displaystring(verbose,'\n%s',['computing pressure according to Pattyn...']);
-	p_g=ComputePressure(m_dh.elements,m_dh.nodes,m_dh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,inputs);
+	p_g=ComputePressure(m_dh.elements,m_dh.nodes,m_dh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
 	
 	if isstokes,
Index: /issm/trunk/src/mex/ComputePressure/ComputePressure.cpp
===================================================================
--- /issm/trunk/src/mex/ComputePressure/ComputePressure.cpp	(revision 3528)
+++ /issm/trunk/src/mex/ComputePressure/ComputePressure.cpp	(revision 3529)
@@ -19,4 +19,6 @@
 	ParameterInputs* inputs=NULL;
 	int      numberofnodes;
+	int      analysis_type;
+	int      sub_analysis_type;
 
 	/* output datasets: */
@@ -36,4 +38,6 @@
 	FetchData(&materials,MATERIALS);
 	FetchParams(&parameters,PARAMETERS);
+	FetchData(&analysis_type,ANALYSIS);
+	FetchData(&sub_analysis_type,SUBANALYSIS);
 
 	/*!Generate internal degree of freedom numbers: */
@@ -45,5 +49,5 @@
 
 	/*!Generate internal degree of freedom numbers: */
-	ComputePressurex(&p_g, elements,nodes,vertices,loads,materials,parameters);
+	ComputePressurex(&p_g, elements,nodes,vertices,loads,materials,parameters,inputs,analysis_type,sub_analysis_type);
 
 	/*write output datasets: */
Index: /issm/trunk/src/mex/ComputePressure/ComputePressure.h
===================================================================
--- /issm/trunk/src/mex/ComputePressure/ComputePressure.h	(revision 3528)
+++ /issm/trunk/src/mex/ComputePressure/ComputePressure.h	(revision 3529)
@@ -24,4 +24,6 @@
 #define PARAMETERS (mxArray*)prhs[5]
 #define INPUTS (mxArray*)prhs[6]
+#define ANALYSIS (mxArray*)prhs[7]
+#define SUBANALYSIS (mxArray*)prhs[8]
 
 /* serial output macros: */
@@ -32,5 +34,5 @@
 #define NLHS  1
 #undef NRHS
-#define NRHS  7
+#define NRHS  9
 
 
Index: /issm/trunk/src/mex/Makefile.am
===================================================================
--- /issm/trunk/src/mex/Makefile.am	(revision 3528)
+++ /issm/trunk/src/mex/Makefile.am	(revision 3529)
@@ -9,4 +9,5 @@
 				BamgConvertMesh\
 				BuildNodeSets\
+				ComputeBasalStress\
 				ComputePressure\
 				ConfigureObjects \
@@ -93,4 +94,7 @@
 					BamgConvertMesh/BamgConvertMesh.h
 
+ComputeBasalStress_SOURCES = ComputeBasalStress/ComputeBasalStress.cpp\
+								  ComputeBasalStress/ComputeBasalStress.h
+
 ComputePressure_SOURCES = ComputePressure/ComputePressure.cpp\
 			  ComputePressure/ComputePressure.h