一种基于特征的CAD体网格模型编辑、优化算法

摘要

本发明公开了一种基于特征的CAD体网格模型编辑、优化算法，包括基于特征的CAD体网格模型编辑和基于特征的CAD四面体网格模型优化；在编辑中，用户可以选择需要保持模型形状特征的一组表面区域，通过对这些区域的平移、旋转等操作，来实现对模型的编辑。在优化时，根据模型变动的区域，首先评估出质量较差的四面体网格单元，对这些单元进行相应的拓扑优化操作，最后再进行网格平滑操作，整个优化过程迭代地进行，最终得到优化后的网格模型数据。本发明可以直接对有限元网格模型进行编辑并优化，可以精确地保持模型原有的特征，完全符合相关工业领域的应用要求。

说明书

技术领域

本发明涉及计算机辅助设计及有限元分析领域，尤其是指一种基于特征的CAD体网格模型编辑、优化算法。

背景技术

计算机辅助设计技术(Computer-Aided Design，CAD)指以计算机作为辅助工具，帮助工程领域的设计人员设计产品的技术。能够有效地加快产品设计分析的速度，提高制造效率及精度，提高经济效益。在计算机辅助设计技术的帮助下，工程技术人员可以实现“无纸化”作业，利用计算机的帮助，高效地完成产品的分析、设计、绘制以及制造等工作。利用计算机辅助设计技术，在节约资源降低设计成本的同时，也极大地加快了产品的设计速度，提高产品质量。美国波音公司从70年代开始就使用CAD技术辅助设计飞机，到90年代生产波音－777客机时，已全部实现无纸化设计和生产，极大地提高了生产效率和数据管理效率。

一般地，产品的生产包含设计、分析与制造几部分。设计人员用CAD软件设计好模型后，由分析人员先把模型网格化，然后用有限元软件进行有限元分析，直到产品性能合格，再去生产。随着软硬件技术的发展，特别是材料技术、3d打印技术的发展，工业4.0和中国制造2025概念的提出，越来越多的工业人员希望把设计与制造这两个独立的过程集成起来。

目前，计算机辅助设计的核心工作是对产品进行几何建模，在利用计算机辅助设计的相关软件对产品设计之后，需要对产品的性能进行验证，检查通过此次设计得到的产品的性能是否可以达到实际应用的要求，这就需要用到计算机辅助工程(Computer-Aided Engineering，CAE)技术。在计算机及辅助工程技术诞生之前，需要对设计的产品研制出具体的样品，并通过实际的物理实验来检验其性能。无法达到要求则需要重新设计、制造、检验，消耗巨大成本的同时，也耗费了很多时间。而借助于计算机辅助工程技术，可以利用计算机来模拟实际的实验环境，对虚拟的样品进行检测，进而在实际的物理实验之前，提前找出大部分的设计缺陷，可以极大地降低产品开发的成本，加快开发速度。

有限元分析技术(Finite Element Analysis，FEA)是计算机辅助工程技术最重要的理论基础。它将连续的区域离散化，生成由特定单元组成的计算模型，便于对其进行数值求解。借助这种技术，工程中大量的仿真问题得以解决。现在，许多基于有限元分析的计算机辅助工程的相关商业软件已经在建筑、机械、航空、航天等多个领域广泛使用，比如ABAQUS、ANSYS、ADINA和Nastran等。

目前，有限元分析技术已经相对比较成熟了，然而，由于历史原因，在计算机辅助设计与有限元分析的集成方面还有一些急需解决的问题。问题主要在于计算机辅助设计与有限元分析系统的数据表示不一致，因此，通过计算机辅助设计的相关软件设计出的几何模型必须要进行一些修改，才能转换为有限元分析系统可以读取的数据。在这个过程中，生成用于合适有限元分析应用的网格往往需要比较长的时间。实际过程中，即使对于模型的修改是局部的，往往也需要重新回到CAD软件中做修改，然后再进行一系列的简化和网格化工作。而在实际的工程设计流程中，往往需要对模型反复的设计与分析，这样就要不停地重复进行模型转换的工作，从而会极大地降低产品开发的效率。

