参考：

感谢原作者们的无私引路和宝贵工作。

前置：
OpenFOAM开发编程基础03 时间类初步 | 𝓐𝓮𝓻𝓸𝓼𝓪𝓷𝓭 (aerosand.cn)

建立本文的项目文件夹并进入

1
2
3
4
// terminal 
cd /home/aerosand/aerosand/ofsp
mkdir 04_mesh
cd 04_mesh

网格类方法

应用准备

1
2
3
4
5
// terminal
foamNewApp 04_01_mesh
cd 04_01_mesh
cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity debug_case
code .

脚本和说明

使用上一篇讨论相同脚本内容即可，略

mesh 是什么

我们先不看 createMesh.H 的源代码，先从应用的需求切入，也就是，我们需要基于具体的 case 建立的计算用的 mesh ，这个 mesh 当然也有自己的数据（包括 point, face 等等），也要有自己的方法（包括返回点列，返回面心等等）。

primitiveMesh

OpenFOAM 提供纯粹基于几何元素构造的 primitiveMesh 类。

通过 Doxygen 查找可以看到 primitiveMesh 不继承自任何其他类，它就是最基础的基类。

1
2
3
// termianl
find $FOAM_SRC -iname primitiveMesh.H
// /usr/lib/openfoam/openfoam2306/src/OpenFOAM/meshed/primitiveMesh/primitiveMesh.H

