Index: ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 18612)
+++ ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 18613)
@@ -64,13 +64,20 @@
 
 	/*parameters: */
 	bool save_results;
+	int  stabilization=1;
 	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+	//femmodel->parameters->FindParam(&stabilization,LevelsetStabilizationEnum);
 
 	/*activate formulation: */
 	femmodel->SetCurrentConfiguration(LevelsetAnalysisEnum);
 
 	if(VerboseSolution()) _printf0_("call computational core:\n");
-	solutionsequence_linear(femmodel);
+	if(stabilization==4){
+		solutionsequence_fct(femmodel);
+	}
+	else{
+		solutionsequence_linear(femmodel);
+	}
 
 	if(save_results){
 		if(VerboseSolution()) _printf0_("   saving results\n");
Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18612)
+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18613)
@@ -203,6 +203,7 @@
 	DamageEvolutionRequestedOutputsEnum,
 	DamageEnum,
 	NewDamageEnum,
+	StressIntensityFactorEnum,
 	MaterialsRhoIceEnum,
 	MaterialsRhoSeawaterEnum,
 	MaterialsRhoFreshwaterEnum,
@@ -399,6 +400,7 @@
 	MeshdeformationSolutionEnum,
 	MeshdeformationAnalysisEnum,
 	LevelsetAnalysisEnum,
+	LevelsetStabilizationEnum,
 	ExtrapolationAnalysisEnum,
 	LsfReinitializationAnalysisEnum,
 	/*}}}*/
Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18612)
+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18613)
@@ -211,6 +211,7 @@
 		case DamageEvolutionRequestedOutputsEnum : return "DamageEvolutionRequestedOutputs";
 		case DamageEnum : return "Damage";
 		case NewDamageEnum : return "NewDamage";
+		case StressIntensityFactorEnum : return "StressIntensityFactor";
 		case MaterialsRhoIceEnum : return "MaterialsRhoIce";
 		case MaterialsRhoSeawaterEnum : return "MaterialsRhoSeawater";
 		case MaterialsRhoFreshwaterEnum : return "MaterialsRhoFreshwater";
@@ -402,6 +403,7 @@
 		case MeshdeformationSolutionEnum : return "MeshdeformationSolution";
 		case MeshdeformationAnalysisEnum : return "MeshdeformationAnalysis";
 		case LevelsetAnalysisEnum : return "LevelsetAnalysis";
+		case LevelsetStabilizationEnum : return "LevelsetStabilization";
 		case ExtrapolationAnalysisEnum : return "ExtrapolationAnalysis";
 		case LsfReinitializationAnalysisEnum : return "LsfReinitializationAnalysis";
 		case ApproximationEnum : return "Approximation";
Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18612)
+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18613)
@@ -214,6 +214,7 @@
 	      else if (strcmp(name,"DamageEvolutionRequestedOutputs")==0) return DamageEvolutionRequestedOutputsEnum;
 	      else if (strcmp(name,"Damage")==0) return DamageEnum;
 	      else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
+	      else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
 	      else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum;
 	      else if (strcmp(name,"MaterialsRhoSeawater")==0) return MaterialsRhoSeawaterEnum;
 	      else if (strcmp(name,"MaterialsRhoFreshwater")==0) return MaterialsRhoFreshwaterEnum;
@@ -258,11 +259,11 @@
 	      else if (strcmp(name,"MasstransportVertexPairing")==0) return MasstransportVertexPairingEnum;
 	      else if (strcmp(name,"MasstransportNumRequestedOutputs")==0) return MasstransportNumRequestedOutputsEnum;
 	      else if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum;