网格模型，是对曲面或空间的一种离散化表示，具体来说，是三维空间中，利用顶点位置以及顶点之间的连接关系对模型的一种表示方式。一般地，要对产品的模型进行性能分析，首先需要对产品进行网格化，得到相应的有限元网格模型。在有限元分析中，四面体网格和六面体网格是常用的两种体网格单元。一般地，由于四面体网格能覆盖整个三维模型，并且表示简单紧凑，而且四面体网格化要比六面体网格化更容易，因此，在实际应用中，四面体网格应用更为广泛，许多商业化的CAE软件都支持四面体网格化。四面体网格被广泛应用于基于物理的变形，模拟仿真，产品质量评估等。如图1所示，在产品设计等过程中，三维模型往往首先通过计算机辅助设计软件生成出来，之后使用有限元分析工具将其离散化为网格模型，进而来评估此次设计后模型的质量，如果评估结果不符合设计要求，则需要重复整个过程。而且在实际的工程设计流程中，设计时模型会经常变化，而在以往的方法中，模型上任何细微的变化，都需要对整个模型重复进行利用计算机辅助设计系统的生成以及有限元分析工具的优化的过程，这是非常费时的。随着社会工业的发展，传统的开发流程明显已不适应社会的发展。

为加快产品开发效率，缩短开发周期，将CAD和CAE两个过程结合起来十分必要，这也是在工业4.0中的必然要求。在这个结合过程中，通常有两个解决方案。其中一个是在CAD软件中引入相应的有限元分析功能，以便在设计后能直接进行有限元分析。目前，一些商业CAD软件已经初步具备部分有限元分析的功能，如SolidWorks、UG/NX等。然而，即使如此，目前这些商用软件也只是把CAD与CAE的数据简单地统一起来。如果模型做了变动，仍然需要对整个模型进行简化和网格化等工作。另一个方法，则是研究直接能编辑有限元网格的方法。

目前，已经存在大量的网格编辑方法。然而，这些方法绝大多数都是针对面网格模型的，在体网格模型的编辑方面，相关的研究工作相对较少。此外，这些方法的适用对象基本都是图形学模型和动画模型，考虑的大多数是编辑前后的光滑性和细节性等视觉上的问题，虽然这些方法对许多模型的变形可以取得较好的效果，但是工业应用上的机械零件CAD模型往往具有许多尖锐的特征(feature)，而且还具有一些编辑语义的局部形状存在，直接应用这些已有的方法，很难保证编辑前后的语义和特征。另外，如果对CAD体网格模型编辑的幅度较大，被编辑区域的网格质量会变差，这就需要研究一种体网格优化算法来对网格进行优化。

发明内容

本发明的目的在于克服现有技术的缺点与不足，提供了一种基于特征的CAD体网格模型编辑、优化算法，该方法可以在编辑的过程中，精确地保持零件模型原有的形状特征，可以实时地对体网格模型的局部形状特征进行编辑，而不必重复进行CAD模型编辑和生成有限元网格等过程，提高了模型设计的效率；另外，该方法能有效、快速地对四面体网格模型进行局部优化，并得到符合有限元分析的体网格模型，保证了编辑后的网格模型能继续满足有限元分析需求。

为实现上述目的，本发明提供的技术方案为：一种基于特征的CAD体网格模型编辑、优化算法，包括以下步骤：

1)基于特征的CAD体网格模型编辑

1.1)读入CAD体网格模型并对外表面区域进行分割，提取相应的外表面区域；

1.2)用户交互式地选取外表面区域，作为要编辑的形状特征的操作区域；

1.3)对用户选定的每一个形状特征的“操作区域”，进行以下计算：

针对当前“操作区域”，建立局部坐标系，该“操作区域”对应的“操作点”作为当前局部坐标系的原点；计算并记录网格所有顶点在当前局部坐标系中的坐标，便于通过操作区域、操作点来对整个网格进行编辑；

1.4)对每个“操作点”，构建优化方程，并利用Eigen求解优化方程，求得每个“操作点”对网格上每个顶点的权重；并对“操作点”的权重进行归一化操作，使得网格上的每个顶点，各个操作点对其的权重总和为1；

1.5)记录用户对每个“操作区域”进行的旋转平移编辑操作，将相应的操作作用在以对应“操作点”为原点的局部坐标系的每个顶点坐标上；将每个“操作点”为原点的局部坐标系下的每个顶点的坐标数据变换到全局坐标系即最初模型所处的坐标系，将每个顶点的坐标数据乘以其对应“操作点”的权重，加和得到最终顶点的坐标数据，这样，每个顶点的坐标都得以更新，模型的编辑操作也已经完成；

2)基于特征的CAD四面体网格模型优化

2.1)如果执行的迭代次数小于迭代阈值，则执行步骤2.2)，否则，执行步骤2.3)；

2.2)对于网格的每一条边，如果其边长大于4/5*target_edge_length，其中target_edge_length是目标边长，则对其执行边分裂操作；如果其边长小于4/5*target_edge_length，那么对其执行边收缩操作；本步骤的两种操作将会对网格模型的拓扑单元进行增加与删除；然后执行步骤2.4)；

2.3)检查网格的每一条边与每一个面，对于符合条件的边与面，对其邻接的四面体单元执行四面体翻转的操作，然后执行步骤2.4)；

