 					derivRho += derivVelRho.w;
 				}
 				else if (fabs(rhoPresMuA.w - rhoPresMuB.w) > 0) { //implies: one of them is solid/boundary, ther other one is solid/boundary of different type or different solid
-//					derivV += DifVelocity_SSI_DEM(dist3, d, rSPH, velMasA, velMasB);
-					real3 dV = DifVelocity_SSI_DEM(dist3, d, rSPH, velMasA, velMasB);
+//					real3 dV = DifVelocity_SSI_DEM(dist3, d, rSPH, velMasA, velMasB);
+					real3 dV = DifVelocity_SSI_Lubrication(dist3, d, rSPH, velMasA, velMasB);
+
 
 					if (rhoPresMuA.w > 0 && rhoPresMuA.w <= numObjectsD.numRigidBodies) { //i.e. rigid
 						uint originalIndex = gridMarkerIndex[index];

	real_ m_eff = (velMasA.w * velMasB.w) / (velMasA.w + velMasB.w);
	real3 force = (/*pow(paramsD.sizeScale, 3) * */kS * l - kD * m_eff * dot(R3(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
}
//--------------------------------------------------------------------------------------------------------------------------------
// collide a particle against all other particles in a given cell

				const real4 & velMasA,
				const real4 & velMasB) {
//printf("** DifVelocity_SSI_DEM\n");
	real_ l = paramsD.MULT_INITSPACE * rSPH - d; //penetration distance
	if (l < 0) {
		return R3(0);
	}
