The idea is similar to Kitonum's but has some advantages.
Note that some details will be omitted.

As it was observed,

f:=a^2+b^2+c^2+a*b+a*c+b*c-(a+b-c)*sqrt(2*a*b+a*c+b*c)-(a+c-b)*sqrt(a*b+2*a*c+b*c)-(b+c-a)*sqrt(a*b+a*c+2*b*c);
is positively homogeneous and symmetric, so we may suppose:
-1 <= a, b, c <= 1.

It is easy to see that c=0 implies a=b=0.
It will be enough to consoder c=1 and c=-1.

1.     c=1

f1:= eval(f, c=1);
plots:-implicitplot(f1, a=-1..1,b=-1..1, grid=[1000,1000]);

So, we have two branches. But they are symmetric because f(a,b,1)=f(b,a,1).
So we may consider only a>0, b<0.
Actually a will be in a1..a2  where a1, a2 could be computed [here at least of one sqrt must be 0].
So, the corresponding b will be h(a) = RootOf(f1, b=-1/5 .. -1/3).
[Note that h(a) can be expressed as a RootOf of a polynomial]

2.    c=-1
There are no solutions.

Conclusion:
The roots are:  a=k*u, b=k*h(u), c=k,  for all u in [a1,a2] and k>=0
and all permutations of these.

Edit.
 

a1:=-1+sqrt(1+RootOf(_Z^8-8*_Z^7-22*_Z^6+228*_Z^5+17*_Z^4-624*_Z^3-352*_Z^2+256*_Z+256, index = 2));

a2:=1;

h:= a -> 
RootOf(
a^16-8*a^15*b-4*a^14*b^2-48*a^13*b^3+154*a^12*b^4+40*a^11*b^5-304*a^10*b^6-56*a^9*b^7+451*a^8*b^8-56*a^7*b^9-304*a^6*b^10+40*a^5*b^11+154*a^4*b^12-48*a^3*b^13-4*a^2*b^14-8*a*b^15+b^16-8*a^15-8*a^14*b-192*a^13*b^2-56*a^12*b^3+96*a^11*b^4+456*a^10*b^5+296*a^9*b^6-560*a^8*b^7-560*a^7*b^8+296*a^6*b^9+456*a^5*b^10+96*a^4*b^11-56*a^3*b^12-192*a^2*b^13-8*a*b^14-8*b^15-4*a^14-192*a^13*b-480*a^12*b^2-1104*a^11*b^3+1564*a^10*b^4+3040*a^9*b^5-404*a^8*b^6-4048*a^7*b^7-404*a^6*b^8+3040*a^5*b^9+1564*a^4*b^10-1104*a^3*b^11-480*a^2*b^12-192*a*b^13-4*b^14-48*a^13-56*a^12*b-1104*a^11*b^2+720*a^10*b^3+4616*a^9*b^4+3464*a^8*b^5-4864*a^7*b^6-4864*a^6*b^7+3464*a^5*b^8+4616*a^4*b^9+720*a^3*b^10-1104*a^2*b^11-56*a*b^12-48*b^13+154*a^12+96*a^11*b+1564*a^10*b^2+4616*a^9*b^3+6008*a^8*b^4-1792*a^7*b^5-7416*a^6*b^6-1792*a^5*b^7+6008*a^4*b^8+4616*a^3*b^9+1564*a^2*b^10+96*a*b^11+154*b^12+40*a^11+456*a^10*b+3040*a^9*b^2+3464*a^8*b^3-1792*a^7*b^4-3896*a^6*b^5-3896*a^5*b^6-1792*a^4*b^7+3464*a^3*b^8+3040*a^2*b^9+456*a*b^10+40*b^11-304*a^10+296*a^9*b-404*a^8*b^2-4864*a^7*b^3-7416*a^6*b^4-3896*a^5*b^5-7416*a^4*b^6-4864*a^3*b^7-404*a^2*b^8+296*a*b^9-304*b^10-56*a^9-560*a^8*b-4048*a^7*b^2-4864*a^6*b^3-1792*a^5*b^4-1792*a^4*b^5-4864*a^3*b^6-4048*a^2*b^7-560*a*b^8-56*b^9+451*a^8-560*a^7*b-404*a^6*b^2+3464*a^5*b^3+6008*a^4*b^4+3464*a^3*b^5-404*a^2*b^6-560*a*b^7+451*b^8-56*a^7+296*a^6*b+3040*a^5*b^2+4616*a^4*b^3+4616*a^3*b^4+3040*a^2*b^5+296*a*b^6-56*b^7-304*a^6+456*a^5*b+1564*a^4*b^2+720*a^3*b^3+1564*a^2*b^4+456*a*b^5-304*b^6+40*a^5+96*a^4*b-1104*a^3*b^2-1104*a^2*b^3+96*a*b^4+40*b^5+154*a^4-56*a^3*b-480*a^2*b^2-56*a*b^3+154*b^4-48*a^3-192*a^2*b-192*a*b^2-48*b^3-4*a^2-8*a*b-4*b^2-8*a-8*b+1,
b, -3/10 .. -27/100)   :

