【求教】forAll() boundaryMesh() boundary() 等的理解

大家新年!

我在学习一段代码时,遇到了一些理解上的问题。将问题现总结如下,并在文后附上源代码:

  1. forAll() 的使用问题

这个函数接收什么参数,其参数类型是什么,该函数返回什么。从我找到的资料来看,该函数是一个宏定义,

forAll(list,i)和for(i=0;i<(list).size();i++)是等价的

那是否能说forAll接收两个参数,一个是列表,一个是label类型的计数器,其作用是对该列表进行遍历?

  1. 关于 mesh.boudaryMesh() 和 mesh.boundary() 的区别

按照1.的推断,结合代码

forAll(mesh.boudaryMesh(),patchI) 及forAll(mesh.boudary(),patchI) 

是否能说mesh.boudaryMesh()和mesh.boundary() 的返回值是列表?
分别是关于什么数据的列表?

产生这个疑问的原因还在于我看到过这样的代码

// Let us define a vector whose values will not change over the course of the program execution.
	const vector originVector(0.05,0.05,0.005);

	// Calculate the distance from the origin to the cell centre furthest away.
	// In Python, this is equivalent to:
	// np.sqrt(np.sum((x0-x)**2))
	// The .value() method is called to convert from a dimensionedScalar to a regular scalar.
	const scalar rFarCell = max( // find the maximum value from all distances
		// compute distance of each cell centre from x0; units of mesh.C() are those of length, as this field
		// describes position in the Cartesian reference frame.
		mag(dimensionedVector("x0",dimLength,originVector)-mesh.C())
		).value(); // convert to dim-less scalar

如果dimensionedVector()产生的是一个带单位的矢量,那么是不是意味着在OF里,一个矢量能和一个由中心坐标矢量组成的列表进行加减运算,即一个矢量能和多个矢量进行加减运算?抑或是我关于mesh.C()返回值为由中心坐标矢量组成的列表的猜想是错误的?

  1. 关于mesh.C()和mesh.Cf()

这两个函数的作用是否返回网格的中心点坐标组成的列表 和 网格内部面中心坐标组成的列表

  1. 关于OF中scalar(标量)和field(场)的理解

两者的区别是否在于:标量在空间计算域的分布是连续的,而场的分布是离散的,在OF中当一个量被指定为场时,只能在网格上某点(比如网格中心或者网格某个面的中心)上获取它,也即必须先划分网格,才能在此基础上定义某个物理量的场。出现这个疑问是因为发现OF中定义了很多场,既有标量场(scalarField),比如压强场p,也有标量(比如nu)。

  1. 作为一个代码初学者,希望各位前辈能给一些建议。

比如我应该去什么网站查找关于forAll()的参数及返回值或者其用法。

  1. 一个小建议。希望论坛的代码区域能自动生成行号,以方便提问时引用行号。

源码如下

#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<<runTime.timeName() << endl;

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

    return 0;
}

非常感谢你关于代码行号的建议。由于论坛源代码还没有比较好的插件,代码行号的显示还在测试中。

1 个赞

forAll() 的定义位于 src/foam/containers/Lists/UList/UList.H,源码如下

#define forAll(list, i) \
    for (Foam::label i=0; i<(list).size(); i++)

因此,只要是列表,都可以使用 forAll 进行遍历。需要注意的是,forAll 只是一个宏,代替的是 for( ) 循环,所以,它并不返回什么参数。而通过这个,可以获得 list 大小和 iiforAll() 之前无须定义。

1 个赞

对于mesh.C() 和 mesh.Cf() ,你的解释是对的。

mesh.C() 返回整个网格各 cell 中心点的坐标。
mesh.Cf() 返回整个网格内部各面的中心坐标。

1 个赞

boundaryMesh() 位于 src/OpenFOAM/meshes/polyMesh/polyMesh.H

        //- Return boundary mesh
        const polyBoundaryMesh& boundaryMesh() const
        {
            return boundary_;
        }

mesh.boundaryMesh() 是边界的列表,例如,其中一项 mesh.boundaryMesh()[0] 的输出是

type            wall;
nFaces          5;
startFace       40;

boundary() 位于 src/finiteVolume/fvMesh/fvMesh.C

const Foam::fvBoundaryMesh& Foam::fvMesh::boundary() const
{
    return boundary_;
}

也是边界的列表,由于类型不同,因此对于数据处理的方式也不太一样。直接输出 boundary()[0] 会报错。

理解的话,从这里来看比较好。Weller 说过[1],

References

  1. CFD Online: FvMesh

Further reading

  1. OpenFOAM Programming Tips | PPT

对于 scalarscalarField 区别:

  • scalar 是一个标量,与网格无关。
  • scalarField 是一组标量,跟网格有对应关系。体心标量场类型为 volScalarField
  • Field 就是场,与网格有关,可以有标量场 volScalarField,也可以有向量场 volVectorField

代码学习,可以看这里 https://cpp.openfoam.org/ 最全了。还有就是网上多搜索下。也可以来这里多多讨论呢。

还可以关注下:sourceflux.dewolfdynamics.com