2.4)执行网格顶点坐标的平滑操作，过程如下：

2.4.1)对网格外表面上的顶点进行平滑，平滑时只根据这个顶点的一邻域中在外表面上的顶点来确定顶点平滑后的坐标；平滑过程中注意保持模型的语义形状特征，位于网格特征线上的顶点不应该进行平滑操作；另外，如果网格的外表面不是平面，那么表面顶点平滑后的位置往往不在外表面上，这就需要在平滑之后，对外表面上的顶点进行一个向原始网格即输入网格的外表面投影的操作；

2.4.2)对网格的内部顶点进行平滑；

2.5)如果执行了前四个步骤后，网格的拓扑元素数目相比于执行前变化不大，或者达到执行次数的上限，停止，否则，执行步骤2.1)。

在步骤1.4)中，所述的操作点权重的计算，其方法为：

定义M＝(V,E,F,H)为具有n个顶点的四面体网格，V,E,F,H分别为该网格的顶点集合、边集合、面集合以及四面体集合，对于顶点i∈M，其笛卡尔坐标被表示为v_i＝(x_i,y_i,z_i)；定义R_j＝(V_j,E_j,F_j),j＝1,...,m为用户选取的所有的外表面操作区域，共m个，P_j,j＝1,...,m为每个操作区域所对应的操作点；w_j是操作点j在网格所有顶点上的权重集合，以向量的形式表示，w_j(v_j)则是操作点j对顶点i的权重；

由于Ls包含了网格的所有拓扑邻接信息，并且，操作点的权重，在当前操作区域所包含的顶点上应该为1，在非当前操作区域包含的顶点上，应该为0；由此，能够将求解操作点权重的问题建模为一个最优化问题，通过对这个最优化问题的求解来求得当前操作点对网格上所有顶点的权重，见下面的公式：

$\underset{w_j,j=1,\mathrm{...},m}{\mathrm{argmin}}\left({||L_s{w}_j||}^2\right)---\left(1\right)$

满足：

$w_j\left(v_j\right)=1,\forall i\in R_j---\left(2\right)$

$w_j\left(v_i\right)=0,\forall i\in (R-R_j)---\left(3\right)$

$\underset{j=1}{\overset{m}{\Sigma }}{w}_j\left(v_i\right)=1,\forall i\in M---\left(4\right)$

$0\le w_j\left(v_i\right)\le 1,\forall i\in M---\left(5\right)$

其中，R为用户选取的所有操作区域的集合，R_j＝(V_j,E_j,F_j),j＝1,...,m的并集，L_s是体网格下的拓扑拉普拉斯算子，它的构建方式与面网格下的类似，即：

${\left(L_s\right)}_{ij}=\left(\begin{array}{c}{d}_i,i=j\\ -1,(i,j)\in E\\ 0,otherwise\end{array}\right)---\left(6\right)$

L_s也被称为网格的拉普拉斯矩阵，能够表示网格的顶点之间的连接关系；d_i表示顶点i的度数，即与顶点i通过边相连的顶点的个数；(i,j)表示网格中由顶点i与顶点j构成的边，E为网格的边集合；通过对此优化问题的求解，能够得到每个操作点在每个顶点上的权重w_j(v_i)；并且，每个顶点上，所有操作点对其的权重和为1，每个操作点在每个顶点上的权重取值范围是[0,1]。

在步骤2.2)中，所述的边收缩方法，在保持形状特征方面，需要注意以下条件：

如果一条边的两个顶点为A、B，一个顶点A位于四面体的内部，而另一个顶点B位于四面体的外表面上，那么则应该拒绝进行边收缩的操作，否则将会使模型的外表面以顶点B为中心凹陷进去；

如果边AB不在网格模型的外表面上，但是边的两个顶点A、B均在网格的外表面上，那么说明顶点A、B分属于模型的两个不同的外表面区域，这种情况下应该拒绝进行边收缩的操作，否则将会使模型的外表面以这两个顶点为中心，向内凹陷成一个点，破坏模型形状的同时，也容易产生非法的拓扑单元；

如果想要进行边收缩的边的一个顶点在特征线的交汇处，那么另一个顶点不管是在网格的内部还是在其表面，都应该拒绝进行边收缩的操作，否则将会严重损毁CAD模型的形状特征。

在步骤2.3)中，所述的外表面四面体翻转操作，其方法为：

外表面的边与两个四面体相邻，直接进行四面体2-2翻转；

外表面的边与三或四个四面体相邻，这时候则能够先对相邻四面体中完全位于网格内部的边即指没有顶点位于网格外表面上的边执行边收缩操作，这样能够将其转换为当前边仅与两个四面体相邻的情况来操作，即利用四面体2-2翻转来处理；

