 			int2 updatePortion = I2(referencePart); //first two component of the referenceArray denote to the fluid and boundary particles
 			thrust::fill(flexIdentifierD.begin() + (updatePortion.x - startFlexMarkers),
 					flexIdentifierD.begin() + (updatePortion.y - startFlexMarkers), flexIdx);
+		}
+	}
+}
+////--------------------------------------------------------------------------------------------------------------------------------
+void Calc_NumNodesMultMarkers_Per_Beam(
+		thrust::device_vector<int> & ANCF_NumNodesMultMarkers_Per_BeamD,
+		const thrust::device_vector<int> & ANCF_NumMarkers_Per_BeamD,
+		const thrust::host_vector<int2> & ANCF_ReferenceArrayNodesOnBeams,
+		int numFlexBodies) {
+	if (numFlexBodies > 0) {
+		for (int flexIdx = 0; flexIdx < numFlexBodies; flexIdx++) {
+			int2 flexPortion = ANCF_ReferenceArrayNodesOnBeams[flexIdx];
+			int numNodes = flexPortion.y - flexPortion.x;
+			ANCF_NumNodesMultMarkers_Per_BeamD[flexIdx] = numNodes * ANCF_NumMarkers_Per_BeamD[flexIdx];
 		}
 	}
 }

	thrust::device_vector<real3> ANCF_VelSlopesD = ANCF_VelSlopes;
	thrust::device_vector<real_> ANCF_Beam_LengthD = ANCF_Beam_Length;
	thrust::device_vector<int2> ANCF_ReferenceArrayNodesOnBeamsD = ANCF_ReferenceArrayNodesOnBeams;
	//*******************
	thrust::device_vector<int> ANCF_NumMarkers_Per_BeamD;
	thrust::device_vector<int> ANCF_NumMarkers_Per_Beam_CumulD;
	thrust::device_vector<int> ANCF_NumNodesMultMarkers_Per_BeamD;
	thrust::device_vector<int> ANCF_NumNodesMultMarkers_Per_Beam_CumulD;

	//*******************


	thrust::device_vector<real3> flexSPH_MeshPos_LRF_D(numFlex_SphMarkers);

	ANCF_VelSlopesD.clear();
	ANCF_Beam_LengthD.clear();
	ANCF_ReferenceArrayNodesOnBeamsD.clear();

	ANCF_NumMarkers_Per_BeamD.clear();
	ANCF_NumMarkers_Per_Beam_CumulD.clear();
	ANCF_NumNodesMultMarkers_Per_BeamD.clear();
	ANCF_NumNodesMultMarkers_Per_Beam_CumulD.clear();



	posRigidCumulativeD.clear();
	velMassRigidD.clear();
