# How to generate random non-overlapping cylinders of fixed radius and height in a fixed cube?

I have seen these two questions (first) (second) which are related to generating non-overlapping cylinders in a cube. I am trying to adapt them to my goal which is to generate random non-overlapping fixed cylinders (radius=10nm and height=20nm) in a cube of (100*100*100 nm).

The number of cylinders would be 20, 40 and 60 for the three cases that I would like to build. This makes the volume fraction in the range of around 0.1, 0.2 and 0.3 for the three cases.

Then at the end I would like to export the coordinates of each cylinders (x1,y1,z1) &(x2,y2,z2) in a file.

I have been using the following reference code as an example and trying to modify this.

    p1.p1 + p2.p2 - 2 p1.p2
((p1.p1 + p2.p2 - 2 p1.p2) /. {p1 -> p1i + dp1 t1, p2 -> p2i + dp2 t2}) // tf // Expand
(* p1i.p1i - 2 p1i.p2i + p1i.(dp1 t1) - 2 p1i.(dp2 t2) + p2i.p2i +
p2i.(dp2 t2) + (dp1 t1).p1i - 2 (dp1 t1).p2i + (dp1 t1).(dp1 t1) -
2 (dp1 t1).(dp2 t2) + (dp2 t2).p2i + (dp2 t2).(dp2 t2) *)

tf[e_] := e /. Dot[z1_ + z2_, z3_ + z4_] ->
Dot[z1, z3] + Dot[z2, z3] + Dot[z1, z4] + Dot[z2, z4]

