Index: /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 17170)
+++ /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 17171)
@@ -15,5 +15,5 @@
 /*}}}*/
 void ExtrapolationAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
-	int    stabilization,finiteelement;
+	int    finiteelement;
 
 	/*Finite element type*/
@@ -83,4 +83,5 @@
 	IssmDouble dlevelset[dim],normal[dim];
 	IssmDouble norm_dlevelset;
+	IssmDouble hx,hy,hz,kappa;
 	IssmDouble* xyz_list = NULL;
 
@@ -101,5 +102,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
+	for(int ig=gauss->begin();ig<gauss->end();ig++){/*{{{*/
 		gauss->GaussPoint(ig);
 
@@ -128,12 +129,16 @@
 					&Ke->values[0],1);
 
+		/* Stabilization *//*{{{*/
 		stabilization=1;
 		if (stabilization==0){/* no stabilization, do nothing*/}
 		else if(stabilization==1){
-			/*Streamline upwinding*/
-			for(row=0;row<dim;row++)
-				for(col=0;col<dim;col++)
-					D[row][col]=h/(2.*1.)*normal[row]*normal[col];
-
+			/* Artificial Diffusion */
+			element->ElementSizes(&hx,&hy,&hz);
+			h=sqrt( pow(hx*normal[0],2) + pow(hy*normal[1],2));
+			kappa=h/2.; 
+			D[0][0]=D_scalar*kappa;
+			D[0][1]=0.;
+			D[1][0]=0.;
+			D[1][1]=D_scalar*kappa;
 			TripleMultiply(Bprime,dim,numnodes,1,
 						&D[0][0],dim,dim,0,
@@ -141,5 +146,16 @@
 						&Ke->values[0],1);
 		}
-	}
+		else if(stabilization==2){
+			/*Streamline upwinding - do not use this for extrapolation: yields oscillating results due to smoothing along normal, not across */
+			for(row=0;row<dim;row++)
+				for(col=0;col<dim;col++)
+					D[row][col]=h/(2.*1.)*normal[row]*normal[col];
+
+			TripleMultiply(Bprime,dim,numnodes,1,
+						&D[0][0],dim,dim,0,
+						Bprime,dim,numnodes,0,
+						&Ke->values[0],1);
+		}/*}}}*/
+	}/*}}}*/
 
 	/*Clean up and return*/
@@ -150,5 +166,5 @@
 	return Ke;
 
-	}/*}}}*/
+}/*}}}*/
 ElementVector* ExtrapolationAnalysis::CreatePVector(Element* element){/*{{{*/