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

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

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

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

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

如何实现映射？

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

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

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

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

某个点从标准立方体对应的自然坐标 ξ⃗=(ξ,η,ζ)\vec{\xi}=(\xi,\ \eta,\ \zeta)ξ=(ξ, η, ζ) 映射到不规则六面体对应的全局坐标 x⃗=(x,y,z)\vec{x}=(x,\ y,\ z)x=(x, y, z) 的映射关系是：
x⃗=∑i=18Ni(ξ,η,ζ)⋅ξi⃗,whereNi(ξ,η,ζ)=18(1+ξiξ)(1+ηiη)(1+ζiζ),vertexes′coordinates(ξi,ηi,ζ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)\vec{x}=\sum_{i=1}^{8}{N_{i}(\xi,\ \eta,\ \zeta)\cdot\vec{\xi_{i}}}, \ \ where\ N_{i}(\xi,\ \eta,\ \zeta)=\frac{1}{8}(1+\xi_{i}\xi)(1+\eta_{i}\eta)(1+\zeta_{i}\zeta), \\ vertexes'\ coordinates\ (\xi_{i},\ \eta_{i},\ \zeta_{i}) = (-1,\ 1,\ 1), or\ (-1,\ -1,\ 1),\\ \hspace{14em}or\ (-1,\ -1,\ -1),\ or\ (-1,\ 1,\ -1),\\ \hspace{14em}or\ (1,\ 1,\ 1),\ or\ (1,\ -1,\ 1),\\ \hspace{14em}or\ (1,\ -1,\ -1),\ or\ (1,\ 1,\ -1)x=i=1∑8Ni(ξ, η, ζ)⋅ξi, where Ni(ξ, η, ζ)=81(1+ξiξ)(1+ηiη)(1+ζiζ),vertexes′ coordinates (ξi, ηi, ζ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)
， N_{i} 是立方体内部点p对于顶点 i 的权重，如果内部点p与顶点 iii 之间距离越近，那么权重 NiN_{i}Ni 就越高。权重满足两个条件是：

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

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

有了坐标映射，就可以知道不规则六面体相对于立方体在三个不同方向上的长度伸缩关系，从而得到不规则六面体的体积。标准立方体的自然坐标方向的微矢量 dξ⃗,dη⃗,dζ⃗\vec{d\xi},\ \vec{d\eta},\ \vec{d\zeta}dξ, dη, dζ 可以分别映射到全局坐标中得到新矢量 dx⃗,dy⃗,dz⃗\vec{dx},\ \vec{dy},\ \vec{dz}dx, dy, dz :
dξ⃗=(dξ,0,0)⟶maptodx⃗=(∂x∂ξdξ,∂y∂ξdξ,∂z∂ξdξ)dη⃗=(0,dη,0)⟶maptody⃗=(∂x∂ηdη,∂y∂ηdη,∂z∂ηdη)dζ⃗=(0,0,dζ)⟶maptodz⃗=(∂x∂ζdζ,∂y∂ζdζ,∂z∂ζdζ)\vec{d\xi}=(d\xi,\ 0,\ 0)\hspace{1em}\mathop{\longrightarrow}^{map\ to}\hspace{1em}\vec{dx}=(\frac{\partial x}{\partial \xi}d\xi,\ \frac{\partial y}{\partial \xi}d\xi,\ \frac{\partial z}{\partial \xi}d\xi) \\ \vec{d\eta}=(0,\ d\eta,\ 0)\hspace{1em}\mathop{\longrightarrow}^{map\ to}\hspace{1em}\vec{dy}=(\frac{\partial x}{\partial \eta}d\eta,\ \frac{\partial y}{\partial \eta}d\eta,\ \frac{\partial z}{\partial \eta}d\eta) \\ \vec{d\zeta}=(0,\ 0,\ d\zeta)\hspace{1em}\mathop{\longrightarrow}^{map\ to}\hspace{1em}\vec{dz}=(\frac{\partial x}{\partial \zeta}d\zeta,\ \frac{\partial y}{\partial \zeta}d\zeta,\ \frac{\partial z}{\partial \zeta}d\zeta) dξ=(dξ, 0, 0)⟶map todx=(∂ξ∂xdξ, ∂ξ∂ydξ, ∂ξ∂zdξ)dη=(0, dη, 0)⟶map tody=(∂η∂xdη, ∂η∂ydη, ∂η∂zdη)dζ=(0, 0, dζ)⟶map todz=(∂ζ∂xdζ, ∂ζ∂ydζ, ∂ζ∂zdζ)
把不规则六面体切分成由上面的微矢量 dx⃗,dy⃗,dz⃗\vec{dx},\ \vec{dy},\ \vec{dz}dx, dy, dz 构成的平行六面体微元，对这些微元进行积分可以得到体积：
volume=∫Ω~(dx⃗×dy⃗)⋅dz⃗=∫Ω∣∂x∂ξdξ∂y∂ξdξ∂z∂ξdξ∂x∂ηdη∂y∂ηdη∂z∂ηdη∂x∂ζdζ∂y∂ζdζ∂z∂ζdζ∣=∫Ω∣∂x∂ξ∂y∂ξ∂z∂ξ∂x∂η∂y∂η∂z∂η∂x∂ζ∂y∂ζ∂z∂ζ∣dξdηdζvolume = \int_{\tilde{\Omega}}{(\vec{dx}\times\vec{dy})\cdot\vec{dz}} = \int_{\Omega}{\left|\begin{matrix} \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{matrix}\right|} = \int_{\Omega}{ \left|\begin{matrix} \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{matrix}\right| d\xi d\eta d\zeta } volume=∫Ω~(dx×dy)⋅dz=∫Ω∣∣∂ξ∂xdξ∂η∂xdη∂ζ∂xdζ∂ξ∂ydξ∂η∂ydη∂ζ∂ydζ∂ξ∂zdξ∂η∂zdη∂ζ∂zdζ∣∣=∫Ω∣∣∂ξ∂x∂η∂x∂ζ∂x∂ξ∂y∂η∂y∂ζ∂y∂ξ∂z∂η∂z∂ζ∂z∣∣dξdηdζ
即关于雅可比矩阵的行列式值的积分。
这个积分会在积分点的增加之下最终收敛到一个精确的结果。
这里做一个测试，随便取一个六面体单元，其八个顶点的坐标分别是：

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

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

