Index: /issm/trunk-jpl/externalpackages/vim/addons/vimrc
===================================================================
--- /issm/trunk-jpl/externalpackages/vim/addons/vimrc	(revision 17226)
+++ /issm/trunk-jpl/externalpackages/vim/addons/vimrc	(revision 17227)
@@ -100,9 +100,13 @@
 
 " save & "make" the current file in all modes
-map <F8> :w <Enter> :make <Enter><Enter>
-map! <F8>  <ESC> :w <Enter> :make <Enter><Enter>
+map <F6> :w <Enter> :make <Enter><Enter>
+map! <F6>  <ESC> :w <Enter> :make <Enter><Enter>
 " make update: nice for longer documents
-map <F7> :w <Enter> :make update <Enter><Enter>
-map! <F7> <ESC> :w <Enter> :make update <Enter><Enter>
+map <F5> :w <Enter> :make update <Enter><Enter>
+map! <F5> <ESC> :w <Enter> :make update <Enter><Enter>
+
+" switch between tabs
+nmap <F7> :tabp <CR>
+nmap <F8> :tabn <CR>
 
 
@@ -141,4 +145,13 @@
   autocmd BufWritePost   *.sh         !chmod +x %
 
+	" Commenting blocks of code.
+	autocmd FileType c,cpp,java,scala let b:comment_leader = '// '
+	autocmd FileType sh,ruby,python   let b:comment_leader = '# '
+	autocmd FileType conf,fstab       let b:comment_leader = '# '
+	autocmd FileType tex              let b:comment_leader = '% '
+	autocmd FileType mail             let b:comment_leader = '> '
+	autocmd FileType vim              let b:comment_leader = '" '
+	noremap <silent> ,cc :<C-B>silent <C-E>s/^/<C-R>=escape(b:comment_leader,'\/')<CR>/<CR>:nohlsearch<CR>
+	noremap <silent> ,cu :<C-B>silent <C-E>s/^\V<C-R>=escape(b:comment_leader,'\/')<CR>//e<CR>:nohlsearch<CR>
 endif " has("autocmd")
 " ----------------------------------------------------------------------}}}
Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 17226)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 17227)
@@ -285,4 +285,6 @@
 					./modules/CreateNodalConstraintsx/CreateNodalConstraintsx.h\
 					./modules/CreateNodalConstraintsx/CreateNodalConstraintsx.cpp\
+					./modules/SetActiveNodesLSMx/SetActiveNodesLSMx.h\
+					./modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp\
 					./modules/UpdateDynamicConstraintsx/UpdateDynamicConstraintsx.h\
 					./modules/UpdateDynamicConstraintsx/UpdateDynamicConstraintsx.cpp\
Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 17226)
+++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 17227)
@@ -33,5 +33,5 @@
 	iomodel->FetchDataToInput(elements,VxEnum);
 	iomodel->FetchDataToInput(elements,VyEnum);
-	
+	iomodel->FetchDataToInput(elements,IceMaskNodeActivationEnum);
 }
 /*}}}*/
@@ -67,19 +67,12 @@
 	bool save_results;
 
-	/* extrapolate */
-	Analysis* analysis = new ExtrapolationAnalysis();
-	femmodel->parameters->SetParam(VxEnum,ExtrapolationVariableEnum);
-	analysis->Core(femmodel);
-	delete analysis;
-
 	/*activate formulation: */
 	femmodel->SetCurrentConfiguration(LevelsetAnalysisEnum);
 
+	if(VerboseSolution()) _printf0_("call computational core:\n");
+	solutionsequence_linear(femmodel);
+
 	/*recover parameters: */
 	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
-
-	if(VerboseSolution()) _printf0_("call computational core:\n");
-	solutionsequence_linear(femmodel);
-
 	if(save_results){
 		if(VerboseSolution()) _printf0_("   saving results\n");
@@ -328,2 +321,56 @@
 }/*}}}*/
 
