在用有限元法或者有限體積法求解流體力學(xué)問題時(shí)，需要先將求解區(qū)域劃分成網(wǎng)格。區(qū)別于在物體表面生成的網(wǎng)格（surface mesh），我們稱這種劃分三維區(qū)域的網(wǎng)格為體網(wǎng)格（volume mesh）。
體網(wǎng)格根據(jù)其單元形狀可以分為四面體網(wǎng)格（tetra-mesh），六面體網(wǎng)格（hexa-mesh），以及四面體或六面體為主的多面體網(wǎng)格（tetra/hexa-dominated mesh）。而根據(jù)其生成方式又可以分為結(jié)構(gòu)化與非結(jié)構(gòu)化網(wǎng)格（stuctured, non-structured），貼體與非貼體網(wǎng)格（conformal, non-conformal），等等。
對(duì)于四面體網(wǎng)格，較為常用的生成算法包括Delaunay法和波前法（advancing front）等。而對(duì)于六面體網(wǎng)格，較為常用的算法有：映射法（mapping），掃掠法（sweeping），以及正交切割單元法（Cartesian cut-cell）。映射法和掃掠法只適用于特定類型的幾何模型；而正交切割單元法具有較強(qiáng)的普適性，只需要提供模型的表面網(wǎng)格就可以自動(dòng)生成六面體為主的多邊形網(wǎng)格（并非純六面體網(wǎng)格）。這里我們主要介紹一下這種方法。
有些體網(wǎng)格生成算法直接從描述幾何模型的參數(shù)方程出發(fā)（例如映射法，掃掠法），而另一類體網(wǎng)格生成算法從模型的表面網(wǎng)格出發(fā)（例如Delaunay法，波前法），正交切割單元法屬于后者。它對(duì)于輸入的面網(wǎng)格有一些要求：
1) 純?nèi)牵╰riangular）：所有單元均為三角形。
2) 滿足水密性條件（water-tight）：任意一條邊恰好屬于兩個(gè)三角形單元。
3) 滿足流形條件（manifold）：任意兩個(gè)三角形要么不相交，要么共享一條邊。
算法最終得到的是一組網(wǎng)格單元（cut-cell element），每個(gè)網(wǎng)格單元由一個(gè)正方體和若干個(gè)切割面（cut-face）共同構(gòu)成。在必要時(shí)可以計(jì)算出正方體被切割后形成的多面體，所有這些多面體共同構(gòu)成模型內(nèi)部空間區(qū)域的一個(gè)剖分。
1) 初始化：
首先根據(jù)面網(wǎng)格的包圍盒（bounding box），以及用戶指定的分辨率（resolution）生
成均勻的正方體單元，作為第一層網(wǎng)格單元的基礎(chǔ)。
2) 切割：
對(duì)于輸入的面網(wǎng)格中的每個(gè)三角形，確定所有與其相交的正方體單元，并計(jì)算出被正方體切割后得到的切割面（cut-face），儲(chǔ)存在相應(yīng)的正方體單元中。
每個(gè)三角形上的所有切割面構(gòu)成這個(gè)三角形的一個(gè)剖分。
3) 細(xì)分：
對(duì)于在模型內(nèi)部而未與表面網(wǎng)格相交的正方體單元，可直接放入體網(wǎng)格當(dāng)中。
對(duì)于與表面網(wǎng)格相交的正方體單元，根據(jù)其存儲(chǔ)的切割面數(shù)量和法向，決定是否將其細(xì)分（refine）。
對(duì)于滿足細(xì)分條件的單元，將其分成大小相等的8個(gè)小正方體單元，并將原本單元中存儲(chǔ)的所有切割面再次切分后，分別存儲(chǔ)到對(duì)應(yīng)的小單元之中。
4) 組裝：
當(dāng)所有正方體單元均不滿足細(xì)分條件，或者達(dá)到預(yù)設(shè)的最大細(xì)分次數(shù)時(shí)，細(xì)分過程終止。
將所有含有切割面的單元放入體網(wǎng)格當(dāng)中。
此時(shí)，體網(wǎng)格中包含了所有完全位于模型內(nèi)部的單元，以及與模型表面相交的單元，它們構(gòu)成模型內(nèi)部空間區(qū)域的一個(gè)剖分。
3. 數(shù)據(jù)結(jié)構(gòu)
為了完整存儲(chǔ)一套體網(wǎng)格的全部數(shù)據(jù)，我們需要記錄三個(gè)列表：節(jié)點(diǎn)坐標(biāo)的列表（point list），切割面的列表（face list），以及網(wǎng)格的單元列表（cell list）。這里我們把節(jié)點(diǎn)和切割面單獨(dú)存成列表，而不是放到網(wǎng)格單元里面，是為了避免重復(fù)記錄相同的數(shù)據(jù)。例如一個(gè)切割面屬于兩個(gè)網(wǎng)格單元，那么我們只需要記錄一份這個(gè)面的數(shù)據(jù)，然后在兩個(gè)網(wǎng)格單元中分別存儲(chǔ)這個(gè)面在列表中的下標(biāo)即可（相當(dāng)于記錄了它的地址）。
切割面
（cut-face）的存儲(chǔ) 切割面是由若干個(gè)節(jié)點(diǎn)構(gòu)成的平面多邊形，它是從輸入
面網(wǎng)格中的某個(gè)三角形被正 方體切割得到的。
對(duì)于每個(gè)切割面，我們只需
要記錄構(gòu)成它
的所有節(jié)點(diǎn)，以及它所在平 面的法向量即可。
切割面（cut-face）的存儲(chǔ)切
割面是由若干個(gè)節(jié)點(diǎn)構(gòu)成的平面多邊形，它是從輸入面網(wǎng)格中的某個(gè)三角形被正 方體切割得到的。
對(duì)于每個(gè)切割面，我們只需要記錄構(gòu)成它的所有節(jié)點(diǎn)，以及它所在平面的法向量即可。
網(wǎng)格單元（cut-cell）的存儲(chǔ)
每個(gè)網(wǎng)格單元由一個(gè)正方體和若干切割面構(gòu)成，在需要時(shí)可以將其轉(zhuǎn)化為多面體單 元。由于自動(dòng)細(xì)分機(jī)制，一個(gè)單元若被細(xì)分，則存儲(chǔ)其 8 個(gè)子單元的地址（下標(biāo)），因 此全部網(wǎng)格單元構(gòu)成一個(gè)八叉樹結(jié)構(gòu)。
如圖，這個(gè)網(wǎng)格單元中包含 4 個(gè)切割面，每個(gè)切割面由面網(wǎng)格中的一個(gè)三 角形和正方體相交得到。它們將正方體切割成兩個(gè)多面體。只要記錄這些 切割面和立方體，就可以在需要的時(shí)候計(jì)算出切割后形成的兩個(gè)多面體。
這里我們重點(diǎn)討論第一步（切割）中用到的計(jì)算方法。細(xì)分和組裝的部分我們這里先不 詳細(xì)討論（如有需要后面可以為此另寫一篇專題）。切割三角形的計(jì)算可分為三步：
1） 對(duì)于每條三角形邊，計(jì)算其與正方體平面的交點(diǎn)，放入節(jié)點(diǎn)列表中；
同時(shí)，將節(jié)點(diǎn) 有序插入到這條邊的節(jié)點(diǎn)列表中。
這里“有序”指的是每個(gè)節(jié)點(diǎn)的前，后節(jié)點(diǎn)為幾何上與其相 鄰的節(jié)點(diǎn)。
在 C++中可以借助 std::map 來實(shí)現(xiàn)有序插入。
2） 將第一步算出的交點(diǎn)兩兩組成混合邊（hybrid-edge）。
對(duì)于每條混合邊，計(jì)算其 與正方體平面的交點(diǎn)，放入節(jié)點(diǎn)列表中；
同時(shí)，將節(jié)點(diǎn)有序插入到這條邊的節(jié)點(diǎn)列表中（這 部分與第一步類似）。
3） 將第一步和第三步生成的節(jié)點(diǎn)組裝成平面多邊形（即切割面），放入切割面列表中；
同時(shí)，將切割面的下標(biāo)記錄在對(duì)應(yīng)的網(wǎng)格單元中。
下面具體討論這些步驟中用到的特殊方法和結(jié)構(gòu)。
混合邊指的是三角形與切割它的平面的交集所在的線段。由于三角形與平面均是三維空 間中的凸集（convex set），因此它們的交集也是凸集，這意味著無論它們相對(duì)位置如何， 交集都只有一條線段（包括退化為一個(gè)點(diǎn)或者是空集的情況）。也就是說，只要我們記錄下每個(gè)與三角形相交的平面上產(chǎn)生的交點(diǎn)，最后就能得到一條在這個(gè)平面上的混合邊。在切割 三角形的計(jì)算中，混合邊與三角形邊有著同等的重要性（統(tǒng)稱為“邊”）。
有兩種退化情況需要注意：第一種是如果三角形的一條邊恰好落在某個(gè)平面上，那么可 以跳過計(jì)算與平面相交的過程，而直接把這條三角形邊作為平面上的混合邊；另一種情況是， 如果某個(gè)平面與三角形的邊相交得到的兩個(gè)交點(diǎn)過于接近（小于浮點(diǎn)數(shù)誤差），或者恰好經(jīng) 過三角形的某個(gè)端點(diǎn)，那么應(yīng)當(dāng)將這個(gè)端點(diǎn)作為這個(gè)平面上退化形式的混合邊。
另外還需要注意的是，在計(jì)算混合邊與正方體平面的交點(diǎn)時(shí)，若不加處理，會(huì)導(dǎo)致重復(fù) 計(jì)算（見下右圖）。因此對(duì)于這種交點(diǎn)，需要記錄其所在的網(wǎng)格線（grid line），確保每條 網(wǎng)格線上只有一個(gè)這樣的節(jié)點(diǎn)。如果計(jì)算新的交點(diǎn)時(shí)發(fā)現(xiàn)對(duì)應(yīng)的網(wǎng)格線上已經(jīng)存在一個(gè)這樣 的節(jié)點(diǎn)，那么應(yīng)將這個(gè)節(jié)點(diǎn)作為交點(diǎn)，而不是生成新的節(jié)點(diǎn)。
在將節(jié)點(diǎn)組裝成平面多邊形的過程中，我們也需要用到一些計(jì)算方法。
1） 首先，我們需要記錄每個(gè)節(jié)點(diǎn)的相鄰節(jié)點(diǎn)（neigh bor），這可以通過按順序遍歷每 條邊（三角邊或者混合邊）上的節(jié)點(diǎn)來得到。
2） 接下來，對(duì)于每個(gè)節(jié)點(diǎn)，又需要將其所有相鄰節(jié)點(diǎn)按照逆時(shí)針方向排序（因?yàn)樗鼈?都在同一個(gè)平面上，所以是可以做到的）。關(guān)于這部分的具體算法在參考文獻(xiàn)[3]中有詳細(xì) 的描述，這里就不再贅述了。
3） 我們將所有未被組裝的節(jié)點(diǎn)有序?qū)Γ窗脒叄琱alf-edge）放在一個(gè)集合當(dāng)中。從 一條半邊出發(fā)，根據(jù)相鄰節(jié)點(diǎn)的順序找到下一個(gè)節(jié)點(diǎn)，直到回到初始節(jié)點(diǎn)，得到一個(gè)完整的 （順時(shí)針）回路構(gòu)成一個(gè)平面多邊形。將這個(gè)過程中訪問到的半邊從集合中移除，直到集合 為空，此時(shí)我們便得到了三角形上所有的平面多邊形。
這里展示一些用我們自己的代碼生成的正交網(wǎng)格結(jié)果。我們將結(jié)果轉(zhuǎn)化為多面體網(wǎng)格的 形式（即計(jì)算出每個(gè)網(wǎng)格單元中切割后的多面體，并將位于模型內(nèi)部的多面體放入一個(gè)列表 中），并將其輸出為 Paraview 可以打開的文件格式。為了更直觀的看到模型內(nèi)部單元的形 狀和分布，我用 Paraview 中的 clip filter 工具切掉了網(wǎng)格中的一些多面體單元。由于還未做 細(xì)分的部分，這里顯示的僅為用均勻背景網(wǎng)格切割得到的結(jié)果。
參考文獻(xiàn)
[1] Owen, and J. Steven. "An Introduction to Automatic Mesh Generation Algorithms - Part I. " (2016).
[2] Aftosmis, M. J. , M. J. Berger, and J. E. Melton. "Adaptive Cartesian Mesh Generation." (2000).
[3] Tao, M., et al. "Mandoline: robust cut-cell generation for arbitrary triangle meshes." ACM Transactions on Graphics 38.6(2019):1-17.
文章來源：仿真學(xué)習(xí)與應(yīng)用