找到 primitiveMesh.H 并打开，我们看一下它提供的构造方法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
...
class primitiveMesh
{
	...
protected:
	...
public:
	...
	//- Construct from components
	primitiveMesh
	(
		const label nPoints,
		const label nInternalFaces,
		const label nFaces,
		const label nCells	
	);
	...
...

可以看到 primitiveMesh 类从基本的点、面、单元构造对象。

点、面、单元由 blockMesh 生成，或者由第三方网格软件生成再转换成 OpenFOAM 格式。无论如何，点、面、单元被生成在 case/constant/polyMesh/ 文件夹下。注意，这里的几何数据只是数据而已，不具备任何的方法。

我们需要将这些离散数据封装到信息流中，以方便后续使用。

OpenFOAM 提供 IOobject 类来封装接入这些数据。

1
2
3
// terminal 
find $FOAM_SRC -iname IOobject.H
// /usr/lib/openfoam/openfoam2306/src/OpenFOAM/db/IOobject/IOobject.H

摘取 IOobject.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
...
class IOobject
:
	public IOobjectOption
{
public:
	...
private:
	...
protected:
	...
public:
	...
	inline IOobject
	(
		const word& name, // 基于的对象名称
		const fileName& instance, // 基于的对象文件名称
		const objectRegistry& registry, // 注册相关，暂时无需了解
		IOobjectOption::readOption rOpt, // 读取选项
		IOobjectOption::writeOption wOpt = IOobjectOption::NO_WRITE, // 写入选项
		bool registerObject = true, // 有默认值，可忽略
		bool globalObject = false // 有默认值，可忽略
		// 暂时了解用法即可
	);
	...
...

polyMesh

在 primitiveMesh 类的基础上，OpenFOAM 又提供派生的 polyMesh 类。polyMesh 类基本上还是几何拓扑的，具有一些几何拓扑的方法。

1
2
// terminal
find $FOAM_SRC -iname polyMesh.H

摘取 polyMesh.H 的一部分如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
...
class polyMesh
: 
	public objectRegistry,
	public primitiveMesh
{
public:
	...
private:
	...
public:
	...
	//- Read construct from IOobject
	explicit polyMesh(const IOobject& io, const bool doInt = true);
	// explicit 关键词禁用了隐式转换，避免发生意外的类型转换
	...
...

使用 polyMesh 类，我们可以构造 mesh 对象，并进行一些网格操作

例如 04_01_mesh.C 中写入

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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#include "fvCFD.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"

    Foam::word regionName(Foam::polyMesh::defaultRegion);

    Foam::polyMesh mesh
    (
        IOobject
        (
            // fvMesh::defaultRegion,
            // 可以直接使用defaultRegion，也可以新建一个regionName
            regionName,
            runTime.timeName(),
            runTime,
            IOobject::MUST_READ
        )
    );

    Info<< "Max cell centre: " << max(mesh.cellCentres()) << endl;
    Info<< "Max cell volumes: " << max(mesh.cellVolumes()) << endl;
    Info<< "Max cell face cetres: " << max(mesh.faceCentres()) << endl;
	// 为了减少输出，每个都套了一个 max 函数

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< nl;
    runTime.printExecutionTime(Info);

    Info<< "End\n" << endl;

    return 0;
}


// ************************************************************************* //

有读者可能会疑问 mesh 的构造为什么基于 time 文件夹而不是 polyMesh 子文件夹，毕竟网格的几何信息并不在 time 文件夹下（其实 time 类的相关定义有向后搜寻，是包括 mesh 信息的）。这里涉及到类的实现源代码，不用深究也不建议深究。

编译运行，结果如下

1
2
3
4
5
6
7
8
9
10
Create time

Max cell centre: (0.0975 0.0975 0.005)
Max cell volumes: 2.5e-07
Max cell face cetres: (0.1 0.1 0.01)

ExecutionTime = 0 s  ClockTime = 0 s

End

fvMesh

OpenFOAM 在 polyMesh 类的基础上，添加了有限体积的方法，进而派生出了 fvMesh 类。

1
2
// terminal
find $FOAM_SRC -iname fvMesh.H

看一下 fvMesh.H 中的代码片段

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
...
class fvMesh
:
    public polyMesh,
    public lduMesh,
    public fvSchemes,
    public surfaceInterpolation,    // needs input from fvSchemes
    public fvSolution,
    public data
{
protected:
	...
public:
	...
	explicit fvMesh(const IOobject& io, const bool doInit=true);
	...
...

使用 fvMesh 类，我们可以构造 mesh 对象，并进行一些操作，包括网格的操作。

例如 04_01_mesh.C 中重新写入

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
27
28
29
30
31
32
33
34
35
36
#include "fvCFD.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"

    Foam::word regionName(Foam::polyMesh::defaultRegion);

    Foam::fvMesh mesh
    (
        IOobject
        (
            regionName,
            runTime.timeName(),
            runTime,
            IOobject::MUST_READ
        )
    );

    Info<< "Max cell centre: " << max(mesh.C()) << endl;
    Info<< "Max cell volumes: " << max(mesh.V()) << endl;
    Info<< "Max cell face cetres: " << max(mesh.Cf()) << endl;


    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< nl;
    runTime.printExecutionTime(Info);

    Info<< "End\n" << endl;

    return 0;
}

编译运行，结果如下

1
2
3
4
5
6
7
8
9
10
11
Create time

Max cell centre: max(C) [0 1 0 0 0 0 0] (0.1 0.1 0.005)
Max cell volumes: max(V) [0 3 0 0 0 0 0] 2.5e-07
Max cell face cetres: max(Cf) [0 1 0 0 0 0 0] (0.1 0.1 0.005)

ExecutionTime = 0 s  ClockTime = 0 s

End


读者可以把 max 函数去掉重新编译查看结果。基于这三个子类和父类的讨论，我们可以想见 createMesh.H 主要是什么代码语句。

createMesh.H

看一下和网格相关的头文件 createMesh.H

1
2
3
// termianl
find $FOAM_SRC -iname createMesh.H
// /usr/lib/openfoam/openfoam2306/src/OpenFOAM/include/createMesh.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
27
28
29
30
31
32
Info<< "Default region: " << fvMesh::defaultRegion << endl;
Info<< "Mesh sudirectory: " << fvMesh::meshSubDir << endl;
// OpenFOAM 预设了默认的域和网格子路径

Foam::autoPtr<Foam::fvMesh> meshPtr(nullptr);
// 新建指向网格对象的指针，初始化为空指针
Foam::word regionName(Foam::polyMesh::defaultRegion);
// 新建区域名称，初始化默认域

Foam::Info << "Create mesh";
Foam::Info << " for time = " << runTime.timeName() << Foam::nl;

meshPtr.reset // 通过meshPtr的reset方法重设指针指向
(
	new Foam::fvMesh // 基于参数构造fvMesh类的指针
	(
		Foam::IOobject
		(
			regionName,
			runTime.timeName(),
			runTime,
			Foam::IOobject::MUST_READ
		),
		false
	)
);
meshPtr().init(true);

Foam::fvMesh& mesh = meshPtr();

Foam::Info << Foam::endl;

如果不使用指针，也可以像前文一样使用 fvMesh 直接创建 mesh 。

主源码

主源码正式内容如下

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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
#include "fvCFD.H"

#include "IOmanip.H" // 输出格式控制

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"

    /*
    * About createMesh.H
    */
    // #include "createMesh.H"
    // 展开 createMesh.H 的主要内容
    Info<< "Default region: " << fvMesh::defaultRegion << endl;
    Info<< "Mesh sudirectory: " << fvMesh::meshSubDir << endl;
    
    Foam::autoPtr<Foam::fvMesh> meshPtr(nullptr);
    Foam::word regionName(Foam::polyMesh::defaultRegion);

    Foam::Info << "Create mesh";
    Foam::Info << " for time = " << runTime.timeName() << Foam::nl;

    meshPtr.reset
    (
        new Foam::fvMesh
        (
            Foam::IOobject
            (
                regionName,
                runTime.timeName(),
                runTime,
                Foam::IOobject::MUST_READ
            ),
            false
        )
    );
    meshPtr().init(true);
    Foam::fvMesh& mesh = meshPtr();
    Foam::Info << Foam::endl;
    // 见上一小节的解释
    

    Info<< "Mesh directory: " << mesh.meshDir() << endl; // methods from polyMesh
	// 输出mesh存放的文件夹路径
	
    const pointField& p = mesh.points(); // mesh对象里保存的点
    Info<< "Number of points: " << p.size() << endl;
    for (int i=0; i<3; ++i) // 显示前3个点
    {
        Info<< "("
            << setf(ios_base::scientific)
            << setw(15)
            << p[i][0]; // 第i点的坐标分量1
        Info<< ", "
            << setf(ios_base::scientific)
            << setw(15)
            << p[i][1]; // 第i点的坐标分量2
        Info<< ", "
            << setf(ios_base::scientific)
            << setw(15)
            << p[i][2] // 第i点的坐标分量3
            << ")" << endl;
    }

    const faceList& f = mesh.faces(); // mesh对象里的面
    Info<< "Number of faces: " << f.size() << endl;
    forAll(f,i) // OpenFOAM提供的遍历方法
    {
        Info<< "("
            << setw(6) << f[i][0] // 组成第i面的第1个点的序号
            << setw(6) << f[i][1] // 组成第i面的第2个点的序号
            << setw(6) << f[i][2] // 组成第i面的第3个点的序号
            << setw(6) << f[i][3] // 组成第i面的第4个点的序号
            << ")" << endl;
    }
    // 可以打开 debug_case/constant/polyMesh/faces 比较

    const labelList& fOwner = mesh.faceOwner();
    Info<< "Number of face owner: " << fOwner.size() << endl;
    forAll(fOwner,i)
    {
        Info<< setw(6) << fOwner[i] << endl;
    }
    // 输出自然列表序号的面的owner的单元序号
	// 例如输出
	// 0                    // 0号面的owner是单元0
	// 0                    // 1号面的owner是单元0
	// 1                    // 2号面的owner是单元1
	// 2                    // 3号面的owner是单元2
	// ...
	// 以此类推

    const labelList& fNeigh = mesh.faceNeighbour();
    Info<< "Number of face neighbour: " << fNeigh.size() << endl;
    forAll(fNeigh,i)
    {
        Info<< setw(6) << fNeigh[i] << endl;
    }
    // 输出自然列表序号的面的neighbour的单元序号
    // 例如
    // 1                   // 0号面的neighbour是单元1
    // 2                   // 1号面的neighbour是单元2
    // 3                   // 2号面的neighbour是单元3
    // 3                   // 3号面的neighbour是单元3
    // 其他的面不是任何单元的neighbour

    const polyBoundaryMesh& bm = mesh.boundaryMesh();
    Info<< "Number of boundary mesh: " << bm.size() << endl;
    forAll(bm,i)
    {
        Info<< "Boundary name: " << bm[i].name()
            << "\tBoundary type: " << bm[i].type()
            << endl;
    }
    // 输出设置的边界


    Info<< nl << endl;
    
    Info<< "Bounding box: " << mesh.bounds() << endl;
    Info<< "Mesh volume: " << sum(mesh.V()).value() << endl;

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< nl;
    runTime.printExecutionTime(Info);

    Info<< "End\n" << endl;

    return 0;
}


// ************************************************************************* //

为了方便结果显示，修改调试算例 debug_case/system/blockMeshDict 中设置的网格划分数量

1
2
3
4
5
6
...
blocks
(
    hex (0 1 2 3 4 5 6 7) (2 2 1) simpleGrading (1 1 1)
);
...

编译运行

结果如下

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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
Create time

Default region: region0
Mesh sudirectory: polyMesh
Create mesh for time = 0

Mesh directory: "polyMesh"
Number of points: 18
(   0.000000e+00,    0.000000e+00,    0.000000e+00)
(   5.000000e-02,    0.000000e+00,    0.000000e+00)
(   1.000000e-01,    0.000000e+00,    0.000000e+00)
Number of faces: 20
(     1     4    13    10)
(     3    12    13     4)
(     4    13    14     5)
(     4     7    16    13)
(     6    15    16     7)
(     7    16    17     8)
(     0     9    12     3)
(     3    12    15     6)
(     2     5    14    11)
(     5     8    17    14)
(     0     1    10     9)
(     1     2    11    10)
(     0     3     4     1)
(     3     6     7     4)
(     1     4     5     2)
(     4     7     8     5)
(     9    10    13    12)
(    12    13    16    15)
(    10    11    14    13)
(    13    14    17    16)
Number of face owner: 20
     0
     0
     1
     2
     2
     3
     0
     2
     1
     3
     0
     1
     0
     2
     1
     3
     0
     2
     1
     3
Number of face neighbour: 4
     1
     2
     3
     3
Number of boundary mesh: 3
Boundary name: movingWall	Boundary type: wall
Boundary name: fixedWalls	Boundary type: wall
Boundary name: frontAndBack	Boundary type: empty


Bounding box: (0.000000e+00 0.000000e+00 0.000000e+00) (1.000000e-01 1.000000e-01 1.000000e-02)
Mesh volume: 1.000000e-04

ExecutionTime = 0.000000e+00 s  ClockTime = 0.000000e+00 s

End


网格类使用示例

让我们尝试使用更多的网格类方法（成员函数）。

应用准备

1
2
3
4
5
// terminal
foamNewApp 04_02_mesh
cd 04_02_mesh
cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity debug_case
code .

脚本和说明

使用上节相同脚本内容即可，略

主源码

主源码如下

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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
#include "fvCFD.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"

