Computing volume from node coordinates on the surface of a non-convex solid

Question

I have a question regarding computation of volume of an irregular shape. I have looked around Mathematica -StackExchange but I could not find a solution which suits me the best. Problem Description: I have the coordinates of the nodes on the surface of my body , which looks like this (image from Abaqus)

In the end I want to compute the Volume of the body. I strated of with ConvexHull and DelaunyMesh to compute my volume. But this was overestimating the volume and the hull lokked like this

After looking around (a lot) on stackexchange, I used the code from Simon Wood.

DataIn = AA[[All, {2, 3, 4}]] points = DeleteDuplicates[DataIn] pts, tetrahedra} = TetGenDelaunay[points] csr[{aa_, bb_, cc_, dd_}] := With[{a = aa - dd, b = bb - dd, c = cc - dd}, Norm[a.a Cross[b, c] + b.b Cross[c, a] + c.c Cross[a, b]]/(2 Norm[a.Cross[b, c]])] radii = csr[pts[[#]]] & /@ tetrahedra alphashape[rmax_] := Pick[tetrahedra, radii, r_ /; r < rmax] faces[tetras_] := Flatten[tetras /. {a_, b_, c_, d_} :> {{a, b, c}, {a, b, d}, {a, c, d}, {b, c, d}}, 1] externalfaces[faces_] := Cases[Tally[Sort /@ faces], {face_, 1} :> face] polys = externalfaces@faces@alphashape[0.5]; Gr = Graphics3D[GraphicsComplex[pts, Polygon@polys], Axes -> True, BoxRatios -> Automatic , Boxed -> True];

I played around a bit with the alphashape rmax value to get a decent image as below

Once I got this fit, I believed it was a cake walk to obtain volume, but I was too optimistic.

I used BoundaryDiscretizeGraphics@Show@DiscretizeGraphics[Gr] to generate a boundary discretization, from which I thought I could create a boundaryregion mesh and get the volume. The DiscretizeGraphics[Gr] seemed to work, but not BoundaryDiscretizeGraphics. It returned the following without an error -

I was able to use RegionMesh directly on DiscretizeGraphics[Gr], but it created 6 surfaces with 2D elements and returned TotalSurfaceArea of 6.3. I expect the volume to be around 1.2. It would be great if someone could help me out here. Unlike the image below, I have solids wich are not symmetric about any axis. All the surfaces have the Periodic Wave (Like the Top and Bottom Surface of this problem). Is there an effecient way to compute the volume of such non-convex solids.

I have uploaded the coordinates data file here and here.

Thank You!

Edit 1

I had also used a form of shoelace method to compute the volume (I took this code from one of the post here, I lost the reference unfortunately)

Disc = DiscretizeGraphics[Gr] Total@With[{coords = MeshCoordinates[Disc]}, MeshCells[Disc, 2] /. Polygon[{a_, b_, c_}] :-> coords[[a]].Cross[coords[[b]], coords[[c]]]]/6

This returned a value 0.0256571

Based on @george2079's suggestion for a similar method here , I obtained the same value but with a nagative sign. So I suppose I have to relook on how I compute the polys

Edit 2 / A cubersome solution

Just for information, the solid has a big Void on the inside like in the image below

I tried importing the entire (~close to 100000 elements) mesh in to mathematica and fill the void with additional mesh, but didn't succed. So I went back to ABAQUS, extracted the elements and nodes on the Surface, renumbered both elements and nodes to start with 1 , and then used this

mesh = ToBoundaryMesh[ToElementMesh["Coordinates" -> NodeCoords, "MeshElements" -> {HexahedronElement[ElementCon]}]]

This generated a Quad boundary mesh with 2 boundary surfaces

rd = MeshRegion[ ToElementMesh[mesh, "RegionHoles" -> None]] gave me a tet mesh and filled the inner region.

Volume@rd gave me a value of 1.09503, which seems quite reasonable. I need to repeat this for a cuboid of known volume and check if the method gives the right value. The tedious job for me here is to reorder/renumber my nodes and elements. If I could bypass this, it would save a lot of time for me.

The unordered nodes/elements are here and the ordered ones are here

Answer 1

This is a partial answer to OP's question that shows how to get quadrilaterals on the outer surface of solid mesh, assuming that we start from Abaqus input file (.inp). The answer uses functions from the package I am developing to help importing different mesh files from FEM software to Mathematica ElementMesh object. (I hope this kind of self-promotion is acceptable here?)

This is a small sample of Abaqus input/mesh to work with.

content = Import["https://pastebin.com/raw/q6uEwaYZ"];
file = Export[FileNameJoin[{NotebookDirectory[], "Abaqus_Hex8_mesh.inp"}],content, "Text"];

Assuming that we have ImportMesh package installed.

Get["ImportMesh`"]
mesh = ImportMesh[file]
(* ElementMesh[{{-5., 5.}, {-5., 5.}, {-5., 5.}}, {HexahedronElement["<" 27 ">"]}] *)

Then we add some displacement to node coordinates to simulate real calculated displacement and create mesh in deformed configuration.

EDIT: At first I added random displacements in each node, but it is easier to test with displacements that would come from homogeneous deformation of the domain. Then we can easily calculate volume of the mesh with other methods.

getDeformedMesh[mesh_ElementMesh, displacements_] := ToElementMesh[
  "Coordinates" -> mesh["Coordinates"] + displacements,
  "MeshElements" -> mesh["MeshElements"],
  (* To avoid error message, because in the original mesh file some 
  nodes are not used in element connectivity specification.*)
  "CheckIncidentsCompletness" -> False
  ]

displacements = {0.5 (#[[3]] + 5), 0, 0.5 (#[[3]] + 5)} & /@ mesh["Coordinates"];

defMesh = getDeformedMesh[mesh, displacements]
(* ElementMesh[{{-5., 10.}, {-5., 5.}, {-5., 10.}}, {HexahedronElement[ "<" 27 ">"]}] *)

defMesh["Wireframe"["MeshElement" -> "MeshElements",
  "MeshElementStyle" -> FaceForm[LightBlue], ImageSize -> 200, 
  Axes -> True]
]

ElementMesh with solid elements can be converted to ElementMesh with only "BoundaryElements". Nice feature here is that we get different element markers for each surface. QuadElement connectivity can be then used with @Henrik's answer.

bmesh = ToBoundaryMesh[defMesh]; 

quads = Join @@ ElementIncidents@bmesh["BoundaryElements"];
pts = bmesh["Coordinates"];

(* Henriks's function *)
Total[getSignedPyramidVolume[Partition[pts[[Flatten[quads]]], 4]]]
(* -1499.91 *)

EDIT: After Henrik has updated his function the value looks reasonable (except the sign is wrong). This is confirmed by alternative calculation of volume as the sum of volume of all solid elements given below. So it seems we have connected the whole workflow for OP's question?

Volume@MeshRegion[defMesh]
(* 1499.91 *)

Answer 2

Once you have stored the vertex positions of the enclosing quad mesh in the variable pts and extracted the quads of the outer boundary in consistent orientation and stored it in the variable quads, you can obtain the enclosed volume (up to sign) by summing the oriented volumes of the pyramids spanned by each quad and the origin. This can be computed with the help of the following function

Block[{PP, P, f},
     PP = Table[Compile`GetElement[P, i, j], {i, 1, 4}, {j, 1, 3}];
     f = {s, t, r} \[Function] 
       Evaluate[((PP[[1]] (1 - s) + s PP[[2]]) (1 - t) + 
           t (PP[[4]] (1 - s) + s PP[[3]])) r];

     getSignedPyramidVolume = With[{
        code = N[Integrate[
            Integrate[
             Integrate[
              Det[D[f[s, t, r], {{s, t, r}, 1}]],
              {r, 0, 1}],
             {s, 0, 1}],
            {t, 0, 1}] // Simplify]
        },
       Compile[{{P, _Real, 2}},
        code,
        CompilationTarget -> "C",
        RuntimeAttributes -> {Listable},
        Parallelization -> True,
        RuntimeOptions -> "Speed"
        ]
       ]
     ]

by

Total[getSignedPyramidVolume[Partition[pts[[Flatten[quads]]], 4]]]

Edit

Fixed ordering of the quad.