Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 17249)
+++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 17250)
@@ -66,10 +66,4 @@
 	/*parameters: */
 	bool save_results;
-
-	/* extrapolate */
-	Analysis* analysis = new ExtrapolationAnalysis();
-	femmodel->parameters->SetParam(VxEnum,ExtrapolationVariableEnum);
-	analysis->Core(femmodel);
-	delete analysis;
 
 	/*activate formulation: */
@@ -165,5 +159,5 @@
 				/* Artificial Diffusion */
 				element->ElementSizes(&hx,&hy,&hz);
-				vel=sqrt(vx*vx + vy*vy );
+				vel=sqrt(vx*vx + vy*vy) + 1e-14;
 				h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) ); //FIXME: is this correct?
 
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17249)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17250)
@@ -1044,6 +1044,6 @@
 void StressbalanceAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
 	/*Default, do nothing*/
-	_printf0_("   Updating active and non-active nodes in Femmodel \n");
-	//SetActiveNodesLSMx(femmodel->elements);
+	_printf0_("   Updating active and non-active nodes for StressbalanceAnalysis \n");
+	SetActiveNodesLSMx(femmodel->elements);
 	return;
 }/*}}}*/
Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/transient_core.cpp	(revision 17249)
+++ /issm/trunk-jpl/src/c/cores/transient_core.cpp	(revision 17250)
@@ -150,7 +150,20 @@
 		if(islevelset){
 			if(VerboseSolution()) _printf0_("   computing movement of ice boundaries\n");
+
+			/* extrapolate required variables */
+			Analysis* extanalysis = new ExtrapolationAnalysis();
+			int vars[2] = {VxEnum, VyEnum};
+			for(int iv=0;iv<2;iv++){
+				femmodel->parameters->SetParam(vars[iv],ExtrapolationVariableEnum); 
+				extanalysis->Core(femmodel);
+			}
+			delete extanalysis;	
+
 			analysis = new LevelsetAnalysis();
 			analysis->Core(femmodel);
 			delete analysis;
+
+			/* update vertices included for next calculation */
+			GetMaskOfIceVerticesLSMx(femmodel);
 		}
 
Index: /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp	(revision 17249)
+++ /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp	(revision 17250)
@@ -16,6 +16,28 @@
 		int         numnodes = element->GetNumberOfNodes();
 		IssmDouble *mask     = xNew<IssmDouble>(numnodes);
-		element->GetInputListOnNodes(&mask[0],IceMaskNodeActivationEnum);
-		
+		// include switch for elements with multiple different sets of nodes
+		switch(element->GetElementType()){
+			case MINIEnum:case TaylorHoodEnum:{
+				Input* input=element->GetInput(IceMaskNodeActivationEnum);
+				if(!input) _error_("Input " << EnumToStringx(IceMaskNodeActivationEnum) << " not found in element");
+
+				/* Start looping on the number of vertices: */
+				Gauss* gauss=element->NewGauss();
+				for(int iv=0;iv<element->NumberofNodesVelocity();iv++){
+					gauss->GaussNode(element->VelocityInterpolation(),iv);
+					input->GetInputValue(&mask[iv],gauss);
+				}
+				for(int iv=0;iv<element->NumberofNodesPressure();iv++){
+					gauss->GaussNode(element->PressureInterpolation(),iv);
+					input->GetInputValue(&mask[element->NumberofNodesVelocity()+iv],gauss);
+				}
+				delete gauss;
+				break;
+			}
+			default:
+				element->GetInputListOnNodes(&mask[0],IceMaskNodeActivationEnum);
+				break;
+		}
+
 		for(int in=0;in<numnodes;in++){
 			Node* node=element->GetNode(in);
@@ -33,10 +55,14 @@
 void GetMaskOfIceVerticesLSMx(FemModel* femmodel){/*{{{*/
 
+	/* Intermediaries */
+	int i;
 	/*Initialize vector with number of vertices*/
 	int numvertices=femmodel->vertices->NumberOfVertices();
 	Vector<IssmDouble>* vec_mask_ice=new Vector<IssmDouble>(numvertices); //vertices that have ice at next time step
-
+	for(i=0;i<numvertices;i++){
+		vec_mask_ice->SetValue(i,0.,INS_VAL);
+	}
 	/*Fill vector with values: */
-	for(int i=0;i<femmodel->elements->Size();i++){
+	for(i=0;i<femmodel->elements->Size();i++){
 		Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
 		SetMaskOfIceElement(vec_mask_ice, element);
