## Randomly packing spheres of fixed radius within a cube

12

9

How can I have Mathematica randomly place spheres in a cube so they won't overlap? The cube is $20 \times 20 \times 20$, and the spheres have a radius of $0.7$.

Interesting problem. Could you explain the context of the problem? Just curiosity? – Dr. belisarius – 2012-10-12T01:37:16.090

How many spheres? – Mark McClure – 2012-10-12T01:38:48.137

I'm trying to solve a packing density problem with mathematica. Eventually, I want to fill the cube to a volume fraction of 1/3 of the spheres. Right now I'm struggling to plot the cube and hard spheres. – Jen – 2012-10-12T01:43:37.233

So you are just trying to plot right now? You can do cube = {Opacity[.3], Cuboid[{0, 0, 0}, {20, 20, 20}]}; spheres = Table[{Black, Sphere[RandomReal[{.35, 20 - .35}, 3], .7]}, {j, 1, 10}]; Graphics3D[{cube, spheres}] – David Slater – 2012-10-12T01:54:50.220

It's easier to pack them than to place them randomly, especially if you want to assume that they cannot sit still, suspended by the air. – DavidC – 2012-10-12T05:44:08.390

Have you seen this, this, and this?

– J. M.'s ennui – 2012-10-12T07:56:58.287

I meant to say "easier to pack them into a lattice structure"... – DavidC – 2012-10-12T09:00:51.903

well he didn't actually say random...simple cubic gets over 1/3 iirc – george2079 – 2012-10-12T11:58:02.697

Experimentally, I can get 287 spheres of radius 0.3125 inches into a cube of 3.937 inches cubed. So, all we do is scale stuff up, right??? This means, I think, that one can get 287 spheres of radius 0.7 inches into a cube of 8.81888 inches cubed. This means--and if I am missing something, I apologize--that one can get 3347 spheres of radius 0.7 into a cube of 20 x 20 x 20. One third of THAT number is: 1115 spheres of radius 0.7 into a 20x20x20 cube. – user62138 – 2018-12-31T19:03:36.917

10

Here's a brute force and soporifically slow method.

spheres = {};
Dynamic[Length[spheres]]


Then create a new random sphere, and if it's not too close to any existing sphere add it to the list:

While[Length[spheres] < 1856, s = RandomReal[{0.7, 20 - 0.7}, 3];
If[And @@ (Norm[# - s] > 1.4 & /@ spheres), AppendTo[spheres, s]]];


According to my calculations, you need 1856 spheres of radius 0.7 to fill 1/3 of the volume of a 20 x 20 x 20 cube.

To display:

cube = {Opacity[0.3], Cuboid[{0, 0, 0}, {20, 20, 20}]};
Graphics3D[{cube, Sphere[#, 0.7] & /@ spheres}, Boxed -> False] As I stated at the beginning, this is very slow and inefficient. The above image is from when I gave up after 1500 spheres. You can Alt. and resume at your leisure. It's entirely possible, what with the randomness, that you will never fit the 1856 spheres in.

1It's a good way to get the OP started. An extension would be to allow for wiggling the locations. One could also treat it as like an electromagnetics problem, for instance — think of the spheres as all having unit charge and reposition the spheres just as they would when another charged sphere is added to the existing bunch, subject to the constraint that none leave the box. – rm -rf – 2012-10-12T03:36:00.937

The alternate type approach is to randomly place spheres so that each exactly touches three others. Much tougher to implement but it will get you ultimately to higher density (50+%). – george2079 – 2012-10-12T04:01:00.810

1A relatively simple scheme could be to select a random sphere in the list, then creating a new sphere that touches it in a random direction. This got me to 1856 spheres in a few minutes. – wxffles – 2012-10-12T04:54:05.383

Why is the standing assumption that the metric space is Euclidean (L2-norm)? – alancalvitti – 2012-10-12T06:22:59.353

7

Here's a brute-force method (based on a bounded random walk) that attempts to generate n nonintersecting spheres of radius r within a cube of side length s, as well as reporting the actual number of spheres it managed to generate after lim unsuccessful iterations:

sphrPak[s_, r_, n_Integer?Positive, lim_Integer?Positive] :=
Module[{sh = s/2, k, p0, pc, sphList},

(* generate initial sphere within box *)
While[
p0 = RandomReal[{-sh, sh}, 3];
And @@ Thread[Map[Min[sh - #, sh + #] &, p0] < r]];
sphList = {p0};

k = 0;
While[Length[sphList] < n && k < lim,
(* center for new sphere chosen in random direction *)
pc = p0 + 2 r Normalize[RandomVariate[NormalDistribution[], 3]];
If[(* is sphere within the cube? *)
PolyhedronData["Cube", "RegionFunction"] @@
ScalingTransform[ConstantArray[1/(s - 2 r), 3]][pc] &&
(* does the sphere not overlap with other spheres? *)
(And @@ Thread[Map[EuclideanDistance[pc, #] &, sphList] >= 2 r]),
k = 0; AppendTo[sphList, pc]; p0 = pc,
(* else *) ++k];
];

If[k == lim, Print[StringForm["Only 1 spheres were generated.", Length[sphList]]]];

Graphics3D[Sphere[sphList, r], Axes -> Automatic,
PlotRange -> {{-sh, sh}, {-sh, sh}, {-sh, sh}}]];

sphrPak[s_, r_, n_Integer] := sphrPak[s, r, n, Quotient[2 n, 3]]


Try it out:

BlockRandom[SeedRandom[1024, Method -> "MersenneTwister"]; (* for reproducibility *)
sphrPak[20, 0.7, 1500]]
Only 725 spheres were generated. Of course, a different attempt might yield a result with more spheres within the cube.

More attempts:

BlockRandom[SeedRandom[4092, Method -> "MersenneTwister"];
sphrPak[20, 0.7, 1500]]
Only 805 spheres were generated. BlockRandom[SeedRandom[2012, Method -> "MersenneTwister"];
sphrPak[20, 0.7, 1500]]
Only 932 spheres were generated. Here's a variation that attempts to do more than one bounded random walk, where the starting point of the $n$th random walk is a randomly chosen point from the $(n-1)$-th random walk:

sphrPak2[s_, r_, n_Integer?Positive, walks_Integer?Positive, its_Integer?Positive] :=
Module[{sh = s/2, j, k, p0, pc, sphList},
While[
p0 = RandomReal[{-sh, sh}, 3];
And @@ Thread[Map[Min[sh - #, sh + #] &, p0] < r]];
sphList = {p0}; j = 0;
While[Length[sphList] < n && j < walks,
k = 0;
While[Length[sphList] < n && k < its,
pc = p0 + 2 r Normalize[RandomVariate[NormalDistribution[], 3]];
If[(PolyhedronData["Cube", "RegionFunction"] @@
ScalingTransform[ConstantArray[1/(s - 2 r), 3]][pc]) &&
(And @@ Thread[Map[EuclideanDistance[pc, #] &, sphList] >= 2 r]),
k = 0; AppendTo[sphList, pc]; p0 = pc,
++k];
];
If[Length[sphList] < n,
p0 = RandomChoice[Map[EuclideanDistance[#, sphList[]] &,
sphList] -> sphList]];
++j;
];
If[j == walks, Print[StringForm["Only 1 spheres were generated.",
Length[sphList]]]];
Graphics3D[Sphere[sphList, r], Axes -> Automatic,
PlotRange -> {{-sh, sh}, {-sh, sh}, {-sh, sh}}]]


Try it out:

BlockRandom[SeedRandom[85, Method -> "Rule50025CA"];
sphrPak2[20, 0.7, 1500, 1, 1000]]
Only 662 spheres were generated. BlockRandom[SeedRandom[85, Method -> "Rule50025CA"];
sphrPak2[20, 0.7, 1500, 20, 1000]]
Only 1335 spheres were generated. Why do you need to check if the center AND the hull are inside the cube? – Dr. belisarius – 2012-10-12T12:14:18.057

1@bel: I'm exploiting the fact that And[] does short-circuit evaluation; if it fails the easy test, then the other tests are not touched. On the other hand, I suppose I can merge those first two tests; wait... – J. M.'s ennui – 2012-10-12T12:19:19.487

6

I will present another brute-force method for densely packing spheres that works by randomly depositing the spheres in successive "layers". The first layer is built by solving the two-dimensional analog of the problem (that is, packing disks in a square). Having constructed the first layer, successive layers are built up by dropping new spheres, checking that the incoming spheres do not intersect with spheres in previous layers (as well as with other spheres in the same layer). For the second layer, the new spheres are compared against spheres in the first layer, while for the third and succeeding layers, the new sphere is tested for intersections with the two previous layers. All this is done up until we hit the top of the cube.

Here is a Mathematica implementation of the idea presented above:

sphrPak[s_, r_, its_Integer: 2000] :=
Module[{j, k, lap, lap2, layer, layers, p0, pc, sphrs, z},

(* pack disks into square *)

(* pick first disk *)
While[
p0 = RandomReal[{r, s - r}, 2];
And @@ Thread[Map[Min[#, s - #] &, p0] < r]];

(* successively pack disks into square *)
layer = {p0}; k = 0;
While[k < its,
pc = RandomReal[{r, s - r}, 2]; (* candidate disk *)
If[And @@ Thread[Map[EuclideanDistance[pc, #] &, layer] >= 2 r],
k = 0; AppendTo[layer, pc];,
++k];
];

(* build and store first layer *)
layers = {PadRight[layer, {Length[layer], 3}, r]};

(* build successive layers *)
j = 2;
While[True,

(* previous layers *)
lap = layers[[j - 1]]; If[j > 2, lap2 = layers[[j - 2]]];

(* find "seed" sphere for layer *)
p0 = {}; k = 0;
While[k < its,
z = RandomReal[{j, j + 1} r];
pc = Append[RandomReal[{r, s - r}, 2], z]; ++k; (* candidate seed *)
If[z <= s - r &&
(And @@ Thread[Map[EuclideanDistance[pc, #] &, lap] >= 2 r]),
If[j > 2,
If[And @@ Thread[Map[EuclideanDistance[pc, #] &, lap2] >= 2 r],
p0 = pc; Break[]],
p0 = pc; Break[]]]];

If[p0 === {}, Break[]]; (* no seed found *)

(* start building new layer *)
layer = {p0}; k = 0;
While[k < its,
z = RandomReal[{j, j + 1} r];
pc = Append[RandomReal[{r, s - r}, 2], z]; (* candidate point *)
If[z <= s - r &&
(And @@ Thread[Map[EuclideanDistance[pc, #] &, layer] >= 2 r]) &&
(And @@ Thread[Map[EuclideanDistance[pc, #] &, lap] >= 2 r]),
If[j > 2,
If[And @@ Thread[Map[EuclideanDistance[pc, #] &, lap2] >= 2 r],
k = 0; AppendTo[layer, pc], ++k],
k = 0; AppendTo[layer, pc]], ++k];
];
AppendTo[layers, layer];
++j];
Print[StringForm["1 spheres generated.", Length[sphrs = Join @@ layers]]];
Graphics3D[Sphere[sphrs, r], Axes -> Automatic, PlotRange -> {{0, s}, {0, s}, {0, s}}]]


Try it out:

BlockRandom[SeedRandom[2012, Method -> "MersenneTwister"]; (* for reproducibility *)
sphrPak[20, 0.7]]
1996 spheres generated. BlockRandom[SeedRandom[1024, Method -> "MersenneTwister"];
sphrPak[20, 0.7]]
2041 spheres generated. being nit picky, if i understand this scheme correctly it seems you are introducing a certain non-randomness in that you have a layered structure with constant z planes. That may be important depending on what you want to do with the result. Perhaps throw in a small random z within each layer. Might need to check two layers below if you do that though. – george2079 – 2012-10-15T12:25:41.570

No @george, the heights at each layer (except the first) are not constant; that is, unless z = RandomReal[{j, j + 1} r] is doing something peculiar that you know and I don't. Yes, I check two layers below as well; that's what the If[j > 2,(* stuff *)] sections of the seeding and layer-building parts are for. – J. M.'s ennui – 2012-10-15T12:30:33.757

sorry..in my head that was outside the loop. Now i'm not sure that seeding the layer is necessary, you could find the mean z for each layer and move up a fixed amount. In the end i expect that seeding is a small part of the time though. – george2079 – 2012-10-15T13:08:21.877

2

So using another random walk-like method (well I tried two which I realized are equivalent) I was able to do this pretty quickly and compactly:

generateRandomAnchoring[] :=

With[{p = RandomReal[{r, m - r}, 3], c = RandomChoice@points,
d = (2 r)},
c + d*Normalize@(c - p)
];
generateRandomWalking[] :=

With[{p = RandomChoice[points], \[Theta] =
RandomReal[2 \[Pi]], \[Phi] = RandomReal[2 \[Pi]], d = (2 r)},
p + {d*Cos[\[Phi]] Sin[\[Theta]], d*Sin[\[Theta]] Sin[\[Phi]],
d*Cos[\[Theta]]}
];
tryPoint[s_] :=
If[AllTrue[s, r < # < m - r &] &&
AllTrue[points, (Norm[# - s] >= 2 r) &], AppendTo[points, s]];


Both of these take a random pregenerated point and attach a point to it at a random position, just in one case it generates a spherical polar coordinate and in the other it just attaches a randomly generated point. They're entirely equivalent as best I can reason out.

Then we'll just initialize the state and update the point generating code from wxffles answer so that it also spits out the time elapsed:

r = .7;
m = 20;
points = {};
tryPoint@RandomReal[{r, m - r}, 3];

With[{now = Now},
Monitor[While[Length[points] < 1825,
tryPoint@generateRandomAnchoring[]
], Column@{Now - now, Length@points}]
]


It runs for about a minute before hitting the 1825 mark. And we'll use the same viewer code as wxffles too:

cube = {Opacity[0.3], Cuboid[{0, 0, 0}, {20, 20, 20}]};
Graphics3D[{cube, Sphere[#, 0.7] & /@ points}, Boxed -> False]


And this is what we see: Then just adjust the number of spheres and timing to your liking and you have a your packed cube.

### Edit

Note that one could also add a weight-tracking scheme for the random choice of points (in the second of the methods). I.e., points that have been chosen and for which the addition failed decrease in weight by some percentage. The likelihood that this will be a significant improvement is low, but it could eke out some performance gains.