数值风云

OpenFOAM开发教程系列:第三讲_认识网格

OFtutorial03_understandingTheMesh

第三讲:理解网格

代码:

#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"

作用:

检查算例结构是否完整,创建时间系统(实例名为:runTime)和fvmesh(实例名为: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()
     << " internal faces in it. Wubalubadubdub!" << nl << endl;

作用:

runTime和mesh是对象(或者类)的实例。

要调用的子函数 .timeName(), .C(), .Cf()在这些类中部署。

很重要的一点是:mesh.C()以及mesh.Cf()返回一个向量场,代表了每个网格或者面的中心。

调用mesh.C().size()方法也可以返回mesh的大小。

注意:mesh.Cf()指internal face(内部网格)

代码:

for (label cellI = 0; cellI < mesh.C().size(); cellI++)
    if (cellI%20 == 0) // only show every twentieth cell not to spam the screen too much
        Info << "Cell " << cellI << " with centre at " << mesh.C()[cellI] << endl;
Info << endl; // spacer

作用:

可以通过C++标准for循环来遍历每一个网格。

代码:

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;

作用:

每个cell都是由面组成的。这些面可能是OpenFOAM中的内部面,边界面或者patch。内部面有一个主网格还有一个邻居网格。

代码:

forAll(mesh.boundaryMesh(), patchI)
    Info << "Patch " << patchI << ": " << mesh.boundary()[patchI].name() << " with "
         << mesh.boundary()[patchI].Cf().size() << " faces. Starts at total face "
         << mesh.boundary()[patchI].start() << endl;
Info << endl;

作用:

边界条件可以通过boundaryMesh对象获取。

实际上,每个边界面的描述都包含在了constant/polyMesh/faces文件中。但是,在这个文件中,内部面被首先定义。

此外,在constant/polyMesh/boundary文件中定义了起始面的ID,表明边界文件的起始位置。OpenFOAM也为循环提供宏定义,这可以减少很多代码量。

代码:

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();
const List<point>& cents = mesh.faceCentres();

forAll(fcs,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;

作用:

对于内部的面来说,.Sf()方法可直接被mesh实例调用。

在计算面的面积的时候有一个更简短的方法.magSf(),这个方法返回面积的数量值。

对于内部面,法向量从owner指向neighbor,且owner的索引比neighbor小。对于边界面,面的法向量指向计算域的外侧。

也可以查看定义面的点。

首先定义面中各种对象的引用。这些变量被定义为常量,因为我们不想去更改他们。

注意:

这些lists代表了面物理上的定义,所以包含了变截面。用户可以通过mesh.boundary()[patchI].Cf().size()mesh.boundary()[patchI].start()方法来判断面试边界面还是内部面。

代码:

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;
}

作用:

在原始的cavity教程中,frontAndBack边界被定义为“empty”类型。这是一种很特殊的边界条件,会引起很多意料之外的问题。比如说.Cf()会尺寸为0.

可能需要检测patch的类型以避免类似的问题。

代码:

word patchName("movingWall");
patchID = mesh.boundaryMesh().findPatchID(patchName);
Info << "Retrieved patch " << patchName << " at index " << patchID << " using its name only." << nl << endl;

作用:

patches也可以通过他们的名称访问。当用户想从dictionary中引用特定的patch时,这可能会被用到。(比如:计算特定的patch的受力)

总体代码:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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"

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()
         << " internal faces in it. Wubalubadubdub!" << nl << endl;

    // It's possible to iterate over every cell in a standard C++ for loop
    for (label cellI = 0; cellI < mesh.C().size(); cellI++)
        if (cellI%20 == 0) // only show every twentieth cell not to spam the screen too much
            Info << "Cell " << cellI << " with centre at " << mesh.C()[cellI] << endl;
    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.
    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;

    // Boundary conditions may be accessed through the boundaryMesh object.
    // In reality, each boundary face is also included in the constant/polyMesh/faces
    // description. But, in that file, the internal faces are defined first.
    // In addition, the constant/polyMesh/boundary file defines the starting faceI
    // indices from which boundary face definitions start.
    // OpenFOAM also provides a macro definition for for loops over all entries
    // in a field or a list, which saves up on the amount of typing.
    forAll(mesh.boundaryMesh(), patchI)
        Info << "Patch " << patchI << ": " << mesh.boundary()[patchI].name() << " with "
             << mesh.boundary()[patchI].Cf().size() << " faces. Starts at total face "
             << mesh.boundary()[patchI].start() << endl;
    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.
    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;

    // For internal faces, method .Sf() can be called directly on the mesh instance.
    // Moreover, there is a shorthand method .magSf() which returns the surface area
    // as a scalar.
    // For internal faces, the normal vector points from the owner to the neighbour
    // and the owner has a smaller cellI index than the neighbour. For boundary faces,
    // the normals always point outside of the domain (they have "imaginary" neighbours
    // which do not exist).

    // It is possible to look at the points making up each face in more detail.
    // First, we define a few shorthands by getting references to the respective
    // objects in the mesh. These are defined as constants since we do not aim to
    // alter the mesh in any way.
    // NOTE: these lists refer to the physical definition of the mesh and thus
    // include boundary faces. Use can be made of the mesh.boundary()[patchI].Cf().size()
    // and mesh.boundary()[patchI].start() methods to check whether the face is internal
    // or lies on a boundary.
    const faceList& fcs = mesh.faces();
    const List<point>& pts = mesh.points();
    const List<point>& cents = mesh.faceCentres();

    forAll(fcs,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;

    // In the original cavity tutorial, on which the test case is based,
    // the frontAndBack boundary is defined as and "empty" type. This is a special
    // BC case which may cause unexpected behaviour as its .Cf() field has size of 0.
    // Type of a patch may be checked to avoid running into this problem if there
    // is a substantial risk that an empty patch type will appear
    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).
    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;
}


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




5 个赞