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

#はじめに

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

OpenFOAMでは、メッシュの各セル・面・点に与えられるデータをGeometricField\<Type, PatchField, GeoMesh>クラスを使って管理しています。
[OpenFOAM C++ Documentation: GeometricField](http://www.openfoam.com/documentation/cpp-guide/html/a01001.html)


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

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

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

```cpp
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を読み取って初期化する場合は次のようにします。

```cpp:GeometricFieldの宣言方法(1)--ディレクトリに作成したファイルを読み取って初期化
volScalarField p
(
    IOobject
    (
        "p",                   // GeometricFieldの名前
        runTime.timeName(),    // 時間(time step)
        mesh,                  // メッシュオブジェクト
        IOobject::MUST_READ,   // 読み込みオプション --> 対応する時間のディレクトリから読み込む
        IOobject::AUTO_WRITE   // 書き込みオプション
    ),
    mesh                       // メッシュオブジェクト
);
```

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

```cpp: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を返し初期化
);
```

```cpp: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)   // ゼロテンソルで初期化
);
```

```cpp: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](http://www.openfoam.com/documentation/cpp-guide/html/a01822.html#affcae20935bdf1726ec1f822228a2330)
GeometricFieldは必ず名づけられて登録されているので、その名前をキーにしてGeometricFieldを検索するメンバ関数が複数用意されています。


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

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

```cpp:GeometricFieldのconst参照
const volScalarField& p = mesh.lookupObject<volScalarField>("p");
```

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

```cpp: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](http://foam.sourceforge.net/docs/cpp/a00876.html)

例えばFoam::pressureInletOutletVelocity::updateCoeffs()では次のように使われています。
[OpenFOAM-2.0.x: pressureInletOutletVelocityFvPatchVectorField.C](https://github.com/OpenFOAM/OpenFOAM-2.0.x/blob/master/src/finiteVolume/fields/fvPatchFields/derived/pressureInletOutletVelocity/pressureInletOutletVelocityFvPatchVectorField.C)

```cpp:fvPatchFieldのconst参照
const fvsPatchField<scalar>& phip =
    patch().lookupPatchField<surfaceScalarField, scalar>(phiName_);
```

# TIPS

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