Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17547)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17548)
@@ -4257,6 +4257,7 @@
 	IssmDouble  dvx[3],dvy[3],dvz[3],B,n;
 	IssmDouble *xyz_list = NULL;
 	IssmDouble  Jdet,r;
+	Gauss*      gauss = NULL;
 
 	parameters->FindParam(&meshtype,MeshTypeEnum);
 	switch(meshtype){
@@ -4308,7 +4309,7 @@
 			sigmapyz_input=element->GetInput(DeviatoricStressyzEnum); _assert_(sigmapyz_input);
 		}
 
-		Gauss* gauss=element->NewGauss(5);
+		gauss=element->NewGauss(5);
 		for(int ig=gauss->begin();ig<gauss->end();ig++){
 			gauss->GaussPoint(ig);
 			element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -4373,30 +4374,128 @@
 						Ke,1);
 
 			/*Create Right hand sides*/
-			for(int ii=0;ii<tnumnodes;ii++) pe_xx[ii] += sigmapxx*tbasis[ii]*gauss->weight*Jdet;
-			for(int ii=0;ii<tnumnodes;ii++) pe_yy[ii] += sigmapyy*tbasis[ii]*gauss->weight*Jdet;
-			for(int ii=0;ii<tnumnodes;ii++) pe_xy[ii] += sigmapxy*tbasis[ii]*gauss->weight*Jdet;
+			for(int ii=0;ii<tnumnodes;ii++) pe_xx[ii] += (r*epsxx+sigmapxx)*tbasis[ii]*gauss->weight*Jdet;
+			for(int ii=0;ii<tnumnodes;ii++) pe_yy[ii] += (r*epsyy+sigmapyy)*tbasis[ii]*gauss->weight*Jdet;
+			for(int ii=0;ii<tnumnodes;ii++) pe_xy[ii] += (r*epsxy+sigmapxy)*tbasis[ii]*gauss->weight*Jdet;
 			if(dim==3){
-				for(int ii=0;ii<tnumnodes;ii++) pe_zz[ii] += sigmapzz*tbasis[ii]*gauss->weight*Jdet;
-				for(int ii=0;ii<tnumnodes;ii++) pe_xz[ii] += sigmapxz*tbasis[ii]*gauss->weight*Jdet;
-				for(int ii=0;ii<tnumnodes;ii++) pe_yz[ii] += sigmapyz*tbasis[ii]*gauss->weight*Jdet;
+				for(int ii=0;ii<tnumnodes;ii++) pe_zz[ii] += (r*epszz+sigmapzz)*tbasis[ii]*gauss->weight*Jdet;
+				for(int ii=0;ii<tnumnodes;ii++) pe_xz[ii] += (r*epsxz+sigmapxz)*tbasis[ii]*gauss->weight*Jdet;
+				for(int ii=0;ii<tnumnodes;ii++) pe_yz[ii] += (r*epsyz+sigmapyz)*tbasis[ii]*gauss->weight*Jdet;
 			}
 		}
 
 		/*Solve the systems*/
-		_error_("S");
+		IssmDouble* d_xx = xNew<IssmDouble>(tnumnodes);
+		IssmDouble* d_yy = xNew<IssmDouble>(tnumnodes);
+		IssmDouble* d_xy = xNew<IssmDouble>(tnumnodes);
+		IssmDouble* d_zz = NULL;
+		IssmDouble* d_xz = NULL;
+		IssmDouble* d_yz = NULL;
+		if(dim==2){
+			_assert_(tnumnodes==3);
+			Matrix3x3Solve(&d_xx[0],Ke,pe_xx);
+			Matrix3x3Solve(&d_yy[0],Ke,pe_yy);
+			Matrix3x3Solve(&d_xy[0],Ke,pe_xy);
+			element->AddInput(StrainRatexxEnum,d_xx,P1DGEnum);
+			element->AddInput(StrainRateyyEnum,d_yy,P1DGEnum);
+			element->AddInput(StrainRatexyEnum,d_xy,P1DGEnum);
+		}
+		else{
+			_assert_(tnumnodes==4);
+			Matrix3x3Solve(&d_xx[0],Ke,pe_xx);
+			Matrix3x3Solve(&d_yy[0],Ke,pe_yy);
+			Matrix3x3Solve(&d_xy[0],Ke,pe_xy);
+			Matrix3x3Solve(&d_zz[0],Ke,pe_zz);
+			Matrix3x3Solve(&d_xz[0],Ke,pe_xz);
+			Matrix3x3Solve(&d_yz[0],Ke,pe_yz);
+			element->AddInput(StrainRatexxEnum,d_xx,P1DGEnum);
+			element->AddInput(StrainRateyyEnum,d_yy,P1DGEnum);
+			element->AddInput(StrainRatexyEnum,d_xy,P1DGEnum);
+			element->AddInput(StrainRatezzEnum,d_zz,P1DGEnum);
+			element->AddInput(StrainRatexzEnum,d_xz,P1DGEnum);
+			element->AddInput(StrainRateyzEnum,d_yz,P1DGEnum);
+		}
 
