MomentOfInertia for Polyhedra incorrectly computed?
Some ranting about the built-in implementation
Given that the inertia tensor is defined by integrals of quadratic polynomials which can be easily evaluated exactly on each simplex of a MeshRegion
, it is absolutely inacceptable that Mathematica uses inappropriate approximation formulae for the integrals.
As a temporary workaround: PolyhedronData["Octahedron", "Region"]
is a BoundaryMeshRegion
. One can discretize it into many small tetrahedra and compute MomentOfInertia
of that tensor.
Itrue = PolyhedronData["Octahedron", "InertiaTensor"]* PolyhedronData["Octahedron", "Volume"];
R = PolyhedronData["Octahedron", "Region"];
S = DiscretizeRegion[R, MaxCellMeasure -> 0.00001];
Max[Abs[MomentOfInertia[R, RegionCentroid[R]] - Itrue]]
Max[Abs[MomentOfInertia[S, RegionCentroid[S]] - Itrue]]
0.000581981
9.66155*10^-9
So, the error is now much smaller, but it is still there (an the computations were done with a ridiculously large amount of effort).
An alternative implementation
So, let's be constructive. A first step towards a more correct implementation could be the following. It uses a 4-point Gauss quadrature on tetrahedra (which is exact on polynomials up to order 2), so it will only work for simplicial regions. Other regions have to be discretized first. The method relies also on the tetrahedra and triangles of MeshRegion
s and BoundaryMeshRegion
s being consistently oriented (so this may be a potential source of bugs).
Quiet[Block[{PP, P, QQ, Q, r, X, A, quadraturepoints, quadratureweights, vol, integrand},
PP = Table[Compile`GetElement[P, i, j], {i, 1, 4}, {j, 1, 3}];
QQ = Table[Compile`GetElement[Q, j], {j, 1, 3}];
r = {X[[1]], X[[2]], X[[3]]} - QQ;
quadraturepoints = Table[Mean[PP] (1 - 1/Sqrt[5]) + PP[[i]]/Sqrt[5], {i, 1, 4}];
vol = 1/3! Det[Table[PP[[i]] - PP[[1]], {i, 2, 4}]];
quadratureweights = vol ConstantArray[1/4, 4];
integrand = X \[Function] Evaluate[r.r IdentityMatrix[3] - TensorProduct[r, r]];
getMomentOfInertia =
With[{code = N[quadratureweights.integrand /@ quadraturepoints]},
Compile[{{P, _Real, 2}, {Q, _Real, 1}},
code,
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
]
];
getMomentOfInertiaB =
With[{code = -N[
quadratureweights.integrand /@ quadraturepoints /. {
Compile`GetElement[P, 4, 1] -> 0,
Compile`GetElement[P, 4, 2] -> 0,
Compile`GetElement[P, 4, 3] -> 0
}
]},
Compile[{{P, _Real, 2}, {Q, _Real, 1}},
code,
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
]
]
]
];
MyMomentOfInertia[R_MeshRegion?Region`Mesh`Utilities`SimplexMeshQ,
p_] := Total@getMomentOfInertia[
Partition[
MeshCoordinates[R][[
Flatten[MeshCells[R, 3, "Multicells" -> True][[1, 1]]]]], 4],
p
];
MyMomentOfInertia[R_MeshRegion, p_] :=
MyMomentOfInertia[
DiscretizeRegion[R, MaxCellMeasure -> ∞], p];
MyMomentOfInertia[R_MeshRegion] :=
MyMomentOfInertia[R, RegionCentroid[R]];
MyMomentOfInertia[
B_BoundaryMeshRegion?Region`Mesh`Utilities`SimplexMeshQ, p_] :=
Total@getMomentOfInertiaB[
Partition[
MeshCoordinates[B][[
Flatten[MeshCells[B, 2, "Multicells" -> True][[1, 1]]]]], 3],
p
];
MyMomentOfInertia[B_BoundaryMeshRegion, p_] :=
MyMomentOfInertia[
DiscretizeRegion[B, MaxCellMeasure -> ∞], p];
MyMomentOfInertia[B_MeshRegion] :=
MyMomentOfInertia[B, RegionCentroid[B]];
Here a usage example:
Itrue = PolyhedronData["Octahedron", "InertiaTensor"] PolyhedronData["Octahedron", "Volume"];
R = PolyhedronData["Octahedron", "Region"];
S = DiscretizeRegion[R, MaxCellMeasure -> ∞];
Max[Abs[MomentOfInertia[R] - Itrue]] // RepeatedTiming
Max[Abs[MomentOfInertia[S] - Itrue]] // RepeatedTiming
Max[Abs[MyMomentOfInertia[R] - Itrue]] // RepeatedTiming
Max[Abs[MyMomentOfInertia[S] - Itrue]] // RepeatedTiming
{0.00164, 0.000581981}
{0.0017, 0.000581981}
{0.0000736, 6.93889*10^-18}
{0.0000736, 2.08167*10^-17}
It is not only more accurate but also faster than the current implementation.
Edit
Support just told me that the development team is informed. Let's await what happens in the upcoming versions...
This appears to be more than just floating point error. The exact computation yields the same result:
Itrue = PolyhedronData["Octahedron", "InertiaTensor"] * PolyhedronData["Octahedron", "Volume"];
ebmr = BoundaryMeshRegion[##, WorkingPrecision -> ∞]& @@ PolyhedronData["Octahedron", "GraphicsComplex"];
Iexact = MomentOfInertia[ebmr]
{{8 Sqrt[2]/243, 0, 0}, {0, 8 Sqrt[2]/243, 0}, {0, 0, 8 Sqrt[2]/243}}
Itrue - Iexact // Simplify
{{1/(1215 Sqrt[2]), 0, 0}, {0, 1/(1215 Sqrt[2]), 0}, {0, 0, 1/(1215 Sqrt[2])}}
N[%]
{{0.000581981, 0., 0.}, {0., 0.000581981, 0.}, {0., 0., 0.000581981}}
One workaround is to integrate manually:
M = {{y^2 + z^2, -x y, -y z}, {-y x, x^2 + z^2, -y z}, {-z x, -z y, x^2 + y^2}};
Integrate[M, {x, y, z} ∈ R]
{{0.0471405, 0., 0.}, {0., 0.0471405, 0.}, {0., 0., 0.0471405}}