OpenFOAMのフィールドオブジェクト参照方法

#はじめに

OpenFOAM/foam-extendのメッシュデータ構造

OpenFOAMでは、メッシュの各セル・面・点に与えられるデータをGeometricField<Type, PatchField, GeoMesh>クラスを使って管理しています。
OpenFOAM C++ Documentation: GeometricField

テンプレート引数
Type	scalar（double）、vector（三次元ベクトル）、symmTensor（二階対称テンソル）、 ...
PatchField	fvPatchField（セルに対応する境界条件）、fvsPatchField（面に対応する境界条件）
GeoMesh	volMesh（セルに対応するメッシュ）、surfaceMesh（面に対応するメッシュ）

OpenFOAMにおいてFieldクラスはListクラス（＝配列）に算術演算を定義したものです。
Fieldクラスに、入出力、データ登録、メッシュタイプ（セルまたは面）、メッシュパッチ（境界条件）を追加したものがGeometricFieldになります。

全ての組み合わせにtypedefで名前が付けられているので、

volVectorField U = ...;         // = GeometricField<vector, fvPatchField, volMesh>
surfaceVectorField gradp = ...; // = GeometricField<vector, fvsPatchField, surfaceMesh>

のように宣言することができます。

GeometricFieldの生成方法

例えばケースディレクトリ構造が以下のようになっている場合を考えます。

case/
  |----0/
  |    |----p
  |    |----U
  |
  |----const/
  |    |----fluidProperties
  |    |----polyMesh/
  |         |----boundary
  |         ... (mesh data)
  |
  |----system/
       |----fvSolution
       |----fvScheme
       |----controlDict

圧力pをvolScalarFieldとして宣言し、case/0/pを読み取って初期化する場合は次のようにします。

GeometricFieldの宣言方法(1)--ディレクトリに作成したファイルを読み取って初期化

volScalarField p
(
    IOobject
    (
        "p",                   // GeometricFieldの名前
        runTime.timeName(),    // 時間(time step)
        mesh,                  // メッシュオブジェクト
        IOobject::MUST_READ,   // 読み込みオプション --> 対応する時間のディレクトリから読み込む
        IOobject::AUTO_WRITE   // 書き込みオプション
    ),
    mesh                       // メッシュオブジェクト
);

他にも以下のような宣言＆初期化方法があります。

GeometricFieldの宣言方法(2)--GeometricFieldを与えて初期化

surfaceScalarField pf
(
    IOobject
    (
        "pf",                  // GeometricFieldの名前
        runTime.timeName(),    // 時間(time step)
        mesh,                  // メッシュオブジェクト
        IOobject::NO_READ,     // 読み込みオプション --> 対応する時間のディレクトリからは読み込まない
        IOobject::AUTO_WRITE   // 書き込みオプション
    ),
    fvc::interpolate(p, "p")   // 既に宣言したvolScalarField pを内挿したsurfaceScalarFieldを返し初期化
);

GeometricFieldの宣言方法(3)--ある値で初期化

volTensorField gradU
(
    IOobject
    (
        "grad(U)",             // GeometricFieldの名前
        runTime.timeName(),    // 時間(time step)
        mesh,                  // メッシュオブジェクト
        IOobject::NO_READ,     // 読み込みオプション
        IOobject::AUTO_WRITE   // 書き込みオプション
    ),
    mesh,                      // メッシュオブジェクト
    dimensionedTensor("zero", dimless, tensor::zero)   // ゼロテンソルで初期化
);

GeometricFieldの宣言方法(4)--tmpクラスオブジェクトを使って初期化

fluidModel fluid(mesh);              // 自作クラスfluidModel（流体に関する情報を取り扱うクラス）
volScalarField rho = fluid.rho();    // fluidModelのメンバ関数rho()、密度を返す。
// rho()はcase/constant/fluidPropertiesから密度を読み込みGeometricFieldを作成しtmpクラスを使って返す。

注意しなければならないことは、GeometricFieldは次元があるので、(3)の場合は必ずDimensionedScalar/DimensionedVector/...で初期化する必要がある点です。

GeometricFieldの参照方法

自作クラスを作っていると、ソルバーで宣言したGeometricFieldを参照して計算する必要がでます。OpenFOAMでは、objectRegistryクラスが生成されたGeometricFieldを管理しており、objectRegistryを経由して全てのGeometricFieldを参照することが可能です。
OpenFOAM C++ Documentation: objectRegistry
GeometricFieldは必ず名づけられて登録されているので、その名前をキーにしてGeometricFieldを検索するメンバ関数が複数用意されています。

fvMeshオブジェクトを経由する場合

例えば圧力フィールドpをあるクラスのメンバー関数内で参照したい場合は、次のようにfvMeshオブジェクトからobjectRegistryにアクセスし参照することができます。
（ほとんど全てのクラスはfvMeshオブジェクトを引数に受け取って生成されるため、プライベート変数としてfvMeshオブジェクトのアドレスを保持しており、このようなアクセスが可能となっています。）
OpenFOAM C++ Documentation: fvMesh

GeometricFieldのconst参照

const volScalarField& p = mesh.lookupObject<volScalarField>("p");

ただしこれはconst参照なのでp自体を変更することはできません。pを変更したい場合はconst_castする必要があります。

const_castによるGeometricFieldの非const参照

volScalarField& p =
    const_cast<volScalarField&>
    (
        mesh.lookupObject<volScalarField>("p")
    );

const_castを多用するのはあまりよくないと思うので、何か他の手を考えるべきでしょう。

fvPatchオブジェクトを経由する場合

fvPatchFieldを自作した場合、同じpatchに属する他の変数のfvPatchFieldを計算に使いたいことがあります。例えば、pフィールドのpatch Aを操作しているときにUフィールドのpatch Aにアクセスするようなときです。この場合、fvPatchのメンバ関数lookupPatchFieldを使用してアクセスすることができます。
OpenFOAM C++ Documentaion: fvPatch

例えばFoam::pressureInletOutletVelocity::updateCoeffs()では次のように使われています。
OpenFOAM-2.0.x: pressureInletOutletVelocityFvPatchVectorField.C

fvPatchFieldのconst参照

const fvsPatchField<scalar>& phip =
    patch().lookupPatchField<surfaceScalarField, scalar>(phiName_);

TIPS

1. Mesh boundaryにアクセスするにはmesh.boundaryMesh()
2. Field boundaryにアクセスするにはfield.boundaryField()