Index: ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 19093)
+++ ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 19094)
@@ -671,8 +671,6 @@
 
 	int         eplflip_lock;
 	int         numnodes      =basalelement->GetNumberOfNodes();
-	IssmDouble  new_active;
-	IssmDouble  recurent;
 	IssmDouble* epl_thickness =xNew<IssmDouble>(numnodes);
 	IssmDouble* old_active    =xNew<IssmDouble>(numnodes);
 	IssmDouble* sedhead       =xNew<IssmDouble>(numnodes);
@@ -699,25 +697,23 @@
 
 	for(i=0;i<numnodes;i++){
 
-		/*Dealing with the initial value to define zigzaging, preceding iteration the first time we see the node and current evolving vector after*/
-		recurence->GetValue(&recurent,basalelement->nodes[i]->Sid());
-		if (recurent>0)vec_mask->GetValue(&old_active[i],basalelement->nodes[i]->Sid());
-		new_active=old_active[i];
-		recurence->SetValue(basalelement->nodes[i]->Sid(),1.,ADD_VAL);
+		/*If node is now closed bring its thickness back to initial*/
+		if (old_active[i]==0.){
+			epl_thickness[i]=init_thick;
+		}
 
 		/*Now starting to look at the activations*/
 		if(residual[i]>0.){
 			vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
-			new_active=1.0;
+			if(old_active[i]==0.)	recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
 		}
 		/*If mask was already one, keep one or colapse*/
 		else if(old_active[i]>0.){
 			vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
-			new_active=1.0;
 			/*If epl thickness gets under colapse thickness, close the layer*/
 			if(epl_thickness[i]<colapse_thick){
 				vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
-				new_active=0.0;
+				recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
 			}
 		}
 		/*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/
@@ -727,22 +723,10 @@
 				/*Increase of the domain is on the downstream node in term of sediment head*/
 				if(sedhead[j] == sedheadmin){
 					vec_mask->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL);
+					if(old_active[i]==0.)	recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
 				}
 			}
 		}
-		/*Now checking for nodes zigzaging between open and close*/
-		if (old_active[i] != new_active){
-			eplzigzag_counter[basalelement->nodes[i]->Lid()] ++;
-			/*If node changed too much of state, fix it to it's last known state*/
-			if(eplzigzag_counter[basalelement->nodes[i]->Lid()]>eplflip_lock & eplflip_lock!=0){
-				vec_mask->SetValue(basalelement->nodes[i]->Sid(),old_active[i],INS_VAL);
-				new_active=old_active[i];
-			}
-		}
-			/*If node is now closed bring its thickness back to initial*/
-		if (new_active==0.){
-			epl_thickness[i]=init_thick;
-		}
 	}
 	basalelement->AddInput(HydrologydcEplThicknessEnum,epl_thickness,basalelement->GetElementType());
 