-	      else if (strcmp(name,"QmuIsdakota")==0) return QmuIsdakotaEnum;
          else stage=3;
    }
    if(stage==3){
-	      if (strcmp(name,"MassFluxSegments")==0) return MassFluxSegmentsEnum;
+	      if (strcmp(name,"QmuIsdakota")==0) return QmuIsdakotaEnum;
+	      else if (strcmp(name,"MassFluxSegments")==0) return MassFluxSegmentsEnum;
 	      else if (strcmp(name,"MassFluxSegmentsPresent")==0) return MassFluxSegmentsPresentEnum;
 	      else if (strcmp(name,"QmuMassFluxSegmentsPresent")==0) return QmuMassFluxSegmentsPresentEnum;
 	      else if (strcmp(name,"QmuNumberofpartitions")==0) return QmuNumberofpartitionsEnum;
@@ -381,11 +382,11 @@
 	      else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum;
 	      else if (strcmp(name,"FlaimAnalysis")==0) return FlaimAnalysisEnum;
 	      else if (strcmp(name,"FlaimSolution")==0) return FlaimSolutionEnum;
-	      else if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum;
          else stage=4;
    }
    if(stage==4){
-	      if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum;
+	      if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum;
+	      else if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum;
 	      else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum;
 	      else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
 	      else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
@@ -411,6 +412,7 @@
 	      else if (strcmp(name,"MeshdeformationSolution")==0) return MeshdeformationSolutionEnum;
 	      else if (strcmp(name,"MeshdeformationAnalysis")==0) return MeshdeformationAnalysisEnum;
 	      else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
+	      else if (strcmp(name,"LevelsetStabilization")==0) return LevelsetStabilizationEnum;
 	      else if (strcmp(name,"ExtrapolationAnalysis")==0) return ExtrapolationAnalysisEnum;
 	      else if (strcmp(name,"LsfReinitializationAnalysis")==0) return LsfReinitializationAnalysisEnum;
 	      else if (strcmp(name,"Approximation")==0) return ApproximationEnum;
@@ -503,12 +505,12 @@
 	      else if (strcmp(name,"Adjointx")==0) return AdjointxEnum;
 	      else if (strcmp(name,"Adjointy")==0) return AdjointyEnum;
 	      else if (strcmp(name,"Adjointz")==0) return AdjointzEnum;
-	      else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum;
-	      else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
+	      if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum;
+	      else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
+	      else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
 	      else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
 	      else if (strcmp(name,"Converged")==0) return ConvergedEnum;
 	      else if (strcmp(name,"Fill")==0) return FillEnum;
@@ -626,12 +628,12 @@
 	      else if (strcmp(name,"J")==0) return JEnum;
 	      else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum;
 	      else if (strcmp(name,"Step")==0) return StepEnum;
-	      else if (strcmp(name,"Time")==0) return TimeEnum;
-	      else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum;
+	      if (strcmp(name,"Time")==0) return TimeEnum;
+	      else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
+	      else if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum;
 	      else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum;
 	      else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum;
 	      else if (strcmp(name,"MassfluxatgateName")==0) return MassfluxatgateNameEnum;
@@ -749,12 +751,12 @@
 	      else if (strcmp(name,"BasalforcingsOceanTurningAngle")==0) return BasalforcingsOceanTurningAngleEnum;
 	      else if (strcmp(name,"BasalforcingsOceanSsh")==0) return BasalforcingsOceanSshEnum;
 	      else if (strcmp(name,"BasalforcingsOceanVx")==0) return BasalforcingsOceanVxEnum;
-	      else if (strcmp(name,"BasalforcingsOceanVy")==0) return BasalforcingsOceanVyEnum;
-	      else if (strcmp(name,"SurfaceforcingsRhoAir")==0) return SurfaceforcingsRhoAirEnum;
          else stage=7;
    }
    if(stage==7){
-	      if (strcmp(name,"SurfaceforcingsAirCoef")==0) return SurfaceforcingsAirCoefEnum;
+	      if (strcmp(name,"BasalforcingsOceanVy")==0) return BasalforcingsOceanVyEnum;
+	      else if (strcmp(name,"SurfaceforcingsRhoAir")==0) return SurfaceforcingsRhoAirEnum;
+	      else if (strcmp(name,"SurfaceforcingsAirCoef")==0) return SurfaceforcingsAirCoefEnum;
 	      else if (strcmp(name,"SurfaceforcingsAirLinDragCoef")==0) return SurfaceforcingsAirLinDragCoefEnum;
 	      else if (strcmp(name,"SurfaceforcingsAirQuadDragCoef")==0) return SurfaceforcingsAirQuadDragCoefEnum;
 	      else if (strcmp(name,"SurfaceforcingsWindVx")==0) return SurfaceforcingsWindVxEnum;
Index: ../trunk-jpl/src/c/shared/Elements/elements.h
===================================================================
--- ../trunk-jpl/src/c/shared/Elements/elements.h	(revision 18612)
+++ ../trunk-jpl/src/c/shared/Elements/elements.h	(revision 18613)
@@ -21,6 +21,7 @@
 				     IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday, 
 					  IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout);
 IssmDouble DrainageFunctionWaterfraction(IssmDouble waterfraction, IssmDouble dt=0.);
