Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 21401)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 21402)
@@ -174,4 +174,80 @@
 void       Penta::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){/*{{{*/
 	_error_("Not supported yet!");
+}
+/*}}}*/
+void       Penta::CalvingRateDev(){/*{{{*/
+
+	if(!this->IsOnBase()) return;
+
+	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;
+	IssmDouble  sigma_vm,sigma_max,epse_2,groundedice;
+
+	/* Get node coordinates and dof list: */
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+	/*Depth average B for stress calculation*/
+	this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
+	this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
+	this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
+
+	/*Retrieve all inputs and parameters we will need*/
+	Input* vx_input = inputs->GetInput(VxAverageEnum); _assert_(vx_input);
+	Input* vy_input = inputs->GetInput(VyAverageEnum); _assert_(vy_input);
+	Input* gr_input = inputs->GetInput(MaskGroundediceLevelsetEnum); _assert_(gr_input);
+	IssmDouble  B   = this->GetMaterialParameter(MaterialsRheologyBbarEnum);
+	IssmDouble  n   = this->GetMaterialParameter(MaterialsRheologyNEnum);
+
+	/* Start looping on the number of vertices: */
+	GaussPenta* gauss=new GaussPenta();
+	for(int iv=0;iv<3;iv++){
+		gauss->GaussVertex(iv);
+
+		/*Get velocity components and thickness*/
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		gr_input->GetInputValue(&groundedice,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]);
+		_assert_(!xIsNan<IssmDouble>(lambda1));
+		_assert_(!xIsNan<IssmDouble>(lambda2));
+
+		/*Process Eigen values (only account for extension)*/
+		lambda1 = max(lambda1,0.);
+		lambda2 = max(lambda2,0.);
+
+		/*Calculate sigma_vm*/
+		epse_2    = 1./2. *(lambda1*lambda1 + lambda2*lambda2);
+		sigma_vm  = sqrt(3.) * B * pow(epse_2,1./(2.*n));
+		sigma_max = 1000.e+3;
+
+		if(groundedice<0) sigma_max=200.e+3;
+
+		/*Assign values*/
+		calvingratex[iv]=vx*sigma_vm/sigma_max;
+		calvingratey[iv]=vy*sigma_vm/sigma_max;
+		calvingrate[iv]=sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
+	}
+
+	/*Add input*/
+	this->inputs->AddInput(new PentaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
+
+	this->InputExtrude(CalvingratexEnum,-1);
+	this->InputExtrude(CalvingrateyEnum,-1);
+	this->InputExtrude(CalvingCalvingrateEnum,-1);
+
+	/*Clean up and return*/
+	delete gauss;
 }
 /*}}}*/
@@ -641,5 +717,5 @@
 
 		/*Create gauss point in the middle of the basal edge*/
-		Gauss* gauss=NewGaussBase(1);
+		Gauss* gauss=NewGauss();
 		gauss->GaussPoint(0);
 
@@ -2213,8 +2289,73 @@
 }
 /*}}}*/