    #include "createMesh.H"

    Info<< "Time                    : " << runTime.timeName() << nl
        << "Number of mesh cells    : " << mesh.C().size() << nl
        << "Number of internal faces: " << mesh.Cf().size() << nl
        << endl;

    for (label cellI=0; cellI<mesh.C().size(); ++cellI) // C++原生方式遍历
    {
        if (cellI%100 == 0)
        {
            Info<< "Cell " << cellI 
                << " with center at " << mesh.C()[cellI] // 返回单元中心坐标
                << endl;
        }
    }
    Info<< nl << endl;
    

    forAll(mesh.owner(),faceI) // OpenFOAM方式遍历
    {
        if (faceI%200 == 0)
        {
            Info<< "Internal face " << faceI
                << " with center at " << mesh.Cf()[faceI] // 返回面中心坐标
                << " with owner cell " << mesh.owner()[faceI] // 返回单元序号
                << " and neighbour cell " << mesh.neighbour()[faceI] // 返回单元序号
                << endl;
        }
    }
    Info<< nl << endl;

    forAll(mesh.boundaryMesh(), patchI)
    {
        Info<< "Patch " << patchI
            << " is " << mesh.boundary()[patchI].name() // 返回边界名称
            << " with " << mesh.boundary()[patchI].Cf().size() << " faces. "
            // 返回边界单元面的数量
            << "Start from face " << mesh.boundary()[patchI].start()
            // 返回边界其实面序号
            << endl;
    }
    Info<< nl << endl;

