 	thrust::copy(ANCF_SlopesVelH2.begin(), ANCF_SlopesVelH2.end(), ANCF_SlopesVelD2.begin());
 	//---------------------------------------------
 
+	flex_NodesForcesD1.clear();
+	flex_NodesForcesD2.clear();
+
 	ANCF_NodesH2.clear();
 	ANCF_SlopesH2.clear();
 	ANCF_NodesVelH2.clear();

}
//------------------------------------------------------------------------------
// Why -= and not += : because M*X2 + K*X = F Therefore, in an explicit method: M*X2 = F - K*X where K*X is our elastic forces (f)
__device__ __host__ inline void Add_f_ToForces(real3 * flex_FSI_NodesForcesD1, real3 * flex_FSI_NodesForcesD2, real_ k, const Elem12 & f, int nodeIdx) {
	flex_FSI_NodesForcesD1[nodeIdx		] += k * R3(f.e00, f.e01, f.e02);
	flex_FSI_NodesForcesD2[nodeIdx		] += k * R3(f.e03, f.e04, f.e05);
	flex_FSI_NodesForcesD1[nodeIdx + 1	] += k * R3(f.e06, f.e07, f.e08);
	flex_FSI_NodesForcesD2[nodeIdx + 1	] += k * R3(f.e09, f.e10, f.e11);
}
////------------------------------------------------------------------------------
//__device__ __host__ inline void MapBeamDataTo_1D_Array(real_* e, const real3 * beamData1, const real3 * beamData2, int2 nodesPortion) { //beamData1: postion, beamData2: slope

	ANCF_NodesVelH.clear();
	ANCF_SlopesVelH.clear();

	flex_NodesForcesH1.clear();
	flex_NodesForcesH2.clear();
	ANCF_ReferenceArrayNodesOnBeamsH.clear();
	ANCF_Beam_LengthH.clear();
	ANCF_IsCantileverH.clear();
