Index: /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp	(revision 22976)
+++ /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp	(revision 22977)
@@ -62,6 +62,5 @@
 	int        domaintype,dim;
 	IssmDouble Jdet,D_scalar;
-	IssmDouble vel,vx,vy,dvxdx,dvydy;
-	IssmDouble dvx[2],dvy[2];
+	IssmDouble vel,vx,vy;
 	IssmDouble* xyz_list = NULL;
 
@@ -80,71 +79,55 @@
 	ElementMatrix* Ke     = element->NewElementMatrix();
 	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
-	IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
+	IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
 	IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
-	IssmDouble    D[2][2] = {0.};
+	IssmDouble     D[2][2] = {0.};
 
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
-	Input* vx_input=element->GetInput(BaseSlopeXEnum); _assert_(vx_input);
-	Input* vy_input=element->GetInput(BaseSlopeYEnum); _assert_(vy_input);
+	//Input* vx_input=element->GetInput(BaseSlopeXEnum); _assert_(vx_input);
+	//Input* vy_input=element->GetInput(BaseSlopeYEnum); _assert_(vy_input);
+	Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
+	Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
 
 	IssmDouble h = element->CharacteristicLength();
 
 	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss=element->NewGauss(2);
+	Gauss* gauss=element->NewGauss(4);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 		gauss->GaussPoint(ig);
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		GetB(B,element,dim,xyz_list,gauss);
+		element->NodalFunctions(basis,gauss);
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 		GetBprime(Bprime,element,dim,xyz_list,gauss);
 
+		D_scalar=gauss->weight*Jdet;
+
+		/*Diffusion */
+		if(sqrt(vx*vx+vy*vy)<10./31536000.){
+			IssmPDouble kappa = 1.e-4;
+			for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					Ke->values[i*numnodes+j] += D_scalar*kappa*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]);
+				}
+			}
+		}
+
+		/*Advection: */
 		vx_input->GetInputValue(&vx,gauss);
 		vy_input->GetInputValue(&vy,gauss);
-		vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
-		vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
-		dvxdx=dvx[0];
-		dvydy=dvy[1];
-		D_scalar=gauss->weight*Jdet;
-
-		/*Advection*/
-		//D[0][0]=D_scalar*dvxdx;
-		//D[0][1]=0.;
-		//D[1][0]=0.;
-		//D[1][1]=D_scalar*dvydy;
-		//TripleMultiply(B,2,numnodes,1,
-		//			&D[0][0],2,2,0,
-		//			B,2,numnodes,0,
-		//			&Ke->values[0],1);
-
-		//vx = 0.;
-		//vy = -.1;
-		if(sqrt(vx*vx+vy*vy)>0.28){
-			/*Advection*/
-			D[0][0]=D_scalar*vx;
-			D[1][1]=D_scalar*vy;
-			TripleMultiply(B,2,numnodes,1,
-						&D[0][0],2,2,0,
-						Bprime,2,numnodes,0,
-						&Ke->values[0],1);
-		}
-		else{
-			/*Diffusion*/
-			D[0][0]=1e5*h/2.0;
-			D[1][1]=1e5*h/2.0;
-			TripleMultiply(Bprime,2,numnodes,1,
-						&D[0][0],2,2,0,
-						Bprime,2,numnodes,0,
-						&Ke->values[0],1);
-		}
-
-		/*Stabilization*/
-		D[0][0]=h/2.0*fabs(vx);
-		D[0][1]=0.;
-		D[1][0]=0.;
-		D[1][1]=h/2.0*fabs(vy);
-		TripleMultiply(Bprime,2,numnodes,1,
+		for(int i=0;i<numnodes;i++){
+			for(int j=0;j<numnodes;j++){
+				Ke->values[i*numnodes+j] += D_scalar*(vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i]);
+			}
+		}
+
+		/*Artifficial diffusivity*/
+		vel=sqrt(vx*vx + vy*vy)+1.e-14;
+		D[0][0]=D_scalar*h/(2.*vel)*fabs(vx*vx);  D[0][1]=D_scalar*h/(2.*vel)*fabs(vx*vy);
+		D[1][0]=D_scalar*h/(2.*vel)*fabs(vy*vx);  D[1][1]=D_scalar*h/(2.*vel)*fabs(vy*vy);
+		TripleMultiply(Bprime,dim,numnodes,1,
 					&D[0][0],2,2,0,
-					Bprime,2,numnodes,0,
+					Bprime,dim,numnodes,0,
 					&Ke->values[0],1);
 	}
@@ -152,5 +135,6 @@
 	/*Clean up and return*/
 	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(dbasis);
 	xDelete<IssmDouble>(Bprime);
 	delete gauss;