Introduction
I'm not really an expert at these functional equations, but it seems that I'm on a bit of a streak (two of my previous three answers come from this tag), and as everybody knows: being in flow is a wonderful thing.
The first thing I did was check Approach0 and the likes for similar questions. I didn't find the exact question anywhere, but I must admit that I found some really similar questions. However, similarity can easily count for nothing, especially with olympiad problems. That's why I'm not going to bother citing those similar problems.
I will definitely try and explain my thought process every step of the way so that others can also understand how to attack similar problems. I will also, briefly, step down a blind alley to illustrate a useful-looking method that I believe will not work for this scenario.
First thoughts and notation
When I hit this problem for the first time, I made some observations that I have here on record.
Let $P(x,y)$ be the assertion $$f(xf(y)+f(x)+y) = xy+f(x)+f(y)$$ for $x,y \in \mathbb R$. I observed that
The right-hand side is symmetric in $x,y$.
The left hand side is $f(\cdot)$ for some $\cdot$. Therefore, every RHS produced by $P(x,y)$ must be in the range of $f$.
On $f(0)$ and the range of $f$
Then I began with the obvious substitutions
$$
P(x,0) : f(xf(0)+f(x)) = f(x)+f(0)
$$
$$
P(0,y) : f(f(0)+y) = f(0)+f(y)
$$
$$
P(x,-f(x)) : f(xf(-f(x))) = -xf(x) + f(x) - f(-f(x))
$$
The third one is useful : using $x=0$, $P(0,-f(0))$ gives $f(-f(0)) = 0$ following cancellations.
On the uniqueness of $f(y)=0$, and the range of $f$
As the OP mentioned, $0$ has the unique preimage $0$.
Indeed, if $f(y)=0$ then $P(y,y)$ gives $y^2=0$ so that $y=0$.
We derived earlier that $f(-f(0)) = 0$. Therefore, it follows that $-f(0)=0$, hence $f(0)=0$ and $0$ is the unique such value by the previous paragraph.
Once we get this, look back at $P(x,0)$ and observe that $f(0)=0$ to see that $f(f(x))=f(x)$ for all $x$.
The range of $f$, and a blind alley
We saw that $f(f(x)) = f(x)$ is true for all $x$. We can be tempted to show that $f$ is surjective now. In this direction, I'll explain what I did, because the commenters above, who very eagerly attempted to solve this question, will be happy to know that I followed them and treated them like my gurus.
Let $x,y$ be arbitrary. Start with $P(f(x),f(y))$ and use $f(f(x))$ to get $$
f(f(x)f(y)+f(x)+f(y)) = f(x)f(y)+f(x)+f(y)
$$
That is, we've proven that under the function $g(x,y)=xy+x+y$, the range of $f$ is closed. This can be used to derive various corollaries from the comments.
Another observation can be made by plugging in $P(-1,-1)$ which gives $f(-1)=-1$. (THIS will be truly crucial)
One can also try to prove the following : $f(f(x)^2)= f(x)^2$ for all $x$ , and $f(-f(x)) = -f(x)$ for all $x$. Thus, we've shown that the range is closed under various operations.
However, this is actually a blind alley. While I cannot assert it with utmost confidence, observe that
One can't attempt substitutions with rational numbers in general: try $f(\frac 12)$, for example, and see where you get. Or maybe try setting $y = \frac 1x$, and you'll have a problem. It seems that you can't get out of the integers.
Even if you get to the rational numbers, there's no way to place any monotonicity conditions on $f$, or continuity conditions, which allow you to move from the rational numbers to the real numbers.
Therefore, attempts to make this work are likely to be futile, though I absolutely, absolutely invite everyone to try.
EDIT: It seems that Sil, in the comments, has found a way out using a clever substitution! I'll have a go and change everything I've said above if it works out.
Hitting a right approach
A right approach , in this case, comes from observing that if we expect, in light of everything we've said before, that $f(x)=x$ for $x \in \mathbb R$, then the function $g(x) = f(x)-x$ should be identically zero for $x \in \mathbb R$.
However, when you form the functional equation for $g(x)$, something that one may not observe in the original equation manifests itself, and we see niceties occur.
Let's do it : let $g(x) = f(x)-x$ for $x \in \mathbb R$. Start with $P(x,y)$ and begin to substitute $g$ for $f$ everywhere, knowing that $f(z) = g(z)+z$.
\begin{align}
f(xf(y)+f(x)+y) &= xy+f(x)+f(y) \\
\implies g(xf(y)+f(x)+y) + xf(y)+f(x)+y &= xy+g(x)+x+g(y)+y \\
\implies g(xg(y)+xy+g(x)+x+y)+xg(y)+xy+g(x)+x+y &= xy+g(x)+x+g(y)+y \\
\implies g(xg(y)+xy+g(x)+x+y)+xg(y) &= g(y) \\
\implies g(xg(y)+xy+g(x)+x+y) &= (1-x)g(y)
\end{align}
We make a crucial observation about $g$ now. Let's call the last identity in the chain above $P^*(x,y)$.
Claim : if there is a $y_0$ such that $g(y_0) \neq 0$, then $g$ is surjective.
Proof : If $g(y_0) \neq 0$, then $g(y_0)$ belongs in the range of $g$ obviously. Let $T \in \mathbb R$ be arbitrary and let $x_0 = 1-\frac{T}{g(y_0)}$. Then, $(1-x_0)g(y_0)=T$. By observing $P^*(x_0,y_0)$, one sees that $T$ is in the range of $g$, as desired.
We will show that $g$ cannot be surjective in the next section, creating the contradiction.
An ideal choice, and a contradiction to $g$ being surjective
From surjectivity, we will now pick the most suitable candidate for the next substitution, leading to a vast simplification.
The key idea is the following : look at the left hand side of $P^*$, and you have $g(xg(y)+g(x)+xy+x+y)$. We want a value of $x$ or $y$ that can perhaps cancel a great number of terms. One such substitution is easily observed: $x=-1$. Recall that $f(-1)=-1$ (Yes, this was important!) so $g(-1) = 0$.
Now $P^*(-1,y)$ gives $$
g(-g(y)+0-y-1+y) = 2g(y) \implies g(-g(y)-1) = 2g(y)
$$
However, we mentioned that $g$ was surjective! Therefore, the above equation actually holds for all $z$ by merely substituting $y = g^{-1}(z)$ above where $y$ is any preimage of $z$.
Thus, we obtain $$
g(-z-1) = 2z
$$
for all $z$. This is a contradiction : take $z=-1$, then $g(0)=-2$ is obtained, which is a contradiction because $f(0)=0$ therefore $g(0)=0$.
Finally, we obtain that $g \equiv 0$, and $f(x) =x$ everywhere.