一种基于有限元和无网格耦合的柔性物体实时切割仿真方法

摘要

本发明涉及一种基于有限元和无网格耦合的柔性物体实时切割仿真方法，主要针对虚拟手术实时交互式软组织形变与切割仿真算法进行研究并实现，该方法包括以下步骤：通过基于材质的多子域划分及体素化、各子域内有限元建模、各子域间支持多材质的无网格建模、无网格区域及相关有限元子域的耦合、切割操作响应。本方法对仿真算法进行并行化设计，借助GPU强大的计算能力，达到实时交互式效率。

说明书

技术领域

本发明属于物理仿真技术领域，具体涉及基于有限元和无网格耦合的柔性物体实时切割 仿真方法。

背景技术

在计算机图形学领域，形变模型最早由Terzopolous等人引入。关于逼真模拟形变物体 的多种方法的很好的总数可以查阅。例如，边界元模型，自适应和多尺度方法，无网格技术 以及有限元方法已经被提出。奥布赖恩和霍金斯描述的基于有限元模拟脆性断裂，后来扩大 到韧性断裂。Nesme等提出了一种考虑到粗元素拥有不同材料属性的复合元公式。

切割形变物体的四面体划分方法有Bielser等人引入。为了减少病态元的数量，Nienhuys 和van der Stappen提出了沿元素表面切割。Cotinet和Forest等人删除了被切割的元素。降低 病型元素的平滑切割被Nienhuys等人通过自适应模型边和面匹配切割面实现。通过将划分 限制为几种细分模式，切割所造成的额外模拟元素的数量可以减少。此方法的多分辨率解法 由Ganovelli等人提出。虚拟节点算法通过沿切割两侧复制仿真元素和重新分配材质成分避 免了病型元素。

Wicke等引入了多面体细分，分割出事的四面体成多面体然后再进一步细分这些元素。 扩展有限元法通过特定的基函数丰富了有限元模型来扑捉仿真元素的间断。Nienhuys和van  der Stappen分别实现了扩展有限元法在虚拟手术和切割2D薄壳。考虑到在粗弹性仿真中的 好材料元素，Sifakis剪裁高分辨率材料边界表面模型替代粗仿真模型。

基于八叉树的液体和气体的物理模拟已经发表。限制和不受限制的八叉树方法都在使 用。为了实现小规模细节的高分辨率，焦点之一是派生自适应差分离散方程。Haber和 Heldmann引入了在受限制八叉树上求解偏微分方程数值解的有限元离散方法。

发明内容

本发明要解决的技术问题是：提供一种对于多材质模型的柔性形变和切割进行物理仿真 方法，并利用GPU硬件的计算能力，提高了计算与绘制效率，该方法主要利用有限元方法 对统一材质进行建模同时利用无网格方法对多材质区域进行建模，通过将两种方法进行耦合 达到了对多材质模型的柔性形变和切割的物理仿真，并且通过对算法的并行优化使得物理仿 真达到交互式效率。

本发明解决上述技术问题的技术方案为：一种基于有限元和无网格耦合的柔性物体实时 切割仿真方法，包括如下步骤：

步骤（1）、通过基于材质的多子域划分将仿真对象划分成具有各项同性的子域；

步骤（2）、各子域依据仿真精度要求进行独立的层次细节体素化生成基于八叉树六面体 网格；

步骤（3）、各子域根据各自对应的六面体网格进行有限元建模；

步骤（4）、在各子域间建立无网格采样区域，建立基于材质距离的平滑距离场；

步骤（5）、在无网格采样区域，基于材质距离场进行无网格方法建模；

步骤（6）、对无网格区域及所相关的有限元子域进行耦合；

步骤（7）、当切割操作发生时，对切割所涉及到的有限元区域进行网格细分并通过克隆 网格体现拓扑变化，并在新生成的有限元区域上重新进行无网格方法建模。

进一步的，步骤（1）包括：读取OBJ模型数据，根据材质边界曲面函数将对模型建立 子域划分投票函数（voting function）。