+/* Update of constraints */
+// void LevelsetAnalysis::UpdateNoIceConstraints(FemModel* femmodel){/*{{{*/
+// 
+// 	IssmDouble* mask_ice = GetMaskOfIce(femmodel->elements, femmodel->nodes); 
+// 
+// 	for (int i=0;i<femmodel->nodes->Size();i++){
+// 		Node* node=dynamic_cast<Node*>(femmodel->nodes->GetObjectByOffset(i));
+// 		if(node->InAnalysis(LevelsetAnalysisEnum)){
+// 			if(mask_ice[node->Sid()]==1.){//FIXME: what should be done with actual spcs to ice model?
+// // 				node->DofInFSet(0); /*remove spc*/ 
+// 			}
+// 			else{
+// 				IssmDouble defval=0.;
+// 				node->ApplyConstraint(1,defval); /*apply spc*/ 
+// 			}
+// 		}
+// 	}
+// }/*}}}*/
+// IssmDouble* LevelsetAnalysis::GetMaskOfIce(Elements* elements, Nodes* nodes){/*{{{*/
+// 
+// 	int                 i;
+// 	IssmDouble*         mask_ice      = NULL;
+// 	Vector<IssmDouble>* vec_mask_ice  = NULL;
+// 	Element*            element       = NULL;
+// 
+// 	/*Initialize vector with number of vertices*/
+// 	IssmDouble numnodes=nodes->NumberOfNodes(LevelsetAnalysisEnum);
+// 	vec_mask_ice=new Vector<IssmDouble>(numnodes); //nodes at ice front that have ice at next time step
+// 	for(i=0;i<numnodes;i++)
+// 		vec_mask_ice[i]=0.;
+// 	/*Fill vector vertices that have no contact to ice: */
+// 	for(i=0;i<elements->Size();i++){
+// 		element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
+// 		SetMaskOfIceElement(vec_mask_ice, element);
+// 	}
+// 
+// 	/*Assemble vector and serialize */
+// 	vec_mask_ice->Assemble();
+// 	mask_ice=vec_mask_ice->ToMPISerial();
+// 	delete vec_mask_ice;
+// 	return mask_ice;
+// 
+// }/*}}}*/
+// void LevelsetAnalysis::SetMaskOfIceElement(Vector<IssmDouble>* vec_mask_ice, Element* element){/*{{{*/
+// 
+// 	/* Intermediaries */
+// 	int numnodes = element->GetNumberOfNodes();
+// 	
+// 	if(element->IsIceInElement()){
+// 		for(int i = 0;i<numnodes;i++){
+// 			vec_mask_ice->SetValue(element->nodes[i]->Sid(),1.,INS_VAL);
+// 		}
+// 	}
+// }/*}}}*/
Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h	(revision 17226)
+++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h	(revision 17227)
@@ -32,4 +32,8 @@
 	void GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
 
+	/* Updating constraints */
+// 	void UpdateNoIceConstraints(FemModel* femmodel);
+// 	IssmDouble* GetMaskOfIce(Elements* elements, Nodes* nodes);
+// 	void SetMaskOfIceElement(Vector<IssmDouble>* vec_mask_ice, Element* element);
 };
 #endif
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17226)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17227)
@@ -189,4 +189,6 @@
 	iomodel->FetchDataToInput(elements,LoadingforceYEnum);
 	iomodel->FetchDataToInput(elements,DamageDEnum);
+	iomodel->FetchDataToInput(elements,IceMaskNodeActivationEnum);
+
 
 	if(iomodel->meshtype==Mesh3DEnum){
@@ -1024,5 +1026,5 @@
 }/*}}}*/
 void StressbalanceAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
-	/*Default, do nothing*/
+	SetActiveNodesLSMx(femmodel->elements);
 	return;
 }/*}}}*/
@@ -1143,4 +1145,8 @@
 	delete Ke1;
 	delete Ke2;
+	
+	// FIXME: TESTING
+// 	_printf0_("\n element ID: " << element->Id() << ";\n");
+// 	Ke->Echo();
 	return Ke;
 }/*}}}*/
@@ -1265,4 +1271,7 @@
 		newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
 		D_scalar=2.*newviscosity*thickness*gauss->weight*Jdet;
+		
+// 		_printf0_("element, node, weight, Jdet: " << element->Id() << "; " << ig << "; " << gauss->weight << "; " << Jdet << "\n");
+
 		for(int i=0;i<3;i++) D[i*3+i]=D_scalar;
 
@@ -1315,4 +1324,8 @@
 	delete pe1;
 	delete pe2;
+
+// 	_printf0_("\n element ID: " << element->Id() << ";\n");
+// 	pe->Echo();
+
 	return pe;
 }/*}}}*/
@@ -1372,5 +1385,5 @@
 
 	/*If no front, return NULL*/
-	if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
+	if(!element->IsIcefront()) return NULL;
 
 	/*Intermediaries*/
@@ -1395,6 +1408,10 @@
 	IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
 	element->GetVerticesCoordinates(&xyz_list);
-	element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
+
+// 	element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
+ 	element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
+
 	element->NormalSection(&normal[0],xyz_list_front);
+// 	element->GetNormalFromLSF(&normal[0]);
 
 	/*Start looping on Gaussian points*/
