# NIntegrate of surface area of intersecting spheres yields zero

I have a bunch of spheres (it's actually diamond cubic structure. The 0.6 radius doesn't matter),

img = {Ball[{-1/2, -1/2, 1/2}, 0.6], Ball[{1/2, 1/2, 1/2}, 0.6],
Ball[{-1/2, 1/2, -1/2}, 0.6], Ball[{1/2, -1/2, -1/2}, 0.6],
Ball[{-1, -1, -1}, 0.6], Ball[{-1, -1, 1}, 0.6],
Ball[{-1, 1, -1}, 0.6], Ball[{-1, 1, 1}, 0.6],
Ball[{1, -1, -1}, 0.6], Ball[{1, -1, 1}, 0.6],
Ball[{1, 1, -1}, 0.6], Ball[{1, 1, 1}, 0.6], Ball[{0, 0, 1}, 0.6],
Ball[{1, 0, 0}, 0.6], Ball[{0, 1, 0}, 0.6], Ball[{-1, 0, 0}, 0.6],
Ball[{0, 0, -1}, 0.6], Ball[{0, -1, 0}, 0.6]};


I am interested in its volume and surface area for any radius inside a unit box. It looks like this (Show[Graphics3D[img], PlotRange -> {{-1, 1}, {-1, 1}, {-1, 1}}])

By doing

