Index: /issm/trunk/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- /issm/trunk/src/c/analyses/LevelsetAnalysis.cpp	(revision 19126)
+++ /issm/trunk/src/c/analyses/LevelsetAnalysis.cpp	(revision 19127)
@@ -11,17 +11,26 @@
 #include "../solutionsequences/solutionsequences.h"
 
-int LevelsetAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
+void LevelsetAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
+	return;
+}
+/*}}}*/
+void LevelsetAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
+	return;
+}/*}}}*/
+void LevelsetAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
+	int finiteelement=P1Enum;
+	if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
+	::CreateNodes(nodes,iomodel,LevelsetAnalysisEnum,finiteelement);
+	iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
+}
+/*}}}*/
+int  LevelsetAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
 	return 1;
 }
 /*}}}*/
-void LevelsetAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
-	return;
-}
-/*}}}*/
 void LevelsetAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
-	int  finiteelement;
 
 	/*Finite element type*/
-	finiteelement = P1Enum;
+	int finiteelement = P1Enum;
 
 	/*Update elements: */
@@ -38,24 +47,43 @@
 	iomodel->FetchDataToInput(elements,VxEnum);
 	iomodel->FetchDataToInput(elements,VyEnum);
-	iomodel->FetchDataToInput(elements,MasstransportCalvingrateEnum);
-}
-/*}}}*/
-void LevelsetAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
-	int finiteelement=P1Enum;
-	if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
-	::CreateNodes(nodes,iomodel,LevelsetAnalysisEnum,finiteelement);
-	iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
-}
-/*}}}*/
-void LevelsetAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
+
+	/*Get calving parameters*/
+	bool iscalving;
+	int  calvinglaw;
+	iomodel->Constant(&iscalving,TransientIscalvingEnum);
+	if(iscalving){
+		iomodel->Constant(&calvinglaw,CalvingLawEnum);
+		iomodel->Constant(&iscalving,TransientIscalvingEnum);
+		switch(calvinglaw){
+			case DefaultCalvingEnum:
+				iomodel->FetchDataToInput(elements,CalvingCalvingrateEnum);
+				iomodel->FetchDataToInput(elements,CalvingMeltingrateEnum);
+				break;
+			case CalvingLevermannEnum:
+				iomodel->FetchDataToInput(elements,CalvinglevermannCoeffEnum);
+				iomodel->FetchDataToInput(elements,CalvinglevermannMeltingrateEnum);
+				break;
+			case CalvingPiEnum:
+				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");
+		}
+	}
+}
+/*}}}*/
+void LevelsetAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
+	parameters->AddObject(iomodel->CopyConstantObject(LevelsetStabilizationEnum));
 	return;
 }
 /*}}}*/
-void LevelsetAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
-	return;
-}/*}}}*/
 
 /*Finite element Analysis*/
-void LevelsetAnalysis::Core(FemModel* femmodel){/*{{{*/
+void           LevelsetAnalysis::Core(FemModel* femmodel){/*{{{*/
 
 	#if !defined(_DEVELOPMENT_)
@@ -64,6 +92,8 @@
 
 	/*parameters: */
+	int  stabilization;
 	bool save_results;
 	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+	femmodel->parameters->FindParam(&stabilization,LevelsetStabilizationEnum);
 
 	/*activate formulation: */
@@ -71,10 +101,15 @@
 
 	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");
-		int outputs = MaskIceLevelsetEnum;
-		femmodel->RequestedOutputsx(&femmodel->results,&outputs,1);
+		int outputs[2] = {MaskIceLevelsetEnum, CalvingCalvingrateEnum};
+		femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
 	}
 }/*}}}*/
@@ -97,11 +132,90 @@
 	_error_("not implemented yet");
 }/*}}}*/
-void LevelsetAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
+void           LevelsetAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
+	 * For node i, Bi can be expressed in the actual coordinate system
+	 * by: 
+	 *       Bi=[ N ]
+	 *          [ N ]
+	 * where N is the finiteelement function for node i.
+	 *
+	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Get nodal functions*/
+	IssmDouble* basis=xNew<IssmDouble>(numnodes);
+	element->NodalFunctions(basis,gauss);
+
+	/*Build B: */
+	for(int i=0;i<numnodes;i++){
+		B[numnodes*0+i] = basis[i];
+		B[numnodes*1+i] = basis[i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(basis);
+}/*}}}*/
+void           LevelsetAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
+	 * For node i, Bi' can be expressed in the actual coordinate system
+	 * by: 
+	 *       Bi_prime=[ dN/dx ]
+	 *                [ dN/dy ]
+	 * where N is the finiteelement function for node i.
+	 *
+	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+	/*Build B': */
+	for(int i=0;i<numnodes;i++){
+		Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
+		Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
+
+}/*}}}*/
+IssmDouble     LevelsetAnalysis::GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1){/*{{{*/
+	// returns distance d of point q to straight going through points s0, s1
+	// d=|a x b|/|b|
+	// with a=q-s0, b=s1-s0
+	
+	/* Intermediaries */
+	const int dim=2;
+	int i;
+	IssmDouble a[dim], b[dim];
+	IssmDouble norm_b;
+
+	for(i=0;i<dim;i++){
+		a[i]=q[i]-s0[i];
+		b[i]=s1[i]-s0[i];
+	}
+	
+	norm_b=0.;
+	for(i=0;i<dim;i++)
+		norm_b+=b[i]*b[i];
+	norm_b=sqrt(norm_b);
+	_assert_(norm_b>0.);
+
+	return fabs(a[0]*b[1]-a[1]*b[0])/norm_b;
+}/*}}}*/
+void           LevelsetAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 	_error_("not implemented yet");
 }/*}}}*/
