# Computing volume of intersection of two regions

I am trying to compute the volume of intersection of the following two regions:

a = 0.857597;
b = 1.653926;
hexagon = Polygon[{{0, (b - a)/2, 1/2}, {(b - a)/2, 0, 1/2},
{1/2, 0, (b - 1)/(2 a)}, {1/2, (b - 1)/2, 0}, {(b - 1)/2, 1/2, 0},
{0, 1/2, (b - 1)/(2 a)}}];
octahedron = ImplicitRegion[Abs[x] + Abs[y] + a Abs[z] <= b/2, {x, y, z}];
region2 = ImplicitRegion[1 >= RegionDistance[hexagon, {x, y, z}], {x, y, z}];


NIntegrate directly doesn't work:

NIntegrate[1, {x, y, z} ∈ RegionIntersection[octahedron, region2]]


It results in a crash after using up the memory (32GB).

I tried to use DiscretizeRegion first:

octd = DiscretizeRegion[octahedron, {{-1, 1}, {-1, 1}, {-1, 1}}];
regd = DiscretizeRegion[region2, {{-1, 2}, {-1, 2}, {-1, 2}}]; (* This takes 40 minutes *)
RegionIntersection[octd, regd]


This returns an error: “BoundaryMeshRegion: The boundary surface is not closed because the edges <<2>> only come from a single face.”

I also tried to discretize the regions using NDSolveFEMToElementMesh.

Needs["NDSolveFEM"];
ToElementMesh[region2, {{-1, 2}, {-1, 2}, {-1, 2}}]


This crashes without using significant memory. Computing finite element mesh on the first region does not crash, but intersecting it with the second region results in a crash without significant memory usage.

octf = ToElementMesh[octahedron, {{-1, 1}, {-1, 1}, {-1, 1}}];
RegionIntersection[octf, regd]


I have reported the issues with ToElementMesh to Wolfram Support.

Is there any workaround?

\$Version (* 12.1.0 for Mac OS X x86 (64-bit) (March 18, 2020) *)

