 __host__ void ChSystemGpuMesh_impl::runTriangleBroadphase() {
     METRICS_PRINTF("Resetting broadphase info!\n");
 
-    packSphereDataPointers();
-
     std::vector<unsigned int, cudallocator<unsigned int>> Triangle_NumSDsTouching;
     std::vector<unsigned int, cudallocator<unsigned int>> Triangle_SDsCompositeOffsets;
 
-    Triangle_NumSDsTouching.resize(meshSoup->nTrianglesInSoup, 0);
-    Triangle_SDsCompositeOffsets.resize(meshSoup->nTrianglesInSoup, 0);
-
-    const int nthreads = CUDA_THREADS_PER_BLOCK;
-    int nblocks = (meshSoup->nTrianglesInSoup + nthreads - 1) / nthreads;
-    triangleSoup_CountSDsTouched<<<nblocks, nthreads>>>(meshSoup, Triangle_NumSDsTouching.data(), gran_params,
-                                                        tri_params);
+    unsigned int numTriangles = meshSoup->nTrianglesInSoup;
+    Triangle_NumSDsTouching.resize(numTriangles, 0);
+    Triangle_SDsCompositeOffsets.resize(numTriangles, 0);
+
+    int nblocks = (numTriangles + CUDA_THREADS_PER_BLOCK - 1) / CUDA_THREADS_PER_BLOCK;
+    triangleSoup_CountSDsTouched<<<nblocks, CUDA_THREADS_PER_BLOCK>>>(meshSoup, Triangle_NumSDsTouching.data(),
+                                                                      gran_params, tri_params);
 
     gpuErrchk(cudaDeviceSynchronize());
     gpuErrchk(cudaPeekAtLastError());
 
-    unsigned int numTriangles = meshSoup->nTrianglesInSoup;
-    unsigned int num_entries = 0;
-
     // do prefix scan
-    {
-        void* d_temp_storage = NULL;
-        size_t temp_storage_bytes = 0;
-        unsigned int* out_ptr = Triangle_SDsCompositeOffsets.data();
-        unsigned int* in_ptr = Triangle_NumSDsTouching.data();
-
-        // copy data into the tmp array
-        gpuErrchk(cudaMemcpy(out_ptr, in_ptr, numTriangles * sizeof(unsigned int), cudaMemcpyDeviceToDevice));
-        cub::DeviceScan::ExclusiveSum(d_temp_storage, temp_storage_bytes, in_ptr, out_ptr, numTriangles);
-
-        gpuErrchk(cudaDeviceSynchronize());
-        gpuErrchk(cudaPeekAtLastError());
-        // Allocate temporary storage
-        gpuErrchk(cudaMalloc(&d_temp_storage, temp_storage_bytes));
-
-        gpuErrchk(cudaDeviceSynchronize());
-        gpuErrchk(cudaPeekAtLastError());
-        // Run exclusive prefix sum
-        cub::DeviceScan::ExclusiveSum(d_temp_storage, temp_storage_bytes, in_ptr, out_ptr, numTriangles);
-
-        gpuErrchk(cudaDeviceSynchronize());
-        gpuErrchk(cudaPeekAtLastError());
-        gpuErrchk(cudaFree(d_temp_storage));
-        num_entries = out_ptr[numTriangles - 1] + in_ptr[numTriangles - 1];
-        // printf("%u entries total!\n", num_entries);
-    }
+    size_t temp_storage_bytes = 0;
+    unsigned int* out_ptr = Triangle_SDsCompositeOffsets.data();
+    unsigned int* in_ptr = Triangle_NumSDsTouching.data();
+
+    // copy data into the tmp array
+    gpuErrchk(cudaMemcpy(out_ptr, in_ptr, numTriangles * sizeof(unsigned int), cudaMemcpyDeviceToDevice));
+    cub::DeviceScan::ExclusiveSum(NULL, temp_storage_bytes, in_ptr, out_ptr, numTriangles);
+    gpuErrchk(cudaDeviceSynchronize());
+
+    // get pointer to device memory; this memory block will be used internally by CUB, for scratch area
+    void* d_scratch_space = (void*)stateOfSolver_resources.pDeviceMemoryScratchSpace(temp_storage_bytes);
+    // Run exclusive prefix sum
+    cub::DeviceScan::ExclusiveSum(d_scratch_space, temp_storage_bytes, in_ptr, out_ptr, numTriangles);
+    gpuErrchk(cudaDeviceSynchronize());
+    unsigned num_entries = out_ptr[numTriangles - 1] + in_ptr[numTriangles - 1];
 
     // for (unsigned int i = 0; i < Triangle_NumSDsTouching.size(); i++) {
     //     printf("Triangle %u touches %u SDs, offset is %u\n", i, Triangle_NumSDsTouching[i],

        if (meshSoup->nTrianglesInSoup != 0 && mesh_collision_enabled) {
            runTriangleBroadphase();
        }
        gpuErrchk(cudaPeekAtLastError());
        gpuErrchk(cudaDeviceSynchronize());

        if (meshSoup->numTriangleFamilies != 0 && mesh_collision_enabled) {
            // TODO please do not use a template here

        }

        runSphereBroadphase();
        elapsedSimTime += (float)(stepSize_SU * TIME_SU2UU);  // Advance current time
    }

