Index: /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 24430)
+++ /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 24431)
@@ -207,58 +207,4 @@
 ElementVector* ExtrapolationAnalysis::CreatePVector(Element* element){/*{{{*/
 	return NULL;
-
-}/*}}}*/
-void           ExtrapolationAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss, int dim){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
-	 * For node i, Bi can be expressed in the actual coordinate system
-	 * by: 
-	 *       Bi=[ N ]
-	 *          [ N ]
-	 * where N is the finiteelement function for node i.
-	 *
-	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Get nodal functions*/
-	IssmDouble* basis=xNew<IssmDouble>(numnodes);
-	element->NodalFunctions(basis,gauss);
-
-	/*Build B: */
-	for(int i=0;i<numnodes;i++)
-		for(int d=0;d<dim;d++)
-			B[numnodes*d+i] = basis[i];
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(basis);
-}/*}}}*/
-void           ExtrapolationAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss, int dim){/*{{{*/
-	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
-	 * For node i, Bi' can be expressed in the actual coordinate system
-	 * by: 
-	 *       Bi_prime=[ dN/dx ]
-	 *                [ dN/dy ]
-	 * where N is the finiteelement function for node i.
-	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Get nodal functions derivatives*/
-	IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
-	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-	/*Build B': */
-	for(int i=0;i<numnodes;i++)
-		for(int d=0;d<dim;d++)
-			Bprime[numnodes*d+i] = dbasis[d*numnodes+i];
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
-
 }/*}}}*/
 void           ExtrapolationAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
@@ -288,5 +234,5 @@
 	}
 }/*}}}*/
-int				ExtrapolationAnalysis::GetExtrapolationCase(Element* element){/*{{{*/
+int            ExtrapolationAnalysis::GetExtrapolationCase(Element* element){/*{{{*/
 
 	/* Get case of extrapolation, depending on domain quality, and extrapolation variable */
Index: /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.h	(revision 24430)
+++ /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.h	(revision 24431)
@@ -26,6 +26,4 @@
 	ElementMatrix* CreateKMatrix(Element* element);
 	ElementVector* CreatePVector(Element* element);
-	void           GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss, int dim);
-	void           GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss, int dim);
 	void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
 	void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 24430)
+++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 24431)
@@ -291,4 +291,5 @@
 	IssmDouble xi,tau;
 	IssmDouble dvx[2],dvy[2];
+	IssmDouble D[4];
 	IssmDouble* xyz_list = NULL;
 
@@ -309,7 +310,4 @@
 	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
 	IssmDouble*		dbasis = xNew<IssmDouble>(dim*numnodes);
-	IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
-	IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
-	IssmDouble*    D      = xNewZeroInit<IssmDouble>(dim*dim);
 
 	/*Retrieve all inputs and parameters*/
