Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17537)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17538)
@@ -3424,4 +3424,5 @@
 	IssmDouble*    vdbasis = xNew<IssmDouble>(dim*vnumnodes);
 	IssmDouble*    tbasis  = xNew<IssmDouble>(tnumnodes);
+	IssmDouble*    D       = xNewZeroInit<IssmDouble>(tausize*tnumnodes*tausize*tnumnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -3530,4 +3531,18 @@
 			}
 		}
+
+		/*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);
 	}
 
@@ -3541,4 +3556,5 @@
 	xDelete<IssmDouble>(Dstar);
 	xDelete<IssmDouble>(d);
+	xDelete<IssmDouble>(D);
 	xDelete<IssmDouble>(tau);
 	return pe;