外表面的边与大于四个四面体相邻，将这条边相邻的四面体单元看作一个整体，提取其与其他四面体相邻的三角面片以及外表面边翻转后的两个三角面片组成一个封闭的区域，利用TETGEN对此封闭区域进行四面体生成，删除原四面体后，用新生成的四面体填充这块区域，实现外表面边翻转的操作。

本发明与现有技术相比，具有如下优点与有益效果：

1、本发明算法可以快速、有效的编辑体网格CAD模型，并保持模型本身的语义形状特征。算法利用有限元网格的特点，把复杂的体拉普拉斯算子进行简化后建立相应的优化方程来编辑网格。由于该优化方程最后可以转化为线性系统的求解，因此，对于数量庞大的体网格模型来说，求解效率比较高，适合用户根据需要对局部形状特征进行实时的交互式的编辑。此外，本算法是通过变动区域的仿射变换来重构体网格顶点的，因此，非常适合应用于修改网格模型上的形状特征的参数，得到局部修改体网格模型的效果。

2、本发明算法通过定义在CAD体网格表面上的用于表示形状特征的“外表面区域”与“特征线”等概念，对位于网格外表面与内部的四面体，分别进行针对性的拓扑优化与几何优化操作。针对边收缩与四面体翻转操作提出了一些策略来保持网格的形状特征，并针对四面体网格外表面上的边，提出了边翻转的解决方案。实验表明，算法在保持模型形状特征的同时，对于低质量的四面体网格优化效果明显，并且可以快速的对CAD网格模型局部编辑后的区域进行体网格优化，使得优化后的网格模型能满足有限元分析的需求。

3、本发明算法可以直接对有限元网格模型进行局部编辑，并对编辑后质量较差的网格区域进行优化，并且整个过程中可以精确的保持模型原有的特征，有效解决实际工程中重复重网格化的问题，加快产品开发的速度，具有较高的实际意义。

附图说明

图1为目前工业产品设计流程示意图。

图2为本发明的基于特征的CAD体网格模型编辑流程图。

图3为本发明的基于特征的CAD四面体网格优化流程图。

图4为本发明的基于特征的CAD四面体网格优化中的边分裂操作示意图。

图5为本发明的基于特征的CAD四面体网格优化中的边收缩操作示意图。

图6为四面体2-3翻转与3-2翻转的示意图。

图7为四面体4-4翻转的示意图。

图8a为模型1编辑前的结果。

图8b为对模型1编辑后的结果。

图9a为模型1优化前内部的结果。

图9b为模型1优化后内部的结果。

图10为模型1优化前后，被优化的局部区域四面体质量分布图。

图11a为优化前的模型2内部的结果。

图11b为模型2优化后内部的结果。

图12为模型2优化前后，被优化的局部区域四面体质量分布图。

图13a为优化前的模型3结果。

图13b为模型3优化后的结果。

图14为模型3优化前后，被优化的局部区域四面体质量分布图。

具体实施方式

下面结合具体实施例对本发明作进一步说明。

本实施例所述的基于特征的CAD体网格模型编辑、优化算法，其具体情况如下：

1.基于特征的CAD体网格模型编辑

我们定义了相应的外表面区域来作为形状特征的几何面。每个外表面区域都是一组相邻的“表面网格单元”(对于四面体网格，这里的“表面网格单元”指的是三角形单元；而对六面体网格，“表面网格单元”则指的是四边形单元)的集合。

在模型的编辑过程中，用户可以指定某些外表面区域作为模型编辑的“操作区域”，位于操作区域上的所有顶点的重心，被称为该“操作区域”的“操作点”。对每个“操作点”构建局部坐标系，用户对不同的操作区域执行平移、旋转等仿射变换，该变换作用在局部坐标系上。通过对各个“操作点”的局部坐标的值，以及“操作点”对每个顶点的权重，来构建出网格每个顶点的坐标，实现模型的编辑操作。

如图2所示，基于特征的CAD体网格模型编辑步骤如下：

1)读入CAD体网格模型并对外表面区域进行分割，提取相应的外表面区域；

2)用户交互式地选取外表面区域，作为要编辑的形状特征的操作区域，在图2的过程2中，选择的两个操作区域用深色显示了出来；

3)对用户选定的每一个形状特征的“操作区域”，进行以下计算：

a)针对当前“操作区域”，建立局部坐标系，该“操作区域”对应的“操作点”作为当前局部坐标系的原点。

b)计算并记录网格所有顶点在当前局部坐标系中的坐标，便于通过操作区域、操作点来对整个网格进行编辑。

4)对每个“操作点”，构建优化方程，并利用Eigen求解优化方程，求得每个“操作点”对网格上每个顶点的权重。并对“操作点”的权重进行归一化操作，使得网格上的每个顶点，各个操作点对其的权重总和为1。

