Index: /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 17425)
+++ /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 17426)
@@ -33,5 +33,7 @@
 void ExtrapolationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
 	int finiteelement=P1Enum;
+	if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
 	::CreateNodes(nodes,iomodel,ExtrapolationAnalysisEnum,finiteelement);
+	iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
 }
 /*}}}*/
@@ -76,37 +78,50 @@
 ElementMatrix* ExtrapolationAnalysis::CreateKMatrix(Element* element){/*{{{*/
 
+	if(!element->IsOnBed()) return NULL;
+	Element* basalelement = element->SpawnBasalElement();
+
 	/*Intermediaries */
-	const int dim = 2;
+	int		  meshtype,dim;
 	int        i,row,col,stabilization;
 	bool	   extrapolatebydiffusion = true;
 	IssmDouble Jdet,D_scalar,h;
-	IssmDouble dlsf[dim],normal[dim];
 	IssmDouble norm_dlsf;
 	IssmDouble hx,hy,hz,kappa;
 	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*    B      = xNew<IssmDouble>(dim*numnodes);
 	IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
-	IssmDouble     D[dim][dim];
+	IssmDouble*		dlsf   = xNew<IssmDouble>(dim);
+	IssmDouble*		normal = xNew<IssmDouble>(dim);
+	IssmDouble*    D		 = xNew<IssmDouble>(dim*dim);
 
 	/*Retrieve all inputs and parameters*/
-	Input* lsf_slopex_input=element->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
-	Input* lsf_slopey_input=element->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
-	element->GetVerticesCoordinates(&xyz_list);
-	h = element->CharacteristicLength();
+	Input* lsf_slopex_input=basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
+	Input* lsf_slopey_input=basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
+	basalelement->GetVerticesCoordinates(&xyz_list);
+	h = basalelement->CharacteristicLength();
 
 	/* 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);
-		GetB(B,element,xyz_list,gauss);
-		GetBprime(Bprime,element,xyz_list,gauss);
+		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		GetB(B,basalelement,xyz_list,gauss);
+		GetBprime(Bprime,basalelement,xyz_list,gauss);
 		
 		D_scalar=gauss->weight*Jdet;
@@ -117,10 +132,10 @@
 				for(col=0;col<dim;col++)
 					if(row==col)
-						D[row][col] = D_scalar;
+						D[row*dim+col] = D_scalar;
 					else
-						D[row][col] = 0.;
+						D[row*dim+col] = 0.;
 
 			TripleMultiply(Bprime,dim,numnodes,1,
-					&D[0][0],dim,dim,0,
+					D,dim,dim,0,
 					Bprime,dim,numnodes,0,
 					&Ke->values[0],1);
@@ -143,9 +158,9 @@
 				for(col=0;col<dim;col++)
 					if(row==col)
-						D[row][col]=D_scalar*normal[row];
+						D[row*dim+col] = D_scalar*normal[row];
 					else
-						D[row][col]=0.;
+						D[row*dim+col] = 0.;
 			TripleMultiply(B,dim,numnodes,1,
-						&D[0][0],dim,dim,0,
+						D,dim,dim,0,
 						Bprime,dim,numnodes,0,
 						&Ke->values[0],1);
@@ -156,13 +171,15 @@
 			else if(stabilization==1){
 				/* Artificial Diffusion */
-				element->ElementSizes(&hx,&hy,&hz);
+				basalelement->ElementSizes(&hx,&hy,&hz);
 				h=sqrt( pow(hx*normal[0],2) + pow(hy*normal[1],2));
 				kappa=h/2.+1.e-14; 
-				D[0][0]=D_scalar*kappa;
-				D[0][1]=0.;
-				D[1][0]=0.;
-				D[1][1]=D_scalar*kappa;
+				for(row=0;row<dim;row++)
+					for(col=0;col<dim;col++)
+						if(row==col)
+							D[row*dim+col] = D_scalar*kappa;
+						else
+							D[row*dim+col] = 0.;
 				TripleMultiply(Bprime,dim,numnodes,1,
-							&D[0][0],dim,dim,0,
+							D,dim,dim,0,
 							Bprime,dim,numnodes,0,
 							&Ke->values[0],1);
@@ -172,8 +189,8 @@
 				for(row=0;row<dim;row++)
 					for(col=0;col<dim;col++)
-						D[row][col]=h/(2.*1.)*normal[row]*normal[col];
+						D[row*dim+col]=h/(2.*1.)*normal[row]*normal[col];
 
 				TripleMultiply(Bprime,dim,numnodes,1,
-							&D[0][0],dim,dim,0,
+							D,dim,dim,0,
 							Bprime,dim,numnodes,0,
 							&Ke->values[0],1);
@@ -181,4 +198,8 @@
 		}
 	}/*}}}*/
+
+	//TESTING
+	_printf_("in element: " << basalelement->Id() << "\n");
+	Ke->Echo();
 
 	/*Clean up and return*/
@@ -186,4 +207,6 @@
 	xDelete<IssmDouble>(B);
 	xDelete<IssmDouble>(Bprime);
+	xDelete<IssmDouble>(dlsf);
+	xDelete<IssmDouble>(normal);
 	delete gauss;
 	return Ke;
@@ -191,4 +214,7 @@
 }/*}}}*/
 ElementVector* ExtrapolationAnalysis::CreatePVector(Element* element){/*{{{*/
+
+	if(!element->IsOnBed()) return NULL;
+	Element* basalelement = element->SpawnBasalElement();
 
 	/*Intermediaries */
@@ -196,8 +222,8 @@
 	
 	/*Fetch number of nodes */
-	int numnodes = element->GetNumberOfNodes();
+	int numnodes = basalelement->GetNumberOfNodes();
 
 	/*Initialize Element vector*/
-	ElementVector* pe = element->NewElementVector();
+	ElementVector* pe = basalelement->NewElementVector();
 	for(i=0;i<numnodes;i++) 
 		pe->values[i]=0.;