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

	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();
}


}
//------------------------------------------------------------------------------
__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,

		const real3 * ANCF_NodesVelD,
		const real3 * ANCF_SlopesVelD,
		const real3 * ANCF_NodesAccD,
		const real3 * ANCF_SlopesAccD,

		const int2 * ANCF_ReferenceArrayNodesOnBeamsD,
		const bool * ANCF_IsCantileverD
	)
{
	uint i = blockIdx.x * blockDim.x + threadIdx.x;
	if (i >= numFlexBodiesD) return;

	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);
}
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
