Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp	(revision 25263)
+++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp	(revision 25264)
@@ -169,6 +169,5 @@
 	ElementMatrix* Ke     = basalelement->NewElementMatrix(NoneApproximationEnum);
 	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
-	IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
-	IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
+	IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
 	IssmDouble*    D      = xNew<IssmDouble>(dim*dim);
 
@@ -189,26 +188,30 @@
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		basalelement->NodalFunctions(basis,gauss);
+		basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
 		vx_input->GetInputValue(&vx,gauss);
 		if(dim==2) vy_input->GetInputValue(&vy,gauss);
 
+		/*Transient term*/
 		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,basalelement,dim,xyz_list,gauss);
-		GetBprime(Bprime,basalelement,dim,xyz_list,gauss);
-
+		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*/
 		D_scalar=dt*gauss->weight*Jdet;
 		for(int i=0;i<dim*dim;i++) D[i]=0.;
-		D[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);
+
+		if(dim==1){
+			/*\phi_i v\cdot\nabla\phi_j*/
+			for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[i]*vx*dbasis[0*numnodes+j];
+		}
+		else{
+			_assert_(dim==2);
+			for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					/*\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]);
+				}
+			}
+		}
 
 		if(stabilization==2){
@@ -241,8 +244,17 @@
 		if(stabilization==1 || stabilization==2){
 			for(int i=0;i<dim*dim;i++) D[i]=D_scalar*D[i];
-			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] += dbasis[0*numnodes+i]*D[0]*dbasis[0*numnodes+j];
+			}
 		}
 	}
@@ -251,6 +263,5 @@
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(basis);
-	xDelete<IssmDouble>(B);
-	xDelete<IssmDouble>(Bprime);
+	xDelete<IssmDouble>(dbasis);
 	xDelete<IssmDouble>(D);
 	delete gauss;
@@ -332,61 +343,4 @@
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	return pe;
-
-}/*}}}*/
-void           FreeSurfaceBaseAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2. 
-	 * 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(1*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 j=0;j<dim;j++){
-			B[numnodes*j+i] = basis[i];
-		}
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(basis);
-}/*}}}*/
-void           FreeSurfaceBaseAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2. 
-	 * 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(2*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 j=0;j<dim;j++){
-			Bprime[numnodes*j+i] = dbasis[j*numnodes+i];
-		}
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
 
 }/*}}}*/
Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.h	(revision 25263)
+++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.h	(revision 25264)
@@ -26,6 +26,4 @@
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
-		void           GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void           GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
 		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/FreeSurfaceTopAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp	(revision 25263)
+++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp	(revision 25264)
@@ -145,6 +145,5 @@
 	ElementMatrix* Ke     = topelement->NewElementMatrix(NoneApproximationEnum);
 	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
-	IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
-	IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
+	IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
 	IssmDouble*    D      = xNew<IssmDouble>(dim*dim);
 
@@ -165,17 +164,14 @@
 		topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		topelement->NodalFunctions(basis,gauss);
+		topelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
 		vx_input->GetInputValue(&vx,gauss);
 		if(dim==2) vy_input->GetInputValue(&vy,gauss);
 
+		/*Transient term*/
 		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,topelement,dim,xyz_list,gauss);
-		GetBprime(Bprime,topelement,dim,xyz_list,gauss);
-
+		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*/
 		D_scalar=dt*gauss->weight*Jdet;
 		for(int i=0;i<dim*dim;i++) D[i]=0.;
@@ -183,8 +179,17 @@
 		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);
+		if(dim==1){
+			/*\phi_i v\cdot\nabla\phi_j*/
+			for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[i]*vx*dbasis[0*numnodes+j];
+		}
+		else{
+			_assert_(dim==2);
+			for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					/*\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]);
+				}
+			}
+		}
 
 		if(stabilization==2){
@@ -218,8 +223,17 @@
 		if(stabilization==1 || stabilization==2){
 			for(int i=0;i<dim*dim;i++) D[i]=D_scalar*D[i];
-			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] += dbasis[0*numnodes+i]*D[0]*dbasis[0*numnodes+j];
+			}
 		}
 	}
@@ -228,6 +242,5 @@
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(basis);
-	xDelete<IssmDouble>(B);
-	xDelete<IssmDouble>(Bprime);
+	xDelete<IssmDouble>(dbasis);
 	xDelete<IssmDouble>(D);
 	delete gauss;
@@ -306,61 +319,4 @@
 
 }/*}}}*/
-void           FreeSurfaceTopAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2. 
-	 * 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(1*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 j=0;j<dim;j++){
-			B[numnodes*j+i] = basis[i];
-		}
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(basis);
-}/*}}}*/
-void           FreeSurfaceTopAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2. 
-	 * 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(2*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 j=0;j<dim;j++){
-			Bprime[numnodes*j+i] = dbasis[j*numnodes+i];
-		}
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
-
-}/*}}}*/
 void           FreeSurfaceTopAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 	   _error_("not implemented yet");
Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.h	(revision 25263)
+++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.h	(revision 25264)
@@ -26,6 +26,4 @@
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
-		void           GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void           GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
 		void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
 		void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);