 						flexParams,
 						float(tStep)/stepEnd,
 						delT);
-
-//		ApplyBoundary(posRadD, rhoPresMuD, numAllMarkers, posRigidD, numRigidBodies, ANCF_NodesD, ANCF_ReferenceArrayNodesOnBeamsD, numFlexBodies);
+		ApplyBoundary(posRadD, rhoPresMuD, numAllMarkers, posRigidD, numRigidBodies, ANCF_NodesD, ANCF_ReferenceArrayNodesOnBeamsD, numFlexBodies);
 		//************
 
 
 
 
-//			/* post_process for Segre-Silberberg */ goes before ApplyBoundary
-//			if(tStep >= 0) {
-//				real2 channelCenter = .5 * R2(paramsH.cMax.y + paramsH.cMin.y, paramsH.cMax.z + paramsH.cMin.z);
-//				FindPassesFromTheEnd(posRigidD, distributionD, numRigidBodies, channelCenter, channelRadius, numberOfSections);
-//			}
+					//			/* post_process for Segre-Silberberg */ goes before ApplyBoundary
+					//			if(tStep >= 0) {
+					//				real2 channelCenter = .5 * R2(paramsH.cMax.y + paramsH.cMin.y, paramsH.cMax.z + paramsH.cMin.z);
+					//				FindPassesFromTheEnd(posRigidD, distributionD, numRigidBodies, channelCenter, channelRadius, numberOfSections);
+					//			}
 
 
 

								ANCF_ReferenceArrayNodesOnBeamsD,
								ANCF_NumMarkers_Per_BeamD,
								ANCF_NumMarkers_Per_Beam_CumulD,
						//		thrust::device_vector<int> & ANCF_NumNodesMultMarkers_Per_BeamD,
								ANCF_NumNodesMultMarkers_Per_Beam_CumulD,

								flexIdentifierD,

		thrust::device_vector<real3> ANCF_SlopesVelD2 = ANCF_SlopesVelD;

		//******** RK2
		ForceSPH(posRadD, velMasD, vel_XSPH_D, rhoPresMuD, bodyIndexD, derivVelRhoD, referenceArray, numAllMarkers, paramsH, 0.5 * delT); //?$ right now, it does not consider paramsH.gravity or other stuff on rigid bodies. they should be applied at rigid body solver
		UpdateFluid(posRadD2, velMasD2, vel_XSPH_D, rhoPresMuD2, derivVelRhoD, referenceArray, 0.5 * delT); //assumes ...D2 is a copy of ...D
		//UpdateBoundary(posRadD2, velMasD2, rhoPresMuD2, derivVelRhoD, referenceArray, 0.5 * delT);		//assumes ...D2 is a copy of ...D
		UpdateRigidBody(posRadD2, velMasD2, posRigidD2, posRadRigidCumulativeD2, velMassRigidD2, qD2, AD1_2, AD2_2, AD3_2, omegaLRF_D2, derivVelRhoD, rigidIdentifierD,
				rigidSPH_MeshPos_LRF_D, referenceArray, jD1, jD2, jInvD1, jInvD2, paramsH, numRigidBodies, startRigidMarkers, numRigid_SphMarkers, float(tStep)/stepEnd, 0.5 * delT);
		UpdateFlexibleBody(posRadD2, velMasD2,
				ANCF_NodesD2, ANCF_SlopesD2, ANCF_NodesVelD2, ANCF_SlopesVelD2,
				derivVelRhoD,
