 	GaussQuadrature GQ;
 	IntitializeGaussQuadrature(GQ);
 	//---------------------------------------------
+#if flexGPU
 	cutilSafeCall( cudaMemcpyToSymbolAsync(flexParamsD, &flexParams, sizeof(ANCF_Params)));
 	cutilSafeCall( cudaMemcpyToSymbolAsync(dTD, &dT, sizeof(real_)));
 	cutilSafeCall( cudaMemcpyToSymbolAsync(numFlexBodiesD, &numFlexBodies, sizeof(int))); //This can be updated in the future, to work with numObjectsD
 	cutilSafeCall( cudaMemcpyToSymbolAsync(GQD, &GQ, sizeof(GaussQuadrature))); //This can be updated in the future, to work with numObjects
-
+#else
+	flexParamsD = flexParams;
+	dTD = dT;
+	numFlexBodiesD = numFlexBodies;
+	GQD = GQ;
+#endif
 	//---------------------------------------------
 	//####### Calculate Forces
 	uint nBlock_FlexBodies;
 	uint nThreads_FlexBodies;
 	computeGridSize(numFlexBodies, 128, nBlock_FlexBodies, nThreads_FlexBodies);
 
-	CalcElasticForces<<<nBlock_FlexBodies, nThreads_FlexBodies>>>(
+	CalcElasticForcesD<<<nBlock_FlexBodies, nThreads_FlexBodies>>>(
 			R3CAST(flex_FSI_NodesForcesD1), R3CAST(flex_FSI_NodesForcesD2),
 			R3CAST(ANCF_NodesD), R3CAST(ANCF_SlopesD), R3CAST(ANCF_NodesVelD), R3CAST(ANCF_SlopesVelD),
 			I2CAST(ANCF_ReferenceArrayNodesOnBeamsD), R1CAST(ANCF_Beam_LengthD));
 	cudaThreadSynchronize();
-	CUT_CHECK_ERROR("Kernel execution failed: CalcElasticForces");
+	CUT_CHECK_ERROR("Kernel execution failed: CalcElasticForcesD");
 
 	int totalNumberOfFlexNodes = flex_FSI_NodesForcesD1.size();
 	thrust::device_vector<real3> ANCF_NodesAccD(totalNumberOfFlexNodes);

//		e[6 * j + 5] = data2.z;
//	}
//}
//------------------------------------------------------------------------------
__device__ __host__ inline void ItegrateInTimeKernel(
		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,

		int2 nodesPortionAdjusted
		)
{
	int numNodesAdjusted = nodesPortionAdjusted.y - nodesPortionAdjusted.x;
	for (int j = 0; j < numNodesAdjusted; j++) {
		int nodeIdx = nodesPortionAdjusted.x + j;
		ANCF_NodesD2[nodeIdx] += dTD * ANCF_NodesVelD[nodeIdx];
		ANCF_SlopesD2[nodeIdx] += dTD * ANCF_SlopesVelD[nodeIdx];

		ANCF_NodesVelD2[nodeIdx] += dTD * ANCF_NodesAccD[nodeIdx];
		ANCF_SlopesVelD2[nodeIdx] += dTD * ANCF_SlopesAccD[nodeIdx];
	}
}

//------------------------------------------------------------------------------
__global__ void CalcElasticForces(
		real3 * flex_FSI_NodesForcesD1,
		real3 * flex_FSI_NodesForcesD2,
		const real3 * ANCF_NodesD,


#include <cstdio>

#define flexGPU true
#if flexGPU
#define flexParamsD flexParamsD
#define dTD dTD
#define numFlexBodiesD numFlexBodiesD
#define GQD GQD
#else
#define flexParamsD flexParamsH
#define dTD dTH
#define numFlexBodiesD numFlexBodiesH
#define GQD GQH
#endif

__constant__ ANCF_Params flexParamsD;
__constant__ real_ dTD;
__constant__ int numFlexBodiesD;
__constant__ GaussQuadrature GQD;

ANCF_Params flexParamsH;
real_ dTH;
int numFlexBodiesH;
GaussQuadrature GQH;

//------------------------------------------------------------------------------
__device__ __host__ inline void Zero12(Elem12 & m) {
