 			Real Wd = W3(d);
 			Real3 velMasB = FETCH(sortedVelMas, j);
 			//			deltaVDenom += mR4(
-			//					velMasB.w / rhoPresMuB.x * mR3(velMasB) * Wd,
-			//					velMasB.w / rhoPresMuB.x * Wd);
+			//					paramsD.markerMass / rhoPresMuB.x * velMasB * Wd,
+			//					paramsD.markerMass / rhoPresMuB.x * Wd);
 			//			deltaVDenom += mR4(
-			//					mR3(velMasB) * Wd,
+			//					velMasB * Wd,
 			//					Wd);
 
 			if (rhoPresMuB.w < -.1) { // only fluid pressure is used to update BCE pressure see Eq 27 of Adami, 2012 paper


	//	//*** Artificial viscosity type 1.3
	//	Real rAB_Dot_GradW = dot(dist3, gradW);
	//	Real3 derivV = -velMasB.w * (rhoPresMuA.y / (rhoPresMuA.x * rhoPresMuA.x) + rhoPresMuB.y / (rhoPresMuB.x *
	// rhoPresMuB.x)) * gradW
	//		+ velMasB.w / (rhoPresMuA.x * rhoPresMuB.x) * 2.0f * paramsD.mu0 * rAB_Dot_GradW / ( d * d +
	// epsilonMutualDistance * paramsD.HSML * paramsD.HSML ) * mR3(velMasA - velMasB);
	//	return mR4(derivV,
	//		rhoPresMuA.x * velMasB.w / rhoPresMuB.x * dot(vel_XSPH_A - vel_XSPH_B, gradW));
}
//--------------------------------------------------------------------------------------------------------------------------------
__device__ inline Real3 DifVelocity_SSI_DEM(const Real3& dist3, const Real& d,
		const Real4& velMasA, const Real4& velMasB) {
	// printf("** DifVelocity_SSI_DEM\n");
	Real l = paramsD.MULT_INITSPACE * paramsD.HSML - d;  // penetration distance
	if (l < 0) {

					+ rhoPresMuB.y / (rhoPresMuB.x * rhoPresMuB.x)) * gradW
			+ paramsD.markerMass * (8.0f * multViscosity) * paramsD.mu0
					* pow(rhoPresMuA.x + rhoPresMuB.x, Real(-2))
					* rAB_Dot_GradW_OverDist * (velMasA - velMasB);
	Real derivRho = rhoPresMuA.x * paramsD.markerMass / rhoPresMuB.x
			* dot(vel_XSPH_A - vel_XSPH_B, gradW);
	//	Real zeta = 0;//.05;//.1;
	//	Real derivRho = rhoPresMuA.x * paramsD.markerMass * invrhoPresMuBx * (dot(vel_XSPH_A - vel_XSPH_B, gradW)
	//			+ zeta * paramsD.HSML * (10 * paramsD.v_Max) * 2 * (rhoPresMuB.x / rhoPresMuA.x - 1) *
	// rAB_Dot_GradW_OverDist
	//			);
