 
 		const int2 * ANCF_ReferenceArrayNodesOnBeamsD,
 		const bool * ANCF_IsCantileverD,
+		real_ dT,
 		uint i
 	)
 {

		IntegrateInTimeKernel(
				ANCF_NodesD2, ANCF_SlopesD2, ANCF_NodesVelD2, ANCF_SlopesVelD2,
				ANCF_NodesVelD, ANCF_SlopesVelD, ANCF_NodesAccD, ANCF_SlopesAccD,
				ANCF_ReferenceArrayNodesOnBeamsD, ANCF_IsCantileverD, i	);
#if !flexGPU
	}
#endif

		Zero12(e_ee);
		// Elastic Force, 1/2: tension force, GQ 5th order. Maybe 4th order is enough as well.
		for (int k = 0; k < 5; k ++) {
			real_ gqPoint = (lE - 0) / 2 * GQ.GQ5_p[k] + (lE + 0) / 2;
			eps_eps_e(e_ee, gqPoint, lE, e);
			SumElem12(f_e, e_ee, flexParams.E * flexParams.A * (lE - 0) / 2 * GQ.GQ5_w[k]);
		}
		// Elastic Force, 2/2: bending force, GQ 3rd order.
		Elem12 k_ke;
		Zero12(k_ke);
		for (int k = 0; k < 3; k ++) {
			real_ gqPoint = (lE - 0) / 2 * GQ.GQ3_p[k] + (lE + 0) / 2;
			kappa_kappa_e(k_ke, gqPoint, lE, e);
			SumElem12(f_e, k_ke, flexParams.E * flexParams.I  * (lE - 0) / 2 * GQ.GQ3_w[k]);
		}
		// Gravitational Foce
		Elem12 f_g;
		Zero12(f_g);
		gravitational_force(f_g, lE, flexParams.rho, flexParams.A, flexParams.gravity);
		SumElem12(f_e, f_g, -1);
		// Add element forces to associated nodes
		Add_f_ToForces(flex_FSI_NodesForcesD1, flex_FSI_NodesForcesD2, -1, f_e, nodeIdx);
	}
}
//------------------------------------------------------------------------------
__GLOBAL__ void CalcElasticForcesD(
		real3 * flex_FSI_NodesForcesD1,
		real3 * flex_FSI_NodesForcesD2,
		const real3 * ANCF_NodesD,
