So-called area coordinates $(x,y,z)$ in the (rectangular & isosceles) triangle are defined as follows: $$ x/3 = \frac{\Delta APC}{\Delta ABC} \quad ; \quad y/3 = \frac{\Delta BPA}{\Delta ABC} \quad ; \quad z/3 = \frac{\Delta CPB}{\Delta ABC} $$ The $x$ and $y$ area-coordinates conform to the common Cartesian coordinates in this case, but the $z$ coordinate is different.
So far so good about the additional condition.
Solve $\,z\,$ therefrom and substitute the result $\,z=3-x-y\,$ into the inequality, giving:
$$
f(x,y) = \frac{xy}{5y^3+4}+\frac{yz}{5z^3+4}+\frac{zx}{5x^3+4} \leqslant \frac13
$$
Not a formal proof is attempted again
by plotting a contour map of the function, as depicted. Levels (`nivo`) of these isolines are defined
(in Delphi Pascal) as:
nivo := min + sqrt(g/grens)*(max-min); { grens = 10 ; g = 1..grens }
The blackness of the isolines is proportional to the (positive) function values; they are almost white
near the minimum and almost black near the maximum values.
Maximum and minimum values of the function are observed to be:
4.87041188016219E-0003 < f < 3.33333101070422E-0001And the maximum is observed to be near the triangle's center of gravity. The $\color{blue}{\mbox{blue}}$ spot is where $\,\left| f(x,y) - 1/3 \right| < 0.001$ .