CN113987856A - 一种基于标架场的复杂多约束结构网格生成方法 - Google Patents

Abstract

本发明公开了针对复杂几何特征约束及尺寸约束的模型的一种基于标架场的复杂多约束结构网格生成方法。现有四边形网格生成方法很难在达到在满足复杂约束的同时存在高鲁棒性的方法并保证较高质量的网格结果。本发明采用的标架场指导网格方法一方面有着较高质量的保证，此外，标架场仅在初始流线的生成扮演着一定作用，后续方法均脱离标架场，避免了标架场的求解不适用网格生成的情况，此方法在测试的多种复杂约束情况下，网格在保证几何约束，并满足尺寸约束的情况下能够以较高质量的结果生成。

Description

一种基于标架场的复杂多约束结构网格生成方法

技术领域

本发明属于计算机辅助设计领域，具体涉及一种针对CAD模型的准结构化四边形网格生成方法。

背景技术

关于结构化四边形网格的生成，大体上有两个主要的领域在进行研究，其一：工程仿真分析，四边形网格相较于三角网格在在有限元和有限体积方法方面有着更加良好的几何支持，其二：计算机图形领域，除了数值仿真方面，四边形网格在曲面建模和纹理中也有很好的用途。

近几十年来，工程分析领域一直致力于四边网格的生成，产生了两大类工业技术，其一：利用一些半自动辅助的方式手动划分四边形区域，然后根据用户需要生成规范的四边形网格填充区域，这样的方式非常耗时，主要用于苛刻的数值模拟，如CFD。其二：有一些全自动的方式如基于标架场的前沿推进法生成直角三角网格，再进行合并细分等操作生成全四边形网格，然而这样的方法存在大量的不规则顶点，缺乏用户控制。

另一方面，计算机图形学领域彻底开发和探索了基于标架场的四边形网格划分技术，核心是自动生成四边形拓扑分区，然而这样的拓扑分区并不一定确保能够提取出来，这样的技术应用于工程分析方面缺乏足够的鲁棒性和稳定性。

本文提出的准四边形网格划分基于标架场划分出拓扑分区，却又不要求其一定是四边形的拓扑分区，在这样的拓扑分区中应用模板方法对每个分区生成网格，鲁棒性得到了一定的保证，经测试，也基本可以满足用户的指定尺寸场需求，且生成的网格奇异点少，质量较高。

发明内容

本发明的目的是针对现有技术的不足，提出一种针对CAD模型的基于有限元标架场计算的复杂多约束四边形网格结构生成方法，根据复杂几何特征约束及尺寸约束情况，对网格中的标架场进行拓扑提取，进行合理的拓扑简化，使用模板法对拓扑生成网格。

本发明具体如下：

步骤1、首先将输入的背景三角网格及对应的尺寸约束与几何特征约束中的尺寸约束转化为相应的几何约束；

1.1对每个三角形面使用marching triangle法提取等值线；

1.2连接提取的等值线，并进行一定的光滑；

1.3过滤数量少的等值线；

1.4再等值线处进行重新网格化，并将等值线作为几何特征约束存储；

步骤2、使用有限元方法计算背景网格的标架场并提取其中的奇异点；

交叉表示：

使用标准单位向量表示u＝(cos(4θ),sin(4θ))＝(u₁,u₂)，其中θ是四个分支之一的角度与局部参考范围。与一般的角度表示θ相反，向量表示u对四边形对称性是不变的，适用于线性有限元插值。由于单位向量的线性组合通常不是单位向量，三角形内的插值并不是严格的十字表示，但是通过投影到单位圆上可以很容易地恢复十字，即对向量表示进行归一化。

有限元离散化：

为了离散交叉场的矢量场表示，在三角剖分的每条边上定义一个交叉，并使用Crouzeix-Raviart插值，如中所示。交叉域问题的未知数是每条边eij处的表示向量场分量u1、u2。由于希望交叉场与表面边界及特征线对齐，因此Dirichlet边界条件在每个边界与特征线处为u＝(cos(0),sin(0))＝(1,0)。

利用热传导方程对标架场进行平滑

为了在域内有一个平滑的交叉场，一个自然的方法是最小化狄利克雷能量：

由于

这个目标函数是非线性的。通过直接求解向量表示的拉普拉斯算子(即

)来最小化狄利克雷能量显然是不合适的，因为这些值可能会在远离边界的地方坍缩为0。这个问题的解决方案是求解一个Ginzburg-Landau非线性方程来惩罚离开交叉流形的值，但该方法对于复杂问题来说过于耗时。一种有效的替代方法是MBO方法，其中交替进行热扩散(方程1)和投影(方程2)步骤。本方法使用递减扩散时间步长，并根据网格大小的考虑选择它们，并且将它们重新分组以允许重新使用由直接线性求解器计算的矩阵分解。

扩散：

投影：

检测奇异点：

对于表面上的每个点x，通过计算沿以x为中心的小闭合圆γ的角度差来检测交叉场奇点。这定义了一个奇点指数：

根据定义，交叉场奇点的索引an为-1时，它对应于四边形网格中索引k＝-1的不规则顶点，即价五。

