Index: ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 19032)
+++ ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 19033)
@@ -67,6 +67,10 @@
 				iomodel->FetchDataToInput(elements,CalvingpiCoeffEnum);
 				iomodel->FetchDataToInput(elements,CalvingpiMeltingrateEnum);
 				break;
+			case CalvingDevEnum:
+				iomodel->FetchDataToInput(elements,CalvingpiCoeffEnum);
+				iomodel->FetchDataToInput(elements,CalvingMeltingrateEnum);
+				break;
 			default:
 				_error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
 		}
@@ -173,18 +177,19 @@
 	Input* meltingrate_input  = NULL;
 
 	/*Load velocities*/
-	if(domaintype==Domain2DhorizontalEnum){
-		vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
-		vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
-	}
-	else{
-		if(dim==1){
+	switch(domaintype){
+		case Domain2DverticalEnum:
 			vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
-		}
-		if(dim==2){
+			break;
+		case Domain2DhorizontalEnum:
+			vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
+			vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
+			break;
+		case Domain3DEnum:
 			vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input);
 			vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input);
-		}
+			break;
+		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
 	}
 
 	/*Load calving inputs*/
@@ -197,37 +202,56 @@
 				meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
 				break;
 			case CalvingLevermannEnum:
-				if(domaintype==Domain2DhorizontalEnum){
-					calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
-					calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
-				}
-				else{
-					if(dim==1){
+				switch(domaintype){
+					case Domain2DverticalEnum:
 						calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
-					}
-					if(dim==2){
+						break;
+					case Domain2DhorizontalEnum:
+						calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
+						calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
+						break;
+					case Domain3DEnum:
 						calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
 						calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
-					}
+						break;
+					default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
 				}
 				meltingrate_input = basalelement->GetInput(CalvinglevermannMeltingrateEnum);     _assert_(meltingrate_input);
 				break;
 			case CalvingPiEnum:
-				if(domaintype==Domain2DhorizontalEnum){
-					calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
-					calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
-				}
-				else{
-					if(dim==1){
+				switch(domaintype){
+					case Domain2DverticalEnum:
 						calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
-					}
-					if(dim==2){
+						break;
+					case Domain2DhorizontalEnum:
+						calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
+						calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
+						break;
+					case Domain3DEnum:
 						calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
 						calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
-					}
+						break;
+					default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
 				}
 				meltingrate_input = basalelement->GetInput(CalvingpiMeltingrateEnum);     _assert_(meltingrate_input);
 				break;
+			case CalvingDevEnum:
+				switch(domaintype){
+					case Domain2DverticalEnum:
+						calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
+						break;
+					case Domain2DhorizontalEnum:
+						calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
+						calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
+						break;
+					case Domain3DEnum:
+						calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
+						calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
+						break;
+					default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+				}
+				meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
+				break;
 			default:
 				_error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
 		}
@@ -254,7 +278,7 @@
 		/* Advection */
 		GetB(B,basalelement,xyz_list,gauss); 
 		GetBprime(Bprime,basalelement,xyz_list,gauss); 
-		vx_input->GetInputValue(&v[0],gauss); // in 3D case, add mesh velocity 
+		vx_input->GetInputValue(&v[0],gauss);
 		vy_input->GetInputValue(&v[1],gauss); 
 
 		/*Get calving speed*/
@@ -272,43 +296,26 @@
 
 					if(norm_dlsf>1.e-10)
 					 for(i=0;i<dim;i++){ 
-					  c[i]=calvingrate*dlsf[i]/norm_dlsf;
-					  m[i]=meltingrate*dlsf[i]/norm_dlsf;
+					  c[i]=calvingrate*dlsf[i]/norm_dlsf; m[i]=meltingrate*dlsf[i]/norm_dlsf;
 					 }
 					else
 					 for(i=0;i<dim;i++){
-						 c[i]=0.;
-						 m[i]=0.;
+						 c[i]=0.; m[i]=0.;
 					 }
-
 					break;
 
 				case CalvingLevermannEnum:
-					calvingratex_input->GetInputValue(&c[0],gauss); // in 3D case, add mesh velocity 
+				case CalvingPiEnum:
+				case CalvingDevEnum:
+					calvingratex_input->GetInputValue(&c[0],gauss);
 					if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
 					meltingrate_input->GetInputValue(&meltingrate,gauss);
-
 					norm_calving=0.;
 					for(i=0;i<dim;i++) norm_calving+=pow(c[i],2);
 					norm_calving=sqrt(norm_calving)+1.e-14;
-
 					for(i=0;i<dim;i++) m[i]=meltingrate*c[i]/norm_calving;
-					
 					break;
 
-				case CalvingPiEnum:
-					calvingratex_input->GetInputValue(&c[0],gauss); // in 3D case, add mesh velocity
-					if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
-					meltingrate_input->GetInputValue(&meltingrate,gauss);
-
-					norm_calving=0.;
-					for(i=0;i<dim;i++) norm_calving+=c[i]*c[i];
-					norm_calving=sqrt(norm_calving)+1.e-14;
-
-					for(i=0;i<dim;i++) m[i]=meltingrate*c[i]/norm_calving;
-
-					break;
-
 				default:
 					_error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
 			}
Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 19032)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 19033)
@@ -1567,7 +1567,7 @@
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		element->NodalFunctions(basis, gauss);
 
