In topology optimization, the treatment of stress constraints for very large scale problems (more than 100 million elements and more than 600 million stress constraints) has so far not been tractable due to the failure of robust agglomeration methods, that is, their inability to accurately handle the locality of the stress constraints. This article presents a three-dimensional design methodology that alleviates this shortcoming using both deterministic and robust problem formulations. The robust formulation, based on the three-field density projection approach, is extended and proved necessary to handle manufacturing uncertainty in three-dimensional stress-constrained problems. Several numerical examples are solved and further postprocessed with body-fitted meshes using commercial software. The numerical investigations demonstrate that: (1) the employed solution approach based on the augmented Lagrangian method is able to handle very large problems, with hundreds of millions of stress constraints; (2) three-dimensional stress-based results are extremely sensitive to slight manufacturing variations; (3) if appropriate interpolation parameters are adopted, voxel-based (fixed grid) models can be used to compute von Mises stresses with excellent accuracy; and (4) in order to ensure manufacturing tolerance in three-dimensional stress-constrained topology optimization, a combination of double filtering and more than three density field realizations may be required.