5)记录用户对每个“操作区域”进行的旋转平移等编辑操作，将相应的操作作用在以对应“操作点”为原点的局部坐标系的每个顶点坐标上。将每个“操作点”为原点的局部坐标系下的每个顶点的坐标数据变换到全局坐标系(最初模型所处的坐标系)，将每个顶点的坐标数据乘以其对应“操作点”的权重，加和得到最终顶点的坐标数据。这样，每个顶点的坐标都得以更新，模型的编辑操作也已经完成。

操作点权重的计算与编辑后的网格顶点坐标重构为这个算法的关键内容，下面对这两个步骤展开介绍。

1)操作点权重的计算

定义M＝(V,E,F,H)为具有n个顶点的四面体网格，V,E,F,H分别为该网格的顶点集合、边集合、面集合以及四面体集合，对于顶点i∈M，其笛卡尔坐标被表示为v_i＝(x_i,y_i,z_i)；定义R_j＝(V_j,E_j,F_j),j＝1,...,m为用户选取的所有的外表面操作区域，共m个，P_j,j＝1,...,m为每个操作区域所对应的操作点；w_j是操作点j在网格所有顶点上的权重集合，以向量的形式表示，w_j(v_j)则是操作点j对顶点i的权重；

由于Ls包含了网格的所有拓扑邻接信息，并且，操作点的权重，在当前操作区域所包含的顶点上应该为1，在非当前操作区域包含的顶点上，应该为0；由此，能够将求解操作点权重的问题建模为一个最优化问题，通过对这个最优化问题的求解来求得当前操作点对网格上所有顶点的权重，见下面的公式：

$\underset{w_j,j=1,\mathrm{...},m}{\mathrm{argmin}}\left({||L_s{w}_j||}^2\right)---\left(1\right)$

满足：

$w_j\left(v_j\right)=1,\forall i\in R_j---\left(2\right)$

$w_j\left(v_i\right)=0,\forall i\in (R-R_j)---\left(3\right)$

$\underset{j=1}{\overset{m}{\Sigma }}{w}_j\left(v_i\right)=1,\forall i\in M---\left(4\right)$

$0\le w_j\left(v_i\right)\le 1,\forall i\in M---\left(5\right)$

其中，R为用户选取的所有操作区域的集合，R_j＝(V_j,E_j,F_j),j＝1,...,m的并集，L_s是体网格下的拓扑拉普拉斯算子，它的构建方式与面网格下的类似，即：

${\left(L_s\right)}_{ij}=\left(\begin{array}{c}{d}_i,i=j\\ -1,(i,j)\in E\\ 0,otherwise\end{array}\right)---\left(6\right)$

L_s也被称为网格的拉普拉斯矩阵，能够表示网格的顶点之间的连接关系；d_i表示顶点i的度数，即与顶点i通过边相连的顶点的个数；(i,j)表示网格中由顶点i与顶点j构成的边，E为网格的边集合；通过对此优化问题的求解，能够得到每个操作点在每个顶点上的权重w_j(v_i)；并且，每个顶点上，所有操作点对其的权重和为1，每个操作点在每个顶点上的权重取值范围是[0,1]。

公式(1)-(5)可以通过最小二乘方法求解，公式(1)作为要求解的目标方程，公式(2)-(5)则作为约束。对算法利用C++实现时，可以通过Eigen库来求解这个优化问题。

2)编辑后的网格顶点坐标重构

利用步骤1)所介绍的优化问题的求解，我们得到了每个操作点对网格所有顶点的权重，通过这些权重，以及计算好的在每个操作点所记录的坐标系下的网格顶点坐标，可以重构出网格所有顶点的坐标数据。见公式(7)，表示操作点j所记录的坐标系下，点i的坐标值。

$v_i=\underset{j=1}{\overset{m}{\Sigma }}{w}_j\left(v_i\right)·{v_i}^{\left(j\right)},\forall i\in M---\left(7\right)$

当用户利用已选择的某个操作区域来对网格进行了某种仿射变换，定义为T_j,j＝1,...,m，可以将此仿射变换作用在对应操作点所记录的网格顶点坐标上，再通过公式(7)求得变换后的网格顶点坐标数据。如公式(8)所示：

$v_i^{\prime }=\underset{j=1}{\overset{m}{\Sigma }}{w}_j\left(v_i\right)·T_j{v}_i^{\left(j\right)},\forall i\in M---\left(8\right)$

