はじめに
下のようなメッシュ細分化を行った際の内容をまとめました。
意地でもヘキサ（ポリヘドラ）メッシュを扱いたいという人向けの内容になっております。 実際はここまでこだわる必要のない場合が殆どだと思いますのでドツボにはまらないようにご注意ください。
refineMesh
OpenFOAMではrefineMeshというユーティリティを使用することで歪みや形状などを保持したままヘキサ型＋ポリヘドラ型（２Dだと四角形＋五角形）にメッシュを細分化することができます。
この設定はsystem/refineMeshDictで以下のように行います。
// 細分化するcellSet set refineSet; // 細分化方法 coordinateSystem global; //coordinateSystem patchLocal; //coordinateSystem fieldBased; // 細分化の方向の定義 globalCoeffs { tan1 (1 0 0); tan2 (0 1 0); } // 定義した方向のうち使用するもの directions ( tan1 tan2 ); // ヘキサメッシュを保持するかどうか useHexTopology true; // 形状を変更してもよいか（useHexTopologyとの併用は不可） geometricCut false; // 細分化過程を出力するかどうか writeMesh false;
この中の細分化方法ではglobalを設定しています。 これは細分化の方向を指定する方法であり上記のような全体が四角形の場合の分割は容易に行うことができます。
円メッシュ細分化時の問題点
しかし例えば全体が円形のメッシュを分割する場合には以下のようにうまく分割できない場合が発生します。
このような場合にはglobalではなくfieldBasedを使用する必要がありますが、これがなかなか曲者でしたのでやりかたをまとめてみました。
※なおヘキサメッシュを3方向に細分化する場合であれば角度によらずほぼうまくいきます。この問題は1方向及び2方向で細分化する場合に発生するようです。
テストケース
今回は例として以下のような細分化を行う手順を紹介していきます。
topoSetによる下準備
分割する箇所をcellCetとして指定する必要があるためtopoSetによりrefineSetという名前を付けます。 また内側（四角形部分）と外側（円周部分）をinnerSet、outerSetとして区別しておきます。
system/topoSetDict
height 0.3; actions ( { name innerSet; type cellSet; action new; source zoneToCell; zone innerZone; } { name innerSet; type cellSet; action delete; source boxToCell; box (-10 -10 $height) (10 10 10); } { name outerSet; type cellSet; action new; source zoneToCell; zone outerZone; } { name outerSet; type cellSet; action delete; source boxToCell; box (-10 -10 $height) (10 10 10); } { name refineSet; type cellSet; action new; source boxToCell; box (-10 -10 -10) (10 10 $height); } );
fieldBasedで読み込むファイル
fieldBasedをどのように設定するかについてはrefineMesh.CのあるディレクトリのrefineMeshDictに以下のように書かれています。
$FOAM_UTILITIES/mesh/manipulation/refineMesh/refineMeshDict
// List of directions to refine, if "fieldBased". Keep in mind that these // fields must be of type "vectorField", not "volVectorField". //directions //( // radialDirectionFieldName // angularDirectionFieldName // heightDirectionFieldName //);
つまりvectorFieldの形式で各セルの分割方向を記述したファイルを用意すればOKとのことです。 volVectorFieldではなくvectorFieldというところがポイントであり両者の違いは後ろの方で説明します。
ファイルの作成
任意のvectorFieldを出力するために以下のことを行いました。
- vectorFieldを作成し保存するユーティリティを作成する
- そのユーティリティを実行する際にcodedFunctionObjectで各セルの細分化方向を指定する
一見ハードルが高そうに見えますがOpenFOAMのチュートリアルにあるものの改造です。
ユーティリティの作成
calcRadiusFieldを参考に以下のようなmakeDirectionというユーティリティを作成しました。
詳しくは後述します。
codedFunctionObjectで各セルの細分化方向を指定する
実際に細分化方向を規定するのはここです。 controlDictに以下のようにfunctionObjectを記述します。
system/controlDict
// 上略 functions { directions { libs ("libutilityFunctionObjects.so"); type coded; name directions; codeInclude #{ #include "topoSet.H" #}; codeEnd #{ // 細分化方法を指定するコード #}; } }
細分化を指定するコードについては後述します。
ユーティリティの実行
上記のユーティリティを実行すると以下のように各セルにおける細分化方向が記述されたファイルが『0』ディレクトリに出力されます。
makeDirection -overwrite
refineMeshDictの編集
最後にrefineMeshDictを以下のように編集します。
set refineSet; // 細分化するcellSet coordinateSystem fieldBased; // 作成した細分化方向 directions ( radialDirection angularDirection // heightDirection // 今回は高さ方向は無し ); useHexTopology true; geometricCut false; writeMesh false;
refineMeshの実行
この状態でrefineMeshを実行することでメッシュを細分化することができます。
refineMesh -overwrite
最後に
ヘキサメッシュでないとモヤモヤするような悪癖のある方は是非試してみてください。
補足
ポリヘドラの表示方法
ポリヘドラメッシュをParaViewで表示する際には以下のVTK Polyhedraにチェックを入れます。
vectorFieldとvolVectorFieldの違い
volVectorFieldはvectorFieldに次元の情報と境界の情報が追加されている形になります。
例）セルが３つのメッシュの場合
volVectorField
FoamFile
{
version 2.0;
format ascii;
class volVectorField; // ファイルのタイプ
location "0";
object angularVolDirection;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 0 0 0 0]; // 次元
internalField nonuniform List<vector> // internalFieldの中身がvectorField
3 // セルの数
(
(0 0 1) // 0セル目のベクトル
(0 0 1) // 1セル目のベクトル
(0 0 1) // 2セル目のベクトル
);
boundaryField // 境界の情報
{
wall
{
type calculated;
value uniform (0 0 0);
}
}
vectorField
// ヘッダ FoamFile { version 2.0; format ascii; class vectorField; // ファイルのタイプ location "0"; object angularDirection; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // 3 // セルの数 ( (0 0 1) // 0セル目のベクトル (0 0 1) // 1セル目のベクトル (0 0 1) // 2セル目のベクトル );
そのためvolVectorFieldを作成して、classの項目を変更して次元と境界情報とinternalField nonuniform List
ユーティリティの説明
vectorIOField radialDirection
(
IOobject
(
"radialDirection",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh.V().size()
);
のようにして『0』ディレクトリに『radialDirection』という名前のファイルを出力します。
ParaViewで表示するために同じ内容のvolVectorFieldも作成します。
volVectorField radialVolDirection
(
IOobject
(
"radialVolDirection",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedVector("r", dimless, Zero)
);
フィールドを作成した後に
runTime.functionObjects().end();
とすることで、controlDict内のcodeEndに書かれたコードを実行できます。
codedFunctionの中身
細分化コードでは以下のようにまず各vectorFieldを読み込みます。
vectorField& radialDirection = mesh().lookupObjectRef<vectorField>("radialDirection"); vectorField& angularDirection = mesh().lookupObjectRef<vectorField>("angularDirection"); vectorField& heightDirection = mesh().lookupObjectRef<vectorField>("heightDirection");
また場所を指定するために以下のようにcellSetを読み込みます。
const word setType = "cellSet"; PtrList<topoSet> sets(2); sets.set ( 0, topoSet::New ( setType, mesh(), "innerSet", IOobject::MUST_READ, IOobject::NO_WRITE ).ptr() );
そして以下のようにすることでcellSet内でセルループを行うことができます。
forAllConstIter(labelHashSet, sets[0], iter) { const label& i = iter.key(); // i = セル番号 // 各セルについてベクトルを指定 }
このループの中では以下のように指定していきます。
innerSet
const label& i = iter.key(); radialDirection[i] = vector( 1.0, 1.0, 0.0); // 北東方向 angularDirection[i] = vector(-1.0, 1.0, 0.0); // 北西方向 heightDirection[i] = vector( 0.0, 0.0, 1.0); // 高さ方向
outerSet：
const label& i = iter.key(); const vector& C = mesh().C()[i]; const scalar& x = C[0]; const scalar& y = C[1]; const scalar& z = C[2]; radialDirection[i] = vector( x, y, 0.0); // 半径方向 angularDirection[i] = vector( -y, x, 0.0); // 円周方向 heightDirection[i] = vector(0.0, 0.0, 1.0); // 高さ方向