    label nIndex(0); // 自己设置一个序号
    forAll(mesh.boundaryMesh(), patchI) // 边界中遍历
    {
        Info<< "Patch " << patchI << nl
            << "\tits face " << nIndex
            << " adjacent to cell " << mesh.boundary()[patchI].patch().faceCells()[nIndex] << nl
            // 返回相邻单元序号
            << "\tits normal vector " << mesh.boundary()[patchI].Sf()[nIndex] << nl
            // 返回边界中面的面法向量
            << "\tits surface area " << mag(mesh.boundary()[patchI].Sf()[nIndex]) << nl
            // 返回面大小
            << endl;
    }
    Info<< nl << endl;


    const faceList& fcs = mesh.faces();
    const pointField& pts = mesh.points(); // 也可以使用 List<point>& 类型
    const List<point>& cents = mesh.faceCentres(); // 所有面心

    forAll(fcs,faceI) // 在所有面中遍历
    {
        if (faceI%200 == 0)
        {
            if (faceI < mesh.Cf().size()) // 如果是内部面
            {
                Info<< "Internal face ";
            } else // 否则是边界面
            {
                forAll(mesh.boundary(),patchI)
                {
                    if ( (mesh.boundary()[patchI].start() <= faceI) &&
                        (faceI < mesh.boundary()[patchI].start() + 
                        mesh.boundary()[patchI].Cf().size()))
                        // 如果在此边界内
                    {
                        Info<< "Face on patch " << patchI << ", faceI ";
                        break;
                    }
                }
            }
            Info<< faceI << " with centre at " << cents[faceI] // 返回面心坐标
                << " has " << fcs[faceI].size() << " vertices: "; // 返回面节点数量
            forAll(fcs[faceI],vertexI) // 在该面遍历
            {
                Info<< " " << pts[fcs[faceI][vertexI]]; // 输出节点坐标
            }
            Info<< endl;
        }
    }
    Info<< nl << endl;

