Computing volume of intersection of two regions

Question

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 NDSolve`FEM`ToElementMesh.

Needs["NDSolve`FEM`"];
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) *)

Answer 1

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["NDSolve`FEM`"]
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}, 
   hexradius];
hexhack = 
  Flatten[{MeshPrimitives[hexagon, 1] /. Line -> Cylinder, 
    MeshPrimitives[hexagon, 0] /. Point -> Ball, cylinderhack}];

Load the link and convert the primitives into open cascade shapes:

Needs["OpenCascadeLink`"]
shapes = OpenCascadeShape /@ hexhack;
union = OpenCascadeShapeUnion[shapes];
oocOcta = OpenCascadeShape[ToBoundaryMesh[octahedron]];
res = OpenCascadeShapeIntersection[union, oocOcta];

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}};

Answer 2

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.