六面体单元的体积计算方法

有限元中最常见的八节点六面体单元，通常来说其形状是一个不规则的六面体，对于不规则的六面体，体积应该怎样计算呢？我摸索到了两个方法，一个可以从理论上收敛到精确解，另一个方法可以快速得到一个相对比较精确的解。

方法1：立方体映射法（不规则六面体与立方体之间的映射）

这种方法的主要思想是把一个不规则六面体映射到一个立方体当中，通过计算不规则六面体与立方体之间的体积比例（分解为三维空间中三个不同方向上的尺寸比例），来得到不规则六面体的体积。

图1：不规则六面体与立方体之间的映射

用OpenGL画出来的不规则六面体

如何实现映射？

所谓映射，那就是点和点之间的一对一映射，不规则六面体内部的任意一个点，都可以唯一地映射到立方体内部的一个点，反之立方体内部的任意一个点也唯一对应不规则六面体中的一个点。怎样实现点和点之间的映射呢，

1. 第一步可以从顶点出发，把六面体的8个顶点和立方体的8个点对应起来
2. 第二步，根据匹配好的8个顶点的位置坐标，可以对内部所有点的位置进行线性插值来实现内部点之间的一一对应。插值方式可以使用形函数(shape function)

图2：自然坐标和全局坐标的映射

形函数插值法表示的坐标映射(线性插值)：

某个点从标准立方体对应的自然坐标 $\vec{\xi }=(\xi ,\eta ,\zeta )$ 映射到不规则六面体对应的全局坐标 $\vec{x}=(x,y,z)$ 的映射关系是：
$$\begin{gathered} \vec{x}=\sum _{i=1}^8{N}_i(\xi ,\eta ,\zeta )\cdot \overrightarrow{{\xi }_i},where{N}_i(\xi ,\eta ,\zeta )=\frac{1}{8}(1+{\xi }_i\xi )(1+{\eta }_i\eta )(1+{\zeta }_i\zeta ), \\ vertexe{s}^{\prime }coordinates({\xi }_i,{\eta }_i,{\zeta }_i)=(-1,1,1),or(-1,-1,1), \\ or(-1,-1,-1),or(-1,1,-1), \\ or(1,1,1),or(1,-1,1), \\ or(1,-1,-1),or(1,1,-1) \end{gathered}$$
， N_{i} 是立方体内部点p对于顶点 i 的权重，如果内部点p与顶点 $i$ 之间距离越近，那么权重 $N_i$ 就越高。权重满足两个条件是：

1. 如果点p和某个顶点A重合，那么顶点A的权重就是1，而其它顶点的权重就是0；
2. 所有顶点的权重加起来等于1
   权重 $N_i$可以表示为点p的自然坐标 $(\xi ,\eta ,\zeta )$ 的函数，称为形函数 $N_i(\xi ,\eta ,\zeta )$ ，函数形式如上面的公式所示；$\overrightarrow{{\xi }_i}$ 是立方体各个顶点 $i$ 的自然坐标。自然坐标的取值范围是 [-1, 1], 而顶点/端点的取值则要么是 1 要么是 -1；

利用坐标映射的体积比例计算体积（雅可比矩阵法）