-void LevelsetAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
+void           LevelsetAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
 	_error_("Not implemented yet");
 }/*}}}*/
-void LevelsetAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
+void           LevelsetAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
 	int domaintype;
@@ -117,64 +231,5 @@
 	}
 }/*}}}*/
-void LevelsetAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
-	/*Default, do nothing*/
-	return;
-}/*}}}*/
-void LevelsetAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
-	 * For node i, Bi can be expressed in the actual coordinate system
-	 * by: 
-	 *       Bi=[ N ]
-	 *          [ N ]
-	 * where N is the finiteelement function for node i.
-	 *
-	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Get nodal functions*/
-	IssmDouble* basis=xNew<IssmDouble>(numnodes);
-	element->NodalFunctions(basis,gauss);
-
-	/*Build B: */
-	for(int i=0;i<numnodes;i++){
-		B[numnodes*0+i] = basis[i];
-		B[numnodes*1+i] = basis[i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(basis);
-}/*}}}*/
-void LevelsetAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
-	 * For node i, Bi' can be expressed in the actual coordinate system
-	 * by: 
-	 *       Bi_prime=[ dN/dx ]
-	 *                [ dN/dy ]
-	 * where N is the finiteelement function for node i.
-	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Get nodal functions derivatives*/
-	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
-	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-	/*Build B': */
-	for(int i=0;i<numnodes;i++){
-		Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
-		Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
-
-}/*}}}*/
-void LevelsetAnalysis::SetDistanceOnIntersectedElements(FemModel* femmodel){/*{{{*/
+void           LevelsetAnalysis::SetDistanceOnIntersectedElements(FemModel* femmodel){/*{{{*/
 
 	/* Intermediaries */
@@ -187,9 +242,9 @@
 
 	Vector<IssmDouble>* vec_dist_zerolevelset = NULL;
-	GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexEnum);
+	GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexPIdEnum);
 	
 	/* set NaN on elements intersected by zero levelset */
 	for(i=0;i<femmodel->elements->Size();i++){
-		element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
 		if(element->IsZeroLevelset(MaskIceLevelsetEnum))
 			for(k=0;k<element->GetNumberOfVertices();k++)
@@ -199,5 +254,5 @@
 	/* set distance on elements intersected by zero levelset */
 	for(i=0;i<femmodel->elements->Size();i++){
-		element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
 		if(element->IsZeroLevelset(MaskIceLevelsetEnum)){
 			SetDistanceToZeroLevelsetElement(vec_dist_zerolevelset, element);
@@ -232,5 +287,5 @@
 	delete dist_zerolevelset;
 }/*}}}*/
-void LevelsetAnalysis::SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_signed_dist, Element* element){/*{{{*/
+void           LevelsetAnalysis::SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_signed_dist, Element* element){/*{{{*/
 
 	if(!element->IsZeroLevelset(MaskIceLevelsetEnum))
@@ -285,26 +340,6 @@
 	xDelete<IssmDouble>(xyz_list_zero);
 }/*}}}*/
-IssmDouble LevelsetAnalysis::GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1){/*{{{*/
-	// returns distance d of point q to straight going through points s0, s1
-	// d=|a x b|/|b|
-	// with a=q-s0, b=s1-s0
-	
-	/* Intermediaries */
-	const int dim=2;
-	int i;
-	IssmDouble a[dim], b[dim];
-	IssmDouble norm_b;
-
-	for(i=0;i<dim;i++){
-		a[i]=q[i]-s0[i];
-		b[i]=s1[i]-s0[i];
-	}
-	
-	norm_b=0.;
-	for(i=0;i<dim;i++)
-		norm_b+=b[i]*b[i];
-	norm_b=sqrt(norm_b);
-	_assert_(norm_b>0.);
-
-	return fabs(a[0]*b[1]-a[1]*b[0])/norm_b;
-}/*}}}*/
+void           LevelsetAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
+	/*Default, do nothing*/
+	return;
+}/*}}}*/