Is this the Tucker-Gergonne-Nagel cubic ? or some other one

For a fixed non-equilateral triangle $ABC$, it appears to be the case that the locus of points $M$ for which the equation $$ xyz=rS $$ is satisfied is a simple closed curve in the interior of the given triangle.

As an example, for a triangle $ABC$ with $a=5,b=6,c=5$, the diagram above shows the triangle together with the locus of qualifying points $M$.

For the above example, using vertices $$ A=(-3,0),;;;B=(0,4),;;;C=(3,0) $$ the locus of qualifying points $M$ is the set of points $(x,y)$ with $1le y le {large{frac{12}{7}}};$such that the equation $$ x^2 = frac { 9(y-1)(4-y)(12-7y) } { 16(7y-3) } $$ is satisfied.

This is not an answer, but it's too long for comment.

Denote $c=AB$, $a=BC$, $c=AB$ and $p=frac{a+b+c}{2}=frac{P_{ABC}}{2}$. It's well-known that if $M$ is the Gergonne point of the triangle $ABC$, then $AB_1=AC_1=p-a$, $BC_1=BA_1=p-b$, $CA_1=CB_1=p-c$ and $S=frac{(p-a)(p-b)(p-c)}{r}$.

Now, suppose that point $M$ inside $triangle ABC$ satisfy condition $S=frac{xyz}{r}$. Then, we have the following system $$ begin{cases} xyz=(p-a)(p-b)(p-c), \ xyz=(b-x)(c-y)(a-z), end{cases} $$ where $(x,y,z)in(0,b)times(0,c)times(0,a)$. (The second equation is the Ceva's theorem)

Moreover, it's easy to see that if $(x,y,z)$ satisfies this system, then corresponding cevians will be concurrent and $S=frac{xyz}{r}$. Thus, we need to determine all solutions of the system above. We know already two of them: $(x,y,z)=(p-a,p-b,p-c)$ and $(x,y,z)=(p-c,p-a,p-b)$.

Probably, it would be better to rewrite this system as follows (here $u=p-a$, $v=p-b$, $w=p-c$; note that $a=(p-b)+(p-c)=v+w$ and similarly $b=w+u$ and $c=u+v$): $$ begin{cases} xyz=uvw, \ xyz=(v+w-x)(w+u-y)(u+v-z), end{cases} text{where}~(x,y,z)in (v+w,w+u,u+v). $$ Exapnding the second equation (along with the first equation) gives $$ begin{cases} xyz=uvw \ sumlimits_{cyc}(w+u)(u+v)x-sumlimits_{cyc}(v+w)yz=sumlimits_{cyc}u^2(v+w) end{cases} $$ Even for $u=v=w=1$ this system, I think, has infinitely many soltions $$ begin{cases} xyz=1, \ 2(x+y+z)-(xy+yz+zx)=3. end{cases} $$