【OpenFOAM用のC++勉強(3)】メッシュ境界の情報取得

2022年12月4日 20:37

こんにちは。
OpenFOAM用のC++の勉強のメモを自分のために残しています。

OpenFOAM v2112

前回の記事で書いたコードの続きに書いていきます。
ファイル構成は以下のようになっています。

.
├── Make
│   ├── files
├── meshOF.C
└── mycavity
    ├── 0
    │   ├── U
    │   └── p
    ├── constant
    │   └── transportProperties
    └── system
        ├── PDRblockMeshDict
        ├── blockMeshDict
        ├── controlDict
        ├── fvSchemes
        └── fvSolution

mycavityは後ほどチュートリアルからコピーしてきます。

meshOF.C

#include "argList.H"
#include "fvMesh.H"
#include "volFields.H"
#include "surfaceFields.H"

using namespace Foam;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:

int main(int argc, char *argv[])
{
    Info << "argc = " << argc << endl;
    Info << "argv = " << argv << endl;

    argList args(argc, argv);
    Info << "commandLine = " << args.commandLine() << endl;
    Info << "rootPath = " << args.rootPath() << endl;
    Info << "globalCaseName = " << args.globalCaseName() << endl;
    Info << "path = " << args.path() << endl;
    Info << "globalPath = " << args.globalPath() << endl;

    Info << "================== Create mesh ==================" << endl;
    #include "createTime.H"
    
	Info << "default region: " << fvMesh::defaultRegion << endl;
	Info << "mesh sub dir: " << fvMesh::meshSubDir << endl;

    fvMesh mesh
    (
        IOobject
        (
            fvMesh::defaultRegion,
            runTime.timeName(),
            runTime,
            IOobject::MUST_READ
        )
    );

    Info << "============= cell volumes ==============="<< endl;
    const DimensionedField< scalar, volMesh >  V = mesh.V();
    Info <<  mesh.V() << endl;

    Info << "============= Face area vector ===============" << endl;
    surfaceVectorField Sf = mesh.Sf();
    Info << Sf << endl;
    
    Info << "============= Face area magnitude ==============="<< endl;
    surfaceScalarField  magSf = mesh.magSf();
    Info << magSf << endl;
    
    Info << "============= Cell centres ===============" << endl;
    volVectorField C = mesh.C();
    Info << C << endl;
    
    Info << "============= Face centres ===============" << endl;
    surfaceVectorField Cf = mesh.Cf();
    Info << Cf << endl;

    Info << "============= mesh boundary name ===============" << endl;
    // Return reference to name.
    Info << mesh.name() << endl; //const word&


    Info << "============= mesh size ===============" << endl;
    // mesh.C() is volVectorField
    // https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1GeometricField.html
    Info << "mesh.C().size() = " << mesh.C().size() << endl;

    Info << "============= mesh.C().internalField() ===============" << endl;
    Info << mesh.C().internalField() << endl;

    Info << "============= mesh boundary patch name ===============" << endl;
    // 	Return reference to boundary mesh.
    forAll(mesh.boundary(),meshI)
    {
        Info << mesh.boundary()[meshI].name() << endl; //const fvBoundaryMesh&
    }
    
    Info << "============= mesh.C().boundaryField() ===============" << endl;
    //Info << mesh.C().boundaryField() << endl;
    forAll(mesh.C().boundaryField(),meshI)
    {
        Info << " === start === " << endl;
        Info << mesh.boundary()[meshI].name() << endl; 
        Info << mesh.C().boundaryField()[meshI] << endl; 
        Info << " === end === " << endl;
    }



    //Info << "============= Face fluxes ===============" << endl;
    //surfaceScalarField phi = mesh.phi();
    //Info << phi << endl;


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

    return 0;
}


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

追加したのは以下の部分です。

    Info << "============= mesh boundary name ===============" << endl;
    // Return reference to name.
    Info << mesh.name() << endl; //const word&

    Info << "============= mesh boundary patch name ===============" << endl;
    // 	Return reference to boundary mesh.
    forAll(mesh.boundary(),meshI)
    {
        Info << mesh.boundary()[meshI].name() << endl; //const fvBoundaryMesh&
    }

    Info << "============= mesh size ===============" << endl;
    // mesh.C() is volVectorField
    // https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1GeometricField.html
    Info << "mesh.C().size() = " << mesh.C().size() << endl;

    Info << "============= mesh.C().internalField() ===============" << endl;
    Info << mesh.C().internalField() << endl;

    Info << "============= mesh.C().boundaryField() ===============" << endl;
    //Info << mesh.C().boundaryField() << endl;
    forAll(mesh.C().boundaryField(),meshI)
    {
        Info << " === start === " << endl;
        Info << mesh.boundary()[meshI].name() << endl; 
        Info << mesh.C().boundaryField()[meshI] << endl; 
        Info << " === end === " << endl;
    }

