 	}
 	thrust::device_vector<int2> dummyNodesFlexIdentify(totalNumberOfFlexNodes);
 	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(),
-			flex_FSI_NodesForces1.begin(), binary_pred_int2, thrust::plus<real3>());
-	(void) thrust::reduce_by_key(flexMapEachMarkerOnAllBeamNodesD.begin(), flexMapEachMarkerOnAllBeamNodesD.end(), flexNodesForcesAllMarkers2.begin(), dummyNodesFlexIdentify.begin(),
-			flex_FSI_NodesForces2.begin(), binary_pred_int2, thrust::plus<real3>());
+
+	//** Superposes the forces of all BCE markers on nodes onto single force on each node. First 3 nodal coordinates.
+	(void) thrust::reduce_by_key(flexMapEachMarkerOnAllBeamNodesD.begin(), flexMapEachMarkerOnAllBeamNodesD.end(), flexNodesForcesAllMarkers1.begin(),
+			dummyNodesFlexIdentify.begin(), flex_FSI_NodesForces1.begin(),
+			binary_pred_int2, thrust::plus<real3>());
+	//** Superposes the forces of all BCE markers on nodes onto single force on each node. Second 3 nodal coordinates.
+	(void) thrust::reduce_by_key(flexMapEachMarkerOnAllBeamNodesD.begin(), flexMapEachMarkerOnAllBeamNodesD.end(), flexNodesForcesAllMarkers2.begin(),
+			dummyNodesFlexIdentify.begin(), flex_FSI_NodesForces2.begin(),
+			binary_pred_int2, thrust::plus<real3>());
 	dummyNodesFlexIdentify.clear();
 	flexNodesForcesAllMarkers1.clear();
 	flexNodesForcesAllMarkers2.clear();
 
-	//**
-
-
+	//################################################### Update nodal coordinates (integrate in time)
+	//** uses the force values on nodal coordinates to update nodal position and slope and velocities ('2' at the end).
+	//*** it uses the current un-updated nodal coordinates (without '2' at the end) to calculate elastic forces.
 	int n = 1;
 	for (int i = 0; i < n; i++) {
 		Update_ANCF_Beam(

				numFlexBodies, flexParams, dT/n
				);
	}


//	 ....
//	 ....
//	 ....
//	 ....
//	//end

//	//TODO: add paramsH.gravity to Flex objects
//	thrust::device_vector<real3> gravityForces3(numRigidBodies);
//	thrust::fill(gravityForces3.begin(), gravityForces3.end(), paramsH.gravity);
//	thrust::transform(totalAccRigid3.begin(), totalAccRigid3.end(), gravityForces3.begin(), totalAccRigid3.begin(), thrust::plus<real3>());
//	gravityForces3.clear();
//	//

	//################################################### update rigid body things
	computeGridSize(numFlex_SphMarkers, 256, nBlocks_numFlex_SphMarkers, nThreads_SphMarkers);
	UpdateFlexMarkersPosition<<<nBlocks_numFlex_SphMarkers, nThreads_SphMarkers>>>(
			R3CAST(posRadD), R4CAST(velMasD),


	real3 absOmega = Calc_ANCF_Point_Omega(ANCF_NodesVelD, ANCF_SlopesVelD, indexOfClosestNode, sE, lE, rX); //interpolation using ANCF beam, cubic hermit equation
		//ff1
	//	velMasD[absMarkerIndex] = R4(beamPointVel + cross(absOmega, dist3), markerMass); //wrong
	velMasD[absMarkerIndex] = R4(beamPointVel, markerMass);
}
////--------------------------------------------------------------------------------------------------------------------------------
void MakeRigidIdentifier(
