Index: /issm/trunk-jpl/src/m/classes/damage.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/damage.m	(revision 16197)
+++ /issm/trunk-jpl/src/m/classes/damage.m	(revision 16198)
@@ -20,5 +20,5 @@
 		penalty_factor      = NaN;
 		
-		%parameters for law 'initial': 
+		%general parameters for evolution law: 
 		stress_threshold    = NaN;
 		c1                  = NaN;
@@ -71,5 +71,5 @@
 			obj.penalty_threshold=0;
 		
-			%damage parameters for 'initial' law of damage propagation
+			%damage evolution parameters 
 			obj.stress_threshold=0;
 			obj.healing=0;
@@ -102,5 +102,5 @@
 			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)');
+					error('Invalid evolution law (md.damage.law) for a damage solution');
 				end
 			else 
@@ -115,5 +115,5 @@
 			fielddisplay(obj,'law','damage law (string) from {''undamaged'',''pralong''}');
 			fielddisplay(obj,'spcdamage','damage constraints (NaN means no constraint)');
-			fielddisplay(obj,'max_damage','maximum damage sustained by ice (0<=max_damage<1)');
+			fielddisplay(obj,'max_damage','maximum possible damage (0<=max_damage<1)');
 			
 			fielddisplay(obj,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG');
Index: /issm/trunk-jpl/src/m/classes/damage.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/damage.py	(revision 16197)
+++ /issm/trunk-jpl/src/m/classes/damage.py	(revision 16198)
@@ -16,14 +16,23 @@
 			
 		#damage: 
-		isdamage = float('NaN')
-		self.D   = float('NaN')
-		self.law   = ''
+		self.D						= float('NaN')
+		self.law						= ''
+		self.spcdamage				= float('NaN')
+		self.max_damage			= float('NaN')
 		
-		#parameters for law 'initial': 
-		self.stress_threshold    = float('NaN')
-		self.c1                  = float('NaN')
-		self.c2                  = float('NaN')
-		self.c3                  = float('NaN')
-		self.c4                  = float('NaN')
+		#numerical
+		stabilization				= float('NaN')
+		penalty_threshold			= float('NaN')
+		maxiter						= float('NaN')
+		penalty_lock				= float('NaN')
+		penalty_factor				= float('NaN')
+
+		#general parameters for evolution law: 
+		self.stress_threshold   = float('NaN')
+		self.c1                 = float('NaN')
+		self.c2                 = float('NaN')
+		self.c3                 = float('NaN')
+		self.c4                 = float('NaN')
+		self.healing				= float('NaN')
 
 		if not len(args):
@@ -38,4 +47,12 @@
 		s+="%s\n" % fielddisplay(self,"D","damage tensor (scalar for now)")
 		s+="%s\n" % fielddisplay(self,"law","damage law (string) from ['undamaged','pralong']")
+		s+="%s\n" % fielddisplay(self,"spcdamage","damage constraints (NaN means no constraint)")
+		s+="%s\n" % fielddisplay(self,"max_damage","maximum possible damage (0<=max_damage<1)")
+
+		s+="%s\n" % fielddisplay(self,"stabilization","0: no, 1: artificial_diffusivity, 2: SUPG")
+		s+="%s\n" % fielddisplay(self,"maxiter","maximum number of non linear iterations")
+		s+="%s\n" % fielddisplay(self,"penalty_lock","stabilize unstable damage constraints that keep zigzagging after n iteration (default is 0, no stabilization)")
+		s+="%s\n" % fielddisplay(self,"penalty_lock","threshold to declare convergence of damage evolution solution (default is 0)")
+		s+="%s\n" % fielddisplay(self,"penalty_lock","scaling exponent (default is 3)")
 
 		if (self.law=='pralong'):
@@ -51,7 +68,25 @@
 
 		#damage parameters: 
-		self.isdamage=0
 		self.D=0
 		self.law='undamaged'
+
+		self.max_damage=1-1e-5 #if damage reaches 1, solve becomes singular, as viscosity becomes nil
+		
+		#Type of stabilization used
+		self.stabilization=2
+			
+		#Maximum number of iterations
+		self.maxiter=100
+
+		#factor used to compute the values of the penalties: kappa=max(stiffness matrix)*10^penalty_factor
+		self.penalty_factor=3
+			
+		#stabilize unstable damage constraints that keep zigzagging after n iteration (default is 0, no stabilization)
+		self.penalty_lock=0
+			
+		#threshold to declare convergence of thermal solution (default is 0)
+		self.penalty_threshold=0
+		
+		#damage evolution parameters 
 		self.stress_threshold=0
 		self.c1=0
@@ -59,12 +94,22 @@
 		self.c3=0
 		self.c4=0
+		self.healing=0
 
 	# }}}
 	def checkconsistency(self,md,solution,analyses):    # {{{
 
-		md = checkfield(md,'damage.D','>=',0,'size',[md.mesh.numberofvertices])
-		if self.isdamage:
-			md = checkfield(md,'damage.law','values',['undamaged','pralong'])
+		md = checkfield(md,'damage.D','>=',0,'<=',self.max_damage,'size',[md.mesh.numberofvertices])
+		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 self.law == 'pralong':
+			md = checkfield(md,'damage.healing','>=',0);
 			md = checkfield(md,'damage.c1','>=',0)
 			md = checkfield(md,'damage.c2','>=',0)
@@ -72,4 +117,7 @@
 			md = checkfield(md,'damage.c4','>=',0)
 			md = checkfield(md,'damage.stress_threshold','>=',0)
+		elif strcmpi(self.law,'undamaged'):
+			if (solution==DamageEvolutionSolutionEnum):
+				raise RuntimeError('Invalid evolution law (md.damage.law) for a damage solution');
 
 		return md
@@ -77,8 +125,15 @@
 	def marshall(self,md,fid):    # {{{
 
-		WriteData(fid,'object',self,'fieldname','isdamage','format','DoubleMat','mattype',1)
 		WriteData(fid,'object',self,'fieldname','D','format','DoubleMat','mattype',1)
-		if self.isdamage:
-			WriteData(fid,'object',self,'fieldname','law','format','String')
+		WriteData(fid,'object',self,'fieldname','law','format','String')
+		WriteData(fid,'object',self,'fieldname','spcdamage','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
+		WriteData(fid,'object',self,'fieldname','max_damage','format','Double');
+
+		WriteData(fid,'object',self,'fieldname','stabilization','format','Integer');
+		WriteData(fid,'object',self,'fieldname','penalty_threshold','format','Integer');
+		WriteData(fid,'object',self,'fieldname','maxiter','format','Integer');
+		WriteData(fid,'object',self,'fieldname','penalty_lock','format','Integer');
+		WriteData(fid,'object',self,'fieldname','penalty_factor','format','Double');
+
 		if self.law=='pralong':
 			WriteData(fid,'object',self,'fieldname','c1','format','Double')
Index: /issm/trunk-jpl/src/m/classes/model/model.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/model/model.py	(revision 16197)
+++ /issm/trunk-jpl/src/m/classes/model/model.py	(revision 16198)
@@ -616,4 +616,5 @@
 		md.masstransport.spcthickness=project3d(md,'vector',md.masstransport.spcthickness,'type','node')
 		md.balancethickness.spcthickness=project3d(md,'vector',md.balancethickness.spcthickness,'type','node')
+		md.damage.spcdamage=project3d(md,'vector',md.damage.spcdamage,'type','node')
 		md.stressbalance.referential=project3d(md,'vector',md.stressbalance.referential,'type','node')
 		md.stressbalance.loadingforce=project3d(md,'vector',md.stressbalance.loadingforce,'type','node')
@@ -633,4 +634,5 @@
 		#damage
 		md.damage.D=project3d(md,'vector',md.damage.D,'type','node')
+
 		#parameters
 		md.geometry.surface=project3d(md,'vector',md.geometry.surface,'type','node')