(* We find the maximum surface area of a brick in the first octant with diagonal corners at (0,0,0) and on the plane 2x+4y+z=4 *)

z=4-2 x-4y
S=Expand[2(x y+x z+y z)]  (* surface area *)
Plot3D[S,{x,0,2},{y,0,1},RegionFunction -> Function[{x,y},4-2 x-4y≥0]]
ContourPlot[S,{x,0,2},{y,0,1},Contours->30,RegionFunction -> Function[{x,y},4-2 x-4y≥0]]

G=Grad[S,{x,y}] (* Find the critical points *)
crit=Solve[G=={0,0},{x,y}]    (* No critical point in the allowed region *)
S/.crit[[1]]

(* Need to consider the three boundary curves *)
f=S/.{x->0,y->t}
Plot[f,{t,0,1}]
cand1=Solve[D[f,t]==0]
f/.cand1[[1]]

g=S/.{x->t,y->0}
Plot[g,{t,0,2}]
cand2=Solve[D[g,t]==0]
g/.cand2[[1]]

yy=Solve[z==0,y]   (* express y in terms of x assuming z=0 *)
h=Expand[S/.yy[[1]]/.{x->t}]
Plot[h,{t,0,2}]
cand3=Solve[D[g,t]==0]
h/.cand3[[1]]

(* The global maximum is 4 when x=1, y=0, z=2 *)
S/.{x->1,y->0}
z/.{x->1,y->0}