                                  marker_Fp.data_d,      // input
                                  marker_volume.data_d,  // input
                                  output.data_d);
-    //    output.copyDeviceToHost();
-    //
-    //    for (int i = 0; i < host_settings.num_mpm_nodes; i++) {
-    //        printf("Nd: %d [%.20f %.20f %.20f]\n", i, output[i * 3 + 0], output[i * 3 + 1], output[i * 3 + 2]);
-    //    }
 
     kMultiplyB<<<CONFIG(size)>>>(input.data_d, old_vel_node_mpm.data_d, node_mass.data_d, output.data_d);
 }

        {
            LOOP_TWO_RING_GPU(  //
                real3 vnew(v_array[current_node * 3 + 0], v_array[current_node * 3 + 1], v_array[current_node * 3 + 2]);
                real3 vold(old_vel_node_mpm[current_node * 3 + 0], old_vel_node_mpm[current_node * 3 + 1],
                           old_vel_node_mpm[current_node * 3 + 2]);
                real3 v0 = vold + vnew;                                   //
                real3 v1 = dN(xi - current_node_location, inv_bin_edge);  //
                delta_F += OuterProduct(v0, v1);                          //
                )
        }
        Mat33 m_FE = marker_Fe[p];

    //    }

    delta_v.resize(host_settings.num_mpm_nodes * 3);
    delta_v = old_vel_node_mpm;

    MPM_BBSolver(rhs, delta_v);

    kIncrementVelocity<<<CONFIG(host_settings.num_mpm_nodes)>>>(delta_v.data_d, old_vel_node_mpm.data_d,
                                                                grid_vel.data_d);
    kUpdateParticleVelocity<<<CONFIG(host_settings.num_mpm_markers)>>>(
        grid_vel.data_d, old_vel_node_mpm.data_d, pos.data_d, vel.data_d, marker_Fe.data_d, marker_Fp.data_d);
    vel.copyDeviceToHost();