	// 一般来说， empty 的边界需要特别注意，有时候需要专门处理
	// 需要找出 empty 的边界的面
    forAll(mesh.boundaryMesh(),patchI) // 遍历组成边界的所有面
    {
        const polyPatch& pp = mesh.boundaryMesh()[patchI];
        if (isA<emptyPolyPatch>(pp)) // 如果面是empty的
        {
            Info<< "Patch " << patchI
                << ": " << mesh.boundary()[patchI].name() // 返回边界名称
                << " is empty."
                << endl;
        }
    }
    Info<< nl << endl;

    word patchName("movingWall");
    label patchID = mesh.boundaryMesh().findPatchID(patchName);
    // 找到符合名称的边界面
    Info<< "Retrived patch " << patchName
        << " at index " << patchID
        << " using its name only."
        << endl;
    Info<< nl << endl;
	

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< nl;
    runTime.printExecutionTime(Info);

    Info<< "End\n" << endl;

    return 0;
}


// ************************************************************************* //

编译运行

结果如下

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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
Create time

Create mesh for time = 0

Time                    : 0
Number of mesh cells    : 400
Number of internal faces: 760

Cell 0 with center at (0.0025 0.0025 0.005)
Cell 100 with center at (0.0025 0.0275 0.005)
Cell 200 with center at (0.0025 0.0525 0.005)
Cell 300 with center at (0.0025 0.0775 0.005)


