Index: /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 17362)
+++ /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 17363)
@@ -57,5 +57,11 @@
 	if(VerboseSolution()) _printf0_("extrapolation of " << EnumToStringx(extvar_enum) << ": call computational core:\n");
 	solutionsequence_linear(femmodel);
-
+	
+	save_results=true;
+	if(save_results){
+		if(VerboseSolution()) _printf0_("   saving results\n");
+		int outputs[2] = {VxEnum,VyEnum};
+		femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
+	}
 
 }/*}}}*/
@@ -84,13 +90,9 @@
 	/*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: */
@@ -100,62 +102,20 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		GetB(B,element,xyz_list,gauss);
 		GetBprime(Bprime,element,xyz_list,gauss);
 
-		/* 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.);
-		for(i=0;i<dim;i++)
-			normal[i]=dlsf[i]/norm_dlsf;
-		
 		D_scalar=gauss->weight*Jdet;
 
 		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[row][col]=((row==col)?D_scalar:0.);
+
+		TripleMultiply(Bprime,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.; 
-			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;