Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 17426)
+++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 17427)
@@ -29,5 +29,5 @@
 		}
 	}
-
+	
 	iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
 	iomodel->FetchDataToInput(elements,VxEnum);
@@ -35,9 +35,14 @@
 	iomodel->FetchDataToInput(elements,MasstransportCalvingrateEnum);
 	
+	if(iomodel->meshtype==Mesh3DEnum){
+		iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
+	}
 }
 /*}}}*/
 void LevelsetAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
 	int finiteelement=P1Enum;
+	if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
 	::CreateNodes(nodes,iomodel,LevelsetAnalysisEnum,finiteelement);
+	iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
 }
 /*}}}*/
@@ -92,6 +97,9 @@
 ElementMatrix* LevelsetAnalysis::CreateKMatrix(Element* element){/*{{{*/
 
+	if(!element->IsOnBed()) return NULL;
+	Element* basalelement = element->SpawnBasalElement();
+
 	/*Intermediaries */
-	int  dim = 2; // solve for LSF in horizontal plane only
+	int  dim, meshtype; // solve for LSF in horizontal plane only
 	int i, row, col;
 	IssmDouble kappa;
@@ -102,9 +110,18 @@
 	IssmDouble* xyz_list = NULL;
 
+	/*Get problem dimension*/
+	basalelement->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DverticalEnum:   dim = 1; break;
+		case Mesh2DhorizontalEnum: dim = 2; break;
+		case Mesh3DEnum:           dim = 2; break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
 	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes    = element->GetNumberOfNodes();
+	int numnodes    = basalelement->GetNumberOfNodes();
 
 	/*Initialize Element vector and other vectors*/
-	ElementMatrix* Ke       = element->NewElementMatrix();
+	ElementMatrix* Ke       = basalelement->NewElementMatrix();
 	IssmDouble*    basis    = xNew<IssmDouble>(numnodes);
 	IssmDouble*    B        = xNew<IssmDouble>(dim*numnodes);
@@ -117,23 +134,37 @@
 
 	/*Retrieve all inputs and parameters*/
-	element->GetVerticesCoordinates(&xyz_list);
-	element->FindParam(&dt,TimesteppingTimeStepEnum);
-	Input* vx_input  = element->GetInput(VxEnum);     _assert_(vx_input);
-	Input* vy_input  = element->GetInput(VyEnum);     _assert_(vy_input);
-	Input* lsf_slopex_input  = element->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
-	Input* lsf_slopey_input  = element->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
-	Input* calvingrate_input  = element->GetInput(MasstransportCalvingrateEnum);     _assert_(calvingrate_input);
+	basalelement->GetVerticesCoordinates(&xyz_list);
+	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
+	Input* vx_input  = NULL;
+	Input* vy_input  = NULL;
+	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: */
-	Gauss* gauss=element->NewGauss(2);
+	Gauss* gauss=basalelement->NewGauss(2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 		gauss->GaussPoint(ig);
 
-		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		D_scalar=gauss->weight*Jdet;
 
 		/* Transient */
 		if(dt!=0.){
-			element->NodalFunctions(basis,gauss);
+			basalelement->NodalFunctions(basis,gauss);
 			TripleMultiply(basis,numnodes,1,0,
 						&D_scalar,1,1,0,
@@ -144,6 +175,6 @@
 
 		/* Advection */
-		GetB(B,element,xyz_list,gauss); 
-		GetBprime(Bprime,element,xyz_list,gauss); 
+		GetB(B,basalelement,xyz_list,gauss); 
+		GetBprime(Bprime,basalelement,xyz_list,gauss); 
 		vx_input->GetInputValue(&v[0],gauss); // in 3D case, add mesh velocity 
 		vy_input->GetInputValue(&v[1],gauss); 
@@ -186,5 +217,5 @@
 			case 1:
 				/* Artificial Diffusion */
-				element->ElementSizes(&hx,&hy,&hz);
+				basalelement->ElementSizes(&hx,&hy,&hz);
 				h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) ); 
 				kappa=h*vel/2.;
@@ -203,5 +234,5 @@
 			case 2:
 				/* Streamline Upwinding */
-				element->ElementSizes(&hx,&hy,&hz);
+				basalelement->ElementSizes(&hx,&hy,&hz);
 				h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) );
 				for(row=0;row<dim;row++) 
@@ -229,12 +260,15 @@
 	xDelete<IssmDouble>(c);
 	xDelete<IssmDouble>(dlsf);
+	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	delete gauss;
 	return Ke;
 }/*}}}*/
 ElementVector* LevelsetAnalysis::CreatePVector(Element* element){/*{{{*/
+	
+	if(!element->IsOnBed()) return NULL;
+	Element* basalelement = element->SpawnBasalElement();
 
 	/*Intermediaries */
-	const int dim = 2;
-	int i, ig, k;
+	int i, ig, meshtype;
 	IssmDouble  Jdet,dt;
 	IssmDouble  lsf;
@@ -242,9 +276,9 @@
 	
 	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
+	int numnodes = basalelement->GetNumberOfNodes();
 
 	/*Initialize Element vector*/
-	ElementVector* pe = element->NewElementVector();
-	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	ElementVector* pe = basalelement->NewElementVector();
+	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
 	
 	if(dt!=0.){
@@ -253,14 +287,14 @@
 
 		/*Retrieve all inputs and parameters*/
-		element->GetVerticesCoordinates(&xyz_list);
-		Input* levelset_input     = element->GetInput(MaskIceLevelsetEnum);                    _assert_(levelset_input);
+		basalelement->GetVerticesCoordinates(&xyz_list);
+		Input* levelset_input     = basalelement->GetInput(MaskIceLevelsetEnum);                    _assert_(levelset_input);
 
 		/* Start  looping on the number of gaussian points: */
-		Gauss* gauss=element->NewGauss(2);
+		Gauss* gauss=basalelement->NewGauss(2);
 		for(ig=gauss->begin();ig<gauss->end();ig++){
 			gauss->GaussPoint(ig);
 
-			element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-			element->NodalFunctions(basis,gauss);
+			basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+			basalelement->NodalFunctions(basis,gauss);
 
 			/* old function value */
@@ -277,4 +311,5 @@
 		for(i=0;i<numnodes;i++) 
 			pe->values[i]=0.; 
+	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	return pe;
 }/*}}}*/