有了坐标映射，就可以知道不规则六面体相对于立方体在三个不同方向上的长度伸缩关系，从而得到不规则六面体的体积。标准立方体的自然坐标方向的微矢量 $\overrightarrow{d\xi },\overrightarrow{d\eta },\overrightarrow{d\zeta }$ 可以分别映射到全局坐标中得到新矢量 $\overrightarrow{dx},\overrightarrow{dy},\overrightarrow{dz}$ :
$$\begin{gathered} \overrightarrow{d\xi }=(d\xi ,0,0){⟶}^{mapto}\overrightarrow{dx}=(\frac{\partial x}{\partial \xi }d\xi ,\frac{\partial y}{\partial \xi }d\xi ,\frac{\partial z}{\partial \xi }d\xi ) \\ \overrightarrow{d\eta }=(0,d\eta ,0){⟶}^{mapto}\overrightarrow{dy}=(\frac{\partial x}{\partial \eta }d\eta ,\frac{\partial y}{\partial \eta }d\eta ,\frac{\partial z}{\partial \eta }d\eta ) \\ \overrightarrow{d\zeta }=(0,0,d\zeta ){⟶}^{mapto}\overrightarrow{dz}=(\frac{\partial x}{\partial \zeta }d\zeta ,\frac{\partial y}{\partial \zeta }d\zeta ,\frac{\partial z}{\partial \zeta }d\zeta ) \end{gathered}$$
把不规则六面体切分成由上面的微矢量 $\overrightarrow{dx},\overrightarrow{dy},\overrightarrow{dz}$ 构成的平行六面体微元，对这些微元进行积分可以得到体积：
$$volume={\int }_{\tilde{\Omega }}(\overrightarrow{dx}\times \overrightarrow{dy})\cdot \overrightarrow{dz}={\int }_{\Omega }\left|\begin{array}{ccc}\frac{\partial x}{\partial \xi }d\xi & \frac{\partial y}{\partial \xi }d\xi & \frac{\partial z}{\partial \xi }d\xi \\ \frac{\partial x}{\partial \eta }d\eta & \frac{\partial y}{\partial \eta }d\eta & \frac{\partial z}{\partial \eta }d\eta \\ \frac{\partial x}{\partial \zeta }d\zeta & \frac{\partial y}{\partial \zeta }d\zeta & \frac{\partial z}{\partial \zeta }d\zeta \end{array}\right|={\int }_{\Omega }\left|\begin{array}{ccc}\frac{\partial x}{\partial \xi }& \frac{\partial y}{\partial \xi }& \frac{\partial z}{\partial \xi }\\ \frac{\partial x}{\partial \eta }& \frac{\partial y}{\partial \eta }& \frac{\partial z}{\partial \eta }\\ \frac{\partial x}{\partial \zeta }& \frac{\partial y}{\partial \zeta }& \frac{\partial z}{\partial \zeta }\end{array}\right|d\xi d\eta d\zeta$$
即关于雅可比矩阵的行列式值的积分。
这个积分会在积分点的增加之下最终收敛到一个精确的结果。
这里做一个测试，随便取一个六面体单元，其八个顶点的坐标分别是：

图3：设置六面体的顶点坐标(随便设置)

由于这里从规则六面体到不规则六面体的mapping使用了trilinear 三线性插值，trilinear多项式关于自然坐标的最高次数达到三次，因此使用二阶高斯积分可以达到积分精度（二阶高斯积分最高可达到3阶多项式的精度）。

方法1的代码

用来计算体积的代码如下所示，其中雅可比矩阵涉及求导的运算，因此这里使用pyTorch的自动微分来得到雅可比矩阵，然后计算行列式并且积分（二阶高斯积分）：

import torch as tch