在离散设置中，提取顶点i周围的定向边缘单环(ei1,...,ein)并计算角度差的总和。在实践中存在三种情况：和为零且顶点是规则的或和为一或负一且顶点为标架场域的奇点。由于对标架场的Crouzeix-Raviart离散化(每条边一个角度)，交叉场奇点可能位于网格的顶点、边或三角形上。当奇点位于边上时，它的两个相邻顶点的索引都不为零，而当奇点位于三角形上时，三个相邻顶点的索引都不为零。

步骤3、在奇异点处及计算出的部分边界点处实行龙格库塔法进行流线的延伸；

3.1、首先在上述点处的一领域及二领域处计算起始点的第一个离散点：

①首先在边界奇异点的相关领域计算出第一个离散点

②在边界奇异点处计算是否存在内部的奇异点，如有执行③，否则跳至④；

③如果奇异点在面内，将此离散点延伸至三角形的边/点上，否则直接作为一个离散点，而后删除此奇异点；

④将上述得到的所有与边界奇异点相关的离散点进行筛选，去除方向接近的离散点

⑤计算剩下的内部奇异点的相关领域的第一个离散点

⑥如果存在奇异点的二领域内有其他奇异点，如有执行⑦，否则跳至⑧，

⑦挑选一个奇异点，将另一奇异点作为此奇异点的一个离散点，如果此离散点在面内，延伸至三角形的边/点，否则直接作为一个离散点，而后删除另一个奇异点；

⑧将上述得到的所有与内部奇异点相关的离散点进行筛选，去除方向接近的离散点

3.2、根据龙格库塔方法将得到的奇异点与相应的第一个离散点进行流线的延伸，终止条件如下所述：

①直到奇异点的一领域，以奇异点终止

②以约束边/边界边的终止处的最近的三角网格点终止；

③同一条流线以同样的方向经过同一个背景网格三角形，终止。

步骤4、简化流线并构造拓扑关系：

4.1对于步骤3中得到的流线结构，检测其流线情况，针对不同的情况进行不同的简化，直到没有一种可简化的情况出现，所有简化情况大致分为三种情形：

①两条流线起止点相同，分别记为A,B，且长度相差较大，改长的流线为A->A,短流线不变

②两条流线起止点相同，且长度相差不大，去除其中一条流线

③两条流线起止点有一个相同记为A，额外的点分别为B,C，不失一般性，设A->B，长于A->C，改长流线为B->C；

4.2根据最终的流线结构构造拓扑分区关。

步骤5、弦折叠的拓扑简化，基于四边形网格对偶结构以及保持网格中奇异点的数量和基本拓扑关系不变从而设计一种四边形网格拓扑简化操作与优化框架，算法执行效率高，简化效果好，可适用于各类四边形网格的简化。

弦、横边、纵边的定义如下：

弦：在四边形区域的拓扑划分中根据每条边所相对的对边的集合构成的四边区域。

横边：弦内从边界开始到另一条边界对边中所有的对边。

纵边：弦内除横边外的其余边(连接横边顶点的边)。

弦的折叠操作：

从模型的对偶表示中删掉一整条弦，意味着删除这根弦所经过的所有四边形面片。具体的删除操作是通过折叠弦经过的所有横边完成，也就是说将弦横边的顶点进行两两合并，对应的纵边上的区域线两两合并。四边形网格中的对偶弦折叠不会影响网格的连通性，并且可以保证折叠得到的网格依然是纯四边形网格。

四边形区域的折叠操作：

四边形区域的折叠就是把弦内部分四边形区域组合而成的多块区域相对的两个顶点进行合并而实现的拓扑区域删除操作。相对于对偶弦的折叠来说，这是个局部性的网格简化操作。可以将将四边区域想象而成两个三角形，要合并的顶点之间的边是公共边，而四边形的折叠可以看成是三角形的边折叠。然而这条公共边由弦的两条纵边进行线性插值合成。由于纵边不是简单一条直线，是由多条区域线组合而成并且边上的点的数量不一致，因此插值的过程还需要对每条区域线取相同点的数量进行拟合，对应计算公式如下：

q_i＝(1-t_i)×a_i+t_i×b_i,q_i+1＝(1-t_i+1)×a_i+1+t_i+1×b_i+1

式中i表示区域线的序号，a_ia_i+1,b_ib_i+1为弦内纵边上的区域线，a_ib_i为弦内的横边，区域线由许多折线构成，j表示区域线上的折线序号,t_i为

x_i,j为纵边区域线上的点

四边形区域的折叠操作影响范围仅仅是于被折叠四边形相连的区域，并且整个操作满足不引入非四边区域的限制。

步骤5包括如下子步骤：

①：在四边形区域的拓扑划分中根据对边划分弦(从边界边一直找对边直至另一条边界边或者从该边出发绕一圈回到该边)

②：将弦按类型分类：

(1)无奇异点存在或者奇异点存在于同一条纵边

(2)同一条横边上存在两个奇异点

(3)两条纵边均存在奇异点并且不在同一条横边上

③：对3类弦进行检查：

(1)如果存在被合并的纵边为边界边或者约束边时，抛弃该弦。

(2)第一类的弦不需进行处理，第二、三类的弦转至步骤四。