+IssmDouble StressIntensityIntegralWeight(IssmDouble depth, IssmDouble water_depth, IssmDouble thickness);
 
 /*Print arrays*/
 void printarray(IssmPDouble* array,int lines,int cols=1);
Index: ../trunk-jpl/src/c/shared/Elements/StressIntensityIntegralWeight.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Elements/StressIntensityIntegralWeight.cpp	(revision 0)
+++ ../trunk-jpl/src/c/shared/Elements/StressIntensityIntegralWeight.cpp	(revision 18613)
@@ -0,0 +1,31 @@
+/* \file StressIntensityIntegralWeight.cpp
+ * \Weight to integrate the stress along the ice flow direction to compute the stress intensity factor : Linear Fracture Mechanics (see "Combining damage and fracture mechanics to model calving",Krug 2014 in the appendix) 
+ */
+
+#include <math.h>
+
+#include "../Numerics/types.h"
+
+IssmDouble StressIntensityIntegralWeight(IssmDouble depth, IssmDouble water_depth, IssmDouble thickness){
+
+	/*output: */
+	IssmDouble beta;
+
+	/*intermediaries: */
+	IssmDouble M1,M2,M3,x,y,d;
+	const double pi = 3.141592653589793;
+	x    = water_depth/thickness;
+	y    = depth;
+	d    = water_depth;
+
+	M1   = 0.0719768-1.513476*x-61.1001*pow(x,2)+1554.95*pow(x,3)-14583.8*pow(x,4)+71590.7*pow(x,5)-205384*pow(x,6)+356469*pow(x,7)-368270*pow(x,8)+208233*pow(x,9)-49544*pow(x,10);
+	//printf("M1 : %g",M1);
+	M2   = 0.246984+6.47583*x+176.456*pow(x,2)-4058.76*pow(x,3)+37303.8*pow(x,4)-181755*pow(x,5)+520551*pow(x,6)-904370*pow(x,7)+936863*pow(x,8)-531940*pow(x,9)+127291*pow(x,10);
+	//printf("M2 : %g",M2);
+	M3   = 0.529659-22.3235*x+532.074*pow(x,2)-5479.53*pow(x,3)+28592.2*pow(x,4)-81388.6*pow(x,5)+128746*pow(x,6)-106246*pow(x,7)+35780.7*pow(x,8);
+	//printf("M3 : %g",M3);
+
+	beta = 2/sqrt(2*pi*(d-y))*(1+M1*sqrt(1-y/d)+M2*(1-y/d)+M3*pow((1-y/d),1.5));
+
+	return beta;
+}
Index: ../trunk-jpl/src/c/Makefile.am
===================================================================
--- ../trunk-jpl/src/c/Makefile.am	(revision 18612)
+++ ../trunk-jpl/src/c/Makefile.am	(revision 18613)
@@ -212,6 +212,7 @@
 					./shared/Sorting/sorting.h\
 					./shared/Elements/elements.h\
 					./shared/Elements/Cuffey.cpp\
