The rheology of highly concentrated monodisperse emulsions is studied by rigorous multidrop numerical simulations for three types of steady macroscopic flow, (i) simple shear ((gamma)over dotx(2), 0 0), (ii) planar extension (PE) ((Gamma)over dotx(1), -(Gamma)over dotx(2), 0) and (iii) mixed ((gamma)over dotx(2), (gamma)over dot chi x(1), 0), where (gamma)over dot and (Gamma)over dot are the deformation rates, and chi is an element of (-1, 1) is the flow parameter, in order to construct and validate a general constitutive model for emulsion flows with arbitrary kinematics. The algorithm is a development of the multipole-accelerated boundary-integral (BI) code of Zinchenko & Davis (J. Fluid Mech., vol. 455, 2002, pp. 21-62). It additionally incorporates periodic boundary conditions for (ii) and (iii) (based on the reproducible lattice dynamics of Kraynik-Reinelt for PE), control of surface overlapping, much more robust controllable surface triangulations for long-time simulations, and more efficient acceleration. The emulsion steady-state viscometric functions (shear viscosity and normal stress differences) for (i) and extensiometric functions (extensional viscosity and stress cross-difference) for (ii) are studied in the range of drop volume fractions c = 0.45-0.55, drop-to-medium viscosity ratios lambda = 0.25-10 and various capillary numbers Ca, with 100-400 drops in a periodic cell and 2000-4000 boundary elements per drop. High surface resolution is important for all three flows at small Ca. Large system size and strains (gamma)over dott of up to several thousand are imperative in some shear-flow simulations to identify the onset of phase transition to a partially ordered state, and evaluate (although still not precisely) the viscometric functions in this state. Below the phase transition point, the shear viscosity versus Ca shows a kinked behaviour, with the local minimum most pronounced at lambda = 1 and c = 0.55. The lambda = 0.25 emulsions flow in a partially ordered manner in a wide range of Ca even when c = 0.45. Increase of lambda to 3-10 shifts the onset of ordering to much smaller Ca, often outside the simulation range. In contrast to simple shear, phase transition is never observed in PE or mixed flow. A generalized five-parameter Oldroyd model with variable coefficients is fitted to our extensiometric and viscometric functions at arbitrary flow intensities (but outside the phase transition range). The model predictions compare very well with precise simulation results for strong mixed flows, chi = 0.25. Time-dependent PE flow is also considered. Ways to overcome the phase transition and drop breakup limitations on constitutive modelling are discussed.