方法1的代码

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

import torch as tchdef eleVol_gaussInteg(nodesXYZ):"""compute the volume of a hexahedron element using Gauss integrationv3------v7/|      /|v0------v4|| |     | || v2----|-v6y ^    |/      |/|    v1------v5--->/    xz   """assert tch.is_tensor(nodesXYZ)assert nodesXYZ.size()[0] == 8 or nodesXYZ.size()[1] == 3ncNode = 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-pointsvolumes = 0.tmp = 1./3.**0.5integPoints = [ [-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 pointsnaturalCoo = \tch.tensor(integPoint, dtype=nodesXYZ.dtype, requires_grad=True)### corresponding weights (shape function) of 8 vertexesweights = 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 pyTorchjacobi = 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 1volumes += float(tch.det(jacobi)) * integWeight[i]return volumes

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

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

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

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

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

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

方法2的代码

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

import torch as tchdef 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: mohanxuanif 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 0triangles.append([xyzs[idx0].tolist(),xyzs[idx1].tolist(),center.tolist()])triangles = tch.tensor(triangles)return sum(tch.det(x) for x in triangles) / 6.

六面体单元的体积计算方法相关推荐

matlab patch 六面体,《有限元基础教程》_【MATLAB算例】4.8.2(1) 基于8节点六面体单元的空间块体分析(Hexahedral3D8Node)...
[MATLAB 算例]4.8.2(1) 基于8节点六面体单元的空间块体分析(Hexahedral3D8Node) 如图4-23所示的一个空间块体,在右端部受两个集中力F 作用,其中的参数为: 1051 ...
WEBGL 体积计算方法
WEBGL 体积土方量计算方法测算效果图代码实现本文主要介绍了基于 Cesium 实现体积计算的方法介绍,前几天做项目时候遇到客户提出加入体积计算的项目要求,琢磨了几天,终于搞出来了,在此将一 ...
六面体体积求解(规则不规则)
规则和不规则六面体求解方法1:立方体映射法(不规则六面体与立方体之间的映射) 这种方法的主要思想是把一个不规则六面体映射到一个立方体当中,通过计算不规则六面体与立方体之间的体积比例(分解为三维空间中 ...
【JY】有限单元分析的常见问题及单元选择
因你精彩即刻关注再次整理了下笔记,在看本文前,可以先看下:[JY]有限元分析的单元类型分享一波~ 我们常用的有限元方法有以下非常需要注意的要点(特别是实体单元的应用):剪切锁死.体积锁死.沙漏模式 ...
FVM in CFD 学习笔记_第6章_有限体积网格
学习自F. Moukalled, L. Mangani, M. Darwish所著The Finite Volume Method in Computational Fluid Dynamics - ...
3 设置网格数的大小_流体仿真中，六面体（Hex）网格的求解效率真的比四面体（Tet）高”很多”么？...
流体仿真中,六面体(Hex)网格与四面体(Tet)网格的争论一直伴随着整个CFD的发展过程,坊间也流传着许许多多关于六面体网格.结构化网格.四面体网格.甚至是Cutcell网格等相关内容的种类繁多的观 ...
流体仿真中，六面体（Hex）网格的求解效率真的比四面体（Tet）高”很多”么？
作者 | 张杨流体仿真中,六面体(Hex)网格与四面体(Tet)网格的争论一直伴随着整个CFD的发展过程,坊间也流传着许许多多关于六面体网格.结构化网格.四面体网格.甚至是Cutcell网格等相关内 ...
解决LS-DYNA中负体积方法
对于承受很大变形的材料,例如泡沫.由于在计算过程中,某个单元可能变得非常扭曲以至于单元的体积出现负值,负体积可能发生在材料未达到失效标准前.LS-DYNA中计算得到负体积(negativevolume ...
c语言质数和合数事例,新人教版三年级数学下册第八单元单元教学计划
第八单元数学广角--搭配 (二) 新知识点: 1.简单事物的排列数. 2.简单事物的组合数. 教学要求: 1.联系学生的生活实际,使学生通过观察.猜测.试验等活动,找出简单事物的排列数和组合数. 2 ...
LS-DYNA常用单元公式选择指南
单元公式 LS-DYNA是一种通用有限元程序,用于分析结构的大变形静力和动力响应,包括结构耦合到流体.主要的解决方法是基于显式动力学.所有的有限元模型,都必然涉及网格划分,ANSYS LS-DYNA ...

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