Index: ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 19387)
+++ ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 19388)
@@ -129,7 +129,7 @@
 	IssmDouble epsf,stress_threshold,eps0;
 	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*/
 	int numnodes = element->GetNumberOfNodes();
@@ -139,7 +139,6 @@
 	element->FindParam(&epsf,DamageC1Enum);
 	element->FindParam(&stress_threshold,DamageStressThresholdEnum);
 	element->FindParam(&domaintype,DomainTypeEnum);
-	element->FindParam(&damagelaw,DamageLawEnum);
 
 	/*Compute stress tensor: */
 	element->ComputeStrainRate();
@@ -201,7 +200,7 @@
 	IssmDouble s_xx,s_xy,s_xz,s_yy,s_yz,s_zz,s1,s2,s3,stmp;
 	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*/
 	int numnodes = element->GetNumberOfNodes();
@@ -214,7 +213,6 @@
 	element->FindParam(&healing,DamageHealingEnum);
 	element->FindParam(&stress_threshold,DamageStressThresholdEnum);
 	element->FindParam(&domaintype,DomainTypeEnum);
-	element->FindParam(&damagelaw,DamageLawEnum);
 
 	/*Get problem dimension*/
 	switch(domaintype){
@@ -290,16 +288,8 @@
 			}
 			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{
 			if(equivstress==1){/* max principal stress */
@@ -311,15 +301,8 @@
 			}
 			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);
 		}
 	}
 	/*Add input*/
@@ -533,9 +516,6 @@
 			this->CreateDamageFInputPralong(element);
 			break;
 		case 2:
-			this->CreateDamageFInputPralong(element);
-			break;
-		case 3:
 			this->CreateDamageFInputExp(element);
 			break;
 		default:
Index: ../trunk-jpl/src/c/modules/Damagex/Damagex.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/Damagex/Damagex.cpp	(revision 0)
+++ ../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: ../trunk-jpl/src/c/modules/Damagex/Damagex.h
===================================================================
--- ../trunk-jpl/src/c/modules/Damagex/Damagex.h	(revision 0)
+++ ../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

Property changes on: ../trunk-jpl/src/c/modules/Damagex
___________________________________________________________________
Added: svn:ignore
   + *.deps
*.dirstamp


Index: ../trunk-jpl/src/c/modules/modules.h
===================================================================
--- ../trunk-jpl/src/c/modules/modules.h	(revision 19387)
+++ ../trunk-jpl/src/c/modules/modules.h	(revision 19388)
@@ -19,6 +19,7 @@
 #include "./ControlInputSetGradientx/ControlInputSetGradientx.h"
 #include "./CreateNodalConstraintsx/CreateNodalConstraintsx.h"
 #include "./CreateJacobianMatrixx/CreateJacobianMatrixx.h"
+#include "./Damagex/Damagex.h"
 #include "./DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h"
 #include "./ExpToLevelSetx/ExpToLevelSetx.h"
 #include "./ElementConnectivityx/ElementConnectivityx.h"
Index: ../trunk-jpl/src/c/Makefile.am
===================================================================
--- ../trunk-jpl/src/c/Makefile.am	(revision 19387)
+++ ../trunk-jpl/src/c/Makefile.am	(revision 19388)
@@ -366,6 +366,7 @@
 endif
 if DAMAGEEVOLUTION
 issm_sources += ./analyses/DamageEvolutionAnalysis.cpp
+issm_sources += ./modules/Damagex/Damagex.cpp
 endif
 if STRESSBALANCE
 issm_sources += ./analyses/StressbalanceAnalysis.cpp
Index: ../trunk-jpl/src/c/cores/damage_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/damage_core.cpp	(revision 19387)
+++ ../trunk-jpl/src/c/cores/damage_core.cpp	(revision 19388)
@@ -27,6 +27,7 @@
 	femmodel->parameters->FindParam(&stabilization,DamageStabilizationEnum);
 
 	if(VerboseSolution()) _printf0_("   computing damage\n");
+	Damagex(femmodel); /* optionally calculate damage analytically first */
 	femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum);
 	if(stabilization==4){
 		solutionsequence_fct(femmodel);