④：根据优先级度量的能量公式对选取的弦进行排序，选择能量最大且大于0的的弦进行拓扑区域简化。如果不存在能量大于0的弦，则结束简化。

优先级度量包括：

使用能量公式选择删除那些狭长且不规整的四边区域，将一个四边区域拓扑剖分结构简化为奇异点最少的最优拓扑结构；能量公式如下：

式中K控制弦能量比值；

⑤：检查弦的类型进行区域简化：

(1)第二类型的弦，直接进行四边形区域的折叠操作；

(2)第三类型的弦将根据奇异点的位置划分为non-zip-patch区域与zip-patch区域(zip-patch为相邻相对角的奇异点之间的区域)

non-zip-patch区域删除离其最近的奇异点的弦内对边纵边和点，然后与该纵边外的四边区域合并zip-patch区域根据原来的横边和纵边线性插值拟合出一条新的纵边，删除原先两条纵边，再将横边根据新点一分为二与弦外边进行合并。

(3)返回①。

步骤6、针对上述得到的区域划分，以背景尺寸长计算每个区域边的分段数，并对此段数进行一定的优化：

由于使用模板法生成网格依赖于每条边的段数设置，因此区域线量化很重要。首先对模型区域线进行网格剖分，然后对模型区域面进行模板生成，受限于曲线离散化，这对于生成三角剖分(或四面体剖分)效果很好，但对于四边形剖分(或六面体剖分)作为四边形拓扑约束(弦)是全局的并通过模型曲线。

对区域线的量化非常简单。使用通过整合全局尺寸场计算的理想边数的整数舍入值，除了四边形区域的区域线，在对边施加相等性。这种非最佳量化会导致不同区域面的最终四边形网格中出现许多必要的偶极子。

步骤6包括如下子步骤：

子步骤1，使用量化计算方式，对所有的区域线结合尺寸场计算理想边数的整数舍入值；所述的区域线上的分段数量化计算方式：

子步骤2，使用四边区域量化的拓扑约束，将一些四边区域的分段数施加相等性；

子步骤3，得到的最终的理想边数对所有的区域线进行二次量化。

区域线上的分段数量化计算

考虑一条区域线C，该区域线C由从a到b变化的t参数化。曲线上根据尺寸图的浮点边数为：

为了得到整数个边n_e，使用一个简单的舍入：

与区域线的第i个内点相关联的参数t_i,i∈[1,n_e-1]是这样的：

这个积分方程可以通过数值积分求解，通过沿区域线添加值直到总和等于

在样本之间进行线性插值。为了对区域线C进行网格划分，通过使用(等式11)计算参数t_i并评估区域线参数化：x_i＝f(t_i)来计算顶点位置x_i,i∈[1,n_e-1]。使用这种方法，根据尺寸图将点很好地放置在曲线上，从较小的特征区域平滑过渡到较粗的区域。

四边形区域的拓扑约束

整合尺寸场没有考虑四边形网格拓扑的特殊性，这样四边形被组织成拓扑和弦(相邻元素的双环)。在这项工作中，调整了四边形区域面相对两侧的边数n_e。考虑一个边界由四条区域线(C₁、C₂、C₃、C₄)组成的区域面，强制对边上的点数相等，即n_e1＝n_e3和n_e2＝n_e4，除非积分值非常不同。当两个相邻的四边形区域面共享一条区域线C_c时，两个面上的n_ec值必须相同。这意味着等式约束在传播。为了解决传播问题，构建了与四边形区域面相关的拓扑弦。弦由拓扑平行的区域线组成，它们都接收相同的固定点数，并通过对先前在每条曲线上计算的理想值求平均值来计算这些点。通过这种简单的传播，弦仅在两个相邻的四边形面共享一条曲线时才传播，该曲线是它们的四个边之一。

步骤7、使用模板法对每个分区进行网格生成，步骤：

1)按照顺/逆时针输入多边形区域角点坐标和相应边的细分数。

2)根据输入的细分数选择符合的拓扑模板，通过计算相应参数确定patch的拓扑结构。

3)计算边界上点的坐标，使用拉普拉斯光顺确定patch的几何位置。

subject to v_i＝ω_i,i∈C

其中，ε为网格的边集，C为边界点集，ω_i∈R²为输入第i个固定边界点坐标。

简化输入：

由于考虑所有可能的输入情况是一项挑战，所以可以考虑将问题简化到一个等效的，更容易解决的子问题，此操作大大减少了需要考虑的情况，以此确保算法具有通用性。

对于输入(l₀,...,l_N-1)，其对于某些k，l_k-1和l_k+1均大于1，定义d＝min(l_k-1,l_k+1)-1，由此输入可以简化为(l₀,...,l_N-1)，其中

可以重复上述的简化操作，直到无法再进行简化为止，从而获得了最大程度简化的输入。

模板选择公式化为ILP：

每个拓扑模板都定义了输入(l₀,...,l_N-1)与参数之间的线性关系，这些参数包含第i条边的填充量p_i、边缘流个数x和y。

上述的通用形式记为：

Ax＝b

其中，A是一个N×M的矩阵，其中M为该模板的参数个数，x表示模板参数的M维向量(待求解)，b是由输入(l₀,...,l_N-1)和模板共同决定的N维向量，并要求

为非负整数，可以将其表达为整数线性规划(ILP)：

argmax c^Tx

subject to Ax＝b

x≥0,x∈Z^M

其中，c代表希望最大化的目标的M维向量，只要该ILP有解，则意味着该模板对于输入而言是可行的。由于M即某个模板的参数个数很少(最多为10个)，因此可以很快求解。由于N≤M，所以上述ILP通常有多解，这就意味着，对于同一个拓扑模板，有多种参数组合满足输入要求，将目标设为最大化

在此目标下，当对应的变量代表边界填充量时，c_j＝1；否则，c_j＝0，此时的拓扑生成包含的奇异点个数是最少的。

给定问题域的区域划分结果以及各个子区域相应边的细分数，各个子区域四边形网格生成使用上述算法求得，然后将其拼接为最终的网格效果。

步骤8、网格光顺：

①规则顶点：给定一个固定的边界，求解Winslow非线性椭圆PDEΔ_xu＝0，其中u(x)是某个坐标计算空间，x是物理空间中的坐标。这样做的好处是两个坐标分量是耦合的，结果四边形形状很好，即使在大的扭曲下也能强制执行一些正交性。通过应用有限差分(FDM)离散化到Winslow方程，可以推导出局部平滑四边形网格中规则顶点的核。假设(x1,...,x8)是有序的围绕规则顶点x的模板的顶点，其新位置由下式给出：

其中α₀＝(x₃-x₇)·(x₃-x₇)，α₁＝(x₁-x₅)·(x₁-x₅)，β＝(x₁-x₅)·(x₃-x₇)

②不规则顶点：使用基于角度的平滑方法，其想法是移动单环平分线上的顶点。作为上面的Winslow内核，它在网格不太受限制的情况下效果很好，但可能会在复杂配置中创建翻转单元。

本发明的实质性效果在于：鲁棒性得到了保证，满足用户的指定尺寸场需求，且生成的网格奇异点少，质量较高。

附图说明

图1为本发明的流程图；

图2-1为一个模型的三角性背景网格图；

图2-2为图2-1模型几何特征线约束图；

图2-3为图2-1模型背景尺寸约束图；

图3为标架场计算结果及奇异点分布示意图；

图4为初始流线结果示意图；

图5为简化流线后的拓扑构造示意图；

图6为简化后的拓扑结构示意图；

图7为根据最终拓扑生成的初始的背景四边形网格示意图；

图8为网格几何优化后的最终四边网格示意图。

具体实施方式

下面结合附图对本发明进一步说明。

如图1所示，针对复杂约束模型的四边形网格生成方法，具体步骤如下：

步骤1、首先将输入的背景三角网格及对应的尺寸约束与几何特征约束中的尺寸约束转化为相应的几何约束；

1.1对每个三角形面使用marching triangle法提取等值线；

1.2连接提取的等值线，并进行一定的光滑；

1.3过滤数量少的等值线；

1.4再等值线处进行重新网格化，并将等值线作为几何特征约束存储；

值得一提的是，此步骤在几何特征约束能够将网格尺寸变化大的地方包括的话，并不需要特地再加这样的等值线，在此例子中，未加入等值线。图2-1为背景三角网格示意图，图2-2为几何特征线约束示意图，图2-3为背景尺寸场约束。

步骤2、使用有限元方法计算背景网格的标架场并提取其中的奇异点；

交叉表示：

使用标准单位向量表示u＝(cos(4θ),sin(4θ))＝(u₁,u₂)，其中θ是四个分支之一的角度与局部参考范围。与一般的角度表示θ相反，向量表示u对四边形对称性是不变的，适用于线性有限元插值。由于单位向量的线性组合通常不是单位向量，三角形内的插值并不是严格的十字表示，但是通过投影到单位圆上可以很容易地恢复十字，即对向量表示进行归一化。由于希望交叉场与表面边界及特征线对齐，因此Dirichlet边界条件在每个边界与特征线处为u＝(cos(0),sin(0))＝(1,0)。

用热方程平滑

为了在域内有一个平滑的交叉场，自然的方法是最小化狄利克雷能量：

在这项工作中，使用递减扩散时间步长，根据网格大小的考虑选择它们，并且将它们重新分组以允许重新使用由直接线性求解器计算的矩阵分解。

扩散：

投影：

检测奇异点：

对于表面上的每个点x，通过计算沿以x为中心的小闭合圆γ的角度差来检测交叉场奇点。这定义了一个奇点指数：

根据定义，交叉场奇点的索引an为-1时，它对应于四边形网格中索引k＝-1的不规则顶点，即价五。

在离散设置中，提取顶点i周围的定向边缘单环(ei1,...,ein)并计算角度差的总和。在实践中存在三种情况：和为零且顶点是规则的或和为一或负一且顶点为标架场域的奇点。由于对标架场的Crouzeix-Raviart离散化(每条边一个角度)，交叉场奇点可能位于网格的顶点、边或三角形上。当奇点位于边上时，它的两个相邻顶点的索引都不为零，而当奇点位于三角形上时，三个相邻顶点的索引都不为零。

如图3,为从图2-1的背景网格及2-2中的几何特征约束计算出的标架场及奇异点信息。

