Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17536)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17537)
@@ -3419,5 +3419,5 @@
 	/*Initialize vectors*/
 	ElementVector* pe      = element->NewElementVector(FSvelocityEnum);
-	IssmDouble*    Dstar   = xNew<IssmDouble>((dim*vnumnodes)*(tausize*tnumnodes));
+	IssmDouble*    Dstar   = xNewZeroInit<IssmDouble>((dim*vnumnodes)*(tausize*tnumnodes));
 	IssmDouble*    tau     = xNew<IssmDouble>(tausize*tnumnodes);
 	IssmDouble*    d       = xNew<IssmDouble>(tausize*tnumnodes);
@@ -3478,8 +3478,56 @@
 
 		/*Create Dstar*/
+		/*In dim = 2
+		 *
+		 *       <----------------- tausize ---------------> x tnumnodes
+		 *       |  gamma_ij^x         0       gamma_ij^y  | ^
+		 * Dij = |                                         | dim
+		 *       |     0          gamma_ij^y   gamma_ij^x  | v
+		 *                                                   x
+		 *                                                   vnumnodes
+		 *
+		 *In dim = 3
+		 *
+		 *       |  gamma_ij^x         0          0         gamma_ij^y  gamma_ij^z      0      |
+		 *       |                                                                             |
+		 * Dij = |     0          gamma_ij^y      0         gamma_ij^x     0        gamma_ij^z |
+		 *       |                                                                             |
+		 *       |     0               0      gamma_ij^z        0       gamma_ij^x  gamma_ij^y |
+		 *
+		 * gamma_ij^x = zeta_j dphi_i/dx
+		 *
+		 * where:
+		 *   - zeta_j is the nodal function for the j^th node of the tensor (P1DG)
+		 *   - phi_i  is the nodal function for the i^th node of the velocity (P2)*/
 		element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
 		element->NodalFunctionsTensor(tbasis,gauss);
-
-		_error_("STOP");
+		if(dim==2){
+			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];
+				}
+			}
+		}
+		else{
+			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];
+				}
+			}
+		}
 	}