Index: ../trunk-jpl/src/c/analyses/SeaiceAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/SeaiceAnalysis.h	(revision 18501)
+++ ../trunk-jpl/src/c/analyses/SeaiceAnalysis.h	(revision 18502)
@@ -31,6 +31,7 @@
 		void UpdateConstraints(FemModel* femmodel);
 
 		/*Sea ice specifics*/
+		void UpdateDamageAndStress(FemModel* femmodel);
 		void CreateCTensor(IssmDouble* C,Element* element,Gauss* gauss);
 		void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
 		void GetM(IssmDouble* M,Element* element,IssmDouble* xyz_list,Gauss* gauss);
Index: ../trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp	(revision 18501)
+++ ../trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp	(revision 18502)
@@ -320,8 +320,49 @@
 	xDelete<int>(doflist);
 }/*}}}*/
 void SeaiceAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
-	/*Default, do nothing*/
-	return;
+
+	/*Intermediaries*/
+	Vector<IssmDouble>* mask        = NULL;
+	IssmDouble*         serial_mask = NULL;
+
+	/*Step 1: update mask of active nodes*/
+	mask=new Vector<IssmDouble>(femmodel->nodes->NumberOfNodes(SeaiceAnalysisEnum));
+	for (int i=0;i<femmodel->elements->Size();i++){
+		Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+
+		/*Get current concentration of element and decide whether it is an active element*/
+		int    numnodes            = element->GetNumberOfNodes();
+		Input* concentration_input = element->GetInput(SeaiceConcentrationEnum); _assert_(concentration_input);
+		if(concentration_input->Max()>0.){
+			for(int in=0;in<numnodes;in++) mask->SetValue(element->nodes[in]->Sid(),1.,INS_VAL);
+		}
+	}
+
+	/*Assemble and serialize*/
+	mask->Assemble();
+	serial_mask=mask->ToMPISerial();
+	delete mask;
+
+	/*Update node activation accordingly*/
+	int counter =0;
+	for(int i=0;i<femmodel->nodes->Size();i++){
+		Node* node=dynamic_cast<Node*>(femmodel->nodes->GetObjectByOffset(i));
+		if(node->InAnalysis(SeaiceAnalysisEnum)){
+			if(serial_mask[node->Sid()]==1.){
+				node->Activate();
+				counter++;
+			}
+			else{
+				node->Deactivate();
+			}
+		}
+	}
+	xDelete<IssmDouble>(serial_mask);
+
+	/*Display number of active nodes*/
+	int sum_counter;
+	ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+	if(VerboseSolution()) _printf0_("   Number of active nodes: "<< sum_counter <<"\n");
 }/*}}}*/
 
 /*Intermediary functions*/
@@ -421,3 +462,122 @@
 	/*Clean-up*/
 	xDelete<IssmDouble>(basis);
 }/*}}}*/
