Volume of HexahedronElement with non-coplanar faces

Question

Bug fixed in Version 12.0

---

I would like to calculate the volume of ElementMesh made of HexahedronElement. Even though "MeshOrder" is 1, elements can have "curved" faces (nodes of the same face are not co-planar). This is a MWE with one element.

Needs["NDSolve`FEM`"]

mesh = ToElementMesh[
  "Coordinates" -> {{0,0,0},{1,0,0},{1,1,0},{0,1,0},{0,0,0.5},{1,0,1},{1,1,1},{0,1,1}},
  "MeshElements" -> {HexahedronElement[{{1, 2, 3, 4, 5, 6, 7, 8}}]}
  ];
mesh["Wireframe"["MeshElementStyle" -> FaceForm[LightBlue],ImageSize -> 200]]

I have found 4 different methods and each gives me a different answer. The most straightforward is to use the "MeshElementMeasure" method of ElementMesh object.

Total@Flatten@mesh["MeshElementMeasure"]
(* 0.916667 *)

I converted the hexahedron to 5 tetrahedra and called the same method on them (function is defined bellow).

Total@Flatten[HexToTetrahedron[mesh]["MeshElementMeasure"]]
(* 0.833333 *)

Then my own implementation of Gauss integration of Jacobian determinant over element (function is defined below).

MeshElementVolume[mesh["Coordinates"], HexahedronElement, 1]
(* 0.875 *)

And finally NIntegrate which also works on ElementMesh objects.

NIntegrate[1., {x, y, z} ∈ mesh]
(* 0.876271 *)

What is "the most" correct way to calculate this volume? I know there are some assumptions involved on how to treat the curved face, but surely there must some common way to do this?

---

Definitions of functions used above:

HexToTetrahedron::type="ElementMesh should contain only hexadedral elements.";

HexToTetrahedron[mesh_ElementMesh]:=Module[{
    nodes,origElms,tetConnect,restructure,newElms
    },
    origElms=mesh["MeshElements"];

    If[Head@First[origElms]=!=HexahedronElement,Message[HexToTetrahedron::type];Return[$Failed]];

    tetConnect={{4, 1, 2, 5},{7, 5, 2, 6},{4, 2, 3, 7},{4, 5, 2, 7},{4, 5, 7, 8}};
    restructure=Function[{hexNodes},Part[hexNodes,#]&/@tetConnect];

    newElms=TetrahedronElement[
        Flatten[restructure/@First@ElementIncidents[origElms],1]
    ];

    ToElementMesh[
        "Coordinates"->mesh["Coordinates"],
        "MeshElements"->{newElms}
    ]
]

---

(* Works for one element only. *)
MeshElementVolume[nodes_List,type_,meshOrder_]:=Block[{
    igCrds=ElementIntegrationPoints[type,meshOrder],
    igWgts=ElementIntegrationWeights[type,meshOrder],
    shapeDerivative=ElementShapeFunctionDerivative[type,meshOrder],
    jacobian,r,s,t
    },

    jacobian=Function[{r,s,t},Det[(shapeDerivative@@{r,s,t}).nodes]];

    (jacobian@@@igCrds).igWgts
]

Answer 1

The volume integral over the standard cube (to my surpise, it's Cuboid[{-1, -1, -1}, {1, 1, 1}]) can be computed analytically like this:

pp = Table[Compile`GetElement[p, i, j], {i, 1, 8}, {j, 1, 3}];
f = {r, s, t} \[Function] Evaluate[
    ElementShapeFunction[HexahedronElement, 1][r, s, t].pp
    ];
det = Det[D[f[r, s, t], {{r, s, t}, 1}]];
coeffrules = CoefficientRules[det, {r, s, t}];
vol = Values[coeffrules].Table[
    Integrate[
     FromCoefficientRules[{c -> 1}, {r, s, t}],
     {r, -1, 1}, {s, -1, 1}, {t, -1, 1}
     ], 
    {c, Keys[coeffrules]}];
cvol = With[{code = vol}, Compile[{{p, _Real, 2}},
    code,
    CompilationTarget -> "C",
    RuntimeAttributes -> {Listable},
    Parallelization -> True
    ]];

Now let's try it:

hexdata = Partition[
   mesh["Coordinates"][[Flatten[mesh["MeshElements"][[1, 1]]]]],
   8
   ];
cvol[hexdata]

{0.875}

Answer 2

Interestingly enough, if you convert the mesh to second order, the correct result comes up. So I made a mess in the 1st order symbolic code.

m2 = MeshOrderAlteration[mesh, 2];
m2["MeshElementMeasure"]
{{0.875`}}

Answer 3

This is fixed in version 12.0

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[
   "Coordinates" -> {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 
      0, 0.5}, {1, 0, 1}, {1, 1, 1}, {0, 1, 1}}, 
   "MeshElements" -> {HexahedronElement[{{1, 2, 3, 4, 5, 6, 7, 
        8}}]}];
Total@Flatten@mesh["MeshElementMeasure"]
0.875`

Answer 4

Please note that the shape of a curved 2-dimensional facet cannot be determined by vertex' positions only. In other words, given 4 points in 3D space, there are infinitely many surfaces enclosed by a "frame" (made by connecting the vertices with straight line segments). Thus the volume occupied by your 3D element is not well defined, as the element itself is not well defined.

That being said, I think different code samples you presented produce a correct result for some 3D element with the vertices specified.

From Grandma' Wikipedia: "The interior surface (or area) of such a polygon is not uniquely defined."

Enjoy your FEM!