• @MarcoB the intersection is not with the hexagon, but it's with region2, a blob around the hexagon where all points in the region are within distance 1 from the hexagon. You can visualize it here: RegionPlot3D[ 1 >= RegionDistance[hexagon, {x, y, z}], {x, -2, 2}, {y, -2, 2}, {z, -2, 2}] Commented Jun 15, 2020 at 16:19
• @flinty Ah! You're right. Deleting previous comment. Commented Jun 15, 2020 at 16:21
• Just so you know, this is still broken in v12.1.1 and continues to crash the kernel :( Commented Jun 23, 2020 at 13:17

Here is an approach based on creating exact regions:

a = Rationalize[0.857597, 10^-16];
b = Rationalize[1.653926, 10^-16];
hexagon =
Polygon[{{0, (b - a)/2, 1/2}, {(b - a)/2, 0, 1/2}, {1/2,
0, (b - 1)/(2 a)}, {1/2, (b - 1)/2, 0}, {(b - 1)/2, 1/2, 0}, {0,
1/2, (b - 1)/(2 a)}}] // Simplify;
octahedron =
ImplicitRegion[Abs[x] + Abs[y] + a Abs[z] <= b/2, {x, y, z}];
rd = RegionDistance[hexagon, {x, y, z}];
region2 = ImplicitRegion[1 >= rd, {x, y, z}];
ri = RegionIntersection[octahedron, region2];


This will run for a few seconds but will return an exact region that we then can mesh.

Needs["NDSolveFEM"]
bounds = {{-1, 1}, {-1, 1}, {-1, 1}};
mesh = ToElementMesh[ri, bounds,
"BoundaryMeshGenerator" -> {"RegionPlot",
"SamplePoints" -> {15, 15, 31}}];
mesh["Wireframe"["MeshElementStyle" -> FaceForm[Green]]]


NIntegrate[1, {x, y, z} \[Element] mesh]
0.871456


I have also tried to make use of the OpenCasadeLink based on the approach given by @flinty.

hexcenter = RegionCentroid[hexagon];
hexnormal =
Normalize[
Cross[hexagon[[1, 1]] - hexcenter, hexagon[[1, 2]] - hexcenter]];
hexradius = Norm[hexcenter - hexagon[[1, 1]]];
cylinderhack =
Cylinder[{hexcenter - hexnormal, hexcenter + hexnormal},
hexhack =
Flatten[{MeshPrimitives[hexagon, 1] /. Line -> Cylinder,
MeshPrimitives[hexagon, 0] /. Point -> Ball, cylinderhack}];


Needs["OpenCascadeLink"]


If you have a better representation of the octahedron, then we'd not need to convert to a boundary element mesh that is then converted to open cascade.

Get the boundary element mesh:

bmesh2 = OpenCascadeShapeSurfaceMeshToBoundaryMesh[res];


However, when we look at the MeshRegion version of the boundary element mesh we will see that there is a very slight elevation at the intersection - it's very hard to see at the top left corner:

MeshRegion[bmesh2]


And that can not be meshed with ToElementMesh - which is not ideal but understandable.

Edit by @YizhenChen:

The following representation of the octahedron gives more accurate answers:

octahedron = ConvexHullMesh[{{b/2, 0, 0}, {-b/2, 0, 0}, {0, b/2, 0},
{0, -b/2, 0}, {0, 0, b/(2 a)}, {0, 0, -b/(2 a)}}];


The cylinderhack given by @flinty is also incorrect, because it results in the "very slight elevation" seen in the figure above. The correct one is:

cylinderhack =
Apply[Prism[{hexagon[[1, #1]] + hexnormal,
hexagon[[1, #2]] + hexnormal, hexagon[[1, #3]] + hexnormal,
hexagon[[1, #1]] - hexnormal, hexagon[[1, #2]] - hexnormal,
hexagon[[1, #3]] - hexnormal}] &, #] & /@ {{1, 2, 3},
{1, 3, 4}, {1, 4, 5}, {1, 5, 6}};

• Thank you. I can replicate the result. However, now I worry about the efficiency of computing the exact region, as in my actual project (not included in the question) I need to compute volume of intersection of more complicated regions. I will try similar techniques on my project (with lower precision), and I hope there could be a more numerical approach, in case the computation requires too much memory. Commented Jun 16, 2020 at 12:22
• @YizhenChen, I think you best bet is to try to make use of the open cascade link as much as possible, as that will create very good regions efficiently. Commented Jun 16, 2020 at 12:29
• Unfortunately open cascade link may not work. For example, the volume of the octahedron is NIntegrate[1, {x, y, z} \[Element] octahedron], which gives 0.879253. But NIntegrate[1, {x, y, z} \[Element] MeshRegion@OpenCascadeShapeSurfaceMeshToBoundaryMesh@OpenCascadeShape@ToBoundaryMesh@octahedron] gives the incorrect result 5.20195. Commented Jun 16, 2020 at 14:24
• The number 5.20195 is actually the surface area. To get the volume I tried NIntegrate[1, {x, y, z} \[Element] MeshRegion@ToElementMesh@OpenCascadeShapeSurfaceMeshToBoundaryMesh@OpenCascadeShape@ToElementMesh@octahedron]. This gives the incorrect result 0.874412. Commented Jun 16, 2020 at 14:35
• Now I have figured it out. MeshRegion@ToElementMesh@octahedron outputs something different from the octahedron; it is a smaller polyhedron with 12 faces, hence the volume and surface area are different. To get the octahedron I shouldn't use ImplicitRegion. The octahedron represented as the convex hull of the vertices gives the correct volume and surface area. I will edit your answer to include the correct representations. Commented Jun 17, 2020 at 7:55

This is not ideal, but it gives an approximate resulting region. I first generate random points on the hexagon and add a random vector on the unit sphere. I take the convex hull of the points which is acceptable because the blob must be convex. Finally I discretize the octahedron and intersect with crudehexagonblob:

crudehexagonblob =
ConvexHullMesh[# + RandomPoint[Sphere[#, 1]] & /@
RandomPoint[hexagon, 40000]];
RegionIntersection[DiscretizeRegion[octahedron], crudehexagonblob]


Sadly convex hull is buggy and if I do 50000 or 20000 points I get an empty region, so I did 40000 and it worked. What a mess.

You could find a way to represent region2 differently. I'm thinking you can put spheres at all vertices and cylinders along all edges and join it to a cylinder at the center. I think this combination of spheres and cylinders is identical to region2:

RegionPlot3D[1 >= RegionDistance[hexagon, {x, y, z}], {x, -2, 2}, {y, -2, 2}, {z, -2, 2}]

hexcenter = RegionCentroid[hexagon];
hexnormal = Normalize[Cross[hexagon[[1, 1]] - hexcenter, hexagon[[1, 2]] - hexcenter]];
hexradius = Norm[hexcenter - hexagon[[1, 1]]];
cylinderhack = Cylinder[{hexcenter - hexnormal, hexcenter + hexnormal}, hexradius];
hexhack = Flatten[{
MeshPrimitives[hexagon, 1] /. Line -> Cylinder,
MeshPrimitives[hexagon, 0] /. Point -> Ball,
cylinderhack}];
Graphics3D[hexhack]



Unfortunately I had to use the same hack with ConvexHullMesh and random points to get a mesh out of the RegionUnion of these combined cylinders and spheres, because if you discretize them individually and RegionUnion them together it fails. Still, this mesh is pretty good:

cvxhm = ConvexHullMesh[RandomPoint[RegionUnion[RegionBoundary /@ hexhack], 40000]]


And disappointingly we can't even intersect this with the octahedron! I welcome any advice to get this to work:

(* unfortunately this fails for me in v12.1 *)
RegionIntersection[
DiscretizeRegion@octahedron,
cvxhm
]
`

Even though it doesn't provide a satisfying answer, I hope I've provided something you or somebody else can build on.

• "What a mess." That are the right words to describe this! Commented Jun 15, 2020 at 17:56