Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 16187)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 16188)
@@ -259,4 +259,5 @@
 					./modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp\
 					./modules/ModelProcessorx/CreateNodes.cpp\
+					./modules/ConstraintsStatex/PengridIsPresent.cpp\
 					./modules/ParseToolkitsOptionsx/ParseToolkitsOptionsx.h\
 					./modules/ParseToolkitsOptionsx/ParseToolkitsOptionsx.cpp\
@@ -409,5 +410,4 @@
 					   ./modules/ModelProcessorx/Melting/CreateLoadsMelting.cpp\
 					   ./modules/ConstraintsStatex/ThermalConstraintsState.cpp\
-					   ./modules/ConstraintsStatex/ThermalIsPresent.cpp\
 					   ./modules/ResetConstraintsx/ThermalConstraintsReset.cpp \
 					   ./analyses/thermal_core.cpp\
@@ -551,6 +551,9 @@
 					./modules/ModelProcessorx/Damage/CreateNodesDamage.cpp \
 					./modules/ModelProcessorx/Damage/CreateConstraintsDamage.cpp\
+					./modules/ConstraintsStatex/DamageConstraintsState.cpp\
+					./modules/ResetConstraintsx/DamageConstraintsReset.cpp \
 					./modules/ModelProcessorx/Damage/CreateParametersDamage.cpp\
-					./modules/ModelProcessorx/Damage/CreateLoadsDamage.cpp
+					./modules/ModelProcessorx/Damage/CreateLoadsDamage.cpp\
+					./solutionsequences/solutionsequence_damage_nonlinear.cpp
 
 #}}}
Index: /issm/trunk-jpl/src/c/analyses/damage_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/damage_core.cpp	(revision 16187)
+++ /issm/trunk-jpl/src/c/analyses/damage_core.cpp	(revision 16188)
@@ -1,4 +1,4 @@
-/*!\file: damage_core.cpp
- * \brief: core of the damage solution 
+/* 
+ * \brief: damgage_core.cpp: core for the damage solution
  */ 
 
@@ -12,63 +12,27 @@
 void damage_core(FemModel* femmodel){
 
-	/*parameters: */
-	bool  save_results;
-	bool  issmbgradients,ispdd,isdelta18o,isFS,isfreesurface;
-	int   solution_type;
-	int  *requested_outputs = NULL;
-	int   numoutputs        = 0;
+	/*intermediary*/
+	bool   save_results;
+	bool   dakota_analysis  = false;
+	int    solution_type;
 
-	/*activate configuration*/
-	femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum);
-	  
-	if(VerboseSolution())_printf_("	call damage core\n\n");
+	if(VerboseSolution()) _printf0_("   computing damage\n");
+	
+	//first recover parameters common to all solutions
+	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+	femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
+	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
 
-	/*recover parameters: */
-	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
-	femmodel->parameters->FindParam(&issmbgradients,SurfaceforcingsIssmbgradientsEnum);
-	femmodel->parameters->FindParam(&ispdd,SurfaceforcingsIspddEnum);
-	femmodel->parameters->FindParam(&isdelta18o,SurfaceforcingsIsdelta18oEnum);
-	femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
-	femmodel->parameters->FindParam(&isfreesurface,MasstransportIsfreesurfaceEnum);
-	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
-	femmodel->parameters->FindParam(&numoutputs,MasstransportNumRequestedOutputsEnum);
-	if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,MasstransportRequestedOutputsEnum);
-
-	if(issmbgradients){
-	  if(VerboseSolution())_printf_("	call smb gradients module\n\n");
-	  SmbGradientsx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
-	}
-	if(ispdd){
-		if(isdelta18o){
-			if(VerboseSolution()) _printf0_("   call Delta18oParametrization module\n");
-			Delta18oParameterizationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
-		} 
-		if(VerboseSolution()) _printf0_("   call positive degree day module\n");
-		PositiveDegreeDayx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+	if(dakota_analysis && solution_type!=TransientSolutionEnum){
+		femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum);
+		ResetConstraintsx(femmodel);
 	}
 
-	if(isFS && isfreesurface){
-		if(VerboseSolution()) _printf0_("   call free surface computational core\n");
-		femmodel->SetCurrentConfiguration(FreeSurfaceBaseAnalysisEnum);
-		solutionsequence_linear(femmodel);
-		femmodel->SetCurrentConfiguration(FreeSurfaceTopAnalysisEnum);
-		solutionsequence_linear(femmodel);
-	}
-	else{
-		if(VerboseSolution()) _printf0_("   call computational core\n");
-		solutionsequence_linear(femmodel);
-	}
+	femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum);
+	solutionsequence_damage_nonlinear(femmodel);
 
 	if(save_results){
 		if(VerboseSolution()) _printf0_("   saving results\n");
-		InputToResultx(femmodel,ThicknessEnum);
-		InputToResultx(femmodel,BedEnum);
-		InputToResultx(femmodel,SurfaceEnum);
-		femmodel->RequestedOutputsx(requested_outputs,numoutputs);
+		InputToResultx(femmodel,DamageDEnum);
 	}
-
-	if(solution_type==MasstransportSolutionEnum)femmodel->RequestedDependentsx();
-
-	/*Free ressources:*/
-	xDelete<int>(requested_outputs);
 }
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16187)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16188)
@@ -54,4 +54,5 @@
 		virtual void   GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue)=0;
 		virtual void   GetInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0;
+		virtual void   GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0;
 
 		virtual IssmDouble SurfaceArea(void)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16187)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16188)
@@ -1422,4 +1422,17 @@
 
 	Input* input=inputs->GetInput(enumtype);
+	if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
+
+	GaussPenta* gauss=new GaussPenta();
+	gauss->GaussVertex(this->GetNodeIndex(node));
+
+	input->GetInputValue(pvalue,gauss);
+	delete gauss;
+}
+/*}}}*/
+/*FUNCTION Penta::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype) {{{*/
+void Penta::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){
+
+	Input* input=this->material->inputs->GetInput(enumtype);
 	if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16187)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16188)
@@ -200,4 +200,5 @@
 		void           GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
 		void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
