Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 17743)
+++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 17744)
@@ -47,4 +47,6 @@
 void DamageEvolutionAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
 
+	int finiteelement = P1Enum;
+
 	/*Update elements: */
 	int counter=0;
@@ -52,5 +54,5 @@
 		if(iomodel->my_elements[i]){
 			Element* element=(Element*)elements->GetObjectByOffset(counter);
-			element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
+			element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement);
 			counter++;
 		}
@@ -73,14 +75,17 @@
 void DamageEvolutionAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
 
+	int finiteelement = P1Enum;
+
 	if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(1,MeshVertexonbaseEnum);
-	::CreateNodes(nodes,iomodel,DamageEvolutionAnalysisEnum,P1Enum);
+	::CreateNodes(nodes,iomodel,DamageEvolutionAnalysisEnum,finiteelement);
 	iomodel->DeleteData(1,MeshVertexonbaseEnum);
 }/*}}}*/
 void DamageEvolutionAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
 
-	int stabilization,finitelelement;
+	int stabilization;
+	int finiteelement = P1Enum;
 	iomodel->Constant(&stabilization,DamageStabilizationEnum);
 
-	IoModelToConstraintsx(constraints,iomodel,DamageSpcdamageEnum,DamageEvolutionAnalysisEnum,P1Enum);
+	IoModelToConstraintsx(constraints,iomodel,DamageSpcdamageEnum,DamageEvolutionAnalysisEnum,finiteelement);
 
 }/*}}}*/
@@ -411,5 +416,10 @@
 
 	/*Get all inputs and parameters*/
-	basalelement->AddBasalInput(DamageDEnum,newdamage,P1Enum);
+	if(domaintype==Domain2DhorizontalEnum){
+		basalelement->AddInput(DamageDbarEnum,newdamage,element->GetElementType());
+	}
+	else{
+		basalelement->AddBasalInput(DamageDEnum,newdamage,element->GetElementType());
+	}
 
 	/*Free ressources:*/
@@ -431,5 +441,5 @@
 	IssmDouble J2s,Chi,Psi,PosPsi,NegPsi;
 	IssmDouble damage,tau_xx,tau_xy,tau_yy;
-	int equivstress;
+	int equivstress,domaintype;
 
 	/*Fetch number of vertices and allocate output*/
@@ -443,4 +453,5 @@
 	element->FindParam(&healing,DamageHealingEnum);
 	element->FindParam(&stress_threshold,DamageStressThresholdEnum);
+	element->FindParam(&domaintype,DomainTypeEnum);
 
 	/*Compute stress tensor: */
@@ -451,5 +462,13 @@
 	Input* tau_xy_input  = element->GetInput(DeviatoricStressxyEnum);     _assert_(tau_xy_input);
 	Input* tau_yy_input  = element->GetInput(DeviatoricStressyyEnum);     _assert_(tau_yy_input);
-	Input* damage_input    = element->GetInput(DamageDEnum); _assert_(damage_input);
+	Input* damage_input = NULL;
+	if(domaintype==Domain2DhorizontalEnum){
+		damage_input = element->GetInput(DamageDbarEnum); 	_assert_(damage_input);
+	}
+	else{
+		damage_input = element->GetInput(DamageDEnum);   _assert_(damage_input);
+	}
+
+
 
 	/*retrieve the desired type of equivalent stress*/