Using sage to plot 3d region Say I want to plot 
$$E_1 = \{(x,y,z)\in\mathbb{R}^3\mid 0\le z\le 1-y, \sqrt{x}\le y\le 1, 0\le x \le 1\}$$
or
$$E_2 = \{(x,y,z)\in\mathbb{R}^3\mid 0\le z\le \sqrt{9-y^2}, {3x}\le y\le 3, 0\le x \le 1\}. $$
Is it possible to accomplish this in Sage? For example, for $E_1$, I could do something like 
x,y,z=var('x,y,z');
F = parametric_plot3d((x,sqrt(x),y), (x,0,1), (y,0,1), color='red');
G = parametric_plot3d((x,y,1-y), (x,0,1), (y,0,1),opacity=.5);
XY = plot3d(0, (x,0,1), (y,0,1),opacity=.5);
XZ = parametric_plot3d((x,0,z), (x,0,1), (z,0,1),opacity=.5);
YZ = parametric_plot3d((0,y,z), (y,0,1), (z,0,1),opacity=.5);
show(F+G+XY+XZ+YZ);
But of course displays each of the surfaces with in the box $[0,1]\times[0,1]\times[0,1]$, not the exact region $E_1$.
 A: Plotting a 3d region defined by inequalities in SageMath
Sage has no built-in function for plotting 3d regions.
The best workaround is probably to parametrise the faces
and plot them with parametric_plot3d.
Start by defining a function to plot a parametric surface,
taking as arguments


*

*a triple xyz of coordinate functions of two variables

*a $u$-range uu

*a $v$-range vv

*an optional opacity argument


and call this function surf:
def surf(xyz, uu, vv, opacity=None):
    return parametric_plot3d(xyz, uu, vv, opacity=opacity)

Below we use an opacity (or "alpha") of 0.3 (ie transparency 0.7)
so that we can see well through the faces and better understand
the volume they bound.
For $E_1$, we parametrise with $u \in [0, 1]$ and $v \in [0, 1]$,
setting $v$ equal to $y$.
sage: uu, vv, op = (0, 1), (0, 1), 0.3
sage: E1_faces = (
....:     [lambda u, v: u*v^2, lambda u, v: v, lambda u, v: 0      ],
....:     [lambda u, v: v^2,   lambda u, v: v, lambda u, v: u*(1-v)],
....:     [lambda u, v: u*v^2, lambda u, v: v, lambda u, v: 1-v    ],
....:     [lambda u, v: 0,     lambda u, v: v, lambda u, v: u*(1-v)])
sage: E1 = add((surf(f, uu, vv, op) for f in E1_faces), Graphics())
sage: E1.show(viewer='threejs', aspect_ratio=1)
Launched html viewer for Graphics3d Object

For $E_2$, we parametrise with $u \in [0, 1]$ and $v \in [0, 3]$,
with $v$ equal to $y$ again.
sage: uu, vv, op = (0, 1), (0, 3), 0.3
sage: E2_faces = (
....:     [lambda u, v: u*v/3, lambda u, v: v, lambda u, v: 0            ],
....:     [lambda u, v: v/3,   lambda u, v: v, lambda u, v: u*sqrt(9-v^2)],
....:     [lambda u, v: u*v/3, lambda u, v: v, lambda u, v: sqrt(9-v^2)  ],
....:     [lambda u, v: 0,     lambda u, v: v, lambda u, v: u*sqrt(9-v^2)])
sage: E2 = add((surf(f, uu, vv, op) for f in E2_faces), Graphics())
sage: E2.show(viewer='threejs', aspect_ratio=1)
Launched html viewer for Graphics3d Object

Online demo: sagecell.
A: This is an old question, but if you do not want to parametrize your surfaces, you can try something like this:
x, y, z = var('x,y,z')

def plot_solid_from_inequalities(constraints, xrange, yrange, zrange,
                                 plot_points=100, epsilon=1e-5, **kwargs):
    """
    Plots a 3d solid given a list of inequalities (constraints).
    The `constraints` argument must contain a list of expressions in the
    variables `x`, `y` and `z`. Each expression will be the left hand side
    of an inequality in the form `g(x,y,z) >= 0`.
    
    If some of the boundary surfaces do not show up, or are incomplete,
    try to increase the value of `epsilon`.
    """

    def region(constraint):
        """
        Create a region bounded by all the constraints except the one
        that is currently being plotted.
        """
        return lambda x, y, z: all(c(x=x, y=y, z=z) >= 0
                                   for c in constraints
                                   if c != constraint)

    plots = [implicit_plot3d(f - epsilon, xrange, yrange, zrange,
                             plot_points=plot_points, region=region(f),
                             **kwargs)
             for f in constraints]

    return(sum(plots))

Then you could plot the two examples by doing this:
constraints = [z, 1 - y - z, y - sqrt(x), 1 - y, x, 1 - x]

plot_solid_from_inequalities(
   constraints, 
   (x,0,1), (y,0,1), (z,0,1), 
   color="yellow", alpha=0.7
)

and

constraints = [z, 1 - y - z, y - sqrt(x), 1 - y, x, 1 - x]

plot_solid_from_inequalities(
   constraints, 
   (x,0,1), (y,0,1), (z,0,1), 
   color="yellow", alpha=0.7
)

The idea is to use implicit plots of the boundary surfaces, and restrict them to the region defined by the inequalities.
One problem with this is that it is very slow.  When I try this on a sage cell server, I have to reduce plot_points to 50, otherwise it times out.  The resulting plot is then quite ugly.  It works nicely on CoCalc, though.
