                                                             marker_volume.data_d,  // input
                                                             node_mass.data_d,      // input
                                                             PolarR.data_d,         // output
+                                                            PolarS.data_d,         // output
                                                             grid_vel.data_d);      // output
 
     kRhs<<<CONFIG(host_settings.num_mpm_nodes)>>>(node_mass.data_d, grid_vel.data_d, rhs.data_d);

        }
    }
}
CUDA_GLOBAL void kSVD(Mat33* marker_Fe_hat, Mat33* PolarR, Mat33* PolarS) {
    const int p = blockIdx.x * blockDim.x + threadIdx.x;
    if (p < device_settings.num_mpm_markers) {
        Mat33 U, V, R, S, W;
        real3 E;
        SVD(marker_Fe_hat[p], U, E, V);
        // Perform polar decomposition F = R*S
        R = MultTranspose(U, V);
        S = V * MultTranspose(Mat33(E), V);

        PolarR[p] = R;
        PolarS[p] = S;
    }
}
CUDA_GLOBAL void kMultiplyA(const real3* sorted_pos,  // input
                            const real* v_array,
                            const real* old_vel_node_mpm,

        marker_Fe_hat[p] = Fe_hat_t * marker_Fe[p];
    }
}

// CUDA_GLOBAL void kSVD(Mat33* marker_Fe_hat, Mat33* PolarR, Mat33* PolarS) {
//    const int p = blockIdx.x * blockDim.x + threadIdx.x;
//    if (p < device_settings.num_mpm_markers) {
//        Mat33 U, V, R, S, W;
//        real3 E;
//        SVD(marker_Fe_hat[p], U, E, V);
//        // Perform polar decomposition F = R*S
//        R = MultTranspose(U, V);
//        S = V * MultTranspose(Mat33(E), V);
//
//        PolarR[p] = R;
//        PolarS[p] = S;
//    }
//}

CUDA_GLOBAL void kApplyForces(const real3* sorted_pos,     // input
                              const Mat33* marker_Fe_hat,  // input
                              const Mat33* marker_Fe,      // input