@@ -1785,5 +1802,5 @@
 
 	/*If no front, return NULL*/
-	if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
+	if(!element->IsIcefront()) return NULL;
 
 	/*Intermediaries*/
@@ -2208,5 +2225,5 @@
 
 	/*If no front, return NULL*/
-	if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
+	if(!element->IsIcefront()) return NULL;
 
 	/*Intermediaries*/
@@ -2974,5 +2991,5 @@
 
 	/*If no front, return NULL*/
-	if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
+	if(!element->IsIcefront()) return NULL;
 
 	/*Intermediaries*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 17226)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 17227)
@@ -220,4 +220,5 @@
 		virtual int    PressureInterpolation()=0;
 		virtual bool   IsZeroLevelset(int levelset_enum)=0;
+		virtual bool   IsIcefront(void)=0;
 		virtual void   ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum)=0;
 		virtual void   GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 17226)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 17227)
@@ -114,4 +114,5 @@
 		int    PressureInterpolation();
 		bool   IsZeroLevelset(int levelset_enum);
+		bool   IsIcefront(void){_error_("not implemented yet");};
 		void   ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
 		void   GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 17226)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 17227)
@@ -159,4 +159,5 @@
 		void        ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");};
 		bool        IsZeroLevelset(int levelset_enum){_error_("not implemented");};
+		bool		IsIcefront(void){_error_("not implemented yet");};
 		void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented");};
 		void		GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 17226)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 17227)
@@ -909,12 +909,28 @@
 		}
 	}
-
+	_assert_(nrfrontnodes==2);
+	
+	/* arrange order of frontnodes such that they are oriented counterclockwise */
+	if((NUMVERTICES+indicesfront[0]-indicesfront[1])%NUMVERTICES!=2){
+		int index=indicesfront[0];
+		indicesfront[0]=indicesfront[1];
+		indicesfront[1]=index;
+	}
 	IssmDouble* xyz_front = xNew<IssmDouble>(3*nrfrontnodes);
 	/* Return nodes */
 	for(i=0;i<nrfrontnodes;i++){
 		for(dir=0;dir<3;dir++){
-			xyz_front[3*i+dir]=xyz_list[indicesfront[i]+dir];
-		}
-	}
+			xyz_front[3*i+dir]=xyz_list[3*indicesfront[i]+dir];
+		}
+	}
+
+// 	for(i=0;i<nrfrontnodes;i++){
+// 		_printf0_("coords frontnode " << i << " :[");
+// 		for(dir=0;dir<3;dir++){
+// 			xyz_front[3*i+dir]=xyz_list[indicesfront[i]+dir];
+// 			_printf0_(xyz_front[3*i+dir] << "; ");
+// 		}
+// 		_printf0_("]\n");
+// 	}
 
 	*pxyz_front=xyz_front;
@@ -922,5 +938,5 @@
 	xDelete<int>(indicesfront);
 }/*}}}*/
-void  Tria::GetNormalFromLSF(IssmDouble *pnormal){/*{{{*/
+void  Tria::GetNormalFromLSF(IssmDouble *normal){/*{{{*/
 
 	/* Intermediaries */
@@ -929,5 +945,4 @@
 	IssmDouble* xyz_list = NULL;
 	IssmDouble  dlevelset[dim], norm_dlevelset;
-	IssmDouble  normal[dim]={0.};
 
 	/*Retrieve all inputs and parameters*/
@@ -936,4 +951,6 @@
 	
 	counter=0;
+	for(i=0;i<dim;i++) normal[i]=0.;
+
 	Gauss* gauss = this->NewGauss(2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
@@ -947,7 +964,5 @@
 	}
 	_assert_(counter>0);
-	for(i=0;i<dim;i++) normal[i]/counter;
-	
-	pnormal=&normal[0];
+	for(i=0;i<dim;i++) normal[i]/=IssmDouble(counter);
 
 	delete gauss;
@@ -1610,4 +1625,5 @@
 				name==MaskGroundediceLevelsetEnum ||
 				name==MaskIceLevelsetEnum ||
+				name==IceMaskNodeActivationEnum ||
 				name==SurfaceSlopeXEnum ||
 				name==SurfaceSlopeYEnum ||
@@ -2586,5 +2602,5 @@
 	GetInputListOnVertices(&ls[0],levelset_enum);
 
-	/*If the level set is awlays <0, there is no ice front here*/
+	/* If levelset function changes sign, or is zero on at least two vertices, there is a zero level set here */
 	iszerols= false;
 	if(IsIceInElement()){
@@ -2597,4 +2613,26 @@
 }
 /*}}}*/
+/*FUNCTION Tria::IsIcefront{{{*/
+bool Tria::IsIcefront(void){
+
+	bool isicefront;
+	int i,nrice;
+	IssmDouble ls[NUMVERTICES];
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
+
+	/* If only one vertex has ice, there is an ice front here */
+	isicefront= false;
+	if(IsIceInElement()){
+		nrice=0;	
+		for(i=0;i<NUMVERTICES;i++)
+			if(ls[i]<0.) nrice++;
+		if(nrice==1) isicefront= true;
+	}
+	return isicefront;
+}
+/*}}}*/
+
 
 #ifdef _HAVE_RESPONSES_
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 17226)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 17227)
@@ -132,6 +132,7 @@
 		void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
 		void	    GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