@@ -333,41 +331,43 @@
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		element->NodalFunctions(basis,gauss);
-
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+		/*Transient term*/
+		D_scalar=gauss->weight*Jdet;
+		for(int i=0;i<numnodes;i++){
+			for(int j=0;j<numnodes;j++){
+				Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j];
+			}
+		}
+
+		/*Advection terms*/
 		vxaverage_input->GetInputValue(&vx,gauss);
 		vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+		D_scalar=dt*gauss->weight*Jdet;
 		if(dim==2){
 			vyaverage_input->GetInputValue(&vy,gauss);
 			vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
-		}
-
-		D_scalar=gauss->weight*Jdet;
-		TripleMultiply(basis,1,numnodes,1,
-					&D_scalar,1,1,0,
-					basis,1,numnodes,0,
-					&Ke->values[0],1);
-
-		GetB(B,element,dim,xyz_list,gauss);
-		GetBprime(Bprime,element,dim,xyz_list,gauss);
-
-		dvxdx=dvx[0];
-		if(dim==2) dvydy=dvy[1];
-		D_scalar=dt*gauss->weight*Jdet;
-
-		D[0*dim+0]=D_scalar*dvxdx;
-		if(dim==2) D[1*dim+1]=D_scalar*dvydy;
-
-		TripleMultiply(B,dim,numnodes,1,
-					D,dim,dim,0,
-					B,dim,numnodes,0,
-					&Ke->values[0],1);
-
-		D[0*dim+0]=D_scalar*vx;
-		if(dim==2) D[1*dim+1]=D_scalar*vy;
-
-		TripleMultiply(B,dim,numnodes,1,
-					D,dim,dim,0,
-					Bprime,dim,numnodes,0,
-					&Ke->values[0],1);
-
+			dvxdx=dvx[0];
+			dvydy=dvy[1];
+			for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					/*\phi_i \phi_j \nabla\cdot v*/
+					Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j]*(dvxdx+dvydy);
+					/*\phi_i v\cdot\nabla\phi_j*/
+					Ke->values[i*numnodes+j] += D_scalar*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j]);
+				}
+			}
+		}
+		else{
+			dvxdx=dvx[0];
+			for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					Ke->values[i*numnodes+j] += D_scalar*dvxdx*dbasis[0*numnodes+i]*dbasis[0*numnodes+j];
+					Ke->values[i*numnodes+j] += D_scalar*vx*dbasis[0*numnodes+j]*basis[i];
+				}
+			}
+		}
+
+		for(int i=0;i<4;i++) D[i]=0.;
 		switch(stabilization){
 			case 0:
@@ -383,5 +383,4 @@
 			case 2:
 				/*Streamline upwinding*/
-				element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 				vxaverage_input->GetInputAverage(&vx);
 				if(dim==1){
@@ -397,5 +396,4 @@
 				/*SUPG*/
 				if(dim!=2) _error_("Stabilization "<<stabilization<<" not supported yet for dim != 2");
-				element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 				vxaverage_input->GetInputAverage(&vx);
 				vyaverage_input->GetInputAverage(&vy);
@@ -418,8 +416,21 @@
 			}
 
-			TripleMultiply(Bprime,dim,numnodes,1,
-						D,dim,dim,0,
-						Bprime,dim,numnodes,0,
-						&Ke->values[0],1);
+			if(dim==2){
+				for(int i=0;i<numnodes;i++){
+					for(int j=0;j<numnodes;j++){
+						Ke->values[i*numnodes+j] += (
+									dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) +
+									dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j]) 
+									);
+					}
+				}
+			}
+			else{
+				for(int i=0;i<numnodes;i++){
+					for(int j=0;j<numnodes;j++){
+						Ke->values[i*numnodes+j] += D_scalar*h/(2.*vel)*dbasis[0*numnodes+i] *D[0]*dbasis[0*numnodes+j];
+					}
+				}
+			}
 		}
 		if(stabilization==2){
@@ -478,7 +489,4 @@
 	xDelete<IssmDouble>(basis);
 	xDelete<IssmDouble>(dbasis);
-	xDelete<IssmDouble>(B);
-	xDelete<IssmDouble>(Bprime);
-	xDelete<IssmDouble>(D);
 	delete gauss;
 	return Ke;
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 24430)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 24431)
@@ -47,4 +47,14 @@
 
 /*Other*/
+bool       Element::AnyFSet(){/*{{{*/
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = this->GetNumberOfNodes();
+
+	for(int i=0;i<numnodes;i++){
+		if(nodes[i]->fsize) return true;
+	}
+	return false;
+}/*}}}*/
 void       Element::ComputeLambdaS(){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 24430)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 24431)
@@ -63,4 +63,5 @@
 		/*bool               AllActive(void);*/
 		/*bool               AnyActive(void);*/
+		bool               AnyFSet(void);
 		void               ComputeLambdaS(void);
 		void               ComputeNewDamage();
Index: /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 24430)
+++ /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 24431)
@@ -43,4 +43,5 @@
 		for (i=0;i<femmodel->elements->Size();i++){
 			element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+			if(!element->AnyFSet()) continue;
 			ElementMatrix* Ke = analysis->CreateKMatrix(element);
 			ElementVector* pe = analysis->CreatePVector(element);
@@ -76,4 +77,5 @@
 	for (i=0;i<femmodel->elements->Size();i++){
 		element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		if(!element->AnyFSet()) continue;
 		ElementMatrix* Ke = analysis->CreateKMatrix(element);
 		ElementVector* pe = analysis->CreatePVector(element);