It is possible to represent (and "see")  only surfaces. You cannot see in the interior of a 3-dimensional set in R^3 (between any interior point and your eye there exist an infinity of other interior points obfuscating the first one).

Today you have problems easier by hand than with Maple :-)

f = (x+y)/(x+y+z), g = (x+z)/(x+y+z), h = (y+z)/(x+y+z)  have the same integral (by symmetry).

==> int(f) = 1/3*int(f+g+h) = 1/3*int(2) = 1/3*2*(4*Pi/3)/8 = Pi/9.

Edit.
A real challenge for any CAS would be to compute the integral of
(x^2017+y^2017)/(x^2017+y^2017+z^2017)
in the first octant of the unit ball. [=Pi/9]

dxdy  :=   y=-sqrt(1-x^2) .. sqrt(1-x^2), x=-1..1 ;

int(sin(x^2) *cos(x^2) , dxdy) =

1/2*int( sin(x^2+y^2) + sin(x^2 - y^2), dxdy ) =

1/2*int( sin(x^2+y^2) , dxdy)  + 0 =

1/2* int( sin(r^2)*r, r=0..1, t=0..2*Pi ) =

1/2*(1 - cos(1))*Pi =

Pi * sin(1/2)^2

It's a very difficult task to compute u3 even at a single point.
For example, u3(0) can be simplified to the triple integral of
f := -(2.533317577*10^7*I)*z*exp((-.64-12.56637062*I)*x^2-(7.853982892*10^7*I)*z^2-(25.13274124*I)*y^2)*y*x*BesselJ(0, 25.13274123*y*x)*BesselJ(0, 62831.85308*z*y)

But f is very oscillatory and to compute its triple integral over [0,1]^3 is almost impossible even using MonteCarlo method.
Probably some special methods will be needed.
Just look at this plot:
plot( eval(Re(f),[x=1/2,y=1/2]),z=0..1 );

P.S. Was this integral invented, or it appears in a concrete problem?

eq:=proc(L::listlist,x::name,y::name,z::name)
  primpart(LinearAlgebra:-Determinant(<Matrix([[x,y,z],L[]])|<1,1,1,1>>))=0
end:

map(eq,L,x,y,z);

TriOnUnitSphere:=proc(A::Vector[column](3),B::Vector[column](3),C::Vector[column](3)) #A,B,C = directions
uses LinearAlgebra,plots,plottools;
local a:=Normalize(A,2),b:=Normalize(B,2),c:=Normalize(C,2),u,v;
display(sphere(),  plot3d(1.01*Normalize(c+v*(a+u*(b-a)-c),2), u=0..1,v=0..1,color=red,style=surface,_rest),
pointplot3d([1.01*a,1.01*b,1.01*c],symbolsize=20,color=blue,symbol=solidsphere)  )
end:

TriOnUnitSphere(<1,2,3>, <-1,1,-2>, <-1,-1,3>);

TriOnUnitSphere(<1,0,0>, <-10,1,0>, <0,1,5>, grid=[200,200]);

with(LinearAlgebra):with(plots):with(plottools):
A1:=<0,0,1>: A2:=<0,3/5,4/5>: A3:=<4/5,3/5,0>: A4:=<5/13,0,12/13>:
seg:= (a,b)->spacecurve(Normalize(a+t*(b-a),2),t=0..1, thickness=3, color=red):
display(sphere([0,0,0],1),seg(A1,A2),seg(A2,A3),seg(A3,A4),seg(A1,A4));

You can use e.g.

binaryvariables = indets(L, specindex(integer,x));

where L is a list containing the constraints (and the function).
Or, you may generate their sequence using seq.
 