+		void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
 		void	         GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);
 		void           GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]);
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16187)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16188)
@@ -176,4 +176,5 @@
 void  Tria::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){
 
+
 	/*Skip if water element*/
 	if(NoIceInElement()) return;
@@ -223,5 +224,10 @@
 			return CreateKMatrixStressbalanceSIA();
 			break;
-		 #endif
+		#endif
+	 	#ifdef _HAVE_DAMAGE_
+		case DamageEvolutionAnalysisEnum:
+			return CreateKMatrixDamageEvolution();
+			break;
+		#endif
 		case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
 			return CreateMassMatrix();
@@ -363,5 +369,5 @@
 	/*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
 	switch(analysis_type){
-#ifdef _HAVE_STRESSBALANCE_
+		#ifdef _HAVE_STRESSBALANCE_
 		case StressbalanceAnalysisEnum:
 			return CreatePVectorStressbalanceSSA();
@@ -370,5 +376,10 @@
 			return CreatePVectorStressbalanceSIA();
 			break;
-#endif
+		#endif
+	 	#ifdef _HAVE_DAMAGE_
+		case DamageEvolutionAnalysisEnum:
+			return CreatePVectorDamageEvolution();
+			break;
+		#endif
 		case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
 			return CreatePVectorSlope();
@@ -377,5 +388,5 @@
 			return CreatePVectorMasstransport();
 			break;
-#ifdef _HAVE_HYDROLOGY_
+		#ifdef _HAVE_HYDROLOGY_
 		case HydrologyShreveAnalysisEnum:
 			return CreatePVectorHydrologyShreve();
@@ -387,6 +398,6 @@
 			return CreatePVectorHydrologyDCEfficient();
 			break;
-#endif
-#ifdef _HAVE_BALANCED_
+		#endif
+		#ifdef _HAVE_BALANCED_
 		case BalancethicknessAnalysisEnum:
 			return CreatePVectorBalancethickness();
@@ -401,6 +412,6 @@
 			return CreatePVectorSmoothedSlopeY();
 			break;
-#endif
-#ifdef _HAVE_CONTROL_
+		#endif
+		#ifdef _HAVE_CONTROL_
 		case AdjointBalancethicknessAnalysisEnum:
 			return CreatePVectorAdjointBalancethickness();
@@ -409,5 +420,5 @@
 			return CreatePVectorAdjointHoriz();
 			break;
-#endif
+	 	#endif
 		default:
 			_error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
@@ -1239,4 +1250,17 @@
 
 	Input* input=inputs->GetInput(enumtype);
+	if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
+
+	GaussTria* gauss=new GaussTria();
+	gauss->GaussVertex(this->GetNodeIndex(node));
+
+	input->GetInputValue(pvalue,gauss);
+	delete gauss;
+}
+/*}}}*/
+/*FUNCTION Tria::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype) {{{*/
+void Tria::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){
+
+	Input* input=this->material->inputs->GetInput(enumtype);
 	if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
 
@@ -1623,4 +1647,9 @@
 		case HydrologyDCEfficientAnalysisEnum:
 			InputUpdateFromSolutionOneDof(solution,EplHeadEnum);
+			break;
+		#endif
+	 	#ifdef _HAVE_DAMAGE_
+		case DamageEvolutionAnalysisEnum:
+			InputUpdateFromSolutionDamageEvolution(solution);
 			break;
 		#endif
@@ -6825,4 +6854,297 @@
 #endif
 
+#ifdef _HAVE_DAMAGE_
+/*FUNCTION Tria::CreateKMatrixDamageEvolution {{{*/
+ElementMatrix* Tria::CreateKMatrixDamageEvolution(void){
+
+	switch(GetElementType()){
+		case P1Enum: case P2Enum:
+			return CreateKMatrixDamageEvolution_CG();
+		case P1DGEnum:
+			_error_("DG not implemented yet!");break;
+		default:
+			_error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");
+	}
+}
+/*}}}*/
+/*FUNCTION Tria::CreateKMatrixDamageEvolution_CG {{{*/
+ElementMatrix* Tria::CreateKMatrixDamageEvolution_CG(void){
+
+	/*Intermediaries */
+	int        stabilization;
+	int        dim;
+	IssmDouble Jdet,D_scalar,dt,h;
+	IssmDouble vel,vx,vy,dvxdx,dvydy;
+	IssmDouble dvx[2],dvy[2];
+	IssmDouble xyz_list[NUMVERTICES][3];
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Initialize Element matrix and vectors*/
+	ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+	IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
+	IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
+	IssmDouble     D[2][2];
+
+	/*Retrieve all inputs and parameters*/
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+	this->parameters->FindParam(&dim,MeshDimensionEnum);
+	this->parameters->FindParam(&stabilization,DamageStabilizationEnum);
+	Input* vxaverage_input=NULL;
+	Input* vyaverage_input=NULL;
+	if(dim==2){
+		vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
+		vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
+	}
+	else{
+		vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
+		vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
+	}
+	h=sqrt(2*this->GetArea());
+
+	/* Start  looping on the number of gaussian points: */
+	GaussTria *gauss=new GaussTria(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+		GetNodalFunctions(basis,gauss);
+
+		vxaverage_input->GetInputValue(&vx,gauss);
+		vyaverage_input->GetInputValue(&vy,gauss);
+		vxaverage_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
+		vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
+
+		D_scalar=gauss->weight*Jdet;
+
+		TripleMultiply(basis,1,numnodes,1,
+					&D_scalar,1,1,0,
+					basis,1,numnodes,0,
+					&Ke->values[0],1);
+		GetBMasstransport(B,&xyz_list[0][0],gauss);
+		GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
+
+		dvxdx=dvx[0];
+		dvydy=dvy[1];
+		D_scalar=dt*gauss->weight*Jdet;
+
+		D[0][0]=D_scalar*dvxdx;
+		D[0][1]=0.;
+		D[1][0]=0.;
+		D[1][1]=D_scalar*dvydy;
+		TripleMultiply(B,2,numnodes,1,
+					&D[0][0],2,2,0,
+					B,2,numnodes,0,
+					&Ke->values[0],1);
+
+		D[0][0]=D_scalar*vx;
+		D[1][1]=D_scalar*vy;
+		TripleMultiply(B,2,numnodes,1,
+					&D[0][0],2,2,0,
+					Bprime,2,numnodes,0,
+					&Ke->values[0],1);
+
+		if(stabilization==2){
+			/*Streamline upwinding*/
+			vel=sqrt(vx*vx+vy*vy)+1.e-8;
+			D[0][0]=h/(2*vel)*vx*vx;
+			D[1][0]=h/(2*vel)*vy*vx;
+			D[0][1]=h/(2*vel)*vx*vy;
+			D[1][1]=h/(2*vel)*vy*vy;
+		}
+		else if(stabilization==1){
+			/*SSA*/
+			vxaverage_input->GetInputAverage(&vx);
+			vyaverage_input->GetInputAverage(&vy);
+			D[0][0]=h/2.0*fabs(vx);
+			D[0][1]=0.;
+			D[1][0]=0.;
+			D[1][1]=h/2.0*fabs(vy);
+		}
+		if(stabilization==1 || stabilization==2){
+			D[0][0]=D_scalar*D[0][0];
+			D[1][0]=D_scalar*D[1][0];
+			D[0][1]=D_scalar*D[0][1];
+			D[1][1]=D_scalar*D[1][1];
+			TripleMultiply(Bprime,2,numnodes,1,
+						&D[0][0],2,2,0,
+						Bprime,2,numnodes,0,
+						&Ke->values[0],1);
+		}
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(Bprime);
+	delete gauss;
+	return Ke;
+}
+/*}}}*/
+/*FUNCTION Tria::CreatePVectorDamageEvolution{{{*/
+ElementVector* Tria::CreatePVectorDamageEvolution(void){
+
+	switch(GetElementType()){
+		case P1Enum: case P2Enum:
+			return CreatePVectorDamageEvolution_CG();
+		case P1DGEnum:
+			_error_("DG not implemented yet");
+		default:
+			_error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");
+	}
+}
+/*}}}*/
+/*FUNCTION Tria::CreatePVectorDamageEvolution_CG {{{*/
+ElementVector* Tria::CreatePVectorDamageEvolution_CG(void){
+
+	/*Intermediaries */
+	IssmDouble Jdet,dt;
+	IssmDouble f,damage;
+	IssmDouble xyz_list[NUMVERTICES][3];
+	Input* damage_input  = NULL;
+	Input* f_input  = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+	damage_input  = this->material->inputs->GetInput(DamageDbarEnum);     _assert_(damage_input);
+	
+	/*retrieve damage evolution forcing function: */
+	f_input=this->DamageEvolutionF();
+
+	/*Initialize forcing function f to 0, do not forget!:*/
+	/* Start  looping on the number of gaussian points: */
+	GaussTria* gauss=new GaussTria(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+		GetNodalFunctions(basis,gauss);
+
+		f_input->GetInputValue(&f,gauss);
+		damage_input->GetInputValue(&damage,gauss);
+
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(damage+dt*f)*basis[i];
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	return pe;
+}
+/*}}}*/
+/*FUNCTION Tria::DamageEvolutionF{{{*/
+Input* Tria::DamageEvolutionF(void){
+
+	/*Intermediaries */
+	IssmDouble c1,c2,c3,healing,stress_threshold;
+	IssmDouble s_xx,s_xy,s_yy;
+	IssmDouble J2s;
+	IssmDouble Xis;
+	IssmDouble Psi;
+	IssmDouble PosPsi;
+	IssmDouble NegPsi;
+	Input* damage_input=NULL;
+	Input* z_input=NULL;
+	Input* sigma_xx_input  = NULL;
+	Input* sigma_xy_input  = NULL;
+	Input* sigma_yy_input  = NULL;
+	GaussTria* gauss=NULL;
+	IssmDouble damage,sigma_xx,sigma_xy,sigma_yy;
+	
+	/*output: */
+	IssmDouble f[NUMVERTICES];
+	Input* f_input=NULL;
+
+	/*retrieve parameters:*/
+	this->parameters->FindParam(&c1,DamageC1Enum);
+	this->parameters->FindParam(&c2,DamageC2Enum);
+	this->parameters->FindParam(&c3,DamageC3Enum);
+	this->parameters->FindParam(&healing,DamageHealingEnum);
+	this->parameters->FindParam(&stress_threshold,DamageStressThresholdEnum);
+
+	/*Compute stress tensor: */
+	this->ComputeStressTensor();
+
+	/*retrieve what we need: */
+	sigma_xx_input  = inputs->GetInput(StressTensorxxEnum);  _assert_(sigma_xx_input);
+	sigma_xy_input  = inputs->GetInput(StressTensorxyEnum);  _assert_(sigma_xy_input);
+	sigma_yy_input  = inputs->GetInput(StressTensoryyEnum);  _assert_(sigma_yy_input);
+	damage_input  = this->material->inputs->GetInput(DamageDbarEnum);        _assert_(damage_input);
+
+	/*Damage evolution z mapping: */
+	gauss=new GaussTria();
+	J2s=0;
+	for (int iv=0;iv<NUMVERTICES;iv++){
+		gauss->GaussVertex(iv);
+		
+		damage_input->GetInputValue(&damage,gauss);
+		sigma_xx_input->GetInputValue(&sigma_xx,gauss);
+		sigma_xy_input->GetInputValue(&sigma_xy,gauss);
+		sigma_yy_input->GetInputValue(&sigma_yy,gauss);
+
+		s_xx=sigma_xx/(1-damage);
+		s_xy=sigma_xy/(1-damage);
+		s_yy=sigma_yy/(1-damage);
+
+		J2s=1.0/sqrt(2.0)*sqrt(pow(s_xx,2)+pow(s_yy,2)+2*pow(s_xy,2));
+		
+		Xis=sqrt(3.0)*J2s;
+
+		Psi=Xis-stress_threshold;
+
+		PosPsi=max(Psi,0.0);
+		NegPsi=max(-Psi,0.0);
+
+		f[iv]= c1* ( pow(PosPsi,c2)   -  healing * pow(NegPsi,c2) )  *   pow((1 - damage),-c3);
+
+	}
+	f_input=new TriaInput(NoneEnum,&f[0],P1Enum);
+	
+	/*Clean up and return*/
+	delete gauss;
+	return f_input;
+}
+/*}}}*/
+/*FUNCTION Tria::InputUpdateFromSolutionDamageEvolution {{{*/
+void  Tria::InputUpdateFromSolutionDamageEvolution(IssmDouble* solution){
+
+	const int    numdof=NDOF1*NUMVERTICES;
+
+	int         i;
+	IssmDouble  values[numdof];
+	int        *doflist = NULL;
+
+	/*Get dof list: */
+	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+
+	/*Use the dof list to index into the solution vector: */
+	for(i=0;i<numdof;i++){
+		values[i]=solution[doflist[i]];
+		/*Check solution*/
+		if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
+	}
+
+	/*Get all inputs and parameters*/
+	this->material->inputs->AddInput(new TriaInput(DamageDbarEnum,values,P1Enum));
+
+	/*Free ressources:*/
+	xDelete<int>(doflist);
+}
+/*}}}*/
+#endif
+
 #ifdef _HAVE_DAKOTA_
 /*FUNCTION Tria::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type);{{{*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16187)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16188)
@@ -218,4 +218,5 @@
 		void           GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
 		void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
+		void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
 		void           GetStrainRate2d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input);
 		void	         InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type);
@@ -276,6 +277,14 @@
 		bool    AnyActive(void);
 		#endif
-		#ifdef _HAVE_BALANCED_
-		#endif
+
+		#ifdef _HAVE_DAMAGE_
+		ElementMatrix* CreateKMatrixDamageEvolution(void);
+		ElementMatrix* CreateKMatrixDamageEvolution_CG(void);
+		ElementVector* CreatePVectorDamageEvolution(void);
+		ElementVector* CreatePVectorDamageEvolution_CG(void);
+		Input* DamageEvolutionF(void);
+		void	       InputUpdateFromSolutionDamageEvolution(IssmDouble* solution);
+		#endif
+
 
 		/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.h	(revision 16187)
+++ /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.h	(revision 16188)
@@ -72,4 +72,5 @@
 		IssmDouble Min(void){_error_("Min not implemented for booleans");};
 		IssmDouble MinAbs(void){_error_("Min not implemented for booleans");};
+		void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
 		void Scale(IssmDouble scale_factor);
 		void ArtificialNoise(IssmDouble min,IssmDouble max){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h	(revision 16187)
+++ /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h	(revision 16188)
@@ -72,4 +72,5 @@
 		void SquareMin(IssmDouble* psquaremin,Parameters* parameters){_error_("not implemented yet");};
 		void ConstrainMin(IssmDouble minimum){_error_("not implemented yet");};
+		void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
 		void Scale(IssmDouble scale_factor){_error_("not implemented yet");};
 		void ArtificialNoise(IssmDouble min,IssmDouble max){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.h	(revision 16187)
+++ /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.h	(revision 16188)
@@ -66,4 +66,5 @@
 		void SquareMin(IssmDouble* psquaremin,Parameters* parameters){_error_("not implemented yet");};
 		void ConstrainMin(IssmDouble minimum){_error_("not implemented yet");};
+		void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
 		void Scale(IssmDouble scale_factor){_error_("not implemented yet");};
 		void ArtificialNoise(IssmDouble min,IssmDouble max){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.h	(revision 16187)
+++ /issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.h	(revision 16188)
@@ -66,4 +66,5 @@
 		void SquareMin(IssmDouble* psquaremin,Parameters* parameters);
 		void ConstrainMin(IssmDouble minimum);
+		void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
 		void Scale(IssmDouble scale_factor);
 		void ArtificialNoise(IssmDouble min,IssmDouble max){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Inputs/Input.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/Input.h	(revision 16187)
+++ /issm/trunk-jpl/src/c/classes/Inputs/Input.h	(revision 16188)
@@ -56,4 +56,5 @@
 		virtual IssmDouble Max(void)=0;
 		virtual IssmDouble Min(void)=0;
+		virtual void Set(IssmDouble setvalue)=0;
 		virtual void   Scale(IssmDouble scale_factor)=0;
 		virtual void   ArtificialNoise(IssmDouble min,IssmDouble max)=0;
Index: /issm/trunk-jpl/src/c/classes/Inputs/IntInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/IntInput.h	(revision 16187)
+++ /issm/trunk-jpl/src/c/classes/Inputs/IntInput.h	(revision 16188)
@@ -68,4 +68,5 @@
 		void SquareMin(IssmDouble* psquaremin,Parameters* parameters);
 		void ConstrainMin(IssmDouble minimum){_error_("not implemented yet");};
+		void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
 		void Scale(IssmDouble scale_factor);
 		void ArtificialNoise(IssmDouble min,IssmDouble max){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h	(revision 16187)
+++ /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h	(revision 16188)
@@ -68,4 +68,5 @@
 		void SquareMin(IssmDouble* psquaremin,Parameters* parameters);
 		void ConstrainMin(IssmDouble minimum);
+		void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
 		void Scale(IssmDouble scale_factor);
 		void ArtificialNoise(IssmDouble min,IssmDouble max){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.h	(revision 16187)
+++ /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.h	(revision 16188)
@@ -71,4 +71,5 @@
 		void SquareMin(IssmDouble* psquaremin,Parameters* parameters);
 		void ConstrainMin(IssmDouble minimum){_error_("not implemented yet");};
+		void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
 		void Scale(IssmDouble scale_factor){_error_("not implemented yet");};
 		void ArtificialNoise(IssmDouble min,IssmDouble max){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp	(revision 16187)
+++ /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp	(revision 16188)
@@ -321,4 +321,11 @@
 }
 /*}}}*/
+/*FUNCTION TriaInput::Set{{{*/
+void TriaInput::Set(IssmDouble setvalue){
+
+	const int numnodes=this->NumberofNodes();
+	for(int i=0;i<numnodes;i++)values[i]=setvalue;
+}
+/*}}}*/
 /*FUNCTION TriaInput::ArtificialNoise{{{*/
 void TriaInput::ArtificialNoise(IssmDouble min,IssmDouble max){
@@ -435,4 +442,37 @@
 }
 /*}}}*/
+/*FUNCTION TriaInput::PointwiseDivide{{{*/
+Input* TriaInput::PointwiseDivide(Input* inputB){
+
+	/*Ouput*/
+	TriaInput* outinput=NULL;
+
+	/*Intermediaries*/
+	TriaInput *xinputB  = NULL;
+	const int   numnodes = this->NumberofNodes();
+
+	/*Check that inputB is of the same type*/
+	if(inputB->ObjectEnum()!=TriaInputEnum)     _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
+	xinputB=(TriaInput*)inputB;
+	if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
+
+	/*Allocate intermediary*/
+	IssmDouble* AdotBvalues=xNew<IssmDouble>(numnodes);
+
+	/*Create point wise division*/
+	for(int i=0;i<numnodes;i++){
+		_assert_(xinputB->values[i]!=0);
+		AdotBvalues[i]=this->values[i]/xinputB->values[i];
+	}
+
+	/*Create new Tria vertex input (copy of current input)*/
+	outinput=new TriaInput(this->enum_type,AdotBvalues,this->element_type);
+
+	/*Return output pointer*/
+	xDelete<IssmDouble>(AdotBvalues);
+	return outinput;
+
+}
+/*}}}*/
 /*FUNCTION TriaInput::Configure{{{*/
 void TriaInput::Configure(Parameters* parameters){
Index: /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h	(revision 16187)
+++ /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h	(revision 16188)
@@ -35,5 +35,5 @@
 		int    InstanceEnum();
 		Input* SpawnTriaInput(int location);
-		Input* PointwiseDivide(Input* inputB){_error_("not implemented yet");};
+		Input* PointwiseDivide(Input* inputB);
 		Input* PointwiseMin(Input* inputB);
 		Input* PointwiseMax(Input* inputB);
@@ -68,4 +68,5 @@
 		void SquareMin(IssmDouble* psquaremin,Parameters* parameters);
 		void ConstrainMin(IssmDouble minimum);
+		void Set(IssmDouble setvalue);
 		void Scale(IssmDouble scale_factor);
 		void ArtificialNoise(IssmDouble min,IssmDouble max);
Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 16187)
+++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 16188)
@@ -241,4 +241,9 @@
 		case HydrologyDCInefficientAnalysisEnum:
 			Ke=PenaltyCreateKMatrixHydrologyDCInefficient(kmax);
+			break;
+		#endif
+		#ifdef _HAVE_DAMAGE_
+		case DamageEvolutionAnalysisEnum:
+			Ke=PenaltyCreateKMatrixDamageEvolution(kmax);
 			break;
 		#endif
@@ -278,4 +283,10 @@
 			break;
 		#endif
+		#ifdef _HAVE_DAMAGE_
+		case DamageEvolutionAnalysisEnum:
+			pe=PenaltyCreatePVectorDamageEvolution(kmax);
+			break;
+		#endif
+
 		default:
 			_error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
@@ -430,4 +441,11 @@
 	}
 	#endif
+	#ifdef _HAVE_DAMAGE_
+	else if (analysis_type==DamageEvolutionAnalysisEnum){
+		ConstraintActivateDamageEvolution(punstable);
+		return;
+	}
+	#endif
+
 	else{
 		_error_("analysis: " << EnumToStringx(analysis_type) << " not supported yet");
@@ -695,4 +713,103 @@
 /*}}}*/
 #endif
+#ifdef _HAVE_DAMAGE_
+/*FUNCTION Pengrid::ConstraintActivateDamageEvolution {{{*/
+void  Pengrid::ConstraintActivateDamageEvolution(int* punstable){
+
+	//   The penalty is stable if it doesn't change during to successive iterations.   
+	IssmDouble max_damage;
+	IssmDouble damage;
+	int        new_active;
+	int        unstable=0;
+	int        penalty_lock;
+
+	/*check that pengrid is not a clone (penalty to be added only once)*/
+	if (node->IsClone()){
+		unstable=0;
+		*punstable=unstable;
+		return;
+	}
+
+	//First recover damage  using the element: */
+	element->GetMaterialInputValue(&damage,node,DamageDbarEnum);
+
+	//Recover our data:
+	parameters->FindParam(&penalty_lock,DamagePenaltyLockEnum);
+	parameters->FindParam(&max_damage,DamageMaxDamageEnum);
+	
+	//Figure out if damage is over max_damage, in which case, this penalty needs to be activated.
+
+	if (damage>max_damage){
+		new_active=1;
+	}
+	else{
+		new_active=0;
+	}
+
+	//Figure out stability of this penalty
+	if (active==new_active){
+		unstable=0;
+	}
+	else{
+		unstable=1;
+		if(penalty_lock)zigzag_counter++;
+	}
+
+	/*If penalty keeps zigzagging more than penalty_lock times: */
+	if(penalty_lock){
+		if(zigzag_counter>penalty_lock){
+			unstable=0;
+			active=1;
+		}
+	}
+
+	//Set penalty flag
+	active=new_active;
+
+	//*Assign output pointers:*/
+	*punstable=unstable;
+}
+/*}}}*/
+/*FUNCTION Pengrid::PenaltyCreateKMatrixDamageEvolution {{{*/
+ElementMatrix* Pengrid::PenaltyCreateKMatrixDamageEvolution(IssmDouble kmax){
+
+	IssmDouble    penalty_factor;
+
+	/*Initialize Element matrix and return if necessary*/
+	if(!this->active) return NULL;
+	ElementMatrix* Ke=new ElementMatrix(&node,NUMVERTICES,this->parameters);
+
+	/*recover parameters: */
+	parameters->FindParam(&penalty_factor,DamagePenaltyFactorEnum);
+
+	Ke->values[0]=kmax*pow(10.,penalty_factor);
+
+	/*Clean up and return*/
+	return Ke;
+}
+/*}}}*/
+/*FUNCTION Pengrid::PenaltyCreatePVectorDamageEvolution {{{*/
+ElementVector* Pengrid::PenaltyCreatePVectorDamageEvolution(IssmDouble kmax){
+
+	IssmDouble penalty_factor;
+	IssmDouble max_damage;
+
+	/*Initialize Element matrix and return if necessary*/
+	if(!this->active) return NULL;
+	ElementVector* pe=new ElementVector(&node,1,this->parameters);
+
+	//Recover our data:
+	parameters->FindParam(&penalty_factor,DamagePenaltyFactorEnum);
+	parameters->FindParam(&max_damage,DamageMaxDamageEnum);
+
+	//right hand side penalizes to max_damage
+	pe->values[0]=kmax*pow(10.,penalty_factor)*max_damage;
+
+	/*Clean up and return*/
+	return pe;
+}
+/*}}}*/
+#endif
+
 /*FUNCTION Pengrid::ResetConstraint {{{*/
 void  Pengrid::ResetConstraint(void){
Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h	(revision 16187)
+++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h	(revision 16188)
@@ -91,4 +91,10 @@
 		void           ConstraintActivateThermal(int* punstable);
 		#endif
+		#ifdef _HAVE_DAMAGE_
+		ElementMatrix* PenaltyCreateKMatrixDamageEvolution(IssmDouble kmax);
+		ElementVector* PenaltyCreatePVectorDamageEvolution(IssmDouble kmax);
+		void           ConstraintActivateDamageEvolution(int* punstable);
+		#endif
+
 		#ifdef _HAVE_HYDROLOGY_
 		ElementMatrix* PenaltyCreateKMatrixHydrologyDCInefficient(IssmDouble kmax);
Index: /issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStateLocal.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStateLocal.h	(revision 16187)
+++ /issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStateLocal.h	(revision 16188)
@@ -6,8 +6,20 @@
 #define _CONSTRAINTSSTATELOCAL_H
 
+
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
 #include "../../classes/classes.h"
 
 /*melting: */
 void  ThermalConstraintsState(Loads* loads, int* pconverged, int* pnum_unstable_constraints,int analysis_type);
+
+/*damage: */
+#ifdef _HAVE_DAMAGE_
+void  DamageConstraintsState(Loads* loads, int* pconverged, int* pnum_unstable_constraints,int analysis_type);
+#endif
 
 /*rifts module: */
Index: /issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStatex.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStatex.cpp	(revision 16187)
+++ /issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStatex.cpp	(revision 16188)
@@ -32,12 +32,25 @@
 	#ifdef _HAVE_RIFTS_
 	if(RiftIsPresent(femmodel->loads,analysis_type)){
+
 		RiftConstraintsState(&converged,&num_unstable_constraints,femmodel->loads,min_mechanical_constraints,analysis_type);
 	}
 	#endif
+
 	#ifdef _HAVE_THERMAL_
-	if(ThermalIsPresent(femmodel->loads,analysis_type)){
+
+	if(PengridIsPresent(femmodel->loads,ThermalAnalysisEnum)){
+
 		ThermalConstraintsState(femmodel->loads,&converged,&num_unstable_constraints,analysis_type);
 	}
 	#endif
+
+	#ifdef _HAVE_DAMAGE_
+
+	if(PengridIsPresent(femmodel->loads,DamageEvolutionAnalysisEnum)){
+
+		DamageConstraintsState(femmodel->loads,&converged,&num_unstable_constraints,analysis_type);
+	}
+	#endif
+
 
 	/*Assign output pointers: */
Index: /issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStatex.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStatex.h	(revision 16187)
+++ /issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStatex.h	(revision 16188)
@@ -9,5 +9,5 @@
 
 /* local prototypes: */
-int  ThermalIsPresent(Loads* loads,int analysis_type);
+int  PengridIsPresent(Loads* loads,int analysis_type);
 int  RiftIsPresent(Loads* loads,int analysis_type);
 void ConstraintsStatex(int* pconverged, int* pnum_unstable_constraints,FemModel* femmodel);
Index: /issm/trunk-jpl/src/c/modules/ConstraintsStatex/DamageConstraintsState.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ConstraintsStatex/DamageConstraintsState.cpp	(revision 16188)
+++ /issm/trunk-jpl/src/c/modules/ConstraintsStatex/DamageConstraintsState.cpp	(revision 16188)
@@ -0,0 +1,52 @@
+/*!\file:  DamageConstraintsState.cpp
+ * \brief  melting rate constraints
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "./ConstraintsStateLocal.h"
+
+void  DamageConstraintsState(Loads* loads, int* pconverged, int* pnum_unstable_constraints,int configuration_type){
+
+	int i;
+
+	int unstable=0;
+	int num_unstable_constraints=0;
+	int converged=0;
+	int sum_num_unstable_constraints=0;
+
+	num_unstable_constraints=0;	
+
+
+	/*Enforce constraints: */
+	for(i=0;i<loads->Size();i++){
+		Object* object=(Object*)loads->GetObjectByOffset(i);
+		Load* load=(Load*)object;
+		if(load->InAnalysis(configuration_type)){
+			if(object->ObjectEnum()==PengridEnum){
+
+				Pengrid* pengrid=(Pengrid*)object;
+
+				pengrid->ConstraintActivate(&unstable);
+				num_unstable_constraints+=unstable;
+			}
+		}
+	}
+
+
+	ISSM_MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+	ISSM_MPI_Bcast(&sum_num_unstable_constraints,1,ISSM_MPI_INT,0,IssmComm::GetComm());                
+	num_unstable_constraints=sum_num_unstable_constraints;
+
+	/*Have we converged? : */
+	if (num_unstable_constraints==0) converged=1;
+	else converged=0;
+
+	/*Assign output pointers: */
+	*pconverged=converged;
+	*pnum_unstable_constraints=num_unstable_constraints;
+}
Index: /issm/trunk-jpl/src/c/modules/ConstraintsStatex/PengridIsPresent.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ConstraintsStatex/PengridIsPresent.cpp	(revision 16188)
+++ /issm/trunk-jpl/src/c/modules/ConstraintsStatex/PengridIsPresent.cpp	(revision 16188)
@@ -0,0 +1,35 @@
+/*!\file:  ThermalIsPresent.cpp
+ * \brief
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "./ConstraintsStateLocal.h"
+
+int PengridIsPresent(Loads* loads,int configuration_type){
+
+	int i;
+	int found=0;
+	int mpi_found=0;
+
+	for(i=0;i<loads->Size();i++){
+		Object* object=(Object*)loads->GetObjectByOffset(i);
+		Load* load=(Load*)object;
+		if(load->InAnalysis(configuration_type)){
+			if (object->ObjectEnum()==PengridEnum){
+				found=1;
+				break;
+			}
+		}
+	}
+
+	ISSM_MPI_Reduce (&found,&mpi_found,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+	ISSM_MPI_Bcast(&mpi_found,1,ISSM_MPI_INT,0,IssmComm::GetComm());                
+	found=mpi_found;
+
+	return found;
+}
Index: sm/trunk-jpl/src/c/modules/ConstraintsStatex/ThermalIsPresent.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ConstraintsStatex/ThermalIsPresent.cpp	(revision 16187)
+++ 	(revision )
@@ -1,35 +1,0 @@
-/*!\file:  ThermalIsPresent.cpp
- * \brief
- */ 
-
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-#include "./ConstraintsStateLocal.h"
-
-int ThermalIsPresent(Loads* loads,int configuration_type){
-
-	int i;
-	int found=0;
-	int mpi_found=0;
-
-	for(i=0;i<loads->Size();i++){
-		Object* object=(Object*)loads->GetObjectByOffset(i);
-		Load* load=(Load*)object;
-		if(load->InAnalysis(configuration_type)){
-			if (object->ObjectEnum()==PengridEnum){
-				found=1;
-				break;
-			}
-		}
-	}
-
-	ISSM_MPI_Reduce (&found,&mpi_found,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&mpi_found,1,ISSM_MPI_INT,0,IssmComm::GetComm());                
-	found=mpi_found;
-
-	return found;
-}
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Damage/CreateLoadsDamage.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Damage/CreateLoadsDamage.cpp	(revision 16187)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Damage/CreateLoadsDamage.cpp	(revision 16188)
@@ -8,6 +8,25 @@
 
 void	CreateLoadsDamage(Loads** ploads, IoModel* iomodel){
+
+	/*Recover pointer: */
+	Loads* loads=*ploads;
+
+	/*create penalties for nodes: no node can have a damage > 1*/
+	iomodel->FetchData(1,DamageSpcdamageEnum);
+	CreateSingleNodeToElementConnectivity(iomodel);
+
+	for(int i=0;i<iomodel->numberofvertices;i++){
+
+		/*keep only this partition's nodes:*/
+		if((iomodel->my_vertices[i]==1)){
+			if (xIsNan<IssmDouble>(iomodel->Data(DamageSpcdamageEnum)[i])){ //No penalty applied on spc nodes!
+				loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,DamageEvolutionAnalysisEnum));
+			}
+		}
+	}
+	iomodel->DeleteData(1,DamageSpcdamageEnum);
+
+	/*Assign output pointer: */
+	*ploads=loads;
 	
-	/*No loads*/
-
 }
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Damage/CreateParametersDamage.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Damage/CreateParametersDamage.cpp	(revision 16187)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Damage/CreateParametersDamage.cpp	(revision 16188)
@@ -21,4 +21,10 @@
 	
 	iomodel->Constant(&law,DamageLawEnum);
+		
+	parameters->AddObject(iomodel->CopyConstantObject(DamagePenaltyThresholdEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(DamagePenaltyLockEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(DamagePenaltyFactorEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(DamageMaxiterEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(DamageMaxDamageEnum));
 
 	/*Retrieve law dependent parameters: */
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Damage/UpdateElementsDamage.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Damage/UpdateElementsDamage.cpp	(revision 16187)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Damage/UpdateElementsDamage.cpp	(revision 16188)
@@ -26,4 +26,6 @@
 	iomodel->FetchDataToInput(elements,VzEnum);
 	iomodel->FetchDataToInput(elements,DamageDEnum);
+	iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
+	iomodel->FetchDataToInput(elements,PressureEnum);
 
 	bool dakota_analysis;
Index: /issm/trunk-jpl/src/c/modules/ResetConstraintsx/DamageConstraintsReset.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ResetConstraintsx/DamageConstraintsReset.cpp	(revision 16188)
+++ /issm/trunk-jpl/src/c/modules/ResetConstraintsx/DamageConstraintsReset.cpp	(revision 16188)
@@ -0,0 +1,29 @@
+/*!\file:  DamageConstraintsState.cpp
+ * \brief  damage constraints reset
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "./ResetConstraintsx.h"
+
+void  DamageConstraintsReset(Loads* loads, int configuration_type){
+
+	int i;
+
+	/*Enforce constraints: */
+	for(i=0;i<loads->Size();i++){
+		Object* object=(Object*)loads->GetObjectByOffset(i);
+		Load* load=(Load*)object;
+		if(load->InAnalysis(configuration_type)){
+			if(object->ObjectEnum()==PengridEnum){
+
+				Pengrid* pengrid=(Pengrid*)object;
+				pengrid->ResetConstraint();
+			}
+		}
+	}
+}
Index: /issm/trunk-jpl/src/c/modules/ResetConstraintsx/ResetConstraintsx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ResetConstraintsx/ResetConstraintsx.cpp	(revision 16187)
+++ /issm/trunk-jpl/src/c/modules/ResetConstraintsx/ResetConstraintsx.cpp	(revision 16188)
@@ -31,7 +31,13 @@
 	#endif
 	#ifdef _HAVE_THERMAL_
-	if(ThermalIsPresent(femmodel->loads,analysis_type)){
+	if(PengridIsPresent(femmodel->loads,ThermalAnalysisEnum)){
 		ThermalConstraintsReset(femmodel->loads,analysis_type);
 	}
 	#endif
+	#ifdef _HAVE_DAMAGE_
+	if(PengridIsPresent(femmodel->loads,DamageEvolutionAnalysisEnum)){
+		DamageConstraintsReset(femmodel->loads,analysis_type);
+	}
+	#endif
+
 }
Index: /issm/trunk-jpl/src/c/modules/ResetConstraintsx/ResetConstraintsx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/ResetConstraintsx/ResetConstraintsx.h	(revision 16187)
+++ /issm/trunk-jpl/src/c/modules/ResetConstraintsx/ResetConstraintsx.h	(revision 16188)
@@ -10,4 +10,5 @@
 /* local prototypes: */
 void ThermalConstraintsReset(Loads* loads, int configuration_type);
+void DamageConstraintsReset(Loads* loads, int configuration_type);
 void ResetConstraintsx(FemModel* femmodel);
 
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 16187)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 16188)
@@ -163,5 +163,10 @@
 	DamageStressThresholdEnum,
 	DamageStabilizationEnum,
+	DamagePenaltyThresholdEnum,
+	DamagePenaltyLockEnum,
+	DamagePenaltyFactorEnum,
+	DamageMaxiterEnum,
 	DamageSpcdamageEnum,
+	DamageMaxDamageEnum,
 	MaterialsRhoIceEnum,
 	MaterialsRhoWaterEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 16187)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 16188)
