Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 18135)
+++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 18136)
@@ -266,5 +266,5 @@
 
 	/*Intermediaries*/
-	int      domaintype;
+	int      domaintype,damagelaw;
 	Element* basalelement;
 	IssmDouble  Jdet,dt;
@@ -285,4 +285,5 @@
 	}
 
+
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = basalelement->GetNumberOfNodes();
@@ -295,5 +296,11 @@
 	basalelement->GetVerticesCoordinates(&xyz_list);
 	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
-	this->CreateDamageFInput(basalelement);
+	basalelement->FindParam(&damagelaw,DamageLawEnum);
+	if(damagelaw==1 | damagelaw==2){
+		this->CreateDamageFInputPralong(basalelement);
+	}
+	else if(damagelaw==3){
+		this->CreateDamageFInputExp(basalelement);
+	}
 	Input* damaged_input = basalelement->GetInput(DamageDEnum); _assert_(damaged_input);
 	Input* damagef_input = basalelement->GetInput(DamageFEnum); _assert_(damagef_input);
@@ -440,5 +447,5 @@
 
 /*Intermediaries*/
-void DamageEvolutionAnalysis::CreateDamageFInput(Element* element){/*{{{*/
+void DamageEvolutionAnalysis::CreateDamageFInputPralong(Element* element){/*{{{*/
 
 	/*Intermediaries */
@@ -527,2 +534,73 @@
 	delete gauss;
 }/*}}}*/
+void DamageEvolutionAnalysis::CreateDamageFInputExp(Element* element){/*{{{*/
+
+	/*Intermediaries */
+	IssmDouble epsf,stress_threshold,eps0;
+	IssmDouble damage,B,n,epseff;
+	IssmDouble eps_xx,eps_yy,eps_xy,eps1,eps2,epstmp;
+	int domaintype,damagelaw;
+
+	/*Fetch number of vertices and allocate output*/
+	int numnodes = element->GetNumberOfNodes();
+	IssmDouble* f   = xNew<IssmDouble>(numnodes);
+
+	/*retrieve parameters:*/
+	element->FindParam(&epsf,DamageC1Enum);
+	element->FindParam(&stress_threshold,DamageStressThresholdEnum);
+	element->FindParam(&domaintype,DomainTypeEnum);
+	element->FindParam(&damagelaw,DamageLawEnum);
+
+	/*Compute stress tensor: */
+	element->ComputeDeviatoricStressTensor();
+
+	/*retrieve what we need: */
+	Input* eps_xx_input  = element->GetInput(StrainRatexxEnum);     _assert_(eps_xx_input);
+	Input* eps_xy_input  = element->GetInput(StrainRatexyEnum);     _assert_(eps_xy_input);
+	Input* eps_yy_input  = element->GetInput(StrainRateyyEnum);     _assert_(eps_yy_input);
+	Input*  n_input=element->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
+	Input* damage_input = NULL;
+	Input* B_input = NULL;
+	if(domaintype==Domain2DhorizontalEnum){
+		damage_input = element->GetInput(DamageDbarEnum); 	_assert_(damage_input);
+		B_input=element->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input);
+	}
+	else{
+		damage_input = element->GetInput(DamageDEnum);   _assert_(damage_input);
+		B_input=element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
+	}
+
+	/*Calculate damage evolution source term: */
+	Gauss* gauss=element->NewGauss();
+	for (int i=0;i<numnodes;i++){
+		gauss->GaussNode(element->GetElementType(),i);
+		
+		eps_xx_input->GetInputValue(&eps_xx,gauss);
+		eps_xy_input->GetInputValue(&eps_xy,gauss);
+		eps_yy_input->GetInputValue(&eps_yy,gauss);
+		B_input->GetInputValue(&B,gauss);
+		n_input->GetInputValue(&n,gauss);
+		damage_input->GetInputValue(&damage,gauss);
+	
+		/*Calculate principal effective strain rates*/
+		eps1=(eps_xx+eps_yy)/2.+sqrt(pow((eps_xx-eps_yy)/2.,2)+pow(eps_xy,2));
+		eps2=(eps_xx+eps_yy)/2.-sqrt(pow((eps_xx-eps_yy)/2.,2)+pow(eps_xy,2));
+		if(fabs(eps2)>fabs(eps1)){epstmp=eps2; eps2=eps1; eps1=epstmp;}
+
+		/*Calculate effective strain rate and threshold strain rate*/
+		epseff=1./sqrt(2.)*sqrt(eps1*eps1-eps1*eps2+eps2*eps2);
+		eps0=pow(stress_threshold/B,n);
+
+		if(epseff>eps0){
+			f[i]=1.-pow(eps0/epseff,1./n)*exp(-(epseff-eps0)/(epsf-eps0))-damage;
+		}
+		else f[i]=0;
+	}
+
+	/*Add input*/
+	element->AddInput(DamageFEnum,f,element->GetElementType());
+	
+	/*Clean up and return*/
+	xDelete<IssmDouble>(f);
+	delete gauss;
+}/*}}}*/
Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h	(revision 18135)
+++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h	(revision 18136)
@@ -34,5 +34,6 @@
 
 		/*Intermediaries*/
-		void CreateDamageFInput(Element* element);
+		void CreateDamageFInputPralong(Element* element);
+		void CreateDamageFInputExp(Element* element);
 };
 #endif