-	    void        GetNormalFromLSF(IssmDouble *pnormal);
+		void        GetNormalFromLSF(IssmDouble *normal);
 		bool        IsZeroLevelset(int levelset_enum);
+		bool        IsIcefront(void);
 
 		#ifdef _HAVE_RESPONSES_
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 17226)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 17227)
@@ -406,7 +406,10 @@
 	/*start module: */
 	if(VerboseModule()) _printf0_("   Updating constraints for time: " << time << "\n");
-
-	// analysis->UpdateConstraints();
-	
+	if(VerboseModule()) _printf0_("   and analysis: " << EnumToStringx(analysis_type) << "\n");
+	Analysis* analysis = EnumToAnalysis(analysis_type);
+	_assert_(analysis);
+// 	analysis->UpdateConstraints(this);
+	delete analysis;
+
 	/*Second, constraints might be time dependent: */
 	SpcNodesx(nodes,constraints,parameters,analysis_type); 
Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/transient_core.cpp	(revision 17226)
+++ /issm/trunk-jpl/src/c/cores/transient_core.cpp	(revision 17227)
@@ -162,7 +162,21 @@
 		if(islevelset){
 			if(VerboseSolution()) _printf0_("   computing movement of ice boundaries\n");
-			analysis = new LevelsetAnalysis();
-			analysis->Core(femmodel);
-			delete analysis;
+
+			/* extrapolate */
+			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;
+
+			LevelsetAnalysis* lsanalysis = new LevelsetAnalysis();
+			/* solve level-set equation */
+			lsanalysis->Core(femmodel);
+			delete lsanalysis;
+
+			/* update vertices included for next calculation */
+			GetMaskOfIceVerticesLSMx(femmodel);
 		}
 
Index: /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 17226)
+++ /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 17227)
@@ -119,4 +119,17 @@
 	//Kfs->AllocationInfo();
 
+	// TESTING
+	IssmDouble* Kffd=Kff->ToSerial();
+	Kff->GetSize(&M,&N);
+	_printf0_(" (Kff stiffness matrix size: "<<M<<" x "<<N<<")\n");
+// 	printarray(Kffd,M,N);
+   for(int row=0;row<M;row++){
+		IssmDouble value = Kffd[row*N + row];
+		if(fabs(value)<1.)
+		 _printf0_("Kff(" << row << "," << row << ")=" << value << "\n");
+	 }
+	delete Kffd;
+	
+
 	/*cleanu up and assign output pointers: */
 	delete analysis;
Index: /issm/trunk-jpl/src/c/modules/modules.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/modules.h	(revision 17226)
+++ /issm/trunk-jpl/src/c/modules/modules.h	(revision 17227)
@@ -31,4 +31,5 @@
 #include "./GetVectorFromControlInputsx/GetVectorFromControlInputsx.h"
 #include "./GiaDeflectionCorex/GiaDeflectionCorex.h"
+#include "./SetActiveNodesLSMx/SetActiveNodesLSMx.h"
 #include "./SetControlInputsFromVectorx/SetControlInputsFromVectorx.h"
 #include "./Gradjx/Gradjx.h"
Index: /issm/trunk-jpl/src/m/dev/devpath.py
===================================================================
--- /issm/trunk-jpl/src/m/dev/devpath.py	(revision 17226)
+++ /issm/trunk-jpl/src/m/dev/devpath.py	(revision 17227)
@@ -25,9 +25,9 @@
 from plotmodel import plotmodel
 
