 	}
 	cudaMemcpyToSymbolAsync(dTD, &dT, sizeof(dT));
 
+	//################################################### make force arrays
 	int numFlBcRigid = 2 + numRigidBodies;
 	int2 totalNumberOfFlexNodes2 = ANCF_ReferenceArrayNodesOnBeamsD[ANCF_ReferenceArrayNodesOnBeamsD.size() - 1];
 	int totalNumberOfFlexNodes = totalNumberOfFlexNodes2.y;
-	int totalNumberOfFlexMultNodes = flexMapEachMarkerOnAllBeamNodesD.size();
+	int totalNumberOfFlexBCEMultNodes = flexMapEachMarkerOnAllBeamNodesD.size();
 
 	thrust::device_vector<real3> flex_FSI_NodesForces1(totalNumberOfFlexNodes);
 	thrust::device_vector<real3> flex_FSI_NodesForces2(totalNumberOfFlexNodes);
 	thrust::fill(flex_FSI_NodesForces1.begin(), flex_FSI_NodesForces1.end(),R3(0));
 	thrust::fill(flex_FSI_NodesForces2.begin(), flex_FSI_NodesForces2.end(),R3(0));
 
-
-						//ff1 for Radu
-					//	for (int i = 0; i < ANCF_Beam_LengthD.size(); i++) {
-					//		int2 nodesInterval = ANCF_ReferenceArrayNodesOnBeamsD[i];
-					//		real_ lBeam = ANCF_Beam_LengthD[i];
-					//		real3 pa, pb;
-					//		pa = ANCF_NodesD[nodesInterval.x];
-					//		pb = ANCF_NodesD[nodesInterval.y - 1];
-					//		real3 sa, sb;
-					//		sa = ANCF_SlopesD[nodesInterval.x];
-					//		sb = ANCF_SlopesD[nodesInterval.y - 1];
-					//		printf("beamNumber %d: beamLength %f, start_pointAndSlope %f %f %f %f %f %f, end_point %f %f %f\n\n", i, lBeam, pa.x, pa.y, pa.z, sa.x, sa.y, sa.z, pb.x, pb.y, pb.z);
-					//
-					//	}
-
-
 	//**
-	thrust::device_vector<real3> flexNodesForcesAllMarkers1(totalNumberOfFlexMultNodes);
-	thrust::device_vector<real3> flexNodesForcesAllMarkers2(totalNumberOfFlexMultNodes);
+	thrust::device_vector<real3> flexNodesForcesAllMarkers1(totalNumberOfFlexBCEMultNodes);
+	thrust::device_vector<real3> flexNodesForcesAllMarkers2(totalNumberOfFlexBCEMultNodes);
 
 	uint nBlocks_numFlex_SphMarkers;
 	uint nThreads_SphMarkers;
 	computeGridSize(numFlex_SphMarkers, 256, nBlocks_numFlex_SphMarkers, nThreads_SphMarkers);
+	//** maps each BCE marker FORCE onto all beam nodes. Does not sum them up, writes them next to each other. Assumes the BCE forces
+	//*** at the beam center line. "flexNodesForcesAllMarkers1" and "flexNodesForcesAllMarkers1" are built.
 	MapForcesOnNodes<<<nBlocks_numFlex_SphMarkers, nThreads_SphMarkers>>>(
 			R3CAST(flexNodesForcesAllMarkers1),
 			R3CAST(flexNodesForcesAllMarkers2),


	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);
		velMasD[absMarkerIndex] = R4(beamPointVel, markerMass);
}
////--------------------------------------------------------------------------------------------------------------------------------
void MakeRigidIdentifier(

				numFlexBodies, flexParams, dT/n
				);
	}
	//################################################### update BCE markers position
	//** new nodal velocity and positions are used to update BCE markers position and velocities. "posRadD", "velMasD" are updated (only the flex portion)
	computeGridSize(numFlex_SphMarkers, 256, nBlocks_numFlex_SphMarkers, nThreads_SphMarkers);
	UpdateFlexMarkersPosition<<<nBlocks_numFlex_SphMarkers, nThreads_SphMarkers>>>(
			R3CAST(posRadD), R4CAST(velMasD),