def eleVol_gaussInteg(nodesXYZ):
    """
        compute the volume of a hexahedron element using Gauss integration
                      v3------v7
                     /|      /|
                    v0------v4|
                    | |     | |
                    | v2----|-v6
             y ^    |/      |/
               |    v1------v5
               --->
              /    x
             z   
    """
    assert tch.is_tensor(nodesXYZ)
    assert nodesXYZ.size()[0] == 8 or nodesXYZ.size()[1] == 3

    ncNode = tch.tensor([  # the natural coordinates of 8 nodes of element
        [-1,  1,  1], [-1, -1,  1], 
        [-1, -1, -1], [-1,  1, -1], 
        [ 1,  1,  1], [ 1, -1,  1], 
        [ 1, -1, -1], [ 1,  1, -1],
    ], dtype=tch.int)

    ### integrate the volume by evenly distributed integration-points
    volumes = 0.
    tmp = 1./3.**0.5
    integPoints = [ 
        [-tmp, -tmp, -tmp],
        [-tmp, -tmp,  tmp],
        [-tmp,  tmp, -tmp],
        [-tmp,  tmp,  tmp],
        [ tmp, -tmp, -tmp],
        [ tmp, -tmp,  tmp],
        [ tmp,  tmp, -tmp],
        [ tmp,  tmp,  tmp],
    ]
    integWeight = [1. for _ in range(8)]
    for i, integPoint in enumerate(integPoints):
        ### natural coordinates of the integration points
        naturalCoo = \
            tch.tensor(integPoint, dtype=nodesXYZ.dtype, requires_grad=True)

        ### corresponding weights (shape function) of 8 vertexes
        weights = tch.tensor([0.125 for _ in range(8)], dtype=nodesXYZ.dtype)
        for nodeId in range(8):
            for dm in range(3):
                weights[nodeId] *= 1. + ncNode[nodeId, dm] * naturalCoo[dm]
        globalCoo = tch.einsum("i, ij -> j", weights, nodesXYZ)

        ### Jacobian matrix obtained by auto grad of pyTorch
        jacobi = tch.tensor([])
        for dm in range(3):
            tuple_ = tch.autograd.grad(globalCoo[dm], naturalCoo, retain_graph=True)
            jacobi = tch.cat((jacobi, tuple_[0]))
        jacobi = jacobi.reshape((-1, 3))

        ### determinent of Jacobian matrix, with Gauss weight of 1
        volumes += float(tch.det(jacobi)) * integWeight[i]
    return volumes

此方法的计算结果是 8.4167，下面将与另一种方法进行比较。

方法2：把外表面分割成三角形之后计算体积

基本思路是把六面体各个面上的四边形分割成三角形，把六面体整体看作是由外围的三角形面片构成的物体，利用我上一篇博客中讲到的三角形面片包围的物体的体积计算方法来得到体积。
需要注意的关键点有两点：

1. 外表面的四边形是空间四边形，四个点不在一个面内，因此线性插值后构成的是一个曲面（不是平面）。为了尽可能用三角形接近这个曲面，可以在四边形的中心取一个中心点（中心点必然落在线性插值得到的曲面上），四边形的4个顶点和这个中心点相连，构成4个三角形
2. 每个三角形的取点顺序需要按照右手定则指向六面体的外侧方向

图5：把外表面分割成三角形之后计算体积

用方法2计算一下不规则六面体的体积，与上面的方法1计算同一个六面体，看看两者数值的差异。方法2计算出来的体积是8.4167，两种方法的相对误差是1.8e-7, 几乎一模一样，说明两种方法都能得到精确结果。

方法2的代码

方法2的代码如下，同样用到了行列式，但这一次是把三角形三个点的空间坐标放一起构成矩阵然后计算行列式

import torch as tch

def eleVol_triangles(nodesXYZ):
    """
        use stl volume methods,
        partition the surface of the element into triangles,
        so that element is surrounded by triangle facets,
        use the method to compute volume of 3D object constructed by triangular facets
    """ # author: mohanxuan
    if not tch.is_tensor(nodesXYZ):
        raise ValueError("not tch.is_tensor(nodesXYZ)")
    if nodesXYZ.size()[0] != 8 or nodesXYZ.size()[1] != 3:
        raise ValueError("nodesXYZ.size()[0] != 8 or nodesXYZ.size()[1] != 3")

    facets = [  # each facet must point to the outside of the element
        [0, 3, 2, 1], [4, 5, 6, 7],
        [0, 4, 7, 3], [1, 2, 6, 5],
        [0, 1, 5, 4], [2, 3, 7, 6],
    ]
    triangles = []
    for facet in facets:
        xyzs = tch.tensor(
            [nodesXYZ[nodeId].tolist() for nodeId in facet]
        )
        center = tch.einsum("ij -> j", xyzs) / len(xyzs)
        for idx0 in range(len(facets[0])):
            idx1 = idx0 + 1 if idx0 + 1 < len(facets[0]) else 0
            triangles.append([
                xyzs[idx0].tolist(),
                xyzs[idx1].tolist(),
                center.tolist()
            ])
    triangles = tch.tensor(triangles)
    return sum(tch.det(x) for x in triangles) / 6.