Index: ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 19093)
+++ ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 19094)
@@ -555,7 +555,7 @@
 	Element* element;
 
 	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++){
Index: ../trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.cpp	(revision 19093)
+++ ../trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.cpp	(revision 19094)
@@ -416,7 +416,7 @@
 	int numvertices=femmodel->vertices->NumberOfVertices();
 
 	Vector<IssmDouble>* vec_dist_zerolevelset = NULL;
-	GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexEnum);
+	GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexPIdEnum);
 	
 	/* set distance on elements intersected by zero levelset */
 	for(i=0;i<femmodel->elements->Size();i++){
Index: ../trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp	(revision 19093)
+++ ../trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp	(revision 19094)
@@ -141,7 +141,7 @@
 	vec_phi=new Vector<IssmDouble>(vertices->NumberOfVertices());
 	for(i=0;i<elements->Size();i++){
 		Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-		element->GetVectorFromInputs(vec_phi,MaskGroundediceLevelsetEnum);
+		element->GetVectorFromInputs(vec_phi,MaskGroundediceLevelsetEnum,VertexPIdEnum);
 	}
 	vec_phi->Assemble();
 	phi=vec_phi->ToMPISerial();
Index: ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp	(revision 19093)
+++ ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp	(revision 19094)
@@ -6,24 +6,25 @@
 #include "../../shared/shared.h"
 #include "../../toolkits/toolkits.h"
 
-void GetVectorFromInputsx( Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type){
+void GetVectorFromInputsx(Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type){
 
 	int i;
 	Vector<IssmDouble>* vector=NULL;
 
-	if(type==VertexEnum){
-
-		/*Allocate vector*/
+	switch(type){
+	case VertexPIdEnum: case VertexSIdEnum:
 		vector=new Vector<IssmDouble>(femmodel->vertices->NumberOfVertices());
-
-		/*Look up in elements*/
-		for(i=0;i<femmodel->elements->Size();i++){
-			Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-			element->GetVectorFromInputs(vector,name);
-		}
+		break;
+	case NodesEnum:case NodeSIdEnum:
+		vector=new Vector<IssmDouble>(femmodel->nodes->NumberOfNodes());
+		break;
+	default:
+			_error_("vector type: " << EnumToStringx(type) << " not supported yet!");
 	}
-	else{
-		_error_("vector type: " << EnumToStringx(type) << " not supported yet!");
+	/*Look up in elements*/
+	for(i=0;i<femmodel->elements->Size();i++){
+		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		element->GetVectorFromInputs(vector,name,type);
 	}
 
 	vector->Assemble();
Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_la.cpp
===================================================================
--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_la.cpp	(revision 19093)
+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_la.cpp	(revision 19094)
@@ -43,8 +43,8 @@
 	Vector<IssmDouble>* vel           = NULL;
 	Vector<IssmDouble>* vel_old       = NULL;
 	Vector<IssmDouble>* vel_old_local = NULL;
-	GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexEnum);
-	GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexEnum);
+	GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexPIdEnum);
+	GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexPIdEnum);
 
 	while(true){
 		count++;
@@ -69,7 +69,7 @@
 
 			/*Update solution*/
 			InputUpdateFromSolutionx(femmodel,ug); delete ug;
-			GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexEnum);
+			GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexPIdEnum);
 			/*Check for convergence*/
 			Vector<IssmDouble>* dvel_local=vel_old_local->Duplicate(); vel_old_local->Copy(dvel_local); dvel_local->AYPX(vel,-1.0);
 			IssmDouble ndu=dvel_local->Norm(NORM_TWO);   delete dvel_local;
@@ -100,7 +100,7 @@
 
 		/*Update solution*/
 		InputUpdateFromSolutionx(femmodel,pug); delete pug;
-		GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexEnum);
+		GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexPIdEnum);
 
 		/*Check for convergence*/
 		Vector<IssmDouble>* dvel=vel_old_local->Duplicate(); vel_old_local->Copy(dvel); dvel->AYPX(vel,-1.0);
Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp
===================================================================
--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp	(revision 19093)
+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp	(revision 19094)
@@ -40,7 +40,7 @@
 	GetSolutionFromInputsx(&ug,femmodel);
 	Vector<IssmDouble>* vx     = NULL;
 	Vector<IssmDouble>* vx_old = NULL;
-	GetVectorFromInputsx(&vx,femmodel,VxEnum,VertexEnum);
+	GetVectorFromInputsx(&vx,femmodel,VxEnum,VertexPIdEnum);
 
 	while(true){
 		count++;
@@ -68,7 +68,7 @@
 		/*Update d and tau accordingly*/
 		analysis->InputUpdateFromSolutionFSXTH_d(  femmodel->elements,femmodel->parameters);
 		analysis->InputUpdateFromSolutionFSXTH_tau(femmodel->elements,femmodel->parameters);
-		GetVectorFromInputsx(&vx,femmodel,VxEnum,VertexEnum);
+		GetVectorFromInputsx(&vx,femmodel,VxEnum,VertexPIdEnum);
 
 		/*Check for convergence*/
 		//Vector<IssmDouble>* dug=ug_old->Duplicate(); ug_old->Copy(dug); dug->AYPX(ug,-1.0);
Index: ../trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 19093)
+++ ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 19094)
@@ -610,8 +610,8 @@
 	IssmDouble         *bed     = NULL;
 
 	/*get vertex vectors for bed and thickness: */