@@ -171,5 +171,10 @@
 		case DamageStressThresholdEnum : return "DamageStressThreshold";
 		case DamageStabilizationEnum : return "DamageStabilization";
+		case DamagePenaltyThresholdEnum : return "DamagePenaltyThreshold";
+		case DamagePenaltyLockEnum : return "DamagePenaltyLock";
+		case DamagePenaltyFactorEnum : return "DamagePenaltyFactor";
+		case DamageMaxiterEnum : return "DamageMaxiter";
 		case DamageSpcdamageEnum : return "DamageSpcdamage";
+		case DamageMaxDamageEnum : return "DamageMaxDamage";
 		case MaterialsRhoIceEnum : return "MaterialsRhoIce";
 		case MaterialsRhoWaterEnum : return "MaterialsRhoWater";
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 16187)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 16188)
@@ -174,5 +174,10 @@
 	      else if (strcmp(name,"DamageStressThreshold")==0) return DamageStressThresholdEnum;
 	      else if (strcmp(name,"DamageStabilization")==0) return DamageStabilizationEnum;
+	      else if (strcmp(name,"DamagePenaltyThreshold")==0) return DamagePenaltyThresholdEnum;
+	      else if (strcmp(name,"DamagePenaltyLock")==0) return DamagePenaltyLockEnum;
+	      else if (strcmp(name,"DamagePenaltyFactor")==0) return DamagePenaltyFactorEnum;
+	      else if (strcmp(name,"DamageMaxiter")==0) return DamageMaxiterEnum;
 	      else if (strcmp(name,"DamageSpcdamage")==0) return DamageSpcdamageEnum;
