7
$\begingroup$

The moment of inertia for Regions seems to be wrongly calculated. Simple example:

PolyhedronData["Octahedron", "InertiaTensor"]*PolyhedronData["Octahedron", "Volume"] // MatrixForm

yields correctly (see also wikipedia)

{{1/(15 Sqrt[2]), 0, 0}, {0, 1/(15 Sqrt[2]), 0}, {0, 0, 1/(15 Sqrt[2])}}

(* N[,7] *)

{{0.04714045, 0, 0}, {0, 0.04714045, 0}, {0, 0, 0.04714045}}

However the function

MomentOfInertia[PolyhedronData["Octahedron", "Region"]]

yields incorrectly

{{0.0465585, 0., 0.}, {0., 0.0465585, 0.}, {0., 0., 0.0465585}}

I don't find a way to improve the accuracy of the Region function. The same accuracy problem persists for more complex shapes loaded from files. This makes the computation of the moment of inertia with Mathematica (tested with version 11.3.0.0) unreliable. Similar problem persist for

MomentOfInertia[PolyhedronData["Tetrahedron", "Region"]]

and

MomentOfInertia[PolyhedronData["Cube", "Region"]]

quits the kernel. Did anyone observe a similar behavior or knows how to increase the accuracy of Region computations?

$\endgroup$
  • $\begingroup$ Welcome to Mma.SE. Thanks for taking the tour. Be sure you have learned about asking and what's on-topic. Always edit if improvable. Showing due diligence, giving brief context, including minimal working examples of code and data in formatted form and so on help us to help you and in most cases inspires great answers. The site depends on participation, as you receive give back: vote and answer questions, keep the site useful, be kind, correct mistakes and share what you have learned. $\endgroup$ – rhermans Aug 21 '18 at 17:15
  • 1
    $\begingroup$ Try using the"ImplicitRegion" property rather than a mesh based one (which is where I believe the approximations are coming from). $\endgroup$ – chuy Aug 21 '18 at 19:42
  • $\begingroup$ @Tobias That (in particular the kernel quits) looks like a severe issue! You should report that to Wolfram Support. $\endgroup$ – Henrik Schumacher Aug 22 '18 at 7:56
9
$\begingroup$

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 MeshRegions and BoundaryMeshRegions 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...

$\endgroup$
  • $\begingroup$ Thanks, this works with fixed floating point precision, also for loaded triangle meshes. There is a small typo in the function MyMomentOfInertia[ B_BoundaryMeshRegion...] you should replace R with B in MeshCoordinates[R] and MeshCells[R,...] $\endgroup$ – Tobias Aug 22 '18 at 18:52
  • $\begingroup$ You're welcome. Also thanks for the hint to the typo! Fixed it. $\endgroup$ – Henrik Schumacher Aug 22 '18 at 19:23
  • 1
    $\begingroup$ Since the input region has inexact coordinates I think it does make sense to return an inexact answer. $\endgroup$ – Chip Hurst Aug 29 '18 at 17:21
2
$\begingroup$

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}}
$\endgroup$

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.