 				real4 velMasB = FETCH(sortedVelMas, j);
 				real4 rhoPresMuB = FETCH(sortedRhoPreMu, j);
 
-				//body force in x direction
-				rhoPresMuB.y = (dist3Alpha.x > 0.5 * paramsD.boxDims.x) ? (rhoPresMuB.y - paramsD.bodyForce4.x*paramsD.boxDims.x) : rhoPresMuB.y;
-				rhoPresMuB.y = (dist3Alpha.x < -0.5 * paramsD.boxDims.x) ? (rhoPresMuB.y + paramsD.bodyForce4.x*paramsD.boxDims.x) : rhoPresMuB.y;
-				//body force in x direction
-				rhoPresMuB.y = (dist3Alpha.y > 0.5 * paramsD.boxDims.y) ? (rhoPresMuB.y - paramsD.bodyForce4.y*paramsD.boxDims.y) : rhoPresMuB.y;
-				rhoPresMuB.y = (dist3Alpha.y < -0.5 * paramsD.boxDims.y) ? (rhoPresMuB.y + paramsD.bodyForce4.y*paramsD.boxDims.y) : rhoPresMuB.y;
-				//body force in x direction
-				rhoPresMuB.y = (dist3Alpha.z > 0.5 * paramsD.boxDims.z) ? (rhoPresMuB.y - paramsD.bodyForce4.z*paramsD.boxDims.z) : rhoPresMuB.y;
-				rhoPresMuB.y = (dist3Alpha.z < -0.5 * paramsD.boxDims.z) ? (rhoPresMuB.y + paramsD.bodyForce4.z*paramsD.boxDims.z) : rhoPresMuB.y;
+				modifyPressure(rhoPresMuB, dist3Alpha);
 
 				if (rhoPresMuA.w < 0  ||  rhoPresMuB.w < 0) {
 					if (rhoPresMuA.w == 0) continue;

	// ff1
//	if (rhoPresMuA.w > 0) printf("force value %f %f %f\n", 1e20*derivV.x, 1e20*derivV.y, 1e20*derivV.z);
	return R4(derivV, derivRho);
}
//--------------------------------------------------------------------------------------------------------------------------------
// collide a particle against all other particles in a given cell

		uint* cellStart,
		uint* cellEnd,
		int numBCE) {
	uint thIndex = blockIdx.x * blockDim.x + threadIdx.x;
	if (thIndex >= numBCE) {
		return;
	}
	//Arman take care of this
	uint BCE_Index = thIndex + min(numObjectsD.startRigidMarkers, numObjectsD.startRigidMarkers); // updatePortionD = [start, end] index of the update portion
	uint sortedIndex = mapOriginalToSorted[BCE_Index]; //index in the sorted array

	// read particle data from sorted arrays
	real3 posRadA = FETCH(sortedPosRad, sortedIndex);
	real4 velMasA = FETCH(sortedVelMas, sortedIndex);
	real4 rhoPreMuA = FETCH(sortedRhoPreMu, sortedIndex);

	// get address in grid
	int3 gridPos = calcGridPos(posRadA);

	real3 devS3 = R3(0);
	real3 volS3 = R3(0);

	// examine neighbouring cells
	for (int x = -1; x <= 1; x++) {
		for (int y = -1; y <= 1; y++) {
			for (int z = -1; z <= 1; z++) {
				stressCell(devS3, volS3, gridPos + I3(x, y, z), sortedIndex, posRadA, velMasA, rhoPreMuA, sortedPosRad, sortedVelMas, sortedRhoPreMu,
						cellStart, cellEnd);
			}
		}
	}
}
//--------------------------------------------------------------------------------------------------------------------------------
//without normalization
