セットの作成

2024年9月8日

はじめに

セルやフェイスの集まりでセットやゾーンを作る方法について。

使用バージョン

OpenFOAM v2406

セットの作成

セットを作成する方法として、setSet を使う方法と topoSet を使う方法がある。簡単に作るには setSet を使えばよい。

setSet は対話形式で操作する。setSet を実行するとコマンド入力待ちになる。

$ setSet
...
readline>

時間ディレクトリがない場合は、"-constant" オプションをつけて実行する。

基本コマンド

• help : 使い方の表示
• list : セットの表示
• quit : 終了

セルセットはつぎのコマンドで作成する。

cellSet <setName> <action>

フェイスセットの場合は faceSet を使う。setName は任意である。

action には以下のものがある。

• list : 表示
• clear : クリア
• invert : セットの選択をさかさまにする (選択されているもの以外が選択される)
• remove : セットの削除
• new <source> : セットの新規作成
• add <source> : 選択の追加
• delete <source> : 選択の削除
• subset <source> : セットから選択する (セットは上書きされる)

source にはたとえば以下のようなものがある。

cellSet

• boxToCell (x1 y1 z1) (x2 y2 z2) : セルをボックスで選択
• cylinderToCell (x1 y1 z1) (x2 y2 z2) radius : セルを円筒で選択 (座標は中心軸のベクトルの端点)
• cellToCell <setName> : セルセットを選択
• faceToCell <setName> <neighbour|owner|any|all> : フェイスセットからセルを選択
• zoneToCell <zoneName>: セルゾーンを選択

faceSet

• boxToFace (x1 y1 z1) (x2 y2 z2) : フェイスをボックスで選択
• faceToFace <setName> : フェイスセットを選択
• patchToFace <patchName> : パッチを選択
• zoneToFace <zoneName>: フェイスゾーンを選択

source のリストがほしい場合は、source に適当な文字列を入れてみればよい。

たとえば、(0, 0, 0) と (0.1, 0.01, 0.01) の 2 点が作る線を対角線とする六面体内にその中心が含まれるセル群からセルセットを作る場合、つぎのようにする。

readline>cellSet c0 new boxToCell (0 0 0) (0.1 0.01 0.01)

セルセット c0 にさらにセルの選択を追加したい場合は、つぎのようにする。

readline>cellSet c0 add boxToCell (0 0.01 0) (0.1 0.02 0.01)

セルセット c0 内からセルを選択して別のセルセットを作る場合は、つぎのようにする。

readline>cellSet c1 new cellToCell c0
readline>cellSet c1 subset boxToCell (0.025 0.005 0) (0.075 0.015 0.01)

フェイスセットの作成も同様である。

readline>faceSet f0 new patchToFace movingWall

ゾーンの作成

セルセットからセルゾーンを作る場合、つぎのようにする。

readline>cellZoneSet c0zone new setToCellZone c0

フェイスの場合は faceZoneSet/setToFaceZone を使う。

バッチ形式での実行

コマンドをバッチ形式で実行することもできる。たとえば "setBatch" という名前のファイルに setSet のコマンドを書き込んでおき、つぎのように実行する。

$ setSet -batch setBatch

辞書形式での設定

topoSet を使えば、同じことを辞書による設定で行うことができる。

たとえば、つぎのように設定する。

system/topoSetDict

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      topoSetDict;
}

actions
(
    {
        name    c0;
        type    cellSet;
        action  new;
        source  boxToCell;
        sourceInfo
        {
            box (0 0 0) (0.1 0.01 0.01);
        }
    }
    {
        name    c0zone;
        type    cellZoneSet;
        action  new;
        source  setToCellZone;
        sourceInfo
        {
            set c0;
        }
    }
);

実行

$ topoSet

source によって設定項目が異なる。

cylinderToCell

    {
        name    c0;
        type    cellSet;
        action  new;
        source  cylinderToCell;
        sourceInfo
        {
            p1 (0 0 0);
            p2 (0 0 0.6);
            radius 0.1;
        }
    }

patchToFace

    {
        name    f0;
        type    faceSet;
        action  new;
        source  patchToFace;
        sourceInfo
        {
            name   freeSurface;
        }
    }

setToFaceZone

    {
        name    f0;
        type    faceZoneSet;
        action  new;
        source  setToFaceZone;
        sourceInfo
        {
            faceSet    f0;
        }
    }

