（计算机应用技术专业论文）任意精度浮点算术在delaunay网格生成算法实现中的应用.pdf

浙江大学硕士学位论文 摘要 摘要 网格生成技术在很多领域都有广泛的应用 如计算机图形学 计算机视觉 可视化 地理信息系统和科学计算 本文主要关注科学计算领域的网格生成技术 按包含内部节点的单元数目是否相同 网格分为结构化和非结构化两类 d e l a u n a y 方法是目前最流行的非结构化网格生成方法之一 d e l a u n a y 网格生成方法牵涉到很多浮点计算 由于计算机浮点数截断误差的 影响 这些计算可能是不精确的 由此可能导致很多意想不到的算法健壮性问题 本文引入任意精度浮点算术 部分解决了此类健壮性问题 本文叙述结构如下 第1 章简单介绍了研究背景和非结构化网格生成技术的研究现状 解释了浮 点运算误差影响d e l a u n a y 网格生成算法健壮性的机理 第2 章详细介绍了软件模拟任意精度浮点算术算法 它适用于符合i e e e 7 5 4 标准的浮点运算部件 它的原理是将浮点算术计算的近似结果与误差部分分别精 确地保存在硬件支持的浮点数中 由多个浮点数共同精确地表示计算结果 算法 的关键是保证误差的精确性 对此本章将给出详细的证明 使用这些算术算法和 自适应技术 我们实现了4 个精确的计算几何谓词 第3 章则将第2 章实现的算术算法和谓词应用到d e l a u n a y 网格生成算法的边 界边恢复和边界面恢复环节 具体地 我们重新设计了线段和三角形单元相交及 求解交点位置的精确算法 以及判断两个共面三角形单元绕向一致性的精确算 法 它们是保证边界边和边界面恢复过程健壮性的关键所在 第4 章则通过实际例子来展示本文算法的效果 首先利用改进后的d e l a u n a y 网格生成算法 成功生成了4 个改进前的算法失效的网格实例 并具体分析了改 进的算法失效的原因 此外 我们展示了网格生成算法应用于具体工程实践中所 获得的网格实例 最后 将算法集成于自主知识产权的c a e 系统h e d p 中 完 成一个夹具的结构振动分析 第5 章总结全文 介绍了d e l a u n a y 方法依然存在的健壮性问题和可能的解决 方案 关键词 d e l a u n a y 三角化 网格生成 健壮性 浮点运算 任意精度浮点算术 运算 边界恢复 浙江大学硕士学位论文 a b s t r a c t m e s hg e n e r a t i o ni sw i d e l ya n ds u c c e s s f u l l yu s e di nm a n ya r e a s s u c ha sc o m p u t e r g r a p h i c s c o m p u t e rv i s i o n v i s u a l i z a t i o n g e o g r a p h i ci n f o r m a t i o ns y s t e m s c i e n t i f i c c o m p u t a t i o na n ds oo n t h et h e s i sf o c u s e so nm e s hg e n e r a t i o nt e c h n o l o g yi ns c i e n t i f i c c o m p u t a t i o na r e a m e s hi sc l a s s i f i e dt os t r u c t u r e do n ea n du n s t r u c t u r e do n ea c c o r d i n g t ow h e t h e rt h en u m b e ro fn o d e si n c i d e n tt oa ni n t e r i o rn o d ei ss i m i l a ro rn o t n l e d e l a u n a ym e t h o di so n eo ft h em o s tp o p u l a ru n s t r u c t u r e dm e s hg e n e r a t i o nm e t h o d s m a n yf l o a t i n g p o i n tc o m p u t a t i o n sa l ei n v o l v e di nt h ed e l a u n a ym e t h o d t h e s e c o m p u t a t i o n sm a yn o tb ee x a c td u et or o u n d o f fe r r o r so ff l o a t i n g p o i n tn u m b e r s c o n s e q u e n t l y m a n yu n e x p e c t e dr o b u s tp r o b l e m sa p p e a r i no r d e rt or e s o l v et h e s e p r o b l e m s as c h e m eb yu s i n ga r b i t r a r yp r e c i s i o nf l o a t i n g p o i n ta r i t h m e t i ci sp r e s e n t e d t h es t r u c t u r eo ft h et h e s i si ss h o w na sf o l l o w s c h a p t e r1i n t r o d u c e st h er e s e a r c hb a c k g r o u n d t h ec u r r e n ts t a t u so fu n s t r u c t u r e d m e s hg e n e r a t i o n a n de x p l a i n sw h yt h er o u n d o f fe r r o r so ff l o a t i n g p o i n tc o m p u t a t i o n s c o u l di n t r o d u c et h er o b u s tp r o b l e m so fd e l a u n a ym e s hg e n e r a t i o n c h a p t e r2i n t r o d u c e st h ea r b i t r a r yp r e c i s i o nf l o a t i n g p o i n ta r i t h m e t i ca l g o r i t h mi n d e t a i l t h ea l g o r i t h mi sf r e eo fr o b u s tp r o b l e m si ft h ef l o a t i n g p o i n tu n i t si np r o c e s s o r s c o m p l y 析t ht h ei e e e 7 5 4f l o a t i n g p o i n ts t a n d a r d t h ei d e ai st oe x a c t l ys t o r ea n a r i t h m e t i cr e s u l t i t ht h ea p p r o x i m a t er e s u l ta n dt h ee x a c tr o u n d o f fe r r o rs e p a r a t e di n f l o a t i n g p o i n tn u m b e r ss u p p o r t e db yh a r d w a r e t h em o s ti m p o r t a n tp a r to ft h e a l g o r i t h mi st og u a r a n t e et h ee x a c t n e s so ft h er o u n d o f fe r r o r s a n dw ew i l lp r o v et h e e x a c t n e s sm a t h e m a t i c a l l y u s i n gt h ea r b i t r a r yp r e c i s i o nf l o a t i n g p o i n ta r i t h m e t i ca n d a d a p t i v et e c h n o l o g y f o u r e x a c tc o m p u t a t i o n a lg e o m e t r yp r e d i c a t e sa r ei m p l e m e n t e d c h a p t e r3a p p l i e st h ea r i t h m e t i ca l g o r i t h ma n dp r e d i c a t e si n t r o d u c e di nc h a p t e r2 i na ni m p l e m e n t a t i o no ft h ed e l a u n a ym e s hg e n e r a t i o nm e t h o d t w oe x a c ta l g o r i t h m s a r ed e s i g n e d o n ei st oj u d g ew h e t h e ras e g m e n ta n dat r i a n g l ei n t e r s e c to rn o t w h e r e t h ei n t e r s e c tp o i n tc o u l db el o c a t e de x a c t l ys i m u l t a n e o u s l y 1 1 1 eo t h e ri st oj u d g e w h e t h e rt h et w oc o p l a n a rt r i a n g l e sh a v es i m i l a ro r i e n t a t i o n s 1 1 1 et w oa l g o r i t h m sa r e k e yt og u a r a n t e er o b u s t n e s so fb o u n d a r ye d g ea n db o u n d a r yf a c er e c o v e r yp r o c e d u r e s i i 浙江大学硕士学位论文 c h a p t e r4p r e s e n t ss o m ee x a m p l e st os h o wt h ee f f e c to ft h ei m p r o v e dd e l a t m a y a l g o r i t h m f o u rm e s he x a m p l e sw h i c hc o u l dn o tb eg e n e r a t e db yt h ef o r m e ra l g o r i t h m a l eg i v e nf i r s t l y a n dt h eu n s u c c e s s f u lr e a s o n sa l ea n a l y z e d s o m eo t h e rm e s h e x a m p l e so b t a i n e df r o mp r a c t i c a ls i m u l a t i o np r o j e c t sa r ep r e s e n t e dt od e m o n s t r a t et h e c a p a b i l i t i e so fo u rm e s h i n ga l g o r i t h m m o r e o v e r b yi n t e g r a t i n g 析t hac a es y s t e m c a l l e dh e d pd e v e l o p e db yo u rl a b c e n t e rf o re n g i n e e r i n ga n ds c i e n t i f i c z h e j i a n g u n i v e r s i t y ac o m p l e t es i m u l a t i o np r o c e s sf o rs t r u c t u r ev i b r a t i o ni sp e r f o r m e d c h a p t e r5c o n c l u d e st h et h e s i s a n dd e p i c t ss o m ed r a w b a c k so ft h ec u r r e n t a l g o r i t h ma n dt h e i rp o s s i b l er e s o l u t i o n s k e y w o r d s d e l a u n a yt r i a n g u l a t i o n m e s hg e n e r a t i o n r o b u s t f l o a t i n g p o i n t c o m p u t a t i o n a r b i t r a r yp r e c i s i o nf l o a t i n g p o i n ta r i t h m e t i c c o m p u t a t i o n a lg e o m e t r y p r e d i c a t e s b o u n d a r yr e c o v e r y 浙江大学硕士学位论文图目录 图目录 图1 1 四叉树算法流程图 3 图1 2 前沿推进法算法流程 4 图1 3 前沿推进法实施过程示意图 5 图1 4 四点共圆时的两种不同的d e l a u n a y 三角化 5 图1 5 插入点时容易出错的情况举例 6 图2 1 浮点数的表示 1 1 图2 2 重叠与不重叠的两个数的定义 一1 l 图2 3e r r o r a b 的最低有效位不会小于 13 图2 4x 与2 y 不重叠的两种情况 1 6 图2 5 算法g r o we x p e n s i o n 数据流向及实现图 1 7 图2 6e x p e n s i o ns u m 算法数据流向及过程图 1 8 图2 7s c a l ee x p e n s i o n 计算过程图 2 1 图2 8 数值6 2 在c c 语言中的单精度表示 2 3 图2 9o r i e n t 2 d 自适应算法实现图 2 7 图2 1 0o r i e n t 3 d 自适应计算实现图 一2 8 图3 1 左图为v o r o n o i 图 右图则是其对应的d e l a u n a y 三角化 2 9 图3 2 二维中的d e l a u n a y 三角形 2 9 图3 3 三维约束d e l a u n a y 网格生成流程图 3 1 图3 4 线段a b 外侧内侧定义 3 2 图3 5 线段p q 与a a b c 不共面且不相交各种情况 3 5 图3 6p q 与a a b c 不共面且相交的各种情况 3 6 图3 7p 在 a b c 平面上时各种相交与不相交情况 一3 7 图3 8a e d 与a a b c 的绕向 3 8 图3 9 构建管道信息的算法流程 一4 0 图3 1 0 管道元与遗失边相交的第1 个交的可能所在的各种位置 4 l 图3 11 构建簇信息的算法流程 4 3 图4 1 实例1 n e ws 4 4 图4 2 实例2 i n t a k ej e ta 4 5 i i i 浙江大学硕士学位论文 图目录 图4 3 实例1 n e ws 网格剖切图 4 5 图4 4 实例2 i n t a k ej e ta 网格剖切图 4 6 图4 5 恢复边界边 1 4 5 3 2 8 时相关两个网格单元情况 4 7 图4 6 边 1 4 5 3 2 8 与三角形 1 3 6 5 3 4 5 3 0 3 的位置 4 8 图4 7 边 1 4 5 3 2 8 与三角形 1 3 1 2 3 4 5 3 1 5 的位置 一4 8 图4 8 边 1 4 5 3 2 8 与三角形 1 3 1 2 3 4 5 3 1 5 的位置 4 8 图4 9 实例3 p a r t l 及其尖角处放大图 4 9 图4 1 0 实例3 p a r t l 网格剖切图 5 0 图4 1 1 实例4 q i u t i 及其尖角处放大图 5 0 图4 1 2 实例4 q i u t i 网格剖切图 5 l 图4 1 3 轴承座实例及其网格剖切图 5 2 图4 1 4 滑件实例及其网格剖切图 5 2 图4 1 5 零件l 实例及其网格剖切图 5 3 图4 1 6 零件2 实例及其网格剖切图 5 3 图4 1 7 各实例网格质量柱状图 5 4 图4 18 夹具周期运动激励信号图 5 6 图4 1 9 导入夹具模型并添加边界条件 5 6 图4 2 0 曲面网格和体网格生成 5 7 图4 2 l 夹具模型网格生成后效果 5 7 图4 2 2 打开求解器并创建边界条件 5 8 图4 2 3 指定边界条件具体意义及模型的材料属性 5 8 图4 2 4 可视化计算结果 5 9 i v 浙江大学硕士学位论文表目录 表目录 表2 1i e e e 7 5 4 中的2 种浮点数据类型 1 1 表2 2i e e e 7 5 4 标准格式与c c 语言中格式对应关系 一2 3 表3 1 各种不相交情况对应o r i e n t 3 d 值 3 3 表3 2 各种相交情况对应o r i e n t 3 d 值 3 3 表3 3p 与a b c 共面时的各种相交与不相交情况的对应o r i e n t 3 d 值 3 4 表4 1 相关边及网格单元顶点坐标 4 7 v 浙江大学研究生学位论文独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作及取得的研究成果 除了文中特别加以标注和致谢的地方外 论文中不包含其他人已经发表或撰写过的研究成 果 也不包含为获得逝 江盘鲎或其他教育机构的学位或证书而使用过的材料 与我一 同工作的同志对本研究所做的任何贡献均已在论文中作了明确的说明并表示谢意 学位论文作者签名 误石叭 签字日期 妒莎年 月7 日 学位论文版权使用授权书 本学位论文作者完全了解逝 江盘堂有权保留并向国家有关部门或机构送交本 论文的复印件和磁盘 允许论文被查阅和借阅 本人授权迸姿盘堂可以将学位论文的 全部或部分内容编入有关数据库进行检索和传播 可以采用影印 缩印或扫描等复制手段 保存 汇编学位论文 保密的学位论文在解密后适用本授权书 学位论文作者签名 劣别i 签字日期 砂8 年6 月 7 e t 翩虢翔挺 签字日期 力谚牮f 月q 日 浙江大学硕士学位论文 第1 章绪论 第1 章绪论 1 1 本文研究内容及背景 很多物理现象可建模为一个偏微分方程系统 科学计算以及工程分析中常用 的偏微分方程的数值分析方法包括有限元法 f i n i t ee l e m e n tm e t h o d 有限体积 法 f i n i t ev o l u m em e t h o d 边界元法 b o u n d a r ye l e m e n tm e t h o d 以及有限差分 法 f i n i t ed i f f e r e n c em e t h o d 这些方法的共通之处在于求解前需将连续几何描 述离散成有限个简单形状的几何形体 这一过程即是网格生成 常用的简单几何 形体包括三角形 四边形 四面体 六面体或棱锥体等 当然 网格生成不仅仅 应用于工程与科学计算中 在计算机图形学 计算机视觉 可视化 地理信息系 统和物理学等领域都有成功的应用i l 书j 按包含内部节点的单元数目是否相同 网格分为结构化和非结构化两类 与 此对应的网格生成方法分别称为结构化网格生成和非结构化网格生成 d e l a u n a y 方法是目前最流行的非结构化网格生成方法之一 评价一个网格生成程序的好坏主要看程序的基本功能 程序的健壮性 生成 的网格质量 运行的速度和效率以及自适应程度 其中程序的健壮性衡量的是网 格生成程序对各种复杂输入条件的适应能力 它是程序能否实用化 进而商品化 的关键性衡量指标 由于计算机浮点数存储位数有限 浮点计算产生超过硬件支持的精度时需要 对结果进行舍入 从而导致不精确的浮点计算结果 d e l a u n a y 网格生成方法中涉 及很多浮点运算 其中有些运算的结果是判断执行不同算法逻辑的依据 以文章 非结构化网格生成及其并行化的若干问题研刭9 j 中介绍的d e l a u n a y 网格生成算法 为例 它包含很多此类运算 如 点插入时判断一点是否在某个网格单元的外接球内 边界边恢复过程中判断某网格单元是否为管道元 以及管道元类型 边界面恢复过程中判断某网格单元是否为簇元 以及簇元类型 可以想象 这类运算结果的不精确性必然会导致程序执行错误的算法逻辑 从而导致很多意想不到的健壮性问题 针对这个问题 本文引入软件模拟任意精 度浮点算术 a r b i t r a r yp r e c i s i o nf l o a t i n g p o i n ta r i t h m e t i c 精确实现此类判断 以 缓解d e l a u n a y 网格生成算法因浮点误差导致的健壮性问题i l 0 1 浙江大学硕士学位论文第1 章绪论 1 2 网格生成算法简介 常用的结构化网格生成算法是映射法 映射方法主要有3 类 超限映射法 代数插值法和基于偏微分方程的方法 映射法几何适应能力通常较差 且很难处 理网格疏密过渡 但在贴体单元的生成和特定物理现象的捕捉上仍然有着自己独 特的优势 因此在如计算流体力学 c f d c o m p u t a t i o n a lf l u i dd y n a m i c s 等特定 研究领域仍然有广泛的应用 目前仍在深入研究之中 由于其良好的边界和约束的适应能力 非结构化网格在很多领域都获得了应 用推广 非结构化网格生成方法主要有栅格法 前沿推进法和d e l a u n a y 方法 1 l 1 2 其中d e l a u n a y 方法因其数学意义明确 执行效率高而被广泛研究和使用 目前二维d e l a u n a y 方法的研究相对成熟 但三维d e l a u n a y 方法的研究远没有二 维问题成熟 其中算法健壮性是关注的焦点之一 l 弘1 7 j 下面将先简单介绍上述非结构化网格生成方法 其中d e l a u n a y 方法还将在第 3 章中重点介绍 1 2 1 栅格法 栅格法基于诸如四叉树和八叉树等空间数据结构分解几何体 由此衍生的典 型方法是四叉树法和八叉树法 分别对应二维和三维区域 s h e p h a r d 和他领导的 研究组首先将四叉树一八叉树法引入到网格生成领域 用于三角形 四面体单元 生成1 1 8 1 9 j 四叉树一八叉树法产生非结构网格的基本做法是 首先用一个立方体 矩形 单元覆盖整个物体区域 然后将其等分为八个 四个 立方体 矩形 然后根据单元是否在内部或边界处决定是否继续划分子立方体 矩形 直到没 有再需要被划分的立方体 矩形 为止 最后剔除掉位于物体外部的立方体 矩 形 并对和物体边界相交的立方体 矩形 进行特殊边界处理 然后将立方体 矩形 切割成四面体 三角形 单元 则可以完成对物体的空间分解 这种方 法尽管效率很高 但在物面附近有可能切割出体积很小的单元 而且它很难保证 边界的形状 最终的网格要加入很多的逻辑关系来确保生成单元的质量和边界的 形状 以四叉树算法为例 图1 1 给出了算法执行的具体流程 2 浙江大学硕士学位论文第l 章绪论 定义个包含整个模型的正方形 取顶点坐标为偶数 方便划分 土 将此正方形作为四叉树的根节点 并标记为p a r t i a l 将其等分为4 个正方形节点 作为此节点的子节点否 1 l 分情况标记这4 个正方形为f u l l p a r t i a l o u t 土 将四叉树标记为p a r t i a l f u l l 的叶 子节点图形分割 形成网格 上 移动边界附近的节点 使其能更加精确的描述模型边界 上 用拉普拉斯光滑化优化内部网格单元 1 l 结束 图1 1 四叉树算法流程图 1 2 2 前沿推进法 前沿推进法有时也称为波前法 前沿推进法对几何边界的适应能力强 最终 浙江大学硕士学位论文第1 章绪论 网格质量较好 是比较受欢迎的一种非结构化的网格生成算法 此算法的健壮性 很依赖于插入节点的算法和生成单元的算法 而算法的执行效率则依赖于数据结 构的设计 一般使用树型结构数据结构 因为树型结构可以提高单元查找算法的 效率 算法生成的网格质量则依赖于边界约束以及插点算法和形成单元等算法的 影响 虽然很多研究者提出了此算法的变形和改进算法 但是这些算法的内核都 是差不多的 其执行流程如下 由边界开始 以一条分界面 分界线 作为前沿 通过插入内部节点和形成网格单元向内部逐渐延伸 直到剖分完毕 这个分界面 分界线 就称为前沿 如图1 3 中的红色区域 此算法的流程如图1 3 对于 二维的情况 算法的主要实施过程如图1 2 所示 图1 2 前沿推进法算法流程 4 浙江人学硕士学位论文第1 章绪论 二 a 图1 3 前沿推进法实施过程示意图 a 初始剖分边界及前沿 b 推进后新插入的节点以及新单元 c 划分完毕之后的网格 1 2 3d e l a u n a y 网格生成法 d e l a u n a y 三角化方法的依据是d i r i c h l e t 在18 5 0 年提出的一种利用己知点集 将平面划分成凸多边形的理论 其基本思想是 给定平面上的点集 a i 则平面 上存在一系列互不重叠的凸多边形 v o r o n o i 多边形 v i 也称v o r o n o i 图 任 意一个v o r o n o i 多边形v k 包含一个给定点a k 并且v k 内的任意一点与a k 的距离 都比与 a i 中其它点的距离小 2 0 2 4 1 这些凸多边形的每一条边都对应两个v o r o n o i 区域及其中心点 a j 中的两个点 我们称之为点偶 将所有的点偶连线就得到 d e l a u n a y 三角化 它是v o r o n o i 图的对偶图 注意 当存在4 点共圆情形时 d e l a u n a y 三角化并不唯一 如图1 4 所示 i i il 图i 4 四点共圆时的两种不同的d e l a u n a y 三角化 二维情形下 任意给定点集的三角化算法中 d e l a u n a y 三角化算法可保证三 浙江大学硕士学位论文 第1 章绪论 角化结果最小内角最大化 此外 d e l a u n a y 算法执行效率比前沿推进法高 这使 得d e l a u n a y 方法在网格生成领域获得广泛的应用 在网格生成领域应用d e l a u n a y 三角化算法需要解决两个基本问题 如何生成区域内部点 如何保证边界约束的完整性 1 3 影响d e l a u n a y 网格生成健壮性关键问题 1 3 1 点插入算法的健壮性 以三维情形为例 在点插入过程中 需要判断新插入的点是否在某个单元所 在的外接球内 这个判断依赖以下一系列浮点运算过程 先求四面体外接球球心 坐标 再通过计算球心与单元顶点的距离求得半径 最后通过比较新插入点与球 心距离的平方和半径的平方确定点是在外接球内 外接球上或外接球外 显然 在四点几乎共面时 由于计算球心坐标时要除以一个很接近与0 的数 那么对应 计算结果误差会较大甚至可能是错误的结果 以致判断出错 在四点几乎共球的 情况下 这个判断同样极易因浮点计算截断误差而出现健壮性问题 例如把本来 没有破坏d e l a u n a y 准则的单元添加到需要被替换的空腔中 或者遗漏破坏 d e l a u n a y 准则的单元 未将其添加到需要替换的空腔中 这些错误都会导致错误 的空腔定义 从而导致点插入过程失败 图1 5 显示了几种在点插入的时候容易 出现错误的情况 a 点 l 乎共线 b 新插入点p 在许多单元外接圆附近 图1 5 插入点时容易出错的情况举例 现有算法通过回退加扰动的机制处理点插入过程中的健壮性问题 实际效果 良好 考虑到任意精度浮点算术的计算复杂度给程序效率带来的巨大影响 因此 在点插入的过程中 我们将使用原来的机制 并未引入任意精度浮点运算来保证 6 浙江大学硕士学位论文第1 苹绪论 此过程的健壮性 2 锄7 1 3 2 边界边和边界面恢复算法的健壮性 如何恢复给定边界是d e l a u n a y 网格生成方法研究的难题之一 特别是三维边 界恢复问题在理论和实践上都存在很多困删9 1 理论上存在不添加辅助点就不能 三角化的多面体 因此健壮的边界恢复算法均需要考虑如何添加辅助点的问题 以文献非结构化网格生成及其并行化的若干问题研究1 9 j 中所设计的边界恢复 算法为例 其中既要判断线面相交情况 同时也需要计算新添加辅助点坐标 即 线面的交点坐标 这两个运算都会因为浮点运算误差问题影响边界恢复算法的 健壮性i 圳 本文研究将表明 第一个问题可通过引入任意精度浮点算术解决 但注意这 是有前提的 即输入数据是没有误差的 遗憾的是 由于利用硬件支持的浮点数 保存辅助点坐标时已经引入截断误差 而且这种误差会因为多次求交运算而累 积 当误差累积到一定程度时 输入点坐标本身的不精确性会引起后续浮点运算 得到的判断结果是不正确的 既使这些运算本身是精确的 利用任意精度方法保 存浮点数可以解决上述问题 但这在实现的时空效率及实现本身的复杂性等方面 都存在很多亟待解决的难题 因此本文只是引入了精确浮点运算解决线面相交判 断等算法的健壮性 由于辅助点坐标误差引起的健壮性问题留待后续工作解决 2 8 3 0 o 1 3 3 对边界点插入及边界恢复算法健壮性问题的解决办法 本节将先简单的介绍有关解决办法的思想 具体的算法会在第2 章和第3 章 中详细介绍和说明 对于1 3 1 和1 3 2 小节中阐述的健壮性问题 究其原因 都是计算机硬件中 浮点运算的截断误差导致的 所以我们将引入任意精度浮点算术算法到这些影响 网格生成程序健壮性的关键环节 来减少或者避免这些运算中由于截断误差导致 的程序健壮性问题 有了任意精度浮点算术算法后 我们就可以使用这些浮点算术算法实现一些 我们需要的一些几何谓词 如判断二维中三点构成三角形的绕向的谓词o r i e n t 2 d 判断一点是否在另外三点所构成的单元的外接圆内的谓词i n c i r c l e 判断三维中四 点构成四面体的绕向的谓词o r i e n t 3 d 以及判断一点是否在另外四点构成的单元 的外接球内的谓词i n s p h e r e 等 本文研究表明 利用这些精确谓词代替基于普通 浮点运算的谓词 可有效提升边界恢复算法的健壮性 7 浙江大学硕士学位论文第1 章绪论 1 4 本文内容组织 本文引入软件模拟任意精度浮点算术算法来解决d e l a u n a y 网格生成算法因浮 点计算截断误差造成的健壮性问题 介绍并实现了相关的软件模拟的任意精度浮 点算术算法 并证明其正确性 使用这些算术算法实现了一些精确的计算几何谓 词 并使用算这些术算法和谓词改善了d e l a u n a y 网格生成算法健壮性 本文内容组织如下 第1 章简单介绍了研究背景和非结构化网格生成技术的研究现状 解释了浮 点运算误差影响d e l a u n a y 网格生成算法健壮性的机理 第2 章详细介绍了软件模拟任意精度浮点算术算法 并使用这些算术算法和 自适应技术实现了4 个精确的计算几何谓词 第3 章则将第2 章介绍的算术算法和谓词应用到d e l a u n a y 网格生成算法的边 界边恢复和边界面恢复环节 第4 章则通过实际例子来展示本文算法的效果 第5 章总结全文 介绍了d e l a u n a y 方法依然存在的健壮性问题和可能的解决 方案 8 浙江大学硕士学位论文第2 章任意精度浮点算术及常用计算几何谓词的实现 第2 章任意精度浮点算术及常用计算几何谓词的实 现 2 1 前言 在网格生成的过程中 由于计算机硬件的浮点截断误差 常常会导致很多问 题 因此 在这些极易受到浮点截断误差影响的过程中 我们需要对浮点运算做 扩展精度的改进 使得判断更加精确 在这方面 j o n a t h a nr i c h a r ds h e w c h u k l 3 l j 做了很多出色的工作 本文的软件模拟的任意精度浮点算术算法是以他的工作作 为参考 本章会介绍两种技术来实现软件模拟任意精度浮点算术算法及谓词 技术一 是一系列关于任意精度浮点算术算法及其正确性证明 这些算法可以使得计算的 结果扩展到几百甚至几千位 i e e e 7 5 4 标准中定义的双精度浮点数也只有5 2 位 技术二则是自适应技术 由于我们最终目的是要利用任意精度的浮点运算实现精 确的计算几何判断谓词o r i e n t 3 d 等 对于这些计算几何谓词 他们并不要求你的 计算结果返回的数值精确 而只需要这个谓词返回的数的符号 s i g n 正确 因 此 如何在保证返回数符号正确的情况下 自适应的提前结束计算将非常有意义 虽然这技术可以使得任意精度浮点运算在大部分时候都比较快 但是在相对误差 较大的情况下 它还是会很慢 软件模拟精确计算算法可以分为许多类 一些精确算法是基于操作整数或者 定点数 而另外一些则是基于操作浮点数 为了表示一个实数 前者是用自定义 的数据结构 这种数据结构可以扩展到任意长度的方法来精确的保存一个实数 而后者则是以一个浮点数作为实数的一个构成部分 c o m p o n e n t 用多个浮点数 来共同表示一个需要用扩展精度保存的实数来达到此目的 随着处理器体系结构 的不断改进和优化 现在很多的处理器对于浮点运算的速度甚至超过的了对于整 数的运算速度 这个趋势同样也反应在任意进度浮点算术算法的发展趋势上 对于大部分的软件模拟任意精度浮点运算是用自己定义的可扩展存储长度的 数据结构来精确保存实数 这样做有一个比较明显的缺点 比如现在做一个加法 2 4 0 0 2 删 如果用这种方法保存计算结果需要8 0 0 位 而如果使用多个数作为一 个数的构成方法 只需要分别将做完加法的最接近的结果2 4 0 0 保存于一个双精度 浮点数 再将计算的误差2 4 0 0 保存在另外一个双精度浮点数中即可 只需1 2 8 位 9 浙江大学硕士学位论文第2 章任意精度浮点算术及常用计算几伺谓词的实现 可以节省很多存储空间 但是用自定义的数据结构也有其自己的优点 因为它表 示实数的方法统一 可以使得各个计算结果都保存同样的一个数据类型中 继续 对这两个数的操作也会相对简单 因为他们都基于同一数据结构 而对于由多个 数组成来表示实数的方法则会可能出现一个结果可能只有2 个浮点数构成 而另 外一个则可能是有6 个甚至更多的数组成 这样使得再对这两个计算结果进行操 作时增加了算法上的难度 在本章中 我们将介绍的是利用浮点运算部件的软件模拟的任意精度浮点算 术算法 并且是用多个浮点数作为一个数的一部分来构成和表示实数 这些算法 将在后面的章节中详细的介绍 不同的硬件有不同的特性 我们这里介绍的这些算法是基于浮点算术运算部 件是操作二进制数 并且必须是以 最接近舍入 的方式进行舍入的 所有符合 i e e e 7 5 4 浮点运算标准的处理器都符合以上两个基本条件 在符合上述两个条件 的情况下 我们将证明算法的正确性1 2 4 最后将通过任意精度浮点算术算法实现几个常用的计算几何谓词 如i n c i r c l e i n s p h e r e o r i e n t 2 d o r i e m 3 d 等 其中判断四个给定点的三维方向的谓词 o f i e m 3 d 在改进d e l a u n a y 网格生成程序中用到的一个非常重要的谓词 同时由于任意精度 的浮点运算的时间运算代价大 我们将利用上面提到的技术二 自适应的方法对 其进行改进 使得其在保证计算结果的正确性上 能尽快的运行 以减少时间上 的运算代价 3 2 3 9 1 2 2 软件模拟任意精度浮点算术运算 2 2 1 背景知识 在介绍软件模拟任意精度浮点算术运算之前 我们将介绍与其相关的许多概 念 以及定义 由于本文只实现在符合i e e e 7 5 4 标准浮点算术运算部件这种硬 件条件下的任意精度浮点算术算法 而并不实现通用的 可在任何硬件条件下都 可以正确执行的算法 所以我们的算法有自身特别的一些条件 这些条件可以简 化算法的设计 2 2 1 1i e e e 7 5 4 标准中的浮点数表示 i e e e 7 5 4 标准中是用图2 1 的方式来表示一个浮点数的 1 位作为这个浮点 数符号位 然后用n 位保存这个浮点数的尾数 而指数则用m 位的带符号整数 保存 一般以补码的方式储存指数 方便对指数进行算术运算 如1 0 1 0 1 0 1 0 1 0 1 这个二进制数在计算机浮点数中是 1 0 1 0 1 0 1 0 1 0 1 2 1 0 来保存这个数的具体值 l o 浙江大学硕士学位论文第2 章任意精度浮点算术及常用计算几何谓词的实现 的 其中尾数为 0 1 0 1 0 1 0 1 0 1 注 由于下面的定义都是基于二进制的数 所 以我们提及的所有的数以及概念都是在二进制的前提下 为叙述方面 在本节的 后续部分都不在特别说明 n n 一1nn 一1 0 图2 1 浮点数的表示 对于i e e e 7 5 4 标准中的浮点数定义了2 种常用的浮点数据类型 单精度 双 精度 详见表2 1 表2 1i e e e 7 5 4 中的2 种浮点数据类型 浮点格式总位数n 对应图2 1 m 对应图2 1 单精度 3 22 3 8 双精度 6 45 21 1 2 2 1 2 软件模拟任意精度浮点算术的相关概念 在这里 我们先介绍一些概念 如果两个数a 1 0 0 1 0 0 1 1 b 1 1 0 0 1 1 1 1 他们中一个数中的最低位1 如a 中1 0 0 1 1 最后这个1 在另外一个数的最高位l 如b 中1 1 0 0 1 最前这个1 之后 即是这个1 代表的数值更大 那我们称这 两个数是重叠的 反之 若一个数的最低位1 在另外一个数的最高位l 之前 那 么这两个数是不重叠的 详见图2 2 a i10 0 10 0 11l a i10 0 0 0 0 0 0i bi1 10 0 1 11 1i b1 1 10 0 1 1 1 1l 数a b 是重叠的数a b 数不重叠的 图2 2 重叠与不重叠的两个数的定义 我们称两个数a 和b 是相邻的 如果它们本身并不重叠 并且2 a 和b 重叠 或者2 b 和a 重叠 比如数a 1 0 0 1 0 0 和b 1 0 它们本身并不重叠 但是a 1 0 0 1 0 0 和2 b 1 0 0 却是重叠的 那么我们称a 和b 是相邻的 同样 数a 1 0 0 0 0 和数b 1 0 0 本身也不重叠 但是不管是a 1 0 0 0 0 和2 b 1 0 0 0 还是2 a 1 0 0 0 0 0 和b 1 0 0 都是不重叠的 那么a 和b 就是不相邻的 对于在2 1 节中提到的两种不同的算法 它们之间最主要的区别就是 用自 定义数据结构软件模拟任意精度浮点算术算法 总是使用不断变化或扩展数据结 构的存储位数 并且用最少的位数来精确储存操作数和计算结果 而使用多个数 浙江大学硕士学位论文 第2 章任意精度浮点算术及常用计算几何谓词的实现 作为构成部分的任意精度浮点算法允许计算结果的误差存在 而误差同样也会被 很快的计算出来 而将计算近似结果以及其误差作为构成一个数的不同部分来储 存 并且这两个数是不重叠的 现在我们介绍关于 最接近舍入 定义 在介绍 最接近舍入 定义之前 我们先介绍另外一个定义 就是一个浮点数最后一位所能表示的数值 u n i ti nl a s t p l a c e 即u l p 举例来说明这个定义 8 个二进制位表示的浮点数a 1 0 1 0 0 1 0 1 b 1 0 1 1 1 0 1 0 那么浮点数a 和b 的u l p 分别是 u l p a 2 u l p a 2 一 那么 一个浮点算术运算部件是 最接近舍入 定义为此运算结果r e s u l t r o n d o f r 与理论上 精确结果r e s u l t 之间的误差i e r r o r l1 2 u l p r e s u l t o 帅d o 曲 举例来说 两个4 位的 二进制浮点数a 0 1 1 1 和b 0 1 0 1 a b 1 0 0 0 1 1 由于计算结果也是保存 在一个4 位的二进制浮点数中 精确截断误差的浮点运算部件保存的结果 r e s u l t m 帅d 0 萨1 0 0 1x 2 5 而e r r o r l 这样定义 最接近舍入 还存在一个歧义 就是当误差i e r r o r l 1 2 u l p r e s u l t r o 帅d o 仃 时 到底怎么取舍 如1 0 1 0 1 这个5 位的 结果要保存到一个4 位的二进制浮点数中时 应该是以1 0 1 1 2 4 还是1 0 1 0 2 4 对应的截断误差e r r o r 分别是一1 和l 这两种表示都符合 最接近舍入 误差的 定义 对于i e e e 7 5 4 标准的浮点算术运算部件 对于此类情况都是按照使得截 断后用于表示浮点数的部分偶数的方式来进行截断的 如1 0 1 0 1 这个例子 截断 后的两种情况前四位分别是1 0 1 0 和1 0 1l 前者为偶数 故取前者的表示方式 即1 o l o 2 4 对应的截断误差e r r o r l 并且本文相关算法也是基于i e e e 7 5 4 浮 点算术运算部件这种模式的 下面 为了叙述方便 我们介绍文章后续部分将使用统一的表示方法 理论 上精确的浮点算术运算我们用a d d s u b 分别表示 而计算机中有截断误差的 浮点算术运算我们分别用 掌表示 而对应a b a b a b 的截断误差分别为 e r r o r a b e r r o r a b e r r o r a