003-OpenFOAM的网格及操作

从网格谈起，简单讲述OpenFOAM中如何植入一般的C++代码，如何使用C++ STL和标准库，以及OpenFOAM的forAll语法

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include <vector>
#include <algorithm>
#include <numeric>

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

    // These two create the time system (instance called runTime) and fvMesh (instance called mesh).
    #include "createTime.H"
    #include "createMesh.H"

    // runTime and mesh are instances of objects (or classes).
    // If you are not familiar with what a class or object is, it is HIGHLY RECOMMENDED you visit this
    // website and only come back once you've read everything about classes, inheritance and polymorphism:
    // http://www.cplusplus.com/doc/tutorial/classes/
    // Note how the next lines call functions .timeName(), .C() and .Cf() implemented in the objects.
    // It is also important to realise that mesh.C() and .Cf() return vector fields denoting centres of each
    // cell and internal face.
    // Calling the mesh.C().size() method therefore yields the total size of the mesh.
    Info << "Hello there, the most recent time folder found is " << runTime.timeName() << nl
         << "The mesh has " << mesh.C().size() << " cells and " << mesh.Cf().size() << nl
         << " internal faces in it. Wubalubadubdub!" << nl << endl;
         /*
            ---runTime.timeName() 返回一个scalar，默认值为最后的时间步长的名字
            ---mesh.C(), mesh.Cf 都是返回一个向量场，前者返回网格的中心，后者返回每个面的中心，但是注意
               .Cf() 只是返回内部网格的面的中心。
         */

    // It's possible to iterate over every cell in a standard C++ for loop
    /*
        OpenFOAM里当然也可以使用标准C++的库及其内部函数，但是要记住，不光需要在头文件中include
        在使用对应的函数及方法时，要使用std::命名空间，否则函数将从Foam::中载入
    */
    std::vector<double> cellX ;
    std::vector<double> cellY ;
    std::vector<double> cellY ;
    for (label cellI = 0; cellI < mesh.C().size(); cellI++)     // label是对C++中 int类的继承
    {
        if (cellI%20 == 0){
            Info << "Cell " << cellI << " with centre at " << mesh.C()[cellI] << endl;
            /* Info << "Cell " << cellI << " 's centre x: " << mesh.C()[cellI][0] << endl;
            Info << "Cell " << cellI << " 's centre y: " << mesh.C()[cellI][1] << endl;
            Info << "Cell " << cellI << " 's centre z: " << mesh.C()[cellI][2] << endl; */
        } // only show every twentieth cell not to spam the screen too much
        cellX.push_back(double(mesh.C()[cellI][0]));
        cellY.push_back(double(mesh.C()[cellI][1]));
        cellZ.push_back(double(mesh.C()[cellI][2]));
    }
        
    std::vector<double> center={0,0,0};
    center[0]=std::accumulate(std::begin(cellX),std::end(cellX),0.0) / cellX.size();
    center[1]=std::accumulate(std::begin(cellY),std::end(cellY),0.0) / cellY.size();
    center[2]=std::accumulate(std::begin(cellZ),std::end(cellZ),0.0) / cellZ.size();
    Info << "The center of the cell is: " << center[0] <<' '<< center[1]<<' '<< center[2] <<nl;
    Info << endl; // spacer

    // Each cell is constructed of faces - these may either be internal or constitute a
    // boundary, or a patch in OpenFOAM terms; internal faces have an owner cell
    // and a neighbour.
    // mesh.owner()与mesh.Cf()的对象都是内部网格的面，mesh.owner()返回一个列表，内容是所有网格的序号，
    // mesh.Cf()返回的列表，列表长度与mesh.owner()相同，元素内容是面中心点的坐标
    // mesh.boundary()[Optional_patchI].Cf()等等返回的内容就是边界列表
    // neighbour()的概念很好理解，就是旁边的那个网格，不一定序号就是+1，因为OpenFOAM编号是一层一层的
    for (label faceI = 0; faceI < mesh.owner().size(); faceI++)
        if (faceI%40 == 0)
            Info << "Internal face " << faceI << " with centre at " << mesh.Cf()[faceI]
                 << " with owner cell " << mesh.owner()[faceI]
                 << " and neighbour " << mesh.neighbour()[faceI] << endl;
    Info << endl;

    // 边界条件可以通过边界网格对象访问。
    // 在现实中，每个边界面也包括在"constant/polyMesh/face"中
    // 描述。但是，在该文件中，内部面首先被定义。
    // 另外，"constant/polyMesh/boundary"文件定义了起始faceI
    // 边界面定义开始的索引。
    // OpenFOAM还为所有条目上的for循环提供了宏定义（forAll）
    // 在字段或列表中，这样可以节省输入量。
    forAll(mesh.boundaryMesh(), patchI) 
    // 等价于 for (label patchI = 0; patchI < mesh.boundaryMesh().size(); patchI++)
        Info << "Patch " << patchI << ": " << mesh.boundary()[patchI].name() << " with "
             << mesh.boundary()[patchI].Cf().size() << " faces. Starts at total face "
             << mesh.boundary()[patchI].start() << endl;  // .start()返回scalar，值为开始的面的编号
    Info << endl;


    // Faces adjacent to boundaries may be accessed as follows.
    // Also, a useful thing to know about a face is its normal vector and face area.
    // 对于内部面，方法. Sf()可以直接在网格实例上调用。
    // 此外，还有一个快捷方法.magSf()，它返回表面面积标量。
    // 对于内面，法向量从所有者指向邻居而主人的cellI指数比邻居的要小。对于边界面，
    // 法线总是指向域外(它们有“想象的”邻居这些都不存在)可以更详细地观察组成每个面的点。
    // 首先，我们通过获取对各自的引用来定义一些简写网格中的对象。这些被定义为常数，因为我们不打算以任何方式改变网格。
    // 注意:这些列表指的是网格的物理定义，因此包括边界面。可以使用mesh.boundary()[patchI].Cf().size()
    // 和mesh.boundary()[patchI].start()方法检查面是否在内部
    // 或者在边界上

    label patchFaceI(0);
    forAll(mesh.boundaryMesh(), patchI)
        Info << "Patch " << patchI << " has its face " << patchFaceI << " adjacent to cell "
             << mesh.boundary()[patchI].patch().faceCells()[patchFaceI]
             << ". It has normal vector " << mesh.boundary()[patchI].Sf()[patchFaceI]
             << " and surface area " << mag(mesh.boundary()[patchI].Sf()[patchFaceI])
             << endl;
    Info << endl;

    const faceList& fcs = mesh.faces();
    const List<point>& pts = mesh.points(); // List 是对vector的继承
    const List<point>& cents = mesh.faceCentres();

    forAll(fcs,faceI)
    // 等价于 for (label faceI = 0; faceI < fcs.size(); faceI++)
        if (faceI%80==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; // exit the forAll loop prematurely
                    }
            }

            Info << faceI << " with centre at " << cents[faceI]
                 << " has " << fcs[faceI].size() << " vertices:";
            forAll(fcs[faceI],vertexI)
                // Note how fcs[faceI] holds the indices of points whose coordinates
                // are stored in the pts list.
                Info << " " << pts[fcs[faceI][vertexI]];
            Info << endl;
        }
    Info << endl;

    // “empty”边界类型是一个特别的情况，可能会导致意外的行为，因为它的.Cf()方法返回的列表大小为0。
    // 可以使用isA<emptyPolyPatch>(polyPatch_to_be_Checked)检查补丁的类型，以避免遇到此问题

    label patchID(0);
    const polyPatch& pp = mesh.boundaryMesh()[patchID];
    if (isA<emptyPolyPatch>(pp))
    {
        // patch patchID is of type "empty".
        Info << "You will not see this." << endl;
    }

    // Patches may also be retrieved from the mesh using their name. This could be
    // useful if the user were to refer to a particular patch from a dictionary
    // (like when you do when calculating forces on a particular patch).
    // 使用patch的名称检索网格时，注意要创建一个word变量保存，findPatchXX接收的变量类型为word，不可以是string
    word patchName("movingWall");
    patchID = mesh.boundaryMesh().findPatchID(patchName);
    Info << "Retrieved patch " << patchName << " at index " << patchID << " using its name only." << nl << endl;

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

    return 0;
}


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