Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 18040)
+++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 18041)
@@ -32,7 +32,7 @@
 
 	/*Retrieve law dependent parameters: */
-	char* law  = NULL;
+	int law;
 	iomodel->Constant(&law,DamageLawEnum);
-	if (strcmp(law,"pralong")==0){
+	if (law>0){
 		parameters->AddObject(iomodel->CopyConstantObject(DamageC1Enum));
 		parameters->AddObject(iomodel->CopyConstantObject(DamageC2Enum));
@@ -43,5 +43,4 @@
 		parameters->AddObject(iomodel->CopyConstantObject(DamageEquivStressEnum));
 	}
-	xDelete<char>(law);
 
 }/*}}}*/
@@ -445,5 +444,5 @@
 	IssmDouble J2s,Chi,Psi,PosPsi,NegPsi;
 	IssmDouble damage,tau_xx,tau_xy,tau_yy;
-	int equivstress,domaintype;
+	int equivstress,domaintype,damagelaw;
 
 	/*Fetch number of vertices and allocate output*/
@@ -458,4 +457,5 @@
 	element->FindParam(&stress_threshold,DamageStressThresholdEnum);
 	element->FindParam(&domaintype,DomainTypeEnum);
+	element->FindParam(&damagelaw,DamageLawEnum);
 
 	/*Compute stress tensor: */
@@ -504,8 +504,15 @@
 		}
 		Psi=Chi-stress_threshold;
-		PosPsi=max(Psi,0.);
 		NegPsi=max(-Chi,0.); /* healing only for compressive stresses */
 
-		f[i]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1.-damage),-c3);
+		if(damagelaw==1){
+			PosPsi=max(Psi,0.);
+			f[i]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
+		}
+		else if(damagelaw==2){
+			PosPsi=max(Psi,1.);
+			f[i]= c1*(pow(log10(PosPsi),c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
+		}
+		else _error_("damage law not supported");
 	}
 
