A methodology is developed to integrate \(p\)-adaptive analysis with shape optimization for general velocity fields. The key contributions lie in an easily implementable shape update strategy and in an adaptive scheme to locally increase the polynomial order based on error indicators. Numerical results show that significant shape change is achievable prior to serious mesh distortion. In addition, need for an \(h - p\) adaptive strategy rather than a pure \(p\)-adaptive strategy is also highlighted.