吾生有涯 学海无涯
析模有界 知识无界
        找回密码

OpenFOAM编程案例｜03 访问网格信息

2022-03-23 分类：OpenFOAM编程 阅读() 评论(0)

本文描述利用fvMesh类获取网格信息的一些基本用法。

OpenFOAM应用过程中需要频繁地与网格进行交互，比如获取网格面信息，获得网格中心坐标等。利用fvMesh类中的相应成员函数很容易实现这类需求。

fvMesh类是一个极为复杂的类，其继承自polyMesh、lduMesh、fvSchemes、surfaceInterpolation、fvSolution及data类，因此不仅可以访问网格，还可以访问与网格有关的各种算法格式。其中对网格信息的访问功能（如访问网格面、得到网格边界信息等）主要继承自polyMesh类。

1 准备文件

与前面的案例一样，本案例只需一个源文件及Make文件夹，利用下面的命令创建文件结构。

cd $FOAM_RUN
mkdir demo3
cd demo3
touch demo3.C
mkdir Make
cd Make
touch files
touch options
cd ..

文件结构如下图所示。

  • 修改files文件的内容
demo3.C
EXE = demo3
  • 修改options文件的内容
EXE_INC = 
        -I$(FOAM_SRC)/finiteVolume/lnInclude 
        -I$(FOAM_SRC)/meshTools/lnInclude
 
EXE_LIBS = 
        -lfiniteVolume 
        -lmeshTools

2 测试文件

本案例涉及到网格操作，因此需要准备网格文件。简单起见，使用一个OpenFOAM案例库中的文件进行测试，如使用cavity。

cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity .
cd cavity
blockMesh

执行blockMesh生成网格后可以查看网格信息，如下图所示。

看到网格数量为400，网格面数为1640。

3 源代码

3.1 读取网格数量

利用fvMesh类的C()及Cf()函数可以获取网格体中心坐标列表及内部面中心列表。

编辑源文件demo3.C，输入下面的代码：

#include "fvCFD.H"
 
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
// 在头文件createMesh.H中创建了一个名为mesh的fvMesh对象
#include "createMesh.H"
 
    // 输出当前网格信息
    // 利用mesh.C()得到网格单元体心坐标列表，mesh.Cf()得到网格内部面中心坐标列表
    Info << "当前时间文件：" << runTime.timeName() << nl
         << "网格数量：" << mesh.C().size() << nl
         << "网格面数量：" << mesh.Cf().size() << nl << nl;
 
    // 可以用循环输出一些网格坐标信息
    for (int cellI = 0; cellI < mesh.C().size(); cellI++)
    {
        // 由于网格数量太多，这里做一些省略,间隔20输出一个
        if (cellI % 20 == 0)
        {
            Info << "Cell[" << cellI << "]的中心坐标为：" << mesh.C()[cellI] << endl;
        }
    }
    Info << endl;
    return 0;
}

编译上面的程序并运行。

注意需要到测试路径（本案例为cavity文件夹）下执行此程序。运行结果如下图所示。

可以看到，程序正确地获取了网格数量与网格面的数量，并输出了网格中心坐标。

3.2 得到owner及neighbour的信息

若想要得到网格面的owner及neighbour网格索引，可以使用mesh.owner()及mesh.neighbour()函数。

在前面的代码中添加下面内容：

    // 可以使用mesh.owner()及mesh.neighbour()得到owner及neightbour的网格编号
    for (label faceI = 0; faceI < mesh.owner().size(); faceI++)
    {
        if (faceI % 100 == 0)
        {
            Info << "内部面[" << faceI << "]的中心坐标为：" << mesh.C()[faceI] << nl
                 << "owner[" << faceI << "]为：" << mesh.owner()[faceI] << nl
                 << "neighbour[" << faceI << "]为：" << mesh.neighbour()[faceI] << nl;
        }
    }
    Info << endl;

编译执行，如下图所示。

3.3 获取边界面信息

利用mesh.boundaryMesh()函数可以获取边界面上的信息。

添加以下代码。

// 利用mesh.boundaryMesh()函数获取边界面信息
// forAll是OpenFOAM定义的循环宏
forAll(mesh.boundaryMesh(), patchI)
{
    Info << "Patch[" << patchI << "]：" << mesh.boundary()[patchI].name() << "包括"
        << mesh.boundary()[patchI].Cf().size() << "个网格面" << nl
        << "起始面为：" << mesh.boundary()[patchI].start() << endl;
}
Info << endl;

编译并运行，结果如下图所示。

除了frontAndBack的网格面数外，其他边界的信息都是没有问题的。至于为什么会有问题，因为该边界类型为empty导致，解决方法在文末。

3.4 边界网格面面积

若想要获取边界面上某个网格的面积，可以使用函数sf()，如下面的代码：

