当前位置:文档之家› 任意形状的三角形网格划分

任意形状的三角形网格划分

第9卷 第5期1997年9月计算机辅助设计与图形学学报J.CAD&CGVol.9,No.5Sep.,1997

任意曲面的三角形网格划分

陈永府 张 华 陈 兴 李德群

(华中理工大学塑性模拟及模具技术国家重点实验室 武汉 430074)

摘要 把曲面分为可展曲面和不可展曲面,对可展曲面用曲面展开算法展成平面,对不可展曲面用曲面分割算法转化成平面片,在平面上运用Delaunay三角划分法进行网格划分,然后把网格节点反映射到曲面上,从而实现任意曲面的三角形网格划分。

关键词 Delaunay三角划分,可展曲面,不可展曲面,曲面展开算法,曲面分割算法。

1996201203收稿,1996211222收到修改稿。本文得到“八五”重点攻关项目“注塑成形的计算机仿真与交互计算”资助。陈永府,1972年生,硕士研究生,研究方向为注塑模CAE的前处理。张 华,博士研究生,研究方向为注塑模CAE。陈 兴,副教授,研究方向为注塑模CAD。李德群,博士生导师,研究方向为模具CAD、CAE、CAM。1 引 言

网格划分是计算机图形学研究的重要内容。目前,有很多种网格划分算法,如拓扑分

解法、节点连接法、映射单元法、基于栅格法等。这些算法在二维网格划分上各有千秋,但

对三维任意曲面的网格划分,成功的例子却很少。俄国数学家Delaunay在1934年就证明

了:对于任意给定的平面点集,有且仅有一种三解剖分方法能够满足“最大2最小角”优化

准则,即所有三角形的最小内角之和最大。Sibson[1]证明了平面任意给定点集的Delauay三角划分具有整体最优化的性质,这就是说,对于任意给定的平面点集,Delaunay三角划

分能够得到整体最优的三角形网格,能尽可能地避免病态三角形的出现。所以,Delaunay三角划分在许多应用领域,尤其是在实体几何造型和有限元网格自动生成等研究领域,受

到广泛的重视。但是传统的Delaunay三角划分不能对三维任意曲面进行网格划分,为了

解决这个问题,本文把曲面分为可展曲面和不可展曲面,分别对可展曲面采用曲面展开算

法,对不可展曲面采用曲面分割算法;将曲面转化为平面,然后在平面上利用Delaunay三

角划分的优化性质进行网格划分,再将网格节点反映射到曲面上(限于篇幅,本文略去网

格节点反映射到曲面上的算法),从而实现任意曲面的三角形网格划分。

2 平面任意多边形域的Delaunay三角划分

Lawson[2]根据“最大2最小角”优化准则,通过“对角线交换”法则实现了二维给定点

集的Delaunay三角划分。Watson[3]首次提出“插入多边形”,并在此基础上,利用“外接

© 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.__________________________________________________________www.paper.edu.cn圆”准则,实现了无约束的Delaunay三角划分。而实际上,网格划分总是有一定的边界约

束条件,本文在Watson算法的基础上,实现了任意多边形域的Delaunay三角划分。

为了和图形系统统一,规定沿面上任一条边前进时,其左侧定义为图形区域,右侧定

义为非图形区域,即外环按逆时针走向,内环按顺时针走向。不失一般性,假定所有的内环

和外环均由直线段组成,因为任意曲线均可用小直线段进行拟合。任意多边形域的De2

