 	thrust::device_vector<uint> bodyIndexD=bodyIndex;
 	//thrust::copy(bodyIndex.begin(), bodyIndex.end(), bodyIndexD.begin());
 	thrust::device_vector<float4> derivVelRhoD(mNSpheres);
-
+	printf("a6 yoho\n");
+		//******************** rigid body some initialization
+	thrust::device_vector<int> rigidIdentifierD(0);
+
+	float rigid_SPH_mass;
+	int numRigid_SphParticles = 0;
 	int startRigidParticle = (I2(referenceArray[1])).y;
-	thrust::device_vector<int> rigidIdentifierD(0);
-	//printf("referenceArray.size() %d\n", referenceArray.size());
-	printf("a6 yoho\n");
 	if (referenceArray.size() > 2) {
 		startRigidParticle = (I2(referenceArray[2])).x;
-		int numRigid_SphParticles = referenceArray[2 + numRigidBodies - 1].y - startRigidParticle;
+		numRigid_SphParticles = referenceArray[2 + numRigidBodies - 1].y - startRigidParticle;
 		rigidIdentifierD.resize(numRigid_SphParticles);
 		for (int rigidSphereA = 0; rigidSphereA < numRigidBodies; rigidSphereA++) {
 			int2 updatePortion = I2(referenceArray[2 + rigidSphereA]); //first two component of the referenceArray denote to the fluid and boundary particles
 			thrust::fill(rigidIdentifierD.begin() + (updatePortion.x - startRigidParticle),
 					rigidIdentifierD.begin() + (updatePortion.y - startRigidParticle), rigidSphereA);
 		}
-	}
+
+		//---
+		float4 typicalRigidSPH = mVelMas[referenceArray[2].x];
+		rigid_SPH_mass = typicalRigidSPH.w;
+	} else {
+		float4 dummyFluid = mVelMas[referenceArray[0].x];
+		rigid_SPH_mass = 100 * dummyFluid.w;
+	}
+	cutilSafeCall( cudaMemcpyToSymbolAsync(rigid_SPH_massD, &rigid_SPH_mass, sizeof(rigid_SPH_mass)));
+	cudaMemcpyToSymbolAsync(startRigidParticleD, &startRigidParticle, sizeof(startRigidParticle)); //can be defined outside of the kernel, and only once
+	cudaMemcpyToSymbolAsync(numRigid_SphParticlesD, &numRigid_SphParticles, sizeof(numRigid_SphParticles)); //can be defined outside of the kernel, and only once
+
 	printf("a7 yoho\n");
 	//******************************************************************************
-	int numRigid_SphParticles = rigidIdentifierD.size();
+
 	thrust::device_vector<float3> rigidSPH_MeshPos_LRF_D(numRigid_SphParticles);
-
 	uint nBlocks_numRigid_SphParticles;
 	uint nThreads_SphParticles;
 	computeGridSize(numRigid_SphParticles, 256, nBlocks_numRigid_SphParticles, nThreads_SphParticles);

	}
}
//--------------------------------------------------------------------------------------------------------------------------------
__global__ void SumSurfaceInteractionForces(float3 * totalBodyForces3, float4 * totalSurfaceInteraction4, float4 * velMassRigidD, float rigid_SPH_mass) {
	uint rigidSphereA = blockIdx.x * blockDim.x + threadIdx.x;
	if (rigidSphereA > numRigidBodiesD) {
		return;
	}
	float4 dummyVelMas = velMassRigidD[rigidSphereA];
	totalBodyForces3[rigidSphereA] = rigid_SPH_mass / dummyVelMas.w * F3(totalSurfaceInteraction4[rigidSphereA]);
}
//--------------------------------------------------------------------------------------------------------------------------------
__global__ void CalcTorqueShare(float3* torqueParticlesD, float4* derivVelRhoD, float4* posRadD, int* rigidIdentifierD, float3* posRigidD) {


	thrust::device_vector<float3> totalBodyForces3(numRigidBodies);
	thrust::fill(totalBodyForces3.begin(), totalBodyForces3.end(), F3(0));
	SumSurfaceInteractionForces<<<nBlocks_numRigid_SphParticles, nThreads_SphParticles>>>(F3CAST(totalBodyForces3), F4CAST(totalSurfaceInteraction4), F4CAST(velMassRigidD));
	cudaThreadSynchronize();
	CUT_CHECK_ERROR("Kernel execution failed: SumSurfaceInteractionForces");
	totalSurfaceInteraction4.clear();
