We study the problem of the nonparametric estimation for the density $π$ of the stationary distribution of a $d$-dimensional stochastic differential equation $(X_t)_{t \in [0, T]}$. From the continuous observation of the sampling path on $[0, T]$, we study the estimation of $π(x)$ as $T$ goes to infinity. For $d\ge2$, we characterize the minimax rate for the $\mathbf{L}^2$-risk in pointwise estimation over a class of anisotropic Hölder functions $π$ with regularity $β= (β_1, ... , β_d)$. For $d \ge 3$, our finding is that, having ordered the smoothness such that $β_1 \le ... \le β_d$, the minimax rate depends on whether $β_2 < β_3$ or $β_2 = β_3$. In the first case, this rate is $(\frac{\log T}{T})^γ$, and in the second case, it is $(\frac{1}{T})^γ$, where $γ$ is an explicit exponent dependent on the dimension and $\barβ_3$, the harmonic mean of smoothness over the $d$ directions after excluding $β_1$ and $β_2$, the smallest ones. We also demonstrate that kernel-based estimators achieve the optimal minimax rate. Furthermore, we propose an adaptive procedure for both $L^2$ integrated and pointwise risk. In the two-dimensional case, we show that kernel density estimators achieve the rate $\frac{\log T}{T}$, which is optimal in the minimax sense. Finally we illustrate the validity of our theoretical findings by proposing numerical results.