步骤3、在奇异点处及计算出的部分边界点处实行龙格库塔法进行流线的延伸；

3.1、首先在上述点处的一领域及二领域处计算起始点的第一个离散点：

①首先在边界奇异点的相关领域计算出第一个离散点

②在边界奇异点处计算是否存在内部的奇异点，如有执行③，否则跳至④；

③如果奇异点在面内，将此离散点延伸至三角形的边/点上，否则直接作为一个离散点，而后删除此奇异点；

④将上述得到的所有与边界奇异点相关的离散点进行筛选，去除方向接近的离散点

⑤计算剩下的内部奇异点的相关领域的第一个离散点

⑥如果存在奇异点的二领域内有其他奇异点，如有执行⑦，否则跳至⑧，

⑦挑选一个奇异点，将另一奇异点作为此奇异点的一个离散点，如果此离散点在面内，延伸至三角形的边/点，否则直接作为一个离散点，而后删除另一个奇异点；

⑧将上述得到的所有与内部奇异点相关的离散点进行筛选，去除方向接近的离散点

3.2、根据龙格库塔方法将得到的奇异点与相应的第一个离散点进行流线的延伸，终止条件如下所述：

①直到奇异点的一领域，以奇异点终止

②以约束边/边界边的终止处的最近的三角网格点终止；

③同一条流线以同样的方向经过同一个背景网格三角形，终止；

图4为初始流线的示意图。

步骤4、简化流线结构，构造拓扑分区

4.1对于步骤3中得到的流线结构，检测其流线情况，针对不同的情况进行不同的简化，直到没有一种可简化的情况出现，所有简化情况大致分为三种情形：

①两条流线起止点相同，分别记为A,B，且长度相差较大，改长的流线为A->A,短流线不变

②两条流线起止点相同，且长度相差不大，去除其中一条流线

③两条流线起止点有一个相同记为A，额外的点分别为B,C，不失一般性，设A->B，长于A->C，改长流线为B->C；

4.2求出流线与流线之间的交点，并构造拓扑关系。

图5为初始拓扑分区的示意图。

步骤5、弦折叠的拓扑简化，基于四边形网格对偶结构以及保持网格中奇异点的数量和基本拓扑关系不变从而设计一种四边形网格拓扑简化操作与优化框架，算法执行效率高，简化效果好，可适用于各类四边形网格的简化。

弦、横边、纵边的定义如下：

弦：在四边形区域的拓扑划分中根据每条边所相对的对边的集合构成的四边区域。

横边：弦内从边界开始到另一条边界对边中所有的对边。

纵边：弦内除横边外的其余边(连接横边顶点的边)。

弦的折叠操作：

从模型的对偶表示中删掉一整条弦，意味着删除这根弦所经过的所有四边形面片。具体的删除操作是通过折叠弦经过的所有横边完成，也就是说将弦横边的顶点进行两两合并，对应的纵边上的区域线两两合并。四边形网格中的对偶弦折叠不会影响网格的连通性，并且可以保证折叠得到的网格依然是纯四边形网格。

四边形区域的折叠操作：

四边形区域的折叠就是把弦内部分四边形区域组合而成的多块区域相对的两个顶点进行合并而实现的拓扑区域删除操作。相对于对偶弦的折叠来说，这是个局部性的网格简化操作。可以将将四边区域想象而成两个三角形，要合并的顶点之间的边是公共边，而四边形的折叠可以看成是三角形的边折叠。然而这条公共边由弦的两条纵边进行线性插值合成。由于纵边不是简单一条直线，是由多条区域线组合而成并且边上的点的数量不一致，因此插值的过程还需要对每条区域线取相同点的数量进行拟合，对应计算公式如下：

q_i＝(1-t_i)×a_i+t_i×b_i,q_i+1＝(1-t_i+1)×a_i+1+t_i+1×b_i+1

式中i表示区域线的序号，a_ia_i+1,b_ib_i+1为弦内纵边上的区域线，a_ib_i为弦内的横边，区域线由许多折线构成，j表示区域线上的折线序号,t_i为

x_i,j为纵边区域线上的点

四边形区域的折叠操作影响范围仅仅是于被折叠四边形相连的区域，并且整个操作满足不引入非四边区域的限制。

优先级度量：

为了提高模型的连接性，保证其几何保真度，使用一个可以智能选择删除元素的算法尤其重要。折叠操作的优先级是通过对原始拓扑区域中的弦进行能量公式计算来排序得到的，尽可能的删除那些狭长且不规整的四边区域，将一个四边区域拓扑剖分结构简化为奇异点最少的最优拓扑结构。能量公式如下：

式中K控制弦能量比值，本算法例子中设置为4。

算法具体步骤如下：

①：在四边形区域的拓扑划分中根据对边划分弦(从边界边一直找对边直至另一条边界边或者从该边出发绕一圈回到该边)

②：将弦按类型分类：

(1)无奇异点存在或者奇异点存在于同一条纵边

(2)同一条横边上存在两个奇异点

(3)两条纵边均存在奇异点并且不在同一条横边上

③：对3类弦进行检查：

(1)如果存在被合并的纵边为边界边或者约束边时，抛弃该弦。

(2)第一类的弦不需进行处理，第二、三类的弦转至步骤四。

