 	Calc_NumNodesMultMarkers_Per_Beam(ANCF_NumNodesMultMarkers_Per_BeamD, ANCF_NumMarkers_Per_BeamD, ANCF_ReferenceArrayNodesOnBeams, numFlexBodies);
 	thrust::exclusive_scan(ANCF_NumNodesMultMarkers_Per_BeamD.begin(), ANCF_NumNodesMultMarkers_Per_BeamD.end(), ANCF_NumNodesMultMarkers_Per_Beam_CumulD());
 
-	Calc_mapEachMarkerOnAllBeamNodes_IdentifierD(ANCF_NumNodesMultMarkers_Per_BeamD, ANCF_NumMarkers_Per_BeamD, ANCF_ReferenceArrayNodesOnBeams, numFlexBodies);
+	Calc_mapEachMarkerOnAllBeamNodes_IdentifierD(flexMapEachMarkerOnAllBeamNodesD, ANCF_NumNodesMultMarkers_Per_Beam_CumulD, ANCF_NumMarkers_Per_BeamD, ANCF_ReferenceArrayNodesOnBeams, numFlexBodies);
 
 	//*******************
 


	int numFlBcRigid = 2 + numRigidBodies;
	int totalNumberOfFlexNodes = ANCF_ReferenceArrayNodesOnBeamsD[ANCF_ReferenceArrayNodesOnBeamsD.size() - 1].y;

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

	thrust::device_vector<real3> flexNodesForcesAllMarkers1();//(totalNumberOfFlexNodes * totalNumberOfFlexMarkers); Size:sum(numNodesOfEachBeam*numSPH_MarkersOfEachBeam)
	thrust::device_vector<real3> flexNodesForcesAllMarkers2();//(totalNumberOfFlexNodes * totalNumberOfFlexMarkers);

	uint nBlocks_numFlex_SphMarkers;
	uint nThreads_SphMarkers;

			R1CAST(flexParametricDistD),
			R4CAST(derivVelRhoD));

	if (flexMapEachMarkerOnAllBeamNodesD.size() != flexNodesForcesAllMarkers1.size()) {
		printf("we have size inconsistency between flex nodesForces and nodesPair identifier");
	}
	thrust::device_vector<int> dummyNodesFlexIdentify(flexMapEachMarkerOnAllBeamNodesD.size());
	thrust::equal_to<int2> binary_pred_int2; //if binary_pred int2 does not work, you have to either add operator == to custom_cutil_math, or you have to map nodes identifiers from int2 to int
	(void) thrust::reduce_by_key(flexMapEachMarkerOnAllBeamNodesD.begin(), flexMapEachMarkerOnAllBeamNodesD.end(), flexNodesForcesAllMarkers1.begin(), dummyNodesFlexIdentify.begin(),
			flexNodesForces1.begin(), binary_pred_int2, thrust::plus<real3>());
	(void) thrust::reduce_by_key(flexMapEachMarkerOnAllBeamNodesD.begin(), flexMapEachMarkerOnAllBeamNodesD.end(), flexNodesForcesAllMarkers2.begin(), dummyNodesFlexIdentify.begin(),
			flexNodesForces2.begin(), binary_pred_int2, thrust::plus<real3>());
	flexNodesForcesAllMarkers1.clear();
	flexNodesForcesAllMarkers2.clear();

//	//TODO: update flex bodies here

//	 ....
//	 ....
//	 ....
