$\newcommand{\la}{\lambda}\newcommand{\R}{\mathbb{R}}$Let $I:=(0,1)$ and $Q=I^n$. Let $\la$ denote the Lebesgue measure over $\R^n$. Let $P_j$ denote the partial derivative of $P$ wrt to its $j$th argument. Because the polynomial $P=P(x_1,\dots,x_n)$ is non-constant, at least one of the $P_j$'s is a nonzero polynomial. Without loss of generality, $P_1$ is a nonzero polynomial.
It follows by the main theorem of algebra and the Tonelli theorem that the Lebesgue measure $\la(N)$ of the (closed) set $$N:=\{(x_1,\dots,x_n)\in \R^n\colon P_1(x_1,\dots,x_n)=0\}$$ of the zeros of the nonzero polynomial $P_1$ is $0$.
Consider now the transformation
$$\R^n\ni(x_1,\dots,x_n)\mapsto g(x_1,\dots,x_n):=(P(x_1,\dots,x_n),x_2,\dots,x_n)\in\R^n, $$
the polynomials
\begin{equation}
p_{(y_1,\dots,y_n)}(x_1):=P(x_1,y_2,\dots,y_n)-y_1,
\end{equation}
the sets
\begin{equation}
Z_{(y_1,\dots,y_n)}:=\{x_1\in \R\colon(x_1,y_2,\dots,y_n)\notin N, p_{(y_1,\dots,y_n)}(x_1)=0\},
\end{equation}
for $(y_1,\dots,y_n)\in\R^n$,
the set
\begin{equation}
Y:=\{(y_1,\dots,y_n)\in\R^n\colon p_{(y_1,\dots,y_n)} \text{ has no multiple roots}\},
\end{equation}
and the sets
\begin{equation}
A_m:=\{(y_1,\dots,y_n)\in Y\colon\text{card}\,Z_{(y_1,\dots,y_n)}=m\}
\end{equation}
for $m=0,\dots,d$, where $\text{card}$ denotes the cardinality and $d$ is the degree of $P$ wrt $x_1$. Note that $(A_0,\dots,A_d)$ is a partition of $Y$, and the sets $A_1,\dots,A_m$ are open.
By the main theorem of algebra and the Tonelli theorem, the open set $Y$ is of full Lebesgue measure.
By the inverse function theorem, for each $m\in\{0,\dots,d\}$, each $(y_1,\dots,y_n)\in A_m$, and each integer $k\in[1,m]$, there is a real-analytic function $X_{m,k}$ defined on a neighborhood $U_{(y_1,\dots,y_n)}$ of the point $(y_1,\dots,y_n)\in A_m$ such that \begin{equation} g(X_{m,k}(z_1,\dots,z_n),z_2,\dots,z_n)=(z_1,\dots,z_n) \end{equation} for all $(z_1,\dots,z_n)\in U_{(y_1,\dots,y_n)}$, and the value of the Jacobian determinant of the map \begin{equation} U_{(y_1,\dots,y_n)}\ni(z_1,\dots,z_n)\mapsto(X_{m,k}(z_1,\dots,z_n),z_2,\dots,z_n) \end{equation} at the point $(z_1,\dots,z_n)$ is \begin{equation} \partial_{z_1}X_{m,k}(z_1,\dots,z_n)=\frac1{P_1(X_{m,k}(z_1,\dots,z_n),z_2,\dots,z_n)}. \end{equation}
So, by the change of variables in multiple integrals, the density (say $f$) of the pushforward $g_\#(\la\downharpoonright Q)$ under $g$ of the measure $\la\downharpoonright Q$ wrt to $\la$ is given by the formula \begin{equation} f(y_1,\dots,y_n)=\sum_{m=1}^d \sum_{k=1}^m\frac{1((y_1,\dots,y_n)\in B_{m,k})}{|P_1(X_{m,k}(y_1,\dots,y_n),y_2,\dots,y_n)|} \end{equation} for $(y_1,\dots,y_n)\in\R^n\setminus g(N)$, where \begin{equation} B_{m,k}:=\{(y_1,\dots,y_n)\in A_m\colon(X_{m,k}(y_1,\dots,y_n),y_2,\dots,y_n)\in Q\}. \end{equation} Note that $\la(g(N))=0$, since $\la(N)=0$ and $g$ is Lipschitz. So, we can let $f:=0$ on $g(N)$.
Now, the density (say $f_1$) of the pushforward $P_\#(\la\downharpoonright Q)$ of the measure $\la\downharpoonright Q$ under $P$ wrt to $\la$ is given by the formula \begin{equation} f_1(y_1)=\int_{\R^{n-1}}dy_2\,\cdots\, dy_n\,f(y_1,\dots,y_n) \end{equation} for $y_1\in\R$.