Index: /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 21400)
+++ /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 21401)
@@ -84,12 +84,12 @@
 
 	/*Intermediaries */
-	int	dim, domaintype, extrapolationcase;
-	int	i,row,col,stabilization;
-	bool	extrapolatebydiffusion;
-	IssmDouble Jdet,D_scalar,h;
-	IssmDouble norm_dlsf;
-	IssmDouble hx,hy,hz,kappa;
-	IssmDouble*	xyz_list = NULL;
-	Element*		workelement=NULL;
+	int	      dim, domaintype, extrapolationcase;
+	int	      i,row,col,stabilization;
+	bool	      extrapolatebydiffusion;
+	IssmDouble  Jdet,D_scalar,h;
+	IssmDouble  norm_dlsf;
+	IssmDouble  hx,hy,hz,kappa;
+	IssmDouble *xyz_list    = NULL;
+	Element    *workelement = NULL;
 
 	/*Get problem case*/
@@ -102,10 +102,12 @@
 		case 1: case 2: case 3: workelement=element; break;
 	}
+
 	/* get extrapolation dimension */
 	workelement->FindParam(&domaintype,DomainTypeEnum);
 	switch(domaintype){
-		case Domain2DverticalEnum: dim=1; break;
+		case Domain2DverticalEnum:   dim=1; break;
 		case Domain2DhorizontalEnum: dim=2; break;
-		case Domain3DEnum: dim=3; break;
+		case Domain3DEnum:           dim=3; break;
+      default: _error_("not supported yet");
 	}
 
@@ -117,7 +119,7 @@
 	IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
 	IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
-	IssmDouble*		D	  = xNew<IssmDouble>(dim*dim);
-	IssmDouble*		dlsf  = xNew<IssmDouble>(dim);
-	IssmDouble*		normal= xNew<IssmDouble>(dim);
+	IssmDouble*		dlsf   = xNew<IssmDouble>(dim);
+	IssmDouble*		normal = xNew<IssmDouble>(dim);
+   IssmDouble*		D	    = xNewZeroInit<IssmDouble>(dim*dim);
 
 	/*Retrieve all inputs and parameters*/
@@ -128,5 +130,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=workelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){/*{{{*/
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
 		gauss->GaussPoint(ig);
 
@@ -139,11 +141,7 @@
 		extrapolatebydiffusion=true;
 		if(extrapolatebydiffusion){
-			/* diffuse values outward */
-			for(row=0;row<dim;row++)
-				for(col=0;col<dim;col++)
-					if(row==col && row<2) //extrapolate only in xy-plane
-						D[row*dim+col] = D_scalar;
-					else
-						D[row*dim+col] = 0.;
+
+			/* diffuse values outward only along the xy-plane*/
+         for(int i=0;i<2;i++) D[i*dim+i] = D_scalar;
 
 			TripleMultiply(Bprime,dim,numnodes,1,
@@ -180,5 +178,5 @@
 						&Ke->values[0],1);
 
-			/* stabilization *//*{{{*/
+			/* stabilization */
 			/* do not use streamline upwinding for extrapolation: it yields oscillating results due to diffusion along normal direction, but none across */
 			stabilization=1;
@@ -200,7 +198,6 @@
 							&Ke->values[0],1);
 			}
-			/*}}}*/
 		}
-	}/*}}}*/
+	}
 
 	/*Clean up and return*/
@@ -236,8 +233,10 @@
 	/*Initialize Element vector*/
 	ElementVector* pe = workelement->NewElementVector();
-	for(int i=0;i<numnodes;i++) 
-		pe->values[i]=0.; 
-
-	if(extrapolationcase==0){workelement->DeleteMaterials(); delete workelement;};
+
+	if(extrapolationcase==0){
+      workelement->DeleteMaterials();
+      delete workelement;
+   }
+
 	return pe;
 }/*}}}*/
@@ -326,12 +325,17 @@
 	int domaintype, extrapolationvariable;
 	int extrapolationcase;
+
 	element->FindParam(&domaintype,DomainTypeEnum);
 	switch(domaintype){
-		case Domain2DverticalEnum: extrapolationcase=0; break;
-		case Domain2DhorizontalEnum: extrapolationcase=1;break;
-		case Domain3DEnum:  
+		case Domain2DverticalEnum:   extrapolationcase=0; break;
+		case Domain2DhorizontalEnum: extrapolationcase=1; break;
+		case Domain3DEnum: 
 			element->FindParam(&extrapolationvariable, ExtrapolationVariableEnum);
-			if(extrapolationvariable==ThicknessEnum) extrapolationcase=2; // scalar fields that are constant along z-axis
-			else extrapolationcase=3; // scalar fields that vary along z-axis
+			if(extrapolationvariable==ThicknessEnum){
+            extrapolationcase=2; // scalar fields that are constant along z-axis
+         }
+			else{
+            extrapolationcase=3; // scalar fields that vary along z-axis
+         }
 			break;
 	}