+	      else if (strcmp(name,"DamageMaxDamage")==0) return DamageMaxDamageEnum;
 	      else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum;
 	      else if (strcmp(name,"MaterialsRhoWater")==0) return MaterialsRhoWaterEnum;
@@ -255,13 +260,13 @@
 	      else if (strcmp(name,"SurfaceforcingsHref")==0) return SurfaceforcingsHrefEnum;
 	      else if (strcmp(name,"SurfaceforcingsSmbref")==0) return SurfaceforcingsSmbrefEnum;
-	      else if (strcmp(name,"SurfaceforcingsBPos")==0) return SurfaceforcingsBPosEnum;
+         else stage=3;
+   }
+   if(stage==3){
+	      if (strcmp(name,"SurfaceforcingsBPos")==0) return SurfaceforcingsBPosEnum;
 	      else if (strcmp(name,"SurfaceforcingsBNeg")==0) return SurfaceforcingsBNegEnum;
 	      else if (strcmp(name,"ThermalIsenthalpy")==0) return ThermalIsenthalpyEnum;
 	      else if (strcmp(name,"ThermalIsdynamicbasalspc")==0) return ThermalIsdynamicbasalspcEnum;
 	      else if (strcmp(name,"ThermalMaxiter")==0) return ThermalMaxiterEnum;
-         else stage=3;
-   }
-   if(stage==3){
-	      if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum;
+	      else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum;
 	      else if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum;
 	      else if (strcmp(name,"ThermalPenaltyThreshold")==0) return ThermalPenaltyThresholdEnum;
@@ -378,13 +383,13 @@
 	      else if (strcmp(name,"IntParam")==0) return IntParamEnum;
 	      else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
-	      else if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
+         else stage=4;
+   }
+   if(stage==4){
+	      if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
 	      else if (strcmp(name,"IceFrontType")==0) return IceFrontTypeEnum;
 	      else if (strcmp(name,"SSA2dIceFront")==0) return SSA2dIceFrontEnum;
 	      else if (strcmp(name,"SSA3dIceFront")==0) return SSA3dIceFrontEnum;
 	      else if (strcmp(name,"Matice")==0) return MaticeEnum;
-         else stage=4;
-   }
-   if(stage==4){
-	      if (strcmp(name,"Matpar")==0) return MatparEnum;
+	      else if (strcmp(name,"Matpar")==0) return MatparEnum;
 	      else if (strcmp(name,"Node")==0) return NodeEnum;
 	      else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum;
@@ -501,13 +506,13 @@
 	      else if (strcmp(name,"BasalFriction")==0) return BasalFrictionEnum;
 	      else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum;
-	      else if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum;
+         else stage=5;
+   }
+   if(stage==5){
+	      if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum;
 	      else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum;
 	      else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
 	      else if (strcmp(name,"StressTensor")==0) return StressTensorEnum;
 	      else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum;
-         else stage=5;
-   }
-   if(stage==5){
-	      if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
+	      else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
 	      else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
 	      else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_damage_nonlinear.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_damage_nonlinear.cpp	(revision 16188)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_damage_nonlinear.cpp	(revision 16188)
@@ -0,0 +1,80 @@
+/*
+ * \brief: solutionsequence_damage_nonlinear.cpp: core of the damage solution 
+ */ 
+
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+
+void solutionsequence_damage_nonlinear(FemModel* femmodel){
+
+	/*solution : */
+	Vector<IssmDouble>* Dg=NULL; 
+	Vector<IssmDouble>* Df=NULL; 
+	Vector<IssmDouble>* Df_old=NULL; 
+	Vector<IssmDouble>* ys=NULL; 
+
+	/*intermediary: */
+	Matrix<IssmDouble>* Kff=NULL;
+	Matrix<IssmDouble>* Kfs=NULL;
+	Vector<IssmDouble>* pf=NULL;
+	Vector<IssmDouble>* df=NULL;
+
+	bool converged;
+	int constraints_converged;
+	int num_unstable_constraints;
+	int count;
+	int damage_penalty_threshold;
+	int damage_maxiter;
+
+	/*parameters:*/
+	int  configuration_type;
+
+	/*Recover parameters: */
+	femmodel->parameters->FindParam(&damage_penalty_threshold,DamagePenaltyThresholdEnum);
+	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+	femmodel->parameters->FindParam(&damage_maxiter,DamageMaxiterEnum);
+
+	count=1;
+	converged=false;
+
+	InputUpdateFromConstantx(femmodel,true,ResetPenaltiesEnum);
+	InputUpdateFromConstantx(femmodel,false,ConvergedEnum);
+	femmodel->UpdateConstraintsx();
+
+	for(;;){
+
+		delete Df_old; Df_old=Df;
+		SystemMatricesx(&Kff, &Kfs, &pf,&df, NULL,femmodel);
+		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
+		Reduceloadx(pf, Kfs, ys); delete Kfs;
+		Solverx(&Df, Kff, pf,Df_old, df, femmodel->parameters);
+		delete Kff;delete pf;delete Dg; delete df;
+		Mergesolutionfromftogx(&Dg, Df,ys,femmodel->nodes,femmodel->parameters); delete ys;
+		InputUpdateFromSolutionx(femmodel,Dg);
+
+		ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel);
+
+		if (!converged){
+			if(VerboseConvergence()) _printf0_("   #unstable constraints = " << num_unstable_constraints << "\n");
+			if (num_unstable_constraints <= damage_penalty_threshold)converged=true;
+			if (count>=damage_maxiter){
+				converged=true;
+				_printf0_("   maximum number of iterations (" << damage_maxiter << ") exceeded\n"); 
+			}
+		}
+		count++;
+
+		InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
+
+		if(converged)break;
+	}
+
+	InputUpdateFromSolutionx(femmodel,Dg);
+
+	/*Free ressources: */
+	delete Dg;
+	delete Df;
+	delete Df_old;
+}
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h	(revision 16187)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h	(revision 16188)
@@ -13,4 +13,5 @@
 
 void solutionsequence_thermal_nonlinear(FemModel* femmodel);
