phis = (fac::interpolate(hUs) & aMesh.Le());

if (pressureFeedback)
{
    pbn = friction->rho()*
        (
          - fac::ndiv(phis, Us)
          + gn*n*h
          - fac::ngrad((pbn&n)*h/(2.*friction->rho()))
        );
}
else
{
    pbn = friction->rho()*
        (
          - fac::ndiv(phis, Us)
          + gn*n*h
        );
}

pb = pbn&n;
//pb = fac::average(fac::interpolate(pb));
pb = Foam::max(pb, dimensionedScalar("pb_min", dimensionSet(1, -1, -2, 0, 0, 0, 0), 0));
