Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 22509)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 22510)
@@ -2336,7 +2336,6 @@
 
 	/*Initialize Element matrix and vectors*/
-	ElementMatrix* Ke = element->NewElementMatrix(HOApproximationEnum);
-	IssmDouble*    B  = xNew<IssmDouble>((dim-1)*numdof);
-	IssmDouble*    D  = xNewZeroInit<IssmDouble>((dim-1)*(dim-1));
+	ElementMatrix* Ke     = element->NewElementMatrix(HOApproximationEnum);
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -2370,12 +2369,22 @@
 		}
 
-		this->GetBHOFriction(B,element,dim,xyz_list_base,gauss);
 		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
-		for(int i=0;i<dim-1;i++) D[i*(dim-1)+i]=alpha2*gauss->weight*Jdet;
-
-		TripleMultiply(B,dim-1,numdof,1,
-					D,dim-1,dim-1,0,
-					B,dim-1,numdof,0,
-					&Ke->values[0],1);
+		element->NodalFunctions(basis,gauss);
+
+		if(dim==3){
+			for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					Ke->values[2*i*2*numnodes+2*j]       += alpha2*gauss->weight*Jdet*basis[i]*basis[j];
+					Ke->values[(2*i+1)*2*numnodes+2*j+1] += alpha2*gauss->weight*Jdet*basis[i]*basis[j];
+				}
+			}
+		}
+		else{
+			for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					Ke->values[i*numnodes+j] += alpha2*gauss->weight*Jdet*basis[i]*basis[j];
+				}
+			}
+		}
 	}
 
@@ -2387,6 +2396,5 @@
 	delete friction;
 	xDelete<IssmDouble>(xyz_list_base);
-	xDelete<IssmDouble>(B);
-	xDelete<IssmDouble>(D);
+	xDelete<IssmDouble>(basis);
 	return Ke;
 }/*}}}*/
@@ -2400,11 +2408,8 @@
 	IssmDouble  viscosity,newviscosity,oldviscosity;
 	IssmDouble  viscosity_overshoot,thickness,Jdet;
-	IssmDouble  D_scalar;
 	IssmDouble *xyz_list = NULL;
 
 	/*Get problem dimension*/
 	element->FindParam(&dim,DomainDimensionEnum);
-	if(dim==2) bsize = 2;
-	else       bsize = 5;
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -2414,7 +2419,5 @@
 	/*Initialize Element matrix and vectors*/
 	ElementMatrix* Ke     = element->NewElementMatrix(HOApproximationEnum);
-	IssmDouble*    B      = xNew<IssmDouble>(bsize*numdof);
-	IssmDouble*    Bprime = xNew<IssmDouble>(bsize*numdof);
-	IssmDouble*    D      = xNewZeroInit<IssmDouble>(bsize*bsize);
+	IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -2436,18 +2439,37 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		this->GetBHO(B,element,dim,xyz_list,gauss);
-		this->GetBHOprime(Bprime,element,dim,xyz_list,gauss);
-
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 		element->material->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
 		element->material->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);
 
 		newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
-		D_scalar=2.*newviscosity*gauss->weight*Jdet;
-		for(int i=0;i<bsize;i++) D[i*bsize+i]=D_scalar;
-
-		TripleMultiply(B,bsize,numdof,1,
-					D,bsize,bsize,0,
-					Bprime,bsize,numdof,0,
-					&Ke->values[0],1);
+
+		if(dim==3){
+			for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					Ke->values[2*i*2*numnodes+2*j] += gauss->weight*Jdet*newviscosity*(
+								4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
+								);
+					Ke->values[2*i*2*numnodes+2*j+1] += gauss->weight*Jdet*newviscosity*(
+								2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
+								);
+					Ke->values[(2*i+1)*2*numnodes+2*j] += gauss->weight*Jdet*newviscosity*(
+								2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
+								);
+					Ke->values[(2*i+1)*2*numnodes+2*j+1] += gauss->weight*Jdet*newviscosity*(
+								dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + 4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
+								);
+				}
+			}
+		}
+		else{
+			for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					Ke->values[i*numnodes+j] += gauss->weight*Jdet*newviscosity*(
+								4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
+								);
+				}
+			}
+		}
 	}
 
@@ -2458,7 +2480,5 @@
 	delete gauss;
 	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(D);
-	xDelete<IssmDouble>(Bprime);
-	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(dbasis);
 	return Ke;
 }/*}}}*/