Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 19387)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 19388)
@@ -367,4 +367,5 @@
 if DAMAGEEVOLUTION
 issm_sources += ./analyses/DamageEvolutionAnalysis.cpp
+issm_sources += ./modules/Damagex/Damagex.cpp
 endif
 if STRESSBALANCE
Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 19387)
+++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 19388)
@@ -130,5 +130,5 @@
 	IssmDouble damage,B,n,epseff;
 	IssmDouble eps_xx,eps_yy,eps_xy,eps1,eps2,epstmp;
-	int domaintype,damagelaw;
+	int domaintype;
 
 	/*Fetch number of vertices and allocate output*/
@@ -140,5 +140,4 @@
 	element->FindParam(&stress_threshold,DamageStressThresholdEnum);
 	element->FindParam(&domaintype,DomainTypeEnum);
-	element->FindParam(&damagelaw,DamageLawEnum);
 
 	/*Compute stress tensor: */
@@ -202,5 +201,5 @@
 	IssmDouble J2s,Chi,Psi,PosPsi,NegPsi;
 	IssmDouble damage,tau_xx,tau_xy,tau_xz,tau_yy,tau_yz,tau_zz,stressMaxPrincipal;
-	int equivstress,domaintype,damagelaw,dim;
+	int equivstress,domaintype,dim;
 
 	/*Fetch number of vertices and allocate output*/
@@ -215,5 +214,4 @@
 	element->FindParam(&stress_threshold,DamageStressThresholdEnum);
 	element->FindParam(&domaintype,DomainTypeEnum);
-	element->FindParam(&damagelaw,DamageLawEnum);
 
 	/*Get problem dimension*/
@@ -291,14 +289,6 @@
 			Psi=Chi-stress_threshold;
 			NegPsi=max(-Chi,0.); /* healing only for compressive stresses */
-
-			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");
+			PosPsi=max(Psi,0.);
+			f[i]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
 		}
 		else{
@@ -312,13 +302,6 @@
 			Psi=Chi-stress_threshold;
 			NegPsi=max(-Chi,0.); /* healing only for compressive stresses */
-			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");
+			PosPsi=max(Psi,0.);
+			f[i]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
 		}
 	}
@@ -534,7 +517,4 @@
 			break;
 		case 2:
-			this->CreateDamageFInputPralong(element);
-			break;
-		case 3:
 			this->CreateDamageFInputExp(element);
 			break;
Index: /issm/trunk-jpl/src/c/cores/damage_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/damage_core.cpp	(revision 19387)
+++ /issm/trunk-jpl/src/c/cores/damage_core.cpp	(revision 19388)
@@ -28,4 +28,5 @@
 
 	if(VerboseSolution()) _printf0_("   computing damage\n");
+	Damagex(femmodel); /* optionally calculate damage analytically first */
 	femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum);
 	if(stabilization==4){
Index: /issm/trunk-jpl/src/c/modules/Damagex/Damagex.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Damagex/Damagex.cpp	(revision 19388)
+++ /issm/trunk-jpl/src/c/modules/Damagex/Damagex.cpp	(revision 19388)
@@ -0,0 +1,29 @@
+/*!\file Damagex
+ * \brief: compute damage
+ */
+
+#include "./Damagex.h"
+#include "../../shared/shared.h"
+#include "../../toolkits/toolkits.h"
+
+void Damagex(FemModel* femmodel){
+
+	/*Recover Damage law Enum*/
+	int damagelaw;
+	femmodel->parameters->FindParam(&damagelaw,DamageLawEnum);
+
+	/*Calculate damage*/
+	switch(damagelaw){
+		case 0:
+			if(VerboseModule()) _printf0_("   computing damage analytically\n");
+			femmodel->ElementOperationx(&Element::ComputeNewDamage);
+			break;
+		case 1:
+		case 2:
+			if(VerboseModule()) _printf0_("   computing damage using source term in advection scheme\n");
+			/* Damage calculated using source term in DamageEvolutionAnalysis */
+			break;
+		default:
+			_error_("Damage law "<<EnumToStringx(damagelaw)<<" not implemented yet");
+	}
+}
Index: /issm/trunk-jpl/src/c/modules/Damagex/Damagex.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/Damagex/Damagex.h	(revision 19388)
+++ /issm/trunk-jpl/src/c/modules/Damagex/Damagex.h	(revision 19388)
@@ -0,0 +1,10 @@
+#ifndef _DAMAGEX_H
+#define _DAMGEX_H
+
+#include "../../classes/classes.h"
+#include "../../analyses/analyses.h"
+
+/* local prototypes: */
+void Damagex(FemModel* femmodel);
+
+#endif
Index: /issm/trunk-jpl/src/c/modules/modules.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/modules.h	(revision 19387)
+++ /issm/trunk-jpl/src/c/modules/modules.h	(revision 19388)
@@ -20,4 +20,5 @@
 #include "./CreateNodalConstraintsx/CreateNodalConstraintsx.h"
 #include "./CreateJacobianMatrixx/CreateJacobianMatrixx.h"
+#include "./Damagex/Damagex.h"
 #include "./DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h"
 #include "./ExpToLevelSetx/ExpToLevelSetx.h"