上次讲到了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
mkdir -p constant/polyMesh
cp $FOAM_TUTORIALS /incompressible/icoFoam/cavity/constant/polyMesh/blockMeshDict constant/polyMesh/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;
}
...
和
这已经是程序能初始化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"
...
const scalarField& V = mesh.V();
forAll(V, i)
{
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();
forAll(C, i)
{
Info << "Center of cell # " << i << ": " << C[i] << endl ;
}
const vectorField& Cf = mesh.Cf();
forAll(Cf, i)
{
Info << "Center of face # " << i << ": " << Cf[i] << endl ;
}
const vectorField& Sf = mesh.Sf();
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：
在介绍剩下的两个文件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
const labelUList& owner = mesh.owner();
const labelUList& neighbour = mesh.neighbour();
forAll(owner, facei)
{
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
forAll(mesh.boundary(), patchi)
{
const labelUList& pOwner = mesh.boundary()[patchi].faceCells();
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
const fvBoundaryMesh& boundary () const
在去查找fvBoundaryMesh类的头文件$FOAM_SRC/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.H发现它是继承自fvPatchList
1
2
3
4
5
6
7
...
class fvBoundaryMesh
:
public fvPatchList
{
...
而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。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
...
const word& name () const
{
return polyPatch_.name();
}
...
virtual label size () const
{
return polyPatch_.size();
}
...
const vectorField& Cf () const
tmp<vectorField> Cn() const ;
const vectorField& Sf () const
...
有时候在写后处理程序时希望通过boundary patch的名字找到boundary patch的编号，这可以通过fvBondaryMesh提供的findPatchID函数实现，1
2
label findPatchID (const word& patchName) const ;
举个例子，我们希望找到HelloWorld case里面的fixedWall对应的boundary patch编号，可以用下面的代码实现：1
2
Info << "Patch fixedWalls's patchID = " << mesh.boundary().findPatchID("fixedWalls" ) << endl ;
程序输出是1
Patch fixedWalls's patchID = 1
