The code below should produce a red cube behind a green cube, in turn behind a blue cube. I intend the matrices tmat and pmat to produce a perspective projection from a viewpoint slightly to the left of the midpoint of the left-hand edge of the frame, and I think 8 times as far away from the screen as the frame is wide.

I have deliberately listed the cubes not in the order of their respective depths, in order to demonstrate that they do not seem to occlude each other correctly, unless I first sort them according to their depths.

Questions:

1) Why is the occlusion incorrect, and how can I fix this ? I seem to have to rely on the "Painter's Algorithm", i.e. drawing nearer objects AFTER further objects, which is well known not to give correct results in general, since objects need not be totally ordered by depth. Other experiments I have done with more objects lead me to believe that it is not a simple problem with surface normals caused by some confusion between left- and right-handed coordinate systems, but perhaps I am wrong about that.

2) Is there a way to draw the edges of the cubes, without also drawing the edges that should be hidden ? I have tried removing the command EdgeForm[None], and using all 27 combinations of $a,b \in \{0, 0.5, 1\}$ with FaceForm[{cubedatum[[2]],Opacity[a]}, {cubedatum[[2]], Opacity[b]}] instead of FaceForm[cubedatum[[2]],None], with Lighting->{{"Directional",White,{{0,0,1},{0,0,0}}}}, Lighting->"Ambient", and Lighting->"Neutral", but nothing worked.

I apologise for asking two questions at a time, but I suspect that their answers will be closely related.

tmat = {{1, 0, 0, 0}, {0, 1, 0, 1/2}, {0, 0, 1, 160/33}, {0, 0, 0, 1}};
pmat = {{8, 0, 13/40, 0}, {0, 8, -1/2, 0}, {0, 0, 0, 0}, {0, 0, -1, 
    8}};
cubeside = 1/8;
gap = 1/20;
cubecentrexs = {0, -gap, gap};
cubecentreys = 2 {0, -gap, gap};
cubecentrezs = {0, -gap, gap};
cubecentres = Transpose[{cubecentrexs, cubecentreys, cubecentrezs}];
cubecolours = {Green, Red, Blue};
(* Don't sort the cubes: *)    
cubedata = Transpose[{cubecentres, cubecolours}];
cubes = Table[{FaceForm[cubedatum[[2]], None], EdgeForm[None], 
    Cuboid[cubedatum[[1]] - (cubeside/2) {1, 1, 1}, 
     cubedatum[[1]] + (cubeside/2) {1, 1, 1}]}, {cubedatum, cubedata}];
picwithoutsorting = 
 Graphics3D[cubes, ViewMatrix -> {tmat, pmat}, Boxed -> False, 
  Lighting -> {{"Directional", White, {{0, 0, 1}, {0, 0, 0}}}}]
(* Sort the cubes by their z-coordinate: *)    
cubedata = 
  Sort[Transpose[{cubecentres, 
     cubecolours}], #1[[1, 3]] <= #2[[1, 3]] &];
cubes = Table[{FaceForm[cubedatum[[2]], None], EdgeForm[None], 
    Cuboid[cubedatum[[1]] - (cubeside/2) {1, 1, 1}, 
     cubedatum[[1]] + (cubeside/2) {1, 1, 1}]}, {cubedatum, cubedata}];
picwithsorting = 
 Graphics3D[cubes, ViewMatrix -> {tmat, pmat}, Boxed -> False, 
  Lighting -> {{"Directional", White, {{0, 0, 1}, {0, 0, 0}}}}]

Output:

This is picwithoutsorting:

And this is picwithsorting:

I have found that using ViewVector, ViewCenter and ViewVertical appropriately instead of ViewMatrix gives very nearly the same geometry but with correct occlusions. However, one disadvantage of this workaround is that coordinates seem to be scaled to fit the bounding box of the scene. Very slight differences in size result which is fine for static images, but when generating frames for an animation, can result in a "wobbling" effect. A possible workaround for which might be to use SphericalRegion->True, but this also potentially introduces a lot of unwanted empty space in the scene.