                               const Mat33* marker_Fp,      // input
                               const real* marker_volume,   // input
                               const real* node_mass,       // input
+                              Mat33* PolarR,               // input
                               real* grid_vel) {
     const int p = blockIdx.x * blockDim.x + threadIdx.x;
     if (p < device_settings.num_mpm_markers) {

gpu_vector<real> marker_volume;
gpu_vector<real> grid_vel, delta_v;
gpu_vector<real> rhs;
gpu_vector<Mat33> marker_Fe, marker_Fe_hat, marker_Fp, marker_delta_F;
gpu_vector<real> old_vel_node_mpm;
gpu_vector<real> ml, mg, mg_p, ml_p;
gpu_vector<real> dot_g_proj_norm;

    kMultiplyA<<<CONFIG(size)>>>(pos.data_d,    // input
                                 input.data_d,  //
                                 old_vel_node_mpm.data_d,
                                 PolarR.data_d,         // input
                                 PolarS.data_d,         // input
                                 marker_Fe.data_d,      // input
                                 marker_volume.data_d,  // input
                                 output.data_d);