+					./shared/Elements/StressIntensityIntegralWeight.cpp\
 					./shared/Elements/Paterson.cpp\
 					./shared/Elements/Arrhenius.cpp\
 					./shared/Elements/LliboutryDuval.cpp \
Index: ../trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 18612)
+++ ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 18613)
@@ -1635,6 +1635,18 @@
 	xDelete<int>(control_type);
 }
 /*}}}*/
+void FemModel::StressIntensityFactorx(){/*{{{*/
+
+	/*Initialize input as 0*/
+	InputUpdateFromConstantx(this,0.,StressIntensityFactorEnum);
+
+	/*Update input for basal element only*/
+	for(int i=0;i<elements->Size();i++){
+		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+		element->StressIntensityFactor();
+	}
+}
+	/*}}}*/
 #ifdef  _HAVE_DAKOTA_
 void FemModel::DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses){/*{{{*/
 
Index: ../trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.h	(revision 18612)
+++ ../trunk-jpl/src/c/classes/Elements/Element.h	(revision 18613)
@@ -218,6 +218,7 @@
 		virtual void   ComputeBasalStress(Vector<IssmDouble>* sigma_b)=0;
 		virtual void   ComputeStressTensor(void)=0;
 		virtual void   ComputeDeviatoricStressTensor(void)=0;
+		virtual void	StressIntensityFactor(void)=0;
 
 		virtual void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finite_element)=0;
 		virtual void   InputExtrude(int input_enum,int start)=0;
Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.h	(revision 18612)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.h	(revision 18613)
@@ -55,6 +55,7 @@
 		void        ComputeStressTensor();
 		void        ComputeDeviatoricStressTensor();
 		void        ComputeSurfaceNormalVelocity();
+		void        StressIntensityFactor(void){_error_("not implemented yet");};
 		void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
 		void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
 		void        ResetHooks();
Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18612)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18613)
@@ -346,6 +346,83 @@
 	delete gauss;
 }
 /*}}}*/