Make/files

meshOF.C

EXE = $(FOAM_USER_APPBIN)/mesh003

Make/options

EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude

EXE_LIBS = \
    -lfiniteVolume

wmakeをすると/home/kamakiri/OpenFOAM/kamakiri-v2012/platforms/linux64Gcc63DPInt32Opt/bin/にmesh003という実行ファイルができます。
wmakeでのエラー内容が多すぎる場合は、ターミナルで確認するのが困難なので以下のようにしてlogファイルに掃き出すと良いでしょう。

 wmake 2>&1 | tee log

前回もやりましたがメッシュ情報の取得にはconstant内の以下の情報が必要です。

constant/polyMesh/boundary
constant/polyMesh/faces
constant/polyMesh/neighbour
constant/polyMesh/owner
constant/polyMesh/points

cp -r /opt/OpenFOAM/OpenFOAM-v2012/tutorials/incompressible/icoFoam/cavity/cavity mycavity

cd mycavity

でフォルダ移動します。
blockMeshを実行してメッシュを生成を行います。

blockMesh

そしてmesh003を実行します。
以下のようにしてlog.mesh003に出力結果を掃き出しておくと良いでしょう。

mesh003 > log.mesh003

メッシュの名前

1つ目はfvMeshクラスにはname()関数があるのでそれを使ってみました。
wordという文字列で返ってきます。

https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1fvMesh.html#a51ae053ef677d6c2d35ce01b46089706

出力結果は、

============= mesh boundary name ===============
region0

となります。
なんかよくわかりませんが、全体の領域の名前の事でしょうか。

メッシュのサイズ

メッシュのサイズとはcellの数のことで以下のようにします。

    Info << "============= mesh size ===============" << endl;
    // mesh.C() is volVectorField
    // https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1GeometricField.html
    Info << "mesh.C().size() = " << mesh.C().size() << endl;

mesh.C()とすることでセルの情報を取り出すことができます。
mesh.C()は前回も見ましたがvolVectorField型で返ってくるので、vectorが持っているsize()関数を使って要素数を出力しています。
結果は以下のようになります。

============= mesh size ===============
mesh.C().size() = 400

メッシュのセルの場を取得

    Info << "============= mesh.C().internalField() ===============" << endl;
    Info << mesh.C().internalField() << endl;

mesh.C()がvolVectorField型なので、volVectorField型が持つ関数を調べることで、どういった関数が使えるのかを知ることができます。
ちなみにvolVectorFieldは以下のように定義しなおしているため

GeometricFieldを調べれば良いですね。

https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1GeometricField.html

結果は以下のようになります。

============= mesh.C().internalField() ===============
dimensions      [0 1 0 0 0 0 0];