rotatedBoxToCell

    {
        name    c0;
        type    cellSet;
        action  new;
        source  rotatedBoxToCell;
        sourceInfo
        {
            origin (0.05 0.0217 0.0);
            i (0.0283 0.0283 0.0);
            j (-0.0283 0.0283 0.0);
            k (0.0 0.0 0.01);
        }
    }

rotatedBoxToCell は回転した箱による選択が可能だが、原点および箱を作る 3 ベクトルを指定するという設定になっている。ベクトルは自分で計算しなくてはならない。Python で計算すると次のようになる。

rotated_box.py

#!/bin/env python

import numpy as np

def rotation_matrix(theta):
    c = np.cos(theta)
    s = np.sin(theta)
    return np.array([
	[c, -s, 0],
	[s, c, 0],
	[0, 0, 1]
    ])

center = np.array([0.05, 0.05, 0])
box_min = np.array([0.03, 0.03, 0])
box_max = np.array([0.07, 0.07, 0.01])

theta = np.radians(45)
R = rotation_matrix(theta)

v1 = np.array([box_max[0] - box_min[0], 0, 0])
v2 = np.array([0, box_max[1] - box_min[1], 0])
v3 = np.array([0, 0, box_max[2] - box_min[2]])

origin = np.dot(R, box_min - center) + center
i = np.dot(R, v1)
j = np.dot(R, v2)
k = np.dot(R, v3)

output_file = "system/rotatedBox"
with open(output_file, "w") as f:
    f.write("rotatedBox\n")
    f.write("{\n")
    f.write(f"    origin ({origin[0]} {origin[1]} {origin[2]});\n")
    f.write(f"    i ({i[0]} {i[1]} {i[2]});\n")
    f.write(f"    j ({j[0]} {j[1]} {j[2]});\n")
    f.write(f"    k ({k[0]} {k[1]} {k[2]});\n")
    f.write("}\n")

このコードは system/rotatedBox に出力するようになっており、topoSetDict で取り込むことができる。

#include "rotatedBox"

actions
(
    {
        name    c0;
        type    cellSet;
        action  new;
        source  rotatedBoxToCell;
        sourceInfo
        {
            $rotatedBox;
        }
    }
);

labelToCell

#include "selectedCells"

actions
(
    {
        name    selectedCellsSet;
        type    cellSet;
        action  new;
        source  labelToCell;
        sourceInfo
        {
            value $selectedCells;
        }
    }
);

system/selectedCells

selectedCells
60
(
109
110
128
...
);

labelToCell でセルを直接指定できる。どうやってセルを知るかであるが、function object の writeCellCentres を使えばよい。

system/controlDict

functions
{
    #includeFunc writeCellCentres
}

実行

$ postProcess -time 0

セル中心データが出てくる。

$ ls 0
C  Cx  Cy  Cz  U  p

これをプログラムで処理すればよい。Python だと次のようになる。

#!/bin/env python

import numpy as np
import os

def get_cell_centers(mesh_dir):
    cell_centers_path = os.path.join(mesh_dir, "C")
    cell_centers = []

    with open(cell_centers_path, "r") as f:
        lines = f.readlines()
        for line in lines:
            if line.startswith("(") and line.endswith(")\n"):
                coords = line[1:-2].split()
                cell_centers.append([float(coord) for coord in coords])

    return np.array(cell_centers)

def rotation_matrix(theta):
    c = np.cos(theta)
    s = np.sin(theta)
    return np.array([
	[c, -s, 0],
	[s, c, 0],
	[0, 0, 1]
    ])

def is_inside_rotated_box(points, box_min, box_max, center, R):
    inv_R = np.linalg.inv(R)
    rotated_points = np.dot(points - center, inv_R.T) + center
    inside = np.all((rotated_points >= box_min) & (rotated_points <= box_max), axis=1)
    return inside

center = np.array([0.05, 0.05, 0])
box_min = np.array([0.03, 0.03, 0])
box_max = np.array([0.07, 0.07, 0.01])

theta = np.radians(45)
R = rotation_matrix(theta)

mesh_dir = "0"
cell_centers = get_cell_centers(mesh_dir)

inside = is_inside_rotated_box(cell_centers, box_min, box_max, center, R)
selected_cells = np.where(inside)[0]

output_file = "system/selectedCells"
with open(output_file, "w") as f:
    f.write(f"selectedCells\n")
    f.write(f"{len(selected_cells)}\n")
    f.write("(\n")
    for cell_id in selected_cells:
        f.write(f"{cell_id}\n")
    f.write(");\n")

この例は、上の rotated_box.py と同じ選択をする例になっている。