f:=(x,y) ->exp(-x^2-y^2):
L:=4: n:=64:
r1:=rand(0..1.0): r:=rand(-1..1.0):
a:=Vector(n,L*r): b:=Vector(n,L*r):
c:=Vector(n,0.3*r):
da:=Vector(n, i->0.5*''r1()''):
P:=proc(aa)
  global n,a,b,c,da,h,L;
  a:=a+eval~(da);
  map[inplace](frem,a,L);
  h:=add(c[i]*f(x+a[i], y+b[i]),i=1..n):
  plot3d(h, x=-L..L,y=-L..L, color="Aquamarine", scaling=constrained, view=-2..2, axes=none, style=surface, lightmodel=light1);
end:
plots:-display(seq(P(aa),aa=1..50),insequence=true);
## or ...
Explore(P(aa), aa=0..1, animate,loop,size=[800,800]);

Download req-int.mw

restart;
f:=r^2 *cos(theta)+r*sin(theta):
r1:= 0.3+0.1*cos(theta):
r2:= 0.5+0.1*cos(theta):
plots:-densityplot(f, theta=0..2*Pi, r=r1 .. r2,  
coords=polar, colorstyle = HUE, style = patchnogrid);

If you need labels:

plots:-display(%, labels=["x","y"]);

c:=Matrix([    # result = 22.613, [4, 2, 5, 3, 6, 1, 4]
[0., 3.60555127546399, 2.00000000000000, 8.06225774829855, 5.09901951359278, 1.41421356237310],
[3.60555127546399, 0., 2.23606797749979, 6.32455532033676, 2.23606797749979, 3.60555127546399],
[2.00000000000000, 2.23606797749979, 0., 8.06225774829855, 3.16227766016838, 1.41421356237310],
[8.06225774829855, 6.32455532033676, 8.06225774829855, 0., 8.06225774829855, 9.00000000000000],
[5.09901951359278, 2.23606797749979, 3.16227766016838, 8.06225774829855, 0., 4.47213595499958],
[1.41421356237310, 3.60555127546399, 1.41421356237310, 9.00000000000000, 4.47213595499958, 0.]]):

n := sqrt(numelems(c)):  N:={seq(1..n)}:
z:=add(add( c[i,j]*x[i,j],j=N minus {i}),i=N): 
conx:=seq( add(x[i,j], i=N minus {j})=1,j=N),  seq( add(x[i,j], j=N minus {i})=1,i=N):
conu:= seq(seq(u[i]-u[j]+n*x[i,j] <= n-1, i=N minus {1,j}), j=N minus {1}):

r := Optimization[LPSolve](z, {conx, conu}, #integervariables={seq(u[i],i=N minus {1})},
                           binaryvariables={ seq(seq(x[i,j],i=N minus {j}),j=N)} );

X:=eval(Matrix(n, symbol=x),{r[2][], seq(x[i,i]=0,i=1..n)}):
ord:=1: k:=1: to n do k:=max[index](X[k]); ord:=ord,k od:'order'=ord;

r := [22.6135858310497, [u[2] = -.999999998967403, u[3] = -2.99999997900383, u[4] = -8.88178419700125*10^(-16), u[5] = -2.00000000722817, u[6] = -3.99999998106902, x[1, 2] = 0, x[1, 3] = 0, x[1, 4] = 0, x[1, 5] = 0, x[1, 6] = 1, x[2, 1] = 0, x[2, 3] = 0, x[2, 4] = 1, x[2, 5] = 0, x[2, 6] = 0, x[3, 1] = 0, x[3, 2] = 0, x[3, 4] = 0, x[3, 5] = 1, x[3, 6] = 0, x[4, 1] = 1, x[4, 2] = 0, x[4, 3] = 0, x[4, 5] = 0, x[4, 6] = 0, x[5, 1] = 0, x[5, 2] = 1, x[5, 3] = 0, x[5, 4] = 0, x[5, 6] = 0, x[6, 1] = 0, x[6, 2] = 0, x[6, 3] = 1, x[6, 4] = 0, x[6, 5] = 0]]

order = (1, 6, 3, 5, 2, 4, 1)