It can in principle be done using RegionPlot3D
. But to get sharp edges, one needs to crank up the number of PlotPoints
. Here is an example with four cubes that doesn't take too long to plot:
n = 4;
insideCube[pt_, l_] := And @@ Thread[Abs[pt] < l]
rotations =
MapThread[
RotationMatrix[{{0, 0, 1}, {Cos[#1] Sin[#2], Sin[#1] Sin[#2],
Cos[#2]}}] &, RandomReal[{0, Pi}, {2, 4}]]
(*
==> {{{0.701529, 0.148251, -0.69705}, {0.148251, 0.926364,
0.346226}, {0.69705, -0.346226, 0.627892}}, {{0.989833,
0.065586, -0.12621}, {0.065586, 0.576913,
0.814168}, {0.12621, -0.814168, 0.566746}}, {{0.999495, -0.0316019,
0.00319107}, {-0.0316019, -0.979314,
0.199866}, {-0.00319107, -0.199866, -0.979818}}, {{0.18502,
0.225404, -0.956536}, {0.225404, 0.937659,
0.264555}, {0.956536, -0.264555, 0.122679}}}
*)
rotatedCubes =
Or @@ MapThread[
insideCube[#1.{x, y, z}, #2] &, {rotations,
RandomReal[{1, 1.1}, n]}]
(*
==> (Abs[0.701529 x + 0.148251 y - 0.69705 z] < 1.02814 &&
Abs[0.148251 x + 0.926364 y + 0.346226 z] < 1.02814 &&
Abs[0.69705 x - 0.346226 y + 0.627892 z] <
1.02814) || (Abs[0.989833 x + 0.065586 y - 0.12621 z] < 1.04828 &&
Abs[0.065586 x + 0.576913 y + 0.814168 z] < 1.04828 &&
Abs[0.12621 x - 0.814168 y + 0.566746 z] <
1.04828) || (Abs[0.999495 x - 0.0316019 y + 0.00319107 z] <
1.0917 && Abs[-0.0316019 x - 0.979314 y + 0.199866 z] < 1.0917 &&
Abs[-0.00319107 x - 0.199866 y - 0.979818 z] <
1.0917) || (Abs[0.18502 x + 0.225404 y - 0.956536 z] < 1.02296 &&
Abs[0.225404 x + 0.937659 y + 0.264555 z] < 1.02296 &&
Abs[0.956536 x - 0.264555 y + 0.122679 z] < 1.02296)
*)
RegionPlot3D[rotatedCubes, {x, -2, 2}, {y, -2, 2}, {z, -2,
2}, PlotPoints -> 120, Mesh -> False]

To cut through the shape and see the inside, you'd have to add another condition to RegionPlot3D
or you can restrict the PlotRange
as in this example:
Show[%, PlotRange -> {{-1.5, 1.5}, {-1.5, 1.5}, {-1.5, 0}}]