// 利用函数faceCells()得到网格面相邻的网格
// 利用函数sf()可以得到网格面的面积向量，该向量的模即为面积
//定义一个索引，仅为测试使用，如果要得到几何面积，则需要将该边界上的所有网格面积求和
label patchFaceI(0);
forAll(mesh.boundaryMesh(), patchI)
{
    Info << "Patch[" << patchI << "]的网格面" << patchFaceI << "与网格["
        << mesh.boundary()[patchI].patch().faceCells()[patchFaceI] <<"]相邻,"
        << "其面积向量：" << mesh.boundary()[patchI].Sf()[patchFaceI] << ",其面积为："
        << mag(mesh.boundary()[patchI].Sf()[patchFaceI])
        << endl;     
}
Info << endl;

程序运行结果如下图所示。

对于内部面，可以直接在mesh实例上调用Sf()方法。若要获取网格面的面积，可以使用函数magSf()，其等同于mag(Sf())，返回一个标量面积值。

对于内部面，法向量从owner指向neightbour，并且owner的索引小于neighbour。对于边界面，法向量总是指向计算域外部。

3.5 网格面的节点信息

可以更详细地查看构成每个网格面的点。首先，我们通过获取对网格中各个对象的引用来定义一些变量。这些变量被定义为const，因为我们不打算以任何方式改变网格。注意：这些列表指的是网格的物理定义，因此包括边界面。可以使用mesh.edge()[patchi].Cf().size()和mesh.edge()[patchi].start()方法检查网格面是位于内部还是位于边界上。

如下面的代码：

// 查看网格细节
// 将变量定义为const，则此变量只能读取，而不能修改
const faceList &fcs = mesh.faces();            // 得到所有的网格面列表
const List &pts = mesh.points();        //得到所有的网格节点列表
const List ¢s = mesh.faceCentres(); //得到网格面心点列表
 
// 遍历索引能被80整除的网格面
forAll(fcs, faceI) if (faceI % 80 == 0)
{
    if (faceI < mesh.Cf().size())
    {
        Info << "内部面";
    }
    else
    {
        forAll(mesh.boundary(), patchI) if ((mesh.boundary()[patchI].start() <= faceI) &&
                                            (faceI < mesh.boundary()[patchI].start() + mesh.boundary()[patchI].Cf().size()))
        {
            Info << "网格面在边界[" << patchI << "]上，网格面[";
            break;
        }
    }
 
    Info << faceI << "]的中心为：" << cents[faceI] << "，拥有" << fcs[faceI].size() << "个节点：";
 
    forAll(fcs[faceI], vertexI)
    {
        Info << " " << pts[fcs[faceI][vertexI]];
    }
    Info << endl;
}
Info << endl;
 
return 0;
}

运行 结果如下图所示。

3.6 一些额外的问题

在本测试案例中，FrontAndBack边界被定义为empty类型。这是一种特殊的边界，在使用上面的函数时会出现一些意想不到的问题，如调用Cf函数时会返回空列表，导致size()函数的值为零。为了避免出现此类问题，在进行Cf()函数调用时，可以先检查边界类型。

如下面的代码：

// 判断边界类型是否为empty类型
forAll(mesh.boundaryMesh(), patchID)
{
    const polyPatch &pp = mesh.boundaryMesh()[patchID];
    if (isA(pp))
    {
        Info << "patch[" << patchID <<"]:" << mesh.boundary()[patchID].name()
            << "的边界类型为empty" << endl;
    }
    else
    {
        Info << "patch[" << patchID << "]:"<< mesh.boundary()[patchID].name()
            <<"的边界类型不是empty" << endl;
    }
}

运行结果如下图所示。

也可以通过patch的名称利用函数findPatchID获取其ID，如下面的代码：

// 利用边界名称得到其ID
word patchName("movingWall");
label patchID = mesh.boundaryMesh().findPatchID(patchName);
Info << "patch name:" << patchName << nl << "patchID:" << patchID << endl;

运行结果如下图所示。

4 完整代码

本案例的完全程序代码如下所示。

