Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17627)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17628)
@@ -3503,9 +3503,9 @@
 			for(int i=0;i<vnumnodes;i++){
 				for(int j=0;j<tnumnodes;j++){
-					Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+0] = tbasis[j]*vdbasis[0*vnumnodes+i];
-					Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+2] = tbasis[j]*vdbasis[1*vnumnodes+i];
-
-					Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+1] = tbasis[j]*vdbasis[1*vnumnodes+i];
-					Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+2] = tbasis[j]*vdbasis[0*vnumnodes+i];
+					Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+0] += gauss->weight*Jdet*tbasis[j]*vdbasis[0*vnumnodes+i];
+					Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+2] += gauss->weight*Jdet*tbasis[j]*vdbasis[1*vnumnodes+i];
+
+					Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+1] += gauss->weight*Jdet*tbasis[j]*vdbasis[1*vnumnodes+i];
+					Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+2] += gauss->weight*Jdet*tbasis[j]*vdbasis[0*vnumnodes+i];
 				}
 			}
@@ -3514,33 +3514,33 @@
 			for(int i=0;i<vnumnodes;i++){
 				for(int j=0;j<tnumnodes;j++){
-					Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+0] = tbasis[j]*vdbasis[0*vnumnodes+i];
-					Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+3] = tbasis[j]*vdbasis[1*vnumnodes+i];
-					Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+4] = tbasis[j]*vdbasis[2*vnumnodes+i];
-
-					Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+1] = tbasis[j]*vdbasis[1*vnumnodes+i];
-					Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+3] = tbasis[j]*vdbasis[0*vnumnodes+i];
-					Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+5] = tbasis[j]*vdbasis[2*vnumnodes+i];
-
-					Dstar[(i*dim+2)*tausize*tnumnodes + j*tausize+2] = tbasis[j]*vdbasis[2*vnumnodes+i];
-					Dstar[(i*dim+2)*tausize*tnumnodes + j*tausize+4] = tbasis[j]*vdbasis[0*vnumnodes+i];
-					Dstar[(i*dim+2)*tausize*tnumnodes + j*tausize+5] = tbasis[j]*vdbasis[1*vnumnodes+i];
+					Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+0] += gauss->weight*Jdet*tbasis[j]*vdbasis[0*vnumnodes+i];
+					Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+3] += gauss->weight*Jdet*tbasis[j]*vdbasis[1*vnumnodes+i];
+					Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+4] += gauss->weight*Jdet*tbasis[j]*vdbasis[2*vnumnodes+i];
+
+					Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+1] += gauss->weight*Jdet*tbasis[j]*vdbasis[1*vnumnodes+i];
+					Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+3] += gauss->weight*Jdet*tbasis[j]*vdbasis[0*vnumnodes+i];
+					Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+5] += gauss->weight*Jdet*tbasis[j]*vdbasis[2*vnumnodes+i];
+
+					Dstar[(i*dim+2)*tausize*tnumnodes + j*tausize+2] += gauss->weight*Jdet*tbasis[j]*vdbasis[2*vnumnodes+i];
+					Dstar[(i*dim+2)*tausize*tnumnodes + j*tausize+4] += gauss->weight*Jdet*tbasis[j]*vdbasis[0*vnumnodes+i];
+					Dstar[(i*dim+2)*tausize*tnumnodes + j*tausize+5] += gauss->weight*Jdet*tbasis[j]*vdbasis[1*vnumnodes+i];
 				}
 			}
 		}
-
-		/*contribution -Dstar tau*/
-		for(i=0;i<tausize*tnumnodes;i++) D[i*(tausize*tnumnodes)+i] = -gauss->weight*Jdet;
-		TripleMultiply(Dstar,dim*tnumnodes,tausize*tnumnodes,0,
-					D,tausize*tnumnodes,tausize*tnumnodes,0,
-					tau,tausize*tnumnodes,1,0,
-					&pe->values[0],1);
-
-		/*contribution + r Dstar d*/
-		for(i=0;i<tausize*tnumnodes;i++) D[i*(tausize*tnumnodes)+i] = +r*gauss->weight*Jdet;
-		TripleMultiply(Dstar,dim*tnumnodes,tausize*tnumnodes,0,
-					D,tausize*tnumnodes,tausize*tnumnodes,0,
-					d,tausize*tnumnodes,1,0,
-					&pe->values[0],1);
-	}
+	}
+
+	/*contribution -Dstar tau*/
+	for(i=0;i<tausize*tnumnodes;i++) D[i*(tausize*tnumnodes)+i] = -1.;
+	TripleMultiply(Dstar,dim*vnumnodes,tausize*tnumnodes,0,
+				D,tausize*tnumnodes,tausize*tnumnodes,0,
+				tau,tausize*tnumnodes,1,0,
+				&pe->values[0],1);
+
+	/*contribution + r Dstar d*/
+	for(i=0;i<tausize*tnumnodes;i++) D[i*(tausize*tnumnodes)+i] = +r;
+	TripleMultiply(Dstar,dim*vnumnodes,tausize*tnumnodes,0,
+				D,tausize*tnumnodes,tausize*tnumnodes,0,
+				d,tausize*tnumnodes,1,0,
+				&pe->values[0],1);
 
 	/*Transform coordinate system*/