+void solutionsequence_damage_nonlinear(FemModel* femmodel);
 void solutionsequence_hydro_nonlinear(FemModel* femmodel);
 void solutionsequence_nonlinear(FemModel* femmodel,bool conserve_loads);
Index: /issm/trunk-jpl/src/m/classes/damage.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/damage.m	(revision 16187)
+++ /issm/trunk-jpl/src/m/classes/damage.m	(revision 16188)
@@ -11,5 +11,12 @@
 		law                 = '';
 		spcdamage           = NaN; 
+		max_damage          = NaN;
+	
+		%numerical
 		stabilization       = NaN;
+		penalty_threshold   = NaN;
+		maxiter             = NaN;
+		penalty_lock        = NaN;
+		penalty_factor      = NaN;
 		
 		%parameters for law 'initial': 
@@ -46,8 +53,22 @@
 			obj.D=0;
 			obj.law='undamaged';
+			
+			obj.max_damage=1-1e-5; %if damage reaches 1, solve becomes singular, as viscosity becomes nil
 		
 			%Type of stabilization used
 			obj.stabilization=2;
-	
+			
+			%Maximum number of iterations
+			obj.maxiter=100;
+
+			%factor used to compute the values of the penalties: kappa=max(stiffness matrix)*10^penalty_factor
+			obj.penalty_factor=3;
+			
+			%stabilize unstable damage constraints that keep zigzagging after n iteration (default is 0, no stabilization)
+			obj.penalty_lock=0;
+			
+			%threshold to declare convergence of thermal solution (default is 0)
+			obj.penalty_threshold=0;
+		
 			%damage parameters for 'initial' law of damage propagation
 			obj.stress_threshold=0;
