---

with(plots):
with(plottools):

f1 := torus([0, 0, 0], 1, 3):

f2 := translate(homothety(rotate(f1, 0, Pi/2, 0), 1/2), 0, 3, 0):
f3 := translate(homothety(rotate(f1, Pi/2, 0, 0), 1/2), 3, 0, 0):
f4 := translate(homothety(rotate(f1, 0, Pi/2, 0), 1/2), 0, -3, 0):
f5 := translate(homothety(rotate(f1, Pi/2, 0, 0), 1/2), -3, 0, 0):

f6 := display(f1, f2, f3, f4, f5, scaling = constrained):

f7 := translate(homothety(rotate(f6, 0, Pi/2, 0), 1/2), 0, 3, 0):
f8 := translate(homothety(rotate(f6, Pi/2, 0, 0), 1/2), 3, 0, 0):
f9 := translate(homothety(rotate(f6, 0, Pi/2, 0), 1/2), 0, -3, 0):
f10 := translate(homothety(rotate(f6, Pi/2, 0, 0), 1/2), -3, 0, 0):

f11 := display(f6, f7, f8, f9, f10, scaling = constrained):

f12 := translate(homothety(rotate(f11, 0, Pi/2, 0), 1/2), 0, 3, 0):
f13 := translate(homothety(rotate(f11, Pi/2, 0, 0), 1/2), 3, 0, 0):
f14 := translate(homothety(rotate(f11, 0, Pi/2, 0), 1/2), 0, -3, 0):
f15 := translate(homothety(rotate(f11, Pi/2, 0, 0), 1/2), -3, 0, 0):

display(f11, f12, f13, f14, f15, style = surface, scaling = constrained):