 		real3 * ANCF_Nodes,
 		real3 * ANCF_Slopes,
 		real3 * ANCF_VelNodes,
-		real3 * ANCF_VelSlopes,
-		int numFlex_SphMarkers) {
-
-	uint index = blockIdx.x * blockDim.x + threadIdx.x;
-	if (index >= numFlex_SphMarkers) {
-		return;
-	}
-	uint flexMarkerIndex = index + startFlexMarkersD; // updatePortionD = [start, end] index of the update portion
+		real3 * ANCF_VelSlopes) {
+
+	uint index = blockIdx.x * blockDim.x + threadIdx.x;
+	if (index >= numFlex_SphMarkersD) {
+		return;
+	}
+	uint absMarkerIndex = index + startFlexMarkersD; // updatePortionD = [start, end] index of the update portion
 	real_ s = parametricDist[index];
 	int flexBodyIndex = flexIdentifierD[index];
 	real_ l = ANCF_Beam_Length[flexBodyIndex];

		SimParams paramsH,
		float fracSimulation,
		real_ dT) {
	int numFlexBodies = ANCF_ReferenceArrayNodesOnBeams.size();
	if (referenceArray.size() < 3 || numFlexBodies < 1) {
		return;
	}
	cudaMemcpyToSymbolAsync(dTD, &dT, sizeof(dT));

	int numFlBcRigid = 2 + numRigidBodies;
	int totalNumberOfFlexNodes = ANCF_ReferenceArrayNodesOnBeams[ANCF_ReferenceArrayNodesOnBeams.size() - 1].y;
	int numFlex_SphMarkers = referenceArray[numFlBcRigid + numFlexBodies - 1].y - referenceArray[numFlBcRigid].x;

	thrust::device_vector<real3> flexNodesForces1();//(totalNumberOfFlexNodes); Size:sum(numNodesOfEachBeam*numSPH_MarkersOfEachBeam)
	thrust::device_vector<real3> flexNodesForces2();

		real3* rigidSPH_MeshPos_LRF_D,
		real3* posRadD,
		int* rigidIdentifierD,
		real3* posRigidD) {
	uint index = blockIdx.x * blockDim.x + threadIdx.x;
	uint rigidMarkerIndex = index + startRigidMarkersD; // updatePortionD = [start, end] index of the update portion
	if (index >= numRigid_SphParticlesD) {
		return;
	}
	real3 dist3 = posRadD[rigidMarkerIndex] - posRigidD[rigidIdentifierD[index]];
	rigidSPH_MeshPos_LRF_D[index] = dist3;
}
//--------------------------------------------------------------------------------------------------------------------------------

__global__ void Populate_FlexSPH_MeshPos_LRF_kernel(
		real3* flexSPH_MeshPos_LRF_D,
		real3 * posRadD,
		int* flexIdentifierD,
		real_* parametricDist,
		real_* ANCF_Beam_Length,
		int* ANCF_NumNodes_Per_Beam,
		real3 * ANCF_Nodes,
		real3 * ANCF_Slopes) {
	uint index = blockIdx.x * blockDim.x + threadIdx.x;
	if (index >= numFlex_SphMarkersD) {
		return;
	}
	uint absMarkerIndex = index + startFlexMarkersD; // updatePortionD = [start, end] index of the update portion
	real_ s = parametricDist[index];
	int flexBodyIndex = flexIdentifierD[index];
	real_ l = ANCF_Beam_Length[flexBodyIndex];
	int nNodes = ANCF_NumNodes_Per_Beam[flexBodyIndex];

	int indexOfClosestNode = int(s / l) * nNodes;
	if (indexOfClosestNode == nNodes) indexOfClosestNode--;

	real3_ beamPointPos = ANCF_Point_Pos(ANCF_Nodes, ANCF_Slopes, indexOfClosestNode, s, l); //interpolation using ANCF beam, cubic hermit equation

	real3_ dist3 = posRadD[absMarkerIndex] - beamPointPos;
	flexSPH_MeshPos_LRF_D[index] = dist3;
}

//--------------------------------------------------------------------------------------------------------------------------------
//the rigid body torque has been calculated in global RF. This kernel maps it to local RF to be appropriate for the formulas
//local torque = T' = A' * T
