 		posRigid.z += (paramsD.cMax.z - paramsD.cMin.z);
 		posRigidD[index] = posRigid;
 		return;
+	}
+}
+//--------------------------------------------------------------------------------------------------------------------------------
+//applies periodic BC along x, for ridid bodies
+__global__ void ApplyPeriodicBoundaryKernel_FlexBodies(real3* ANCF_NodesD, int2* ANCF_ReferenceArrayNodesOnBeamsD) {
+	uint index = blockIdx.x * blockDim.x + threadIdx.x;
+	if (index >= numFlexBodiesD) {
+		return;
+	}
+	int2 nodesInterval = ANCF_ReferenceArrayNodesOnBeamsD[index];
+	bool flagX = true, flagY = true, flagZ = true;
+	for (int i = nodesInterval.x; i < nodesInterval.y; i++) {
+		real3 nodePos = ANCF_NodesD[i];
+		if (nodePos.x <= paramsD.cMax.x && nodePos.x >= paramsD.cMin.x) flagX = false;
+		if (nodePos.y <= paramsD.cMax.y && nodePos.y >= paramsD.cMin.y) flagY = false;
+		if (nodePos.z <= paramsD.cMax.z && nodePos.z >= paramsD.cMin.z) flagZ = false;
+	}
+	for (int i = nodesInterval.x; i < nodesInterval.y; i++) {
+		real3 nodePos = ANCF_NodesD[i];
+		if (flagX) nodePos.x = fmod(nodePos.x - paramsD.cMin.x, paramsD.boxDims.x) + paramsD.cMin.x;
+		if (flagY) nodePos.y = fmod(nodePos.y - paramsD.cMin.y, paramsD.boxDims.y) + paramsD.cMin.y;
+		if (flagZ) nodePos.z = fmod(nodePos.z - paramsD.cMin.z, paramsD.boxDims.z) + paramsD.cMin.z;
 	}
 }
 //--------------------------------------------------------------------------------------------------------------------------------

	int nNodes = nodesInterval.y - nodesInterval.x;
	real_ lE = lBeam / (nNodes - 1); //Element length
	real_ sE = fmod(sOverBeam, lE);
	real3 beamPointPos = Calc_ANCF_Point_Pos(ANCF_NodesD, ANCF_SlopesD, indexOfClosestNode, sE, lE); //interpolation using ANCF beam, cubic hermit equation

//		real3 pa = ANCF_NodesD[indexOfClosestNode];
//		real3 pb = ANCF_NodesD[indexOfClosestNode + 1];

		SimParams paramsH,
		const ANCF_Params & flexParams) {
	//****************************** bin size adjustement and contact detection stuff *****************************
	int3 SIDE = I3(int((paramsH.cMax.x - paramsH.cMin.x) / paramsH.binSize0 + .1), int((paramsH.cMax.y - paramsH.cMin.y) / paramsH.binSize0 + .1),
			int((paramsH.cMax.z - paramsH.cMin.z) / paramsH.binSize0 + .1));
	real_ mBinSize = paramsH.binSize0; //Best solution in that case may be to change cMax or cMin such that periodic sides be a multiple of binSize

	printf("SIDE: %d, %d, %d\n", SIDE.x, SIDE.y, SIDE.z);
	//**********************************************************************************************************
	paramsH.gridSize = SIDE;
	//paramsH.numCells = SIDE.x * SIDE.y * SIDE.z;
	paramsH.worldOrigin = paramsH.cMin;
	paramsH.cellSize = R3(mBinSize, mBinSize, mBinSize);
	paramsH.boxDims = paramsH.cMax - paramsH.cMin;
	printf("boxDims: %f, %f, %f\n", paramsH.boxDims.x, paramsH.boxDims.y, paramsH.boxDims.z);

	setParameters(&paramsH); 														// sets paramsD in SDKCollisionSystem