Map[# /. z_ :> (z /. t1 -> 1) t1^Count[z, t1, 4] &,
Map[# /. z_ :> (z /. t2 -> 1) t2^Count[z, t2, 4] &, %]]
(* t1^2 dp1.dp1 - 2 t1 t2 dp1.dp2 + t1 dp1.p1i - 2 t1 dp1.p2i + t2^2 dp2.dp2 +
t2 dp2.p2i + t1 p1i.dp1 - 2 t2 p1i.dp2 + p1i.p1i - 2 p1i.p2i + t2 p2i.dp2 + p2i.p2i *)

FullSimplify[Solve[{D[%, t1] == 0, D[%, t2] == 0}, {t1, t2}],
TransformationFunctions -> {Automatic, tf1}]
(* {{t1 -> (dp2.dp1 (-dp2.p1i + dp2.p2i) +
dp2.dp2 (p1i.dp1 - p2i.dp1))/((dp1.dp2)^2 - dp1.dp1 dp2.dp2),
t2 -> (dp1.dp1 (-dp2.p1i + dp2.p2i) +
dp2.dp1 (p1i.dp1 - p2i.dp1))/((dp1.dp2)^2 - dp1.dp1 dp2.dp2)}} *)

tf1[e_] := e /. Dot[z1_, z2_] -> Dot[z2, z1]

int[cyl1_, cyl2_] :=
Module[{p1i = cyl1[[1, 1]], dp1 = cyl1[[1, 2]] - cyl1[[1, 1]], r1 = cyl1[],
p2i = cyl2[[1, 1]], dp2 = cyl2[[1, 2]] - cyl2[[1, 1]], r2 = cyl2[],
loc, t1, t2}, loc = {t1 -> (dp2.dp1 (-dp2.p1i + dp2.p2i) +
dp2.dp2 (p1i.dp1 - p2i.dp1))/((dp1.dp2)^2 - dp1.dp1 dp2.dp2),
t2 -> (dp1.dp1 (-dp2.p1i + dp2.p2i) +
dp2.dp1 (p1i.dp1 - p2i.dp1))/((dp1.dp2)^2 - dp1.dp1 dp2.dp2)};
(-r2/Norm[dp2] < (t1 /. loc) < 1 + r2/Norm[dp2]) &&
(-r1/Norm[dp1] < (t2 /. loc) < 1 + r1/Norm[dp1]) &&
((t1^2 dp1.dp1 - 2 t1 t2 dp1.dp2 + t1 dp1.p1i - 2 t1 dp1.p2i + t2^2 dp2.dp2 +
t2 dp2.p2i + t1 p1i.dp1 - 2 t2 p1i.dp2 + p1i.p1i - 2 p1i.p2i + t2 p2i.dp2 +
p2i.p2i) /. loc) < (r1 + r2)^2]

cylinders = Table[{RandomReal[{-100, 100}, {2, 3}], RandomReal}, {100}];
nint = ParallelTable[Or @@ Table[int[cylinders[[i]], cylinders[[j]]], {j, i + 1, 100}],
{i, 100}];
c = Cases[{cylinders, nint} // Transpose, {z_, False} -> z, {1}];
Graphics3D[{EdgeForm[None], Directive[Opacity@RandomReal[{.4, .9}], Hue[RandomReal[]]],
Cylinder[First@#, Last@#]} & /@ c, Boxed -> False, ImageSize -> 800]


• @bbgodfrey Could you check this one? It's almost the same one you did before. Aug 4 '17 at 12:45
• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Aug 6 '17 at 2:02

Edited to Improve Definition of cylinders

This problem differs from the two referenced in the question in three important ways. The cylinders are precisely confined within a box, the filling factor is much larger, and the aspect ratio of the cylinders is essentially unity. With perfect ordering, 125 cylinders can fit into the box. As soon as only a few are placed randomly, the maximum number of remaining cylinders that fit drops precipitously. To take these considerations into account, int (which determines whether two cylinders are non-overlapping) was modified from the referenced answers.

int[cyl1_, cyl2_] := Module[{p1i = cyl1[[1, 1]], dp1 = cyl1[[1, 2]] - cyl1[[1, 1]],
r1 = cyl1[], p2i = cyl2[[1, 1]], dp2 = cyl2[[1, 2]] - cyl2[[1, 1]],
r2 = cyl2[], loc, t1, t2},
loc = {t1 -> Clip[(dp2.dp1 (-dp2.p1i + dp2.p2i) + dp2.dp2 (p1i.dp1 - p2i.dp1))/
((dp1.dp2)^2 - dp1.dp1 dp2.dp2), {- r1/Norm[dp1], 1 + r1/Norm[dp1]}],
t2 -> Clip[(dp1.dp1 (-dp2.p1i + dp2.p2i) + dp2.dp1 (p1i.dp1 - p2i.dp1))/
((dp1.dp2)^2 - dp1.dp1 dp2.dp2), {- r2/Norm[dp2], 1 + r2/Norm[dp2]}]};
((t1^2 dp1.dp1 - 2 t1 t2 dp1.dp2 + t1 dp1.p1i - 2 t1 dp1.p2i + t2^2 dp2.dp2 +
t2 dp2.p2i + t1 p1i.dp1 - 2 t2 p1i.dp2 + p1i.p1i - 2 p1i.p2i + t2 p2i.dp2 +
p2i.p2i) /. loc) > 2 (r1 + r2)^2]


Note that this function is only approximate and occasionally determines that two cylinders overlap when they do not. A precise determination could be obtained from

Volume@RegionIntersection[Cylinder @@ cylinders[[i]], Cylinder @@ cylinders[[j]]] == 0


but is orders of magnitude slower than int. The number of falsely excluded cylinders is not large. In addition, a second function is introduced to identify and eliminate cylinders that extend slightly beyond the 100 x 100 x 100 box. It is exact and quite fast.

inb[cyl_] := Module[{cos, d, rl = Last[cyl], hlf = Norm[cyl[[1, 1]] - cyl[[1, 2]]]/2},
And @@ Map[(cos = (cyl[[1, 1, #]] - cyl[[1, 2, #]])/(2 hlf);
d = hlf Abs[cos] + rl Sqrt[1 - cos^2];
(-50 + d < (cyl[[1, 1, #]] + cyl[[1, 2, #]])/2 < 50 - d)) &, {1, 2, 3}]]


The code to determine the non-overlapping cylinders then is,

(SeedRandom; ict = 2000000; rl = 10; hlf = 10; n = 60;
cylinders = Table[s = RandomReal[{-50 + rl, 50 - rl}, {3}]; zz = RandomReal[{-hlf, hlf}];
phi = RandomReal[{-Pi, Pi}];
{{s - {Sqrt[hlf^2 - zz^2] Cos[phi], Sqrt[hlf^2 - zz^2] Sin[phi], zz},
s + {Sqrt[hlf^2 - zz^2] Cos[phi], Sqrt[hlf^2 - zz^2] Sin[phi], zz}}, rl}, ict];
cylinders = Pick[cylinders, inb /@ cylinders];
c = {cylinders[]};
Do[If[(And @@ Table[int[cylinders[[i]], c[[j]]], {j, Length[c]}]) && (Length[c] < n),
AppendTo[c, cylinders[[i]]]], {i, 2, Length[cylinders]}];
Length[c]) // AbsoluteTiming

(* {2412.55, 16} *)

Graphics3D[{EdgeForm[None], Directive[Opacity@RandomReal[{.4, .9}], Hue[RandomReal[]]],
Cylinder @@ #} & /@ c, ImageSize -> 800, PlotRange -> {{-50, 50}, {-50, 50}, {-50, 50}}] Although some cylinders in the figure may appear to be overlap, rotating the figure within Mathematica shows that they are not. Significantly increasing the number of randomly placed cylinders will be difficult unless rl is decreased.

Result for rl == hlf == 7.

Reducing rl and hlf to 7 nearly triples the number of non-overlapping cylinders obtained, although slowly.

(* {6881.68, 45} *) Result for 2:1 aspect ratio

At the request of the OP in a comment below, the code above has been modified modestly to accommodate independent values for the cylinder radius, rl, and half length, hlf. Choosing

ict = 10000; rl = 5; hlf = 10; n = 60


and the code otherwise as above yields

(* {29.1495, 42} *) At the OP's request, the parameter n has been added to the code above. It specifies the maximum number of non-overlapping cylinders to be generated. Of course, fewer may be generated, depending on the values of the other parameters, ict, rl, and hlf.
• I found and fixed a typo in inb. The simplest way to export the cylinder data as text is Export["C:\\Temp\\cyls.txt", c], where the first argument is the location to which the data is to be exported. By the way, I now know how to make int exact, which will permit a few more cylinders to be packed into the box. It will take me a while to implement the code. Please do not forget to upvote and accept my answer, if it meets your needs. Aug 6 '17 at 15:54