+void       Penta::StressIntensityFactor(){/*{{{*/
+
+	/* Check if we are on the base */
+	if(!IsOnBase()) return;
+
+	IssmDouble  ki[6]={0.};
+	IssmDouble  const_grav=9.81;
+	IssmDouble  rho_ice=900;
+	IssmDouble  rho_water=1000;
+	IssmDouble  Jdet[3];
+	IssmDouble  pressure,vx,vy,vel,deviaxx,deviaxy,deviayy,water_depth,prof,stress_xx,thickness;
+
+	Penta* penta=this;
+	for(;;){
+	
+		IssmDouble  xyz_list[NUMVERTICES][3];
+		/* Get node coordinates and dof list: */
+		::GetVerticesCoordinates(&xyz_list[0][0],penta->vertices,NUMVERTICES);
+
+		///*Compute the Jacobian for the vertical integration*/
+		Jdet[0]=(xyz_list[3][2]-xyz_list[0][2])*0.5;
+		Jdet[1]=(xyz_list[4][2]-xyz_list[1][2])*0.5;
+		Jdet[2]=(xyz_list[5][2]-xyz_list[2][2])*0.5;
+	
+		/*Retrieve all inputs we will need*/
+		Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
+		Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
+		Input* vel_input=inputs->GetInput(VelEnum);                                _assert_(vel_input);
+		Input* pressure_input=inputs->GetInput(PressureEnum);                      _assert_(pressure_input);
+		Input* deviaxx_input=inputs->GetInput(DeviatoricStressxxEnum);             _assert_(deviaxx_input);
+		Input* deviaxy_input=inputs->GetInput(DeviatoricStressxyEnum);             _assert_(deviaxy_input);
+		Input* deviayy_input=inputs->GetInput(DeviatoricStressyyEnum);             _assert_(deviayy_input);
+		Input* surface_input=inputs->GetInput(SurfaceEnum);								_assert_(surface_input);
+		Input* thickness_input=inputs->GetInput(ThicknessEnum);							_assert_(thickness_input);
+		
+		/* Start looping on the number of 2D vertices: */
+		for(int ig=0;ig<3;ig++){
+			GaussPenta* gauss=new GaussPenta(ig,3+ig,11);
+			for (int iv=gauss->begin();iv<gauss->end();iv++){
+				gauss->GaussPoint(iv);
+
+				/* Get the value we need*/
+				pressure_input->GetInputValue(&pressure,gauss);
+				vx_input->GetInputValue(&vx,gauss);
+				vy_input->GetInputValue(&vy,gauss);
+				vel_input->GetInputValue(&vel,gauss);
+				deviaxx_input->GetInputValue(&deviaxx,gauss);
+				deviaxy_input->GetInputValue(&deviaxy,gauss);
+				deviayy_input->GetInputValue(&deviayy,gauss);
+				surface_input->GetInputValue(&water_depth,gauss);
+				thickness_input->GetInputValue(&thickness,gauss);
+				prof=water_depth-penta->GetZcoord(&xyz_list[0][0],gauss);
+
+				/*stress_xx= Deviatoric stress along the ice flow direction plus cryostatic pressure */
+				stress_xx=(vx*vx*(deviaxx)+vy*vy*(deviayy)+2*vy*vx*deviaxy)/(vel*vel+1.e-6);
+
+				if(prof<water_depth&prof<thickness){
+					/* Compute the local stress intensity factor*/ 
+					ki[ig]+=Jdet[ig]*gauss->weight*stress_xx*StressIntensityIntegralWeight(prof,min(water_depth,thickness),thickness);
+				}
+			}
+			delete gauss;
+		}
+			
+		/*Stop if we have reached the surface/base*/
+		if(penta->IsOnSurface()) break;
+		
+		/*get upper Penta*/
+		penta=penta->GetUpperPenta();
+		_assert_(penta->Id()!=this->id);
+	}
+
+	/*Add input*/
+	this->inputs->AddInput(new PentaInput(StressIntensityFactorEnum,&ki[0],P1Enum));
+	this->InputExtrude(StressIntensityFactorEnum,-1);
+}
+/*}}}*/
 void       Penta::Configure(Elements* elementsin, Loads* loadsin, Nodes* nodesin,Vertices* verticesin, Materials* materialsin, Parameters* parametersin){/*{{{*/
 
 	int analysis_counter;
Index: ../trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.h	(revision 18612)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.h	(revision 18613)
@@ -56,6 +56,7 @@
 		void   ComputeSigmaNN(){_error_("not implemented yet");};
 		void   ComputeStressTensor();
 		void   ComputeDeviatoricStressTensor();
+		void   StressIntensityFactor();
 		void   Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
 		void   ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
 		int    FiniteElement(void);
Index: ../trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Seg.h	(revision 18612)
+++ ../trunk-jpl/src/c/classes/Elements/Seg.h	(revision 18613)
@@ -56,6 +56,7 @@
 		void        ComputeSigmaNN(){_error_("not implemented yet");};
 		void        ComputeStressTensor(){_error_("not implemented yet");};
 		void        ComputeDeviatoricStressTensor(){_error_("not implemented yet");};
+		void        StressIntensityFactor(void){_error_("not implemented yet");};
 		void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
 		void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
 		void        ResetHooks(){_error_("not implemented yet");};
Index: ../trunk-jpl/src/c/classes/Elements/Tetra.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 18612)
+++ ../trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 18613)
@@ -56,6 +56,7 @@
 		void        ComputeSigmaNN(){_error_("not implemented yet");};
 		void        ComputeStressTensor(){_error_("not implemented yet");};
 		void        ComputeDeviatoricStressTensor(){_error_("not implemented yet");};