+void SeaiceAnalysis::UpdateDamageAndStress(FemModel* femmodel){/*{{{*/
+	/* The damage variable is updated as a function of the actual elastic deformation
+	 * In both cases, a Coulombic enveloppe is used, define by the cohesion C, tan(phi) and tract_coef.
+	 * In both cases, a maximal compressive strength is fixed at compr_max
+	 * The enveloppe is defined in N/m^2.
+	 * The coeficients of the internal stress are then multiplied by the ice thickness to be used in the vertical integrated momentiom equation.
+	 */
+
+	/* Mohr-Coulomb criterion calculation
+	 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	 *                                                                   
+	 *                            sigma_s                                
+	 *        Mohr.Coulomb branch     |                                  
+	 *                     \          |                                  
+	 *                      * |       |                                  
+	 *                        *       |   cohesion (C=cfix+calea)        
+	 *                        | *     |  /                               
+	 *                        |   *   | /    tract                       
+	 *                        |     * |/    /                            
+	 *             -comp_max  |       *    /                             
+	 *                      \ |       | * /                              
+	 *                       \|      0| | *                              
+	 *             -------------------------*------------ sigma_n        
+	 *                        |       | | *                              
+	 *                        |       | *                                
+	 *                        |       *                                  
+	 *                        |     * |                                  
+	 *                        |   *   |                                  
+	 *                        | *     |                                  
+	 *                        *       |                                  
+	 *                      * |       |                                  
+	 *                    *           |                                  
+	 *                                |                                  
+	 *                                |                                  
+	 *                                                                   
+	 * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	 */
+
+	/*Intermediaties*/
+	IssmDouble sigma_xx,sigma_yy,sigma_xy,sigma_s,sigma_n,sigma_target;
+	IssmDouble compression_coef,traction_coef,cohesion,tan_phi;
+	IssmDouble traction,compression_max;
+	IssmDouble damage_test,damage_new,damage,tmp;
+
+	/*Loop over the elements of this partition and update accordingly*/
+	for(int i=0;i<femmodel->elements->Size();i++){
+		Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+
+		/*Get Mohr-Coulomb parameters*/
+		compression_coef = element->GetMaterialParameter(MaterialsCompressionCoefEnum);
+		traction_coef    = element->GetMaterialParameter(MaterialsTractionCoefEnum);
+		cohesion         = element->GetMaterialParameter(MaterialsCohesionEnum); //C
+		tan_phi          = element->GetMaterialParameter(MaterialsInternalFrictionCoefEnum); //mu = tan(phi)
+
+		/*Get current stress state*/
+		Input*   sigma_xx_input = element->GetInput(StressTensorxxEnum);
+		Input*   sigma_yy_input = element->GetInput(StressTensoryyEnum);
+		Input*   sigma_xy_input = element->GetInput(StressTensorxyEnum);
+		if(sigma_xx_input) sigma_xx_input->GetInputAverage(&sigma_xx);
+		else               sigma_xx = 0.;
+		if(sigma_yy_input) sigma_yy_input->GetInputAverage(&sigma_yy);
+		else               sigma_yy = 0.;
+		if(sigma_xy_input) sigma_xy_input->GetInputAverage(&sigma_xy);
+		else               sigma_xy = 0.;
+
+		/* Compute the invariants of the elastic deformation and instantaneous deformation rate */
+		sigma_s=sqrt(pow((sigma_xx-sigma_yy)/2.,2)+pow(sigma_xy,2));
+		sigma_n=(sigma_xx+sigma_yy)/2.;
+
+		/* same maximal tensile strength */
+		traction=traction_coef*cohesion/tan_phi;
+		compression_max=compression_coef*cohesion;
+
+		/* estimate the internal constraints using the current elastic deformation */
+		Input* damage_input    = element->GetInput(MaterialsDamageEnum); _assert_(damage_input);
+		Input* damagenew_input = element->GetInput(MaterialsDamageEnum); _assert_(damagenew_input);//FIXME
+		damage_input->GetInputAverage(&damage);
+		damagenew_input->GetInputAverage(&damage_new);
+		damage_test = damage;
+		if(sigma_n>traction || sigma_n<-compression_max){
+			if(sigma_n>traction){
+				sigma_target=traction;
+			}
+			else{
+				sigma_target=-compression_max;
+			}
+
+			tmp=1.-sigma_target/sigma_n*(1.-damage);
+			if(tmp<1. && tmp>damage_new){
+				damage_test=((damage_test>tmp)?(damage_test):(tmp)); /* max(damage_test,tmp); */
+			}
+		}
+		if(sigma_s>cohesion-sigma_n*tan_phi){
+			tmp=1.0-cohesion/(sigma_s+sigma_n*tan_phi)*(1-damage);
+			if(tmp<1. && tmp>damage_new){
+				damage_test=((damage_test>tmp)?(damage_test):(tmp)); /*max(damage_test,tmp); */
+			}
+		}
+
+		/* The damage variable is changed */
+		damage_new=damage_test;
+		element->AddInput(MaterialsDamageEnum,&damage_test,P0Enum);//FIXME
+
+		/* Recompute the internal stress*/
+		if(damage<1.){
+			sigma_xx = (1.-damage_new)/(1.-damage)*sigma_xx;
+			sigma_yy = (1.-damage_new)/(1.-damage)*sigma_yy;
+			sigma_xy = (1.-damage_new)/(1.-damage)*sigma_xy;
+		}
+		else{
+			sigma_xx = 0.;
+			sigma_yy = 0.;
+			sigma_xy = 0.;
+		}
+		element->AddInput(StressTensorxxEnum,&sigma_xx,P0Enum);
+		element->AddInput(StressTensoryyEnum,&sigma_yy,P0Enum);
+		element->AddInput(StressTensorxyEnum,&sigma_xy,P0Enum);
+	}
+}/*}}}*/
Index: ../trunk-jpl/src/c/cores/seaice_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/seaice_core.cpp	(revision 18501)
+++ ../trunk-jpl/src/c/cores/seaice_core.cpp	(revision 18502)
@@ -17,7 +17,7 @@
 	/*recover parameters: */
 	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
 
-	if(VerboseSolution()) _printf0_("call computational core:\n");
+	/*Launch solution sequence for the only analysis we are interested in*/
 	femmodel->SetCurrentConfiguration(SeaiceAnalysisEnum);
 	solutionsequence_linear(femmodel);
 
