An extended 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 of density and velocity 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. The numerical solution methodology, a combination of nonlinear Successive Over-Relaxation(SOR) and Simulated Annealing(SA), is also discussed. Comparison of prediction with measured carbon poloidal and toroidal rotation velocities in a co-injected and a counter-injected H-mode discharges in DIII-D [J. Luxon, Nucl. Fusion 42, 614 (2002)] indicates agreement to within <10% except in the very edge in the co-injected discharge.