+		/*Update tau accordingly*/
+		IssmDouble* tau_xx = xNew<IssmDouble>(tnumnodes);
+		IssmDouble* tau_yy = xNew<IssmDouble>(tnumnodes);
+		IssmDouble* tau_xy = xNew<IssmDouble>(tnumnodes);
+		IssmDouble* tau_zz = NULL;
+		IssmDouble* tau_xz = NULL;
+		IssmDouble* tau_yz = NULL;
+		Input* epsxx_input=element->GetInput(StrainRatexxEnum); _assert_(epsxx_input);
+		Input* epsyy_input=element->GetInput(StrainRateyyEnum); _assert_(epsyy_input);
+		Input* epsxy_input=element->GetInput(StrainRatexyEnum); _assert_(epsxy_input);
+		Input* epszz_input=NULL; Input* epsxz_input=NULL; Input* epsyz_input=NULL;
+		if(dim==3){
+			epszz_input=element->GetInput(StrainRatezzEnum); _assert_(epszz_input);
+			epsxz_input=element->GetInput(StrainRatexzEnum); _assert_(epsxz_input);
+			epsyz_input=element->GetInput(StrainRateyzEnum); _assert_(epsyz_input);
+		}
+		Gauss* gauss = element->NewGauss();
+		for(int ig=0;ig<tnumnodes;ig++){
+			gauss->GaussNode(P1DGEnum,ig);
+
+			/*Get D(u)*/
+			vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+			vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+			if(dim==3){
+				vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
+			}
+			epsxx = dvx[0];
+			epsyy = dvy[1];
+			epsxy = 0.5*(dvx[1] + dvy[0]); 
+			if(dim==3){
+				epszz = dvz[2]; 
+				epsxz = 0.5*(dvx[2] + dvz[0]);
+				epsyz = 0.5*(dvy[2] + dvz[1]);
+			}
+
+			/*Get tau^(n-1) from inputs*/
+			sigmapxx_input->GetInputValue(&sigmapxx,gauss);
+			sigmapyy_input->GetInputValue(&sigmapyy,gauss);
+			sigmapxy_input->GetInputValue(&sigmapxy,gauss);
+			if(dim==3){
+				sigmapzz_input->GetInputValue(&sigmapzz,gauss);
+				sigmapxz_input->GetInputValue(&sigmapxz,gauss);
+				sigmapyz_input->GetInputValue(&sigmapyz,gauss);
+			}
+
+			/*Get d and update tau accordingly*/
+			tau_xx[ig] = sigmapxx + r*(epsxx - d_xx[ig]);
+			tau_yy[ig] = sigmapyy + r*(epsyy - d_yy[ig]);
+			tau_xy[ig] = sigmapxy + r*(epsxy - d_xy[ig]);
+			if(dim==3){
+				tau_zz[ig] = sigmapzz + r*(epszz - d_zz[ig]);
+				tau_xz[ig] = sigmapxz + r*(epsxz - d_xz[ig]);
+				tau_yz[ig] = sigmapyz + r*(epsyz - d_yz[ig]);
+			}
+		}
+
+		/*Add inputs*/
+		element->AddInput(StrainRatexxEnum,tau_xx,P1DGEnum);
+		element->AddInput(StrainRateyyEnum,tau_yy,P1DGEnum);
+		element->AddInput(StrainRatexyEnum,tau_xy,P1DGEnum);
+		if(dim==3){
+			element->AddInput(StrainRatezzEnum,tau_zz,P1DGEnum);
+			element->AddInput(StrainRatexzEnum,tau_xz,P1DGEnum);
+			element->AddInput(StrainRateyzEnum,tau_yz,P1DGEnum);
+		}
+
 		/*Clean up and */
