Abstract: Some of the most interesting refraction properties of phononic crystals are revealed by examining the anti-plane shear waves in doubly periodic elastic composites with unit cells containing rectangular and/or elliptical multi-inclusions. The corresponding band structure, group velocity, and energy–flux vector are calculated using a powerful mixed variational method that accurately and efficiently yields all the field quantities over multiple frequency pass-bands. The background matrix and the inclusions can be anisotropic, each having distinct elastic moduli and mass densities. Equifrequency contours and energy–flux vectors are readily calculated as functions of the wave-vector components. By superimposing the energy–flux vectors on equifrequency contours in the plane of the wave-vector components, and supplementing this with a three-dimensional graph of the corresponding frequency surface, a wealth of information is extracted essentially at a glance. This way it is shown that a composite with even a simple square unit cell containing a central circular inclusion can display negative or positive energy and phase velocity refractions, or simply performs a harmonic vibration (standing wave), depending on the frequency and the wave-vector. Moreover, that the same composite when interfaced with a suitable homogeneous solid can display: (1) negative refraction with negative phase velocity refraction; (2) negative refraction with positive phase velocity refraction; (3) positive refraction with negative phase velocity refraction; (4) positive refraction with positive phase velocity refraction; or even (5) complete reflection with no energy transmission, depending on the frequency, and direction and the wavelength of the plane-wave that is incident from the homogeneous solid to the interface. For elliptical and rectangular inclusion geometries, analytical expressions are given for the key calculation quantities. Expressions for displacement, velocity, linear momentum, strain, and stress components, as well as the energy–flux and group velocity components are given in series form. The general results are illustrated for rectangular unit cells, one with two and the other with four inclusions, although any number of inclusions can be considered. The energy–flux and the accompanying phase velocity refractions at an interface with a homogeneous solid are demonstrated. Finally, by comparing the results of the present solution method with those obtained using the Rayleigh quotient and, for the layered case, with the exact solutions, the remarkable accuracy and the convergence rate of the present solution method are demonstrated. Graphical abstract: Anomalous refractive characteristics of phononic crystals are revealed using anti-plane shear waves in doubly periodic elastic composites. It is shown that negative or positive refraction can be accompanied by negative or positive phase velocity refraction, and/or the crystal can have complete refraction with no energy transmission, all depending on the wavelength and frequency of the plane wave that is incident from a homogeneous solid to its interface with the composite. A powerful computational tool is proposed which applies to one-, two-, three-dimensional phononic crystals yielding results with great accuracy. [Figure not available: see fulltext.]