Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18281)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18282)
@@ -88,14 +88,4 @@
       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: */
@@ -121,10 +111,21 @@
       damage_input->GetInputValue(&D,gauss);
 
+		/* Compute threshold strain rate from threshold stress */
+		eps_0=pow(sigma_0/B,n); 
+		_assert_(eps_f>eps_0);
+
 		/* 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 */
-		
+		/* provide a reasonable initial guess */
+		if(D==0){
+			k1=eps_0;
+		}
+		else{
+			k1=exp(n*eps_0/(eps_f-eps_0)-n*log(1.-D)+log(eps_0)); /* initial guess */
+		}
+	
+		counter=0;
 		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));
+			k2=k1-(k1+(eps_f-eps_0)/n*log(k1)-eps_0+(eps_f-eps_0)*(log(1.-D)-1./n*log(eps_0)))/(1.+(eps_f-eps_0)/n/k1);
 
       	if( fabs(k2-k1)/(fabs(k2))<threshold ){
@@ -132,5 +133,5 @@
       	}
       	else{
-      	   k1=k2;
+				k1=k2;
       	   counter++;
       	}
@@ -139,10 +140,6 @@
    	}
 
-		/* 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));
+			newD[i]=1.-pow(eps_0/eps_eff,1./n)*exp(-(eps_eff-eps_0)/(eps_f-eps_0));
 		}
 		else newD[i]=D;