④：根据优先级度量的能量公式对选取的弦进行排序，选择能量最大且大于0的的弦进行拓扑区域简化。如果不存在能量大于0的弦，则结束简化。

⑤：检查弦的类型进行区域简化：

(1)第二类型的弦，直接进行四边形区域的折叠操作

(2)第三类型的弦将根据奇异点的位置划分为non-zip-patch区域与zip-patch区域(zip-patch为相邻相对角的奇异点之间的区域)

non-zip-patch区域删除离其最近的奇异点的弦内对边纵边和点，然后与该纵边外的四边区域合并zip-patch区域根据原来的横边和纵边线性插值拟合出一条新的纵边，删除原先两条纵边，再将横边根据新点一分为二与弦外边进行合并。

(3)返回①

图6为简化后的拓扑分区示意图。

步骤6、针对上述得到的区域划分，以背景尺寸长计算每个区域边的分段数，并对此段数进行一定的优化：

由于使用模板法生成网格依赖于每条边的段数设置，因此区域线量化很重要。首先对模型区域线进行网格剖分，然后对模型区域面进行模板生成，受限于曲线离散化，这对于生成三角剖分(或四面体剖分)效果很好，但对于四边形剖分(或六面体剖分)作为四边形拓扑约束(弦)是全局的并通过模型曲线。例如，想象一个简单的矩形区域，根据连续尺寸场，其四边的理想边缘数为[9.4,4.1,10.6,3.9]。通过整数舍入的简单量化将分别选择[9,4,11,4]边。虽然存在具有这种量化的四边形网格，但它不是结构化的，因为它必须包括一对+1/–1不规则顶点，以便从一侧的9个边过渡到另一侧的11个边。

对区域线的量化非常简单。使用通过整合全局尺寸场计算的理想边数的整数舍入值，除了四边形区域的区域线，在对边施加相等性。这种非最佳量化会导致不同区域面的最终四边形网格中出现许多必要的偶极子(一对+1/–1不规则顶点)。

区域线上的分段数量化计算

考虑一条区域线C，该区域线C由从a到b变化的t参数化。曲线上根据尺寸图的浮点边数为：

其中dl＝||x'(t)||dt。为了得到整数个边n_e，使用一个简单的舍入：

与区域线的第i个内点相关联的参数t_i,i∈[1,n_e-1]是这样的：

这个积分方程可以通过数值积分求解，通过沿区域线添加值直到总和等于

在样本之间进行线性插值。为了对区域线C进行网格划分，通过使用(等式11)计算参数t_i并评估区域线参数化：x_i＝f(t_i)来计算顶点位置x_i,i∈[1,n_e-1]。使用这种方法，根据尺寸图将点很好地放置在曲线上，从较小的特征区域平滑过渡到较粗的区域。

四边形区域的拓扑约束

整合尺寸场没有考虑四边形网格拓扑的特殊性，这样四边形被组织成拓扑和弦(相邻元素的双环)。在这项工作中，调整了四边形区域面相对两侧的边数n_e。考虑一个边界由四条区域线(C₁、C₂、C₃、C₄)组成的区域面，强制对边上的点数相等，即n_e1＝n_e3和n_e2＝n_e4，除非积分值非常不同。当两个相邻的四边形区域面共享一条区域线C_c时，两个面上的n_ec值必须相同。这意味着等式约束在传播。为了解决传播问题，构建了与四边形区域面相关的拓扑弦。弦由拓扑平行的区域线组成，它们都接收相同的固定点数，通过对先前在每条曲线上计算的理想值求平均值来计算这些点。通过这种简单的传播，弦仅在两个相邻的四边形面共享一条曲线时才传播，该曲线是它们的四个边之一。

步骤7、使用模板法对每个分区进行网格生成

步骤：

1)按照顺/逆时针输入多边形区域角点坐标和相应边的细分数。

2)根据输入的细分数选择符合的拓扑模板，通过计算相应参数确定patch的拓扑结构。

3)计算边界上点的坐标，使用拉普拉斯光顺确定patch的几何位置。

subject to v_i＝ω_i,i∈C

其中，ε为网格的边集，C为边界点集，ω_i∈R²为输入第i个固定边界点坐标。

简化输入

由于考虑所有可能的输入情况是一项挑战，所以可考虑将问题简化到一个等效的，更容易解决的子问题，此操作大大减少了需要考虑的情况，以此确保算法具有通用性。

对于输入(l₀,...,l_N-1)，其对于某些k，l_k-1和l_k+1均大于1，定义d＝min(l_k-1,l_k+1)-1，由此输入可以简化为(l₀,...,l_N-1)，其中

可以重复上述的简化操作，直到无法再进行简化为止，从而获得了最大程度简化的输入。

模板选择公式化为ILP

每个拓扑模板都定义了输入(l₀,...,l_N-1)与参数之间的线性关系，这些参数包含第i条边的填充量p_i、边缘流个数x和y。

上述的通用形式记为：

Ax＝b

其中，A是一个N×M的矩阵，其中M为该模板的参数个数，x表示模板参数的M维向量(待求解)，b是由输入(l₀,...,l_N-1)和模板共同决定的N维向量，并要求

