# 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}] – flinty Jun 15 at 16:19
• @flinty Ah! You're right. Deleting previous comment. – MarcoB Jun 15 at 16:21
• Just so you know, this is still broken in v12.1.1 and continues to crash the kernel :( – flinty Jun 23 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. – Yizhen Chen Jun 16 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. – user21 Jun 16 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. – Yizhen Chen Jun 16 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. – Yizhen Chen Jun 16 at 14:35
• @YizhenChen What are the expected surface area and the expected volume for the octahedron? – user21 Jun 17 at 5:30

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! – Henrik Schumacher Jun 15 at 17:56