 					  // sure!
 	Real kD = 40; // 20;//40.0;//20.0; //420.0;				//damping coef. // 40 is good don't change it.
 	Real3 n = dist3 / d;  // unit vector B to A
-	Real m_eff = (velMasA.w * velMasB.w) / (velMasA.w + velMasB.w);
+	Real m_eff = 0.5 * paramsD.markerMass; //(mA * mB) / (mA + mB);
 	Real3 force = (/*pow(paramsD.sizeScale, Real(3)) * */kS * l
-			- kD * m_eff * dot(mR3(velMasA - velMasB), n)) * n; // relative velocity at contact is simply assumed as the relative vel of the centers. If you are
+			- kD * m_eff * dot(velMasA - velMasB, n)) * n; // relative velocity at contact is simply assumed as the relative vel of the centers. If you are
 																// updating the rotation, this should be modified.
-	return force / velMasA.w;  // return dV/dT same as SPH
+	return force / paramsD.markerMass;  // return dV/dT same as SPH
 }
 //--------------------------------------------------------------------------------------------------------------------------------
 __device__ inline Real3 DifVelocity_SSI_Lubrication(const Real3& dist3,
-		const Real& d, const Real4& velMasA, const Real4& velMasB) {
+		const Real& d, const Real3& velMasA, const Real3& velMasB) {
 	// printf("** DifVelocity_SSI_Lubrication\n");
 	Real Delta_c = paramsD.HSML;
 	Real s = d - paramsD.MULT_INITSPACE * paramsD.HSML;


	//	//*** Artificial viscosity type 1.2
	//	Real nu = 22.8f * paramsD.mu0 / 2.0f / (rhoPresMuA.x * rhoPresMuB.x);
	//	Real3 derivV = -velMasB.w * (
	//		rhoPresMuA.y / (rhoPresMuA.x * rhoPresMuA.x) + rhoPresMuB.y / (rhoPresMuB.x * rhoPresMuB.x)
	//		- nu * vAB_Dot_rAB / ( d * d + epsilonMutualDistance * paramsD.HSML * paramsD.HSML )
	//		) * gradW;
	//	return mR4(derivV,
	//		rhoPresMuA.x * velMasB.w / rhoPresMuB.x * dot(vel_XSPH_A - vel_XSPH_B, gradW));

	//*** Artificial viscosity type 2
	Real rAB_Dot_GradW = dot(dist3, gradW);

	// write new velocity back to original unsorted location
	//  uint originalIndex = gridMarkerIndex[index];

	// Real newDensity = densityShare + paramsD.markerMass * W3(0); //?$ include the particle in its summation as well
	// if (rhoPreMuA.w < -.1) { rhoPreMuA.x = newDensity; }
	// rhoPreMuA.y = Eos(rhoPreMuA.x, rhoPreMuA.w);
	//   oldRhoPreMu[originalIndex] = rhoPreMuA;
