 				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;

		thrust::device_vector<uint>  & cellEnd,
		uint numAllMarkers,
		uint numCells) {
	//#if USE_TEX
	//    cutilSafeCall(cudaBindTexture(0, oldPosTex, sortedPosRad, numAllMarkers*sizeof(Real4)));
	//    cutilSafeCall(cudaBindTexture(0, oldVelTex, sortedVelMas, numAllMarkers*sizeof(Real4)));
	//    cutilSafeCall(cudaBindTexture(0, cellStartTex, cellStart, numCells*sizeof(uint)));
	//    cutilSafeCall(cudaBindTexture(0, cellEndTex, cellEnd, numCells*sizeof(uint)));    
	//#endif

	// thread per particle
	uint numThreads, numBlocks;


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