 				Real d = length(dist3);
 				if (d > RESOLUTION_LENGTH_MULT * paramsD.HSML) continue;
 				Real4 rhoPresMuB = FETCH(sortedRhoPreMu, j);
-				if (!( rhoPresMuA.w <0 && rhoPresMuB.w < 0 )) continue;//# A and B must be fluid, accoring to colagrossi (2003), the other phase (i.e. rigid) should not be considered)
+				if (rhoPresMuB.w > -.1 ) continue;//# B must be fluid (A was checked originally and it is fluid at this point), accoring to colagrossi (2003), the other phase (i.e. rigid) should not be considered)
 				Real multRho = 2.0f / (rhoPresMuA.x + rhoPresMuB.x);
 				Real4 velMasB = FETCH(sortedVelMas, j);
 				deltaV += velMasB.w * mR3(velMasB - velMasA) * W3(d) * multRho;

	if (index >= numAllMarkers) return;

	// read particle data from sorted arrays
	Real3 posRadA = FETCH(sortedPosRad, index);
	Real4 velMasA = FETCH(sortedVelMas, index);
	Real4 rhoPreMuA = FETCH(sortedRhoPreMu, index);
	Real3 deltaV = mR3(0);

	// get address in grid


	cudaThreadSynchronize();
	CUT_CHECK_ERROR("Kernel execution failed: newVel_XSPH_D");
}
//--------------------------------------------------------------------------------------------------------------------------------
void CalcBCE_Stresses(
