 				real3 vr = R3(velMasB - velMasA);
 				real3 gradW = GradW(dist3);
 
+				// Randles and Libersky, 1996
+				devS3 += -paramsD.mu0 * velMasB.w / rhoPresMuB.x *
+						R3(vr.x * gradW.y + vr.y * gradW.x, vr.x * gradW.z + vr.z * gradW.x, vr.y * gradW.z + vr.z * gradW.y);
+				volS3 += -paramsD.mu0 * velMasB.w / rhoPresMuB.x * 4.0 / 3.0 *
+						R3(vr.x * gradW.x, vr.y * gradW.y, vr.z * gradW.z);
 			}
 		}
 	}

		uint* cellStart,
		uint* cellEnd,
		int numBCE) {
	uint thIndex = blockIdx.x * blockDim.x + threadIdx.x;
	if (thIndex >= numBCE) {
		return;
	}
	//Arman take care of this
	uint BCE_Index = thIndex + min(numObjectsD.startRigidMarkers, numObjectsD.startRigidMarkers); // updatePortionD = [start, end] index of the update portion
	uint sortedIndex = mapOriginalToSorted[BCE_Index]; //index in the sorted array

	// read particle data from sorted arrays

			}
		}
	}

	devStressD[index] = devS3;
	volStressD[index] = volS3;
}
//--------------------------------------------------------------------------------------------------------------------------------
//without normalization
