, 2019

For $x = y = z = 1$ we have $f(f(1)) = 1$, and hence for $x = y = z = f(1)$ we have
$$3f(f(1)) = 3f(1)^2 \Rightarrow f(1)^2 = 1 \Rightarrow f(1) = 1.$$
For $y = z = 1$, we have for every $x \in \mathbb{R}_+$:
$$\begin{aligned} f(f(x)) + f(x) &= 2x \\ f(f(x)) &= 2x - f(x) \end{aligned} \qquad (1)$$
For $z = 1$ from the given relation we get:
$$f(xf(y)) + f(y) + f(f(x)) = xy + y + x$$
And from (1) we find:
$$f(xf(y)) = xy + y - x + f(x) - f(y), \qquad (2)$$
for all $x, y \in \mathbb{R}_+$. Putting to (2) the $f(x)$ in the place of $x$ we have, from (1), for all $x, y \in \mathbb{R}_+$:
$$f(f(x)f(y)) = f(x)y + y + 2x - 2f(x) - f(y) \qquad (3)$$
and interchanging $x$ and $y$ we get the relation
$$f(f(y)f(x)) = f(y)x + x + 2y - 2f(y) - f(x) \qquad (4)$$
From (3) and (4) we have, for all $x, y \in \mathbb{R}_+$:
$$(f(x) - 1)(y - 1) = (f(y) - 1)(x - 1), \qquad (5)$$
From which for $y = 2$ we find
$$f(x) = cx - c + 1 = c(x - 1) + 1 \qquad (6)$$
where $c = f(2) - 1$. By substitution to the given equation we have
$$c^2(xy + yz + zx - x - y - z) + c(x + y + z) - 3c + 3 = xy + yz + zx,$$
from which for $x = y = z$ we find
$$3(c^2 - 1)x^2 + 3c(1-c)x + 3(1-c) = 0 \\ \Leftrightarrow c^2 - 1 = 0,\ c(1-c) = 0,\ 1-c = 0 \Leftrightarrow c = 1.$$
$c = 1$, and hence: $f(x) = x$, for every $x > 0$.