-		thickness_input->GetInputValue(&thickness,gauss);
+		thickness_input->GetInputValue(&thickness,gauss); _assert_(thickness>0);
 		surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
 
 		for(int i=0;i<numnodes;i++){
Index: ../trunk-jpl/src/c/cores/transient_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/transient_core.cpp	(revision 19032)
+++ ../trunk-jpl/src/c/cores/transient_core.cpp	(revision 19033)
@@ -102,19 +102,31 @@
 		}
 
 		if(islevelset){
-			if(iscalving && calvinglaw==CalvingLevermannEnum){
-				if(VerboseSolution()) _printf0_("   computing calving rate\n");
-				femmodel->StrainRateparallelx();
-				femmodel->StrainRateperpendicularx();
-				femmodel->CalvingRateLevermannx();
+			if(iscalving){
+				switch(calvinglaw){
+					case DefaultCalvingEnum:
+						break;
+					case CalvingLevermannEnum:
+						if(VerboseSolution()) _printf0_("   computing Levermann's calving rate\n");
+						femmodel->StrainRateparallelx();
+						femmodel->StrainRateperpendicularx();
+						femmodel->CalvingRateLevermannx();
+						break;
+					case CalvingPiEnum:
+						if(VerboseSolution()) _printf0_("   computing Pi's calving rate\n");
+						femmodel->StrainRateparallelx();
+						femmodel->DeviatoricStressx();
+						femmodel->CalvingRatePix();
+						break;
+					case CalvingDevEnum:
+						femmodel->CalvingRateDevx();
+						femmodel->ElementOperationx(&Element::CalvingRateDev);
+						break;
+					default:
+						_error_("Caving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
+				}
 			}
-			if(iscalving && calvinglaw==CalvingPiEnum){
-				if(VerboseSolution()) _printf0_("   computing calving rate\n");
-				femmodel->StrainRateparallelx();
-				femmodel->DeviatoricStressx();
-				femmodel->CalvingRatePix();
-			}
-			if(VerboseSolution()) _printf0_("   computing movement of ice boundaries\n");
+			if(VerboseSolution()) _printf0_("   computing levelset transport\n");
 			/* smoothen slope of lsf for computation of normal on ice domain*/
 			levelsetfunctionslope_core(femmodel);
 
Index: ../trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 19032)
+++ ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 19033)
@@ -671,6 +671,15 @@
 
 }
 /*}}}*/
+void FemModel::ElementOperationx(void (Element::*function)(void)){ /*{{{*/
+
+	for(int i=0;i<elements->Size();i++){
+		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+		(element->*function)();
+	}
+
+}
+/*}}}*/
 void FemModel::Responsex(IssmDouble* responses,const char* response_descriptor){/*{{{*/
 
 	int response_descriptor_enum=StringToEnumx(response_descriptor);
@@ -1742,6 +1751,14 @@
 	}
 }
 /*}}}*/