Internal face 0 with center at (0.005 0.0025 0.005) with owner cell 0 and neighbour cell 1
Internal face 200 with center at (0.0125 0.03 0.005) with owner cell 102 and neighbour cell 122
Internal face 400 with center at (0.03 0.0525 0.005) with owner cell 205 and neighbour cell 206
Internal face 600 with center at (0.0375 0.08 0.005) with owner cell 307 and neighbour cell 327


Patch 0 is movingWall with 20 faces. Start from face 760
Patch 1 is fixedWalls with 60 faces. Start from face 780
Patch 2 is frontAndBack with 0 faces. Start from face 840


Patch 0
	its face 0 adjacent to cell 380
	its normal vector (0 5e-05 0)
	its surface area 5e-05

Patch 1
	its face 0 adjacent to cell 0
	its normal vector (-5e-05 0 0)
	its surface area 5e-05

Patch 2
	its face 0 adjacent to cell 0
	its normal vector (0 0 -2.5e-05)
	its surface area 2.5e-05



Internal face 0 with centre at (0.005 0.0025 0.005) has 4 vertices:  (0.005 0 0) (0.005 0.005 0) (0.005 0.005 0.01) (0.005 0 0.01)
Internal face 200 with centre at (0.0125 0.03 0.005) has 4 vertices:  (0.01 0.03 0) (0.01 0.03 0.01) (0.015 0.03 0.01) (0.015 0.03 0)
Internal face 400 with centre at (0.03 0.0525 0.005) has 4 vertices:  (0.03 0.05 0) (0.03 0.055 0) (0.03 0.055 0.01) (0.03 0.05 0.01)
Internal face 600 with centre at (0.0375 0.08 0.005) has 4 vertices:  (0.035 0.08 0) (0.035 0.08 0.01) (0.04 0.08 0.01) (0.04 0.08 0)
Face on patch 1, faceI 800 with centre at (0.1 0.0025 0.005) has 4 vertices:  (0.1 0 0) (0.1 0.005 0) (0.1 0.005 0.01) (0.1 0 0.01)
1000 with centre at (0.0425 0.0025 0) has 4 vertices:  (0.04 0 0) (0.04 0.005 0) (0.045 0.005 0) (0.045 0 0)
1200 with centre at (0.0925 0.0025 0) has 4 vertices:  (0.09 0 0) (0.09 0.005 0) (0.095 0.005 0) (0.095 0 0)
1400 with centre at (0.0425 0.0025 0.01) has 4 vertices:  (0.04 0 0.01) (0.045 0 0.01) (0.045 0.005 0.01) (0.04 0.005 0.01)
1600 with centre at (0.0925 0.0025 0.01) has 4 vertices:  (0.09 0 0.01) (0.095 0 0.01) (0.095 0.005 0.01) (0.09 0.005 0.01)


Patch 2: frontAndBack is empty.


Retrived patch movingWall at index 0 using its name only.



ExecutionTime = 0 s  ClockTime = 0 s

End


计算网格体积比

我们基于前面网格方法的讨论来尝试一个计算网格参数的应用，也许以后会用在求解器开发上。

代码来自 Tom Smith 讲义《programming with OpenFOAM》（见文首参考链接）

应用准备

1
2
3
4
// terminal 
foamNewApp 04_03_mesh
cd 04_03_mesh
cp -r $FOAM_TUTORIALS/incompressible/simpleFoam/pitzDaily debug_case

脚本和说明

以后非特别说明，不再赘述。

主源码

