 		thrust::device_vector<real3> & posRadD,
 		thrust::device_vector<real4> & velMasD,
 		const int numRigidBodies,
+		const int numFlexBodies,
+		const int numFlex_SphMarkers,
 		thrust::device_vector<real3> & ANCF_Nodes,
 		thrust::device_vector<real3> & ANCF_Slopes,
 		thrust::device_vector<real3> & ANCF_VelNodes,


	real3_ beamPointVel = ANCF_Point_Vel(ANCF_Nodes, ANCF_Slopes, ANCF_VelNodes, ANCF_VelSlopes, indexOfClosestNode, s, l); //interpolation using ANCF beam, cubic hermit equation

	real3_ r =




	printf("before first kernel\n");
	Populate_RigidSPH_MeshPos_LRF_kernel<<<nBlocks_numRigid_SphMarkers, nThreads_SphMarkers>>>(R3CAST(rigidSPH_MeshPos_LRF_D), R3CAST(posRadD), I1CAST(rigidIdentifierD), R3CAST(posRigidD), startRigidMarkers, numRigid_SphMarkers);
	cudaThreadSynchronize();
	CUT_CHECK_ERROR("Kernel execution failed: CalcTorqueShare");	printf("after first kernel\n");

	//******************************************************************************
	//******************** flex body some initialization

	int numFlBcRigid = 2 + numRigidBodies;
	int numFlexBodies = ANCF_Beam_Length.size();
//	int totalNumberOfFlexNodes = ANCF_ReferenceArrayNodesOnBeams[ANCF_ReferenceArrayNodesOnBeams.size() - 1].y;

	int startFlexMarkers = (referenceArray[numFlBcRigid-1]).y;
	int numFlex_SphMarkers = referenceArray[numFlBcRigid + numFlexBodies - 1].y - startFlexMarkers;
	cudaMemcpyToSymbolAsync(startFlexMarkersD, &startFlexMarkers, sizeof(startFlexMarkers)); //can be defined outside of the kernel, and only once
	cudaMemcpyToSymbolAsync(numFlex_SphMarkersD, &numFlex_SphMarkers, sizeof(numFlex_SphMarkers)); //can be defined outside of the kernel, and only once

		//******************************************************************************
	thrust::device_vector<real3> flexSPH_MeshPos_LRF_D(numFlex_SphMarkers);
	uint nBlocks_numFlex_SphMarkers;
	uint nThreads_SphMarkers;
	computeGridSize(numFlex_SphMarkers, 256, nBlocks_numFlex_SphMarkers, nThreads_SphMarkers);
	printf("before first kernel\n");
	Populate_FlexSPH_MeshPos_LRF_kernel<<<nBlocks_numFlex_SphMarkers, nThreads_SphMarkers>>>//(R3CAST(rigidSPH_MeshPos_LRF_D), R3CAST(posRadD), I1CAST(rigidIdentifierD), R3CAST(posRigidD), startRigidMarkers, numRigid_SphMarkers);
	cudaThreadSynchronize();
	CUT_CHECK_ERROR("Kernel execution failed: CalcTorqueShare");	printf("after first kernel\n");

	//******************************************************************************
	thrust::device_vector<real4> qD1 = mQuatRot;
	thrust::device_vector<real3> AD1(numRigidBodies);