+void FemModel::CalvingRateDevx(){/*{{{*/
+
+	for(int i=0;i<elements->Size();i++){
+		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+		element->CalvingRateDev();
+	}
+}
+/*}}}*/
 void FemModel::StrainRateparallelx(){/*{{{*/
 
 	for(int i=0;i<elements->Size();i++){
Index: ../trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.h	(revision 19032)
+++ ../trunk-jpl/src/c/classes/Elements/Element.h	(revision 19033)
@@ -170,6 +170,7 @@
 		virtual void       AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0;
 		virtual void	    CalvingRateLevermann(void)=0;
 		virtual void       CalvingRatePi(void)=0;
+		virtual void       CalvingRateDev(void){_error_("not implemented yet");};
 		virtual IssmDouble CharacteristicLength(void)=0;
 		virtual void       ComputeBasalStress(Vector<IssmDouble>* sigma_b)=0;
 		virtual void       ComputeDeviatoricStressTensor(void)=0;
Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 19032)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 19033)
@@ -327,6 +327,63 @@
 
 }
 /*}}}*/
+void       Tria::CalvingRateDev(){/*{{{*/
+
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
+	IssmDouble  calvingratex[NUMVERTICES];
+	IssmDouble  calvingratey[NUMVERTICES];
+	IssmDouble  calvingrate[NUMVERTICES];
+	IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
+
+	/* Get node coordinates and dof list: */
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+	/*Retrieve all inputs and parameters we will need*/
+	Input* vx_input=inputs->GetInput(VxEnum);        _assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum);        _assert_(vy_input);
+
+	/* Start looping on the number of vertices: */
+	GaussTria* gauss=new GaussTria();
+	for(int iv=0;iv<NUMVERTICES;iv++){
+		gauss->GaussVertex(iv);
+
+		/*Get velocity components and thickness*/
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		vel=sqrt(vx*vx+vy*vy)+1.e-14;
+
+		/*Compute strain rate and viscosity: */
+		this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+
+		/*Get Eigen values*/
+		Matrix2x2Eigen(&lambda1,&lambda2,&ex,&ey,epsilon[0],epsilon[2],epsilon[1]);
+
+		/*Process Eigen values*/
+		lambda1>0? lambda1 = pow(lambda1,.3) : lambda1=0.;
+		lambda2>0? lambda2 = pow(lambda2,.3) : lambda2=0.;
+		lambda1 = lambda1*5.e-2;
+		lambda2 = lambda2*5.e-2;
+		_assert_(!xIsNan<IssmDouble>(lambda1));
+		_assert_(!xIsNan<IssmDouble>(lambda2));
+
+		/*Assign values*/
+		//calvingratex[iv]=ex*lambda1 - ey*lambda2;
+		//calvingratey[iv]=ey*lambda1 + ex*lambda2;
+		calvingratex[iv]=vx/vel*(lambda1 + lambda2);
+		calvingratey[iv]=vy/vel*(lambda1 + lambda2);
+		calvingrate[iv]=sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
+	}
+
+	/*Add input*/
+	this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
+
+	/*Clean up and return*/
+	delete gauss;
+}
+/*}}}*/
 IssmDouble Tria::CharacteristicLength(void){/*{{{*/
 
 	return sqrt(2*this->GetArea());
Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.h	(revision 19032)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.h	(revision 19033)
@@ -51,13 +51,14 @@
 		/*Element virtual functions definitions: {{{*/
 		void        AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
 		void			CalvingRateLevermann();
+		void			CalvingRatePi();
+		void			CalvingRateDev();
 		IssmDouble  CharacteristicLength(void);
 		void        ComputeBasalStress(Vector<IssmDouble>* sigma_b);
 		void        ComputeDeviatoricStressTensor();
 		void        ComputeSigmaNN();
 		void        ComputeStressTensor();
 		void        ComputeSurfaceNormalVelocity();
-		void			CalvingRatePi();
 		void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
 		void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index);
 		void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
Index: ../trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.h	(revision 19032)
+++ ../trunk-jpl/src/c/classes/FemModel.h	(revision 19033)
@@ -61,6 +61,7 @@
 		void SetCurrentConfiguration(int configuration_type,int analysis_type);
 
 		/*Modules*/ 
+		void ElementOperationx(void (Element::*function)(void));
 		void GetInputLocalMinMaxOnNodesx(IssmDouble** pmin,IssmDouble** pmax,IssmDouble* ug);
 		void MassFluxx(IssmDouble* presponse);
 		void MaxAbsVxx(IssmDouble* presponse);
@@ -87,6 +88,7 @@
 		void DeviatoricStressx();
 		void CalvingRateLevermannx();
 		void CalvingRatePix();
+		void CalvingRateDevx();
 		#ifdef  _HAVE_DAKOTA_
 		void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses);
 		#endif