v_i'表示进行网格模型编辑之后的网格顶点坐标数据，在实际的实验过程中，获得了来自于用户的对某个操作区域的仿射变换相关参数之后，将该变换作用在对应操作点的局部坐标系的网格顶点数据上，得到编辑后的局部坐标。然后，将局部坐标转换到全局坐标空间下，应用公式(8)即可获得编辑后的网格坐标数据，完成体网格模型编辑操作。

2.基于特征的CAD四面体网格模型优化

在对基于特征的体网格编辑进行介绍时，引入了体网格“外表面区域”的概念。通过对读入的CAD体网格模型进行预处理，获得网格的外表面区域集合。外表面区域可以作为表示网格表面形状特征的几何面，那么网格表面上不同外表面区域交界上的边就是辨识CAD网格形状特征的关键，这里将位于两个外表面区域交界上的边集合组成的空间中的线段集，称为“特征线”。本优化就是通过对特征线上的边与点进行特殊的处理，来实现对CAD体网格模型语义形状特征的保持。

如图3所示，基于特征的CAD四面体网格模型优化步骤如下：

1)如果优化算法目前执行的迭代次数小于迭代阈值(iteration_threshold)，则执行步骤2)，否则，执行步骤3)。

2)对于网格的每一条边，如果其边长大于4/5*target_edge_length，target_edge_length是目标边长，则对其执行边分裂(split)操作；如果其边长小于4/3*target_edge_length，那么对其执行边收缩(collapse)操作。本步骤的两种操作将会对网格模型的拓扑单元(顶点，线，面，四面体)进行增加与删除。这里的4/5与4/3两个参数的选取借鉴自Botsch等[Botsch M,Kobbelt L.A remeshing approach to multiresolution modeling[A].Proceedings of the 2004 Eurographics/ACM SIGGRAPH symposium on Geometry processing[C].ACM,2004:185-192.]所提出的三角网格的优化算法。执行步骤4)。

3)检查网格的每一条边与每一个面(仅在执行四面体2-3翻转时，以面为检查单位)，对于符合条件的边与面，对其邻接的四面体单元执行四面体翻转(flip)的操作。四面体翻转的具体条件以及需注意的问题将在下文中阐明。执行步骤4)。

4)执行网格顶点坐标的平滑(smooth)操作，过程如下：

a)对网格外表面上的顶点进行平滑，平滑时只根据这个顶点的一邻域中在外表面上的顶点来确定顶点平滑后的坐标。平滑过程中注意保持模型的语义形状特征，位于网格特征线上的顶点不应该进行平滑操作。另外，如果网格的外表面不是平面，那么表面顶点平滑后的位置往往不在外表面上，这就需要在平滑之后，对外表面上的顶点进行一个向原始网格(即输入网格)的外表面投影的操作。

b)对网格的内部顶点进行平滑。

5)如果算法执行了前四个步骤后，网格的拓扑元素数目相比于执行前变化很小，或者达到算法执行次数的上限，算法停止。否则，执行步骤1)。

优化算法的流程如上所述，目标边长的可以由用户指定，也可以按照以下公式求得，定义M＝(V,E,F,H)为具有n个顶点的四面体网格,V,E,F,H分别为该网格的顶点集合、边集合、面集合以及四面体集合，则目标边长的计算方法如下：

H_Q＝{h|h∈H^S_quality(h)＞τ}(9)

E_Q＝{e|e∈E_tetrahedro_n(h)^h∈H_Q}(10)

$t\arg et\_edge\_length=\frac{{\Sigma }_{e\in E_Q}{S}_{length}\left(e\right)}{\left|E_Q\right|}---\left(11\right)$

S_quality为计算四面体单元质量的函数，τ为四面体质量阈值，H_Q为网格中质量大于阈值的四面体集合。E_tetrahedron函数用来获得四面体的六条边的集合。E_Q为H_Q中所有四面体包含的边的集合。利用公式(11)计算出E_Q中所有边的平均长度，获得目标边长。S_length为计算边长度的函数，|E_Q|为E_Q集合中元素的数目。步骤2)、3)是对四面体网格进行拓扑优化操作，步骤4)是对网格模型进行几何平滑优化。步骤2)是通过边分裂与边收缩这两项拓扑操作来改变网格的四面体单元密度，使网格的边长快速趋近于目标边长。步骤3)的目的则是提高网格的质量，包括四面体网格外表面的规则程度以及内部四面体单元的质量。一般情况下，步骤2)执行3～8次后，网格的边长会与目标边长相差很小，可以不再进行边分裂与边收缩的操作，转而执行四面体翻转。所以，迭代阈值默认设为8，也可以由用户设定。对于网格局部压缩的编辑操作，步骤2)需要执行的次数较少，对于局部的拉伸的操作，需要执行步骤2)的次数比较多。步骤5)中，算法的执行次数的上限同样可由用户设定，默认设置为20。

