CN116894282B - 空间点集与多连通网格区域拓扑关系的识别方法及系统 - Google Patents
空间点集与多连通网格区域拓扑关系的识别方法及系统 Download PDFInfo
- Publication number
- CN116894282B CN116894282B CN202311148126.9A CN202311148126A CN116894282B CN 116894282 B CN116894282 B CN 116894282B CN 202311148126 A CN202311148126 A CN 202311148126A CN 116894282 B CN116894282 B CN 116894282B
- Authority
- CN
- China
- Prior art keywords
- triangle
- point
- ray
- eps
- sequence
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 71
- 238000004458 analytical method Methods 0.000 claims abstract description 12
- 101100311458 Lytechinus variegatus SUM-1 gene Proteins 0.000 claims description 42
- 238000012216 screening Methods 0.000 claims description 27
- 230000008569 process Effects 0.000 claims description 26
- 239000002775 capsule Substances 0.000 claims description 25
- 239000007787 solid Substances 0.000 claims description 15
- 238000004422 calculation algorithm Methods 0.000 claims description 14
- 238000012790 confirmation Methods 0.000 claims description 13
- 238000004364 calculation method Methods 0.000 claims description 9
- 238000010276 construction Methods 0.000 claims description 9
- 230000005540 biological transmission Effects 0.000 claims description 7
- 238000007667 floating Methods 0.000 claims description 7
- 239000012634 fragment Substances 0.000 claims description 4
- 238000012546 transfer Methods 0.000 claims description 4
- 230000005484 gravity Effects 0.000 claims description 3
- 239000012530 fluid Substances 0.000 abstract description 4
- 238000010586 diagram Methods 0.000 description 7
- 238000004891 communication Methods 0.000 description 5
- 238000005192 partition Methods 0.000 description 5
- 241001522296 Erithacus rubecula Species 0.000 description 4
- 238000003491 array Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 238000010200 validation analysis Methods 0.000 description 2
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000006870 function Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000011895 specific detection Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Geometry (AREA)
- Computer Hardware Design (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Computational Mathematics (AREA)
- Algebra (AREA)
- Computing Systems (AREA)
- Fluid Mechanics (AREA)
- Mathematical Physics (AREA)
- Image Generation (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明涉及计算流体力学、计算几何学与高性能计算的交叉技术领域，公开了空间点集与多连通网格区域拓扑关系的识别方法及系统，该方法，通过辅助三角形的构造，对点集内元素与多连通网格区域的相对位置进行识别；以及，通过射线与辅助三角形最近交点进行属性分析。本发明解决了现有技术存在的可靠性低、鲁棒性差等问题。
Description
技术领域
本发明涉及计算流体力学、计算几何学与高性能计算的交叉技术领域，具体是空间点集与多连通网格区域拓扑关系的识别方法及系统。
背景技术
当前，在计算流体力学领域，多部件重叠网格系统的并行装配已成为多体相对运动情形下流场非定常模拟的关键手段之一。它基于一定的挖洞准则，实现计算单元、插值单元以及洞内单元的分类。在众多的挖洞准则中，需要优先识别空间点与多连通网格区域的拓扑关系，排除封闭物体壁面包围的洞内点元素，从而在并行化查询运算中降低跨处理器的数据传输总量。图1是一个三维空间部件网格的截面，在外部人工边界面与内部固壁边界面之间，是基于四面体、六面体、三棱柱以及金字塔四种标准单元的无缝填充。图2结合图1的边界条件展示了来自其他部件的网格点与当前多连通网格区域之间三种典型的拓扑关系，即点（P）在网格区域上、点（Q）在固壁边界内与点（R）在人工边界外。对上述关系，研究人员大多采用计算几何中的射线法加以判断，其基本思路是：从当前点向任意方向引出一条射线，如果射线与多面体边界面的交点总数为奇数，那么该点在多面体内部，否则在多面体外部。因为射线法适用于空间点与复连通域拓扑关系的识别，且不受区域凹凸性的约束，所以在工程中得到了广泛的应用。然而，如图3所示，该方法需要同时判断射线与边界面多边形的顶点、棱边和面元素内部相交三种情形，它们对应的数值公差难以统一，加之相邻多边形共享部分顶点与棱边的情形，从而使交点总数的奇偶性判断存在不确定性。这一固有缺陷导致射线法的鲁棒性问题至今未能得到彻底的解决，对于建立具有高可靠性的重叠网格计算模块来说也是非常不利的。
发明内容
为克服现有技术的不足，本发明提供了空间点集与多连通网格区域拓扑关系的识别方法及系统，解决现有技术存在的可靠性低、鲁棒性差等问题。
本发明解决上述问题所采用的技术方案是：
空间点集与多连通网格区域拓扑关系的识别方法，通过辅助三角形的构造，对点集内元素与多连通网格区域的相对位置进行识别；以及，通过射线与辅助三角形最近交点进行属性分析。
作为一种优选的技术方案，包括以下步骤：
S1，三角剖分序列构建：构造固壁边界与人工边界上的三角剖分序列；
S2，三角形序列收集与共享：通过全局规约实现三角形序列的收集与共享；
S3，二叉树检索结构建立：并行建立三角形投影序列的二叉树检索结构；
S4，查询点集初筛与余集构造：使用射线法对查询点集进行初筛且构造余集；
S5，拓扑关系确认：对余集元素与网格区域的拓扑关系进行确认；
S6，结果合并：合并点集初筛与余集确认结果至同一序列中。
作为一种优选的技术方案，S1包括以下步骤：
S11，对网格区域进行分解，将部件网格及其边界面切割成若干片段，若干片段均匀分布在多个处理器上；
S12，每个处理器对人工边界片段与固壁边界片段上面的多边形元素进行遍历，建立三角形剖分序列，并将三角形剖分序列存储在C++标准模板库的向量容器中。
作为一种优选的技术方案，S2包括以下步骤：
S21，利用合信息传递接口MPI的全局规约操作对各处理器上的三角形序列的数据进行收集与共享；
S22，如果记编号为i的处理器PR[i]上的局部三角形序列为TS[i]，元素总量为SIZE[i]，那么经过MPI的规约操作，统计出P[0]至P[M-1]处理器上三角形元素总量为TOTAL=SIZE[0]+...+SIZE[M-1]；
S23，在每个处理器上，以TOTAL为总长度为整体三角形序列Sigma开辟动态内存，再次使用MPI规约操作，使Sigma实现完整的边界面三角形信息存储，实现Sigma=TS[0]+...+TS[M-1]。
作为一种优选的技术方案，S3包括以下步骤：
S31，生成平面上投影三角形序列：设定以空间点集中每个元素为起点的射线统一指向笛卡尔坐标系的x轴正向、y轴正向或z轴正向之一；假设所有的射线皆沿着z轴正向，那么在x-o-y平面上，整体三角形序列Sigma将投影得到新的长度相同的三角形序列Omega；其中，如果Sigma中的某个三角形元素sigma[i]与坐标轴z平行，那么投影后的元素omega[i]将退化成x-o-y平面上的线段，o为坐标系原点；
S32，计算投影三角形包围盒序列的平面范围参数：给定x-o-y平面上的投影三角形omega[i]，计算omega[i]三个顶点之x坐标的最小值xmin[i]、x坐标的最大值xmax[i]，y坐标的最小值ymin[i]以及y坐标的最大值ymax[i]，从而得到当前投影三角形omega[i]的平面包围盒box[i]；box[i]被视为四维Euclid空间中的一个点，其坐标为(xmin[i], xmax[i], ymin[i], ymax[i])，由此便得到与投影三角形序列Omega绑定的四维Euclid空间点列Box={box[i]|i=0,...,TOTAL-1}，其中，TOTAL为投影三角形包围盒序列的元素总量；
S33，计算三角形投影序列的二叉树检索结构之根结点参数：依次遍历四维Euclid空间点列Box的各个元素，计算出各个坐标分量的上界、下界，上界用长度为4的数组upper存储，下界用长度为4的数组lower存储，公式如下：
lower[0]=min{xmin[i]|i=0,...,TOTAL-1},
lower[1]=min{xmax[i]|i=0,...,TOTAL-1},
lower[2]=min{ymin[i]|i=0,...,TOTAL-1},
lower[3]=min{ymax[i]|i=0,...,TOTAL-1},
upper[0]=max{xmin[i]|i=0,...,TOTAL-1},
upper[1]=max{xmax[i]|i=0,...,TOTAL-1},
upper[2]=max{ymin[i]|i=0,...,TOTAL-1},
upper[3]=max{ymax[i]|i=0,...,TOTAL-1}；
定义一个四维Euclid空间中的长方体区域
Domain=(lower[0]-eps,upper[0]+eps)(lower[1]-eps,upper[1]+eps)/>(lower[2]-eps,upper[2]+eps)/>(lower[3]-eps,upper[3]+eps)；其中，/>为集合直积符号，eps为计算机双精度浮点数的误差分辨率，eps的引入可使Domain完全覆盖有限空间点列Box的全体，长方体区域Domain为二叉树检索结构的根结点，Box的首元素box[0]与根结点绑定，且对根结点初始区域进行对半划分使得分割平面与lower[0]和upper[0]均平行；
S34，依次插入孩子结点并生成完整的四维空间二叉树检索结构：依次遍历四维Euclid空间点列Box的剩余TOTAL-1个元素{box[1],...,box[TOTAL-1]}，根据剩余TOTAL-1个元素的坐标确定在当前二叉树上的结点位置，并在当前结点上完成子区域的对半划分；根据交替方向原则，二叉树结点层数每增加一层，空间划分方向改变一次，第m层的空间分割平面与lower[m%4]和upper[m%4]平行；其中，m%s表示整数m除以整数s的余数。
作为一种优选的技术方案，S4中，初筛过程中并不直接判断射线与给定三角形的相对关系，而是在给定三角形的基础上构造与之共面的两个辅助三角形，以此筛查从三角形内部穿过的射线，避免射线与给定三角形顶点或棱边的相交逻辑计算。
作为一种优选的技术方案，S4包括以下步骤：
S41，初始化射线与边界面的两个交点计数器：设置射线与固壁边界交点计数器solid_counter初值为0，射线与人工边界交点计数器artificial_counter初值为0；
S42，结合二叉树检索结构筛查给定射线的三角形包围盒集合：对于由点空间点P出发的射线PQ，设其在x-o-y平面上的投影坐标为(xp,yp)，遍历步骤S3所建立的二叉树，快速检索到与PQ所在直线相交的包围盒子集，与PQ所在直线相交的包围盒子集记为SubBoxSet(xp,yp)；
S43，计算包围盒相关的空间三角形与射线相交的几何参数：在步骤S3基础上，任意一个空间三角形皆与一个平面包围盒绑定；对于给定的射线PQ，设与SubBoxSet(xp,yp)相关的空间三角形子集为SubTriangleSet(xp,yp)；设三角形ABC为SubTriangleSet(xp,yp)的一个元素，三个顶点ABC的空间坐标分别为(xa,ya,za)、(xb,yb,zb)与(xc,yc,zc)，那么三角形ABC的重心G的坐标(xg,yg,zg)表述为xg=(xa+xb+xc)/3,yg=(ya+yb+yc)/3,zg=(za+zb+zc)/3；
取一个无量纲小量tau,根据如下公式计算辅助三角形DEF三个顶点D、E、F的空间坐标值(xd,yd,zd)、(xe,ye,ze)与(xf,yf,zf)：
xd=xg+(1-tau)*(xa-xg), yd=yg+(1-tau)*(ya-yg), zd=zg+(1-tau)*(za-zg)
xe=xg+(1-tau)*(xb-xg), ye=yg+(1-tau)*(yb-yg), ze=zg+(1-tau)*(zb-zg)
xf=xg+(1-tau)*(xc-xg), yf=yg+(1-tau)*(yc-yg), zf=zg+(1-tau)*(zc-zg)
同时，根据如下公式计算辅助三角形RST三个顶点R、S、T的空间坐标值(xr,yr,zr)、(xs,ys,zs)与(xt,yt,zt)：
xr=xg+(1+tau)*(xa-xg), yr=yg+(1+tau)*(ya-yg), zr=zg+(1+tau)*(za-zg)
xs=xg+(1+tau)*(xb-xg), ys=yg+(1+tau)*(yb-yg), zs=zg+(1+tau)*(zb-zg)
xt=xg+(1+tau)*(xc-xg), yt=yg+(1+tau)*(yc-yg), zt=zg+(1+tau)*(zc-zg)
设射线PQ的单位法向量为，沿着射线方向，通过参数c3将它与三角形ABC所在平面交点Q的坐标(xq,yq,zq)表达为：
xq=xp+nx*c3, yq=yp+ny*c3, zq=zp+nz*c3
而在三角形DEF所在平面上，Q的坐标(xq,yq,zq)通过参数c0,c1,c2表达为：
xq=xd*c0+xe*c1+xf*c2, yq=yd*c0+ye*c1+yf*c2, zq=zd*c0+ze*c1+zf*c2
其中c1+c2+c3=1，消去xq,yq,zq，得到线性代数方程组如下：
①
使用的射线PQ穿过辅助三角形DEF内点的逻辑为
C1=(方程组①有唯一解)&&(c0>0)&&(c1>0)&&(c2>0)&&(c3>0)
同理，射线PQ与辅助三角形RST相关的线性方程组为
②
使用的射线PQ穿过辅助三角形RST内点的逻辑为
C2=(方程组②有唯一解)&&(g1>0)&&(g2>0)&&(g3>0)&&(g4>0)
如果存在三角形使射线PQ满足条件(!C1)&&C2，那么在初筛步骤中停止点P与多连通网格区域之间的拓扑关系后将其插入当前处理器的空间点余集RS中；否则执行如下逻辑：如果当前三角形位于人工边界，那么artificial_counter累加1，否则solid_counter累加1；
S44，根据交点奇偶性判断空间点元素与网格区域之间的拓扑关系：结合步骤S43的结果，如果artificial_counter为奇数，那么P点位于人工边界内部，否则位于人工边界外部；如果solid_counter为奇数；那么P点位于固壁边界内部，否则位于固壁边界外部；如果P位于固壁边界外部且位于人工边界内部，那么P位于网格覆盖的区域上。
作为一种优选的技术方案，S5中，对余集元素与多连通网格区域的拓扑关系进行跨处理器识别，对于没有被分布式网格区域上任意体单元覆盖的空间点P，如果距离P点最近的交点位于固壁上，那么P位于固壁边界内部；如果P点最近的交点位于人工边界上或者交点集合为空，那么P位于人工边界外部。
作为一种优选的技术方案，S5包括以下步骤：
S51，在各处理器上独立构造体单元的空间包围盒序列：网格区域上的体单元类型仅限定为四面体、六面体、三棱柱与金字塔四种基本类型；在编号为i的处理器PR[i]上，设子网格块上的体单元序列为Kappa；对Kappa的第j号元素kappa[j]，计算其所有顶点之x坐标的最小值xmin[j]、x坐标的最大值xmax[j]、y坐标的最小值ymin[j]以及y坐标的最大值ymax[j]，z坐标的最小值zmin[j]以及z坐标的最大值zmax[j]，从而得到当前体单元kappa[j]的空间包围盒capsule[j]；capsule[j]被视为六维Euclid空间中的一个点，其坐标为(xmin[j], xmax[j], ymin[j], ymax[j], zmin[j], zmax[j])；得到与体单元序列Kappa绑定的六维Euclid空间点列Capsule ={capsule[j]|j=0,...,SUM-1}；其中，SUM为当前处理器上子网格体单元总量；
S52，计算体单元序列的二叉树检索结构之根结点参数：依次遍历六维Euclid空间点列Capsule的各个元素，计算出各个坐标分量的上界、下界，上界用长度为6的数组upper存储，下界用长度为6的数组lower存储，公式如下：
lower[0]=min{xmin[j]|j=0,...,SUM-1},
lower[1]=min{xmax[j]|j=0,...,SUM-1},
lower[2]=min{ymin[j]|j=0,...,SUM-1},
lower[3]=min{ymax[j]|j=0,...,SUM-1},
lower[4]=min{zmin[j]|j=0,...,SUM-1},
lower[5]=min{zmax[j]|j=0,...,SUM-1},
upper[0]=max{xmin[j]|j=0,...,SUM-1},
upper[1]=max{xmax[j]|j=0,...,SUM-1},
upper[2]=max{ymin[j]|j=0,...,SUM-1},
upper[3]=max{ymax[j]|j=0,...,SUM-1},
upper[4]=max{zmin[j]|j=0,...,SUM-1},
upper[5]=max{zmax[j]|j=0,...,SUM-1}.
定义一个六维Euclid空间中的长方体区域
Range=(lower[0]-eps,upper[0]+eps)(lower[1]-eps,upper[1]+eps)/>(lower[2]-eps,upper[2]+eps)/>(lower[3]-eps,upper[3]+eps)/>(lower[4]-eps,upper[4]+eps)/>(lower[5]-eps,upper[5]+eps)；其中，/>为集合直积符号，eps为计算机双精度浮点数的误差分辨率，eps的引入使Range完全覆盖有限空间点列Capsule的全体；长方体区域Range为二叉树检索结构的根结点，Capsule的首元素capsule[0]与根结点绑定，且对根结点初始空间进行划分，分割平面与lower[0]和upper[0]均平行）；
S53，依次插入孩子结点并生成完整的六维空间二叉树检索结构：依次遍历六维Euclid空间点列Capsule的剩余SUM-1个元素{capsule[1],...,capsule[SUM-1]}，根据它们的坐标确定在当前二叉树上的结点位置，并在当前结点上完成空间划分；根据交替方向原则，二叉树结点层数每增加一层，空间划分方向改变一次；第m层的空间分割平面与lower[m%6]和upper[m%6]平行；其中，m%s表示整数m除以整数s的余数；
S54，跨处理器识别被体单元覆盖的余集元素：首先通过MPI规约操作对每个处理器上的余集元素进行收集，然后对其并集进行共享；在这个跨处理器的数据传输过程中，每个余集元素所在的原始处理器编号以及在余集中的编号被记录下来，以便于查询结果的逆向返回与更新；结合步骤S52与步骤S53中建立的六维空间二叉树结构快速识别覆盖给定余集元素P的空间包围盒子集，进而通过Newton迭代法深度判断P点与相关体单元的覆盖关系；当整体余集元素在六维空间二叉树结构上完成分布式检索后，返回那些被体单元覆盖的元素信息；各处理器上的余集RS[0],RS[1],...,RS[M-1]中的未定元素仅存在两种分类：位于固壁边界内部或位于人工边界外部；
S55，使用射线投影算法对余集中未定元素进行分析：在编号为i的处理器PR[i]上，设P为相应余集RS[i]中的一个未定元素；依次执行步骤S42与步骤S43中方程②的求解过程，使用的射线PQ穿过辅助三角形RST内点的逻辑C2=(方程组②有唯一解)&&(g1>0)&&(g2>0)&&(g3>0)&&(g4>0)；如果射线不与任何三角形相交，那么P位于人工边界外侧；否则，取距离P点距离最近的交点，记为Q，P的位置分为如下两种情况：如果Q点所属的辅助三角形位于人工边界，那么P位于人工边界外部；如果Q点所属的辅助三角形位于固壁边界，那么P位于固壁边界内部。
空间点集与多连通网格区域拓扑关系的识别系统，用于实现所述的空间点集与多连通网格区域拓扑关系的识别方法，包括依次相连的以下模块：
三角剖分序列构建模块：用以，构造固壁边界与人工边界上的三角剖分序列；
三角形序列收集与共享模块：用以，通过全局规约实现三角形序列的收集与共享；
二叉树检索结构建立模块：用以，并行建立三角形投影序列的二叉树检索结构；
查询点集初筛与余集构造模块：用以，使用射线法对查询点集进行初筛且构造余集；
拓扑关系确认模块：用以，对余集元素与网格区域的拓扑关系进行确认；
结果合并模块：用以，合并点集初筛与余集确认结果至同一序列中。
本发明相比于现有技术，具有以下有益效果：
（1）本发明在重叠网格隐式装配的并行实现过程中，点集初筛与余集确认模型能够快速识别封闭固壁边界包围的洞内节点，大幅降低跨处理器查询过程中的数据传输总量，从整体上提升算法的并行执行速度；
（2）本发明辅助三角形的构造与应用，避免了传统射线法与多边形棱边、顶点相交时复杂的公差分析，在分布式内存管理模式下更容易实现；在交点总数奇偶性判断基础上，增加的投影长度比较过程提升了射线分析方法在重叠网格挖洞时的可靠性。
附图说明
图1是本发明三维空间部件网格截面图；
图2是本发明空间点与多连通网格区域的相对位置示意图；
图3是本发明射线与空间多边形相交情形分类示意图；
图4是本发明所述的空间点集与多连通网格区域拓扑关系的识别方法的步骤示意图；
图5是本发明分布式内存环境下点集初筛和余集生成示意图；
图6是本发明三角形抽象类型的数据表达示意图；
图7是本发明辅助三角形的构造及其与射线求交示意图；
图8是本发明分布式内存环境下空间点集之余集确认过程；
图9是本发明实施例Robin直升机混合重叠网格系统示意图；
图10是本发明实施例Robin直升机混合重叠网格截面（挖洞前，y=0）；
图11是本发明实施例Robin直升机混合重叠网格截面（挖洞后，y=0）。
具体实施方式
下面结合实施例及附图，对本发明作进一步的详细说明，但本发明的实施方式不限于此。
实施例1
如图1至图11所示，本发明属于计算流体力学、计算几何学与高性能计算的交叉领域，具体涉及一种在分布式内存管理模式下，空间点集与多连通网格区域拓扑关系的快速识别方法。
本发明公开了一种空间点集与多连通网格区域拓扑关系的识别方法，根据执行顺序，它包含点集初筛与余集确认两个过程。在点集初筛阶段，通过辅助三角形的构造，避免了传统射线法在多边形顶点与棱线附近复杂的公差分析，对点集内绝大多数元素与多连通网格区域的相对位置进行了快速识别。在余集确认阶段，通过射线与辅助三角形最近交点的属性分析，保证了这些空间点与网格区域相对位置识别算法的准确性。
在重叠网格隐式并行装配算法中，本发明可以大幅降低跨处理器查询过程的信息传输总量，从整体上提升算法的并行执行速度。
针对上述问题，本发明结合重叠网格装配算法的具体应用场景，提供一种空间点集与多连通网格区域拓扑关系的识别技术。它依次执行点集初筛与余集确认两个过程。因为从测度分析的角度来看，如图3所示，由P点发出的任意一条射线，它与边界面多边形顶点或棱边相交的概率远小于与多边形内点相交的概率，所以在点集初筛阶段，通过相对误差控制算法，排除从多边形顶点或棱边附近穿过的射线，仅分析严格从多边形内部（远离多边形边界）穿过的射线，根据交点总数的奇偶性判断射线起点（属于初始空间点集）与多连通网格区域拓扑关系。点集初筛算法完成后，初始点集中大多数网格点与当前多连通域的相对位置已经得到精确判断，其它元素将在后续的余集确认过程中加以处理。图3中，P是空间中给定的一个待查询点，PA、PB、PC、PQ分别为以P为起点的四条射线。其中，A位于空间多边形的棱线上，B位于空间多边形的顶点，C为空间多边形内点，Q位于空间多边形外，PQ不与空间多边形相交。仍以图1与图2为例，本发明结合人工边界与固壁边界之间的网格建立二叉树检索结构，经过快速查询判断，如果余集中的点未被网格上的体单元覆盖，那么该点位置只有两种可能，即位于人工边界外侧或者固壁边界内侧。此时，如果从当前点发出的一条射线最先穿越固壁边界，那么它位于固壁边界内侧。反之，如果该射线最先穿越人工边界，或者不与任何边界相交，那么它位于人工边界外侧。在分布式内存管理模式下，每个处理器上针对子网格块独立生成二叉树检索结构，而边界几何信息通过全局规约操作实现整体共享，这样的通信成本在高性能计算平台上是可以快速实现的。因为本发明项目的点集初筛过程已经大幅降低了跨处理器查询数据的总量，所以后续的余集确认过程的时间成本也得到显著缩减。
为实现上述目的，根据图4与图5所示，本发明项目采用的技术解决方案为：一种空间点集与多连通网格区域拓扑关系的识别方法。在分布式内存环境下，它的并行实现过程包含如下6个基本步骤：
步骤S1：构造固壁边界与人工边界上的三角剖分序列。当区域分解完成后，图1标记的部件网格及其边界面被切割成若干片段，均匀分布在多个处理器上。在当前步骤中，每个处理器对人工边界片段与固壁边界片段上面的多边形元素进行遍历，建立三角形剖分序列，并存储在C++标准模板库的向量（vector）容器中。图6是本发明中三角形（MYTriangle）抽象类型的数据表达示意图，它内核包含长度为9的双精度（double）浮点型数组（coordinate），用来存储三个顶点的坐标序列；同时，使用整型（int）变量type记录三角形元素所在边界的属性值（BC_TYPE）；X0、Y0、Z0分别表示顶点0的x轴坐标、y轴坐标、z轴坐标，X1、Y1、Z1分别表示顶点1的y轴坐标、y轴坐标、z轴坐标，X2、Y2、Z2分别表示顶点2的x轴坐标、y轴坐标、z轴坐标。一般地，如果多边形元素顶点N总数为3，那么直接新建一个三角形对象triangle并为其开辟内存空间，记录coordinate与type数据后，将triangle之地址值添加至当前处理器上的三角形向量容器collector中；如果多边形元素顶点总数N大于3，那么通过对角线切割的形式将其分解为N-2个三角形triangle[0], triangle[1],...,triangle[N-3]，依次为每个三角形开辟内存空间并记录coordinate与type数据，然后把每个三角形之地址值添加至当前处理器上的三角形向量容器collector中。
步骤S2：通过全局规约实现三角形序列的收集与共享。因为本发明项目构造的并行算法需要把边界面视为一个整体后进行运算，所以当前步骤结合信息传递接口MPI的全局规约操作对各处理器上的三角形序列的数据进行收集与共享。在图5中，如果记编号为i的处理器PR[i]上的局部三角形序列为TS[i]，元素总量为SIZE[i]，那么经过MPI的规约操作，可以统计出P[0]至P[M-1]处理器上三角形元素总量为TOTAL=SIZE[0]+...+SIZE[M-1]。在每个处理器上，以TOTAL为总长度为整体三角形序列Sigma开辟动态内存，再次使用MPI规约操作，使Sigma实现完整的边界面三角形信息存储，即Sigma=TS[0]+...+TS[M-1]。
步骤S3：并行建立三角形投影序列的二叉树检索结构。因为在后续的步骤中，需要快速检索与给定射线相交的三角形子集，所以这里建立二叉树检索结构以实现这个目的。它依次执行三角形序列投影、平面包围盒序列计算、根结点计算以及四维Euclid空间二叉树的生成四个过程。
进一步地，步骤S3的具体操作步骤可以分解为：
步骤S31：平面上投影三角形序列的生成。为简化计算过程且降低计算成本，本发明项目设定以空间点集中每个元素为起点的射线统一指向笛卡尔坐标系的x轴正向、y轴正向或z轴正向之一。这里不妨假设所有的射线皆沿着z轴正向（选取其它坐标轴方向的情形与之完全类似），那么在x-o-y（o为坐标系原点）平面上，整体三角形序列Sigma将投影得到新的长度相同的三角形序列（记为Omega）。这里需要特别指出的是，如果Sigma中的某个三角形元素sigma[i]与坐标轴z平行，那么投影后的元素omega[i]将退化成x-o-y平面上的线段，但是，本发明项目仍用前述三角形类型对其进行描述，这并不影响后续逻辑运算的正确性。
步骤S32：计算投影三角形包围盒序列的平面范围参数。给定x-o-y平面上的投影三角形omega[i]，计算其三个顶点之x坐标的最小值xmin[i]、x坐标的最大值xmax[i]，y坐标的最小值ymin[i]以及y坐标的最大值ymax[i]，从而得到当前投影三角形omega[i]的平面包围盒box[i]。在本发明项目中，box[i]被视为四维Euclid空间中的一个点，其坐标为(xmin[i], xmax[i], ymin[i], ymax[i])。由此便得到与投影三角形序列Omega绑定的四维Euclid空间点列Box={box[i]|i=0,...,TOTAL-1}，TOTAL为投影三角形包围盒序列的元素总量。
步骤S33：计算三角形投影序列的二叉树检索结构之根结点参数。依次遍历四维Euclid空间点列Box的各个元素，计算出各个坐标分量的上下界（分别用长度为4的数组upper与lower存储），即
lower[0]=min{xmin[i]|i=0,...,TOTAL-1},
lower[1]=min{xmax[i]|i=0,...,TOTAL-1},
lower[2]=min{ymin[i]|i=0,...,TOTAL-1},
lower[3]=min{ymax[i]|i=0,...,TOTAL-1},
upper[0]=max{xmin[i]|i=0,...,TOTAL-1},
upper[1]=max{xmax[i]|i=0,...,TOTAL-1},
upper[2]=max{ymin[i]|i=0,...,TOTAL-1},
upper[3]=max{ymax[i]|i=0,...,TOTAL-1}.
此时，定义一个四维Euclid空间中的长方体区域
Domain=(lower[0]-eps,upper[0]+eps)(lower[1]-eps,upper[1]+eps)/>(lower[2]-eps,upper[2]+eps)/>(lower[3]-eps,upper[3]+eps)，其中/>为集合直积符号，eps为计算机双精度浮点数的误差分辨率，它的引入可使Domain完全覆盖有限空间点列Box的全体。在当前步骤中，长方体区域Domain即为二叉树检索结构的根结点。Box的首元素box[0]与根结点绑定，且对根结点初始区域进行对半划分（分割平面与lower[0]和upper[0]平行）。
步骤S34：依次插入孩子结点并生成完整的四维空间二叉树检索结构。依次遍历四维Euclid空间点列Box的剩余TOTAL-1个元素{box[1],...,box[TOTAL-1]}，根据它们的坐标确定在当前二叉树上的结点位置，并在当前结点上完成子区域的对半划分。根据交替方向原则，二叉树结点层数每增加一层，空间划分方向改变一次。因为从根节点（第0层）开始分割平面与lower[0]和upper[0]平行，所以第m层的空间分割平面与lower[m%4]和upper[m%4]平行（在C++语言中m%s运算表示整数m除以整数s的余数）。
步骤S4：使用射线法对查询点集进行初筛且构造余集。与传统射线法的应用方式不同，初筛过程中并不直接判断射线与给定三角形的相对关系，而是在给定三角形的基础上构造与之共面的两个辅助三角形，以此筛查那些严格从三角形内部穿过的射线，避免射线与给定三角形顶点或棱边的相交逻辑计算。这里承接步骤S31的假设，即所有从待检测点发出的射线皆沿着z轴正向。
进一步地，步骤S4对每一条给定射线的具体检测步骤可以分解为：
步骤S41：初始化射线与边界面的两个交点计数器。设置射线与固壁边界交点计数器solid_counter初值为0，与人工边界交点计数器artificial_counter初值为0。
步骤S42：结合二叉树检索结构筛查给定射线的三角形包围盒集合。对于由点空间点P出发的射线PQ，设其在x-o-y平面上的投影坐标为(xp,yp)，那么，遍历步骤S3所建立的二叉树，即可快速检索到与PQ所在直线相交的包围盒子集，记为SubBoxSet(xp,yp)，SubBoxSet(xp,yp)可以为空集。
步骤S43：计算包围盒相关的空间三角形与射线相交的几何参数。在步骤S3基础上，任意一个空间三角形皆与一个平面包围盒绑定。对于给定的射线PQ，设与SubBoxSet(xp,yp)相关的空间三角形子集为SubTriangleSet(xp,yp)。如图7所示，三角形ABC为SubTriangleSet(xp,yp)的一个元素，三个顶点ABC的空间坐标分别为(xa,ya,za)、(xb,yb,zb)与(xc,yc,zc)，那么它的重心G的坐标(xg,yg,zg)可以表述为
xg=(xa+xb+xc)/3,yg=(ya+yb+yc)/3,zg=(za+zb+zc)/3
取一个无量纲小量tau（本发明项目的后续实施例中取tau=0.0001）,根据如下公式计算辅助三角形DEF（内含于三角形ABC的虚线三角形）三个顶点D、E、F的空间坐标值(xd,yd,zd)、(xe,ye,ze)与(xf,yf,zf)：
xd=xg+(1-tau)*(xa-xg), yd=yg+(1-tau)*(ya-yg), zd=zg+(1-tau)*(za-zg)
xe=xg+(1-tau)*(xb-xg), ye=yg+(1-tau)*(yb-yg), ze=zg+(1-tau)*(zb-zg)
xf=xg+(1-tau)*(xc-xg), yf=yg+(1-tau)*(yc-yg), zf=zg+(1-tau)*(zc-zg)
同时，根据如下公式计算辅助三角形RST（包围三角形ABC的虚线三角形）三个顶点R、S、T的空间坐标值(xr,yr,zr)、(xs,ys,zs)与(xt,yt,zt)：
xr=xg+(1+tau)*(xa-xg), yr=yg+(1+tau)*(ya-yg), zr=zg+(1+tau)*(za-zg)
xs=xg+(1+tau)*(xb-xg), ys=yg+(1+tau)*(yb-yg), zs=zg+(1+tau)*(zb-zg)
xt=xg+(1+tau)*(xc-xg), yt=yg+(1+tau)*(yc-yg), zt=zg+(1+tau)*(zc-zg)
设射线PQ的单位法向量为，沿着射线方向，通过参数c3将它与三角形ABC所在平面交点Q的坐标(xq,yq,zq)表达为：
xq=xp+nx*c3, yq=yp+ny*c3, zq=zp+nz*c3
而在三角形DEF所在平面上，Q的坐标(xq,yq,zq)通过参数c0,c1,c2表达为：
xq=xd*c0+xe*c1+xf*c2, yq=yd*c0+ye*c1+yf*c2, zq=zd*c0+ze*c1+zf*c2
其中c1+c2+c3=1，消去xq,yq,zq，得到线性代数方程组如下：
①
本发明项目使用的射线PQ穿过辅助三角形DEF内点的逻辑为
C1=(方程组①有唯一解)&&(c0>0)&&(c1>0)&&(c2>0)&&(c3>0)
同理，射线PQ与辅助三角形RST相关的线性方程组为
②
本发明项目使用的射线PQ穿过辅助三角形RST内点的逻辑为
C2=(方程组②有唯一解)&&(g1>0)&&(g2>0)&&(g3>0)&&(g4>0)
如果存在三角形使射线PQ满足条件(!C1)&&C2，那么在初筛步骤中停止点P与多连通网格区域之间的拓扑关系后将其插入当前处理器的空间点余集RS中；否则执行如下逻辑：如果当前三角形位于人工边界，那么artificial_counter累加1，否则solid_counter累加1。
步骤S44：根据交点奇偶性判断空间点元素与网格区域之间的拓扑关系。结合步骤S43的结果，如果artificial_counter为奇数，那么P点位于人工边界内部，否则位于人工边界外部。如果solid_counter为奇数，那么P点位于固壁边界内部，否则位于固壁边界外部。进一步地，如果P位于固壁边界外部且位于人工边界内部，那么P位于网格覆盖的区域上。
步骤S5：对余集元素与网格区域的拓扑关系进行确认。经过初始筛查，余集元素总量已经远小于初始空间点集的元素总量。从当前步骤开始，本发明项目需要对余集元素与多连通网格区域的拓扑关系进行跨处理器识别。对于没有被分布式网格区域上任意体单元覆盖的空间点P，如果距离P点最近的交点位于固壁上，那么P位于固壁边界内部；如果P点最近的交点位于人工边界上或者交点集合为空，那么P位于人工边界外部。
进一步地，步骤S5的具体步骤可以分解为：
步骤S51：在各处理器上独立构造体单元的空间包围盒序列。在本发明项目中，网格区域上的体单元类型仅限定为四面体、六面体、三棱柱与金字塔四种基本类型。在编号为i的处理器PR[i]上，设子网格块上的体单元序列为Kappa。对Kappa的第j号元素kappa[j]，计算其所有顶点之x坐标的最小值xmin[j]、x坐标的最大值xmax[j]，y坐标的最小值ymin[j]以及y坐标的最大值ymax[j]，z坐标的最小值zmin[j]以及z坐标的最大值zmax[j]，从而得到当前体单元kappa[j]的空间包围盒capsule[j]。在本发明项目中，capsule[j]被视为六维Euclid空间中的一个点，其坐标为(xmin[j], xmax[j], ymin[j], ymax[j], zmin[j], zmax[j])。由此便得到与体单元序列Kappa绑定的六维Euclid空间点列Capsule ={capsule[j]|j=0,...,SUM-1}。SUM为当前处理器上子网格体单元总量。
步骤S52：计算体单元序列的二叉树检索结构之根结点参数。依次遍历六维Euclid空间点列Capsule的各个元素，计算出各个坐标分量的上下界（分别用长度为6的数组upper与lower存储），即
lower[0]=min{xmin[j]|j=0,...,SUM-1},
lower[1]=min{xmax[j]|j=0,...,SUM-1},
lower[2]=min{ymin[j]|j=0,...,SUM-1},
lower[3]=min{ymax[j]|j=0,...,SUM-1},
lower[4]=min{zmin[j]|j=0,...,SUM-1},
lower[5]=min{zmax[j]|j=0,...,SUM-1},
upper[0]=max{xmin[j]|j=0,...,SUM-1},
upper[1]=max{xmax[j]|j=0,...,SUM-1},
upper[2]=max{ymin[j]|j=0,...,SUM-1},
upper[3]=max{ymax[j]|j=0,...,SUM-1},
upper[4]=max{zmin[j]|j=0,...,SUM-1},
upper[5]=max{zmax[j]|j=0,...,SUM-1}.
此时，定义一个六维Euclid空间中的长方体区域
Range=(lower[0]-eps,upper[0]+eps)(lower[1]-eps,upper[1]+eps)/>(lower[2]-eps,upper[2]+eps)/>(lower[3]-eps,upper[3]+eps)/>(lower[4]-eps,upper[4]+eps)/>(lower[5]-eps,upper[5]+eps)，其中/>为集合直积符号，eps为计算机双精度浮点数的误差分辨率，它的引入可使Range完全覆盖有限空间点列Capsule的全体。在当前步骤中，长方体区域Range即为二叉树检索结构的根结点。Capsule的首元素capsule[0]与根结点绑定，且对根结点初始空间进行划分（分割平面与lower[0]和upper[0]平行）。
步骤S53：依次插入孩子结点并生成完整的六维空间二叉树检索结构。依次遍历六维Euclid空间点列Capsule的剩余SUM-1个元素{capsule[1],...,capsule[SUM-1]}，根据它们的坐标确定在当前二叉树上的结点位置，并在当前结点上完成空间划分。根据交替方向原则，二叉树结点层数每增加一层，空间划分方向改变一次。因为从根节点（第0层）开始分割平面与lower[0]和upper[0]平行，所以第m层的空间分割平面与lower[m%6]和upper[m%6]平行（在C++语言中m%s运算表示整数m除以整数s的余数）。
步骤S54：跨处理器识别被体单元覆盖的余集元素。图8展示了本发明项目在分布式内存环境下对空间点集之余集进行确认的过程。首先通过MPI规约操作对每个处理器上的余集元素进行收集，然后对其并集进行共享。在这个跨处理器的数据传输过程中，每个余集元素所在的原始处理器编号以及在余集中的编号被记录下来，以便于查询结果的逆向返回与更新。结合步骤S52与步骤S53中建立的六维空间二叉树结构快速识别覆盖给定余集元素P的空间包围盒子集，进而通过Newton迭代法深度判断P点与相关体单元的覆盖关系。当整体余集元素在六维空间二叉树结构上完成分布式检索后，返回那些被体单元覆盖的元素信息，至此，各处理器上的余集RS[0],RS[1],...,RS[M-1]中的未定元素仅存在两种分类，即位于固壁边界内部或位于人工边界外部。
步骤S55：使用射线投影算法对余集中未定元素进行分析。这里仍然承接步骤S31的假设，即所有从待检测点发出的射线皆沿着z轴正向。在编号为i的处理器PR[i]上，设P为相应余集RS[i]中的一个未定元素。依次执行步骤S42与步骤S43中方程②的求解过程，这里仅考虑给定三角形的外侧辅助三角形与射线相交的结果，即使用的射线PQ穿过辅助三角形RST内点的逻辑C2=(方程组②有唯一解)&&(g1>0)&&(g2>0)&&(g3>0)&&(g4>0)。如果射线不与任何三角形相交，那么P位于人工边界外侧。否则，取距离P点距离最近的交点，记为Q，分两种情况讨论：如果Q点所属的辅助三角形位于人工边界，那么P位于人工边界外部。如果Q点所属的辅助三角形位于固壁边界，那么P位于固壁边界内部。
步骤S6：合并点集初筛与余集确认结果至同一序列中。根据本发明项目的点集初筛结果和余集确认结果，记录原始空间点集所有元素与多连通网格区域拓扑关系。
本发明在重叠网格隐式装配的并行实现过程中，点集初筛与余集确认模型能够快速识别封闭固壁边界包围的洞内节点，大幅降低跨处理器查询过程中的数据传输总量，从整体上提升算法的并行执行速度；
本发明辅助三角形的构造与应用，避免了传统射线法与多边形棱边、顶点相交时复杂的公差分析，在分布式内存管理模式下更容易实现；在交点总数奇偶性判断基础上，增加的投影长度比较过程提升了射线分析方法在重叠网格挖洞时的可靠性。
实施例2
如图1至图11所示，作为实施例1的进一步优化，在实施例1的基础上，本实施例还包括以下技术特征：
为使本领域的普通技术人员能够更好地理解本发明的技术方案，下面以国家数值风洞工程中的通用CFD软件——风雷为例，对本发明的技术方案做进一步的描述。原有的风雷软件具备动态部件网格的重叠装配功能，从而能够进行直升机悬停与前飞状态下的非定常流场模拟。然而，伴随着网格规模的增加，重叠装配模块在并行通信方面的时间成本与空间成本迅速增加。结合具体实例分析发现，部件壁面边界包围的大量网格点是在参与了大量集合通信后被识别为洞内点的。图9展示了Robin直升机的混合重叠网格系统，其背景网格为块自适应型笛卡尔网格（网格点总量650万，体单元总量623万），包围旋翼与机身贴体网格为各向异性的结构网格（网格点总量1334万，体单元总量1306万）。图10为截面y=0处的网格重叠结构，经过程序统计，各个固壁边界包围的背景网格点总量为237万，约占背景网格总量的36.5%。采用本发明项目的技术后，在点集初筛阶段，位于固壁边界内部的背景网格点有98%得到了快速识别，只有剩余的2%参与了跨处理器的数据传输，并在余集确认中得到了精确识别。同时，重叠装配过程所需的CPU时间降低了20%。图11截面y=0处的最终挖洞结果，它充分验证了本发明项目算法的准确性。实施例表明，本发明的应用保证了几何分析结果的正确性，同时降低了重叠网格隐式并行装配算法的通信成本。
如上所述，可较好地实现本发明。
本说明书中所有实施例公开的所有特征，或隐含公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合和/或扩展、替换。
以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，依据本发明的技术实质，在本发明的精神和原则之内，对以上实施例所作的任何简单的修改、等同替换与改进等，均仍属于本发明技术方案的保护范围之内。
Claims (2)
1.空间点集与多连通网格区域拓扑关系的识别方法，其特征在于，通过辅助三角形的构造，对点集内元素与多连通网格区域的相对位置进行识别；以及，通过射线与辅助三角形最近交点进行属性分析；
包括以下步骤：
S1，三角剖分序列构建：构造固壁边界与人工边界上的三角剖分序列；
S2，三角形序列收集与共享：通过全局规约实现三角形序列的收集与共享；
S3，二叉树检索结构建立：并行建立三角形投影序列的二叉树检索结构；
S4，查询点集初筛与余集构造：使用射线法对查询点集进行初筛且构造余集；
S5，拓扑关系确认：对余集元素与网格区域的拓扑关系进行确认；
S6，结果合并：合并点集初筛与余集确认结果至同一序列中；
S1包括以下步骤：
S11，对网格区域进行分解，将部件网格及其边界面切割成若干片段，若干片段均匀分布在多个处理器上；
S12，每个处理器对人工边界片段与固壁边界片段上面的多边形元素进行遍历，建立三角形剖分序列，并将三角形剖分序列存储在C++标准模板库的向量容器中；
S2包括以下步骤：
S21，利用合信息传递接口MPI的全局规约操作对各处理器上的三角形序列的数据进行收集与共享；
S22，如果记编号为i的处理器PR[i]上的局部三角形序列为TS[i]，元素总量为SIZE[i]，那么经过MPI的规约操作，统计出P[0]至P[M-1]处理器上三角形元素总量为TOTAL=SIZE[0]+...+SIZE[M-1]；
S23，在每个处理器上，以TOTAL为总长度为整体三角形序列Sigma开辟动态内存，再次使用MPI规约操作，使Sigma实现完整的边界面三角形信息存储，实现Sigma=TS[0]+...+TS[M-1]；
S3包括以下步骤：
S31，生成平面上投影三角形序列：设定以空间点集中每个元素为起点的射线统一指向笛卡尔坐标系的x轴正向、y轴正向或z轴正向之一；假设所有的射线皆沿着z轴正向，那么在x-o-y平面上，整体三角形序列Sigma将投影得到新的长度相同的三角形序列Omega；其中，如果Sigma中的某个三角形元素sigma[i]与坐标轴z平行，那么投影后的元素omega[i]将退化成x-o-y平面上的线段，o为坐标系原点；
S32，计算投影三角形包围盒序列的平面范围参数：给定x-o-y平面上的投影三角形omega[i]，计算omega[i]三个顶点之x坐标的最小值xmin[i]、x坐标的最大值xmax[i]，y坐标的最小值ymin[i]以及y坐标的最大值ymax[i]，从而得到当前投影三角形omega[i]的平面包围盒box[i]；box[i]被视为四维Euclid空间中的一个点，其坐标为(xmin[i], xmax[i], ymin[i], ymax[i])，由此便得到与投影三角形序列Omega绑定的四维Euclid空间点列Box={box[i]|i=0,...,TOTAL-1}，其中，TOTAL为投影三角形包围盒序列的元素总量；
S33，计算三角形投影序列的二叉树检索结构之根结点参数：依次遍历四维Euclid空间点列Box的各个元素，计算出各个坐标分量的上界、下界，上界用长度为4的数组upper存储，下界用长度为4的数组lower存储，公式如下：
lower[0]=min{xmin[i]|i=0,...,TOTAL-1},
lower[1]=min{xmax[i]|i=0,...,TOTAL-1},
lower[2]=min{ymin[i]|i=0,...,TOTAL-1},
lower[3]=min{ymax[i]|i=0,...,TOTAL-1},
upper[0]=max{xmin[i]|i=0,...,TOTAL-1},
upper[1]=max{xmax[i]|i=0,...,TOTAL-1},
upper[2]=max{ymin[i]|i=0,...,TOTAL-1},
upper[3]=max{ymax[i]|i=0,...,TOTAL-1}；
定义一个四维Euclid空间中的长方体区域
Domain=(lower[0]-eps,upper[0]+eps)(lower[1]-eps,upper[1]+eps)/>(lower[2]-eps,upper[2]+eps)/>(lower[3]-eps,upper[3]+eps)；其中，/>为集合直积符号，eps为计算机双精度浮点数的误差分辨率，eps的引入可使Domain完全覆盖有限空间点列Box的全体，长方体区域Domain为二叉树检索结构的根结点，Box的首元素box[0]与根结点绑定，且对根结点初始区域进行对半划分使得分割平面与lower[0]和upper[0]均平行；
S34，依次插入孩子结点并生成完整的四维空间二叉树检索结构：依次遍历四维Euclid空间点列Box的剩余TOTAL-1个元素{box[1],...,box[TOTAL-1]}，根据剩余TOTAL-1个元素的坐标确定在当前二叉树上的结点位置，并在当前结点上完成子区域的对半划分；根据交替方向原则，二叉树结点层数每增加一层，空间划分方向改变一次，第m层的空间分割平面与lower[m%4]和upper[m%4]平行；其中，m%s表示整数m除以整数s的余数；
S4中，初筛过程中并不直接判断射线与给定三角形的相对关系，而是在给定三角形的基础上构造与之共面的两个辅助三角形，以此筛查从三角形内部穿过的射线，避免射线与给定三角形顶点或棱边的相交逻辑计算；
S4包括以下步骤：
S41，初始化射线与边界面的两个交点计数器：设置射线与固壁边界交点计数器solid_counter初值为0，射线与人工边界交点计数器artificial_counter初值为0；
S42，结合二叉树检索结构筛查给定射线的三角形包围盒集合：对于由点空间点P出发的射线PQ，设其在x-o-y平面上的投影坐标为(xp,yp)，遍历步骤S3所建立的二叉树，快速检索到与PQ所在直线相交的包围盒子集，与PQ所在直线相交的包围盒子集记为SubBoxSet(xp,yp)；
S43，计算包围盒相关的空间三角形与射线相交的几何参数：在步骤S3基础上，任意一个空间三角形皆与一个平面包围盒绑定；对于给定的射线PQ，设与SubBoxSet(xp,yp)相关的空间三角形子集为SubTriangleSet(xp,yp)；设三角形ABC为SubTriangleSet(xp,yp)的一个元素，三个顶点ABC的空间坐标分别为(xa,ya,za)、(xb,yb,zb)与(xc,yc,zc)，那么三角形ABC的重心G的坐标(xg,yg,zg)表述为xg=(xa+xb+xc)/3,yg=(ya+yb+yc)/3,zg=(za+zb+zc)/3；
取一个无量纲小量tau,根据如下公式计算辅助三角形DEF三个顶点D、E、F的空间坐标值(xd,yd,zd)、(xe,ye,ze)与(xf,yf,zf)：
xd=xg+(1-tau)*(xa-xg), yd=yg+(1-tau)*(ya-yg), zd=zg+(1-tau)*(za-zg)
xe=xg+(1-tau)*(xb-xg), ye=yg+(1-tau)*(yb-yg), ze=zg+(1-tau)*(zb-zg)
xf=xg+(1-tau)*(xc-xg), yf=yg+(1-tau)*(yc-yg), zf=zg+(1-tau)*(zc-zg)
同时，根据如下公式计算辅助三角形RST三个顶点R、S、T的空间坐标值(xr,yr,zr)、(xs,ys,zs)与(xt,yt,zt)：
xr=xg+(1+tau)*(xa-xg), yr=yg+(1+tau)*(ya-yg), zr=zg+(1+tau)*(za-zg)
xs=xg+(1+tau)*(xb-xg), ys=yg+(1+tau)*(yb-yg), zs=zg+(1+tau)*(zb-zg)
xt=xg+(1+tau)*(xc-xg), yt=yg+(1+tau)*(yc-yg), zt=zg+(1+tau)*(zc-zg)
设射线PQ的单位法向量为，沿着射线方向，通过参数c3将它与三角形ABC所在平面交点Q的坐标(xq,yq,zq)表达为：
xq=xp+nx*c3, yq=yp+ny*c3, zq=zp+nz*c3；
而在三角形DEF所在平面上，Q的坐标(xq,yq,zq)通过参数c0,c1,c2表达为：
xq=xd*c0+xe*c1+xf*c2, yq=yd*c0+ye*c1+yf*c2, zq=zd*c0+ze*c1+zf*c2
其中c1+c2+c3=1，消去xq,yq,zq，得到线性代数方程组如下：
①
使用的射线PQ穿过辅助三角形DEF内点的逻辑为
C1=(方程组①有唯一解)&&(c0>0)&&(c1>0)&&(c2>0)&&(c3>0)
同理，射线PQ与辅助三角形RST相关的线性方程组为
②
使用的射线PQ穿过辅助三角形RST内点的逻辑为
C2=(方程组②有唯一解)&&(g1>0)&&(g2>0)&&(g3>0)&&(g4>0)
如果存在三角形使射线PQ满足条件(!C1)&&C2，那么在初筛步骤中停止点P与多连通网格区域之间的拓扑关系后将其插入当前处理器的空间点余集RS中；否则执行如下逻辑：如果当前三角形位于人工边界，那么artificial_counter累加1，否则solid_counter累加1；
S44，根据交点奇偶性判断空间点元素与网格区域之间的拓扑关系：结合步骤S43的结果，如果artificial_counter为奇数，那么P点位于人工边界内部，否则位于人工边界外部；如果solid_counter为奇数；那么P点位于固壁边界内部，否则位于固壁边界外部；如果P位于固壁边界外部且位于人工边界内部，那么P位于网格覆盖的区域上；
S5中，对余集元素与多连通网格区域的拓扑关系进行跨处理器识别，对于没有被分布式网格区域上任意体单元覆盖的空间点P，如果距离P点最近的交点位于固壁上，那么P位于固壁边界内部；如果P点最近的交点位于人工边界上或者交点集合为空，那么P位于人工边界外部；
S5包括以下步骤：
S51，在各处理器上独立构造体单元的空间包围盒序列：网格区域上的体单元类型仅限定为四面体、六面体、三棱柱与金字塔四种基本类型；在编号为i的处理器PR[i]上，设子网格块上的体单元序列为Kappa；对Kappa的第j号元素kappa[j]，计算其所有顶点之x坐标的最小值xmin[j]、x坐标的最大值xmax[j]、y坐标的最小值ymin[j]以及y坐标的最大值ymax[j]，z坐标的最小值zmin[j]以及z坐标的最大值zmax[j]，从而得到当前体单元kappa[j]的空间包围盒capsule[j]；capsule[j]被视为六维Euclid空间中的一个点，其坐标为(xmin[j], xmax[j], ymin[j], ymax[j], zmin[j], zmax[j])；得到与体单元序列Kappa绑定的六维Euclid空间点列Capsule ={capsule[j]|j=0,...,SUM-1}；其中，SUM为当前处理器上子网格体单元总量；
S52，计算体单元序列的二叉树检索结构之根结点参数：依次遍历六维Euclid空间点列Capsule的各个元素，计算出各个坐标分量的上界、下界，上界用长度为6的数组upper存储，下界用长度为6的数组lower存储，公式如下：
lower[0]=min{xmin[j]|j=0,...,SUM-1},
lower[1]=min{xmax[j]|j=0,...,SUM-1},
lower[2]=min{ymin[j]|j=0,...,SUM-1},
lower[3]=min{ymax[j]|j=0,...,SUM-1},
lower[4]=min{zmin[j]|j=0,...,SUM-1},
lower[5]=min{zmax[j]|j=0,...,SUM-1},
upper[0]=max{xmin[j]|j=0,...,SUM-1},
upper[1]=max{xmax[j]|j=0,...,SUM-1},
upper[2]=max{ymin[j]|j=0,...,SUM-1},
upper[3]=max{ymax[j]|j=0,...,SUM-1},
upper[4]=max{zmin[j]|j=0,...,SUM-1},
upper[5]=max{zmax[j]|j=0,...,SUM-1}.
定义一个六维Euclid空间中的长方体区域
Range=(lower[0]-eps,upper[0]+eps)(lower[1]-eps,upper[1]+eps)/>(lower[2]-eps,upper[2]+eps)/>(lower[3]-eps,upper[3]+eps)/>(lower[4]-eps,upper[4]+eps)/>(lower[5]-eps,upper[5]+eps)；其中，/>为集合直积符号，eps为计算机双精度浮点数的误差分辨率，eps的引入使Range完全覆盖有限空间点列Capsule的全体；长方体区域Range为二叉树检索结构的根结点，Capsule的首元素capsule[0]与根结点绑定，且对根结点初始空间进行划分，分割平面与lower[0]和upper[0]均平行）；
S53，依次插入孩子结点并生成完整的六维空间二叉树检索结构：依次遍历六维Euclid空间点列Capsule的剩余SUM-1个元素{capsule[1],...,capsule[SUM-1]}，根据它们的坐标确定在当前二叉树上的结点位置，并在当前结点上完成空间划分；根据交替方向原则，二叉树结点层数每增加一层，空间划分方向改变一次；第m层的空间分割平面与lower[m%6]和upper[m%6]平行；其中，m%s表示整数m除以整数s的余数；
S54，跨处理器识别被体单元覆盖的余集元素：首先通过MPI规约操作对每个处理器上的余集元素进行收集，然后对其并集进行共享；在这个跨处理器的数据传输过程中，每个余集元素所在的原始处理器编号以及在余集中的编号被记录下来，以便于查询结果的逆向返回与更新；结合步骤S52与步骤S53中建立的六维空间二叉树结构快速识别覆盖给定余集元素P的空间包围盒子集，进而通过Newton迭代法深度判断P点与相关体单元的覆盖关系；当整体余集元素在六维空间二叉树结构上完成分布式检索后，返回那些被体单元覆盖的元素信息；各处理器上的余集RS[0],RS[1],...,RS[M-1]中的未定元素仅存在两种分类：位于固壁边界内部或位于人工边界外部；
S55，使用射线投影算法对余集中未定元素进行分析：在编号为i的处理器PR[i]上，设P为相应余集RS[i]中的一个未定元素；依次执行步骤S42与步骤S43中方程②的求解过程，使用的射线PQ穿过辅助三角形RST内点的逻辑C2=(方程组②有唯一解)&&(g1>0)&&(g2>0)&&(g3>0)&&(g4>0)；如果射线不与任何三角形相交，那么P位于人工边界外侧；否则，取距离P点距离最近的交点，记为Q，P的位置分为如下两种情况：如果Q点所属的辅助三角形位于人工边界，那么P位于人工边界外部；如果Q点所属的辅助三角形位于固壁边界，那么P位于固壁边界内部。
2.空间点集与多连通网格区域拓扑关系的识别系统，其特征在于，用于实现权利要求1所述的空间点集与多连通网格区域拓扑关系的识别方法，包括依次相连的以下模块：
三角剖分序列构建模块：用以，构造固壁边界与人工边界上的三角剖分序列；
三角形序列收集与共享模块：用以，通过全局规约实现三角形序列的收集与共享；
二叉树检索结构建立模块：用以，并行建立三角形投影序列的二叉树检索结构；
查询点集初筛与余集构造模块：用以，使用射线法对查询点集进行初筛且构造余集；
拓扑关系确认模块：用以，对余集元素与网格区域的拓扑关系进行确认；
结果合并模块：用以，合并点集初筛与余集确认结果至同一序列中。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311148126.9A CN116894282B (zh) | 2023-09-07 | 2023-09-07 | 空间点集与多连通网格区域拓扑关系的识别方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311148126.9A CN116894282B (zh) | 2023-09-07 | 2023-09-07 | 空间点集与多连通网格区域拓扑关系的识别方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116894282A CN116894282A (zh) | 2023-10-17 |
| CN116894282B true CN116894282B (zh) | 2023-11-24 |
Family
ID=88311077
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311148126.9A Active CN116894282B (zh) | 2023-09-07 | 2023-09-07 | 空间点集与多连通网格区域拓扑关系的识别方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116894282B (zh) |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6373489B1 (en) * | 1999-01-12 | 2002-04-16 | Schlumberger Technology Corporation | Scalable visualization for interactive geometry modeling |
| CN101303414A (zh) * | 2008-05-22 | 2008-11-12 | 北京航空航天大学 | 一种基于水平集的地层面及地质体生成方法 |
| CN102004922A (zh) * | 2010-12-01 | 2011-04-06 | 南京大学 | 基于骨架特征的高分辨率遥感影像飞机提取方法 |
| CN107038322A (zh) * | 2017-05-26 | 2017-08-11 | 哈尔滨工程大学 | 一种任意形状放射源的辐射剂量仿真方法 |
| CN112992294A (zh) * | 2021-04-19 | 2021-06-18 | 中国空气动力研究与发展中心计算空气动力研究所 | 多孔介质lbm计算网格生成方法 |
| CN113255251A (zh) * | 2021-07-14 | 2021-08-13 | 中国空气动力研究与发展中心低速空气动力研究所 | 一种真实感冰型的渲染方法 |
| CN113689556A (zh) * | 2021-10-25 | 2021-11-23 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种块自适应型笛卡尔网格快速图映射方法及系统 |
| CN114494650A (zh) * | 2022-04-06 | 2022-05-13 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种分布式非结构网格跨处理器面对接方法及系统 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US11188738B2 (en) * | 2017-02-08 | 2021-11-30 | The Research Foundation For The State University Of New York | System and method associated with progressive spatial analysis of prodigious 3D data including complex structures |
-
2023
- 2023-09-07 CN CN202311148126.9A patent/CN116894282B/zh active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6373489B1 (en) * | 1999-01-12 | 2002-04-16 | Schlumberger Technology Corporation | Scalable visualization for interactive geometry modeling |
| CN101303414A (zh) * | 2008-05-22 | 2008-11-12 | 北京航空航天大学 | 一种基于水平集的地层面及地质体生成方法 |
| CN102004922A (zh) * | 2010-12-01 | 2011-04-06 | 南京大学 | 基于骨架特征的高分辨率遥感影像飞机提取方法 |
| CN107038322A (zh) * | 2017-05-26 | 2017-08-11 | 哈尔滨工程大学 | 一种任意形状放射源的辐射剂量仿真方法 |
| CN112992294A (zh) * | 2021-04-19 | 2021-06-18 | 中国空气动力研究与发展中心计算空气动力研究所 | 多孔介质lbm计算网格生成方法 |
| CN113255251A (zh) * | 2021-07-14 | 2021-08-13 | 中国空气动力研究与发展中心低速空气动力研究所 | 一种真实感冰型的渲染方法 |
| CN113689556A (zh) * | 2021-10-25 | 2021-11-23 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种块自适应型笛卡尔网格快速图映射方法及系统 |
| CN114494650A (zh) * | 2022-04-06 | 2022-05-13 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种分布式非结构网格跨处理器面对接方法及系统 |
Non-Patent Citations (1)
| Title |
|---|
| 并行化非结构重叠网格隐式装配技术;常兴华 等;《航空学报》;第39卷(第06期);第1-11页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116894282A (zh) | 2023-10-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Lin et al. | Collision and proximity queries | |
| Tang et al. | Interactive continuous collision detection between deformable models using connectivity-based culling | |
| Markowsky et al. | Fleshing out wire frames | |
| Krishnan et al. | Spherical shell: A higher order bounding volume for fast proximity queries | |
| Romney et al. | An efficient system for geometric assembly sequence generation and evaluation | |
| Wang et al. | A survey of 3D solid reconstruction from 2D projection line drawings | |
| Kallmann | Shortest Paths with Arbitrary Clearance from Navigation Meshes. | |
| JPH0756678B2 (ja) | 対話形形状モデリングシステム | |
| Lin et al. | A syntactic approach to 3-D object representation | |
| CN116894282B (zh) | 空间点集与多连通网格区域拓扑关系的识别方法及系统 | |
| Sulaiman et al. | Bounding volume hierarchies for collision detection | |
| CN115564925B (zh) | 基于B-rep模型和笛卡尔网格切片的网格生成方法 | |
| Horvat et al. | Ray-casting point-in-polyhedron test | |
| CN101609565A (zh) | 基于L-Rep模型的三维实体布尔运算方法 | |
| Oaks et al. | Polyhedral Mesh Generation. | |
| Ulyanov et al. | Interactive vizualization of constructive solid geometry scenes on graphic processors | |
| Barki et al. | Contributing vertices-based Minkowski sum of a nonconvex--convex pair of polyhedra | |
| Daum et al. | Boundary representation-based implementation of spatial BIM queries | |
| Byrne et al. | Applications of the VOLA format for 3D data knowledge discovery | |
| Choe et al. | Reduction of LiDAR point cloud maps for localization of resource-constrained robotic systems | |
| CN117422810B (zh) | 结构与参数引导的室内要素规则化与关系推理方法及终端 | |
| Borrmann et al. | An iterative, octree-based algorithm for distance computation between polyhedra with complex surfaces | |
| Veltkamp | Hierarchical approximation and localization | |
| Krishnan et al. | Interactive boundary computation of boolean combinations of sculptured solids | |
| Schinko et al. | Vertex Climax: Converting Geometry into a Non-nanifold Midsurface. |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |