规则和不规则六面体求解

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

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

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

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

如何实现映射？

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

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

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

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

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

即关于雅可比矩阵的行列式值的积分。

这个积分会在积分点的增加之下最终收敛到一个精确的结果。

这里做一个测试，随便取一个六面体单元，其八个顶点的坐标分别是：

图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有限体积网格,用Matlab实现简单有限体积求解器
用于瞬态对流扩散PDE的简单而通用的FVM求解器一个简单的有限体积工具这段代码是化学/石油工程师开发一种简单的工具来求解对流扩散方程的一般形式的结果: αgeneral / equationt + ...
排水沟槽开挖土方的计算方法（平行相似梯形组成的六面体体积分割计算方法）
排水沟槽是一种利用重力敷设排水管道需要开挖的沟槽.该沟槽的特点是: 1.沟槽的底宽从起点到终点保持不变(管径不变的情况下): 2.沟槽的边坡放坡边坡1:m保持不变."1:m"表示 ...
C语言表达式的求解规则,C语言实现整数四则运算表达式的计算
一.问题重述 [问题描述] 从标准输入中读入一个整数算术运算表达式,如5 - 1 * 2 * 3 + 12 / 2 / 2 = .计算表达式结果,并输出. 要求: 1.表达式运算符只有+.-.*./ ...
3 设置网格数的大小_流体仿真中，六面体（Hex）网格的求解效率真的比四面体（Tet）高”很多”么？...
流体仿真中,六面体(Hex)网格与四面体(Tet)网格的争论一直伴随着整个CFD的发展过程,坊间也流传着许许多多关于六面体网格.结构化网格.四面体网格.甚至是Cutcell网格等相关内容的种类繁多的观 ...
流体仿真中，六面体（Hex）网格的求解效率真的比四面体（Tet）高”很多”么？
作者 | 张杨流体仿真中,六面体(Hex)网格与四面体(Tet)网格的争论一直伴随着整个CFD的发展过程,坊间也流传着许许多多关于六面体网格.结构化网格.四面体网格.甚至是Cutcell网格等相关内 ...
构建规则格网进行体积计算
构建规则格网进行体积计算 1.构建规则格网 1.1生成所有格网点 2.计算体积 2.1计算凸包所包含的所有格网点 2.2插值计算凸包内格网点的高程 2.3计算体积总体步骤: 生成凸包多边形构建规则 ...
python求正方体体积_「高中数学」简单几何体的面积与体积相关知识点整理+例题...
一.知识要点 (一)圆柱.圆锥.圆台的侧面积将侧面沿母线展开在平面上,则其侧面展开图的面积即为侧面面积. 1.圆柱的侧面展开图--矩形圆柱的侧面积 2.圆锥的侧面展开图--扇形圆锥的侧面积 3. ...
puzzle（0112）不规则数独、变种数独
目录一,不规则数独 1,规则 2,对称性含正方形的不规则数独中心对称的不规则数独其他对称不规则数独 3,计算机求解二,变种数独 1,规则 2,对角线数独 3,窗口数独 4,超级窗口数独 5, ...

六面体体积求解(规则不规则)