本算法的关键步骤为拓扑优化操作与几何优化操作，下面展开详细介绍：

1)拓扑优化操作

算法执行的拓扑优化操作包括边分裂[Muthukrishnan S N,Shiakolas P S,Nambiar R V,et al.Simple Algorithm for Adaptive Refinement of Three-dimensional Finite Element Tetrahedral Meshes[J].AIAA Journal.1995,33(5):928-932.]，边收缩[De Cougny H L,Shephard M S.Parallel Refinement and Coarsening of Tetrahedral Meshes[J].International Journal for Numerical Methods in Engineering.1999,46(7):1101-1125.]，四面体翻转三类操作。四面体翻转由包括位于网格表面的四面体翻转与网格内部的四面体翻转。内部的四面体翻转[Klingner B M,Shewchuk J R.Aggressive tetrahedral mesh improvement[A].Proceedings of the 16th international meshing roundtable[C].Springer Berlin Heidelberg,2008:3-23.]主要包括四面体的“3-2翻转”(3-2flip)，“2-3翻转”(2-3flip)以及“4-4翻转”(4-4flip)。

a)边分裂

边分裂操作如图4所示，在执行时不会破坏网格原本的形状特征，也不会产生非法的四面体单元，没有太多需要注意的地方。

b)边收缩

边收缩操作如图5所示，在执行时要注意一下几点：

如果一条边的两个顶点为A、B，一个顶点A位于四面体的内部，而另一个顶点B位于四面体的外表面上，那么则应该拒绝进行边收缩的操作。否则将会使模型的外表面以顶点B为中心凹陷进去。

如果边AB不在网格模型的外表面上，但是边的两个顶点A、B均在网格的外表面上，那么说明顶点A、B分属于模型的两个不同的外表面区域，这种情况下应该拒绝进行边收缩的操作，否则将会使模型的外表面以这两个顶点为中心，向内凹陷成一个点，破坏模型形状的同时，也容易产生非法的拓扑单元。

如果想要进行边收缩的边的一个顶点在特征线的交汇处，那么另一个顶点不管是在网格的内部还是在其表面，都应该拒绝进行边收缩的操作，否则将会严重损毁CAD模型的形状特征。

假设要收缩的边为AB，集合S_a为与点A相邻并且不与边AB相邻的四面体集合，集合S_b为与点B相邻并且不与边AB相邻的四面体集合。如果集合S_a与集合S_b中的四面体两两之间有面相邻的情况，那么不能对AB执行边收缩的操作，否则会出现非法的四面体单元。

由于执行边收缩后，除了删除这条边相邻的一圈四面体之外，这条边两个顶点会被合并为一个顶点，这会引起附近的区域中未被删除的四面体的一个顶点位置发生变化。如果在变化后出现四面体翻转甚至交叉的情况，应该拒绝执行边收缩的操作。

c)外表面四面体翻转

外表面四面体翻转是指对位于网格外表面上的边的邻接四面体进行的翻转操作。目的是为了调整位于网格外表面上的顶点的度数，使得每个顶点的“外表面度数”(与当前顶点相邻的其他外表面顶点的数目)趋近于6，从而尽可能使本算法的四面体网格优化算法优化后的外表面网格与单纯的对外表面网格执行三角网格优化算法的结果一致。因此，这里的外表面四面体翻转前后的网格外表面的拓扑变化与现存的三角网格的边翻转操作实现的效果完全一致，这里也将这类操作称为“边翻转”。

不过由于体网格外表面上的边所邻接的四面体的数目并不固定，所以需要根据边所邻接的四面体的数目来对边翻转操作的实现进行分类，这也是本算法发明的创新点之一。

外表面的边与两个四面体相邻，这种情况也被称为四面体的“2-2翻转”(2-2 flip)，在这个操作的过程中，要翻转的边邻接的四面体所占有的空间体积没有被改变，而是被新的四面体单元所填充。两个四面体所共用的面以及共用的位于外表面上的边均被翻转，所以有些人也习惯将这里的“边翻转”称为“面翻转”。

外表面的边与三或四个四面体相邻，这时候则可以先对相邻四面体中完全位于网格内部的边(指没有顶点位于网格外表面上的边)执行边收缩操作，这样可以将其转换为当前边仅与两个四面体相邻的情况来操作。

外表面的边与大于四个四面体相邻，这种情况下如果再采取上一种情况提到的策略，连续的进行边收缩会对周边的四面体单元造成较大的修改，很可能产生更多低质量四面体，因此本算法对这种情况提出了新的策略。将这条边相邻的四面体单元看作一个整体，提取其与其他四面体相邻的三角面片以及外表面边翻转后的两个三角面片组成一个封闭的区域。利用TETGEN对此封闭区域进行四面体生成，删除原四面体后，用新生成的四面体填充这块区域，实现外表面边翻转的操作。