+void       Penta::ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble d,xn,yn;
+
+	/*Get current levelset and vertex coordinates*/
+	IssmDouble ls[NUMVERTICES];
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
+	InputDuplicate(MaskIceLevelsetEnum,PressureEnum);
+
+	/*Get distance from list of segments and reset ls*/
+	for(int j=0;j<NUMVERTICES;j++){
+		IssmDouble dmin = 1.e+50;
+		for(int i=0;i<numsegments;i++){
+			IssmDouble x = xyz_list[j][0];
+			IssmDouble y = xyz_list[j][1];
+			IssmDouble l2 = (segments[4*i+2]-segments[4*i+0])*(segments[4*i+2]-segments[4*i+0]) + (segments[4*i+3]-segments[4*i+1])*(segments[4*i+3]-segments[4*i+1]);
+
+			/*Segment has a length of 0*/
+			if(l2==0.){
+				d = (x-segments[4*i+0])*(x-segments[4*i+0])+(y-segments[4*i+1])*(y-segments[4*i+1]);
+				if(d<dmin) dmin = d;
+				continue;
+			}
+
+			/*Consider the line extending the segment, parameterized as v + t (w - v).
+			 *We find projection of point p onto the line.
+			 *It falls where t = [(p-v) . (w-v)] / |w-v|^2*/
+			IssmDouble t = ((x-segments[4*i+0])*(segments[4*i+2]-segments[4*i+0]) + (y-segments[4*i+1])*(segments[4*i+3]-segments[4*i+1]))/l2;
+			if(t < 0.0){
+				// Beyond the 'v' end of the segment
+				d = (x-segments[4*i+0])*(x-segments[4*i+0])+(y-segments[4*i+1])*(y-segments[4*i+1]);
+			}
+			else if (t > 1.0){
+				// Beyond the 'w' end of the segment
+				d = (x-segments[4*i+2])*(x-segments[4*i+2])+(y-segments[4*i+3])*(y-segments[4*i+3]);
+			}
+			else{
+				// Projection falls on the segment
+				xn = segments[4*i+0] + t * (segments[4*i+2] - segments[4*i+0]);
+				yn = segments[4*i+1] + t * (segments[4*i+3] - segments[4*i+1]);
+				d = (x-xn)*(x-xn)+(y-yn)*(y-yn);
+			}
+
+			if(d<dmin) dmin = d;
+		}
+
+		/*Update signed distance*/
+		dmin = sqrt(dmin);
+		if(dmin>10000) dmin=10000;
+		if(ls[j]>0){
+			ls[j] = dmin;
+		}
+		else{
+			ls[j] = - dmin;
+		}
+	}
+
+	/*Update Levelset*/
+	this->inputs->AddInput(new PentaInput(MaskIceLevelsetEnum,&ls[0],P1Enum));
+}
+/*}}}*/
 void       Penta::SetClone(int* minranks){/*{{{*/
 
 	_error_("not implemented yet");
 }
+
 /*}}}*/
 void       Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 21401)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 21402)
@@ -48,4 +48,5 @@
 		void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
 		void           AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
+		void           CalvingRateDev();
 		void           CalvingRateLevermann();
 		IssmDouble     CharacteristicLength(void){_error_("not implemented yet");};
@@ -146,4 +147,5 @@
 		void           ResetFSBasalBoundaryCondition(void);
 		void           ResetHooks();
+		void           ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments);
 		void	         SetClone(int* minranks);
 		void           SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 21401)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 21402)
@@ -1734,4 +1734,8 @@
 	int num_procs = IssmComm::GetSize();
 
+	/*Get domain type (2d or 3d)*/
+	int domaintype;
+	this->parameters->FindParam(&domaintype,DomainTypeEnum);
+	
 	/*1: go throug all elements of this partition and figure out how many
 	 * segments we have (corresopnding to levelset = 0)*/
@@ -1739,5 +1743,8 @@
 	for(int i=0;i<elements->Size();i++){
 		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
-		element->WriteLevelsetSegment(segments);
+		if(!element->IsOnBase()) continue;
+		Element* basalelement = element->SpawnBasalElement();
+		basalelement->WriteLevelsetSegment(segments);
+		if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	}
 
@@ -1771,5 +1778,20 @@
 	for(int i=0;i<elements->Size();i++){
 		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+		if(!element->IsOnBase()) continue;
 		element->ResetLevelsetFromSegmentlist(allsegmentlist,numseg);
+	}
+
+
+	/*Extrude if necessary*/
+	int elementtype;
+	this->parameters->FindParam(&elementtype,MeshElementtypeEnum);
+	if(elementtype==PentaEnum){
+		InputExtrudex(this,MaskIceLevelsetEnum,-1);
+	}
+	else if(elementtype==TriaEnum){
+		/*no need to extrude*/
+	}
+	else{
+		_error_("not implemented yet");
 	}