Computing volume of intersection of two regions
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.
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}};