-#c = get_ipython().config
-#c.InteractiveShellApp.exec_lines = []
-#c.InteractiveShellApp.exec_lines.append('%load_ext autoreload')
-#c.InteractiveShellApp.exec_lines.append('%autoreload 2')
-#c.InteractiveShellApp.exec_lines.append('print "Warning: disable autoreload in startup.py to improve performance." ')
+c = get_ipython().config
+c.InteractiveShellApp.exec_lines = []
+c.InteractiveShellApp.exec_lines.append('%load_ext autoreload')
+c.InteractiveShellApp.exec_lines.append('%autoreload 2')
+c.InteractiveShellApp.exec_lines.append('print "Warning: disable autoreload in startup.py to improve performance." ')
 
 print("\n  ISSM development path correctly loaded\n\n")
Index: /issm/trunk-jpl/src/m/plot/processmesh.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/processmesh.py	(revision 17226)
+++ /issm/trunk-jpl/src/m/plot/processmesh.py	(revision 17227)
@@ -14,6 +14,6 @@
 	if md.mesh.numberofvertices==0:
 		raise ValueError('processmesh error: mesh is empty')
-	if md.mesh.numberofvertices==md.mesh.numberofelements:
-		raise ValueError('processmesh error: the number of elements is the same as the number of nodes')
+# 	if md.mesh.numberofvertices==md.mesh.numberofelements:
+# 		raise ValueError('processmesh error: the number of elements is the same as the number of nodes')
 
 	if len(data)==0 or not isinstance(data,dict):
Index: /issm/trunk-jpl/test/NightlyRun/test336.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test336.py	(revision 17226)
+++ /issm/trunk-jpl/test/NightlyRun/test336.py	(revision 17227)
@@ -8,5 +8,5 @@
 from MatlabFuncs import *
 
-md=triangle(model(),'../Exp/Square.exp',150000.)
+md=triangle(model(),'../Exp/Square.exp',10000.)
 md=setmask(md,'','')
 md=parameterize(md,'../Par/SquareSheetConstrained.py') 
@@ -20,5 +20,5 @@
 md.transient.isgroundingline=False
 md.transient.isgia=False
-md.transient.islevelset=True
+md.transient.islevelset=False
 
 # init levelset function
@@ -27,9 +27,28 @@
 xmin=min(md.mesh.x)
 xmax=max(md.mesh.x)
-xmed=(xmax+xmin)/2
-ymed=(ymax+ymin)/2
-distx=numpy.absolute(md.mesh.x.reshape(-1,1)-xmed)
-disty=numpy.absolute(md.mesh.y.reshape(-1,1)-ymed)
-md.mask.ice_levelset=numpy.maximum(distx,disty)-1.e5
+xmed=(xmax+xmin)/2.
+ymed=(ymax+ymin)/2.
+md.mask.ice_levelset=numpy.sqrt(numpy.power(md.mesh.x-xmed,2.)+numpy.power(md.mesh.y-ymed,2.)) - (xmax-xmin)/3.
+
+# set spcs
+mask=1.*numpy.ones((md.mesh.numberofvertices,1))
+nrverts=md.mesh.elements.shape[1]
+for i in range(0,md.mesh.numberofelements):
+	elt=numpy.copy(md.mesh.elements[i])
+	elt-=1
+	isiceinelement=False
+	for iv in range(0,nrverts):
+		if(md.mask.ice_levelset[elt[iv]]<=0.):
+			isiceinelement=True
+	if(isiceinelement):
+		for iv in range(0,nrverts):
+			mask[elt[iv]]=2.
+
+v=0.
+for i in range(0,md.mesh.numberofvertices):
+	if(mask[i]==1.):
+		md.stressbalance.spcvx[i]=v
+		md.stressbalance.spcvy[i]=v
+		md.stressbalance.spcvz[i]=v
 
 md=solve(md,TransientSolutionEnum())
@@ -37,5 +56,5 @@
 #Fields and tolerances to track changes
 field_names     =['Vx','Vy','Vel','Pressure','MaskIceLevelset']
-field_tolerances=[1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13]
+field_tolerances=[1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13]#,1e-13,1e-13]
 field_values=[\
 	md.results.TransientSolution[0].Vx,\
@@ -43,9 +62,9 @@
 	md.results.TransientSolution[0].Vel,\
 	md.results.TransientSolution[0].Pressure,\
-	md.results.TransientSolution[0].MaskIceLevelset,\
+# 	md.results.TransientSolution[0].MaskIceLevelset,\
 	md.results.TransientSolution[1].Vx,\
 	md.results.TransientSolution[1].Vy,\
 	md.results.TransientSolution[1].Vel,\
 	md.results.TransientSolution[1].Pressure,\
-	md.results.TransientSolution[1].MaskIceLevelset,\
+# 	md.results.TransientSolution[1].MaskIceLevelset,\
 	]