为非负整数，可以将其表达为整数线性规划(ILP)：

argmax c^Tx

subject to Ax＝b

x≥0,x∈Z^M

其中，c代表希望最大化的目标的M维向量，只要该ILP有解，则意味着该模板对于输入而言是可行的。由于M即某个模板的参数个数很少(最多为10个)，因此可以很快求解。由于N≤M，所以上述ILP通常有多解，这就意味着，对于同一个拓扑模板，有多种参数组合满足输入要求，将目标设为最大化

在此目标下，当对应的变量代表边界填充量时，c_j＝1；否则，c_j＝0，此时的拓扑生成包含的奇异点个数是最少的。

给定问题域的区域划分结果以及各个子区域相应边的细分数，各个子区域四边形网格生成使用上述算法求得，然后将其拼接为最终的网格效果。

图7为根据步骤7所计算的分区点进行模板法生成网格的结果示意图。

步骤8、网格光顺

①规则顶点：给定一个固定的边界，求解Winslow非线性椭圆PDEΔ_xu＝0，其中u(x)是某个坐标计算空间，x是物理空间中的坐标。这样做的好处是两个坐标分量是耦合的，结果四边形形状很好，即使在大的扭曲下也能强制执行一些正交性。通过应用有限差分(FDM)离散化到Winslow方程，可以推导出局部平滑四边形网格中规则顶点的核。假设(x1,...,x8)是有序的围绕规则顶点x的模板的顶点，其新位置由下式给出：

其中α₀＝(x₃-x₇)·(x₃-x₇)，α₁＝(x₁-x₅)·(x₁-x₅)，β＝(x₁-x₅)·(x₃-x₇)

②不规则顶点：使用基于角度的平滑方法，其想法是移动单环平分线上的顶点。作为上面的Winslow内核，它在网格不太受限制的情况下效果很好，但可能会在复杂配置中创建翻转单元。

图8为对步骤7的网格结果进行几何优化的最终网格示意图。

综上，本发明通过将标架场提取简单拓扑与模板生成网格相结合，得到了满足复杂几何特征约束及尺寸约束的高精度、高质量以及有效性要求的四边形网格，从而可以应用到高层次CAD、CAE模型的网格化研究设计中。

Claims (1)

1.一种基于标架场的复杂多约束结构网格生成方法，其特征在于，包括以下步骤：

步骤1、将输入的背景三角网格及对应的尺寸约束与几何特征约束中的尺寸约束转化为相应的几何约束；

1.1对每个三角形面使用marching triangle法提取等值线；

1.2连接提取的等值线，并进行一定的光滑；

1.3过滤数量少的等值线；

1.4再等值线处进行重新网格化，并将等值线作为几何特征约束存储；

步骤2、使用有限元方法计算背景网格的标架场并提取其中的奇异点；

步骤3、在奇异点处及计算出的部分边界点处实行龙格库塔法进行流线的延伸；

3.1、在奇异点和计算出的部分边界点处的一领域及二领域处计算起始点的第一个离散点：

①在边界奇异点的相关领域计算出第一个离散点；

②在边界奇异点处计算是否存在内部的奇异点，如有执行③，否则跳至④；

③如果奇异点在面内，将此离散点延伸至三角形的边/点上，否则直接作为一个离散点，而后删除此奇异点；

④将上述得到的所有与边界奇异点相关的离散点进行筛选，去除方向接近的离散点；

⑤计算剩下的内部奇异点的相关领域的第一个离散点；

⑥如果存在奇异点的二领域内有其他奇异点，如有执行⑦，否则跳至⑧；

⑦挑选一个奇异点，将另一奇异点作为此奇异点的一个离散点，如果此离散点在面内，延伸至三角形的边/点，否则直接作为一个离散点，而后删除另一个奇异点；

⑧将上述得到的所有与内部奇异点相关的离散点进行筛选，去除方向接近的离散点；

3.2、根据龙格库塔方法将得到的奇异点与相应的第一个离散点进行流线的延伸，终止条件如下所述：

①直到奇异点的一领域，以奇异点终止；

②以约束边/边界边的终止处的最近的三角网格点终止；

③同一条流线以同样的方向经过同一个背景网格三角形，终止；

步骤4、简化流线结构，构造初始拓扑；

4.1对于步骤3中得到的流线结构，检测其流线情况，进行简化，直到没有一种可简化的情况出现，简化情况包括：

①两条流线起止点相同，分别记为A,B，且长度相差较大，改长的流线为A->A,短流线不变；

②两条流线起止点相同，且长度相差不大，去除其中一条流线；

③两条流线起止点有一个相同记为A，额外的点分别为B,C，不失一般性，设A->B，长于A->C，改长流线为B->C；

4.2求出流线与流线之间的交点，并构造拓扑关系；

步骤5、弦折叠的拓扑简化，基于四边形网格对偶结构以及保持网格中奇异点的数量和基本拓扑关系不变从而设计一种四边形网格拓扑简化操作与优化框架，适用于各类四边形网格的简化；

弦、横边、纵边的定义如下：

弦：在四边形区域的拓扑划分中根据每条边所相对的对边的集合构成的四边区域；

横边：弦内从边界开始到另一条边界对边中所有的对边；