To plot the solution use:

plot3d([ [k*u, k*h(u), k], [k*u, k, k*h(u)], [k*h(u), k*u,  k],
[k*h(u), k, k*u], [k, k*h(u), k*u], [k, k*u, k*h(u)] ],  u=a1..a2, k=0..10);

A strange fact (which I don't understand) is that in Maple there are both left and right versions of spherical coordinates!

The left one is for plot3d:
[u,v,w] -> [u*cos(v)*sin(w), u*sin(v)*sin(w), u*cos(w)]

and the right one is in VectorCalculus:
[u,v,w] -> [u*cos(w)*sin(v), u*sin(w)*sin(v), u*cos(v)].

Of course, as Carl has shown they can be redefined by the user, but it is confusing.

restart;
S:= rhs~(solve({-x-y+3*z =a, -x+2*y = b, 3*x-2*y-z = c}, [x,y,z])[]):
A:=eval(S, [a=1,b=0,c=0]):
B:=eval(S, [a=0,b=1,c=0]):
C:=eval(S, [a=0,b=0,c=1]):
OO:=[0,0,0]:
plots:-polygonplot3d([[OO,A,B], [OO,A,C], [OO,B,C]]);

In your example, you can do this only if the rank of the system is r=19 and the rank corresponding to the unknowns Z1,...,Z19 is also r=19.
In this case you can solve wrt Z1,...,Z19, e.g.
solve( sys, {seq(Z||i,i=1..r)} );

Edit:
I have interpreted your matrix as representing a piecewise linear (discontinuous) function (as your plot shows).

I think that this is intentional. After all, we can see 10 as a binary number.

b:=convert(0.01,binary);
    .1010001111 10^-6
op(b);
    1010001111, -16
``(op(1,b))*2^(``(op(2,b)));

Use one of the following:

plot3d(p1, -1 .. 1, -1 .. 1);

expr:=piecewise(x^2+y^2 <= 1, x*y-y^2, 0);
plot3d(expr, x = -1 .. 1, y = -1 .. 1);

plot3d('p1(x,y)', x = -1 .. 1, y = -1 .. 1);

The command to plug n=1 into the expression S is:
eval(S, n=1);

The command solve is to find the root(s) of the equation S=0, and the syntax solve(S, n=1) is wrong.
You should read (or at least browse) the Maple User Manual.

No need for assumptions, but instead:

- Don't use square brackets [ ... ],  they are list delimiters.
- Don't use inert Int , or use value(%) after that.
- (optional) use Maple (1D) input.

int( (x/L)^3 - 1, x=0..L );

           
                              

Yes, arrays of plots refuse to be exported (by context menu or by plottools:-exportplot).
Workaround:

restart;
p:=seq(plot(x^k,x=-1..1),k=1..4):
plotsetup(png,plotoutput="d:/temp/ex-array.png");
plots:-display(<p[1],p[2];p[3],p[4]>);
plotsetup(default);

subs(_Z9 = n, rhs(symbolic));
(of course the numeric suffix may differ).

 

When coords=cylindrical,  r  is ploted as a function of theta and z.
Even in your worksheet the notation is misleading; you use plot3d( f(r), theta=..., r=...)
but actually r should be z; it does not matter, of course, but your intentions are then not clear.