OpenFOAM里面的网格对象（fvMesh）

上次讲到了Time对象，并介绍了我们的程序如何跟system/controlDict文件关联。这次我们介绍程序如何读取网格、网格文件在OpenFOAM里面如何表示、以及如何获取网格信息。

程序如何读取网格

上次介绍setTime.H时提到过，程序读取网格是通过在main函数开头用包含createMesh.H文件来实现的。 现在我们改动我们的HelloWorld程序，包含这个头文件，如下

1
2
3
4
5
6
int main(int argc, char* argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    ...

createMesh.H文件位于$FOAM_SRC/OpenFOAM/include/。 我们可以看看其内容是

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
// createMesh.H
// ~~~~~~~~~~~~
    Foam::Info
        << "Create mesh for time = "
        << runTime.timeName() << Foam::nl << Foam::endl;
    Foam::fvMesh mesh
    (
        Foam::IOobject
        (
            Foam::fvMesh::defaultRegion,
            runTime.timeName(),
            runTime,
            Foam::IOobject::MUST_READ
        )
    );

可以发现其实是创建了一个fvMesh对象，叫做mesh。 创建这个fvMesh对象的时候程序读取了constant/polyMesh/目录下面的
points、faces、owner、neighbour和boundary这5个文件。这5个文件就是OpenFOAM存储网格的文件。它们可以用blockMesh生成，也可以用前处理工具从其他格式的网格文件转换得到（比如用fluentMeshtoFoam、gmshToFoam等等）。 有了这个对象，我们就可以从这个对象获得网格的所有信息（如cell center的坐标、cell的体积等等），这些信息在创建fvMesh对象时，其构造函数就以自动计算好了。

下面举个例子，看看如何获得每个cell 的体积。 简单起见，首先我们准备好一个简单的网格，直接从官方教程的cavity case里面复制blockMeshDict到我们之间创建的HelloWorld case过来：

1
2
3
cd $FOAM_RUN/HelloWorld    #进入之前我们创建的helloWorld case目录
mkdir -p constant/polyMesh #先创建这个目录
cp $FOAM_TUTORIALS/incompressible/icoFoam/cavity/constant/polyMesh/blockMeshDict  constant/polyMesh/blockMeshDict #复制blockMeshDict

修改constant/polyMesh/blockMeshDict文件，把cell个数改为为（2 2 1），也就是X、Y方向各两个cell，Z方向一个cell， 并且把convertToMeters改成1。 如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
...
convertToMeters 1; //改动
vertices
(
    (0 0 0)
    (1 0 0)
    (1 1 0)
    (0 1 0)
    (0 0 0.1)
    (1 0 0.1)
    (1 1 0.1)
    (0 1 0.1)
);
blocks
(
    hex (0 1 2 3 4 5 6 7) (2 2 1) simpleGrading (1 1 1) //改动
);
...

使用blockMesh生成网格。现在constant/polyMesh目录下面就有那5个表示网格的文件了。

接下来还要准备system/fvSchemes文件和system/fvSolutions文件, 即使它们里面是空的，因为fvMesh的构造函数必须读取这两个文件。这里我们用的这两个文件内容(除去头部和尾部)分别是：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
$ cat system/fvSchemes
...
gradSchemes
{
    default         none;
}
divSchemes
{
    default         none;
}
laplacianSchemes
{
    default         none;
}
...

和

1
2
3
solvers
{
}

这已经是程序能初始化mesh所需要的最少的信息了，少了任何一条都不行。

现在我们在HelloWorld程序里面增加如下几行代码，来输出每个cell的体积：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#include "fvCFD.H"
#include "MyLib.H"
int main(int argc, char* argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    ...
    // mesh info
    const scalarField& V = mesh.V(); // cell volume
    forAll(V, i) // forAll是OpenFOAM提供的对容器类循环的宏
    {
        Info << "Volume of cell # " << i << ":  " << V[i] << endl;
    }
    ...
    return 0;
}

这里的const scalarField& V = mesh.V()就是获取mesh对象的所有cell体积，给一个scalarFiled类型的引用。scalrField其实是通过typedef Field<scalar>得到的, 这里初略的理解为一个标量元素的数组，并且可以用[]操作符取元素就够了。（Field<>只是OpenFOAM众多层次的容器类的一层）。

重新wmake编译，运行输出如下

1
2
3
4
5
...
Volume of cell # 0:  0.025
Volume of cell # 1:  0.025
Volume of cell # 2:  0.025
Volume of cell # 3:  0.025

可以发现每个cell的体积都是(1/2)*(1/2)*0.1 = 0.025。

类似的，还可以获取cell center坐标、face center坐标、face的法向量和面积等几何信息，如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
...
   const vectorField& C = mesh.C(); // cell center coordinate
   forAll(C, i)
   {
       Info << "Center of cell # " << i << ":  " << C[i] << endl;
   }
   const vectorField& Cf = mesh.Cf(); // face center coordinate
   forAll(Cf, i)
   {
       Info << "Center of face # " << i << ":  " << Cf[i] << endl;
   }
   const vectorField& Sf = mesh.Sf(); // face normal
   forAll(Sf, i)
   {
       Info << "Normal of internal face # " << i << ":  " << Sf[i] << endl;
   }
   ...

下面这个表格列出了获取网格的一些基本几何信息的函数：

我们暂时先在这里停住，在最后一节再演示如何获取网格的其他信息（数据信息和拓扑信息）。 接下来我们看看任意非结构网格在OpenFOAM内部是如何表示的。

网格在OpenFOAM里面如何表示

我们开始以刚才的2*2*1的网格为例来说明。 constant/polyMesh目录下面的5个文件内容分布是：

  • points : 记录所有的网格顶点，每个点都是由三个分量组成的（三维空间坐标）。 本例中一共有3*3*2 = 18个点。
  • faces : 记录所有的face和组成每个face的顶点编号（就是出现在points里面所有顶点中的次序）。对于本例，一共有20个face，包含了内部face （internal face）和位于边界上的face（boundary face）。
  • boundary: 记录每个boundary patch的信息，包括patch的名字，类型、所包含的face的个数n，以及这些face在所有face中的开始位置（也就是此patch中第一个face出现在face记录的face列表里面的次序）。 每个boundary patch中的face编号在faces中是连续的。 faces 中开头记录的face是internal face。下图展示的是之前例子中的boundary和face：


这里有4个internal face和3个boundary patch，分别是， movingWall（2个face）， fixedWall（6个face）， frontAandBack（8个face）。

在介绍剩下的两个文件owner和neighbour之前，有必要介绍OpenFOAM网格里面的face的owner和neighbour的概念。 无论网格多么复杂，有一个基本特点或者说是要求是： 除了boundary face外的每个face （也就是internal face） 被且仅被两个cell共享。这两个cell分别称为这个face的owner cell和neighbour cell。 规定这个这个face的法向量从owner指向neighbour。 boundary face也有owner cell，但是没有neighbour cell，因为boundary face的法向量指向了计算区域外部，那里没有cell了。现在再看owner和neighbour文件就简单了：

  • owner文件记录了一个整型数组，第i个元素表示第i个face(包括boundary face)的owner cell的在所有cell中的编号(cell ID)，数组大小等于所有face的个数（本例是20）。
  • neighbour文件记录了一个整型数组，第i个元素表示第i个internal face的neighbour cell的在所有cell中的编号(cell ID)，数组大小等于所有internal face的个数（本例是4）。

获取网格信息

之前我们介绍了如何通过V()、C()、Cf()等函数获取一些基本的几何信息，现在我们看如何怎样获取一个face的owner cellID和neighbour cellID。直接看代码：

1
2
3
4
5
6
7
8
9
10
11
// find owner cell and neighbour cell of internal faces
const labelUList& owner = mesh.owner(); // labelUList是一个存贮label的容器。
const labelUList& neighbour = mesh.neighbour();
// loop over ONLY internal face to find owner and neighbour
forAll(owner, facei) //注意这里的owner只存了internal face的owner cellID
{
    const label own = owner[facei];
    const label nei = neighbour[facei];
    Info << "face # " << facei << "'s own: " << own
        << ", nei: " << nei << endl;
}

这段代码输出的是：

1
2
3
4
face # 0's own: 0, nei: 1
face # 1's own: 0, nei: 2
face # 2's own: 1, nei: 3
face # 3's own: 2, nei: 3

现在我们可以验证一下到底是不是这样的， 通过刚刚介绍的C()、 Cf()、Sf()函数输出的位置等信息，我们可以找到每个cell和face的编号以及每个face的法向量如下：


对比刚刚输出的结果可以发现确实是这样。

回到刚才插入的代码，可以发现 fvMesh.owner()和 fvMesh.neighbour()函数可以获取与每个internal face的owner cellID和neighbour cellID。注意这里的owner()返回的数组不包含boundary face的cell ID，这也可以从上例中forAll(owner, facei)只循环4次（4个internal face）可以看出。

那要怎么知道每个boundary face的owner cellID （也就跟这个bounary face相邻的那个cell的编号）呢，毕竟这个信息在处理边界时必须知道。 下面这段代码就是答案

1
2
3
4
5
6
7
8
9
10
11
// loop over boundary faces, except empty boundaries
forAll(mesh.boundary(), patchi) //这里对每个bounary patch循环，mesh.boundary()返回boundary patch的数组，不包含empty类型的boundary patch。
{
    const labelUList& pOwner = mesh.boundary()[patchi].faceCells(); // faceCell函数
    forAll(pOwner, facei)
    {
        Info << "Boundary patch #" << patchi
            << "'s face " << "#" << facei << "'s owner cell : "
            <<  pOwner[facei] << endl;
    }
}

这里是打印出了每个boundary patch的每个face的onwer cellID，所以有两层forAll循环。 boundary patch的所有face的owner cellID数组是通过faceCells()函数获取的。
这段代码的输出是：

1
2
3
4
5
6
7
8
Boundary patch #0's face #0's owner cell : 2
Boundary patch #0's face #1's owner cell : 3
Boundary patch #1's face #0's owner cell : 0
Boundary patch #1's face #1's owner cell : 2
Boundary patch #1's face #2's owner cell : 1
Boundary patch #1's face #3's owner cell : 3
Boundary patch #1's face #4's owner cell : 0
Boundary patch #1's face #5's owner cell : 1

这里的0号bounary patch是movingWall，1号boundary patch是fixedWall。可以看到这里并没有处理因为计算1D/2D问题而声明的empty类型的boundary patch （frontAndBack），事实上，我在写CFD程序时也确实不需要处理empty patch，因为从一个empty patch的某个face上流入的必然等于从其相对那个face上流出的，其合贡献为0。

除了faceCell()函数外，还可以获取boundary patch的更多信息，比如name，start（该patch第一个face中在所有的face中的位置）。更多信息可以查看
$FOAM_SRC/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H文件。

至于为什么是这个文件，我们可以浏览源代码顺藤摸瓜搞清楚这里面的boundary patch和mesh的关系以及其中涉及到的类和相应的文件。 很多时候浏览源代码也可以获得很多信息。

首先这里的mesh.boundary()返回的对象类类型是fvBoundaryMesh, 这可以从fvMesh的头文件
$FOAM_SRC/finiteVolume/fvMesh/fvMesh.H中boundary()的函数声明看到：

1
2
//- Return reference to boundary mesh
const fvBoundaryMesh& boundary() const;

在去查找fvBoundaryMesh类的头文件$FOAM_SRC/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.H发现它是继承自fvPatchList

1
2
3
4
5
6
7
...
class fvBoundaryMesh
:
    public fvPatchList
{
    // Private data
...

而fvPatchList的声明是通过typedef实现的, 在$FOAM_SRC/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatchList.H可以看到

1
typedef PtrList<fvPatch> fvPatchList;

可以看到fvPatchList就是一个fvPatch的指针数组（PtrList）。现在去看看fvPatch类的头文件，在$FOAM_SRC/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H。
这里的fvPatch指的就是我们说的boundary patch。 其Access类型的函数可以返回诸多我们感兴趣的数据，比如

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
...
         //- Return name
         const word& name() const
         {
             return polyPatch_.name();
         }
         ...
         //- Return size
         virtual label size() const
         {
             return polyPatch_.size();
         }
         ...
     // Access functions for geometrical data
         //- Return face centres
         const vectorField& Cf() const;
         //- Return neighbour cell centres
         tmp<vectorField> Cn() const;
         //- Return face area vectors
         const vectorField& Sf() const;
         ...

有时候在写后处理程序时希望通过boundary patch的名字找到boundary patch的编号，这可以通过fvBondaryMesh提供的findPatchID函数实现，

1
2
//- Find patch index given a name
label findPatchID(const word& patchName) const;

举个例子，我们希望找到HelloWorld case里面的fixedWall对应的boundary patch编号，可以用下面的代码实现：

1
2
// find patch ID by name
Info << "Patch fixedWalls's patchID = " << mesh.boundary().findPatchID("fixedWalls") << endl;

程序输出是

1
Patch fixedWalls's patchID = 1