以上，对四面体网格边翻转的实现方式进行了阐述，不过在执行边翻转之前，还需要检测当前边是否应该进行边翻转，本算法按照以下的方式进行检查。

如果进行边翻转的边位于网格的特征线上，则拒绝对其进行边翻转的操作。否则，通过比较翻转前后的附近顶点的外表面度数来决定是否执行边翻转。以假设网格外表面上的边为AB，其两个顶点分别是A和B，如果要对其进行翻转，那么翻转后的边的两个顶点为C和D。设函数valence(v)可以获得顶点v的外表面度数，由于四面体网格模型的外表面必定是闭合的三角网格，那么每个位于外表面上的顶点的最佳外表面度数应该是6。如下式所示：

deviation_pre＝|valence(A)-6|+|valence(B)-6|+|valence(C)-6|+|valence(D)+6|deviation_post＝|valence(A)-1-6|+|valence(B)-1-6|+|valence(C)+1-6|+|valence(D)+1-6|(12)

在上式中，deviation_pre和deviation_post分别表示在进行边翻转前后，四个顶点的外表面度数与目标外表面度数(即6)之间的偏差总和。如果deviation_pre的值大于deviation_post，那么应该对这条边执行边翻转操作。

d)内部四面体翻转

图6为四面体2-3翻转与3-2翻转操作，图7展示了四面体4-4翻转操作。这些都属于内部四面体翻转的过程，这些操作不会破坏网格模型的形状特征，只要再四面体翻转的操作后，局部区域的四面体质量得到提高，并且没有产生非法的四面体单元，即可执行相应的四面体翻转操作。

2)几何优化操作

对于网格外表面上没有位于网格特征线上某一顶点，设其坐标为v，n为此顶点的法向量(可以通过对此顶点相邻的外表面三角面片的法向量加和再归一化得到)，p为此顶点邻接的外表面顶点坐标的质心，那么此顶点新的坐标值v'可以通过以下的公式得到，该公式借鉴自Botsch等[Botsch M,Kobbelt L.A remeshing approach to multiresolution modeling[A].Proceedings of the 2004 Eurographics/ACM SIGGRAPH symposium on Geometry processing[C].ACM,2004:185-192.]所提出的三角网格的优化算法。

v'＝p+nn^T(v-p)(13)

顶点的法向量n在这里的作用是保证平滑操作对顶点坐标的修改后，新的坐标会落在原顶点法向量的切平面上。

如果网格模型的外表面不是平面，在用以上公式平滑之后，外表面顶点新的位置很可能会偏离原模型的外表面，那么需要进行一个向原外表面投影的操作。投影操作的过程比较简单，对平滑后的外表面顶点，计算其法向量与原始网格的外表面的交点，此交点即为投影后的顶点坐标。

对于处于网格特征线上的顶点，为了简化操作，本算法的选择是在平滑操作时跳过这些顶点，在实际的实验中观察发现，这并不会对网格外表面的质量造成太大的影响。

内部顶点的平滑则更加简单一些。对于网格内部的某一顶点，设其坐标值为v，p为此顶点邻接的所有顶点坐标的质心，那么此顶点新的坐标值v'可以通过以下的公式得到：

v'＝(1-ω)v+ωp(14)

这里ω的取值介于0与1之间，本算法对其取值是0.75。

综上所述，本发明提供了基于特征的体网格编辑算法，图8a展示了该算法编辑前的网格模型。图8b展示了算法编辑后的网格模型结果。通过实验结果可以看出，算法可以快速、有效的编辑体网格CAD模型，并保持模型本身的语义形状特征。非常适合应用于修改网格模型上的形状特征的参数，得到局部修改体网格模型的效果。

本发明提供了基于特征的CAD四面体网格优化算法。图9a、图9b展示了模型1优化前后的网格内部结果。图11a、图11b分别展示了模型2优化前后的网格内部结果。图13a、图13b分别为模型3优化前后的网格外表面结果。图10、图12、图14分别展示了模型1、模型2、模型3优化前后，被优化的局部区域四面体质量分布图，图中的用两种颜色的柱状图来表示优化前后的质量分布，第一种颜色表示优化前的质量分布，第二种颜色表示优化后的质量分布。实验表明，算法在保持模型形状特征的同时，对于低质量的四面体网格优化效果明显，并且可以快速的对CAD网格模型局部编辑后的区域进行体网格优化，使得优化后的网格模型能满足有限元分析的需求。

以上所述之实施例子只为本发明之较佳实施例，并非以此限制本发明的实施范围，故凡依本发明之形状、原理所作的变化，均应涵盖在本发明的保护范围内。