 			I1CAST(ANCF_NumNodesMultMarkers_Per_Beam_Cumul),
 			R1CAST(parametricDist),
 			R4CAST(derivVelRhoD),
-			markerMass,
-			numFlex_SphMarkers);
+			markerMass);
 
 	if (nodesAndFlexPairIdentifier.size() != flexNodesForcesAllMarkers1.size()) {
 		printf("we have size inconsistency between flex nodesForces and nodesPair identifier");

	derivVelRhoD.clear();
	rigidIdentifierD.clear();
	rigidSPH_MeshPos_LRF_D.clear();
	qD1.clear();
	AD1.clear();
	AD2.clear();

		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
