Index: ../trunk-jpl/src/c/analyses/SeaiceAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/SeaiceAnalysis.h	(revision 18504)
+++ ../trunk-jpl/src/c/analyses/SeaiceAnalysis.h	(revision 18505)
@@ -35,5 +35,6 @@
 		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);
+		void PostProcess(FemModel* femmodel);
 };
 #endif
Index: ../trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp	(revision 18504)
+++ ../trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp	(revision 18505)
@@ -140,7 +140,7 @@
 		this->CreateCTensor(&C[0][0],element,gauss);
 		for(int i=0;i<3;i++){
 			for(int j=0;j<3;j++){
-				C[i][j] = dt*time_ratio*gauss->weight*Jdet*C[i][j];
+				C[i][j] = dt*thickness*time_ratio*gauss->weight*Jdet*C[i][j];
 			}
 		}
 		TripleMultiply(B,3,numdof,1,
@@ -319,8 +319,7 @@
 void SeaiceAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
 	/*Intermediaries*/
-	int*        doflist=NULL;
-	IssmDouble* xyz_list=NULL;
+	int* doflist=NULL;
 
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = element->GetNumberOfNodes();
@@ -357,7 +356,6 @@
 	xDelete<IssmDouble>(vy);
 	xDelete<IssmDouble>(vx);
 	xDelete<IssmDouble>(values);
-	xDelete<IssmDouble>(xyz_list);
 	xDelete<int>(doflist);
 }/*}}}*/
 void SeaiceAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
@@ -415,7 +413,7 @@
 
 	/*Intermediaries*/
 	IssmDouble D[3][3];
-	IssmDouble damage,thickness,concentration;
+	IssmDouble damage,concentration;
 
 	/*Get Material parameters*/
 	IssmDouble poisson          = element->GetMaterialParameter(MaterialsPoissonEnum);
@@ -424,10 +422,8 @@
 
 	/*Get damage input at this location*/
 	Input* damage_input        = element->GetInput(DamageEnum);              _assert_(damage_input);
-	Input* thickness_input     = element->GetInput(SeaiceThicknessEnum);     _assert_(thickness_input);
 	Input* concentration_input = element->GetInput(SeaiceConcentrationEnum); _assert_(concentration_input);
 	damage_input->GetInputValue(&damage,gauss);
-	thickness_input->GetInputValue(&thickness,gauss);
 	concentration_input->GetInputValue(&concentration,gauss);
 
 	/*Build D tensor once for all*/
@@ -447,7 +443,7 @@
 	/*Let's build C = hE*D */
 	for(int i=0;i<3;i++){
 		for(int j=0;j<3;j++){
-			C[i*3+j] = thickness*E*D[i][j];
+			C[i*3+j] = E*D[i][j];
 		}
 	}
 }/*}}}*/
@@ -622,3 +618,151 @@
 		element->AddInput(StressTensorxyEnum,&sigma_xy,P0Enum);
 	}
 }/*}}}*/
