A neoclassical rotation theory (poloidal and toroidal) is developed from the fluid moment equations, using the Braginskii decomposition of the viscosity tensor extended to generalized curvilinear geometry and a neoclassical calculation of the parallel viscosity coefficient interpolated over collision regimes. Important poloidal dependences are calculated using the Miller equilibrium flux surface geometry representation, which takes into account elongation, triangularity, flux surface compression/expansion and the Shafranov shift. The resulting set of eight (for a two-ion-species plasma model) coupled nonlinear equations for the flux surface averaged poloidal and toroidal rotation velocities and for the up–down and in–out density asymmetries for both ion species are solved numerically. Comparison of prediction with measured carbon poloidal and toroidal rotation velocities in a co-injected and a counter-injected H-mode discharge in DIII-D (Luxon J. 2002 Nucl. Fusion 42 614) indicates agreement to within <10% except in the very edge (ρ > 0.90) in the co-injected discharge.