 ChFsiForceExplicitSPH::~ChFsiForceExplicitSPH() {}
 
 //--------------------------------------------------------------------------------------------------------------------------------
-
 void ChFsiForceExplicitSPH::Initialize() {
     ChFsiForce::Initialize();
     cudaMemcpyToSymbolAsync(paramsD, paramsH.get(), sizeof(SimParams));

        }
    }

    if (paramsD.elastic_SPH) {
        if (sum_w_i < paramsD.C_Wi) {
            derivVelRho.w = -1.0;
        } else {
            derivVelRho.w = 1.0;
        }
    }

    if (!paramsD.elastic_SPH) {
        Real nu = paramsD.mu0 / paramsD.rho0;
        Real dvxdt = -preGra.x / rhoPresMuA.x +
                     (velxLap.x + velxGra.x * velxLap.y + velxGra.y * velxLap.z + velxGra.z * velxLap.w) * nu;
        Real dvydt = -preGra.y / rhoPresMuA.x +
                     (velyLap.x + velyGra.x * velyLap.y + velyGra.y * velyLap.z + velyGra.z * velyLap.w) * nu;
        Real dvzdt = -preGra.z / rhoPresMuA.x +
                     (velzLap.x + velzGra.x * velzLap.y + velzGra.y * velzLap.z + velzGra.z * velzLap.w) * nu;
        Real drhodt = -paramsD.rho0 * (velxGra.x + velyGra.y + velzGra.z);

        Real Det_G = (Gi[0] * Gi[4] * Gi[8] - Gi[0] * Gi[5] * Gi[7] - Gi[1] * Gi[3] * Gi[8] + Gi[1] * Gi[5] * Gi[6] +
                      Gi[2] * Gi[3] * Gi[7] - Gi[2] * Gi[4] * Gi[6]);
        Real Det_L = (Li[0] * Li[4] * Li[8] - Li[0] * Li[5] * Li[7] - Li[1] * Li[3] * Li[8] + Li[1] * Li[5] * Li[6] +
                      Li[2] * Li[3] * Li[7] - Li[2] * Li[4] * Li[6]);
        if (rhoPresMuA.w > -1.5 && rhoPresMuA.w < -0.5) {
            if (Det_G > 0.9 && Det_G < 1.1 && Det_L > 0.9 && Det_L < 1.1 && sum_w_i > 0.9) {
                derivVelRho = mR4(dvxdt, dvydt, dvzdt, drhodt);
            }
        }
    }


    return;
}

//--------------------------------------------------------------------------------------------------------------------------------
__global__ void CopySortedToOriginal_XSPH_D(Real3* sortedXSPH,
                                            Real3* originalXSPH,
                                            uint* gridMarkerIndex,
