Assume linear model $$Y = X \beta + \epsilon \\ \epsilon \sim \mathbb{N}(0, \Omega), $$ where $\Omega$ is a known covariance matrix. The GLS estimator for $\beta$ is well-known: $$\widehat{\beta}_{GLS} = (X^T \Omega^{-1} X)^{-1}X^T \Omega^{-1} y$$ Now if we further assume $X$ to be a square matrix and full-rank, we can infer an updated formula for $\widehat{\beta}_{GLS}$:
$$\widehat{\beta}_{GLS} = X^{-1} \Omega X^{-T}X^T \Omega^{-1} y = X^{-1}y$$
And it's just an OLS estimator that does not in any way take the known covariance $\Omega$ into an account.
I naturally would expect $\widehat{\beta}_{GLS}$ to be corrected for $\Omega$ as it is "hardcoded" into the equation of a linear model, but instead, it is the same as $\widehat{\beta}_{OLS}$; the residuals produced by this estimator are zeros. I could argue that there was a parametric multiplier in front of the covariance matrix $\sigma^2 \Omega$ and then conclude that $\sigma = 0$, but, well, there is none. In a sense, I see how it might be ok, as residuals are exactly at the mean of $\epsilon$ and therefore the likelihood is high.
However, I don't think I really understand the intuition behind it and it looks "broken" to me. Is there anything that can be done about it?