进一步的，步骤（2）包括：建立模型整体的AABB包围盒作为基于八叉树的六面体网 格根节点，如果当前六面体空间内包含模型元素：顶点、三角面片和材质边界曲面，且没有 达到最高层次精度，则对该六面体进行剖分并对新生成的8个子六面体重复上述操作；当所 有六面体网格均符合上述条件时，停止操作，并对八叉树的叶子节点上的六面体元素使用其 中心点作为输入参数调用步骤（1）中的子域划分投票函数确定其所属的子域；该步骤（2） 最终生成了一个反应模型表面和材质分界特征的层次化八叉树六面体网格。

进一步的，步骤（3）包括：为了提高计算效率，采用线性基，并生成六面体有限元形 函数以及形函数的三维空间导数，并以此构建形函数矩阵和应变矩阵；根据当前域所属的材 质属性构建材质矩阵，所述的材质属性为杨氏模量和泊松比，依据有限元公式组织局部刚度 矩阵、质量矩阵和施力向量。

进一步的，步骤（4）包括：对包含材质边界的六面体网格取其二环邻域并将所涉及的 六面体网格作为无网格方法的积分背景网格，随机均匀分配采样点，并为每一个高斯积分点 建立局部影响域内采样点间的材质距离场。

进一步的，步骤（5）包括：根据步骤（4）所得到的距离场生成针对高斯积分点的移动 最小二乘近似（Moving Least Square Approximation）的权函数，并得到满足移动最小二乘近 似的位移函数的形函数，从而得到了对多材质区域的无网格方法建模。

进一步的，步骤（6）包括：找出同属于无网格区域和有限元区域的六面体网格，将步 骤（3）和步骤（5）所得到的位移函数的形函数通过斜坡函数（Ramp Function）进行联立耦 合形成耦合区域的形函数，并根据所生成的形函数生成局部刚度矩阵、质量矩阵和施力向量。

进一步的，步骤（7）包括：将切割工具建模为有三角面片组成的曲面，当切割操作发 生时，通过六面体网格与三角面片做碰撞检测，能够得到被切割的网格，对于粗层次的网格 首先要对其进行网格细分以达到仿真精度，对于被切割的六面体进行克隆操作以满足切割所 造成的拓扑变化，最后对受影响的区域进行步骤（3）和步骤（4），更新仿真方程。

本发明的原理在于：

基于有限元和无网格耦合的柔性物体实时切割仿真方法，包括如下步骤：

1）基于材质的多子域划分及网格化

步骤（1）、读取OBJ模型数据，根据材质边界曲面函数将对模型建立子域划分投票函数 （voting function）。对于一个网格元素而言，当它只满足一个投票函数时，它被认为是一个 内部元素（Interior Element），当它满足多于一个投票函数时，它被认为是边界元素（Boundary  Element）,否则被认为是外部元素（Exterior Element）不参与运算。

步骤（2）、建立模型整体的AABB包围盒作为基于八叉树的六面体网格根节点，如果当 前六面体空间内包含模型元素（顶点、三角面片、材质边界曲面）且没有达到最高层次精度， 则对该六面体进行剖分并对新生成的8个子六面体重复上述操作，在剖分的过程中，为了防 止出现悬挂点（Hanging Node）,需要保证每一个八叉树叶子节点与其邻接节点的层次精度至 多相差一个级别，否则也要进行剖分。当所有六面体网格均符合上述条件时，停止操作，并 对八叉树的叶子节点上的六面体元素使用其中心点作为输入参数调用步骤（1）中的子域划 分投票函数确定其所属的子域。该步骤最终生成了一个反应模型表面和材质分界特征的层次 化八叉树六面体网格。

2）基于网格的物理建模

