# Volume of HexahedronElement with non-coplanar faces

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["NDSolveFEM"]

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"];

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
]

• That looks a bit like a short coming in the hex mesh measure computation. Jun 14 '18 at 9:54
• It's split into tets but the assumtion is that the faces are planar. I'll need to think about a more general way to do it. This will take a bit of thought. MeshRegion does warn about this: MeshRegion[{{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}}, Hexahedron[{{1, 2, 3, 4, 5, 6, 7, 8}}]] Sorry about that. Jun 15 '18 at 12:20
• Would you mind if I change the title to '... with non coplanar faces' ? Dec 13 '18 at 13:00
• @Pinti I don't think that this would break any theoretical assumptions, at least if you employ tri-linear parameterization of quads (the approximation order would remain the same). It makes it just harder to discretize a given geometry. However, tri-quadratic parameterization should be able to approximate the boundary of the domain by an additional order and you would loose this additional order if you require the boundary faces of hexes to be planar. Dec 13 '18 at 14:01
• This is an example that would not work. You'd then need to use something like this :ToElementMesh[Cuboid[{-1, -1, -1}, {1, 1, 1}], MaxCellMeasure -> {"Volume" -> 0.005}, "MeshElementType" -> TetrahedronElement]; Dec 13 '18 at 14:07

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[CompileGetElement[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}

• I think this is essentially the same method as in my third example. Thanks for showing me a faster code, maybe I can make it work for all element types. Jun 14 '18 at 11:05
• You can find the base coordinates of the mother element with: Needs["NDSolveFEM"] order = 1; eleType = HexahedronElement; MeshElementBaseCoordinates[eleType, order] Dec 13 '18 at 9:38

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}}


This is fixed in version 12.0

Needs["NDSolveFEM"]
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
`

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."