Okay, I could not resist to give it another try. This is somewhat related to my research anyway. I use a new answer since this is going to be quite a different approach and there might still something to learn from the old approach.
This time, we write a heavily optimized CompiledFunction, called CapsuleIntersectingQ that (hopefully correctly) detects wether two capsules with given centerlines and given radii intersect. We do that by discussing a convex quadratic function f (see below) in two variables that describes quadratic distances of points on the infinite lines through the centerlines of the capsules. Actually, we need the minimizers of this function on the unit square, so we need to discuss several cases (9 cases to be precise, depending on wether the minimum lies in the interior, on an edge, or in a corner of the unit square). The generic case can solved by applying Newton's method symbolically (only one step is needed, since f is quadratic). I am not 100% sure wether I handle the boundary cases correctly, so please report if you find any intersecting capsules where they should not be or any other oddities.
(* a CompiledFunction that is supposed to test pairs of capsules for intersection *)
CapsuleIntersectingQ =
Quiet[Block[{pp, qq, f, Df, DDf, x, y, S, T, S0, S1, T0, T1, s, t, p1, p2, p3, p4, p5, p6, q1, q2, q3, q4, q5, q6},
pp = {{p1, p2, p3}, {p4, p5, p6}};
qq = {{q1, q2, q3}, {q4, q5, q6}};
f = {x, y} \[Function] Evaluate[Total[({(1 - x), x}.pp - {(1 - y), y}.qq)^2]];
Df = {x, y} \[Function] Evaluate[D[f[x, y], {{x, y}, 1}]];
DDf = {x, y} \[Function] Evaluate[D[f[x, y], {{x, y}, 2}]];
{S, T} = LinearSolve[DDf[0, 0], -Df[0, 0]];
S0 = -Df[0, 0][[1]]/DDf[0, 0][[1, 1]];
S1 = -Df[0, 1][[1]]/DDf[0, 1][[1, 1]];
T0 = -Df[0, 0][[2]]/DDf[0, 0][[2, 2]];
T1 = -Df[1, 0][[2]]/DDf[1, 0][[2, 2]];
With[{
sint = N@S, tint = N@T,
F00 = N@f[0, 0], F01 = N@f[0, 1], F10 = N@f[1, 0],
F11 = N@f[1, 1], Fst = N@f[S, T],
s0 = N@S0, Fs0 = N@f[S0, 0],
s1 = N@S1, Fs1 = N@f[S1, 1],
t0 = N@T0, F0t = N@f[0, T0],
t1 = N@T1, F1t = N@f[1, T1]
},
Compile[{{p, _Real, 2}, {r1, _Real}, {q, _Real, 2}, {r2, _Real}},
Block[{s, t, c, p1, p2, p3, p4, p5, p6, q1, q2, q3, q4, q5, q6},
p1 = Compile`GetElement[p, 1, 1];
p2 = Compile`GetElement[p, 1, 2];
p3 = Compile`GetElement[p, 1, 3];
p4 = Compile`GetElement[p, 2, 1];
p5 = Compile`GetElement[p, 2, 2];
p6 = Compile`GetElement[p, 2, 3];
q1 = Compile`GetElement[q, 1, 1];
q2 = Compile`GetElement[q, 1, 2];
q3 = Compile`GetElement[q, 1, 3];
q4 = Compile`GetElement[q, 2, 1];
q5 = Compile`GetElement[q, 2, 2];
q6 = Compile`GetElement[q, 2, 3];
c = Min[{F00, F01, F10, F11}];
If[c <= (r1 + r2)^2,
True,
s = sint; t = tint;
If[0. <= s <= 1. && 0. <= t <= 1.,
c = Fst,
If[0. <= s0 <= 1., c = Min[c, Fs0]];
If[0. <= s1 <= 1., c = Min[c, Fs1]];
If[0. <= t0 <= 1., c = Min[c, F0t]];
If[0. <= t1 <= 1., c = Min[c, F1t]];
];
c <= (r1 + r2)^2
]
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True
]
]
]];
For completeness, here is our method to produce random rotation matrices from random vectors.
getRotationMatrix =
Compile[{{u, _Real, 1}},
Block[{uu, r, cosr, sinr}, uu = u[[1]]^2 + u[[2]]^2 + u[[3]]^2;
r = Sqrt[uu];
cosr = Cos[r];
sinr = Sin[r];
{{(u[[1]]^2 + cosr u[[2]]^2 + cosr u[[3]]^2)/
uu, (u[[1]] u[[2]] - cosr u[[1]] u[[2]] - u[[3]] r sinr)/
uu, (u[[1]] u[[3]] - cosr u[[1]] u[[3]] + u[[2]] r sinr)/
uu}, {(u[[1]] u[[2]] - cosr u[[1]] u[[2]] + u[[3]] r sinr)/
uu, (cosr u[[1]]^2 + u[[2]]^2 + cosr u[[3]]^2)/
uu, (u[[2]] u[[3]] - cosr u[[2]] u[[3]] - u[[1]] r sinr)/
uu}, {(u[[1]] u[[3]] - cosr u[[1]] u[[3]] - u[[2]] r sinr)/
uu, (u[[2]] u[[3]] - cosr u[[2]] u[[3]] + u[[1]] r sinr)/
uu, (cosr u[[1]]^2 + cosr u[[2]]^2 + u[[3]]^2)/uu}}],
CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True];
And here is the main program. Note that I remove everything related to RegionDistance.
(*number of capsules to stuff into the box*)
n = 10000;
(*maximal number of trials to stuff a random capsule into the box*)
maxiters = 1000;
(*edge length of the box*)
l = 100.;
(*length of capsules*)
L = 15;
(*radius of capsules*)
r = 3.;
(*prototypical centerline of capsules; new capsules are generated by rotating and translating it randomly*)
p0 = Developer`ToPackedArray[{-L/2 {1., 0., 0.}, L/2 {1., 0., 0.}}];
(*a bag for rule them all*)
ptbag = Internal`Bag[{}];
(*initialize the distance function so that the first capsule will be feasible*)
distToCenterlines = Function[x, 2 l, Listable];
(*a counter for the capsules*)
ncapsules = 0;
(*the box (for graphical purposes only*)
cube = Cuboid[-l/2 {1., 1., 1.}, l/2 {1., 1., 1.}];
(*the box where the centerlines have to be contained in*)
safetycube = Cuboid[-(l/2 - r) {1., 1., 1.}, (l/2 - r) {1., 1., 1.}];
(*generate a function that test if a point is contained in safetycube*)
insafetycubeQ = RegionMember[safetycube];
(*this is the amount of random data to generate at once*)
m = 1000;
(*c is a counter keeping track of the amount of random data (matlist and shiftlist below) that has been used already;initialize the counter so that a refresh will be kicked off immediately in the main loop*)
c = m;
(*use a normal distribution to generate random rotation vectors*)
distro = MultinormalDistribution[{0., 0., 0.}, DiagonalMatrix[{1., 1., 1.}]];
Dynamic[{j, iter}]
(*finally,the main loop*)
Do[iter = 0;
(*Boolean that is supposed to be equal to False if the new capsule is feasible*)
b = True;
(*generate random capsules until you find a feasible one*)
While[b,
++iter;
If[iter >= maxiters, Break[]];
(*check if the randomness reservoir is not empty*)
If[c >= m,(*if empty,refill the reservoir*)
matlist = getRotationMatrix[RandomVariate[distro, m]];
shiftlist = RandomReal[{r - l/2, l/2 - r}, {m, 3}];
c = 0;
];
c++;
(*generate the random capsule*)
p = p0.matlist[[c]] + ConstantArray[shiftlist[[c]], 2];
(*perform the cheap test against the box first*)
b = ! And @@ (insafetycubeQ[p]);
If[! b && ncapsules > 0,
centerlines = ArrayReshape[Internal`BagPart[ptbag, All], {ncapsules, 2, 3}];
b = Or @@ CapsuleIntersectingQ[centerlines, r, p, r]
];
];
If[iter >= maxiters,
Print["I tried my best but with ", maxiters, " trials, I was able to stuff only ", ncapsules, " capsules into the box."];
Break[]
];
(*yes,we have one more capsule...*)
ncapsules += 1;
(*stuff the centerline of new capsule into the bag*)
Internal`StuffBag[ptbag, Flatten[p], 6];
, {j, 1, n}]
This time, I let it run until nothing seems to fit in any more. Here is the result:
Graphics3D[{
FaceForm[{Darker@Blue, Opacity[.1]}], Specularity[White, 30],
cube, safetycube,
FaceForm[{Darker@Orange, Opacity[1.0]}], Specularity[White, 30],
Map[CapsuleShape[#, r] &, centerlines]
},
Lighting -> "Neutral", PlotRange -> Transpose[List @@ cube],
Boxed -> False, SphericalRegion -> True
]
CapsuleShape[]in Mathematica. – J. M.'s missing motivation Nov 08 '17 at 17:31