High-order spectral/hp simulations are of practical interest since they are characterized by small diffusion and dispersion errors and converge exponentially for smooth flows. However, they exhibit high-frequency oscillations (Gibbs phenomenon) in the vicinity of discontinuities such as shocks. To eliminate these oscillations and stabilize the solution, they require the addition of high-order numerical dissipation terms. Essentially, this means that the extra degrees of freedom, afforded by using a high-order approximation, are wasted at the discontinuities. All this suggests that a more efficient approach is the use of p-adaptation, where we use a low-order approximation near discontinuities and high-order ones elsewhere to take advantage of the spectral convergence of the method in the regions of smooth flow. Moreover, we propose to control the resolution in the vicinity of discontinuities, while keeping the polynomial order low, with a combination of h-adaptation via mesh refinement and r-adaptation via mesh deformation, using a variational approach. Our ultimate goal is to combine the best properties of all three mesh adaptation techniques to devise an optimal h-r-p adaptation strategy. We will explore different combinations of the adaptation techniques and assess their individual cost, implications on the solver and overall cost by simulating flows for simple test geometries (a jet intake and a NACA0012 airfoil). Finally, we will present a proposal for devising an optimal h-r-p adaptation capable of accurately resolving the flow structures starting from a uniform under-resolved mesh. The proof-of-concept strategies and the individual mesh modification techniques are implemented and available as part of the open-source spectral/hp element framework Nektar++ and its high-order mesh generator NekMesh.