+void SeaiceAnalysis::PostProcess(FemModel* femmodel){/*{{{*/
+
+	/*Intermediaties*/
+	IssmDouble *xyz_list = NULL;
+	IssmDouble *xyz_list_old = NULL;
+	IssmDouble *xyz_list_new = NULL;
+	IssmDouble  dt,meshx,meshy,vx,vy,eps_xx,eps_yy,eps_zz,area_old,area_new;
+	IssmDouble  time_relaxation_stress,time_relaxation_damage,damage,concentration,thickness;
+	IssmDouble  sigma_dot_xx,sigma_dot_yy,sigma_dot_xy;
+	IssmDouble  sigma_pred_xx,sigma_pred_yy,sigma_pred_xy;
+	IssmDouble  sigma_xx,sigma_yy,sigma_xy;
+	IssmDouble  min_c,min_h,max_h;
+	IssmDouble  epsilon[3]; /*eps_xx,eps_yy,eps_xy*/
+	IssmDouble  C[3][3];
+
+	/*Fetch the parameters once for all*/
+	femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+	femmodel->parameters->FindParam(&min_c,SeaiceMinConcentrationEnum);
+	femmodel->parameters->FindParam(&min_h,SeaiceMinThicknessEnum);
+	femmodel->parameters->FindParam(&max_h,SeaiceMaxThicknessEnum);
+
+	/*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 some inputs needed for the update*/
+		element->GetVerticesCoordinates(&xyz_list);
+		time_relaxation_stress = element->GetMaterialParameter(MaterialsTimeRelaxationStressEnum);
+		time_relaxation_damage = element->GetMaterialParameter(MaterialsTimeRelaxationDamageEnum);
+		Input* meshx_input = element->GetInput(MeshXEnum); _assert_(meshx_input);
+		Input* meshy_input = element->GetInput(MeshYEnum); _assert_(meshy_input);
+		Input* vx_input    = element->GetInput(VxEnum);    _assert_(vx_input);
+		Input* vy_input    = element->GetInput(VyEnum);    _assert_(vy_input);
+
+		/*Preallocate future inputs*/
+		int numvertices=element->GetNumberOfVertices();
+		IssmDouble* meshx_new_list = xNew<IssmDouble>(numvertices);
+		IssmDouble* meshy_new_list = xNew<IssmDouble>(numvertices);
+		xyz_list_old = xNew<IssmDouble>(numvertices*3);
+		xyz_list_new = xNew<IssmDouble>(numvertices*3);
+
+		/*1. update vertex positions (using a Gauss object for convenience)*/
+		Gauss* gauss=element->NewGauss();
+		for (int iv=0;iv<numvertices;iv++){
+			gauss->GaussVertex(iv);
+
+			meshx_input->GetInputValue(&meshx,gauss);
+			meshy_input->GetInputValue(&meshy,gauss);
+			vx_input->GetInputValue(&vx,gauss);
+			vy_input->GetInputValue(&vy,gauss);
+
+			meshx_new_list[iv] = meshx + vx*dt;
+			meshy_new_list[iv] = meshy + vy*dt;
+
+			xyz_list_old[iv*3+0] = meshx;
+			xyz_list_old[iv*3+1] = meshy;
+			xyz_list_old[iv*3+2] = 0.;
+
+			xyz_list_new[iv*3+0] = meshx_new_list[iv];
+			xyz_list_new[iv*3+1] = meshy_new_list[iv];
+			xyz_list_new[iv*3+2] = 0.;
+		}
+		element->AddInput(MeshXEnum,meshx_new_list,P1Enum);
+		element->AddInput(MeshYEnum,meshy_new_list,P1Enum);
+
+		/*Now we are going to use a point in the center of the element, because we are P0 for the stress*/
+		gauss->GaussNode(P0Enum,0);
+
+		/*Calculate sigma_dot*/
+		sigma_dot_xx=sigma_dot_yy=sigma_dot_xy=0.;
+		this->CreateCTensor(&C[0][0],element,gauss);
+		element->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
+		for(int j=0;j<3;j++){
+			sigma_dot_xx += C[0][j]*epsilon[0];
+			sigma_dot_yy += C[1][j]*epsilon[1];
+			sigma_dot_xy += C[2][j]*epsilon[2];
+		}
+
+		/*Get predicted stress state*/
+		Input* sigma_xx_input = element->GetInput(StressTensorPredictorxxEnum); _assert_(sigma_xx_input);
+		Input* sigma_yy_input = element->GetInput(StressTensorPredictoryyEnum); _assert_(sigma_yy_input);
+		Input* sigma_xy_input = element->GetInput(StressTensorPredictorxyEnum); _assert_(sigma_xy_input);
+		sigma_xx_input->GetInputAverage(&sigma_xx);
+		sigma_yy_input->GetInputAverage(&sigma_yy);
+		sigma_xy_input->GetInputAverage(&sigma_xy);
+
+		/*Calculate new stress and push to element*/
+		sigma_xx = sigma_xx + dt*(sigma_dot_xx - sigma_xx/time_relaxation_stress);
+		sigma_yy = sigma_yy + dt*(sigma_dot_yy - sigma_yy/time_relaxation_stress);
+		sigma_xy = sigma_xy + dt*(sigma_dot_xy - sigma_xy/time_relaxation_stress);
+		element->AddInput(StressTensorxxEnum,&sigma_xx,P0Enum);
+		element->AddInput(StressTensoryyEnum,&sigma_yy,P0Enum);
+		element->AddInput(StressTensorxyEnum,&sigma_xy,P0Enum);
+
+		/*Update Damage According*/
+		Input* damage_input = element->GetInput(DamageEnum); _assert_(damage_input);
+		damage_input->GetInputAverage(&damage);
+		damage = damage*(1.-dt/time_relaxation_damage);
+		element->AddInput(DamageEnum,&damage,P0Enum);
+
+		/*Prepare new predictor*/
+		this->CreateCTensor(&C[0][0],element,gauss); /*C is updated now that there is a new damage in inputs*/
+		for(int j=0;j<3;j++){
+			sigma_pred_xx += C[0][j]*epsilon[0];
+			sigma_pred_yy += C[1][j]*epsilon[1];
+			sigma_pred_xy += C[2][j]*epsilon[2];
+		}
+		element->AddInput(StressTensorPredictorxxEnum,&sigma_pred_xx,P0Enum);
+		element->AddInput(StressTensorPredictoryyEnum,&sigma_pred_yy,P0Enum);
+		element->AddInput(StressTensorPredictorxyEnum,&sigma_pred_xy,P0Enum);
+
+		/*Calculate Old and new area (FIXME: for now we assume trianlges...)*/
+		element->JacobianDeterminant(&area_old,xyz_list_old,gauss);
+		element->JacobianDeterminant(&area_new,xyz_list_new,gauss);
+		area_old = area_old/SQRT3;
+		area_new = area_new/SQRT3;
+
+		/*Update ice thickness and concentration using element distortion*/
+		Input* thickness_input     = element->GetInput(SeaiceThicknessEnum);     _assert_(thickness_input);
+		Input* concentration_input = element->GetInput(SeaiceConcentrationEnum); _assert_(concentration_input);
+		thickness_input->GetInputAverage(&thickness);
+		concentration_input->GetInputAverage(&concentration);
+		IssmDouble ice_area   = concentration*area_old;
+		IssmDouble ice_volume = thickness*area_old;
+		if(concentration>min_c){
+			concentration = ice_area  /area_new;
+			thickness     = ice_volume/area_new;
+
+			/* lower bounds */
+			concentration = ((concentration>min_c)?(concentration):(min_c));
+			thickness     = ((thickness    >min_h)?(thickness    ):(min_h));
+
+			/* upper bounds (only for the concentration) */
+			concentration    = ((concentration<1.   )?(concentration):(1));
+			thickness        = ((thickness    <max_h)?(thickness    ):(max_h));
+		}
+		element->AddInput(SeaiceConcentrationEnum,&concentration,P0Enum);
+		element->AddInput(SeaiceThicknessEnum,&thickness,P0Enum);
+
+		/*Clean up*/
+		xDelete<IssmDouble>(meshx_new_list);
+		xDelete<IssmDouble>(meshy_new_list);
+		xDelete<IssmDouble>(xyz_list);
+		xDelete<IssmDouble>(xyz_list_old);
+		xDelete<IssmDouble>(xyz_list_new);
+	}
+
+}/*}}}*/
Index: ../trunk-jpl/src/c/cores/seaice_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/seaice_core.cpp	(revision 18504)
+++ ../trunk-jpl/src/c/cores/seaice_core.cpp	(revision 18505)
@@ -12,21 +12,26 @@
 void seaice_core(FemModel* femmodel){
 
 	/*parameters: */
-	bool save_results;
+	bool           save_results;
+	SeaiceAnalysis analysis = SeaiceAnalysis();
 
 	/*recover parameters: */
 	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
 
 	/*Launch solution sequence for the only analysis we are interested in*/
 	femmodel->SetCurrentConfiguration(SeaiceAnalysisEnum);
-	SeaiceAnalysis* analysis = new SeaiceAnalysis();
-	analysis->UpdateDamageAndStress(femmodel);
+	analysis.UpdateDamageAndStress(femmodel);
 	solutionsequence_linear(femmodel);
+	analysis.PostProcess(femmodel);
 
 	if(save_results){
 		if(VerboseSolution()) _printf0_("   saving results\n");
-		int outputs[6] = {VxEnum,VyEnum,VelEnum,SeaiceConcentrationEnum,SeaiceThicknessEnum,DamageEnum};
-		femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],6);
+		int outputs[14] = {VxEnum,VyEnum,VelEnum,
+			SeaiceConcentrationEnum,SeaiceThicknessEnum,DamageEnum,
+			StressTensorPredictorxxEnum,StressTensorPredictoryyEnum,StressTensorPredictorxyEnum,
+			StressTensorxxEnum,StressTensoryyEnum,StressTensorxyEnum,
+			MeshXEnum,MeshYEnum};
+		femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],14);
 	}
 
 }
