 
 }
 //------------------------------------------------------------------------------
-__global__ void SolveAndIntegrateInTimeD(
+__global__ void SolveForAccD(
+		real3 * ANCF_NodesAccD,
+		real3 * ANCF_SlopesAccD,
+
+		const real3 * flex_FSI_NodesForcesD1,
+		const real3 * flex_FSI_NodesForcesD2,
+
+		const int2 * ANCF_ReferenceArrayNodesOnBeamsD,
+		const real_ * ANCF_Beam_LengthD,
+		const bool * ANCF_IsCantileverD
+	)
+{
+	uint i = blockIdx.x * blockDim.x + threadIdx.x;
+	if (i >= numFlexBodiesD) return;
+
+	int2 nodesPortion = ANCF_ReferenceArrayNodesOnBeamsD[i];
+	int numNodes = nodesPortion.y - nodesPortion.x;
+	int numElements = numNodes - 1;
+	real_ lBeam = ANCF_Beam_LengthD[i];
+	real_ lE = lBeam / numElements;
+	bool isCantilever = ANCF_IsCantileverD[i];
+	real_ mult = (1/(flexParamsD.rho * flexParamsD.A * lE));
+
+	min_vec(
+			ANCF_NodesAccD,
+			ANCF_SlopesAccD,
+			flex_FSI_NodesForcesD1,
+			flex_FSI_NodesForcesD2,
+			mult,
+			lE,
+			nodesPortion,
+			isCantilever);
+}
+//------------------------------------------------------------------------------
+__global__ void IntegrateInTimeD(
 		real3 * ANCF_NodesD2,
 		real3 * ANCF_SlopesD2,
 		real3 * ANCF_NodesVelD2,
 		real3 * ANCF_SlopesVelD2,
 
-		real3 * ANCF_NodesAccD,
-		real3 * ANCF_SlopesAccD,
-
-		real3 * flex_FSI_NodesForcesD1,
-		real3 * flex_FSI_NodesForcesD2,
 		const real3 * ANCF_NodesVelD,
 		const real3 * ANCF_SlopesVelD,
+		const real3 * ANCF_NodesAccD,
+		const real3 * ANCF_SlopesAccD,
 
 		const int2 * ANCF_ReferenceArrayNodesOnBeamsD,
-		const real_ * ANCF_Beam_LengthD,
 		const bool * ANCF_IsCantileverD
 	)
 {
 	uint i = blockIdx.x * blockDim.x + threadIdx.x;
 	if (i >= numFlexBodiesD) return;
 
-	SolveAndIntegrateInTimeKernel(
-				ANCF_NodesD2, ANCF_SlopesD2, ANCF_NodesVelD2, ANCF_SlopesVelD2,
-				ANCF_NodesAccD, ANCF_SlopesAccD,
-				flex_FSI_NodesForcesD1, flex_FSI_NodesForcesD2, ANCF_NodesVelD, ANCF_SlopesVelD,
-				ANCF_ReferenceArrayNodesOnBeamsD, ANCF_Beam_LengthD, ANCF_IsCantileverD, i);
-}
-//------------------------------------------------------------------------------
-void SolveAndIntegrateInTimeH(
-		real3 * ANCF_NodesD2,
-		real3 * ANCF_SlopesD2,
-		real3 * ANCF_NodesVelD2,
-		real3 * ANCF_SlopesVelD2,
-
-		real3 * ANCF_NodesAccD,
-		real3 * ANCF_SlopesAccD,
-
-		real3 * flex_FSI_NodesForcesD1,
-		real3 * flex_FSI_NodesForcesD2,
-		const real3 * ANCF_NodesVelD,
-		const real3 * ANCF_SlopesVelD,
-
-		const int2 * ANCF_ReferenceArrayNodesOnBeamsD,
-		const real_ * ANCF_Beam_LengthD,
-		const bool * ANCF_IsCantileverD
-	)
-{
-	for (int i = 0; i < numFlexBodiesD; i++) {
-		SolveAndIntegrateInTimeKernel(
-					ANCF_NodesD2, ANCF_SlopesD2, ANCF_NodesVelD2, ANCF_SlopesVelD2,
-					ANCF_NodesAccD, ANCF_SlopesAccD,
-					flex_FSI_NodesForcesD1, flex_FSI_NodesForcesD2, ANCF_NodesVelD, ANCF_SlopesVelD,
-					ANCF_ReferenceArrayNodesOnBeamsD, ANCF_Beam_LengthD, ANCF_IsCantileverD, i);
+	int2 nodesPortion = ANCF_ReferenceArrayNodesOnBeamsD[i];
+	int2 nodesPortionAdjusted2 = nodesPortion;
+	bool isCantilever = ANCF_IsCantileverD[i];
+	if (isCantilever) {
+		nodesPortionAdjusted2.x = nodesPortionAdjusted2.x + 1;
 	}
+
+	ItegrateInTimeKernel(
+		ANCF_NodesD2, ANCF_SlopesD2, ANCF_NodesVelD2, ANCF_SlopesVelD2,
+		ANCF_NodesVelD, ANCF_SlopesVelD, ANCF_NodesAccD, ANCF_SlopesAccD,
+		nodesPortionAdjusted2);
 }
 //------------------------------------------------------------------------------
 //------------------------------------------------------------------------------

	thrust::device_vector<real3> ANCF_SlopesAccD(totalNumberOfFlexNodes);
	thrust::fill(ANCF_NodesAccD.begin(), ANCF_NodesAccD.end(),R3(0));
	thrust::fill(ANCF_SlopesAccD.begin(), ANCF_SlopesAccD.end(),R3(0));
	SolveAndIntegrateInTimeD<<<nBlock_FlexBodies, nThreads_FlexBodies>>>(
			R3CAST(ANCF_NodesD2), R3CAST(ANCF_SlopesD2), R3CAST(ANCF_NodesVelD2), R3CAST(ANCF_SlopesVelD2),
			R3CAST(ANCF_NodesAccD), R3CAST(ANCF_SlopesAccD),
			R3CAST(flex_FSI_NodesForcesD1), R3CAST(flex_FSI_NodesForcesD2), R3CAST(ANCF_NodesVelD), R3CAST(ANCF_SlopesVelD),
			I2CAST(ANCF_ReferenceArrayNodesOnBeamsD), R1CAST(ANCF_Beam_LengthD), BCAST(ANCF_IsCantileverD));
	cudaThreadSynchronize();
	CUT_CHECK_ERROR("Kernel execution failed: SolveAndIntegrateInTime");
	ANCF_NodesAccD.clear();
	ANCF_SlopesAccD.clear();
}

	}
}
//------------------------------------------------------------------------------
__device__ __host__ inline void ItegrateInTimeKernel(
		real3 * ANCF_NodesD2,
		real3 * ANCF_SlopesD2,
		real3 * ANCF_NodesVelD2,
