            host_settings.bin_edge, host_settings.num_mpm_nodes, host_settings.num_mpm_markers);
 }
 //
-void MPM_Initialize(MPM_Settings& settings, std::vector<real3>& positions) {
-    host_settings = settings;
-
-    cudaCheck(cudaMalloc(&lower_bound, sizeof(real3)));
-    cudaCheck(cudaMalloc(&upper_bound, sizeof(real3)));
-
-    pos.data_h = positions;
-    pos.copyHostToDevice();
-
-    cudaCheck(cudaMemcpyToSymbolAsync(device_settings, &host_settings, sizeof(MPM_Settings)));
-
-    MPM_ComputeBounds();
-    marker_volume.resize(host_settings.num_mpm_markers);
-    node_mass.resize(host_settings.num_mpm_nodes);
-    node_mass = 0;
-
-    cudaCheck(cudaPeekAtLastError());
-    cudaCheck(cudaDeviceSynchronize());
-
-    kRasterize<<<CONFIG(host_settings.num_mpm_markers)>>>(pos.data_d,         // input
-                                                          node_mass.data_d);  // output
-
-    cudaCheck(cudaPeekAtLastError());
-    cudaCheck(cudaDeviceSynchronize());
-
-    kComputeParticleVolumes<<<CONFIG(host_settings.num_mpm_markers)>>>(pos.data_d,        // input
-                                                                       node_mass.data_d,  // output
-                                                                       marker_volume.data_d);
-
-    cudaCheck(cudaPeekAtLastError());
-    cudaCheck(cudaDeviceSynchronize());
-}
+
 void Multiply(gpu_vector<real>& input, gpu_vector<real>& output, gpu_vector<real>& r) {
     //    int size = input.size();
     //    kMultiplyA<<<CONFIG(size)>>>(sorted_pos,  // input


void MPM_Solve(MPM_Settings& settings, std::vector<real3>& positions, std::vector<real3>& velocities) {
    host_settings = settings;
    MPM_ComputeBounds();
    //
    //    // ========================================================================================
    //    kRasterize<<<CONFIG(num_mpm_markers)>>>(pos.data_d,        // input
    //                                            vel.data_d,        // input
    //                                            node_mass.data_d,  // output
    //                                            grid_vel.data_d    // output
    //                                            );
    //
    //    kNormalizeWeights<<<CONFIG(num_mpm_nodes)>>>(node_mass.data_d,  // output
    //                                                 grid_vel.data_d);
    //
    //    old_vel_node_mpm = grid_vel;
    //
    //    kRhs<<<CONFIG(num_mpm_markers)>>>(pos.data_d,            // input
    //                                      marker_Fe_hat.data_d,  // input
    //                                      marker_Fe.data_d,      // input

        LOOP_TWO_RING_GPU(                                              //
            real weight = N(xi - current_node_location, inv_bin_edge);  //
            particle_density += grid_mass[current_node] * weight;       //
            )
        // Inverse density to remove division
        particle_density = (bin_edge * bin_edge * bin_edge) / particle_density;
        volume[p] = device_settings.mass * particle_density;
    }
}