NIntegrate[
Boole[
Or@@((Norm[{x,y,z}-#]<=0.5) &/@ Cases[img,Ball[pos_,_]->pos])
], {x,-1,1},{y,-1,1},{z,-1,1}]


I can get the volume when $r=.5$.

Then I tried (<= changed to ==)

NIntegrate[
Boole[
Or@@((Norm[{x,y,z}-#]==0.5) &/@ Cases[img,Ball[pos_,_]->pos])
], {x,-1,1},{y,-1,1},{z,-1,1}]


in an attempt to get surface area. But it gives zero. Where was I wrong?

## Edit (2015-06-14)

I ran the methods of two of the answers below over $radius=\sqrt{2}/2$. I kept getting:

Unable to compute the area of region \
RegionIntersection[Cuboid[{-1,-1,-1},{1,1,1}],RegionDifference[\
RegionDifference[RegionDifference[RegionDifference[RegionDifference[\
RegionDifference[RegionDifference[RegionDifference[<<2>>],Ball[<<2>>]]\
,Ball[{<<3>>},Power[<<2>>]]],Ball[{1,0,0},1/Sqrt[2]]],Ball[{0,1,0},1/\
Sqrt[2]]],Ball[{-1,0,0},1/Sqrt[2]]],Ball[{0,0,-1},1/Sqrt[2]]],Ball[{0,\
-1,0},1/Sqrt[2]]]]. >>


They both work fine for 0.5. The code is simply as follows.

imgballs = {Ball[{-1/2, -1/2, 1/2}, 0.5], Ball[{1/2, 1/2, 1/2}, 0.5],
Ball[{-1/2, 1/2, -1/2}, 0.5], Ball[{1/2, -1/2, -1/2}, 0.5],
Ball[{-1, -1, -1}, 0.5], Ball[{-1, -1, 1}, 0.5],
Ball[{-1, 1, -1}, 0.5], Ball[{-1, 1, 1}, 0.5],
Ball[{1, -1, -1}, 0.5], Ball[{1, -1, 1}, 0.5],
Ball[{1, 1, -1}, 0.5], Ball[{1, 1, 1}, 0.5], Ball[{0, 0, 1}, 0.5],
Ball[{1, 0, 0}, 0.5], Ball[{0, 1, 0}, 0.5], Ball[{-1, 0, 0}, 0.5],
Ball[{0, 0, -1}, 0.5], Ball[{0, -1, 0}, 0.5]};

(* Michael Seifert's method *)
imgspheres = imgballs /. Ball[x___, 0.5] -> Sphere[x, Sqrt[2]/2];
imgballs2 = imgballs /. Ball[x___, y_] -> Ball[x, Sqrt[2]/2];
reglist = Table[RegionDifference[
RegionIntersection[imgspheres[[i]], Cuboid[{-1, -1, -1}, {1, 1, 1}]],
RegionUnion[Drop[imgballs2, {i}]]], {i, 1, Length[imgspheres]}];
Area /@ reglist
Total[%]

(* george2079's method *)
MapIndexed[(Area@
RegionIntersection[
Fold[RegionDifference, #, Delete[imgballs2, #2]],
Cuboid[{-1, -1, -1}, {1, 1, 1}]]) &, imgspheres]
Total@%


Any idea to fix this? I'm actually interested in this radius.

• The area calculation fails because you are depending on the practically zero probability of NIntegrate's initial sampling happening to exactly hit a surface. Try to formulate it as an integral over the surface of each sphere, with a Boole test if the surface element is in the box. – george2079 Jun 10 '15 at 16:55
• @george2079 Do integration for each sphere and add them up? How can I deal with the intersections then? – Richard Cox Jun 10 '15 at 17:28
• The intersections are just lines (or points), and so they have zero area. Imagine, for example, a finite surface in the xy-plane and one in the xz-plane, intersecting at a line segment along the x-axis. The combined surface area of the complex is just the sum of the individual planar surfaces. – Michael Seifert Jun 10 '15 at 17:48
• @MichaelSeifert Yes, the intersections are lines. What I meant is the parts that overlap. For instance, look at the $1/8$ sphere on the upper left corner. A large part of its area is inside another sphere. I don't wanna count that. – Richard Cox Jun 10 '15 at 18:40
• Your last integral is set up as the triple integral over a volume (cube) of a function that is zero almost everywhere. It is not the integral you seek and its value is actually zero. (To answer "Where was [I] wrong?") – Michael E2 Jun 10 '15 at 20:59

If you're running Mathematica 10 or above, there's a dead simple method using the new Area function:

newimg = img /. Ball[x___, 0.6] -> Sphere[x, 0.5]
Total[Map[Area[RegionIntersection[#, Cuboid[{-1, -1, -1}, {1, 1, 1}]]] &, newimg]]

(* 25.1327 *)


Note a few things here:

1. Mathematica considers Ball to be a 3D object; if we want its surface area we have to use Sphere instead.

2. Mathematica is much better at dealing with the individual spheres and their intersections with the cube than it is at dealing with the whole complex at once. In principle you could also use RegionUnion to define the entire surface as the union of these Spheres, but when I tried that approach Mathematica just gave up.

3. Similarly, you could try using ImplicitRegion, almost as you did above, to define the complex, and then feed that into Area. However, I tried this just now and Mathematica is still thinking about it. (EDIT: I aborted the evaluation after about an hour with no results.) The above code runs in just a few seconds on my machine.

EDIT: It was pointed out in the comments that the OP may have wanted to exclude any of the spheres that lie within another sphere. This can be done as well, though it's a little trickier. Basically, we need to take the RegionDifference between each Sphere and all of the other Balls:

imgballs = {Ball[{-1/2, -1/2, 1/2}, 0.5], Ball[{1/2, 1/2, 1/2}, 0.5],
Ball[{-1/2, 1/2, -1/2}, 0.5], Ball[{1/2, -1/2, -1/2}, 0.5],
Ball[{-1, -1, -1}, 0.5], Ball[{-1, -1, 1}, 0.5],
Ball[{-1, 1, -1}, 0.5], Ball[{-1, 1, 1}, 0.5],
Ball[{1, -1, -1}, 0.5], Ball[{1, -1, 1}, 0.5],
Ball[{1, 1, -1}, 0.5], Ball[{1, 1, 1}, 0.5], Ball[{0, 0, 1}, 0.5],
Ball[{1, 0, 0}, 0.5], Ball[{0, 1, 0}, 0.5], Ball[{-1, 0, 0}, 0.5],
Ball[{0, 0, -1}, 0.5], Ball[{0, -1, 0}, 0.5]};
imgspheres = imgballs /. Ball[x___, 0.5] -> Sphere[x, 0.5];
reglist = Table[RegionDifference[
RegionIntersection[imgspheres[[i]], Cuboid[{-1, -1, -1}, {1, 1, 1}]],
RegionUnion[Drop[imgballs, {i}]]], {i, 1, Length[imgspheres]}];
Area /@ reglist
Total[%]

(* {2.29981, 2.29981, 2.29981, 2.29981, 0.392699, 0.182252, 0.182252,
0.392699, 0.182252, 0.392699, 0.392699, 0.182252, 1.1499, 1.1499, \
1.1499, 1.1499, 1.1499, 1.1499} *)

(* 18.3984 *)


Just for fun, let's try to visualize this region. I had to do some tricky shenanigans here to keep DiscretizeRegion happy; specifically, I found that if I took the differences with all of the other spheres, Mathematica couldn't keep track of which differences made a difference (i.e., which other spheres actually intersected a given sphere.) So I defined the adjmatrix below to pick out the balls that intersected a given sphere. I then had to make sure that I was subtracting something from each sphere, hence the ugly If statement that subtracts a "dummy" ball if there is no "real" ball intersecting a given sphere.

adjmatrix = Outer[((4 \[Pi] 0.5^3/3) > Volume[RegionIntersection[#1, #2]] >
0) &, imgballs, imgballs];
reglist =
Table[RegionDifference[RegionIntersection[imgspheres[[i]], Cuboid[{-1, -1, -1}, {1, 1, 1}]],
If[Pick[imgballs, adjmatrix[[i]]] != {}, RegionUnion[Pick[imgballs, adjmatrix[[i]]]], Sphere[{5, 5, 5}, 0.1]]],
{i, 1, Length[imgspheres]}]
Show[DiscretizeRegion /@ reglist, ImageSize -> Large]


The area of this new reglist is the same as before, as one can easily check.

• Thank you for the wonderful solution! I actually tried ImplicitRegion before posting the question, and you are right it can't work out. – Richard Cox Jun 10 '15 at 21:26
• Doesn't this method also count the area of a sphere lies within another sphere, which @Fred said he doesn't want to count? – Rahul Jun 11 '15 at 19:53
• Hmm, good point. The original code he tried to use would have had that same problem as well, though, (it would have counted all points that were the proper distance from one of the centers, regardless of whether they were inside a different sphere), so I'm a little confused. – Michael Seifert Jun 11 '15 at 19:58
• I've now modified my code to address @Fred's desire to only count surfaces that lie outside any of the other spheres. – Michael Seifert Jun 11 '15 at 21:29
• Using the numerical value doesn't seem to help either; Mathematica can't calculate the areas of some of the spheres. It can calculate the areas of some of them, though, particularly the ones entirely contained in the unit cell. I would focus on calculating the areas of one of those; if you can do that, then you can multiply the result by 8 (the number of atoms per unit cell) to get the result you need. You'll need to let the spheres go outside the cuboid for this radius value, though. – Michael Seifert Jun 15 '15 at 14:07

The old fashion way:

 r = .5;
NIntegrate[ r^2 Sin[phi] Boole[And @@ ((-1 < # < 1 ) & /@
(#[[1]] +
CoordinateTransform[
"Spherical" -> "Cartesian", {r, theta, phi}]))],
{theta, 0, 2 Pi} , {phi, 0, Pi} ] & /@ img // Total


25.1327

Edit: this is excluding areas that are inside any other sphere:

 centers = img[[All, 1]];
NIntegrate[
p = #[[1]] +
CoordinateTransform["Spherical" -> "Cartesian", {r, theta, phi}];
r^2 Sin[phi] Boole[ Min[(Norm[p - #]) & /@
Complement[centers, {#[[1]]}] ] > r &&
(And @@ ((-1 < # < 1) & /@ p))],
{theta, 0, 2 Pi}, {phi, 0, Pi}] & /@ img // Total


It is running very slowly and throwing convergence warnings however..

...finally finished with result:

19.181

This figure is generated using EvaluationMonitor to capture all the sample points that are used to compute the area. (gives an idea of why it is so slow )

## edit

a variant on Michael's approach..

 MapIndexed[(Area@RegionIntersection[
Fold[ RegionDifference, #  , Delete[ imgballs, #2] ],
Cuboid[{-1, -1, -1}, {1, 1, 1}]]) &, imgspheres]
Total@%


18.3984

• You get the same result. But like @Rahul said above, doesn't this count the overlapping part? Sorry I am new to Mathematica and the code is kinda hard for me to grasp. – Richard Cox Jun 11 '15 at 20:21
• After doing some research, I start to appreciate your code. Thanks! – Richard Cox Jun 11 '15 at 23:52
• I was running into a new problem when computing radius equal to sqrt(2)/2. Could you please take a look at my update above? Thanks!! – Richard Cox Jun 15 '15 at 2:20
• Cant test it , but I think you need to replace Delete[ imgballs, #2] with something like Select[ imgballs , Function[{ball}, Norm[ball[[1]]-#[[1]]]<2 r]]. If that works change MapIndexed to Map as well. – george2079 Jun 15 '15 at 14:13
• make that  0 < Norm[ball[[1]]-#[[1]]]<2 r  (still untested, no v10 here..) – george2079 Jun 15 '15 at 16:15

Another version 10 approach (barring errors on my part) to address volume as well as surface area. The (limitations of region) discretizations render only approximations, cf george2079 better numeric integration and the neat approach of Michael Seifert...both of which I have upvoted.

cb = Cuboid[{-1, -1, -1}, {1, 1, 1}];
i = RegionUnion @@ (RegionIntersection[cb, #] & /@ img);
t = RegionUnion @@ img;
o = RegionUnion @@ (RegionDifference[#, cb] & /@ img);
in = DiscretizeRegion[i];
tr = DiscretizeRegion[t];
out = DiscretizeRegion[o];
a0 = Total[Area /@ MeshPrimitives[BoundaryMesh[tr], 2]];
a1 = Total[Area /@ MeshPrimitives[BoundaryMesh[in], 2]];
a2 = Total[Area /@ MeshPrimitives[BoundaryMesh[out], 2]];
sa=a1 - (a1 + a2 - a0)/2
vol=Volume[in]


The surface area (needed to subtract the area of the cut sides): 17.1081, and volume: 5.42381.

The following renderings take some time (but I present for fun...and perhaps to correct my misconceptions):

rp0 = RegionPlot3D[t, PlotPoints -> 60, Boxed -> False, Axes -> False,
Background -> Black];
rp1 = RegionPlot3D[i, PlotPoints -> 100, Boxed -> False,
Axes -> False, Background -> Black];
rp2 = RegionPlot3D[o, PlotPoints -> 40, Boxed -> False, Axes -> False,
Background -> Black];


then

TabView[{"Full object" -> rp0, "Restricted object" -> rp1,
"Complement" -> rp2}]


• So does the 18.3984 contains cut sides? I don't want these cut sides actually. I'm confused again. – Richard Cox Jun 12 '15 at 17:17
• You can easily and efficiently compute the area using Area@RegionBoundary@tr etc. – RunnyKine Jun 12 '15 at 22:01
• @RunnyKine thank you for educating me in this regard. I really appreciate it. – ubpdqn Jun 12 '15 at 22:04
• @Fred the count is the polygons that make up the mesh discretization. – ubpdqn Jun 12 '15 at 22:06