Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 25226)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 25227)
@@ -206,7 +206,6 @@
 	/*Initialize Element vector*/
 	ElementMatrix* Ke     = basalelement->NewElementMatrix();
-	IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
 	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
-	IssmDouble     D[2][2]={0.};
+	IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -222,4 +221,6 @@
 		gauss           ->GaussPoint(ig);
 		basalelement    ->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+		basalelement->NodalFunctions(basis,gauss);
 
 		epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input,epl_head_input,base_input);
@@ -230,11 +231,9 @@
 		//D_scalar=gauss->weight*Jdet;
 		if(dt!=0.) D_scalar=D_scalar*dt;
-		D[0][0]=D_scalar;
-		D[1][1]=D_scalar;
-		GetB(B,basalelement,xyz_list,gauss);
-		TripleMultiply(B,2,numnodes,1,
-					&D[0][0],2,2,0,
-					B,2,numnodes,0,
-					&Ke->values[0],1);
+		for(int i=0;i<numnodes;i++){
+			for(int j=0;j<numnodes;j++){
+				Ke->values[i*numnodes+j] += D_scalar*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]);
+			}
+		}
 
 		/*Transient*/
@@ -243,8 +242,5 @@
 			D_scalar=epl_storing*gauss->weight*Jdet;
 			//D_scalar=(epl_storing/epl_transmitivity)*gauss->weight*Jdet;
-			TripleMultiply(basis,numnodes,1,0,
-						&D_scalar,1,1,0,
-						basis,1,numnodes,0,
-						&Ke->values[0],1);
+			for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
 
 			/*Transfer EPL part*/
@@ -252,8 +248,5 @@
 			D_scalar=dt*transfer*gauss->weight*Jdet;
 			//D_scalar=dt*(transfer/epl_transmitivity)*gauss->weight*Jdet;
-			TripleMultiply(basis,numnodes,1,0,
-										 &D_scalar,1,1,0,
-										 basis,1,numnodes,0,
-										 &Ke->values[0],1);
+			for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
 		}
 	}
@@ -262,5 +255,5 @@
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(basis);
-	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(dbasis);
 	delete gauss;
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
@@ -405,31 +398,4 @@
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	return pe;
-}/*}}}*/
-void HydrologyDCEfficientAnalysis::GetB(IssmDouble* B,Element* element,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=[ 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>(2*numnodes);
-	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-	/*Build B: */
-	for(int i=0;i<numnodes;i++){
-		B[numnodes*0+i] = dbasis[0*numnodes+i];
-		B[numnodes*1+i] = dbasis[1*numnodes+i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
 }/*}}}*/
 void HydrologyDCEfficientAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h	(revision 25226)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h	(revision 25227)
@@ -29,5 +29,4 @@
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
-		void GetB(IssmDouble* B,Element* element,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/HydrologyDCInefficientAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 25226)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 25227)
@@ -228,7 +228,6 @@
 	/*Initialize Element vector*/
 	ElementMatrix* Ke     = basalelement->NewElementMatrix();
-	IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
 	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
-	IssmDouble     D[2][2]= {0.};
+	IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -252,4 +251,6 @@
 		gauss          -> GaussPoint(ig);
 		basalelement   -> JacobianDeterminant(&Jdet,xyz_list,gauss);
+		basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+		basalelement->NodalFunctions(basis,gauss);
 
 		sediment_transmitivity = SedimentTransmitivity(basalelement,gauss,sed_head_input,base_input,SedTrans_input);
@@ -260,21 +261,15 @@
 		//D_scalar=gauss->weight*Jdet;
 		if(dt!=0.) D_scalar=D_scalar*dt;
-		D[0][0]=D_scalar;
-		D[1][1]=D_scalar;
-		GetB(B,basalelement,xyz_list,gauss);
-		TripleMultiply(B,2,numnodes,1,
-									 &D[0][0],2,2,0,
-									 B,2,numnodes,0,
-									 &Ke->values[0],1);
+		for(int i=0;i<numnodes;i++){
+			for(int j=0;j<numnodes;j++){
+				Ke->values[i*numnodes+j] += D_scalar*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]);
+			}
+		}
 
 		/*Transient*/
 		if(dt!=0.){
-			basalelement->NodalFunctions(&basis[0],gauss);
 			D_scalar=sediment_storing*gauss->weight*Jdet;
 			//D_scalar=(sediment_storing/sediment_transmitivity)*gauss->weight*Jdet;
-			TripleMultiply(basis,numnodes,1,0,
-										 &D_scalar,1,1,0,
-										 basis,1,numnodes,0,
-										 &Ke->values[0],1);
+			for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
 
 			/*Transfer EPL part*/
@@ -282,11 +277,7 @@
 				if(active_element){
 					transfer=GetHydrologyKMatrixTransfer(basalelement);
-					basalelement->NodalFunctions(&basis[0],gauss);
 					D_scalar=dt*transfer*gauss->weight*Jdet;
 					//D_scalar=dt*(transfer/sediment_transmitivity)*gauss->weight*Jdet;
-					TripleMultiply(basis,numnodes,1,0,
-												 &D_scalar,1,1,0,
-												 basis,1,numnodes,0,
-												 &Ke->values[0],1);
+					for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
 				}
 			}
@@ -295,6 +286,6 @@
 	/*Clean up and return*/
 	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(B);
 	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(dbasis);
 	delete gauss;
 	if(domaintype!=Domain2DhorizontalEnum){
@@ -473,31 +464,4 @@
 	}
 	return pe;
-}/*}}}*/
-void HydrologyDCInefficientAnalysis::GetB(IssmDouble* B,Element* element,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=[ 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>(2*numnodes);
-	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-	/*Build B: */
-	for(int i=0;i<numnodes;i++){
-		B[numnodes*0+i] = dbasis[0*numnodes+i];
-		B[numnodes*1+i] = dbasis[1*numnodes+i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
 }/*}}}*/
 void HydrologyDCInefficientAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h	(revision 25226)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h	(revision 25227)
@@ -33,5 +33,4 @@
 
 		/*Intermediaries*/
-		void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
 		IssmDouble SedimentStoring(Element* element, Gauss* gauss, Input2* sed_head_input, Input2* base_input);
 		IssmDouble SedimentTransmitivity(Element* element,Gauss* gauss,Input2* sed_head_input, Input2* base_input,Input2* SedTrans_input);