Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 18276)
+++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 18277)
@@ -548,5 +548,5 @@
 
 	/*Compute stress tensor: */
-	element->ComputeDeviatoricStressTensor();
+	element->ComputeStrainRate();
 
 	/*retrieve what we need: */
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 18276)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 18277)
@@ -89,4 +89,5 @@
 	int     numoutputs;
 	char**  requestedoutputs = NULL;
+	int     materials_type;
 
 	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSIAEnum));
@@ -114,4 +115,10 @@
 		parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianRlambdaEnum));
 		parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianThetaEnum));
+	}
+
+	iomodel->Constant(&materials_type,MaterialsEnum);
+	if(materials_type==MatdamageiceEnum){
+		parameters->AddObject(iomodel->CopyConstantObject(DamageC1Enum));
+		parameters->AddObject(iomodel->CopyConstantObject(DamageStressThresholdEnum));
 	}
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18276)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18277)
@@ -46,11 +46,16 @@
 	IssmDouble  eps_xx,eps_xy,eps_yy,eps_xz,eps_yz,eps_zz,eps_eff;
 	IssmDouble  epsmin=1.e-27;
-	int         dim;
-
-	/*Retrieve all inputs we will be needing: */
-	/* TODO: retrieve parameters eps_0 and eps_f and input DamageD(bar?) */
+	IssmDouble  eps_0,eps_f,sigma_0,B,D,n;
+	int         dim,counter=0;
+	IssmDouble  k1,k2,threshold=1.e-12;
+
+	/* Retrieve parameters */
 	this->GetVerticesCoordinates(&xyz_list);
 	this->ComputeStrainRate();
 	parameters->FindParam(&dim,DomainDimensionEnum);
+	parameters->FindParam(&eps_f,DamageC1Enum);
+	parameters->FindParam(&sigma_0,DamageStressThresholdEnum);
+
+	/* Retrieve inputs */
 	Input* eps_xx_input=this->GetInput(StrainRatexxEnum); _assert_(eps_xx_input);
 	Input* eps_yy_input=this->GetInput(StrainRateyyEnum); _assert_(eps_yy_input);
@@ -65,9 +70,37 @@
 	}
 
-	/* Start looping on the number of vertices: */
+	/* Fetch number of nodes and allocate output*/
+   int numnodes = this->GetNumberOfNodes();
+   IssmDouble* newD = xNew<IssmDouble>(numnodes);
+
+	/* Retrieve domain-dependent inputs */
+	Input* n_input=this->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
+   Input* damage_input = NULL;
+   Input* B_input = NULL;
+	int domaintype;
+   parameters->FindParam(&domaintype,DomainTypeEnum);
+	if(domaintype==Domain2DhorizontalEnum){
+      damage_input = this->GetInput(DamageDbarEnum);  _assert_(damage_input);
+      B_input=this->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input);
+   }
+   else{
+      damage_input = this->GetInput(DamageDEnum);   _assert_(damage_input);
+      B_input=this->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
+   }
+	/* Get B and D inputs appropriate for the type of domain */
+	/*n=this->material->GetN();*/
+   /*if(domaintype==Domain2DhorizontalEnum){
+		B=this->material->GetBbar();
+		D=this->material->GetDbar();
+   }
+	else{
+		B=this->material->GetB();
+		D=this->material->GetD();
+	}*/
+	
+	/* Start looping on the number of nodes: */
 	Gauss* gauss=this->NewGauss();
-	int numvertices = this->GetNumberOfVertices();
-	for (int iv=0;iv<numvertices;iv++){
-		gauss->GaussVertex(iv);
+	for (int i=0;i<numnodes;i++){
+		gauss->GaussNode(this->GetElementType(),i);
 
 		eps_xx_input->GetInputValue(&eps_xx,gauss);
@@ -82,14 +115,50 @@
 
 		/* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
-		eps_eff=sqrt(eps_xx*eps_xx+eps_yy*eps_yy+eps_xy*eps_xy+eps_xz*eps_xz+eps_yz*eps_yz+eps_xx*eps_xx*epsmin*epsmin);
-
-		/*TODO: compute kappa from initial D, then compute new D */
-
-	}
-
-	/* TODO: add newdamage input to DamageEnum and NewDamageEnum */
+		eps_eff=sqrt(eps_xx*eps_xx+eps_yy*eps_yy+eps_xy*eps_xy+eps_xz*eps_xz+eps_yz*eps_yz+eps_xx*eps_yy+epsmin*epsmin);
+
+		B_input->GetInputValue(&B,gauss);
+      n_input->GetInputValue(&n,gauss);
+      damage_input->GetInputValue(&D,gauss);
+
+		/* Compute kappa (k) from pre-existing level of damage, using Newton-Raphson */
+		k1=exp(n*eps_0/(eps_f-eps_0)+n*log(1.-D)-log(eps_0)); /* initial guess */
+		
+		while(true){
+      	/*Newton step k2=k1-f(k1)/f'(k1) */
+      	k2=k1-(1.-D-pow(eps_0/k1,1./n)*exp(-(k1-eps_0)/(eps_f-eps_0)))/(pow(eps_0/k1,1./n)*exp(-(k1-eps_0)/(eps_f-eps_0))*(-1./(eps_f-eps_0)-1./n/k1));
+
+      	if( fabs(k2-k1)/(fabs(k2))<threshold ){
+         	break;
+      	}
+      	else{
+      	   k1=k2;
+      	   counter++;
+      	}
+
+      	if(counter>50) break;
+   	}
+
+		/* Compute threshold strain rate from threshold stress */
+		eps_0=pow(sigma_0/B,n); 
+		_assert_(eps_f>eps_0);
+
+		if(eps_eff>k2){
+			newD[i]=1.-pow(eps_0/k2,1./n)*exp(-(k2-eps_0)/(eps_f-eps_0));
+		}
+		else newD[i]=D;
+	}
+
+	/* Add new damage input to DamageEnum and NewDamageEnum */
+	this->AddInput(NewDamageEnum,newD,this->GetElementType());
+	if(domaintype==Domain2DhorizontalEnum){
+		this->AddInput(DamageDbarEnum,newD,this->GetElementType());
+	}
+	else{
+		this->AddInput(DamageDEnum,newD,this->GetElementType());
+	}
 
 	/*Clean up and return*/
 	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(newD);
 	delete gauss;