In many parametric problems governed by partial differential equations, like, for example, aerodynamic shape optimization, a key challenge is the efficient evaluation of high-fidelity (HF) CFD models across a design space. Traditional linear reduction methods, built from proper orthogonal decomposition (POD), do not maintain accuracy in the presence of convected shock waves \cite{Lucia}, \cite{Iollo2014}, \cite{Bernard2018}. This is due to the slow decay of the Kolmogorov n-width in parametric shock problems, which limits their applicability in complex scenarios like transonic aircraft optimization. To address this, recent studies propose nonlinear interpolation methods using Lagrangian or Eulerian mappings, which enable the transport of shock features across parametric variations. When combined with a Convex Displacement Interpolation (CDI) framework \cite{Iollo22}, \cite{Cucchiara2024}, these mappings allow to predict the position of the coherent structures that linear approaches cannot capture efficiently. However, existing approaches to compute mappings often lack generality for application to complex 3D geometries and fail to enforce boundary condition compliance efficiently. This work introduces a general method to compute diffeomorphic mappings applicable to arbitrary 3D geometries. These mappings are designed to respect domain boundaries while enabling shock-aligned interpolation for parametric CFD models. We first established the necessary mathematical framework and then implemented a dedicated mapping algorithm suited for complex domains. Its distinctive feature is the optimization of a differential flow aligning distinctive features of the solutions. Finally, we validate the method using CFD simulations of the ONERA M6 wing at varying angles of attack. The results, benchmarked against high-fidelity simulations, demonstrate that the proposed approach maintains good predictive accuracy in cases involving displaced shocks. The prediction was then used in an industrial CFD code as an initialization for the next high-fidelity simulation. Combined with mesh adaptation methods, the prediction improves both the accuracy and computational time of the HF model for the evaluation of new parameters. This framework can improve the feasibility of reduced-order nonlinear modeling for industrial aerodynamic design, offering both accuracy and efficiency.