Index: /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 17406)
+++ /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 17407)
@@ -57,6 +57,6 @@
 	if(VerboseSolution()) _printf0_("extrapolation of " << EnumToStringx(extvar_enum) << ": call computational core:\n");
 	solutionsequence_linear(femmodel);
-	
-	save_results=true;
+
+	save_results=false;
 	if(save_results){
 		if(VerboseSolution()) _printf0_("   saving results\n");
@@ -79,4 +79,5 @@
 	const int dim = 2;
 	int        i,row,col,stabilization;
+	bool	   extrapolatebydiffusion = true;
 	IssmDouble Jdet,D_scalar,h;
 	IssmDouble dlsf[dim],normal[dim];
@@ -90,9 +91,13 @@
 	/*Initialize Element vector and other vectors*/
 	ElementMatrix* Ke     = element->NewElementMatrix();
+	IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
 	IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
 	IssmDouble     D[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();
 
 	/* Start  looping on the number of gaussian points: */
@@ -102,20 +107,82 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		GetB(B,element,xyz_list,gauss);
 		GetBprime(Bprime,element,xyz_list,gauss);
-
+		
 		D_scalar=gauss->weight*Jdet;
 
-		for(row=0;row<dim;row++)
-			for(col=0;col<dim;col++)
-				D[row][col]=((row==col)?D_scalar:0.);
-
-		TripleMultiply(Bprime,dim,numnodes,1,
+		if(extrapolatebydiffusion){
+			/* diffuse values outward */
+			for(row=0;row<dim;row++)
+				for(col=0;col<dim;col++)
+					if(row==col)
+						D[row][col] = D_scalar;
+					else
+						D[row][col] = 0.;
+
+			TripleMultiply(Bprime,dim,numnodes,1,
 					&D[0][0],dim,dim,0,
 					Bprime,dim,numnodes,0,
 					&Ke->values[0],1);
+		}
+		else{
+			/* extrapolate values along normal */
+			/* Get normal on ice boundary */
+			lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
+			lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
+			norm_dlsf=0.;
+			for(i=0;i<dim;i++)	norm_dlsf+=dlsf[i]*dlsf[i]; 
+			norm_dlsf=sqrt(norm_dlsf); _assert_(norm_dlsf>0.);
+
+			if(norm_dlsf>0.)
+				for(i=0;i<dim;i++)	normal[i]=dlsf[i]/norm_dlsf;
+			else
+				for(i=0;i<dim;i++)	normal[i]=0.;
+			
+			for(row=0;row<dim;row++)
+				for(col=0;col<dim;col++)
+					if(row==col)
+						D[row][col]=D_scalar*normal[row];
+					else
+						D[row][col]=0.;
+			TripleMultiply(B,dim,numnodes,1,
+						&D[0][0],dim,dim,0,
+						Bprime,dim,numnodes,0,
+						&Ke->values[0],1);
+
+			/* Stabilization *//*{{{*/
+			stabilization=1;
+			if (stabilization==0){/* no stabilization, do nothing*/}
+			else if(stabilization==1){
+				/* Artificial Diffusion */
+				element->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;
+				TripleMultiply(Bprime,dim,numnodes,1,
+							&D[0][0],dim,dim,0,
+							Bprime,dim,numnodes,0,
+							&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*/
 	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(B);
 	xDelete<IssmDouble>(Bprime);
 	delete gauss;