+		void        StressIntensityFactor(void){_error_("not implemented yet");};
 		void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
 		void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
 		void        ResetHooks();
Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18612)
+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18613)
@@ -1104,6 +1104,7 @@
 				name==EplHeadEnum ||
 				name==SedimentHeadOldEnum ||
 				name==EplHeadOldEnum ||
+				name==StressIntensityFactorEnum ||
 				name==HydrologydcEplThicknessOldEnum ||
 				name==HydrologydcEplInitialThicknessEnum ||
 				name==HydrologydcEplThicknessEnum ||
@@ -1202,6 +1203,10 @@
 				this->ComputeNewDamage();
 				input=this->inputs->GetInput(output_enum);
 				break;
+			case StressIntensityFactorEnum:
+				this->StressIntensityFactor();
+				input=this->inputs->GetInput(output_enum);
+				break;
 			default:
 				_error_("input "<<EnumToStringx(output_enum)<<" not found in element");
 		}
Index: ../trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.h	(revision 18612)
+++ ../trunk-jpl/src/c/classes/FemModel.h	(revision 18613)
@@ -80,6 +80,7 @@
 		void IceVolumeAboveFloatationx(IssmDouble* pV);
 		void ElementResponsex(IssmDouble* presponse,int response_enum);
 		void BalancethicknessMisfitx(IssmDouble* pV);
+		void StressIntensityFactorx();
 		#ifdef  _HAVE_DAKOTA_
 		void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses);
 		#endif
Index: ../trunk-jpl/src/m/enum/StressIntensityFactorEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/StressIntensityFactorEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/StressIntensityFactorEnum.m	(revision 18613)
@@ -0,0 +1,11 @@
+function macro=StressIntensityFactorEnum()
+%STRESSINTENSITYFACTORENUM - Enum of StressIntensityFactor
+%
+%   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=StressIntensityFactorEnum()
+
+macro=StringToEnum('StressIntensityFactor');
Index: ../trunk-jpl/src/m/enum/EnumDefinitions.py
===================================================================
--- ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18612)
+++ ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18613)
@@ -203,6 +203,7 @@
 def DamageEvolutionRequestedOutputsEnum(): return StringToEnum("DamageEvolutionRequestedOutputs")[0]
 def DamageEnum(): return StringToEnum("Damage")[0]
 def NewDamageEnum(): return StringToEnum("NewDamage")[0]
+def StressIntensityFactorEnum(): return StringToEnum("StressIntensityFactor")[0]
 def MaterialsRhoIceEnum(): return StringToEnum("MaterialsRhoIce")[0]
 def MaterialsRhoSeawaterEnum(): return StringToEnum("MaterialsRhoSeawater")[0]
 def MaterialsRhoFreshwaterEnum(): return StringToEnum("MaterialsRhoFreshwater")[0]
@@ -394,6 +395,7 @@
 def MeshdeformationSolutionEnum(): return StringToEnum("MeshdeformationSolution")[0]
 def MeshdeformationAnalysisEnum(): return StringToEnum("MeshdeformationAnalysis")[0]
 def LevelsetAnalysisEnum(): return StringToEnum("LevelsetAnalysis")[0]
+def LevelsetStabilizationEnum(): return StringToEnum("LevelsetStabilization")[0]
 def ExtrapolationAnalysisEnum(): return StringToEnum("ExtrapolationAnalysis")[0]
 def LsfReinitializationAnalysisEnum(): return StringToEnum("LsfReinitializationAnalysis")[0]
 def ApproximationEnum(): return StringToEnum("Approximation")[0]
Index: ../trunk-jpl/src/m/enum/LevelsetStabilizationEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/LevelsetStabilizationEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/LevelsetStabilizationEnum.m	(revision 18613)
@@ -0,0 +1,11 @@
+function macro=LevelsetStabilizationEnum()
+%LEVELSETSTABILIZATIONENUM - Enum of LevelsetStabilization
+%
+%   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=LevelsetStabilizationEnum()
+
+macro=StringToEnum('LevelsetStabilization');