internalField   nonuniform List<vector> 
400
(
(0.0025 0.0025 0.005)
(0.0075 0.0025 0.005)
(0.0125 0.0025 0.005)
(0.0175 0.0025 0.005)
(0.0225 0.0025 0.005)

境界面の情報の取得

次に境界面の情報を取得します。

    Info << "============= mesh boundary patch name ===============" << endl;
    // 	Return reference to boundary mesh.
    forAll(mesh.boundary(),meshI)
    {
        Info << mesh.boundary()[meshI].name() << endl; //const fvBoundaryMesh&
    }

mesh.boundary()はFoam::fvBoundaryMeshを返し、それがfvPatchListを継承して作れているので名前からリストであること予測できます。

ですのでforAllでループさせながらmesh.boundary()[meshI].name()とすると境界面の名前を取得できます。
結果は以下のようになります。

============= mesh boundary patch name ===============
movingWall
fixedWalls
frontAndBack

境界面の名前を取得できました。

境界面の座標を取得

境界面を取得できたので境界面のセル内の座標を取得します

    Info << "============= mesh.C().boundaryField() ===============" << endl;
    //Info << mesh.C().boundaryField() << endl;
    forAll(mesh.C().boundaryField(),meshI)
    {
        Info << " === start === " << endl;
        Info << mesh.boundary()[meshI].name() << endl; 
        Info << mesh.C().boundaryField()[meshI] << endl; 
        Info << " === end === " << endl;
    }

mesh.C()がセルの座標でvolVectorFieldを返します。
volVectorFieldはGeometricFieldを定義しなおしたセル用のベクトルなので、実態はGeometricFieldとなります。
GeometricField内の関数boundaryField()を使うと境界の情報を取得することができるということです。

結果は以下のようになります。

============= mesh.C().boundaryField() ===============
 === start === 
movingWall
type            sliced;
value           nonuniform List<vector> 
20
(
(0.0025 0.1 0.005)
(0.0075 0.1 0.005)
(0.0125 0.1 0.005)
(0.0175 0.1 0.005)

（省略）

 === end === 
 === start === 
fixedWalls
type            sliced;
value           nonuniform List<vector> 
60
(
(0 0.0025 0.005)
(0 0.0075 0.005)
(0 0.0125 0.005)

（省略）

 === end === 
 === start === 
frontAndBack
type            sliced;
value           nonuniform List<vector> 0();

 === end ===

メッシュ境界の情報の取得

次は境界面に対して情報を取得してみます。ファイル構成は以下のようになっています。

.
├── Make
│   ├── files
│   └── options
├── log
├── meshOF.C
└── mycavity
├── 0
│   ├── U
│   ├── banana
│   ├── banana_org
│   └── p
├── constant
│   ├── polyMesh
│   │   ├── boundary
│   │   ├── faces
│   │   ├── neighbour
│   │   ├── owner
│   │   └── points
│   └── transportProperties
└── system
├── PDRblockMeshDict
├── blockMeshDict
├── controlDict
├── fvSchemes
└── fvSolution

mainOF.C

#include "argList.H"
#include "fvMesh.H"
#include "volFields.H"
#include "surfaceFields.H"

using namespace Foam;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:

int main(int argc, char *argv[])
{
    Info << "argc = " << argc << endl;
    Info << "argv = " << argv << endl;

    argList args(argc, argv);
    Info << "commandLine = " << args.commandLine() << endl;
    Info << "rootPath = " << args.rootPath() << endl;
    Info << "globalCaseName = " << args.globalCaseName() << endl;
    Info << "path = " << args.path() << endl;
    Info << "globalPath = " << args.globalPath() << endl;

    Info << "================== Create mesh ==================" << endl;
    #include "createTime.H"
    
	Info << "default region: " << fvMesh::defaultRegion << endl;
	Info << "mesh sub dir: " << fvMesh::meshSubDir << endl;

    fvMesh mesh
    (
        IOobject
        (
            fvMesh::defaultRegion,
            runTime.timeName(),
            runTime,
            IOobject::MUST_READ
        )
    );

    Info << "============= cell volumes ==============="<< endl;
    const DimensionedField< scalar, volMesh >  V = mesh.V();
    Info <<  mesh.V() << endl;

    Info << "============= Face area vector ===============" << endl;
    surfaceVectorField Sf = mesh.Sf();
    Info << Sf << endl;
    
    Info << "============= Face area magnitude ==============="<< endl;
    surfaceScalarField  magSf = mesh.magSf();
    Info << magSf << endl;
    
    Info << "============= Cell centres ===============" << endl;
    volVectorField C = mesh.C();
    Info << C << endl;
    
    Info << "============= Face centres ===============" << endl;
    surfaceVectorField Cf = mesh.Cf();
    Info << Cf << endl;

    Info << "============= mesh boundary name ===============" << endl;
    // Return reference to name.
    Info << mesh.name() << endl; //const word&


    Info << "============= mesh size ===============" << endl;
    // mesh.C() is volVectorField
    // https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1GeometricField.html
    Info << "mesh.C().size() = " << mesh.C().size() << endl;

    Info << "============= mesh.C().internalField() ===============" << endl;
    Info << mesh.C().internalField() << endl;

    Info << "============= mesh boundary patch name ===============" << endl;
    // 	Return reference to boundary mesh.
    forAll(mesh.boundary(),meshI)
    {
        Info << mesh.boundary()[meshI].name() << endl; //const fvBoundaryMesh&
    }
    
    Info << "============= mesh.C().boundaryField() ===============" << endl;
    //Info << mesh.C().boundaryField() << endl;
    forAll(mesh.C().boundaryField(),meshI)
    {
        Info << " === start === " << endl;
        Info << mesh.boundary()[meshI].name() << endl; 
        Info << mesh.C().boundaryField()[meshI] << endl; 
        Info << " === end === " << endl;
    }

    Info << "============= mesh.boundaryMesh()===============" << endl;
    // https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1polyMesh.html
    // fvMesh class inherited from polyMesh
    // boundaryMesh() is a member function of polyMesh
    // const polyBoundaryMesh & 	boundaryMesh () const
    Info << mesh.boundaryMesh() << endl;
    Info << "boundary names = " << mesh.boundaryMesh().names() << endl;
    Info << "boundary types = " << mesh.boundaryMesh().types() << endl;
    Info << "boundary patchSizes = " << mesh.boundaryMesh().patchSizes() << endl;   
    Info << "boundary physicalTypes = " << mesh.boundaryMesh().physicalTypes() << endl;

    Info << "============= mesh.boundaryMesh()[patchI] ===============" << endl;
    forAll(mesh.boundaryMesh(), patchI)
    {
        Info << "boundary name = " << mesh.boundaryMesh()[patchI].name() << endl;
        Info << "boundary type = " << mesh.boundaryMesh()[patchI].typeName << endl;
    }

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

    return 0;
}

追加したのは以下の部分です。

    Info << "============= mesh.boundaryMesh()===============" << endl;
    // https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1polyMesh.html
    // fvMesh class inherited from polyMesh
    // boundaryMesh() is a member function of polyMesh
    // const polyBoundaryMesh & 	boundaryMesh () const
    Info << mesh.boundaryMesh() << endl;
    Info << "boundary names = " << mesh.boundaryMesh().names() << endl;
    Info << "boundary types = " << mesh.boundaryMesh().types() << endl;
    Info << "boundary patchSizes = " << mesh.boundaryMesh().patchSizes() << endl;   
    Info << "boundary physicalTypes = " << mesh.boundaryMesh().physicalTypes() << endl;

    Info << "============= mesh.boundaryMesh()[patchI] ===============" << endl;
    forAll(mesh.boundaryMesh(), patchI)
    {
        Info << "boundary name = " << mesh.boundaryMesh()[patchI].name() << endl;
        Info << "boundary type = " << mesh.boundaryMesh()[patchI].typeName << endl;
    }

Make/files

meshOF.C

EXE = $(FOAM_USER_APPBIN)/mesh004

Make/options

EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude

EXE_LIBS = \
    -lfiniteVolume

以下のコマンドでコンパイル。

 wmake 2>&1 | tee log

以下のコマンドでmycavityのケースファイルに移動します。

cd mycavity

meshは上でfvMeshクラスをインスタンス化した際に付けた名前です。
※チュートリアル

cp -r /opt/OpenFOAM/OpenFOAM-v2012/tutorials/incompressible/icoFoam/cavity/cavity mycavity

そしてblockMeshを実行しておきます。

blockMesh

これはメッシュ情報には以下のファイルが必要だからです。

constant/polyMesh/boundary
constant/polyMesh/faces
constant/polyMesh/neighbour
constant/polyMesh/owner
constant/polyMesh/points

では、mycavityの中で以下のmesh004を実行してます。
実行結果は長いのでlog.mesh004に掃き出します。

mesh004 > log.mesh004

log.mesh004

argc = 1
argv = 1
/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : _7bdb509494-20201222 OPENFOAM=2012
Arch   : "LSB;label=32;scalar=64"
Exec   : mesh004
Date   : Dec 06 2022
Time   : 20:11:33
Host   : DESKTOP-CT961AV
PID    : 1679
I/O    : uncollated
Case   : /mnt/d/openfoam/20221011_wolf_benchmark/OpenFOAM_Training/OF9/101programming/myWork/mesh004/mycavity
nProcs : 1
trapFpe: Floating point exception trapping enabled (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 5, maxFileModificationPolls 20)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
commandLine = "mesh004"
rootPath = "/mnt/d/openfoam/20221011_wolf_benchmark/OpenFOAM_Training/OF9/101programming/myWork/mesh004"
globalCaseName = "mycavity"
path = "/mnt/d/openfoam/20221011_wolf_benchmark/OpenFOAM_Training/OF9/101programming/myWork/mesh004/mycavity"
globalPath = "/mnt/d/openfoam/20221011_wolf_benchmark/OpenFOAM_Training/OF9/101programming/myWork/mesh004/mycavity"
================== Create mesh ==================
Create time

default region: region0
mesh sub dir: polyMesh
============= cell volumes ===============
dimensions      [0 3 0 0 0 0 0];


value           nonuniform List<scalar> 
400
(
2.5e-07
2.5e-07

（省略）
============= Face area vector ===============
dimensions      [0 2 0 0 0 0 0];
oriented        oriented;


internalField   nonuniform List<vector> 
760
(
(5e-05 0 0)
(0 5e-05 0)
（省略）
============= Face area magnitude ===============
dimensions      [0 2 0 0 0 0 0];


internalField   nonuniform List<scalar> 
760
(
5e-05
（省略）
============= Cell centres ===============
dimensions      [0 1 0 0 0 0 0];


internalField   nonuniform List<vector> 
400
(
(0.0025 0.0025 0.005)

（省略）
============= Face centres ===============
dimensions      [0 1 0 0 0 0 0];


internalField   nonuniform List<vector> 
760
(
(0.005 0.0025 0.005)

（省略）


============= mesh boundary name ===============
region0
============= mesh size ===============
mesh.C().size() = 400
============= mesh.C().internalField() ===============
dimensions      [0 1 0 0 0 0 0];


internalField   nonuniform List<vector> 
400
(
(0.0025 0.0025 0.005)


（省略）
============= mesh boundary patch name ===============
movingWall
fixedWalls
frontAndBack
============= mesh.C().boundaryField() ===============
 === start === 
movingWall
type            sliced;
value           nonuniform List<vector> 
20
(
(0.0025 0.1 0.005)
(0.0075 0.1 0.005)

（省略）


============= mesh.boundaryMesh()===============
3
(
    movingWall
    {
        type            wall;
        inGroups        1(wall);
        nFaces          20;
        startFace       760;
    }
    fixedWalls
    {
        type            wall;
        inGroups        1(wall);
        nFaces          60;
        startFace       780;
    }
    frontAndBack
    {
        type            empty;
        inGroups        1(empty);
        nFaces          800;
        startFace       840;
    }
)
boundary names = 3(movingWall fixedWalls frontAndBack)
boundary types = 3(wall wall empty)
boundary patchSizes = 3(20 60 800)
boundary physicalTypes = 3(  )
============= mesh.boundaryMesh()[patchI] ===============
boundary name = movingWall
boundary type = patch
boundary name = fixedWalls
boundary type = patch
boundary name = frontAndBack
boundary type = patch

以下の部分が今回追加した内容の出力結果です。

============= mesh.boundaryMesh()===============
3
(
    movingWall
    {
        type            wall;
        inGroups        1(wall);
        nFaces          20;
        startFace       760;
    }
    fixedWalls
    {
        type            wall;
        inGroups        1(wall);
        nFaces          60;
        startFace       780;
    }
    frontAndBack
    {
        type            empty;
        inGroups        1(empty);
        nFaces          800;
        startFace       840;
    }
)
boundary names = 3(movingWall fixedWalls frontAndBack)
boundary types = 3(wall wall empty)
boundary patchSizes = 3(20 60 800)
boundary physicalTypes = 3(  )
============= mesh.boundaryMesh()[patchI] ===============
boundary name = movingWall
boundary type = patch
boundary name = fixedWalls
boundary type = patch
boundary name = frontAndBack
boundary type = patch

解説をします。

fvMeshはメッシュ情報の他に、方程式のスキームを扱うfvSchemsや、マトリクス解法を指定するfvSolutionなどのも継承しているクラスです。

その中のpolyMesh.Hを除くとpolyMeshの中にboundaryMesh()関数があり、これが境界情報を扱う関数となっています。

boundaryMesh()関数の結果はpolyBoundaryMeshクラスを返しているので、polyBoudaryMeshのメンバ関数を調べると取り扱いがわかります。

各境界面の情報を取得

各境界面には以下のようにforAllでループを回してひとつひとつ取り出しています。

    Info << "============= mesh.boundaryMesh()[patchI] ===============" << endl;
    forAll(mesh.boundaryMesh(), patchI)
    {
        Info << "boundary name = " << mesh.boundaryMesh()[patchI].name() << endl;
        Info << "boundary type = " << mesh.boundaryMesh()[patchI].typeName << endl;
    }

出力結果は以下の部分です。

============= mesh.boundaryMesh()[patchI] ===============
boundary name = movingWall
boundary type = patch
boundary name = fixedWalls
boundary type = patch
boundary name = frontAndBack
boundary type = patch

ちなみにOpenFOAMのコードの中で出てくるforAllはC++の文法というわけではなく、OpenFOAMが独自で定義したforループです。

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

と定義しています。

参考にしているOpenFOAMの辞書です。
Foundation版をベースにしていますがわからないことはこちらで調べています。

引き続きC++の勉強を続けます。

Twitter➡@t_kun_kamakiri
ブログ➡宇宙に入ったカマキリ（物理ブログ）