3
$\begingroup$

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[[2]], 
        p2i = cyl2[[1, 1]], dp2 = cyl2[[1, 2]] - cyl2[[1, 1]], r2 = cyl2[[2]], 
        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[5]}, {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]

Any heads up would be helpful.

Improve this question
$\endgroup$
2
  • $\begingroup$ @bbgodfrey Could you check this one? It's almost the same one you did before. $\endgroup$
    – skonda2
    Aug 4, 2017 at 12:45
  • $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Aug 6, 2017 at 2:02

1 Answer 1

Reset to default
4
$\begingroup$

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[[2]], p2i = cyl2[[1, 1]], dp2 = cyl2[[1, 2]] - cyl2[[1, 1]], 
    r2 = cyl2[[2]], 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[1492]; 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[[1]]};
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}}]

enter image description here

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} *)

enter image description here

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} *)

enter image description here

Minor Addition

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.

Improve this answer
$\endgroup$
10
  • $\begingroup$ I am having a hard time putting together the whole code. Some error pops out when I attach everything. The pictures you generated are pretty much what I am looking for. Do you know, how export the co-ordinates in a text file (x1,y1,z1) & (x2,y2,z2)? $\endgroup$
    – skonda2
    Aug 6, 2017 at 15:11
  • $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Aug 6, 2017 at 15:54
  • $\begingroup$ What does (* {29.8132, 78} *) mean? I know 78 is the number of cylinders. What does 29.81 mean, is it time? Also, ict value of 5000, do you think it should be 500000? $\endgroup$
    – skonda2
    Aug 6, 2017 at 20:23
  • $\begingroup$ The code runs now, however I can see multiple overlapping cylinders in mathematica. Is it something to do with the int or inb function? In addition, Is there a way we could incorporate a input parameter "n", which decides the number of cylinders inside the cube? $\endgroup$
    – skonda2
    Aug 9, 2017 at 1:02
  • $\begingroup$ I am not looking for high packing fractions. So, a reasonable "n" would be useful. At the end, I would have a code that could generate random non-overlapping cylinders with specified aspect ratio and specified number. Could you show how to add n? $\endgroup$
    – skonda2
    Aug 9, 2017 at 12:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.