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

Here is another approach you might find useful: https://mathematica.stackexchange.com/a/86277/9490 — Jason B., Mar 21 '18 at 14:54
Wait... so you only have the nodes, but no connectivity information? — J. M.'s missing motivation, Mar 21 '18 at 15:24
@JasonB. I had already tried the code from RunnyKine, but this as well seems to give me the area and not the volume. I belive RegionBoundary does not give out the volume. — JHN_28, Mar 21 '18 at 15:30
@J.M. Unfortunately I have only the node coordinates in the dat file. But, I could figure out a way to get the connectivity of QUAD elements on the surface. — JHN_28, Mar 21 '18 at 15:36
How accurate does it have to be? If not too much you might try a by-hand Monte Carlo approach using random points in a surrounding rectangular prism and taking ration(points inside)/(total points). — Daniel Lichtblau, Mar 21 '18 at 16:12
@DanielLichtblau Well up to at least second digit. The initial volume of my solid before deformation was 1 cubic-units, and I don't expect the volume to chage drasticically (< 2 cubic-units). I did look in to Monte Carlo methods, but decided againts it. There are no points inside, all the points I have are on the surface. — JHN_28, Mar 21 '18 at 16:19
If you have the boundary faces you can use methods from here. — Daniel Lichtblau, Mar 21 '18 at 16:26
@Hem_28 If you can provide the Abaqus input file (.inp) and the file with calculated displacements, then it should be easy to compute the volume? — Pinti, Mar 21 '18 at 17:14
@Pinti If I understood you correctly, you are asking me to consider all the elements/nodes used in Abaqus and then generate a mesh in Mathematica and then compute the volume? The problem is, there are lot of vacious regions with in the solid. This would not generate a mesh for entire region with in the solid and would need detecting the vacious space and filling it with elements. And I have to do this operation multiple times and this would be time consuming. Or did you mean something else? — JHN_28, Mar 21 '18 at 17:36
@DanielLichtblau This looks way to complicated for an amature like me :D Anyway I am still wondering how this would give me the Volume of the Solid. Probably I missed something there. I will take a closer look in to the code later. — JHN_28, Mar 21 '18 at 17:52
There must be a way to export the whole boundary mesh (vertex coordinates and index quads). Within Mathematica, we can select that single connected component that makes up the outer boundary. — Henrik Schumacher, Mar 21 '18 at 22:42
if you can extract a (triangulated) surface you can use this https://mathematica.stackexchange.com/a/26015/2079 (convexity not required) — george2079, Mar 21 '18 at 23:22
...and no there is not a straightforward way to get abaqus to give you a surface-only mesh for a solid model... that maybe should be a separate question, though i'm not sure this is the place for it. — george2079, Mar 21 '18 at 23:29
@george2079 A FEM tool that is not able to compute and export the boundary of the computation domain? Rather unlikely (or completely useless). — Henrik Schumacher, Mar 22 '18 at 08:32
@Hem_28 Even having the full solid mesh, including vertex positions and the integer 8-tuples of the hexahedra would suffice to get a grip on that. — Henrik Schumacher, Mar 22 '18 at 08:34
@george2079 So I had tried this numerical computation before, I forgot to mention that in my initial post. I have edited my post with the formula which I had used. It returned the volume as 0.0256571. I suppose the faces generated by my code are not consistently oriented. — JHN_28, Mar 22 '18 at 10:11
@HenrikSchumacher I will check your method. I need to reorganize by vertex (node) index, and consistently update the element connectivity. Currently my node numbers and element numbers are highly random, and not sequential. — JHN_28, Mar 22 '18 at 10:24
@Hem_28 The ordering of the vertices doesn't matter at all as long as you can also export the connectivity information of the element. I do not know how Abaqus calls this data. For Hexahedra, it should be a list of 8-tuples of integers. That has to be exported along with the vertex positions. — Henrik Schumacher, Mar 22 '18 at 10:29
@HenrikSchumacher yes it is a list of 8-tuples of integers. What I meant was, the node numbers on the surfaces are not continuous.(eq. 1001, 1110, 1200, 1350,...), so the element connectivity I read in would correspond to these numbers. As of now I could not find a method where I could add the index for a certain coordinate. If I am able to directly read the coordinates with it's index number (node number), then the element connectivity I have would be valid. I'll dig around a little to see if I can assign index to coordinates. — JHN_28, Mar 22 '18 at 11:01
@Hem_28 Really: No need for reordering. I'll do that for you. — Henrik Schumacher, Mar 22 '18 at 16:04
@Ah If I could by pass the reordering, it would save a great ddeal of work. I have uploaded a file on my actual post. Would be great if there is an alternate method than my cumbersome solution :) — JHN_28, Mar 22 '18 at 16:36
Just in case of... Do you know ? you can extract le volume of each element in Abaqus through command: *ELEMENT OUTPUT EVOL At desired step and frame, you can sum all current element volume to obtain total volume. NicoLM — NicoLM, Nov 03 '20 at 22:15
The FEMAddOns Package now has the ImportMesh package absorbed. Install with ResourceFunction["FEMAddOnsInstall"][]. — user21, Nov 04 '20 at 07:03

Pinti · Answer 1 · 2018-03-22T19:07:44.230

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 *)

Hey I am working on it now. I used your package quickly, but I get a invalid syntax and this error ImportMesh`Private`Abaqus`UpTo[4] is not a valid Span specification. A Span specification should be 1, 2, or 3 integers separated by ;;. (Any of the integers can be omitted or replaced with All.). I have edited my input deck manually, this could be due to that. But, I am comparing it to your example input file and checking if I could rectify the error. — JHN_28, Mar 22 '18 at 12:28
Sorry, I have used a new Span syntax, updated in MMA version 11.2. I will try to figure out a workaround for slightly older versions. — Pinti, Mar 22 '18 at 12:34
Development version on Github is now compatible with MMA 11 and higher. — Pinti, Mar 22 '18 at 13:02
Worth noting you can also just load the package from the web like: Get["https://raw.githubusercontent.com/c3m-labs/ImportMesh/master/ImportMesh.wl"] — b3m2a1, Mar 26 '18 at 08:24
The FEMAddOns Package now has the ImportMesh package absorbed. Install with ResourceFunction["FEMAddOnsInstall"][]. — user21, Nov 04 '20 at 07:03

Henrik Schumacher · Answer 2 · 2018-03-22T18:31:54.550

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.

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

2 Answers2