Volumes and moments of inertia of meshes in 3D

Question

I came across puzzling results for certain 3D meshes, which I reduced to a minimal case. I made a 1-2-3 cuboid in Blender, exported it to this OBJ file, then imported it to Mathematica, which generates a plausible result:

(rawBlenderCuboid = 
 Import["https://filebin.ca/3TcNMJskChjW/one_two_three_cuboid.obj"]) // InputForm

MeshRegion[{{1., -3., -2.}, {1., -3., 2.}, {-1., -3., 2.}, {-1., -3., 
 -2.0000009536743164}, {1., 3., -1.9999990463256836}, {-1., 3., -2.}, 
 {-1., 3., 1.9999990463256836}, {0.9999989867210388, 3., 
 2.0000009536743164}}, {Polygon[{{1, 3, 4}, {5, 7, 8}, {1, 8, 2}, {2, 7, 
   3}, {3, 6, 4}, {5, 4, 6}, {3, 1, 2}, {7, 5, 6}, {8, 1, 5}, {7, 2, 8}, 
   {6, 3, 7}, {4, 5, 1}}]}]

This has no defects according to FindMeshDefects. It has the expected area of $88=2*(2*4+4*6+6*2)$:

Area[rawBlenderCuboid]

Its volume, however, is zero

Volume[rawBlenderCuboid]

0

it should be $48$, as is found for the bounding region:

Volume[BoundingRegion[rawBlenderCuboid]]

48.

Its moment of inertia is inexplicable

MomentOfInertia[rawBlenderCuboid]

{{541.333, 7.56979*10^-6, -0.0000177622}, 
 {7.56979*10^-6, 242.667, -0.0000184774}, 
 {-0.0000177622, -0.0000184774, 421.333}}

A hand calculation with unit density leads one to expect $(208, 80, 160)$ for the principal components of moment of inertia, which is indeed found for the bounding region:

MomentOfInertia[BoundingRegion[rawBlenderCuboid]]

{{208., -7.10543*10^-15, -5.32907*10^-15}, 
 {-7.10543*10^-15, 80.0001, 1.15463*10^-14}, 
 {-5.32907*10^-15, 1.15463*10^-14, 160.}}

Bounding regions seem always to be cuboids, so they are not acceptable for my real examples, which are not convex. Therefore, the rather nice solutions from this SE post about convex hulls would not serve.

A clue comes from the fact that Mathematica does not find the mesh to be solid:

SolidRegionQ[rawBlenderCuboid]

False

whereas the bounding region is reckoned as solid:

SolidRegionQ[BoundingRegion[rawBlenderCuboid]]

True

Both the bounding region and the raw region are both valid regions according to MeshRegionQ, and they're both embedded in 3 dimensions:

RegionEmbeddingDimension[rawBlenderCuboid]

3

RegionEmbeddingDimension[BoundingRegion[rawBlenderCuboid]]

3

so it's not easy to see the more profound differences between them.

These examples from Wolfram's curated collections also have zero volume and are not solid.

I realize from this SE post that I can fit tetrahedra and get to a volume that way, or that I could implement integrals inside polyhedra using a source like this one. Either approach is considerable work, and I wonder whether I am just missing some really simple answer already done in Mathematica. I also think that the approaches documented here for parametric and procedurally generated meshes would not be of direct help. This is my first foray into the new mesh functions, so I am a long way from mastery.

Answer 1

When Mathematica gives correct surface area and volume of zero, what you have is 2-manifold embedded in 3D space, not a solid. It is like the difference between a Sphere and Ball.

Area[Sphere[]]

4 π

Volume[Sphere[]]

0

Area[Ball[]]

∞

Volume[Ball[]]

(4 π)/3

Area[RegionBoundary[Ball[]]]

4 π

Update

box = 
  MeshRegion[
    {{1., -3., -2.}, {1., -3., 2.}, {-1., -3., 2.}, {-1., -3., -2.}, 
     {1., 3., -2.}, {-1., 3., -2.}, {-1., 3., 2.}, {1., 3., 2.}},
    {Polygon[
      {{1, 3, 4}, {5, 7, 8}, {1, 8, 2}, {2, 7, 3}, {3, 6, 4}, {5, 4, 6}, 
       {3, 1, 2}, {7, 5, 6}, {8, 1, 5}, {7, 2, 8}, {6, 3, 7}, {4, 5, 1}}]}]

 RegionDimension @ box

2

This confirms you have a hollow box; a solid brick would give 3.

2nd Update

You can generate a brick from a box with DelaunayMesh.

brick = DelaunayMesh @ MeshCoordinates[box]

RegionDimension @ brick

3

Volume[brick]

48

MomentOfInertia[brick] // Chop

{{208., 0, 0}, {0, 80., 0}, {0, 0, 160.}}

Answer 2

If the mesh is a triangle mesh, the "Stokes theorem" can be applied in the following way:

M = Import["https://filebin.ca/3TcNMJskChjW/one_two_three_cuboid.obj"];
With[{triangles = Partition[MeshCoordinates[M][[Flatten[MeshCells[M, 2][[All, 1]]]]], 3]},
 Total[Det /@ triangles]/6.
]
(* 48. *)

But the easiest way will be to transform your mesh into a BoundaryMeshRegion.

B = BoundaryMeshRegion[MeshCoordinates[M], MeshCells[M, 2]];
Volume[B]
MomentOfInertia[B]
(* 48. *)
(* {{208., 3.74847*10^-6, -3.57628*10^-6}, {3.74847*10^-6,80., -6.49028*10^-6}, {-3.57628*10^-6, -6.49028*10^-6, 160.}} *)

Answer 3

You can import OBJ files directly as BoundaryMeshRegion's

In[5]:= RegionDimension@Import["https://filebin.ca/3TcNMJskChjW/one_two_three_cuboid.obj",
  "BoundaryMeshRegion"]
Out[5]= 3

This is documented on ref/format/OBJ.