launay三角划分算法如下:令多边形的边数为N,第k条边的起点序号为Lk1,终点序号为Lk2。(1)若N=3,则该多边形本身是一个三角形,划分结束;否则令k=1,转入(2);(2)令k=k+1,若Lk2在有向线段L11L12之左,转入(3),否则转入(2);(3)判断当前多边形其余各边是否与线段L11Lk2或L12Lk2相交,若是,转入(2),否则转入(4);(4)保存节点Lk2至候选节点链表中,若k=N,转入(5),否则转入(2);(5)从候选节点链表中找出与节点L11、L12形成外接圆半径最小的节点L02,则节点L11、L12和L02必然构成一个有效的Delaunay三角形,同时对三角形作如下修正:①若线段L11L02及L12L02均不是当前多边形的边界线段,则令N=N+1,L12=L02,LN1=L02,LN2=L12,转入(1);②若线段L11L02(或L12L02)是当前多边形的第m条边,而线段L12L02(或L11L02)不是当前多边形的边,则令L11=L02(或L12=L02),Lm1=LN1,Lm2=LN2,N=N-1,转入(1);③若线段L11L02和L12L02分别是当前多边形的第m条边和第n条边,则将线段L11L12、第m条边和第n条边从当前多边形中删除,N=N-3,转入(1)。多边形域由Delaunay三角划分后所生成的网格大多数都是细长型的,如图1所示,远不能满足几何造型和有限元分析的要求。为此,还必须对其进一步细化,在边界上和图

形区域内部插入适当数量的新节点,以提高网格的质量。图2是经过细化和Laplace光滑

处理以后的三角形网格。

图1 Delaunay

三角划分图2 细化和Laplace光滑后的网格

3 展开可展曲面

一条平面曲线沿着垂直于该平面的某一直线平行移动而形成的曲面,属于可展曲面,把平面曲线称为该曲面的轮廓线,而把沿其平行移动的直线称为该曲面的主轴。

如图3所示,曲线ABCD为轮廓线,LL1为曲面主轴。这里的轮廓线与曲面边界是有

区别的,轮廓线沿主轴平移形成可展曲面,而可展曲面经过“修剪”后的边界才是这里所指

的“曲面边界”。如图4所示,虚线为轮廓线所形成的可展曲面,实线为曲面边界。有时轮7935期陈永府等:任意曲面的三角形网格划分

© 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.中国科技论文在线__________________________________________________________www.paper.edu.cn廓线可以与部分边界重合。

图3

 可展曲面的轮廓线和主轴图4 曲面边界示意图

为叙述方便,假设主轴与Y轴平行(实际上,可为任意直线),轮廓线和曲面边界均由

直线和圆弧组成。设有N段轮廓线和M条边界,且设SSk1、SSk2、SSk3表示第k段轮廓线的起点、终点和中间任一点(直线时,无第三点,圆弧时有第三点,以下同);PSk1、PSk2、

PSk3是SSk1、SSk2、SSk3转化到平面后所对应的点;SEk1、SEk2、SEk3是第k条边界的点;

PEk1、PEk2、PEk3是SEk1、SEk2、SEk3转化到平面后所对应的点;Length(AB)表示A、B两

点之间轮廓线的长度,若A、B点不在轮廓线上,则表示沿主轴向轮廓线投影后所对应的

长度。则曲面展开算法如下:(1)k=1,PSk1->x=SSk1->x,PSk1->y=SSk1->y,PSk1->z=010。(2)PSk2->x=SSk1->x+Length(SSk1SSk2),PSk2->y=SSk2->y,PSk2->z=010。(3)若kx=PSk2->x,PS(k+1)1->y=PSk2->y,PS(k+1)1->z=010;k=k+1,转入(2);否则,k=1,转入(4)。(4)若k>M,则结束;否则,若第k条边界是直线,则转入(5);若是圆弧,则转入(6)。(5)按照(7)将SEk1、SEk2转化成PEk1、PEk2,且PEk1PEk2为直线;k=k+1,转入(4)。(6)按照(7)将SEk1、SEk2、SEk3转化成PEk1、PEk2、PEk3,且若PEk1->y、PEk2->y与PEk3->y相等,则PEk1PEk2PEk3为直线;否则,PEk1PEk2PEk3为圆弧;k=k+1,转入(4)。(7)若SEki(i=1,2,3)对应于第m段轮廓线,则PEki->x=SSm1->x+Length(SSm1SEki),PEki->y=SEki->y,PEki->z=010。将可展曲面展成平面以后,用任意多边形域的Delaunay三角划分算法对其进行网格

