Creating random configurations of spherocylinders or cylinders

Question

About the setting:

We have a 3D simulation box with side $l$ and our catesian coordinate system is set with its origin at the centre of the box. We have a number $N$ of spherocylinders of aspect ratio $L/D=20,$ where $D$ is our basic unit of length, namely, the base diameter of the cylinder (or the spherical caps) and $L$ the length of the cylinder (without the caps). We can take an example of a box with $l=100.$

Question:

The aim is to generate random configurations for the $N$ spherocylinders, such that none of them are in overlap, so not even touching, nor are they in contact with the walls of the box. This is useful for studying packing properties of such anisotropic objects, in contrast to the well-known sphere packing studies. I am curious to find out whether in Mathematica, there are existing libraries that can help one get a head start for implementing a routine that would generate valid configurations. A configuration is simply a list with the coordinates of the centres of the particles and their respective orientation coordinates, so 6 numbers per particle.

Answer 1

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
 ]

Answer 2

Any approach to this will depend on efficiently deciding if two capsules overlap. You could do this using the built-in RegionDisjoint

capsulesOverlapRegion[cap1 : CapsuleShape[{p1a_List, p1b_List}, d1_], 
  cap2 : CapsuleShape[{p2a_List, p2b_List}, d2_]] := ! 
  RegionDisjoint[
   BoundaryDiscretizeRegion[cap1],
   BoundaryDiscretizeRegion[cap2]
   ]

but this isn't terribly fast.

SeedRandom[42];
{c1, c2} = {CapsuleShape[RandomReal[1, {2, 3}], .25], 
   CapsuleShape[RandomReal[1, {2, 3}], .25]};

capsulesOverlapRegion[c1, c2] // RepeatedTiming
(* {0.033, True} *)

Instead, I will use the compiled function described in this post, which measures the minimum distance between two line segments. If this distance is greater than the sum of the diameters, the capsules do not intersect. Here I'll drop the CapsuleShape wrapper, and represent a given capsule as a list of {point1, point2, diameter}

distSegToSegComp = << "https://git.io/vDJMK";
capsulesOverlap[CapsuleShape[{p1a_List, p1b_List}, d1_], 
  CapsuleShape[{p2a_List, p2b_List}, d2_]] := (d1 + d2) > 
  distSegToSegComp[{p1a, p1b}, {p2a, p2b}]

capsulesOverlap[c1, c2] // RepeatedTiming
(* {9.*10^-6, True} *)

Now you can use a loop, generating a random capsule at each step and rejecting it if it overlaps with another capsule.

With a speedy test, we can just use a While loop to get the required number of shapes. Here is an example that generates 300 pseudorandomly oriented capsules, with length 0.3 and diameter 0.05, all within a box extending from {-1,-1,-1} to {1,1,1}.

length = 0.3;
diameter = 0.05;
nCapsuleGoal = 300;
maxAttempts = 3000;

boxLowerPoint = {-1, -1, -1};
boxUpperPoint = {1, 1, 1};
box = BoundaryDiscretizeRegion[
   Cuboid[diameter + boxLowerPoint, boxUpperPoint - diameter]];
pointInBox = 
  And @@ Thread[(diameter + boxLowerPoint) < # < (-diameter + 
        boxUpperPoint)] &;

distSegToSegComp = << "https://git.io/vDJMK";
capsulesOverlap[{p1a_List, p1b_List, d1_}, {p2a_List, p2b_List, 
   d2_}] := (d1 + d2) > distSegToSegComp[{p1a, p1b}, {p2a, p2b}]

capsuleInBox[{p1_List, p2_List, d_}] := 
  pointInBox[p1] && pointInBox[p2];

\[ScriptCapitalD] = 
  MatrixPropertyDistribution[r.{0, 0, length}, 
   r \[Distributed] CircularRealMatrixDistribution[3]];

capsules = {};
nCapsules = 0;
nAttempts = 0;
randomCapsule[] := 
  With[{start = RandomPoint[box]}, {start, 
    start + RandomVariate[\[ScriptCapitalD]], diameter}];

While[nCapsules < nCapsuleGoal && nAttempts < maxAttempts, 
  capsule = randomCapsule[];
  If[capsuleInBox[capsule] && 
    AllTrue[capsules, ! capsulesOverlap[capsule, #] &], 
   AppendTo[capsules, capsule];
   nCapsules++;];
  nAttempts++;];

visualize the results with

Graphics3D[{CapForm[Round], Tube[{#1, #2}, #3] & @@@ capsules, Blue, 
  Opacity[0.1], Cuboid[boxLowerPoint, boxUpperPoint]}, Boxed -> False]

Answer 3

Disclaimer

In the meantime, user929304 found out that my collision test is incorrect. I posted a new answer with a rather new approach.

Strategy

I use a Bag to store the lines between the two sphere centers of each existing capsule. When generating a new capsule I use

1. RegionMember to check wether a new capsule fits into the box and
2. RegionDistance of a MeshRegion containing all centerlines of existing cpasules to check for collisions between the new capsule and the existent capsules.

It is advantageous to pack all existing capsules into a single MeshRegion and to use its RegionDistanceFunction since internally, Mathematica uses a k-d-tree (or a similar data structure): This avoids computing the distance from a new capsule to each preexisting capsule.

Here is the code along with several comments:

(*number of capsules to stuff into the box*) 
n = 100; 
(*maximal number of trials to stuff a random capsule into the box*)
maxiters = 1000;
(*edge length of the box*)
l = 10.;
(*length of capsules*) 
L = 0.5; 
(*radius of capsules*) 
r = 0.125;
(*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[{}];
(*generating edge combinatorics in advance safes much time later*)
edges = Transpose[{2 Range[n] - 1, 2 Range[n]}]; 
(*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];
(*the following is a 1000 times faster substitute for RotationMatrix; sorry for the mess*)
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
   ];

(*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;
 (*a 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, b = Min[distToCenterlines[p]] <= 2. r];
 ];
 (*when the while loop exits, p represents a feasible capsule or we used maxiters trials*)
 If[iter >= maxiters,
  Print["I tried my best but I could only stuff ", ncapsules, 
   " into the box."];
  Break[]
  ];
 (*yes,we have one more capsule...*)
 ncapsules += 1;
 (*stuff the centerline end of new capsule into the bag*)
 Internal`StuffBag[ptbag, Flatten[p], 6];
 (*recompute distance function to existing centerlines*)
 distToCenterlines = RegionDistance[
   MeshRegion[
    Partition[Internal`BagPart[ptbag, All], 3],
    Line[edges[[1 ;; ncapsules]]]
    ]
   ];
 , {j, 1, n}]

Extract the information as desired from by the OP:

capsulecenterlines = ArrayReshape[Internal`BagPart[ptbag, All], {ncapsules, 2, 3}];
capsulecenters =  0.5 (capsulecenterlines[[All, 2]] + capsulecenterlines[[All, 1]]);
capsuledirections = (capsulecenterlines[[All, 2]] - capsulecenterlines[[All, 1]])/L;

These can be exported, e.g., to CSV files...

Export[FileNameJoin[{$HomeDirectory, "capsulecenterdata.csv"}], capsulecenters]
Export[FileNameJoin[{$HomeDirectory, "capsuledirectiondata.csv"}], capsuledirections]

... and also reimported (with check for correctness):

Import[FileNameJoin[{$HomeDirectory, "capsulecenterdata.csv"}], "Data"] == capsulecenters
Import[FileNameJoin[{$HomeDirectory, "capsuledirectiondata.csv"}], "Data"] == capsuledirections

(* True *)
(* True *)

And here is the result in graphical form:

Graphics3D[{
  FaceForm[{Darker@Blue, Opacity[.1]}], Specularity[White, 30],
  cube, safetycube,
  FaceForm[{Darker@Orange, Opacity[1.0]}], Specularity[White, 30], 
  Map[CapsuleShape[#, r] &, capsulecenterlines]
  },
 Lighting -> "Neutral",
 PlotRange -> Transpose[List @@ cube],
 Boxed -> False,
 SphericalRegion -> True
 ]

Suggestions for improvement:

I think it would be possible to speed things up if one could get a grip directly on the k-d-Tree. In the current implementation, the tree is rebuilt completely after a new capsule is added. But there has to be cheaper a method to create the new k-d-Tree by modifying the old one.

I am not 100% sure if the generated rotation matrices a uniformly distributed over rotation group $\operatorname{SO}(3)$. On the other hand, the generated distribution is rotation invariant, isn't it? So we sample from the normalized Haar measure, aren't we?

Edit

In the meantime, I made several improvements:

• The greatest effect came from defining the edge combinatorics right in the beginning so that a lot of messy copying, packing and unpacking of arrays can be skipped.

• Generate random data in larger chunks since this is a bit faster. That has also some effect.

• Tried to use an Internal`Bag for collecting the endpoints of the center lines. The effect however is almost not measurable.

• Replaced RotationMatrix with a much quicker compiled function.