Skip to main content

graphics3d - 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[[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.



Answer



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.


Comments

Post a Comment

Popular posts from this blog

plotting - Filling between two spheres in SphericalPlot3D

Manipulate[ SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, Mesh -> None, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], {n, 0, 1}] I cant' seem to be able to make a filling between two spheres. I've already tried the obvious Filling -> {1 -> {2}} but Mathematica doesn't seem to like that option. Is there any easy way around this or ... Answer There is no built-in filling in SphericalPlot3D . One option is to use ParametricPlot3D to draw the surfaces between the two shells: Manipulate[ Show[SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], ParametricPlot3D[{ r {Sin[t] Cos[1.5 Pi], Sin[t] Sin[1.5 Pi], Cos[t]}, r {Sin[t] Cos[0 Pi], Sin[t] Sin[0 Pi], Cos[t]}}, {r, 1, 2 - n}, {t, 0, Pi}, PlotStyle -> Yellow, Mesh -> {2, 15}]], {n, 0, 1}]
Post a Comment
Read more

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]
Post a Comment
Read more

plotting - Adding a thick curve to a regionplot

Suppose we have the following simple RegionPlot: f[x_] := 1 - x^2 g[x_] := 1 - 0.5 x^2 RegionPlot[{y < f[x], f[x] < y < g[x], y > g[x]}, {x, 0, 2}, {y, 0, 2}] Now I'm trying to change the curve defined by $y=g[x]$ into a thick black curve, while leaving all other boundaries in the plot unchanged. I've tried adding the region $y=g[x]$ and playing with the plotstyle, which didn't work, and I've tried BoundaryStyle, which changed all the boundaries in the plot. Now I'm kinda out of ideas... Any help would be appreciated! Answer With f[x_] := 1 - x^2 g[x_] := 1 - 0.5 x^2 You can use Epilog to add the thick line: RegionPlot[{y < f[x], f[x] < y < g[x], y > g[x]}, {x, 0, 2}, {y, 0, 2}, PlotPoints -> 50, Epilog -> (Plot[g[x], {x, 0, 2}, PlotStyle -> {Black, Thick}][[1]]), PlotStyle -> {Directive[Yellow, Opacity[0.4]], Directive[Pink, Opacity[0.4]],
Post a Comment
Read more