Index: /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 17428)
+++ /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 17429)
@@ -29,7 +29,7 @@
 		}
 	}
-	
 	if(iomodel->meshtype==Mesh3DEnum){
 		iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
+		iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
 	}
 }
@@ -86,5 +86,5 @@
 
 	/*Intermediaries */
-	int		  meshtype,dim;
+	int		   meshtype,dim;
 	int        i,row,col,stabilization;
 	bool	   extrapolatebydiffusion = true;
@@ -110,7 +110,7 @@
 	IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
 	IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
-	IssmDouble*		dlsf   = xNew<IssmDouble>(dim);
-	IssmDouble*		normal = xNew<IssmDouble>(dim);
-	IssmDouble*    D		 = xNew<IssmDouble>(dim*dim);
+	IssmDouble*		D	  = xNew<IssmDouble>(dim*dim);
+	IssmDouble*		dlsf  = xNew<IssmDouble>(dim);
+	IssmDouble*		normal= xNew<IssmDouble>(dim);
 
 	/*Retrieve all inputs and parameters*/
@@ -162,7 +162,7 @@
 				for(col=0;col<dim;col++)
 					if(row==col)
-						D[row*dim+col] = D_scalar*normal[row];
+						D[row*dim+col]=D_scalar*normal[row];
 					else
-						D[row*dim+col] = 0.;
+						D[row*dim+col]=0.;
 			TripleMultiply(B,dim,numnodes,1,
 						D,dim,dim,0,
@@ -181,7 +181,7 @@
 					for(col=0;col<dim;col++)
 						if(row==col)
-							D[row*dim+col] = D_scalar*kappa;
+							D[row*dim+col]=D_scalar*kappa;
 						else
-							D[row*dim+col] = 0.;
+							D[row*dim+col]=0.;
 				TripleMultiply(Bprime,dim,numnodes,1,
 							D,dim,dim,0,
@@ -207,8 +207,9 @@
 	xDelete<IssmDouble>(B);
 	xDelete<IssmDouble>(Bprime);
+	xDelete<IssmDouble>(D);
 	xDelete<IssmDouble>(dlsf);
 	xDelete<IssmDouble>(normal);
+	delete gauss;
 	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
-	delete gauss;
 	return Ke;
 
@@ -224,5 +225,4 @@
 	/*Fetch number of nodes */
 	int numnodes = basalelement->GetNumberOfNodes();
-	basalelement->FindParam(&meshtype,MeshTypeEnum);
 
 	/*Initialize Element vector*/
@@ -230,4 +230,6 @@
 	for(i=0;i<numnodes;i++) 
 		pe->values[i]=0.; 
+
+	basalelement->FindParam(&meshtype,MeshTypeEnum);
 	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	return pe;
@@ -326,12 +328,14 @@
 		node=element->GetNode(in);
 		levelset_input->GetInputValue(&phi,gauss);
-		if(phi<=0.){
-			/* if ice, set dirichlet BC */
-			extvar_input->GetInputValue(&value,gauss);
-			node->ApplyConstraint(1,value);
-		}
-		else {
-			/* no ice, set no spc */
-			node->DofInFSet(0); 
+		if(node->IsActive()){
+			if(phi<=0.){
+				/* if ice, set dirichlet BC */
+				extvar_input->GetInputValue(&value,gauss);
+				node->ApplyConstraint(1,value);
+			}
+			else {
+				/* no ice, set no spc */
+				node->DofInFSet(0); 
+			}
 		}
 	}
Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 17428)
+++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 17429)
@@ -37,4 +37,5 @@
 	if(iomodel->meshtype==Mesh3DEnum){
 		iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
+		iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
 	}
 }
@@ -101,5 +102,5 @@
 
 	/*Intermediaries */
-	int  dim, meshtype; // solve for LSF in horizontal plane only
+	int  dim, meshtype;
 	int i, row, col;
 	IssmDouble kappa;
@@ -136,23 +137,23 @@
 	basalelement->GetVerticesCoordinates(&xyz_list);
 	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
-	Input* vx_input  = NULL;
-	Input* vy_input  = NULL;
+	Input* vx_input=NULL;
+	Input* vy_input=NULL;
+	if(meshtype==Mesh2DhorizontalEnum){
+		vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
+		vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
+	}
+	else{
+		if(dim==1){
+			vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
+		}
+		if(dim==2){
+			vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input);
+			vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input);
+		}
+	}
+
 	Input* lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
 	Input* lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
 	Input* calvingrate_input  = basalelement->GetInput(MasstransportCalvingrateEnum);     _assert_(calvingrate_input);
-
-	if(meshtype==Mesh2DhorizontalEnum){
-		vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
-		vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
-	}
-	else{
-		if(dim==1){
-			vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
-		}
-		if(dim==2){
-			vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input);
-			vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input);
-		}
-	}
 	
 	/* Start  looping on the number of gaussian points: */
@@ -260,6 +261,6 @@
 	xDelete<IssmDouble>(c);
 	xDelete<IssmDouble>(dlsf);
+	delete gauss;
 	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
-	delete gauss;
 	return Ke;
 }/*}}}*/
@@ -306,4 +307,6 @@
 		xDelete<IssmDouble>(xyz_list);
 		xDelete<IssmDouble>(basis);
+		basalelement->FindParam(&meshtype,MeshTypeEnum);
+		if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 		delete gauss;
 	}
@@ -311,5 +314,4 @@
 		for(i=0;i<numnodes;i++) 
 			pe->values[i]=0.; 
-	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	return pe;
 }/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Node.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Node.cpp	(revision 17428)
+++ /issm/trunk-jpl/src/c/classes/Node.cpp	(revision 17429)
@@ -97,5 +97,7 @@
 				analysis_enum==HydrologyDCInefficientAnalysisEnum ||
 				analysis_enum==DamageEvolutionAnalysisEnum || 
-				analysis_enum==HydrologyDCEfficientAnalysisEnum
+				analysis_enum==HydrologyDCEfficientAnalysisEnum ||
+				analysis_enum==LevelsetAnalysisEnum ||
+				analysis_enum==ExtrapolationAnalysisEnum
 				){
 		if(iomodel->meshtype==Mesh3DEnum || iomodel->meshtype==Mesh2DverticalEnum){