Index: issm/trunk/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- issm/trunk/src/c/analyses/LevelsetAnalysis.cpp	(revision 22758)
+++ issm/trunk/src/c/analyses/LevelsetAnalysis.cpp	(revision 23189)
@@ -12,5 +12,7 @@
 
 void LevelsetAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
-	IoModelToConstraintsx(constraints,iomodel,"md.levelset.spclevelset",LevelsetAnalysisEnum,P1Enum);
+	int finiteelement;
+	iomodel->FindConstant(&finiteelement,"md.levelset.fe");
+	IoModelToConstraintsx(constraints,iomodel,"md.levelset.spclevelset",LevelsetAnalysisEnum,finiteelement);
 }
 /*}}}*/
@@ -19,5 +21,6 @@
 }/*}}}*/
 void LevelsetAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
-	int finiteelement=P1Enum;
+	int finiteelement;
+	iomodel->FindConstant(&finiteelement,"md.levelset.fe");
 	if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
 	::CreateNodes(nodes,iomodel,LevelsetAnalysisEnum,finiteelement);
@@ -32,5 +35,6 @@
 
 	/*Finite element type*/
-	int finiteelement = P1Enum;
+	int finiteelement;
+	iomodel->FindConstant(&finiteelement,"md.levelset.fe");
 
 	/*Update elements: */
@@ -43,5 +47,5 @@
 		}
 	}
-	
+
 	iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
 	iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum);
@@ -177,5 +181,5 @@
 
 	/*Calving threshold*/
-	
+
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes    = basalelement->GetNumberOfNodes();
@@ -326,5 +330,5 @@
 				 }
 				break;
-				
+
 			case CalvingLevermannEnum:
 				calvingratex_input->GetInputValue(&c[0],gauss);
@@ -357,5 +361,5 @@
 				 }
 				break;
-			
+
 			case CalvingHabEnum:
 				lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
@@ -378,5 +382,5 @@
 				 }
 				break;
-			
+
 			case CalvingCrevasseDepthEnum:
 				lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
@@ -385,5 +389,5 @@
 
 				if(groundedice<0) meltingrate = 0.;
-				
+
 				norm_dlsf=0.;
 				for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
@@ -516,5 +520,5 @@
 }/*}}}*/
 ElementVector* LevelsetAnalysis::CreatePVector(Element* element){/*{{{*/
-	
+
 	if(!element->IsOnBase()) return NULL;
 	Element* basalelement = element->SpawnBasalElement();
@@ -525,5 +529,5 @@
 	IssmDouble  lsf;
 	IssmDouble* xyz_list = NULL;
-	
+
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = basalelement->GetNumberOfNodes();
@@ -532,5 +536,5 @@
 	ElementVector* pe = basalelement->NewElementVector();
 	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
-	
+
 	if(dt!=0.){
 		/*Initialize basis vector*/
@@ -623,5 +627,5 @@
 	// d=|a x b|/|b|
 	// with a=q-s0, b=s1-s0
-	
+
 	/* Intermediaries */
 	const int dim=2;
@@ -634,5 +638,5 @@
 		b[i]=s1[i]-s0[i];
 	}
-	
+
 	norm_b=0.;
 	for(i=0;i<dim;i++)
@@ -662,57 +666,4 @@
 		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
 	}
-}/*}}}*/
-void           LevelsetAnalysis::SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_signed_dist, Element* element){/*{{{*/
-
-	if(!element->IsZeroLevelset(MaskIceLevelsetEnum))
-		return;
-
-	/* Intermediaries */
-	const int dim=3;
-	int i,d;
-	int numvertices=element->GetNumberOfVertices();
-	IssmDouble s0[dim], s1[dim], v[dim];
-	IssmDouble dist,lsf_old;
-
-	IssmDouble* lsf = xNew<IssmDouble>(numvertices);
-	IssmDouble* sign_lsf = xNew<IssmDouble>(numvertices);
-	IssmDouble* signed_dist = xNew<IssmDouble>(numvertices);
-	IssmDouble* xyz_list = NULL;
-	IssmDouble* xyz_list_zero = NULL;
-
-	/* retrieve inputs and parameters */
-	element->GetVerticesCoordinates(&xyz_list);
-	element->GetInputListOnVertices(lsf,MaskIceLevelsetEnum);
-
-	/* get sign of levelset function */
-	for(i=0;i<numvertices;i++)
-		sign_lsf[i]=(lsf[i]>=0.?1.:-1.);
-
-	element->ZeroLevelsetCoordinates(&xyz_list_zero, xyz_list, MaskIceLevelsetEnum);
-	for(d=0;d<dim;d++){
-		s0[d]=xyz_list_zero[0+d];
-		s1[d]=xyz_list_zero[3+d];
-	}
-
-	/* get signed_distance of vertices to zero levelset straight */
-	for(i=0;i<numvertices;i++){
-		for(d=0;d<dim;d++)
-			v[d]=xyz_list[3*i+d];
-		dist=GetDistanceToStraight(&v[0],&s0[0],&s1[0]);
-		signed_dist[i]=sign_lsf[i]*dist;
-	}
-	
-	/* insert signed_distance into vec_signed_dist, if computed distance is smaller */
-	for(i=0;i<numvertices;i++){
-		vec_signed_dist->GetValue(&lsf_old, element->vertices[i]->Sid());
-		if(xIsNan<IssmDouble>(lsf_old) || fabs(signed_dist[i])<fabs(lsf_old))
-			vec_signed_dist->SetValue(element->vertices[i]->Sid(),signed_dist[i],INS_VAL);
-	}
-
-	xDelete<IssmDouble>(lsf);
-	xDelete<IssmDouble>(sign_lsf);
-	xDelete<IssmDouble>(signed_dist);
-	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(xyz_list_zero);
 }/*}}}*/
 void           LevelsetAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
@@ -724,5 +675,5 @@
 	IssmDouble  bed,water_depth;
 	IssmDouble  levelset;
-	
+
 	femmodel->parameters->FindParam(&calvinglaw,CalvingLawEnum);
 
@@ -756,5 +707,5 @@
 		}
 	}
-	
+
 	if(calvinglaw==CalvingHabEnum){
 
@@ -762,9 +713,9 @@
 		femmodel->elements->InputDuplicate(MaskIceLevelsetEnum, DistanceToCalvingfrontEnum);
 		femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum);
-		
+
 		/*Loop over all elements of this partition*/
 		for(int i=0;i<femmodel->elements->Size();i++){
 			Element* element  = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-			
+
 			rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
 			rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
@@ -797,11 +748,11 @@
 		}
 	}
-	
+
 	if(calvinglaw==CalvingCrevasseDepthEnum){
-	
+
 		int                 nflipped,local_nflipped;
 		Vector<IssmDouble>* vec_constraint_nodes = NULL;
 		IssmDouble* constraint_nodes = NULL;
-		
+
 		/*Get the DistanceToCalvingfront*/
 		femmodel->elements->InputDuplicate(MaskIceLevelsetEnum,DistanceToCalvingfrontEnum);
@@ -895,5 +846,5 @@
 				gauss->GaussNode(element->GetElementType(),in);
 				Node* node=element->GetNode(in);
-				
+
 				if(constraint_nodes[node->Sid()]>0.){
 					node->ApplyConstraint(0,+1.);