主源码如下

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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
#include "fvCFD.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"

    #include "createMesh.H"

    const scalar c = mesh.C().size();
    Info<< "Number of cells = " << c << endl;
    // 网格单元数量

    const labelListList& neighbour = mesh.cellCells();
    // 返回邻单元为新列表

    List<scalar> ratios(0); // 声明一个scalar列表且全初始化为0
    scalar volumeRatio = 0;

    forAll (neighbour, cellI) // 遍历邻单元，其实就是所有单元
    {
        List<label> n = neighbour[cellI];
        // 对每个单元，获得它的邻单元到列表 n

        const scalar cellVolume = mesh.V()[cellI];
        // 本单元体积

        forAll (n, i) // 对于某个本单元，遍历它的邻单元
        {
            label neighbourIndex = n[i];
            scalar neighbourVolume = mesh.V()[neighbourIndex];
            // 邻单元体积

            if (neighbourVolume >= cellVolume) // 体积比大于 1 的情况
            {
                volumeRatio = neighbourVolume / cellVolume;
                ratios.append(volumeRatio);
                // 体积比存入 ratios 列表
                // 方便但是效率低
            }
        }
    }

    Info<< "Maximum volume ratio = " << max(ratios) << endl;
    // 最大体积比


    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< nl;
    runTime.printExecutionTime(Info);

    Info<< "End\n" << endl;

    return 0;
}

编译运行，结果为

1
2
3
4
5
6
7
8
9
10
Create time

Create mesh for time = 0

Number of cells = 12225
Maximum volume ratio = 1.49908

ExecutionTime = 0.13 s  ClockTime = 0 s

End

代码改进

List 类的 append 方法虽然方便，但是效率很低，因为每一次向 List 添加数据的时候，都要创建一个 n+1 大小的新 List ，将老数据拷贝过去，再加入新数据。在数据量很大的时候，每次添加新数据都有新建再复制的过程，效率非常低。

我们考虑到对于要处理的网格，因为网格确定，所以要计算的体积比总量也是确定的。可以提前声明一个固定大小的 List，分配好内存，后续只需要向其添加数据就行。

改进后主源码如下

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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#include "fvCFD.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"

    #include "createMesh.H"

    const scalar c = mesh.C().size();
    Info<< "Number of cells = " << c << endl;

    const labelListList& neighbour = mesh.cellCells();

	/*
	* 代码改进如下 
	*/
    label len = mesh.Cf().size(); // 需要的 ratio 的 List 大小
    scalar initial = 0;
    List<scalar> ratios(len, initial); // 创建固定大小的 List 
    label counter = 0; // ratios 的序号
    scalar volumeRatio = 0;

    forAll (neighbour, cellI) // 遍历邻单元
    {
        List<label> n = neighbour[cellI];

        const scalar cellVolume = mesh.V()[cellI];

        forAll (n, i)
        {
            label neighbourIndex = n[i];
            scalar neighbourVolume = mesh.V()[neighbourIndex];

            if (neighbourVolume >= cellVolume)
            {
                volumeRatio = neighbourVolume / cellVolume;
                ratios[counter] = volumeRatio; // 存储体积比
                counter += 1;
            }
        }
    }

    Info<< "Maximum volume ratio = " << max(ratios) << endl;

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< nl;
    runTime.printExecutionTime(Info);

    Info<< "End\n" << endl;

    return 0;
}

编译运行结果如下

1
2
3
4
5
6
7
8
9
10
Create time

Create mesh for time = 0

Number of cells = 12225
Maximum volume ratio = 1.49908

ExecutionTime = 0.06 s  ClockTime = 0 s

End

可以很明显看到运行时间从 0.13s 减少到了 0.06s，效率大大提高。

代码再改进

如果我们只是要获得最大体积比，那么其实我们不需要保存所有的体积比，保存体积比的最大值就够了。

