![見出し画像](https://assets.st-note.com/production/uploads/images/92600257/rectangle_large_type_2_cbf1583d29260f91c8e2761d07ccea52.jpeg?width=1200)
【OpenFOAM用のC++勉強(3)】メッシュ境界の情報取得
こんにちは。
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://assets.st-note.com/img/1670151422356-rDHiBp08wz.png?width=1200)
出力結果は、
============= 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は以下のように定義しなおしているため
![](https://assets.st-note.com/img/1670151970312-qhoNhMrrMt.png?width=1200)
GeometricFieldを調べれば良いですね。
![](https://assets.st-note.com/img/1670152031507-yZRRK2oyq8.png?width=1200)
結果は以下のようになります。
============= 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を継承して作れているので名前からリストであること予測できます。
![](https://assets.st-note.com/img/1670152802002-xvKtNRb1aV.png?width=1200)
ですので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などのも継承しているクラスです。
![](https://assets.st-note.com/img/1670326707832-99of1AddnX.png?width=1200)
その中のpolyMesh.Hを除くとpolyMeshの中にboundaryMesh()関数があり、これが境界情報を扱う関数となっています。
![](https://assets.st-note.com/img/1670326708046-3aoSWjiPDc.png?width=1200)
boundaryMesh()関数の結果はpolyBoundaryMeshクラスを返しているので、polyBoudaryMeshのメンバ関数を調べると取り扱いがわかります。
![](https://assets.st-note.com/img/1670326708339-GmpRPuDB9u.png?width=1200)
各境界面の情報を取得
各境界面には以下のように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ループです。
![](https://assets.st-note.com/img/1671457362172-WMt9qiqeN1.png?width=1200)
#define forAll(list, i) \
for (Foam::label i=0; i<(list).size(); ++i)
と定義しています。
参考にしているOpenFOAMの辞書です。
Foundation版をベースにしていますがわからないことはこちらで調べています。
引き続きC++の勉強を続けます。
Twitter➡@t_kun_kamakiri
ブログ➡宇宙に入ったカマキリ(物理ブログ)
この記事が気に入ったらサポートをしてみませんか?