Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I would like to be able to visualize the empty space in a Graphics3D plot. For example, if I plot hexagonally close packed spheres using something like:

Show[Graphics3D[Sphere[#, Sqrt[2]/2] & /@ Select[Tuples[{-1, 0, 1}, 3],
 Mod[Total[#], 2] == 0 &], Boxed -> False]]

Mathematica graphics

Then how would I show the space between the spheres? That is, I would like the space occupied by the spheres to be empty and the space between them to be filled in.

share|improve this question
    
The empty space between them. I just found out it is often called an inverse opal structure. See, for example, https://duckduckgo.com/?q=inverse+opal+structure –  Gary L. Gray May 25 at 15:48
    
I appreciate how the answers might do what the OP wants, but the question seems to make no sense. The spheres as closed surfaces divide the space into regions. Using even/odd crossing one can identify two regions, one of which is the space outside all the spheres. Its visualization would look like a solid block with holes inside it, which might be shown with opacity < 1. What does "between" the spheres mean -- inside the convex hull? –  Michael E2 May 26 at 3:45
add comment

3 Answers 3

You may try this

eq = And @@ (Total[({x, y, z} - #)^2] > 1/2 & /@ 
    Select[Tuples[{-1, 0, 1}, 3], Mod[Total[#], 2] == 0 &])
RegionPlot3D[eq, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, Mesh -> None, 
 PlotPoints -> 150]

enter image description here

Notice there are small holes at points of contact between the spheres.

You can also "bound" by a sphere instead of a cuboid, with

eq = And @@ 
  Prepend[Total[({x, y, z} - #)^2] > 1/2 & /@ 
    Select[Tuples[{-1, 0, 1}, 3], Mod[Total[#], 2] == 0 &], 
   x^2 + y^2 + z^2 < 1]

enter image description here

share|improve this answer
    
This is very helpful. Thank you! I am going to play around with it a bit, but it does essentially what I need. –  Gary L. Gray May 25 at 16:59
add comment

Using some tricks (for speed and show-off :)

b = Normal@LatticeData["FaceCenteredCubic", "Basis"]
l1 = Flatten[Table[{i, j, k}.b, {i, -1, 1}, {j, -1, 1}, {k, -1, 1}], 2];
f = Nearest[l1];

RegionPlot3D[
 Length@f[{x, y, z}, {1, 1/N@Sqrt@2}] < 1,
 {x, -1, 1}, {y, -1, 1}, {z, -1, 1},
 Mesh -> None, PlotStyle -> Directive[Yellow, Opacity[0.5]], 
 PlotPoints -> 100]

Mathematica graphics

share|improve this answer
add comment

I thought a cute idea would be to take the rectangle formed by four points of contact between the spheres and blow it out like a bubble:

enter image description here

It seemed a simple idea, but the code to put together the surface of the space bounded by the spheres and the planes through the points of contact was a little longer than I expected.

Graphics3D @ space[]

Mathematica graphics

The edges of the "triangles" and "squares" with arcs for sides lie on planes. We can solve for these and sort the points (some of the edges will be reversed, and it's easier just to sort them).

Code dump

Because I wanted to do the animation, I began saving bits of data. The function space`setup sets up the all the data the first time space[] is called. There is an option PlotPoints that controls how finely the spherical segment is subdivided into polygons. The data for the various properties may accessed with calls of the form space[prop, PlotPoints -> plotPoints].

ClearAll[space, space`setup];
Options[space] = {PlotPoints -> 15};
space["Properties"] = {"fullsurface", "allpoints", "spherepoints", 
   "rectangle", "rectpoints", "allspheres"};

space[property_String: "fullsurface", opts : OptionsPattern[]] /; 
   MemberQ[space["Properties"], property] && 
    space`setup[OptionValue[space, opts, PlotPoints]] := 
  space[property, PlotPoints -> OptionValue[PlotPoints]];

space`setup[Automatic] := space`setup[PlotPoints /. Options[space]];
space`setup[plotPoints_] := 
  Module[{get, vertices, pairs, edges, rectpoints, p, q, r, s, spacexf},
   get[prop_] := space[prop, PlotPoints -> plotPoints];
   vertices = DeleteCases[
     Select[Tuples[{-1, 0, 1}, 3], Mod[Total[#], 2] == 0 &],
     {0, 0, 0}];
   pairs = Subsets[vertices, {2}];
   edges = Pick[pairs, #, Min[#]] &@ N[EuclideanDistance @@@ pairs];
   {p, q, r, s} = # - First@vertices & /@ 
     Mean /@ Select[edges, MemberQ[#, First@vertices] &][[{1, 2, 4, 3}]];
   rectpoints = 
     N @ Flatten[
          Table[(1 - u) p + u q + (1 - v) p + v s - p,
           {u, 0, 1, 1/(plotPoints - 1)}, {v, 0, 1, 1/(plotPoints - 1)}],
          1];
   space["spherepoints", PlotPoints -> plotPoints] = 
    First@vertices + # & /@ (Normalize /@ rectpoints/Sqrt[2]);
   space["rectpoints", PlotPoints -> plotPoints] = 
    First@vertices + # & /@ rectpoints;

   spacexf = Flatten[Table[
      RotationMatrix[r2, {0, 1, 0}].RotationMatrix[r1, {1, 0, 0}],
      {r1, 0, Pi, Pi/2}, {r2, 0, 3 Pi/2, Pi/2}],
     1];
   space["allpoints", PlotPoints -> plotPoints] = Flatten[
     Transpose[
      spacexf.Transpose@get["spherepoints"],
      {1, 3, 2}], 1];
   space["allnormals", PlotPoints -> plotPoints] = Flatten[
     Transpose[
      spacexf.Transpose[-(First@vertices - # & /@ (get["spherepoints"]))],
      {1, 3, 2}], 1];

   space["allspheres", PlotPoints -> plotPoints] = Polygon[Flatten[
      Table[plotPoints^2 k + {
         {i + 1 + plotPoints*(j - 1),  i + plotPoints*(j - 1),  i + plotPoints*j},
         {i + 1 + plotPoints*j,    i + 1 + plotPoints*(j - 1),  i + plotPoints*j}},
       {k, 0, 11}, {i, plotPoints - 1}, {j, plotPoints - 1}],
      3]];
   space["alltriangles", PlotPoints -> plotPoints] = Polygon[
     Map[
      Function[pat, 
       With[{proj = Orthogonalize[Prepend[IdentityMatrix[3], pat]][[2 ;; 3]]}, 
        SortBy[Flatten @ Position[get["allpoints"], _?(#.pat >= 2 &)], 
         ArcTan @@ (proj.get["allpoints"][[#]]) &]]
       ],
      Tuples[{-1, 1}, 3]
      ],
     VertexNormals -> None];
   space["allsquares", PlotPoints -> plotPoints] = Polygon[
     Map[
      Function[pat, 
       SortBy[Flatten @ Position[get["allpoints"], pat], 
        ArcTan @@ 
          get["allpoints"][[#, 
            Flatten @ Position[pat, Except[1. | -1.], Heads -> False]]] &]],
      (Join[#, -#] &@N@IdentityMatrix[3] /. 0. -> _)
      ],
     VertexNormals -> None];
   space["fullsurface", PlotPoints -> plotPoints] = 
    GraphicsComplex[get["allpoints"],
     {EdgeForm[],
      get["allspheres"],
      get["alltriangles"],
      get["allsquares"]
      },
     VertexNormals -> get["allnormals"]
     ];

   space`setup[plotPoints] = True  (* Unset to reset *)
   ];

For the animation:

vertices = DeleteCases[
   Select[Tuples[{-1, 0, 1}, 3], Mod[Total[#], 2] == 0 &],
   {0, 0, 0}];
With[{plotPoints = 15},
  anim = Table[
    Graphics3D[{EdgeForm[],
      {GraphicsComplex[(1 - t) space["rectpoints"] + t space["spherepoints"],
        space["allspheres"][[All, 1 ;; 2*(plotPoints-1)^2]]
        ],
       Red, Point[vertices]}
      }, ViewPoint -> {-3, 1, 0}/2],
    {t, 0, 1, 1/15}
    ]
  ];

Export["Example.gif", anim]
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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