【OpenFOAM用のC++勉強(5)】境界面における速度場と圧力場の取得
こんにちは。
OpenFOAM用のC++の勉強のメモを自分のために残しています。
OpenFOAM v2112
前回の復習（指定した座標での速度場の取得）
前回は指定した座標での速度場や圧力場の値を取得しました。
１．ベクトルとして座標を定義します。
vector U_pos(0.05, 0.05, 0);２．座標からセルNo.を取得します。
abel cellUNo = mesh.findCell(U_pos);３．セルNo.から値の取得
label cellUNo = mesh.findCell(U_pos);
Info << "U" << U_pos << " = " << U[cellUNo] << endl;
Info << "U" << U_pos << " = " << U.ref()[cellUNo] << endl;結果は以下のようになります。
U(0.05 0.05 0) = (-0.196068 0.0263898 0)
U(0.05 0.05 0) = (-0.196068 0.0263898 0)境界面の情報を取得
全体のコードを示します。
#include "argList.H"
#include "fvMesh.H"
#include "volFields.H"
#include "surfaceFields.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:
int main(int argc, char *argv[])
{
Info << "================== Create mesh ==================" << endl;
argList args(argc, argv);
#include "createTime.H"
// mesh
fvMesh mesh
(
IOobject
(
fvMesh::defaultRegion,
runTime.timeName(),
runTime,
IOobject::MUST_READ
)
);
Info << "runTime.timeName() = " << runTime.timeName() << endl;
Info << "============= volScalarField p ===============" << endl;
// volScalarField typedef=> GeometricField => DimensionedField => Field
// https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1DimensionedField.html
word t_ = "0.5";
volScalarField p
(
IOobject
(
"p",
// runTime.timeName(), // word
t_,
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info << t_ << endl;
vector P_pos(0.05, 0.05, 0);
Info << "get p value at P_pos = " << P_pos << endl;
// mesh is fvMesh class.
// fvMesh => polyMesh
/*
Foam::label Foam::polyMesh::findCell
(
const point& p,
const cellDecomposition decompMode
)
*/
label cellPNo = mesh.findCell(P_pos);
Info << "p" << P_pos << " = " << p[cellPNo] << endl;
Info << "p" << P_pos << " = " << p.ref()[cellPNo] << endl; //https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1GeometricField.html#a77a3ea1ce7e2adc04d27301292b095ae
Info << "============= volVectorField U ===============" << endl;
volVectorField U
(
IOobject
(
"U",
// runTime.timeName(), // word
t_,
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
vector U_pos(0.05, 0.05, 0);
Info << "get U value at U_pos = " << U_pos << endl;
label cellUNo = mesh.findCell(U_pos);
Info << "U" << U_pos << " = " << U[cellUNo] << endl;
Info << "U" << U_pos << " = " << U.ref()[cellUNo] << endl;
Info << "============= Boundary p value ===============" << endl;
Info << "mesh.boundaryMesh()" << endl;
Info << mesh.boundaryMesh() << endl;
Info << "============================" << endl;
forAll(mesh.boundaryMesh(), patchI)
{
Info << "mesh.boundaryMesh()[" << patchI << "] = \n" <<
mesh.boundaryMesh()[patchI] << endl;
}
// https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1fvBoundaryMesh.html#ac89fe0aab9fe12bd4f0d3807495fe72b
label inletPressID = mesh.boundaryMesh().findPatchID("movingWall");
Info << "inletvelID = " << inletPressID << endl;
const fvPatchScalarField& PinletList = p.boundaryField()[inletPressID];
forAll(PinletList, fi)
{
Info << "p[" << fi << "] =" << PinletList[fi] << endl;
}
Info << "============= Boundary U value ===============" << endl;
label inletvelID = mesh.boundaryMesh().findPatchID("movingWall");
Info << "inletvelID = " << inletvelID << endl;
const fvPatchVectorField& UinletList = U.boundaryField()[inletvelID];
forAll(UinletList, fi)
{
Info << "U[" << fi << "] =" << UinletList[fi] << endl;
}
Info << "============= p of cell value near boundary ===============" << endl;
// fvPatchScalarField = fvPatchField<scalar>
// fvpatch have path()
// https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1fvPatch.html#acb76222e9692fae6c93d7a49125229f7
const labelList& PfaceCell =PinletList.patch().faceCells();
forAll(PfaceCell, ci)
{
Info << "cellNo = " << PfaceCell[ci] << " : pFaceCell[" << ci << "] =" << p[PfaceCell[ci]] << endl;
}
Info << "============= U of cell value near boundary ===============" << endl;
const labelList& UfaceCell =UinletList.patch().faceCells();
forAll(UfaceCell, ci)
{
Info << "cellNo = " << UfaceCell[ci] << " : UFaceCell[" << ci << "] =" << U[UfaceCell[ci]] << endl;
}
return 0;
}Make/files
meshOF.C
EXE = $(FOAM_USER_APPBIN)/mesh005Make/options
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lfiniteVolume以下のコマンドでコンパイル。
wmake 2>&1 | tee log※チュートリアルは以下からコピーしてくる。
cp -r /opt/OpenFOAM/OpenFOAM-v2012/tutorials/incompressible/icoFoam/cavity/cavity mycavity以下のコマンドでmycavityのケースファイルに移動します。
cd mycavityそしてblockMeshを実行しておきます。
blockMeshこれはメッシュ情報には以下のファイルが必要だからです。
constant/polyMesh/boundary
constant/polyMesh/faces
constant/polyMesh/neighbour
constant/polyMesh/owner
constant/polyMesh/points
そして流速Uや圧力pの情報を取得するには計算を実行しておく必要があります。
icoFoamでは、mycavityの中で以下のmesh005を実行してます。
実行結果は長いのでlog.mesh005に掃き出します。
mesh005 > log.meesh005log.mesh007
================== Create mesh ==================
/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / 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 : mesh007
Date : Dec 11 2022
Time : 21:47:31
Host : DESKTOP-CT961AV
PID : 4762
I/O : uncollated
Case : /mnt/d/openfoam/20221011_wolf_benchmark/OpenFOAM_Training/OF9/101programming/myWork/mesh007/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
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time
runTime.timeName() = 0
============= volScalarField p ===============
0.5
get p value at P_pos = (0.05 0.05 0)
p(0.05 0.05 0) = -0.00133453
p(0.05 0.05 0) = -0.00133453
============= volVectorField U ===============
get U value at U_pos = (0.05 0.05 0)
U(0.05 0.05 0) = (-0.196068 0.0263898 0)
U(0.05 0.05 0) = (-0.196068 0.0263898 0)
============= Boundary p value ===============
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;
}
)
============================
mesh.boundaryMesh()[0] =
type wall;
inGroups 1(wall);
nFaces 20;
startFace 760;
mesh.boundaryMesh()[1] =
type wall;
inGroups 1(wall);
nFaces 60;
startFace 780;
mesh.boundaryMesh()[2] =
type empty;
inGroups 1(empty);
nFaces 800;
startFace 840;
inletPressID = 0
p.boundaryField()[0][5] = -0.240143
p[0] = -4.36666
p[1] = -2.04463
p[2] = -0.973538
p[3] = -0.563177
p[4] = -0.357842
p[5] = -0.240143
p[6] = -0.170115
p[7] = -0.123779
p[8] = -0.087477
p[9] = -0.0535059
p[10] = -0.0166499
p[11] = 0.0276653
p[12] = 0.0846122
p[13] = 0.161734
p[14] = 0.272416
p[15] = 0.442103
p[16] = 0.718223
p[17] = 1.22571
p[18] = 2.41104
p[19] = 4.84854
============= Boundary U value ===============
inletvelID = 0
U[0] = (1 0 0)
U[1] = (1 0 0)
U[2] = (1 0 0)
U[3] = (1 0 0)
U[4] = (1 0 0)
U[5] = (1 0 0)
U[6] = (1 0 0)
U[7] = (1 0 0)
U[8] = (1 0 0)
U[9] = (1 0 0)
U[10] = (1 0 0)
U[11] = (1 0 0)
U[12] = (1 0 0)
U[13] = (1 0 0)
U[14] = (1 0 0)
U[15] = (1 0 0)
U[16] = (1 0 0)
U[17] = (1 0 0)
U[18] = (1 0 0)
U[19] = (1 0 0)
============= p of cell value near boundary ===============
PinletPatch = movingWall
cellNo = 380 : pFaceCell[0] = -4.36666
cellNo = 381 : pFaceCell[1] = -2.04463
cellNo = 382 : pFaceCell[2] = -0.973538
cellNo = 383 : pFaceCell[3] = -0.563177
cellNo = 384 : pFaceCell[4] = -0.357842
cellNo = 385 : pFaceCell[5] = -0.240143
cellNo = 386 : pFaceCell[6] = -0.170115
cellNo = 387 : pFaceCell[7] = -0.123779
cellNo = 388 : pFaceCell[8] = -0.087477
cellNo = 389 : pFaceCell[9] = -0.0535059
cellNo = 390 : pFaceCell[10] = -0.0166499
cellNo = 391 : pFaceCell[11] = 0.0276653
cellNo = 392 : pFaceCell[12] = 0.0846122
cellNo = 393 : pFaceCell[13] = 0.161734
cellNo = 394 : pFaceCell[14] = 0.272416
cellNo = 395 : pFaceCell[15] = 0.442103
cellNo = 396 : pFaceCell[16] = 0.718223
cellNo = 397 : pFaceCell[17] = 1.22571
cellNo = 398 : pFaceCell[18] = 2.41104
cellNo = 399 : pFaceCell[19] = 4.84854
============= U of cell value near boundary ===============
cellNo = 380 : UFaceCell[0] = (0.299604 0.146091 0)
cellNo = 381 : UFaceCell[1] = (0.396272 0.127503 0)
cellNo = 382 : UFaceCell[2] = (0.55742 0.0760895 0)
cellNo = 383 : UFaceCell[3] = (0.674487 0.0437233 0)
cellNo = 384 : UFaceCell[4] = (0.743889 0.0255721 0)
cellNo = 385 : UFaceCell[5] = (0.787154 0.0152422 0)
cellNo = 386 : UFaceCell[6] = (0.815486 0.00921237 0)
cellNo = 387 : UFaceCell[7] = (0.833956 0.00550944 0)
cellNo = 388 : UFaceCell[8] = (0.845386 0.00302719 0)
cellNo = 389 : UFaceCell[9] = (0.851372 0.00114935 0)
cellNo = 390 : UFaceCell[10] = (0.852667 -0.000499394 0)
cellNo = 391 : UFaceCell[11] = (0.849383 -0.00221509 0)
cellNo = 392 : UFaceCell[12] = (0.841 -0.00433859 0)
cellNo = 393 : UFaceCell[13] = (0.826141 -0.00742411 0)
cellNo = 394 : UFaceCell[14] = (0.802078 -0.0125323 0)
cellNo = 395 : UFaceCell[15] = (0.763621 -0.0217729 0)
cellNo = 396 : UFaceCell[16] = (0.699108 -0.0392829 0)
cellNo = 397 : UFaceCell[17] = (0.583952 -0.0726601 0)
cellNo = 398 : UFaceCell[18] = (0.415827 -0.127868 0)
cellNo = 399 : UFaceCell[19] = (0.308819 -0.149468 0)以下コードの解説です。
Info << "============= Boundary p value ===============" << endl;からコードを追加しています。
mesh.boundaryMesh()はpolyMeshクラスの中のメンバ関数boundaryMeshで以下のように、polyBoundaryMeshを返します。
const polyBoundaryMesh & boundaryMesh()const
Info << "============= Boundary p value ===============" << endl;
Info << "mesh.boundaryMesh()" << endl;
Info << mesh.boundaryMesh() << endl;これにより結果は境界面の情報をリストとして取得することができます。
============= Boundary p value ===============
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;
}
)for文でリストからひとつひとつ要素を取得することもできます。
Info << "============= Boundary p value ===============" << endl;
forAll(mesh.boundaryMesh(), patchI)
{
Info << "mesh.boundaryMesh()[" << patchI << "] = \n" <<
mesh.boundaryMesh()[patchI] << endl;
}境界面の名前から境界面の速度と圧力の値を取得
１．境界面の名前からパッチのIDを取得
fvMeshはpolyMeshを継承しており、polyMeshのメンバ関数であるboundaryMesh()はpolyBoundaryMeshを返しています。polyBoundaryMeshはfindPatchIDというメンバ関数があり、これが境界面パッチのIDを取得する関数です。
label inletPressID = mesh.boundaryMesh().findPatchID("movingWall");
Info << "inletPressID = " << inletPressID << endl;inletPressID = 0とパッチ面のIDを取得できました。
２．pはvolScalarField型（GeometricFieldのスカラー版がvolScalarField型）なので、GeometricFieldのメンバ関数boundaryFIeldを使うことで速度や圧力の値を取得できます。
const fvPatchScalarField& PinletList = p.boundaryField()[inletPressID];
forAll(PinletList, fi)
{
Info << "p[" << fi << "] =" << PinletList[fi] << endl;
}結果は以下のようになります。
inletPressID = 0
p[0] =-4.36666
p[1] =-2.04463
p[2] =-0.973538
p[3] =-0.563177
p[4] =-0.357842ちなみにp.boundaryField()だけだと以下のように境界面の定義のリストが出力されます。
p.boundaryField()
3
(
type zeroGradient;
type zeroGradient;
type empty;
)以下のようにp.boundaryField()[境界面ID]から続けて[5]のようにindexを入れると値の取得ができます。
Info << "p.boundaryField()[0][5] = " << p.boundaryField()[0][5] << endl;結果は以下のようになります。
p.boundaryField()[0][5] = -0.240143同様に速度についても以下のように境界面での値を取得できます。
Info << "============= Boundary U value ===============" << endl;
label inletvelID = mesh.boundaryMesh().findPatchID("movingWall");
Info << "inletvelID = " << inletvelID << endl;
const fvPatchVectorField& UinletList = U.boundaryField()[inletvelID];
forAll(UinletList, fi)
{
Info << "U[" << fi << "] =" << UinletList[fi] << endl;
}境界面近傍のセルの値での速度と圧力の取得
fvPatchVectorFieldやfvPatchVectorFieldは以下のように、
fvPatchScalarField = fvPatchField<scalar>
fvPatchVectorField = fvPatchField<vector>
と定義されているので実態はfvPatch型です。
Info << "============= p of cell value near boundary ===============" << endl;
// fvPatchScalarField = fvPatchField<scalar>
// fvpatch have path()
// https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1fvPatch.html#acb76222e9692fae6c93d7a49125229f7
const fvPatch& PinletPatch = PinletList.patch();
Info << "PinletPatch = " << PinletPatch.name() << endl;
const labelList& PfaceCell =PinletList.patch().faceCells();
forAll(PfaceCell, ci)
{
Info << "cellNo = " << PfaceCell[ci] << " : pFaceCell[" << ci << "] = " << p[PfaceCell[ci]] << endl;
}
Info << "============= U of cell value near boundary ===============" << endl;
const labelList& UfaceCell =UinletList.patch().faceCells();
forAll(UfaceCell, ci)
{
Info << "cellNo = " << UfaceCell[ci] << " : UFaceCell[" << ci << "] = " << U[UfaceCell[ci]] << endl;
}const fvPatch& PinletPatch = PinletList.patch()とするとfvPatch型を返すのでメンバ関数findCellsにより値のリストを返します。
返ってきたのがリストなのでfor文でひとつひとつ要素を取り出します。
forAll(PfaceCell, ci)
{
Info << "cellNo = " << PfaceCell[ci] << " : pFaceCell[" << ci << "] = " << p[PfaceCell[ci]] << endl;
}
Info << "============= U of cell value near boundary ===============" << endl;
const labelList& UfaceCell =UinletList.patch().faceCells();
forAll(UfaceCell, ci)
{
Info << "cellNo = " << UfaceCell[ci] << " : UFaceCell[" << ci << "] = " << U[UfaceCell[ci]] << endl;
}結果は以下のようになります。
============= p of cell value near boundary ===============
PinletPatch = movingWall
cellNo = 380 : pFaceCell[0] = -4.36666
cellNo = 381 : pFaceCell[1] = -2.04463
cellNo = 382 : pFaceCell[2] = -0.973538
cellNo = 383 : pFaceCell[3] = -0.563177
cellNo = 384 : pFaceCell[4] = -0.357842
cellNo = 385 : pFaceCell[5] = -0.240143
cellNo = 386 : pFaceCell[6] = -0.170115
cellNo = 387 : pFaceCell[7] = -0.123779
cellNo = 388 : pFaceCell[8] = -0.087477
cellNo = 389 : pFaceCell[9] = -0.0535059
cellNo = 390 : pFaceCell[10] = -0.0166499
cellNo = 391 : pFaceCell[11] = 0.0276653
cellNo = 392 : pFaceCell[12] = 0.0846122
cellNo = 393 : pFaceCell[13] = 0.161734
cellNo = 394 : pFaceCell[14] = 0.272416
cellNo = 395 : pFaceCell[15] = 0.442103
cellNo = 396 : pFaceCell[16] = 0.718223
cellNo = 397 : pFaceCell[17] = 1.22571
cellNo = 398 : pFaceCell[18] = 2.41104
cellNo = 399 : pFaceCell[19] = 4.84854
============= U of cell value near boundary ===============
cellNo = 380 : UFaceCell[0] = (0.299604 0.146091 0)
cellNo = 381 : UFaceCell[1] = (0.396272 0.127503 0)
cellNo = 382 : UFaceCell[2] = (0.55742 0.0760895 0)
cellNo = 383 : UFaceCell[3] = (0.674487 0.0437233 0)
cellNo = 384 : UFaceCell[4] = (0.743889 0.0255721 0)
cellNo = 385 : UFaceCell[5] = (0.787154 0.0152422 0)
cellNo = 386 : UFaceCell[6] = (0.815486 0.00921237 0)
cellNo = 387 : UFaceCell[7] = (0.833956 0.00550944 0)
cellNo = 388 : UFaceCell[8] = (0.845386 0.00302719 0)
cellNo = 389 : UFaceCell[9] = (0.851372 0.00114935 0)
cellNo = 390 : UFaceCell[10] = (0.852667 -0.000499394 0)
cellNo = 391 : UFaceCell[11] = (0.849383 -0.00221509 0)
cellNo = 392 : UFaceCell[12] = (0.841 -0.00433859 0)
cellNo = 393 : UFaceCell[13] = (0.826141 -0.00742411 0)
cellNo = 394 : UFaceCell[14] = (0.802078 -0.0125323 0)
cellNo = 395 : UFaceCell[15] = (0.763621 -0.0217729 0)
cellNo = 396 : UFaceCell[16] = (0.699108 -0.0392829 0)
cellNo = 397 : UFaceCell[17] = (0.583952 -0.0726601 0)
cellNo = 398 : UFaceCell[18] = (0.415827 -0.127868 0)
cellNo = 399 : UFaceCell[19] = (0.308819 -0.149468 0)ちなみに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
ブログ➡宇宙に入ったカマキリ（物理ブログ）
この記事が気に入ったらサポートをしてみませんか？