步骤（3）、为了提高计算效率，本方法采用线性基（p＝[1,x,y,z]），并生成六面体有限 元形函数以及形函数的三维空间导数，并以此构建形函数矩阵和应变矩阵。根据当前域所属 的材质属性（杨氏模量和泊松比）构建材质矩阵。依据有限元公式组织局部刚度矩阵、质量 矩阵和施力向量。由于步骤（2）所生成的六面体网格形状是相同的只是尺寸不同，本方法 只需要计算上述矩阵一次，通过乘以不同的放缩因子得到不同尺寸的对应矩阵。

步骤（4）、对包含材质边界的六面体网格取其二环邻域并将所涉及的六面体网格作为无 网格方法的积分背景网格，随机均匀分配采样点，并为每一个高斯积分点建立局部影响域内 采样点间的材质距离场。距离场测量采用测地线方法，利用单元最短路算法对支持域内采样 点进行计算，为了体现多材质特性，本方法提出一种包含材质属性的两点间距离度量函数， 用此函数可以得到支持多材质的平滑距离场。

步骤（5）、根据步骤（4）所得到的距离场生成针对高斯积分点的移动最小二乘近似 （Moving Least Square Approximation）的权函数，并得到满足移动最小二乘近似的位移函数 的形函数，从而得到了对多材质区域的无网格方法建模。

步骤（6）、找出同属于无网格区域和有限元区域的六面体网格，将步骤（3）和步骤（5） 所得到的位移函数的形函数通过斜坡函数（Ramp Function）进行联立耦合形成耦合区域的形 函数，并根据所生成的形函数生成局部刚度矩阵、质量矩阵和施力向量。将步骤（3）和步 骤（6）所得到的局部矩阵根据全局自由度发布到全局，由于有限元方法和无网格方法具有 紧支特性，所以得到的全局矩阵具有稀疏特性，利用预条件共轭梯度方法（Preconditioned  Conjugate Gradient Method）求解，得到了所有采样点的位移向量。为了进行动态仿真，需要 采用NewMark隐式积分方法，利用所得到的位移向量计算采样点的速度向量和加速度向量。

3）切割操作响应

步骤（7）、本方法将切割工具建模为有三角面片组成的曲面，当切割操作发生时，通过 六面体网格与三角面片做碰撞检测，可以得到被切割的网格，对于粗层次的网格首先要对其 进行网格细分以达到仿真精度，对于被切割的六面体进行克隆操作以满足切割所造成的拓扑 变化，最后对受影响的区域进行步骤（3）和步骤（4），更新全局方程组。对于切割造成的 模型表面的几何变化采用Delaunay三角化进行更新。

本发明与现有技术相比的有点在于：

本发明能够对多材质柔性物体的变形和切割进行实时物理仿真。本发明主要有两点贡 献：第一，提出了一种基于材质距离的无网格建模方法，该方法可以有效地对多材质和几何 形状复杂的物体进行仿真。第二，给出了一种基于多尺度的局部有限元与全局无网格耦合方 法，该方法可以对多材质柔性物体进行材质区域划分，根据区域的重要性进行多尺度建模， 并且支持实时形变和交互式切割操作。

附图说明

图1为算法整体流程图；

图2为模型初始化流程图；其中，（a）初始模型；（b）材质分区后模型；（c）基于八叉 树体素化网格；

图3为层次化六面体网格结构剖面图；

图4为测地线计算距离场示意图；

图5为平滑距离场示意图；

图6为切割响应操作示意图；其中，（a）初始网格；（b）发生切割的细分网格和克隆元；

图7为表面三角网格分割算法示意图；其中，（a）分割前表面三角网格；（b）分割后表 面三角网格；

图8为多材质形变物体仿真输出；其中，（a）纹理输出；（b）仿真网格输出；

图9为多材质形变物体切割仿真输出；其中，（a）纹理输出；（b）表面三角网格输出。

具体实施方式

图1给出了基于有限元和无网格耦合的柔性物体实时切割仿真过程的总体处理流程，下 面结合其他附图及具体实施方式进一步说明本发明。

本发明提供一种基于有限元和无网格耦合的柔性物体实时切割仿真方法，主要步骤介绍 如下：

步骤（1）根据输入的OBJ模型（见图2（a））和材质边界方程建立子域划分投票函数 （voting function，v_i(x,y,z),i＝0..m,m为材质边界个数）。其中投票函数以空间点坐标为输入 参数，当空间坐标包含在所表示的材质域内返回1，否则返回0。图2（b）表示了当前模型 自上而下拥有四种材质，以不同颜色表示不同的材质属性。

步骤（2）建立模型整体的AABB包围盒作为基于八叉树的六面体网格根节点，如果当 前六面体空间内包含模型元素（顶点、三角面片、材质边界曲面）且没有达到最高层次精度， 则对该六面体进行剖分并对新生成的8个子六面体重复上述操作，在剖分的过程中，为了防 止出现悬挂点（Hanging Node）,需要保证每一个八叉树叶子节点与其邻接节点的层次精度至 多相差一个级别，否则也要进行剖分。当所有六面体网格均符合上述条件时，停止操作，并 对八叉树的叶子节点上的六面体元素使用其中心点作为输入参数调用步骤（1）中的子域划 分投票函数确定其所属的子域。该步骤最终生成了一个反应模型表面和材质分界特征的层次 化八叉树六面体网格。图2（c）表示了全局六面体网格，图3表示了模型内部的层次化六面 体网格。

步骤（3）本方法采用线性基（p＝[1,x,y,z]），并生成六面体有限元形函数（公式1）以 及形函数的三维空间导数（公式2），其中ζ,η,μ为局部坐标，f为施加的外力，并以此构建 形函数矩阵和应变矩阵。

$N_1=\frac{1}{8}(1-\zeta )(1-\eta )(1-\mu ),N_2=\frac{1}{8}(1+\zeta )(1-\eta )(1-\mu )$

$N_3=\frac{1}{8}(1+\zeta )(1+\eta )(1-\mu ),N_4=\frac{1}{8}(1-\zeta )(1+\eta )(1-\mu )$

$N_5=\frac{1}{8}(1-\zeta )(1-\eta )(1+\mu ),N_6=\frac{1}{8}(1+\zeta )(1-\eta )(1+\mu )---\left(1\right)$

$N_7=\frac{1}{8}(1+\zeta )(1+\eta )(1+\mu ),N_8=\frac{1}{8}(1-\zeta )(1+\eta )(1+\mu )$

N＝[N₁,N₂,N₃,N₄,N₅,N₆,N₇,N₈]

$B=\left(\begin{array}{ccc}\partial /\partial x& 0& 0\\ 0& \partial /\partial y& 0\\ 0& 0& \partial /\partial z\\ \partial /\partial y& \partial /\partial x& 0\\ 0& \partial /\partial z& \partial /\partial y\\ \partial /\partial z& 0& \partial /\partial x\end{array}\right)·\left(\begin{array}{ccc}N& 0& 0\\ 0& N& 0\\ 0& 0& N\end{array}\right)---\left(2\right)$

根据当前域所属的材质属性（杨氏模量和泊松比）构建材质矩阵。依据有限元公式组织 局部刚度矩阵（公式3）、质量矩阵（公式4）和施力向量（公式5）。

$k_e={\int }_{{\Omega }_e}{B}^T\mathrm{CBd\Omega }---\left(3\right)$

$m_e={\int }_{{\Omega }_e}\rho N^T\mathrm{Nd\Omega }---\left(4\right)$

$f_e={\int }_{{\Gamma }_e}{N}^T\mathrm{fd\Gamma }---\left(5\right)$

由于步骤（2）所生成的六面体网格形状是相同的只是尺寸不同，本方法只需要计算上 述矩阵一次，通过乘以不同的放缩因子得到不同尺寸的对应矩阵。

步骤（4）对包含材质边界的六面体网格取其二环邻域并将所涉及的六面体网格作为无 网格方法的积分背景网格，随机均匀分配采样点，并为每一个高斯积分点建立局部影响域内 采样点间的材质距离场。距离场测量采用测地线方法，利用单元最短路算法对支持域内采样 点进行计算（见图4）。

