6
$\begingroup$

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) *)
$\endgroup$
  • 1
    $\begingroup$ @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}] $\endgroup$ – flinty Jun 15 at 16:19
  • 1
    $\begingroup$ @flinty Ah! You're right. Deleting previous comment. $\endgroup$ – MarcoB Jun 15 at 16:21
  • $\begingroup$ Just so you know, this is still broken in v12.1.1 and continues to crash the kernel :( $\endgroup$ – flinty Jun 23 at 13:17
5
$\begingroup$

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]]]

enter image description here

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]

enter image description here

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}};
| improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Yizhen Chen Jun 16 at 12:22
  • $\begingroup$ @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. $\endgroup$ – user21 Jun 16 at 12:29
  • $\begingroup$ 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. $\endgroup$ – Yizhen Chen Jun 16 at 14:24
  • $\begingroup$ 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. $\endgroup$ – Yizhen Chen Jun 16 at 14:35
  • $\begingroup$ @YizhenChen What are the expected surface area and the expected volume for the octahedron? $\endgroup$ – user21 Jun 17 at 5:30
6
$\begingroup$

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]

hexagon blob plots

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]]

approx mesh blob hexagon

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.

| improve this answer | |
$\endgroup$
  • 3
    $\begingroup$ "What a mess." That are the right words to describe this! $\endgroup$ – Henrik Schumacher Jun 15 at 17:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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