-	GetVectorFromInputsx(&surface  ,this, SurfaceEnum,VertexEnum);
-	GetVectorFromInputsx(&bed      ,this, BaseEnum,    VertexEnum);
+	GetVectorFromInputsx(&surface  ,this, SurfaceEnum,VertexPIdEnum);
+	GetVectorFromInputsx(&bed      ,this, BaseEnum,   VertexPIdEnum);
 
 	/*Allocate vector*/
 	Vector<IssmDouble> *vx=new Vector<IssmDouble>(vertices->NumberOfVertices());
@@ -1830,7 +1830,7 @@
 		if(flag==ScaledEnum){
 
 			/*this response was scaled. pick up the response from the inputs: */
-			GetVectorFromInputsx(&vertex_response,this, StringToEnumx(root),VertexEnum);
+			GetVectorFromInputsx(&vertex_response,this, StringToEnumx(root),VertexPIdEnum);
 
 			/*Now, average it onto the partition nodes: */
 			AverageOntoPartitionx(&qmu_response,elements,nodes,vertices,loads,materials,parameters,vertex_response);
@@ -1913,25 +1913,45 @@
 
 	Vector<IssmDouble>* mask							= NULL;
 	Vector<IssmDouble>* recurence  				= NULL;
+	Vector<IssmDouble>* active						= NULL;
 	IssmDouble*         serial_mask				= NULL;
-	Vector<IssmDouble>* active						= NULL;
+	IssmDouble*         serial_rec  			= NULL;
 	IssmDouble*         serial_active			= NULL;
+	IssmDouble*         old_active        = NULL;
 	int*                eplzigzag_counter =	NULL;
-		
+	int                 eplflip_lock;
+	
 	HydrologyDCEfficientAnalysis* effanalysis =  new HydrologyDCEfficientAnalysis();
 	HydrologyDCInefficientAnalysis* inefanalysis =  new HydrologyDCInefficientAnalysis();
+
 	/*Step 1: update mask, the mask might be extended by residual and/or using downstream sediment head*/
 	mask=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
 	recurence=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
 	this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum); 
+	this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum); 
+	GetVectorFromInputsx(&old_active,this,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum);
+	
 	for (int i=0;i<elements->Size();i++){
 		Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
 		effanalysis->HydrologyEPLGetMask(mask,recurence,eplzigzag_counter,element);
 	}
+	/*check for changes and increment zigzag counter, change the mask if necessary*/
+	recurence->Assemble();
+	serial_rec=recurence->ToMPISerial();
+	for (int i=0;i<nodes->Size();i++){
+		Node* node=xDynamicCast<Node*>(nodes->GetObjectByOffset(i));
+		if(serial_rec[node->Sid()]==1.)eplzigzag_counter[node->Lid()] ++;
+		if(eplzigzag_counter[node->Lid()]>eplflip_lock & eplflip_lock!=0){
+			mask->SetValue(node->Sid(),old_active[node->Sid()],INS_VAL);
+		}
+	}
+
+	
 	this->parameters->SetParam(eplzigzag_counter,this->nodes->Size(),EplZigZagCounterEnum);
 	/*Assemble and serialize*/
 	mask->Assemble();
-	serial_mask=mask->ToMPISerial();
+	serial_mask=mask->ToMPISerial();	
+	
 	xDelete<int>(eplzigzag_counter);
 	delete mask;
 	delete recurence;
Index: ../trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.h	(revision 19093)
+++ ../trunk-jpl/src/c/classes/Elements/Element.h	(revision 19094)
@@ -95,7 +95,7 @@
 		void               GetNodesLidList(int* lidlist);
 		void               GetNodesSidList(int* sidlist);
 		void               GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);
-		void               GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum);
+		void               GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum, int type);
 		void	             GetVertexPidList(int* pidlist);
 		void               GetVerticesConnectivityList(int* connectivitylist);
 		void               GetVerticesCoordinates(IssmDouble** xyz_list);
Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 19093)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 19094)
@@ -1942,7 +1942,6 @@
 		break;
 
 	case NodesEnum:
-
 		/*Get number of nodes and dof list: */
 		numnodes = this->NumberofNodes(this->element_type);
 		values   = xNew<IssmDouble>(numnodes);
@@ -1956,7 +1955,6 @@
 		break;
 
 	case NodeSIdEnum:
-
 		/*Get number of nodes and dof list: */
 		numnodes = this->NumberofNodes(this->element_type);
 		values   = xNew<IssmDouble>(numnodes);
Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 19093)
+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 19094)
@@ -861,20 +861,71 @@
 	*phi=4.*epsilon_eff*epsilon_eff*viscosity;
 }
 /*}}}*/
-void       Element::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){/*{{{*/
+/* void       Element::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){/\*{{{*\/ */
 
+/* 	/\*Fetch number vertices for this element and allocate arrays*\/ */
+/* 	int numvertices = this->GetNumberOfVertices(); */
+/* 	int*        vertexpidlist = xNew<int>(numvertices); */
+/* 	IssmDouble* values        = xNew<IssmDouble>(numvertices); */
+
+/* 	/\*Fill in values*\/ */
+/* 	this->GetVertexPidList(vertexpidlist); */
+/* 	this->GetInputListOnVertices(values,input_enum); */
+/* 	vector->SetValues(numvertices,vertexpidlist,values,INS_VAL); */
+
+/* 	/\*Clean up*\/ */
+/* 	xDelete<int>(vertexpidlist); */
+/* 	xDelete<IssmDouble>(values); */
+
+/* } */
+/* /\*}}}*\/ */
+void       Element::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum,int type){/*{{{*/
+
 	/*Fetch number vertices for this element and allocate arrays*/
-	int numvertices = this->GetNumberOfVertices();
-	int*        vertexpidlist = xNew<int>(numvertices);
-	IssmDouble* values        = xNew<IssmDouble>(numvertices);
+	int         numvertices = this->GetNumberOfVertices();
+	int         numnodes    = this->GetNumberOfNodes();
+	int*        doflist     = NULL;
+	IssmDouble* values      = NULL;
 
-	/*Fill in values*/
-	this->GetVertexPidList(vertexpidlist);
-	this->GetInputListOnVertices(values,input_enum);
-	vector->SetValues(numvertices,vertexpidlist,values,INS_VAL);
-
+	switch(type){
+	case VertexPIdEnum:
+		doflist = xNew<int>(numvertices);
+		values = xNew<IssmDouble>(numvertices);
+		/*Fill in values*/
+		this->GetVertexPidList(doflist);
+		this->GetInputListOnVertices(values,input_enum);
+		vector->SetValues(numvertices,doflist,values,INS_VAL);
+		break;
+	case VertexSIdEnum:
+		doflist = xNew<int>(numvertices);
+		values = xNew<IssmDouble>(numvertices);
+		/*Fill in values*/
+		this->GetVerticesSidList(doflist);
+		this->GetInputListOnVertices(values,input_enum);
+		vector->SetValues(numvertices,doflist,values,INS_VAL);
+		break;
+	case NodesEnum:
+		doflist = xNew<int>(numnodes);
+		values = xNew<IssmDouble>(numnodes);
+		/*Fill in values*/
+		this->GetInputListOnNodes(values,input_enum);
+		this->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+		vector->SetValues(numnodes,doflist,values,INS_VAL);
+		break;
+	case NodeSIdEnum:
+		doflist = xNew<int>(numnodes);
+		values = xNew<IssmDouble>(numnodes);
+		/*Fill in values*/
+		this->GetNodesSidList(doflist);
+		this->GetInputListOnNodes(values,input_enum);
+		vector->SetValues(numnodes,doflist,values,INS_VAL);
+		break;
+	default:
+		_error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
+	}
+	
 	/*Clean up*/
-	xDelete<int>(vertexpidlist);
+	xDelete<int>(doflist);
 	xDelete<IssmDouble>(values);
 
 }