$\mathrm{dist}=\frac{2\left|\right|P_1-P_2\left|\right|}{{\mathrm{stiff}}_1+{\mathrm{stiff}}_2}---\left(6\right)$

为了体现多材质特性，本方法提出一种包含材质属性的两点间距离度量函数，用此函数 可以得到支持多材质的平滑距离场（公式6，其中P_i表示六面体中点坐标，stiff_i表示六面体 的材质信息）。图5表示了步骤（4）得到的当前支持域内的平滑距离场。

步骤（5）根据步骤（4）所得到的距离场生成针对高斯积分点的移动最小二乘近似（Moving  Least Square Approximation）的权函数（公式7），并得到满足移动最小二乘近似的位移函数 的形函数（公式8，9，10），从而得到了对多材质区域的无网格方法建模（公式3，4，5）。

$w\left(r\right)=\left(\begin{array}{cc}2/3-4r^2+4r^3& r\le 1/2\\ 4/3-4r+4r^2-4r^3/3& 1/2<r\le 1\\ 0& r>1\end{array}\right)---\left(7\right)$

$A\left(x\right)=\underset{I=1}{\overset{N}{\Sigma }}{w}_I\left(x\right)p\left(x_I\right)p^T\left(x_I\right)---\left(8\right)$

B(x)＝[w₁(x)p(x₁),w₂(x)p(x₂)…w_N(x)p(x_N)]   （9）

$\Phi (x,\bar{x})=p^T\left(\bar{x}\right)A^{-1}\left(x\right)B\left(x\right)$

步骤（6）找出同属于无网格区域和有限元区域的六面体网格，将步骤（3）和步骤（5） 所得到的位移函数的形函数通过斜坡函数（公式11）进行联立耦合形成耦合区域的形函数（公 式12，Φ为无网格形函数，N为有限元形函数），并根据所生成的形函数生成局部刚度矩阵、 质量矩阵和施力向量。

$R\left(x\right)=\underset{I}{\Sigma }{\Phi }_I\left(x\right),x_I\in {\Gamma }_{\mathrm{EFG}}---\left(11\right)$

${\tilde{N}}_I=\left(\begin{array}{cc}[1-R]N_I+R{\Phi }_I& X_I\in {\Omega }_I\\ R{\Phi }_I& X_I\notin {\Omega }_I\end{array}\right)---\left(12\right)$

$M\stackrel{··}{u}+C\stackrel{·}{u}+\mathrm{Ku}=f---\left(13\right)$

将步骤（3）和步骤（6）所得到的局部矩阵根据全局自由度发布到全局，由于有限元方 法和无网格方法具有紧支特性，所以得到的全局矩阵具有稀疏特性。将矩阵带入拉格朗日运 动方程（公式13，M是质量矩阵，C是阻尼矩阵，K是刚度矩阵，f是施力向量），利用预 条件共轭梯度方法（Preconditioned Conjugate Gradient Method）求解，得到了所有采样点的 位移向量。为了进行动态仿真，需要采用NewMark隐式积分方法，利用所得到的位移向量 计算采样点的速度向量和加速度向量。

步骤（7）将切割工具建模为有三角面片组成的曲面，当切割操作发生时，通过六面体 网格与三角面片做碰撞检测，可以得到被切割的网格，对于粗层次的网格首先要对其进行网 格细分以达到仿真精度，对于被切割的六面体进行克隆操作以满足切割所造成的拓扑变化 （见图6），最后对受影响的区域进行步骤（3）和步骤（4），更新全局方程组。对于切割造 成的模型表面的几何变化采用Delaunay三角化进行更新（见图7）。

本方法的实现使用的软件平台为Microsoft visual studio 2010与OpenGL，使用了CUDA 来加速并行算法的计算效率。硬件平台为3.4GHz Inter(R) Core(TM) i7-2600 CPU、4GB内存 以及NVIDIA GeForce GTX570 GPU。方法效果图如图8，9所示。