@@ -61,8 +82,15 @@
 		function md = checkconsistency(obj,md,solution,analyses) % {{{
 			
-			md = checkfield(md,'damage.D','>=0',0,'size',[md.mesh.numberofvertices 1]);
+			md = checkfield(md,'damage.D','>=',0,'<=',obj.max_damage,'size',[md.mesh.numberofvertices 1]);
+			md = checkfield(md,'damage.max_damage','<',1,'>=',0);
 			md = checkfield(md,'damage.law','values',{'undamaged','pralong'});
 			md = checkfield(md,'damage.spcdamage','forcing',1);
+			
 			md = checkfield(md,'damage.stabilization','numel',[1],'values',[0 1 2]);
+			md = checkfield(md,'damage.maxiter','>=0',0);
+			md = checkfield(md,'damage.penalty_factor','>=0',0);
+			md = checkfield(md,'damage.penalty_lock','>=0',0);
+			md = checkfield(md,'damage.penalty_threshold','>=0',0);
+
 			if strcmpi(obj.law,'pralong'),
 				md = checkfield(md,'damage.healing','>=',0);
@@ -72,4 +100,10 @@
 				md = checkfield(md,'damage.c4','>=',0);
 				md = checkfield(md,'damage.stress_threshold','>=',0);
+			elseif strcmpi(obj.law,'undamaged'),
+				if (solution==DamageEvolutionSolutionEnum),
+					error('you are trying to evolve an undamaged material. Choose a valid evolution law (md.damage.law)');
+				end
+			else 
+				error('invalid damage evolution law');
 			end
 
@@ -81,5 +115,12 @@
 			fielddisplay(obj,'law','damage law (string) from {''undamaged'',''pralong''}');
 			fielddisplay(obj,'spcdamage','damage constraints (NaN means no constraint) [m]');
+			fielddisplay(obj,'max_damage','maximum damage sustained by ice (0<=max_damage<1)');
+			
 			fielddisplay(obj,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG');
+			fielddisplay(obj,'maxiter','maximum number of non linear iterations');
+			fielddisplay(obj,'penalty_lock','stabilize unstable damage constraints that keep zigzagging after n iteration (default is 0, no stabilization)');
+			fielddisplay(obj,'penalty_threshold','threshold to declare convergence of damage evolution solution (default is 0)');
+			fielddisplay(obj,'penalty_factor','scaling exponent (default is 3)');
+
 			if strcmpi(obj.law,'pralong'),
 				fielddisplay(obj,'c1','damage parameter 1');
@@ -96,6 +137,13 @@
 			WriteData(fid,'object',obj,'fieldname','D','format','DoubleMat','mattype',1);
 			WriteData(fid,'object',obj,'fieldname','law','format','String');
+			WriteData(fid,'object',obj,'fieldname','spcdamage','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
+			WriteData(fid,'object',obj,'fieldname','max_damage','format','Double');
+
 			WriteData(fid,'object',obj,'fieldname','stabilization','format','Integer');
-			WriteData(fid,'object',obj,'fieldname','spcdamage','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
+			WriteData(fid,'object',obj,'fieldname','penalty_threshold','format','Integer');
+			WriteData(fid,'object',obj,'fieldname','maxiter','format','Integer');
+			WriteData(fid,'object',obj,'fieldname','penalty_lock','format','Integer');
+			WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double');
+	
 			if strcmpi(obj.law,'pralong'),
 				WriteData(fid,'object',obj,'fieldname','c1','format','Double');
Index: /issm/trunk-jpl/src/m/classes/thermal.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/thermal.m	(revision 16187)
+++ /issm/trunk-jpl/src/m/classes/thermal.m	(revision 16188)
@@ -67,4 +67,5 @@
 			fielddisplay(obj,'penalty_lock','stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization)');
 			fielddisplay(obj,'penalty_threshold','threshold to declare convergence of thermal solution (default is 0)');
+			fielddisplay(obj,'penalty_factor','scaling exponent (default is 3)');
 			fielddisplay(obj,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)');
             		fielddisplay(obj,'isdynamicbasalspc',['enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)']);
Index: /issm/trunk-jpl/src/m/enum/DamageMaxDamageEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/DamageMaxDamageEnum.m	(revision 16188)
+++ /issm/trunk-jpl/src/m/enum/DamageMaxDamageEnum.m	(revision 16188)
@@ -0,0 +1,11 @@
+function macro=DamageMaxDamageEnum()
+%DAMAGEMAXDAMAGEENUM - Enum of DamageMaxDamage
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=DamageMaxDamageEnum()
+
+macro=StringToEnum('DamageMaxDamage');
Index: /issm/trunk-jpl/src/m/enum/DamageMaxiterEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/DamageMaxiterEnum.m	(revision 16188)
+++ /issm/trunk-jpl/src/m/enum/DamageMaxiterEnum.m	(revision 16188)
@@ -0,0 +1,11 @@
+function macro=DamageMaxiterEnum()
+%DAMAGEMAXITERENUM - Enum of DamageMaxiter
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=DamageMaxiterEnum()
+
+macro=StringToEnum('DamageMaxiter');
Index: /issm/trunk-jpl/src/m/enum/DamagePenaltyFactorEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/DamagePenaltyFactorEnum.m	(revision 16188)
+++ /issm/trunk-jpl/src/m/enum/DamagePenaltyFactorEnum.m	(revision 16188)
@@ -0,0 +1,11 @@
+function macro=DamagePenaltyFactorEnum()
+%DAMAGEPENALTYFACTORENUM - Enum of DamagePenaltyFactor
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=DamagePenaltyFactorEnum()
+
+macro=StringToEnum('DamagePenaltyFactor');
Index: /issm/trunk-jpl/src/m/enum/DamagePenaltyLockEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/DamagePenaltyLockEnum.m	(revision 16188)
+++ /issm/trunk-jpl/src/m/enum/DamagePenaltyLockEnum.m	(revision 16188)
@@ -0,0 +1,11 @@
+function macro=DamagePenaltyLockEnum()
+%DAMAGEPENALTYLOCKENUM - Enum of DamagePenaltyLock
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=DamagePenaltyLockEnum()
+
+macro=StringToEnum('DamagePenaltyLock');
Index: /issm/trunk-jpl/src/m/enum/DamagePenaltyThresholdEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/DamagePenaltyThresholdEnum.m	(revision 16188)
+++ /issm/trunk-jpl/src/m/enum/DamagePenaltyThresholdEnum.m	(revision 16188)
@@ -0,0 +1,11 @@
+function macro=DamagePenaltyThresholdEnum()
+%DAMAGEPENALTYTHRESHOLDENUM - Enum of DamagePenaltyThreshold
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=DamagePenaltyThresholdEnum()
+
+macro=StringToEnum('DamagePenaltyThreshold');
Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py
===================================================================
--- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 16187)
+++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 16188)
@@ -163,5 +163,10 @@
 def DamageStressThresholdEnum(): return StringToEnum("DamageStressThreshold")[0]
 def DamageStabilizationEnum(): return StringToEnum("DamageStabilization")[0]
+def DamagePenaltyThresholdEnum(): return StringToEnum("DamagePenaltyThreshold")[0]
+def DamagePenaltyLockEnum(): return StringToEnum("DamagePenaltyLock")[0]
+def DamagePenaltyFactorEnum(): return StringToEnum("DamagePenaltyFactor")[0]
+def DamageMaxiterEnum(): return StringToEnum("DamageMaxiter")[0]
 def DamageSpcdamageEnum(): return StringToEnum("DamageSpcdamage")[0]
+def DamageMaxDamageEnum(): return StringToEnum("DamageMaxDamage")[0]
 def MaterialsRhoIceEnum(): return StringToEnum("MaterialsRhoIce")[0]
 def MaterialsRhoWaterEnum(): return StringToEnum("MaterialsRhoWater")[0]
Index: /issm/trunk-jpl/test/NightlyRun/test271.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test271.m	(revision 16187)
+++ /issm/trunk-jpl/test/NightlyRun/test271.m	(revision 16188)
@@ -4,4 +4,5 @@
 md.damage.D=0.5*ones(md.mesh.numberofvertices,1);
 md.damage.spcdamage=NaN*ones(md.mesh.numberofvertices,1);
+md.damage.law='pralong';
 
 %boundary conditions for damage, to be put in SquareShelf.par
@@ -11,5 +12,5 @@
 md.damage.spcdamage(pos)=0;
 
-md.verbose=verbose('11111111');
+md.verbose=verbose('solution',true);
 
 md=setflowequation(md,'SSA','all');