东莞市弘盛再生资源回收有限公司
栏目分类
热点资讯
新闻资讯

你的位置:东莞市弘盛再生资源回收有限公司 > 新闻资讯 > 图片除了frontAndBack的网格面数外

图片除了frontAndBack的网格面数外

发布日期:2024-06-25 10:51    点击次数:128

图片除了frontAndBack的网格面数外

公司简介

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

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

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

图片

1 准备文件

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

cd $FOAM_RUNmkdir demo3cd demo3touch demo3.Cmkdir Makecd Maketouch filestouch optionscd ..

文件结构如下图所示。

图片

修改files文件的内容
demo3.CEXE = 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 cavityblockMesh

执行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()函数可以获取边界面上的信息。

添加以下代码。

荔蒲县垂经干果有限公司 0,
福州爱泰进出口有限公司 0.55) 0px 2px 10px;">// 利用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<point> &pts = mesh.points();        //得到所有的网格节点列表const List<point> &cents = 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<emptyPolyPatch>(pp))    {        Info << "patch[" << patchID <<"]:" << mesh.boundary()[patchID].name()            << "的边界类型为empty" << endl;    }    else    {        Info << "patch[" << patchID << "]:"<< mesh.boundary()[patchID].name()            <<"的边界类型不是empty" << endl;    }}

运行结果如下图所示。

图片

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

// 利用边界名称得到其IDword 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<point> &pts = mesh.points();        //得到所有的网格节点列表    const List<point> &cents = 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<emptyPolyPatch>(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;} 

(本文完毕)

本站仅提供存储服务,所有内容均由用户发布,如发现有害或侵权内容,请点击举报。

友情链接:

Powered by 东莞市弘盛再生资源回收有限公司 @2013-2022 RSS地图 HTML地图

Copyright 365建站 © 2013-2024 SSWL 版权所有