Index: /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 21418)
+++ /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 21419)
@@ -84,7 +84,7 @@
 
 	/*Intermediaries */
+	bool	      extrapolatebydiffusion = true;
 	int	      dim, domaintype, extrapolationcase;
 	int	      i,row,col,stabilization;
-	bool	      extrapolatebydiffusion;
 	IssmDouble  Jdet,D_scalar,h;
 	IssmDouble  norm_dlsf;
@@ -128,4 +128,9 @@
 	workelement->GetVerticesCoordinates(&xyz_list);
 
+	if(element->ObjectEnum()==PentaEnum){
+		for(i=0;i<3;i++) xyz_list[3*i+2] = 0.;
+		for(i=3;i<6;i++) xyz_list[3*i+2] = 1.;
+	}
+
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=workelement->NewGauss(2);
@@ -139,5 +144,4 @@
 		D_scalar=gauss->weight*Jdet;
 
-		extrapolatebydiffusion=true;
 		if(extrapolatebydiffusion){
 
@@ -214,30 +218,6 @@
 }/*}}}*/
 ElementVector* ExtrapolationAnalysis::CreatePVector(Element* element){/*{{{*/
-
-	/*Intermediaries */
-	Element* workelement=NULL;
-
-	/*Get problem dimension*/
-	int extrapolationcase=GetExtrapolationCase(element);
-	switch(extrapolationcase){
-		case 0: 
-			if(!element->IsOnBase()) return NULL; 
-			workelement = element->SpawnBasalElement(); 
-			break;
-		case 1: case 2: case 3: workelement=element; break;
-	}
-
-	/*Fetch number of nodes */
-	int numnodes = workelement->GetNumberOfNodes();
-
-	/*Initialize Element vector*/
-	ElementVector* pe = workelement->NewElementVector();
-
-	if(extrapolationcase==0){
-      workelement->DeleteMaterials();
-      delete workelement;
-   }
-
-	return pe;
+	return NULL;
+
 }/*}}}*/
 void           ExtrapolationAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss, int dim){/*{{{*/