主源码如下

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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#include "fvCFD.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"

    #include "createMesh.H"

    const scalar c = mesh.C().size();
    Info<< "Number of cells = " << c << endl;

    const labelListList& neighbour = mesh.cellCells();

	/*
	* 代码改进如下 
	*/
    scalar volumeRatio = 0.0;
    scalar currentRatio = 0.0;

    forAll (neighbour, cellI)
    {
        List<label> n = neighbour[cellI];

        const scalar cellVolume = mesh.V()[cellI];

        forAll (n, i)
        {
            label neighbourIndex = n[i];
            scalar neighbourVolume = mesh.V()[neighbourIndex];

            if (neighbourVolume >= cellVolume)
            {
                volumeRatio = neighbourVolume / cellVolume;

                if (volumeRatio > currentRatio) // 只存储最大体积比
                {
                    currentRatio = volumeRatio;
                }
            }
        }
    }

    Info<< "Maximum volume ratio = " << currentRatio << endl;

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< nl;
    runTime.printExecutionTime(Info);

    Info<< "End\n" << endl;

    return 0;
}

编译运行结果如下

1
2
3
4
5
6
7
8
9
10
Create time

Create mesh for time = 0

Number of cells = 12225
Maximum volume ratio = 1.49908

ExecutionTime = 0.05 s  ClockTime = 0 s

End

可以看到运行时间进一步减少，效率又提高了。

代码最终版

仅获得体积比还不够，我们可能希望知道测试算例的网格有多少超过了指定的体积比准则。

我们通过 OpenFOAM 的字典来接入体积比准测这一参数。

为测试算例提供字典 /userApp/debug_case/system/volumeRatioDict ，内容如下

我们约定

  • *Properties 文件放在 /userApp/constant 文件夹下
  • *Dict 文件放在 /userApp/system 文件夹下
1
2
3
4
5
6
7
8
9
10
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

maxRatio 1.2;

主源码如下

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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
#include "fvCFD.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"

    #include "createMesh.H"

    IOdictionary volumeRatioDict // 建立字典对象
    (
        IOobject
        (
            "volumeRatioDict", // 文件名称
            runTime.system(), // 文件路径
            mesh, // 基于 mesh 构造
            IOobject::MUST_READ, // 必读
            IOobject::NO_WRITE // 只读不写
        )
    );

    const scalar c = mesh.C().size();
    Info<< "Number of cells = " << c << endl;

    const labelListList& neighbour = mesh.cellCells();

    scalar volumeRatio = 0.0;
    scalar currentRatio = 0.0;

    label nFail = 0; // 统计超出准测的数量

    scalar maxRatio(readScalar(volumeRatioDict.lookup("maxRatio")));
	// 从字典读入

    forAll (neighbour, cellI)
    {
        List<label> n = neighbour[cellI];

        const scalar cellVolume = mesh.V()[cellI];

        forAll (n, i)
        {
            label neighbourIndex = n[i];
            scalar neighbourVolume = mesh.V()[neighbourIndex];

            if (neighbourVolume >= cellVolume)
            {
                volumeRatio = neighbourVolume / cellVolume;

                if (volumeRatio > currentRatio)
                {
                    currentRatio = volumeRatio;
                }

                if (volumeRatio > maxRatio) // 统计超出准测的数量
                {
                    nFail += 1;
                }
            }
        }
    }

    Info<< "Maximum volume ratio = " << currentRatio << nl
        << "Number of cell volume ratios exceeding " << maxRatio
        << " = " << nFail << nl
        << endl;

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< nl;
    runTime.printExecutionTime(Info);

    Info<< "End\n" << endl;

    return 0;
}

编译运行，结果如下

1
2
3
4
5
6
7
8
9
10
11
12
Create time

Create mesh for time = 0

Number of cells = 12225
Maximum volume ratio = 1.49908
Number of cell volume ratios exceeding 1.2 = 532


ExecutionTime = 0.06 s  ClockTime = 1 s

End

小结

通过前面的几篇文章，OpenFOAM 应用的常见头文件 setRootCase.H ，createTime.H 和 createMesh.H 已经尽可能在不过度深入的情况下进行了讨论。此外，我们也尝试讨论了一个基于离散网格的完整应用，不断改进的过程对我们以后开发自己的求解器也有一些启发。

除了主函参数、时间和离散网格，我们怎么创建以及计算控制方程方程中的物理量呢？