+		xDelete<IssmDouble>(xyz_list);
 		xDelete<IssmDouble>(tbasis);
 		xDelete<IssmDouble>(Ke);
-		xDelete<IssmDouble>(pe_xx);
-		xDelete<IssmDouble>(pe_yy);
-		xDelete<IssmDouble>(pe_zz);
-		xDelete<IssmDouble>(pe_xy);
-		xDelete<IssmDouble>(pe_xz);
-		xDelete<IssmDouble>(pe_yz);
+		xDelete<IssmDouble>(pe_xx); xDelete<IssmDouble>(d_xx);   xDelete<IssmDouble>(tau_xx);
+		xDelete<IssmDouble>(pe_yy); xDelete<IssmDouble>(d_yy);   xDelete<IssmDouble>(tau_yy);
+		xDelete<IssmDouble>(pe_zz); xDelete<IssmDouble>(d_zz);   xDelete<IssmDouble>(tau_zz);
+		xDelete<IssmDouble>(pe_xy); xDelete<IssmDouble>(d_xy);   xDelete<IssmDouble>(tau_xy);
+		xDelete<IssmDouble>(pe_xz); xDelete<IssmDouble>(d_xz);   xDelete<IssmDouble>(tau_xz);
+		xDelete<IssmDouble>(pe_yz); xDelete<IssmDouble>(d_yz);   xDelete<IssmDouble>(tau_yz);
 	}
 
+	delete gauss;
+
 }/*}}}*/
 
 /*Coupling (Tiling)*/
Index: ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp	(revision 17547)
+++ ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp	(revision 17548)
@@ -366,3 +366,16 @@
 	Ainv[8]=(A[0]*A[4]-A[1]*A[3])/det; /* = (a*e-b*d)/det */
 
 }/*}}}*/
+/*FUNCTION Matrix3x3Solve(IssmDouble* X,IssmDouble* A,IssmDouble* B) {{{*/
+void Matrix3x3Solve(IssmDouble* X,IssmDouble* A,IssmDouble* B){
+
+	IssmDouble Ainv[3][3];
+
+	Matrix3x3Invert(&Ainv[0][0],A);
+	for(int i=0;i<3;i++) X[i]=Ainv[i][0]*B[0] + Ainv[i][1]*B[1] + Ainv[i][2]*B[2];
+
+}/*}}}*/
+/*FUNCTION Matrix4x4Solve(IssmDouble* X,IssmDouble* A,IssmDouble *B) {{{*/
+void Matrix4x4Solve(IssmDouble* X,IssmDouble* A,IssmDouble *B){
+	_error_("STOP");
+}/*}}}*/
Index: ../trunk-jpl/src/c/shared/Matrix/matrix.h
===================================================================
--- ../trunk-jpl/src/c/shared/Matrix/matrix.h	(revision 17547)
+++ ../trunk-jpl/src/c/shared/Matrix/matrix.h	(revision 17548)
@@ -14,5 +14,7 @@
 void Matrix2x2Determinant(IssmDouble* Adet,IssmDouble* A);
 void Matrix3x3Invert(IssmDouble* Ainv, IssmDouble* A);
 void Matrix3x3Determinant(IssmDouble* Adet,IssmDouble* A);
+void Matrix3x3Solve(IssmDouble* X,IssmDouble* A,IssmDouble* B);
+void Matrix4x4Solve(IssmDouble* X,IssmDouble* A,IssmDouble* B);
 
 #endif //ifndef _MATRIXUTILS_H_