No numerics hacks here; this really computes the volume symbolically. It is a bit tedious and demands some tricks which may appear more obvious in this answer than they would really be on the first try.
Don't expect this code to magically be useful for other such problems; it definitely has its fragile parts.
The idea of the code below is to split the Steinmetz(-like?) solid into easier pieces for volume calculation; specifically involving just one Cylinder and a handful of HalfSpaces each. As we later see, this is actually a necessity to work around orientation problems in Mathematica. Splitting could be visualised with an exploded view drawing:
In order to accomplish this, we have to find an accurate definition for each subregion. Cylinder-cylinder surface intersections on the surface of the solid (parts of ellipses around the origin with various orientations), and points connecting those different parts (points where different pairs of intersections meet) help. These ellipses naturally lie on planes, and define HalfSpaces which can be used to constrain a subregion, and graph cycles of points on the surface of the same Cylinder define which cylinder is used for each of them as a "cap".
CylindricalDecomposition with "Components" is used to guarantee that each individual curve expression is a connected component (there may be several such components on each ellipse), and thus can be used to figure out which points lie on which continuous curve.
(* Unit-length cylinder direction vectors. *)
directions = Append[
Table[RotationTransform[a, {0, 0, 1}][{1, 1, 1}/Sqrt[3]],
{a, 0, 3 Pi/2, Pi/2}],
{0, 0, 1}];
(* Used often. *)
simplify = FullSimplify[#, Element[x | y | z, Reals]] &;
(* For extra warnings if things don't go as planned on comparisons. *)
eqSimplify = (If[! BooleanQ[#],
Echo[#, "Didn't simplify to a Boolean:"]]; #) &@*FullSimplify;
(* Helper function which returns the two ellipses on
surface intersection of two cylinders passing through
the origin, based on their unit-length direction. *)
ClearAll[cylinderRings];
cylinderRings[v1_List, v2_List] :=
With[{transform =
Thread[
{x, y, z} ->
{Normalize@Cross[v1, v2],
Sqrt[(1 - v1 . v2)/2] Normalize[v1 + v2],
Sqrt[(1 + v1 . v2)/2] Normalize@Cross[Cross[v1, v2], v1 + v2]} .
{x, y, z}]},
ImplicitRegion[simplify[# /. transform], {x, y, z}] & /@
{x^2 + y^2 == 1 && z == 0, x^2 + z^2 == 1 && y == 0}]
(* All expressions for separate dimensional components of
intersections of two cylinders on the surface of the solid. *)
curves =
With[
{solid = RegionIntersection @@ (Cylinder[2 {-#, #}] & /@ directions)},
(* Find topologically separate components of each ring on the surface. *)
(CylindricalDecomposition[
RegionMember[RegionIntersection[solid, #], {x, y, z}],
{x, y, z}, "Components"] & /@
cylinderRings @@ #) & /@
Subsets[directions, {2}] //
simplify // Flatten //
(* Filter out dimensionless solutions. *)
Select[RegionDimension[ImplicitRegion[#, {x, y, z}]] != 0 &]];
(* All points which lie on two separate curves. *)
points = SolveValues[#, {x, y, z}] & /@ Subsets[curves, {2}] //
FullSimplify // Flatten[#, 1] & // DeleteDuplicates;
(* The surface graph; vertices are points above, edges the curves. *)
graph =
(* Create a list of points on each curve. *)
Select[points,
eqSimplify@*RegionMember[ImplicitRegion[#, {x, y, z}]]] & /@
curves //
(* Find the shortest geometric ordering of these points as a
line per each curve. This is a bit of a hack... *)
Map[
First@TakeSmallestBy[Permutations[#],
RegionMeasure[Line[#], 1] &, 1] &] //
(* Create graph edges on basis of these line segments. *)
Map[UndirectedEdge @@@ Partition[#, 2, 1] &] // Flatten // Graph
This graph corresponds to the following (observe the matching number of edges connected to different vertices):
With this graph we can find out subgraphs of points for each Cylinder:
(* Compute per-cylinder subgraphs on the solid. *)
subgraphs =
Subgraph[graph,
(* Select subgraphs with vertices on each cylinder surface. *)
Select[VertexList[graph],
eqSimplify@*RegionMember[
(* Hack around deficiency in RegionBoundary;
implicit cylinder surface (boundary) regions. *)
RotationTransform[{{0, 0, 1}, #}]@
ImplicitRegion[x^2 + y^2 == 1, {x, y, z}]]]] & /@
directions
Now each subregion can be obtained with (region) intersections of aforementioned primitives. This code has a trick up its sleeve: technically one should care about handedness of cycles which form the definition of HalfSpaces constraining the region, but since these Cylinders are in symmetric regarding the origin it's not really necessary - wrong handedness results a symmetric mirror image of the subregion and has the same volume.
In order to pacify Volume computation we reorient the subregion so that the Cylinder is always oriented towards the $z$ axis by rotating HalfSpaces; otherwise Mathematica seems to have trouble succeeding in this task.
(* Compute the solid volume by summing up subregion volumes for
intersection of half-spaces and one cylinder each. *)
Parallelize@MapThread[
Function[{dir, subgraph},
Volume[
RegionIntersection[
(* Important hack: avoid arbitrarily oriented cylinders.
Use cylinder oriented towards the z axis instead
and rotate half-spaces to match. Without this,
Mathematica stumbles and fails to compute exact volumes;
this is probably an interaction between internal
CylindricalDecomposition sub-region result and
Integrate over regions. *)
Append[
(* Half-spaces' normal vectors are computed from
vertices as spanning vectors on each graph edge. *)
RotationTransform[{dir, {0, 0, 1}}]@
HalfSpace[Cross @@ #, {0, 0, 0}] & /@ #,
Cylinder[{{0, 0, -2}, {0, 0, 2}}]]]] & /@
(* The volume computation is performed for each 4-cycle on
their corresponding cylinder surfaces. *)
FindCycle[subgraph, {4}, All]],
{directions, subgraphs}] //
Flatten // Total // FullSimplify
$$-\frac{4}{3} \left(-2-10 \sqrt{2}+\sqrt{6 \left(19+6 \sqrt{2}\right)}\right)$$
or
ResourceFunction["RadicalDenest"][%]
$$-\frac{4}{3} \left(-2-10 \sqrt{2}+6 \sqrt{3}+\sqrt{6}\right)$$
The volume computation takes a while.
This matches numerically acquired results:
N[%, 50]
(* 4.4004547140460115048732334911411402560985863336674 *)
The above volume computation can be sped up in this case (but not necessarily even on rather similar cases) by implementing it as an integration in cylindrical coordinates over the region in the cylinder oriented on the $z$ axis:
(* Create rotated intersections of HalfSpaces like before. *)
MapThread[Function[{dir, subgraph},
RotationTransform[{dir, {0, 0, 1}}]@
RegionIntersection[
HalfSpace[Cross @@ #, {0, 0, 0}] & /@ #] & /@
FindCycle[subgraph, {4}, All]],
{directions, subgraphs}] //
Flatten //
ParallelMap[
(* Convert Cartesian coordinates to cylindrical before integration. *)
Activate@IntegrateChangeVariables[
Inactive[Integrate][
(* Integrate over the region inside the z-oriented cylinder. *)
Boole[RegionMember[#, {x, y, z}]],
Element[{x, y, z}, Cylinder[2 {{0, 0, -1}, {0, 0, 1}}, 1]]],
{r, \[Theta], zz}, "Cartesian" -> "Cylindrical"] &, #] & //
Total // FullSimplify
No matter what the coordinate system, I haven't had success on getting integration results from cylinders on their alignment as they are actually present in the solid. If that would be easily achievable all the work with piecing up the solid would be unnecessary, but it would appear to be less trivial an endeavour.
Code for Steinmetz solid visualisations featured above:
Show[
With[{solid =
RegionIntersection @@ (Cylinder[2 {-#, #}] & /@ directions)},
BoundaryDiscretizeRegion[solid, RegionBounds[solid],
MaxCellMeasure -> {1 -> 0.02}]],
Graphics3D[
{Thick, Black,
ScalingTransform[{1.001, 1.001, 1.001}] /@
Flatten@Parallelize[
MeshPrimitives[
DiscretizeRegion[ImplicitRegion[#, {x, y, z}],
Method -> "Semialgebraic"], {1}] & /@
curves],
Sphere[#, 0.03] & /@ points}]]
Note Method -> "Semialgebraic" - which is used to discretize implicit single-dimensional regions (that's lines and other curves) in 3D embedding. It's slow but reliable in this task unlike other methods.
MapThread[
Function[{dir, subgraph},
RegionIntersection[
BooleanRegion[
BooleanCountingFunction[{{0, Length@#}}, Length@#],
HalfSpace[Cross @@ #, {0, 0, 0}] & /@ #],
Cylinder[{-dir, dir}],
HalfSpace[-Total[VertexList@#], {0, 0, 0}]] & /@
FindCycle[subgraph, {4}, All]],
{directions, subgraphs}] // Flatten //
Map[BoundaryDiscretizeRegion[#, 1.3 {{-1, 1}, {-1, 1}, {-1, 1}},
MaxCellMeasure -> {1 -> 0.0075}] &] //
Map[TranslationTransform[
Normalize[RegionCentroid[#]]/2][#] &] // Show
Unlike the volume computing code this visualisation needs to care about handedness of subregion boundaries (clockwise/counter-clockwise). This is handled by BooleanRegion to first include both, and then accepting only the side where surface points actually lie with another half-space intersection. Since discretisation is just numerics the more convoluted implicit regions don't really matter much here.
This "exploded view" is important for visually inspecting that all pieces are actually accounted for in the volume computation. Recognising if complicated equations actually are exactly zero is a complicated problem, and can lead to Selects on basis of RegionMember results to drop matches on this indecisiveness. (This is also the reason for couple eqSimplify@*RegionMember constructs in the main code.)
Some additional comments on CylindricalDecomposition, more commonly know as cylindrical algebraic decomposition:
It has nothing to do with Cylinders in the sense of this question; it's just a coincidental name of a method. In general it's a very useful tool in real algebraic geometry, easing mechanised symbolic computation on semialgebraic sets (which correspond with what people often consider in "constructive solid geometry"), at least in low enough dimensions, polynomial orders and coefficient complexities.
Regarding my commentary on Mathematica failing to compute Volumes when Cylinders are not oriented "nicely", or in this case on a coordinate axis: CylindricalDecomposition is used, also to my knowledge internally, to split regions into subregions which are easier to handle by other functions. This splitting produces sub-region splits on coordinate axis directions. I tried to figure out what Mathematica finds hard in this case, and got hints that it fails to compute some volume integrals with some cylindrical decomposition generated subregion variations. How these regions are split depend on the orientation of the region, although it's only a difference in rotation...