划分,然后把网格节点反映射到曲面上去,从而实现可展曲面的三角形网格划分。

图5 曲面分割示意图4 分割不可展曲面

由于双三次Bezier曲面应用最为广泛,所以本文只讨论双三次Bezier曲面的分割[6]。

411 双三次Bezier曲面的分割设双三次Bezier曲面P(u,w)被等参数线u=us,w=ws分成4个曲面片:Q(t,v)、

R(t,v)、G(t,v)和W(t,v),其中t,v为分割后的曲面参数,如图5所示。893计算机辅助设计与图形学学报1997年

© 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.中国科技论文在线__________________________________________________________www.paper.edu.cn下面求解曲面片Q(t,v)的控制点阵BQ,P(u,w)为

P(u,w)=u3 u2 u 1MBBBMTBw3 w2 w 1T

首先进行参数变换,令t=u󰃗us,v=w󰃗ws,则

P(u,w)=u3 u2 u 1MBBBMTBw3 w2 w 1T

=(tus)3 (tus)2 (tus) 1MBBBMTB(vws)3 (vws)2 (vws) 1T=t3 t2 t 1󰃖diag(us3 us2 us 1)󰃖MBBBMTB󰃖diag(ws3 ws2 ws 1)󰃖v3 v2 v 1T=t3 t2 t 1󰃖MB1000(1-us)us00

(1-us)22(1-us)usus20

(1-us)33(1-us)2us3(1-us)us2us3󰃖BB

󰃖1(1-ws)(1-ws)2(1-ws)3

0ws2(1-ws)ws3(1-ws)2ws00w2s3(1-ws)ws2

000ws3MTBv3 v2 v 1T

由于曲面P(u,w)中的0≤u≤us,0≤w≤ws部分即为曲面片Q(t,v)(0≤t≤1,0≤v≤1),所以,曲面片Q(t,v)的控制点阵BQ为

BQ=1000(1-us)us00

(1-us)22(1-us)usus20

(1-us)33(1-us)2us3(1-us)us2us3

󰃖BB󰃖1(1-ws)(1-ws)2(1-ws)3

0ws2(1-ws)ws3(1-ws)2ws00ws23(1-ws)ws2

000ws3

同理,可求得曲面片R(t,v)、G(t,v)和W(t,v)的控制点阵BR、BG和BW。

412 曲面片的平面度计算设P(u,w)(Α≤u≤Β,Ν≤w≤Ω)是参数曲面,则d(P(u,w),[Α,Β,Ν,Ω])为P(u,w)(Α≤u≤Β,Ν≤w≤Ω)的高度[7],且

d(P(u,w),[Α,Β,Ν,Ω])=supΑ≤u≤ΒΝ≤w≤Ω(inf(data)Χi≥0,i=1,2,3,4Χ1+Χ2+Χ3+Χ4=1)

data=P(u,w)-[Χ1(P(Α,Ν)+∃)+Χ2(P(Β,Ν)-∃)

+Χ3(P(Β,Ω)+∃)+Χ4(P(Α,Ω)-∃)]∃=(P(Α,Ν)+P(Β,Ω)-P(Α,Ω)-P(Β,Ν))󰃗4

上式中,之所以选(P(Α

,Ν)+∃),(P(Β,Ν)-∃),(P(Β,Ω)+∃),(P(Α,Ω)-∃)4个点,是因为这4个点总是共面的[7]。

曲面片的平面度定义为曲面的高度除以曲面片的面积,曲面片的面积以4个角点构9935期陈永府等:任意曲面的三角形网格划分

© 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.中国科技论文在线__________________________________________________________www.paper.edu.cn

相关主题