During the force calculation at the Lagrangian nodes using the formula

F=2(rho)(Ud-Unof)/(dt)

the non-forced velocity (Unof) is interpolated from the surrounding fluid domain (both inside and outside the immersed boundary) and the Ud is the desired velocity. My confusion is where do i get the density (rho) from. Do i interpolate it from the outside or from both sides or do i take an average density of the whole domain?