#include "fvCFD.H"
 
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
// 在头文件createMesh.H中创建了一个名为mesh的fvMesh对象
#include "createMesh.H"
 
    // 输出当前网格信息
    // 利用mesh.C()得到网格单元体心坐标列表，mesh.Cf()得到网格内部面中心坐标列表
    Info << "当前时间文件：" << runTime.timeName() << nl
         << "网格数量：" << mesh.C().size() << nl
         << "网格面数量：" << mesh.Cf().size() << nl << nl;
 
    // 可以用循环输出一些网格坐标信息
    for (int cellI = 0; cellI < mesh.C().size(); cellI++)
    {
        // 由于网格数量太多，这里做一些省略,间隔20输出一个
        if (cellI % 20 == 0)
        {
            Info << "Cell[" << cellI << "]的中心坐标为：" << mesh.C()[cellI] << endl;
        }
    }
    Info << endl;
 
    // 可以使用mesh.owner()及mesh.neighbour()得到owner及neightbour的网格中心坐标列表
    for (label faceI = 0; faceI < mesh.owner().size(); faceI++)
    {
        if (faceI % 100 == 0)
        {
            Info << "内部面[" << faceI << "]的中心坐标为：" << mesh.C()[faceI] << nl
                 << "owner[" << faceI << "]为：" << mesh.owner()[faceI] << nl
                 << "neighbour[" << faceI << "]为：" << mesh.neighbour()[faceI] << nl;
        }
    }
    Info << endl;
 
    // 利用mesh.boundaryMesh()函数获取边界面信息
    // forAll是OpenFOAM定义的循环宏
    forAll(mesh.boundaryMesh(), patchI)
    {
        Info << "Patch[" << patchI << "]：" << mesh.boundary()[patchI].name() << "包括"
             << mesh.boundary()[patchI].Cf().size() << "个网格面" << nl
             << "起始面为：" << mesh.boundary()[patchI].start() << endl;
    }
    Info << endl;
 
    // 利用函数faceCells()得到网格面相邻的网格
    // 利用函数sf()可以得到网格面的面积向量，该向量的模即为面积
    //定义一个索引，仅为测试使用，如果要得到几何面积，则需要将该边界上的所有网格面积求和
    label patchFaceI(0);
    forAll(mesh.boundaryMesh(), patchI)
    {
        Info << "Patch[" << patchI << "]的网格面" << patchFaceI << "与网格["
             << mesh.boundary()[patchI].patch().faceCells()[patchFaceI] << "]相邻,"
             << "其面积向量：" << mesh.boundary()[patchI].Sf()[patchFaceI] << ",其面积为："
             << mag(mesh.boundary()[patchI].Sf()[patchFaceI])
             << endl;
    }
    Info << endl;
 
    // 查看网格细节
    // 将变量定义为const，则此变量只能读取，而不能修改
    const faceList &fcs = mesh.faces();            // 得到所有的网格面列表
    const List &pts = mesh.points();        //得到所有的网格节点列表
    const List ¢s = mesh.faceCentres(); //得到网格面心点列表
 
    // 遍历索引能被80整除的网格面
    forAll(fcs, faceI) if (faceI % 80 == 0)
    {
        if (faceI < mesh.Cf().size())
        {
            Info << "内部面";
        }
        else
        {
            forAll(mesh.boundary(), patchI) if ((mesh.boundary()[patchI].start() <= faceI) &&
                                                (faceI < mesh.boundary()[patchI].start() + mesh.boundary()[patchI].Cf().size()))
            {
                Info << "网格面在边界[" << patchI << "]上，网格面[";
                break;
            }
        }
 
        Info << faceI << "]的中心为：" << cents[faceI] << "，拥有" << fcs[faceI].size() << "个节点：";
 
        forAll(fcs[faceI], vertexI)
        {
            Info << " " << pts[fcs[faceI][vertexI]];
        }
        Info << endl;
    }
    Info << endl;
 
    // 判断边界类型是否为empty类型
    forAll(mesh.boundaryMesh(), patchID)
    {
        const polyPatch &pp = mesh.boundaryMesh()[patchID];
        if (isA(pp))
        {
            Info << "patch[" << patchID << "]:" << mesh.boundary()[patchID].name()
                 << "的边界类型为empty" << endl;
        }
        else
        {
            Info << "patch[" << patchID << "]:" << mesh.boundary()[patchID].name()
                 << "的边界类型不是empty" << endl;
        }
    }
    Info << endl;
 
    // 利用边界名称得到其ID
    word patchName("movingWall");
    label patchID = mesh.boundaryMesh().findPatchID(patchName);
    Info << "patch name:" << patchName << nl << "patchID:" << patchID << endl;
 
    Info << endl;
 
    return 0;
}

（本文完毕）

本篇文章来源于微信公众号: CFD之道

赞(3) 打赏
版权声明：未经允许，请勿随意用于商业用途。
文章名称：《OpenFOAM编程案例｜03 访问网格信息》
文章链接：https://www.topcfd.cn/18679/
本站资源仅供个人学习交流，请于下载后24小时内删除，不允许用于商业用途，否则法律问题自行承担。
分享到

相关推荐

说两句 抢沙发

评论前必须登录！

 

最新文章

2024 年 12 月
一 二 三 四 五 六 日
 12
3456789
10111213141516
17181920212223
24252627282930

觉得文章有用就打赏一下文章作者吧

非常感谢你的打赏，我们将继续给力更多优质内容，让我们一起创建更加美好的网络世界！

支付宝扫一扫

微信扫一扫

登录

找回密码

注册