纵边：弦内除横边外的其余边(连接横边顶点的边)；

弦的折叠操作：将弦横边的顶点进行两两合并，对应的纵边上的区域线两两合并；

四边形区域的折叠操作：

四边形区域的折叠就是把弦内部分四边形区域组合而成的多块区域相对的两个顶点进行合并而实现的拓扑区域删除操作；将四边区域想象而成两个三角形，要合并的顶点之间的边是公共边，公共边由弦的两条纵边进行线性插值合成，对应计算公式如下：

q_i＝(1-t_i)×a_i+t_i×b_i,q_i+1＝(1-t_i+1)×a_i+1+t_i+1×b_i+1

式中i表示区域线的序号，a_ia_i+1,b_ib_i+1为弦内纵边上的区域线，a_ib_i为弦内的横边，区域线由许多折线构成，j表示区域线上的折线序号,t_i为

x_i,j为纵边区域线上的点；

步骤5包括如下子步骤：

①：在四边形区域的拓扑划分中根据对边划分弦；

②：将弦按类型分类：

(1)无奇异点存在或者奇异点存在于同一条纵边；

(2)同一条横边上存在两个奇异点；

(3)两条纵边均存在奇异点并且不在同一条横边上；

③：对②中的3类弦进行检查：

(1)判断存在被合并的纵边为边界边或者约束边时，抛弃该弦；

(2)第一类的弦不需进行处理，第二、三类的弦转至④；

④：根据优先级度量的能量公式对选取的弦进行排序，选择能量最大且大于0的的弦进行拓扑区域简化；如果不存在能量大于0的弦，则结束简化；

优先级度量包括：

使用能量公式选择删除那些狭长且不规整的四边区域，将一个四边区域拓扑剖分结构简化为奇异点最少的最优拓扑结构；能量公式如下：

式中K控制弦能量比值；

⑤：检查弦的类型进行区域简化：

(1)第二类型的弦，进行四边形区域的折叠操作；

(2)第三类型的弦将根据奇异点的位置划分为non-zip-patch区域与zip-patch区域；

non-zip-patch区域删除离其最近的奇异点的弦内对边纵边和点，然后与该纵边外的四边区域合并；

zip-patch区域根据原来的横边和纵边线性插值拟合出一条新的纵边，删除原先两条纵边，再将横边根据新点一分为二与弦外边进行合并；

(3)返回①；

步骤6、针对上述得到的区域划分，以背景尺寸长计算每个区域边的分段数，并对此段数进行优化：

步骤6包括如下子步骤：

子步骤1，使用量化计算方式，对所有的区域线结合尺寸场计算理想边数的整数舍入值；所述的区域线上的分段数量化计算方式：

设区域线C，该区域线C由从a到b变化的t参数化；曲线上根据尺寸图的浮点边数为：

其中dl＝||x'(t)||dt；为了得到整数个边n_e，使用一个简单的舍入：

子步骤2，使用四边区域量化的拓扑约束，将一些四边区域的分段数施加相等性；

所述的四边形区域的拓扑约束：

调整四边形区域面相对两侧的边数n_e；通过对先前在每条曲线上计算的理想值求平均值来计算这些点；通过这种简单的传播，弦仅在两个相邻的四边形面共享一条曲线时才传播，该曲线是它们的四个边之一；

子步骤3，得到的最终的理想边数对所有的区域线进行二次量化；

与区域线的第i个内点相关联的参数t_i,i∈[1,n_e-1]为：

通过沿区域线添加值直到总和等于

在样本之间进行线性插值；为了对区域线C进行网格划分，使用公式11计算参数t_i并评估区域线参数化：x_i＝f(t_i)计算顶点位置x_i,i∈[1,n_e-1]；使用这种方法，根据尺寸图将点很好地放置在曲线上，从较小的特征区域平滑过渡到较粗的区域；

步骤7、使用模板法对每个分区进行网格生成；

步骤：

1)按照顺/逆时针输入多边形区域角点坐标和相应边的细分数；

2)根据输入的细分数选择符合的拓扑模板，通过计算相应参数确定patch的拓扑结构；

3)计算边界上点的坐标，使用拉普拉斯光顺确定patch的几何位置；拉普拉斯光顺的公式如下：

subject to v_i＝ω_i,i∈C

其中，ε为网格的边集，C为边界点集，ω_i∈R²为输入第i个固定边界点坐标；

步骤8、网格光顺

①规则顶点：给定一个固定的边界，求解Winslow非线性椭圆PDEΔ_xu＝0，其中u(x)是某个坐标计算空间，x是物理空间中的坐标；通过应用有限差分(FDM)离散化到Winslow方程，推导出局部平滑四边形网格中规则顶点的核；设(x1,...,x8)是有序的围绕规则顶点x的模板的顶点，其新位置由下式给出：

其中α₀＝(x₃-x₇)·(x₃-x₇)，α₁＝(x₁-x₅)·(x₁-x₅)，β＝(x₁-x₅)·(x₃-x₇)

②不规则顶点：使用基于角度的平滑方法，其想法是移动单环平分线上的顶点；作为上面的Winslow内核，它在网格不太受限制的情况下效果很好，但可能会在复杂配置中创建翻转单元。

Images