How to plot Clebsch surface with those 27 lines?

Question

I can get its Algebraic Equation by this code:

Entity["Surface", "ClebschDiagonalCubic"][EntityProperty["Surface", "AlgebraicEquation"]]

Then get its graphics:

ContourPlot3D[81 (x^3 + y^3 + z^3) - 9 (x^2 + y^2 + z^2) - 
   189 (x^2 y + x^2 z + x y^2 + x z^2 + y^2 z + y z^2) + 54 x y z - 
   9 (x + y + z) + 126 (x y + x z + y z) - 1 == 0, {x, -1, 1}, {y, -1,
   1}, {z, -1, 1}, Boxed -> False, BoxRatios -> {1, 1, 1}, 
 AxesOrigin -> {0, 0, 0}, AxesStyle -> {Red, Green, Blue}, Mesh -> None]

As we know, there are 27 lines in its surface:

But I don't know how to draw these lines onto the surface

Answer 1

EDIT:

With help from @cvgmt, both fixing a wrong term in the equation(!) and adding constraints for finding lines:

With[{expr =
   1 + 54 x y z - 9 (x + y + z) + 126 (x y + x z + y z) - 
     9 (x^2 + y^2 + z^2) -
     189 (x^2 y + x y^2 + x^2 z + y^2 z + x z^2 + y z^2) + 
     81 (x^3 + y^3 + z^3) == 0},
 (* Expression over a parametric infinite line. *)
 expr /. Thread[{x, y, z} -> {a, b, c} + t {u, v, w}] //
    (* Find lines which lie on the surface. *)
    Solve[
      (* Expression must be true for all t,
      that is, it must lie on the surface everywhere. *)
      Resolve[ForAll[t, #]] &&
       (* Constrain solutions:
       line offset and line direction vectors must be at right angle,
       and line direction vector must be unit length. *)
       {a, b, c} . {u, v, w} == 0 && u^2 + v^2 + w^2 == 1 &&
       (* Prevent finding lines with "direction" negated. *)
       ((w > 0) || (w == 0 && v > 0) || (w == 0 && v == 0 && u > 0)),
      {a, b, c, u, v, w}, Reals] & //
   (* Line segments inside the plotted area. *)
   RegionIntersection[
      InfiniteLine[{a, b, c}, {u, v, w}], Ball[]] /. N[#] & //
  (* Create animated graphics. *)
  Animate[
    Show[
     (* Plot the surface. *)
     ContourPlot3D[expr, Element[{x, y, z}, Ball[]],
      PlotPoints -> 100, ContourStyle -> Opacity[9/10], Boxed -> False,
      BoxRatios -> Automatic, Axes -> None, Mesh -> None,
      RegionBoundaryStyle -> None,
      SphericalRegion -> True, ViewAngle -> 10 Degree,
      ViewPoint -> Dynamic[{5 Sin[a], 5 Cos[a], 5/2}]],
     (* Plot line segments. *)
     Graphics3D[{Thick, CapForm["Butt"], #}]],
    {a, 0, 2 Pi}] &]

---

Original answer:

This is more of an extended comment - but anyway, it would seem that Mathematica has at least some success finding those lines on the surface:

With[{eqn = -1 + 54 x y z - 9 (x + y + z) + 126 (x y + x z + y z) - 
    9 (x^2 + y^2 + z^2) - 
    189 (x^2 y + x y^2 + x^2 z + y^2 z + x z^2 + y z^2) + 
    81 (x^3 + y^3 + z^3)},
 CoefficientList[eqn /. Thread[{x, y, z} -> {a, b, c} + t {u, v, w}], 
        t] & // Map[# == 0 &] //
     Solve[
       And @@ # && {a, b, c} . {u, v, w} == 0 && 
        u^2 + v^2 + w^2 == 1, {a, b, c, u, v, w}, Reals] & // N //
   RegionIntersection[InfiniteLine[{a, b, c}, {u, v, w}], 
      Cuboid[{-1, -1, -1}, {1, 1, 1}]] /. # & //
  ListAnimate@
    Table[
     Show[ContourPlot3D[eqn == 0, {x, -1, 1}, {y, -1, 1}, {z, -1, 1},
       PlotPoints -> 50, ContourStyle -> Opacity[9/10], 
       Boxed -> False, BoxRatios -> Automatic, Axes -> None, 
       Mesh -> None, SphericalRegion -> True, 
       ViewPoint -> {5 Sin[a], 5 Cos[a], 5/2}],
      Graphics3D[{Thick, #}]],
     {a, 0, 2 Pi - Pi/64, Pi/64}] &]

Answer 2

Provide a more concise solution:

expr = 1 + 54 x y z - 9 (x + y + z) + 126 (x y + y z + z x) - 
   9 (x^2 + y^2 + z^2) - 
   189 (x^2 y + x y^2 + y^2 z + y z^2 + z^2 x + z x^2) + 
   81 (x^3 + y^3 + z^3);
sols = Solve[{Resolve[
     ForAll[t, (expr /. 
         Thread[{x, y, z} -> {a, b, c} + t {u, v, w}]) == 0]], {a, b, 
       c} . {u, v, w} == 0, u^2 + v^2 + w^2 == 1, 
    w > 0 || w == 0  && v > 0 || w == 0 && v == 0 && u > 0}, {a, b, c,
     u, v, w}, Reals];
Show[ContourPlot3D[expr == 0, {x, -2, 2}, {y, -2, 2}, {z, -2, 2}, 
  Mesh -> None, PlotPoints -> 60, MaxRecursion -> 4, 
  ContourStyle -> Cyan], 
 Graphics3D[{White, InfiniteLine[{a, b, c}, {u, v, w}] /. sols}], 
 Boxed -> False, Axes -> False, ViewPoint -> Front]

Explanation

1. We can get the lines in the surface by the parametric equation of a line to combine the original equation to get eqns.
2. Since the original equation is $0$, We can let the equation of the substituted line be constant $0$
3. I assumed that the line was translated through the origin to the surface:

So it has to satisfy:{a,b,c}.{u,v,w}==0

4. Since {u,v,w} is a vector, we can assume that it's distributed along the circumference of the unit circle. So it meet:u^2+v^2+w^2==1

5. Now, one last problem, we solve the equation this way and we get $54$ solutions. Because MMA can't tell the difference between $-vec$ and $vec$ being the same vector. We have to limit it with some conditions:

    ((w > 0)(*At an acute Angle to the w axis*)
      || (w == 0 && v > 0)(*Perpendicular to the w axis but at an acute Angle to the v axis*)
      || (w == 0 && v == 0 && u > 0)(*Perpendicular to both the w and v axes,but at an acute Angle to the u axis*))
   

---

But I have an additional puzzle: Why can I only solve three lines on a cubic surface $x^3 + 3 y^3 + z^3-2 x^2 + 5 x y - x + 7=0$? Why can't I get $27$ too as this statement?