Capytaine: a Python-based distribution of Nemoh
https://ancell.in/capytaine/latest/

訳 今井康貴（imaiy@cc.saga-u.ac.jp）

1. インストール

CapytaineはWindows、MacOS [1]、Linuxで利用可能です。

[1] arm64アーキテクチャ(Apple M1)の最新情報については、https://github.com/capytaine/capytaine/issues/190 を参照してください。

CapytaineはPython 3.6以降を必要とします。したがって、現在サポートされているすべてのバージョンのPythonと互換性があります。

1.1 Condaを用いる場合

Capytaineをインストールする最も簡単な方法は、Condaで利用できるプリコンパイルパッケージを用いることです。
1) Anacondaディストリビューション、またはその軽量版であるMinicondaとMiniforgeをダウンロードし、インストールしてください。
2) Condaがインストールされたら、専用の環境を作りたいと思うかもしれません。Capytaineのパッケージはconda-forgeチャンネルで公開されており、以下のコマンドでインストールできます。
conda install -c conda-forge capytaine
必要な依存関係が自動的にインストールさます。

オプションの依存関係

オプションの依存関係は、手動でインストールします。これらはあると便利ですが、Capytaine の主な機能には必要ありません。

名前 インストールコマンドの例 使用方法
matplotlib conda install matplotlib ドキュメントとクックブックのいくつかの例で使用されています。
vtk conda install -c conda-forge vtk 3Dビジュアライゼーションに使用します。
joblib conda install joblib 並列処理に使用します。
quadpy pip install quadpy 高次求積用(実験的)


1.2 Pipを用いたインストール

このパッケージはPyPIで入手可能ですが、ソース配布のみです。つまり、パッケージをインストールするためには、Fortranコンパイラ[2]が必要です。もしそうなら、Capytaineを以下のようにインストールできます。

pip install capytaine

もし、コンパイラをインストールできない場合は、代わりにCondaを使うことをお勧めします。

[2] 例えば、UbuntuやDebianの場合: sudo apt install gfortran. macOSの場合、例えばこの手順を参照してください。
https://github.com/capytaine/capytaine/issues/115#issuecomment-1143987636

2. チュートリアル

2.1 主なコンセプト

メッシュ
平均的な位置にある浮体のメッシュ。Meshクラスのインスタンスとして格納される。
メッシュは、頂点のリスト（頂点は実数値座標の三つ組）と面のリスト（面は頂点のインデックスの四つ組）として定義されます。デフォルトでは、面は四角形であると仮定される。三角形の面は、2つの同じ頂点を持つ四角形としてサポートされます。
Mesh クラスは、頂点と面から計算されたデータ（面の法線、面の中心、面の面積など）も保存します。

DOF
	自由度（Dof）は、浮体の平均位置を中心とした小さな動きや変形を定義する。メッシュの各面の中心にベクトルとして格納されます。
radiating_dof はボディの実際の動きを表し、influcted_dof は（一般化された）力の成分を表します。
注
Galerkin境界要素法の分野の数学者にとっては、自由度の概念は異なる意味を持つかもしれない（Galerkin分解の基底関数）。ここでは、自由度は浮体の物理的な自由度であり、典型的には剛体の並進と回転である。

浮体（FloatingBody）
FloatingBodyは、主にMeshといくつかの自由度の再集合である。
ボディの自由度は名前（例：Heave）で参照される。これらは定義された順番に並んでいるはずですが、コード上は厳密に保証されていません。インデックスでアクセスするよりも、名前でアクセスする方が好ましいでしょう。
メッシュとDOF以外に、質量や重心位置などの物理情報をFloatingBodyインスタンスに格納できる。この情報は、いくつかの特定の動作にのみ必要であり（静水圧を参照）、多くの場合、指定しないままにしておくことができます。

LinearPotentialFlowProblem（リニアポテンシャルフロー問題）
問題は、いくつかのパラメタの集まりです：FloatingBody、波の角周波数，水深, 水密度, 重力加速度.
クラス LinearPotentialFlowProblem は，2つの子クラスを持つ．
RadiationProblem (放射している物体名も必要)
DiffractionProblem (入射波面の角度も必要)
ほとんどのパラメタはオプションです。それらが提供されない場合、デフォルト値が使用されます（問題の設定ページを参照してください）。

ソルバ
コードの核となる部分です。LinearPotentialFlowProblemを入力として受け取り、LinearPotentialFlowResultを返すsolve()メソッドを持っています。
グリーン関数を計算するクラスと行列を構築するクラスを呼び出します。詳しくは，「Resolution」を参照してください．

LinearPotentialFlowResult（リニアポテンシャルフローリザルト）
結果を格納するクラスは，問題を格納するクラスと似ていますが，
radiation問題ではresult.added_massesやresult.radiation_dampings，
diffraction問題ではresult.forceなどの補足データがあります．
力は、自由度の名前と値を関連付けた辞書として格納されています。
値は、この自由度に沿った力の積分値です。
例えば、直交座標の剛体にかかる力ベクトルの成分を取り出すには、Surge、Sway、Heaveに関する力の値を確認します。

2.2 ステップバイステップの例

ipythonのような対話型Pythonコンソールを起動します。Capytaineの主要な機能は全て読み込むことができます。
from capytaine import *.
CapytaineはPythonのloggingモジュールを使っていることに注意してください。次に、オプションとして、以下のコマンドでloggingモジュールを初期化することで、コードからフィードバックを得ることができます。
import logging
logging.basicConfig(level=logging.INFO)
INFOをDEBUGに置き換えて、ソルバ内部で起こっているすべての情報を得ることができます。一方、レベルをWARNINGに設定すると、重要な警告だけがソルバから出力されます（これはloggingモジュールが設定されていない場合のデフォルトの動作です）。

メッシュのロード
このチュートリアルでは、Capytaineに含まれるメッシュジェネレータの1つを使い、簡単な幾何学的形状を作成します。
sphere = Sphere(radius=1.0, center=(0, 0, -2), name="my buoy")
また、ドキュメントの「メッシュと浮体」のセクションにあるように、ユーザーは様々なファイルフォーマットからメッシュをインポートできる。メッシュは、Meshオブジェクトとして保存されます。例えば、次のような構文で頂点の座標、面の中心、面の法線ベクトルにアクセスできる。
sphere.mesh.vertices[:10] # 最初の10個の頂点。
sphere.mesh.faces_centers[5] # 6番目の面の中心（Pythonの配列は0から始まります）.
sphere.mesh.faces_normals[5] # 6番目の面の法線ベクトル．
vtkがインストールされている場合、以下のようにしてメッシュを3D表示できる。
sphere.show()

DOFの定義
diffractionやradiationの問題を解く前に、物体の自由度（dofs）を定義する必要があります。これは、いくつかの方法で行うことができます。
手動で行う方法: ベクトルのリストを定義し、各ベクトルが面の中心でのボディの変位となるようにします。以下の例は、x方向の剛体運動の最も単純な例です。

sphere.dofs['Surge'] = [(1, 0, 0) for face in sphere.mesh.faces].

剛体の並進と回転を定義するためにヘルパー関数が用意されています。例えば
z方向の動きには、FloatingBody.add_translation_dofを使うことができます。
これは、"Surge"、"Sway"、"Heave "などのいくつかのdof名を認識できる。
sphere.add_translation_dof(name="Heave")。

FloatingBody.add_rotation_dof と
FloatingBody.add_all_rigid_body_dofs のドキュメントを参照してください。

自由度は、dofs辞書に格納されています。あるボディの dofs の名前にアクセスするには、たとえば次のようにします。
print(sphere.dofs.keys())
# dict_keys(['Surge', 'Heave'])

静水力学
Capytaineは与えられたメッシュに対して静水圧の計算を直接行うことができます。与えられたFloatingBodyの体積、没水面積、水面面積、浮力中心、メタセンター半径、高さ、静水力学、内部質量などのパラメタを得ることができます。
コードにボディに関する情報を与えてみましょう。
sphere.center_of_mass = (0, 0, -2)
sphere.rotation_center = (0, 0, -2)
回転中心は回転DOFを定義するために使用される点である。
(静水圧法の現在の制限により、剛体のDOF（ここではサージとヒーブ）があるとすぐに、それが回転DOFでなくても回転中心の定義が必要です)。

各静水圧パラメタは専用のメソッドで計算できる:
print(sphere.volume)
# 3.82267415555807
print(sphere.center_of_buoyancy)
# [-3.58784373e-17 -2.59455034e-17 -2.00000000e+00]
print(sphere.compute_hydrostatic_stiffness())
# <xarray.DataArray 'hydrostatic_stiffness' (influenced_dof: 2, radiating_dof: 2)>
# array([[0.00000000e+00, 0.00000000e+00],
#[0.00000000e+00, 2.38246922e-13]])
# Coordinates:
# * influenced_dof (influenced_dof) <U5 'Surge' 'Heave'
# * radiating_dof (radiating_dof) <U5 'Surge' 'Heave'
print(sphere.compute_rigid_body_inertia())
# <xarray.DataArray 'inertia_matrix' (influenced_dof: 2, radiating_dof: 2)>
# array([[3822.67415556,0.],
#[ 0., 3822.67415556]])
# Coordinates:
# * influenced_dof (influenced_dof) <U5 'Surge' 'Heave'
# * radiating_dof (radiating_dof) <U5 'Surge' 'Heave'

ここでの行列は2x2行列です。球体に対して2つのDOFしか定義していないからです。
また、すべての静水圧パラメタを計算し、パラメタと値のdictを返すcompute_hydrostaticsメソッドを使用できる。
Hydrostatics = sphere.compute_hydrostatics()

注意
静水圧パラメタを計算する前に、body.keep_immersed_part()を使用してメッシュを切り出したい場合があります。球体は完全に浸漬されているので、ここでは必要ありません。クロッピングはcompute_hydrostatics()関数に含まれています。

線形ポテンシャル流問題の定義
球のヒーブに対する放射問題を定義してみましょう。
from numpy import infty
problem = RadiationProblem(body=sphere, radiating_dof="Heave", omega=1.0, sea_bottom=-infty, g=9.81, rho=1000)
引数radiating_dofはbody引数で与えられた浮体のdofのうちの一つの名前でなければならない。波の角周波数はω=1rad/s，水深は無限大、重力加速度は9.81m/s2、水の密度は1000kg/m3とした。
これらの最後のパラメタは実際にはオプションである。デフォルトの値を使っているので、次のように放射問題を定義できる。
problem = RadiationProblem(body=sphere, radiating_dof="Heave", omega=1.0) と定義することもできます。
さらにいくつかのパラメタは自動的に計算されます, たとえば:
print(problem.wavenumber)
# 0.1019367991845056
print(problem.period)
# 6.283185307179586

問題の解法
BEMソルバを初期化します。
solver = BEMSolver()
ソルバの設定はこの時点でも可能ですが、このチュートリアルではデフォルトの設定を使用します。それでは、先ほど定義した問題を解いてみましょう。
result = solver.solve(problem)
solve() メソッドは結果オブジェクトを返します。solve() メソッドは結果オブジェクトを返します。結果オブジェクトは問題から得たすべてのデータを含んでいます。
print(result.omega)
# 1.0
print(result.body.name)
# "my buoy"
print(result.radiating_dof)
# "Heave"
print(result.period)
# 6.283185307179586
もちろん、出力データも保存されています。放射の問題を解決したので、付加質量と放射減衰にアクセスできるようになりました。
print(result.added_masses)
# {'Surge': 9.154531598110083e-06, 'Heave': 2207.8423200090374}.
added_masses辞書には、ボディの「影響を受けたDOF」それぞれにかかる結果の力が格納されています。
この例では、放射DOFはヒーブであり、そのx方向反力(result.add_masses['Surge'])は，heave方向の反力(result.add_masses['Heave'])に対して無視できるほど小さい。
print(result.radiation_dampings)
# {'Surge': -5.792518686098536e-07, 'Heave': 13.62318484050783}


結果を配列にまとめる
サージに対する付加質量と放射減衰を計算する。
other_problem = RadiationProblem(body=sphere, radiating_dof="Surge", omega=1.0)
other_result = solver.solve(other_problem)。
この2回目の解法は1回目の解法よりも高速になることに注意してください。
ソルバはこの体の中間データを保存しており、この他の問題を解くためにそれらを再利用します。
結果は以下のようにまとめることができる:
dataset = assemble_dataset([result, other_result])
新しいオブジェクトはxarrayパッケージのNetCDFライクなデータセットである。
これは、結果オブジェクトから追加された質量と放射減衰を整理して保存しています。
この例では、基本的に2つの2x2マトリックスです。この行列は、例えば次のようにアクセスできる。
dataset['added_mass'].sel(radiating_dof=["Surge", "Heave"], influenced_dof=["Surge", "Heave"], omega=1.0)
おそらく、各テストを個別に定義することなく、広範囲のパラメタについて問題を解決したいことでしょう。これはソルバのfill_datasetメソッドで行うことができます。問題の設定を参照してください。

3. メッシュと浮体

3.1 Meshmagick によるメッシュのインポート
既存のメッシュファイルを使用して新しい船体を作成するには、次の構文を使用します。
from capytaine import FloatingBody
body=FloatingBody.from_file('path/to/mesh.dat', file_format='nemoh')
上記の例では、Nemoh のメッシュフォーマットを使用しています。
Meshmagick のおかげで、他の多くのメッシュ形式をインポートできる。ファイルフォーマットはオプションの file_format 引数で指定します。フォーマットが指定されていない場合は、ファイルの拡張子から推測します。
body=FloatingBody.from_file('path/to/mesh.msh') # gmshファイル
現在読み込みに対応している形式は、以下の表のとおりです (Meshmagick のドキュメントを参考にしました)。
ファイルの拡張子，ソフトウェア，キーワード，追加機能
.mar	NEMOH [1]	NEMOH, MAR	対称性
.nem	NEMOH [1]	nemoh_mesh, nem
.gdf	WAMIT [2]	wamit, gdf
.inp	DIODORE [3]	diodore-inp, inp
.DAT	DIODORE[3]	DIODORE-dat
.hst	HYDROSTAR [4]	HYDROSTAR,HST, 対称性
.nat	natural, nat
.msh	GMSH 2 [5]	gmsh, msh
.rad	RADIOSS	RAD,RADIOSS
.stl	stl
.vtu	PARAVIEW [6]	vtu
.vtp	PARAVIEW [6]	vtp
.vtk	PARAVIEW [6]	paraview-legacy, vtk
.tec	TECPLOT [7]	tecplot, tec
.med	SALOME [8]	med, salome


[1] (1,2)NEMOHは、Ecole Centrale de Nantes (LHEEA)で開発されたオープンソースの船舶挙動解析ソフトウエアです。
[2]WAMITは、WAMIT社で開発された船舶航行用BEMソフトウエアです。
[3] (1,2)DIODOREはPRINCIPIA社によって開発された船舶のBEMソフトウェアです。
[4] HYDROSTAR は、BUREAU VERITAS 社が開発した船舶航行用 BEM ソフトウエアです。
[5] GMSH は、C. Geuzaine と J.-F. Remacle によって開発されたオープンソースのメッシュ作成ソフトウェアです。Remacleによって開発されたオープンソースのメッシングソフトウェアです。ファイルフォーマットVer.4は現在未対応です。
[6] (1,2,3)PARAVIEWは、Kitware社によって開発されたオープンソースの可視化ソフトウェアです。
[7] TECPLOTは、Tecplot社が開発した可視化ソフトウェアです
[8] SALOME-MECA は、EDF-R&D が開発したオープンソースの計算力学ソフトウェアです。

メッシュファイルを読み込む際に、すべてのメタデータが考慮されるわけではありません。例えば、対称性は.marと.hstファイルフォーマットでのみ考慮されます。Githubでissueを作成し、改善点をご指摘ください。

3.2 Meshioでメッシュをインポートする

メッシュは meshio ライブラリを使用してインポートすることもできます。上記のMeshmagickメッシュリーダとは異なり、このライブラリはCapytaineにパッケージされていないため、個別にインストールする必要があります。
pip install meshio
meshio メッシュオブジェクトは FloatingBody.from_meshio メソッドで読み込むことができます。
import meshio
mesh=meshio.read("myfile.stl")
body=FloatingBody.from_meshio(mesh, name="My floating body")

この機能は、pygmshがメッシュをmeshioと同じフォーマットで返すので、pygmshを使ってメッシュを生成できる。以下は、pygmshを使ったメッシュ生成の例です（これも別途インストールする必要があります）。

import pygmsh
offset=1e-2
T1=0.16
T2=0.37
r1=0.88
r2=0.35
with pygmsh.occ.Geometry() as geom:
    cyl=geom.add_cylinder([0, 0, 0], [0, 0, -T1],  r1)
    cone=geom.add_cone([0, 0, -T1], [0, 0, -T2], r1, r2)
    geom.translate(cyl, [0, 0, offset])
    geom.translate(cone, [0, 0, offset])
    geom.boolean_union([cyl, cone])
    mesh=geom.generate_mesh(dim=2)
body=FloatingBody.from_meshio(mesh)


3.3 表示とアニメーション

VTK（インストールされている場合）を使って船体を3Dで表示するには、showメソッドを使用する。
body.show()
animateメソッドを使うと、船体に与えられた動きを可視化できる。
anim=body.animate(motion={"Heave": 0.1, "Surge": 0.1j}, loop_duration=1.0)
anim.run()
上記の例では、HeaveとSurgeの線形結合のインタラクティブなアニメーションを提示します。
Jupyter notebookはアニメーションの（非インタラクティブな）ビデオを含めることもできます。
anim.embed_in_notebook(camera_position=(-1.0, -1.0, 1.0), resolution=(400, 300))

3.4 幾何学的変形

既存の船体やメッシュを変換するための関数がいくつか用意されています。
ほとんどの変換方法は2つのバージョンで存在します。
一つは不定形動詞（translate, rotate, ...）の名前で、インプレースの変換です。
もうひとつは過去分詞（translated, rotated, ...）で、同じ変換ですが新しいオブジェクトを返します。
ほとんどの場合、性能は重要ではなく、後者の方法の方がコードのデバッグが若干容易になるため、好まれるはずです。
以下に、利用可能なメソッドのほとんどをリストアップします。これらのメソッドはメッシュにも浮体にも適用することができ、その場合、自由度も変換されます。

# TRANSLATIONS
body.translated_x(10.0)
body.translated_y(10.0)
body.translated_z(10.0)
body.translated([10.0, 5.0, 2.0])
# Translation such that point_a would become equal to point_b
body.translated_point_to_point(point_a=[5, 6, 7], point_b=[4, 3, 2])
# ROTATIONS
body.rotated_x(3.14/5)  # Rotation of pi/5 around the Ox axis
body.rotated_y(3.14/5)  # Rotation of pi/5 around the Oy axis
body.rotated_z(3.14/5)  # Rotation of pi/5 around the Oz axis
# Rotation of pi/5 around an arbitrary axis.
from capytaine import Axis
my_axis=Axis(vector=[1, 1, 1], point=[3, 4, 5])
body.rotated(axis=my_axis, angle=3.14/5)
# Rotation around a point such that vec1 would become equal to vec2
body.rotated_around_center_to_align_vector(
    center=(0, 0, 0),
    vec1=(1, 4, 7),
    vec2=(9, 2, 1)
)
# REFLECTIONS
from capytaine import Plane
body.mirrored(Plane(normal=[1, 2, 1], point=[0, 4, 5]))
上記の方法はすべてPlaneとAxisオブジェクトに適用することもできます。


3.5 結合

メッシュと船体は+演算子で結合できる。
both_bodies=body_1+body_2

演算は連結可能であり、(body_1+body_2)+body_3 は body_1+(body_2+body_3) と等価である。また、通常は関係ない内部的な詳細まで可換である。しかし、2体以上の場合は、代わりにjoin_bodiesメソッドを使うことが推奨される。
all_bodies=body_1.join_bodies(body_2, body_3, body_4)

dofsを持つ2つの浮体がマージされたとき、結果の船体は個々の船体のdofを継承し、新しい名前 body_name__dof_name を持つ。たとえば
body_1.add_translation_dof(name="heave")
body_2.add_translation_dof(name="Heave")
both_bodies=body_1+body_2
assert 'body_1__Heave' in both_bodies.dofs
assert 'body_2__Heave' in both_bodies.dofs

3.6 クリッピング

メッシュや船体はclipメソッドやclippedメソッドでクリッピングできる。幾何学的な変換は、前者がインプレースであるのに対し、後者は新しいオブジェクトを返します。これらのメソッドは引数として Plane オブジェクトをとります。平面は、それに属する点と法線ベクトルによって定義されます。
xOy_Plane=Plane(point=(0, 0, 0), normal=(0, 0, 1))
clipped_body=body.clipped(xOy_Plane)

Planeの法線ベクトルの向きによって、メッシュのどの部分が返されるかが決まるので注意。higher_part=body.clipped(Plane(point=(0, 0, 0), normal=(0, 0, -1)))
lower_part=body.clipped(Plane(point=(0, 0, 0),normal=(0、0、1)))
# body=lower_part+higher_part

keep_immersed_part メソッドは、船体をz=0とz=-hの２水平面でクリップする (デフォルトではインプレース)．
clipped_body=body.keep_immersed_part(sea_bottom=-10, inplace=False)


3.7 重心と回転のDOF

船体の重心は、center_of_mass属性に3つの要素からなるベクトルを代入することで定義できる。
body.center_of_mass=np.array([0.0, -1.0, -1.0])

質量中心はいくつかの流体力計算で使用されます。
回転自由度の定義を除いて、流体力係数には必要ない。
回転自由度を定義するとき、コードは rotation_center, center_of_mass, *geometric_center という属性を探し（この順番）、それらを使って回転軸を定義する。
いずれも定義されていない場合、回転は領域原点(0, 0, 0)を中心に定義する。


3.8 積分の求積の定義

警告
この機能は実験的なものです。現時点では、四角形パネルのみがサポートされる。
BEM問題を解く際、グリーン関数はメッシュ上で積分される必要があります。デフォルトでは、積分はパネルの中心での値を取り、その面積を乗算することで近似されます。より正確な積分を行うには、高次の求積を定義できる。
この機能は、求積法の計算に外部パッケージのquadpyを使用しています。インストールは
pip install quadpy

次に、利用可能な求積法の中から一つを選び、compute_quadratureメソッドに渡します。
from quadpy.quadrilateral import stroud_c2_7_2
body.mesh.compute_quadrature(method=stroud_c2_7_2())
これは必要なときに自動的に使用されます。

警告
メッシュの変形(マージ、クリッピングなど)は求積をリセットする可能性があります。最終的なメッシュに対してのみ計算してください。


4. 静水力学

Capytaineは、与えられたFloatingBodyの静水圧パラメタを計算する。

注
ここで、没水した船体上の関数の積分は、面中心でのデータ関数の和にそれぞれの面積を乗じたものとして近似される。
(1)
ここで、i は面番号、(xi, yi, zi) は i 番目の面領域である。

4.1 静水力学パラメタ

Capytaineでは、各静水圧パラメタに対して個別のメソッドが用意されています。
それらのいくつかは、船体の質量中心を定義する必要があるかもしれません。
これは、以下の例のように属性 center_of_mass を設定することによって行う。

注意
個々の静水圧パラメタを計算する前に、必ず船体水上部を切り取って（クロップして）没水部分のみを残すこと。

import capytaine as cpt
import numpy as np

rigid_sphere = cpt.Sphere(radius=1.0, center=(0, 0, 0), name="my buoy")
rigid_sphere.center_of_mass = (0, 0, -0.3)
rigid_sphere.add_all_rigid_body_dofs()
rigid_sphere.keep_immersed_part()

print("Volume:", rigid_sphere.volume)
print("Center of buoyancy:", rigid_sphere.center_of_buoyancy)
print("Wet surface area:", rigid_sphere.wet_surface_area)
print("Displaced mass:", rigid_sphere.disp_mass(rho=1025))
print("Waterplane center:", rigid_sphere.waterplane_center)
print("Waterplane area:", rigid_sphere.waterplane_area)
print("Metacentric parameters:",
    rigid_sphere.transversal_metacentric_radius,
    rigid_sphere.longitudinal_metacentric_radius,
    rigid_sphere.transversal_metacentric_height,
    rigid_sphere.longitudinal_metacentric_height)


4.2 静水剛性

浮体の静水力学を計算する式は次の通りである。
(2)
ここで
n は表面法線
Vi=unx+vny+wnz はDOFベクトル
Di=∇・Vi は DOF の発散度である。

compute_hydrostatic_stiffness()メソッドは静水力学を計算し、
(DOF count x DOF count) 2D マトリクスを xarray.DataArray として返す。

rigid_sphere.hydrostatic_stiffness = rigid_sphere.compute_hydrostatic_stiffness()

注
剛体DOFは、上記の正確な式で評価できる。しかし, 一般に, 任意のDOFsの発散は不明である。ユーザは Di の値を compute_hydrostatics の入力として渡すことができる。そうでない場合、コードは発散がゼロであると仮定する

elastic_sphere = cpt.Sphere(radius=1.0, center=(0, 0, 0), name="my buoy")
elastic_sphere.center_of_mass = (0, 0, -0.3)
elastic_sphere.keep_immersed_part()
elastic_sphere.dofs["elongate_in_z"] = np.array([(0, 0, z) for (x, y, z) in elastic_sphere.mesh.faces_centers])
dofs_divergence = {"elongate_in_z": np.ones(elastic_sphere.mesh.nb_faces)}
density = 1000.0
gravity = 9.81
elongate_in_z_hs = elastic_sphere.compute_hydrostatic_stiffness(divergence=dofs_divergence, rho=density, g=gravity)
analytical_hs = - density * gravity * (4 * elastic_sphere.volume * elastic_sphere.center_of_buoyancy[2])
print(np.isclose(elongate_in_z_hs.values[0, 0], analytical_hs))
# True

質量が指定されない場合（上記の例のように）、物体は浮力平衡状態にあると仮定される。その質量は、変位した水の体積の質量となる。

非中立の浮体は部分的に実装されています（単一の剛体の場合のみ）。
この場合、FloatingBodyのmass属性を設定することで、船体の質量を与える。

rigid_sphere.mass = 0.8 * rigid_sphere.disp_mass(rho=1000)
rigid_sphere.hydrostatic_stiffness = rigid_sphere.compute_hydrostatic_stiffness()

4.3 慣性行列

compute_rigid_body_inertia() メソッドは、6つの剛体係数を持つ物体の6×6慣性行列を計算する。その他の自由度の慣性係数は、デフォルトでNaNで埋められます。
rigid_sphere.inertia_matrix = elastic_sphere.compute_rigid_body_inertia()

静水力学に関しては、質量属性が指定されていない限り、質量は水の変位質量であると仮定される。
カスタムの行列を提供する。Capytaine で計算されたデータとの一貫性を保つために、dof 名をラベルとして xarray.DataArray でラップすることが推奨される。
elastic_sphere.inertia_matrix = elastic_sphere.add_dofs_labels_to_matrix(np.array([[1000.0]]]))

4.4 全ての静水圧パラメタの計算

各静水圧パラメタを個別に計算する代わりに、compute_hydrostaticsは全ての静水圧パラメタを含むdictを返します。

注意
compute_hydrostatics を使用するために keep_immersed_part を適用する必要はない。

hydrostatics = rigid_sphere.compute_hydrostatics()
print(hydrostatics.keys())
# dict_keys(['g', 'rho', 'center_of_mass', 'wet_surface_area', 'disp_volumes',
# 'disp_volume', 'disp_mass', 'center_ of_buoyancy', 'waterplane_center',
# 'waterplane_area', 'transversal_metacentric_radius',
# 'longitudinal_metacentric_radius' , 'transversal_metacentric_height',
# 'longitudinal_metacentric_height', 'hydrostatic_stiffness',
# 'length_overall', 'breadt h_overall', 'depth', 'draught',
# 'length_at_waterline', 'breadth_at_waterline',
# 'length_overall_submerged', 'breadth_overall_submerged', 'inertia_matrix'])

5. 問題の設定
5.1 xarrayを使ったテスト行列の定義

解くべきポテンシャル流れ問題のパラメタを設定する最短の方法は、xarrayデータセットを構築することです。そして、関数 fill_dataset は自動的にテスト行列で定義されたすべてのシミュレーションを行います。

import numpy as np
import xarray as xr
import capytaine as cpt
body = ...  # Set up the body and its dofs here.
test_matrix = xr.Dataset(coords={
    'omega': np.linspace(0.1, 4, 40),
    'wave_direction': [0, np.pi/2],
    'radiating_dof': list(body.dofs),
    'water_depth': [np.infty],
})
dataset = cpt.BEMSolver().fill_dataset(test_matrix, body)

fill_dataset の第一引数は、実行するケースを指定する xarray.Dataset である。
データセットにパラメタが指定されていない場合は、デフォルト値が使用される（次のセクションを参照）。
第二引数は、単一の FloatingBody か FloatingBody のリストである。後者の場合、リスト中の各ボディは独立して調べられ、出力データセットは1次元多く含む。
fill_dataset は、入力と同じ座標で、追加出力データが入った xarray.Dataset を返す。
座標θがテスト行列に追加された場合、このコードはこれらのθの値に対してKochin関数を計算する。

5.2 LinearPotentialFlowProblem クラス

解くべき問題をより細かく制御するために、LinearPotentialFlowProblemのリストを定義しておくとよい。これは例えば次のように定義される。

problem = cpt.RadiationProblem(body=my_body, omega=1.0, radiating_dof="Heave")
other_problem = cpt.DiffractionProblem(body=my_body, omega=1.0, wave_direction=pi/2)

body 以外はすべてオプションのパラメタです. 以下の表は、それらの定義とデフォルト値です。

パラメタ 説明 (単位) デフォルト値
free_surface 自由曲面の位置 [1] (m) 0(m)
sea_bottom 海底の位置 (m) -Infininity(m)
omega 周波数 (rad/s) 1.0(rad/s)
g 重力加速度(m/s2) 9.81(m/s2)
rho 水の密度(kg/m3) 1000(kg/m3)
wave_direction 波の到来方向 (回折の場合のみ) 0(rad) [2] (1)
radiating_dof 放射するdofの名前 (放射の場合のみ) 最初に見つかったもの

[1] 自由曲面の位置はz=0とz=+infinityの2つしか認められません. 前者は線形ポテンシャル流の通常の場合です. 後者は自由表面を持たない無限ポテンシャル流の領域における物体に対応する.
[2] 波の方向が 0 rad の場合は, x=-infinity から x 軸方向に伝播する波が tox=+infinity に対応する.

警告
Nemoh などの他のソフトウェアと異なり、Capytaine の波動方向はラジアン単位で表現される。

波高は暗黙で1mと仮定される。すべての計算は線形であるため、任意の波高や動きの振幅は、その結果に希望する値を乗じることで取得できます。
以下の属性は、与えられた問題に対して自動的に計算される。

パラメタ 説明 (単位)
深度 水深 h(m)
波数 波数 k(1/m)
波長 波長 lambda=2pi/k (m)
周期 波の周期 T=2pi/omega(s)
無次元数_ω 2*omega^2*h/g (-)
無次元_波数 kh(-)

として取得することができる:
problem.wavenumber
problem.period
# ...


5.3 Nemoh.calのレガシーパラメタファイル

Nemohのレガシーパラメタファイルは、専用の関数で読み込むことができます。
from capytaine.io.legacy import_cal_file
list_of_problems = import_cal_file("path/to/Nemoh.cal")

この関数は、LinearPotentialFlowProblemsのリストを返す。

警告：この機能は実験的なものです。ファイル中のいくつかの設定（自由曲面計算やKochin関数など）は今のところ無視される. 以下のNemoh.calの例を参照してください。

--- 環境 ---
1000 ! RHO ! KG/M**3 ! 流体密度
9.81 ! G ! M/S**2 ! 重力
0 ! 深さ ! M ! 水深
0 0 ! 波浪計測点 ! この線は無視される !
--- 浮遊物の説明
1 ! ボディ数
--- ボディ1
円柱.dat ! メッシュファイル名
540 300 ! メッシュの頂点数、パネル数 ! この行は無視される !
6 ! 自由度数
1 1 0 0 0 0 0 ! サージ
1 0 1 0 0 0 0 ! 揺れ
1 0 0 1 0 0 0 ! ヒーブ
2 1 0 0 0 0 -7.5 ! CdGまわりのロール
2 0 1 0 0 0 -7.5 ! CdGまわりのピッチ
2 0 0 1 0 0 -7.5 ! CdGまわりのヨー
6 ! 一般化力の数
1 1 0 0 0 0 0 ! x方向の力 ! この行は無視される !
1 0 1 0 0 0 0 ! y方向の力 ! この行は無視される !
1 0 0 1 0 0 0 ! z方向の力 ! このラインは無視される !
2 1 0 0 0 0 -7.5 ! CdG x軸まわりのモーメント ! この行は無視される
2 0 1 0 0 0 -7.5 ! CdG y軸まわりのモーメント ! この行は無視される．
2 0 0 1 0 0 -7.5 ! CdG z軸まわりのモーメント ! この行は無視される．
0 ! 付加情報の行数 ! 追加情報は全て無視される !
--- 解く荷重ケース
2 0.1 2.0 ! 波数、最小値、最大値 (rad/s)
1 0.0 0.0 ! 波の方向数、最小、最大 (度)
--- 後処理
0 0.1 10 ! IRF ! この行は無視される !
0 ! 圧力表示 ! この行は無視される !
0 0 180 ! コチン関数 ! この行は無視される !
0 10 100 100 ! 自由曲面の標高 ! この行は無視される !

5.4 コマンドラインインタフェース

警告：この機能は実験的なものです。

Capytaineにはコマンドラインのcapytaineが付属しており、次のように使うことができます。
$ capytaine /ディレクトリ/parameter.cal
引数として渡されたパラメタファイル(Nemoh.calフォーマット)を読み込み、レガシーtecplot出力ファイルをディレクトリpath/to/directory/results/に書き出します。

警告：すでに結果ファイルが存在する場合は、上書きされる。

コマンドに引数が与えられない場合、コードはカレントディレクトリにあるNemoh.calファイルを探します。

6. 解法
6.1 ソルバの設定

BEMSolverクラスは、初期化時に2つの引数（キーワードのみ）を取ります。
from capytaine import BEMSolver
solver = BEMSolver(green_function=..., engine=...)
この2つのオブジェクトについて詳しく説明します。

グリーン関数

Green関数を評価するためのクラスで、AbstractGreenFunctionから派生しています。現在のバージョンでは2つ用意されています．

Delhommeau (デフォルト)
Nemoh で実装されている方法 (参照: [Del87] および [Del89]). 利用可能なオプションの詳細については、ドキュメントを参照してください。

XieDelhommeau
上記の手法の変形で, 自由表面付近でより正確な計算が可能です ([X18] を参照). Delhommeauと同じオプションが使えます。

上級ユーザはグリーン関数を評価するために独自のクラスを書くことができます. クックブックの例を参照してください.

エンジン

MatrixEngineから派生した，相互作用行列を構築するためのクラスです．現バージョンでは2つ用意されています．

BasicMatrixEngine (デフォルト)
Nemohにかなり似たシンプルなエンジンです。このエンジンは、ほとんど最適化されずに完全な行列を構築します。解法を速くするために，反射対称性のみを使用できる．
このオブジェクトは，以下のオプションで初期化できる．
matrix_cache_size (デフォルト: 1)
ソルバは、最後に計算された相互作用行列をメモリに保持します。この設定は、保存される古い行列の数を制御します。0に設定すると、コードのRAM使用量が減りますが、計算時間が長くなる可能性があります。

linear_solver (デフォルト: 'direct')
このオプションは、BEM問題の解法に使用される線形システムのソルバを設定するために使用されます。文字列を渡すと、コードは定義済みのソルバのいずれかを使用するようになります。LU分解を用いた直接法ソルバである「direct」と反復法ソルバである「gmres」の2つが利用可能である。
前者はより堅牢で、計算時間が予測しやすいため、デフォルトで使用されています（バージョン1.4以降）。上級ユーザは、ソルバを gmres に変更したいと思うかもしれません。これは多くの状況で速く、他の状況では完全に失敗します。
また、行列とベクトルを引数にとり、ベクトルを返す任意の関数をソルバに与えることもできます。

import numpy as np
def my_linear_solver(A, b):
        """A dumb solver for testing."""
        return np.linalg.inv(A) @ b
my_bem_solver = cpt.BEMSolver(
   engine=BasicMatrixEngine(linear_solver=my_linear_solver)
   )

このオプションは、例えば反復法ソルバにカスタムプリコンディショニングを適用するために使用できる。

HierarchicalToeplitzMatrixEngine
メッシュの階層構造を利用して、階層的な影響行列を構築する実験エンジンです。
このオブジェクトは以下のオプションで初期化できる。
matrix_cache_size (デフォルト: 1)
上記と同じ。

ACA_distance および ACA_tol
低ランク行列の精度を設定するために使用される，Adaptive Cross Approximation (ACA) のパラメタ．

6.2 レガシーインターフェース

NemohクラスはCapytaineのバージョン1.0におけるメインのソルバクラスであった。後方互換性のため、現在のバージョンでも利用可能です。現在はBEMSolverのサブクラスで、常にDelhommeauのグリーン関数を使用し、バージョン1.0と同じ引数を受け付けます。
BEMSolverの使用を推奨します。

6.3 問題の解法

ソルバが初期化されると、solve()メソッドで問題を解くことができます。
result = solver.solve(problem, keep_details=False)

オプション引数keep_details（デフォルト値：True）は、ソースとポテンシャルの分布を結果オブジェクトに保存するかどうかを制御します。これらのデータは、Kochin関数の計算や自由表面高度の再構成のようないくつかのポスト処理に必要です。しかし、物体にかかる力だけに興味がある場合は、メモリスペースを節約するためにこれらを破棄できる。

問題のリストは次のようにして最適な順序で一度に解くことができる。
list_of_results = solver.solve_all(list_of_problems, keep_details=False)


6.4 並列化

Capytaineは2種類の並列化機能を備えています。
	joblib	OpenMP
シングル(solve)
バッチ処理(solve_all, fill_dataset)

OpenMPによる単一問題解決

単一問題を解く場合、行列の構築と線形代数演算（インストールに応じてBLASまたはMKLを使用）をOpenMPで並列化できる。この機能はデフォルトでインストールされ、オンになっています。使用するスレッド数は環境変数OMP_NUM_THREADSとMKL_NUM_THREADS（通常condaで配布されているIntelのMKLライブラリを使用した線形代数の場合）で制御することが可能です。この環境変数はPythonインタプリタが起動する前に設定する必要があることに注意してください。また、スレッド数を動的に変更したい場合は、threadpoolctl ライブラリで行うこともできます(GH47 も参照)。

joblibによる一括解決

いくつかの独立した問題を解決するとき、それらを並行して解決できる。この機能(バージョン1.4の新機能)は、オプションの依存関係であるjoblibがインストールされている必要があります。メソッド solve_all() と fill_dataset() はオプションのキーワード引数 n_jobs を取り、バッチ処理中に並列で実行するジョブの数を制御します。joblibはユーザーからのフィードバック（ロギングやエラー報告）を妨げる可能性があるため、現在デフォルトで無効になっています。

n_jobs=1（デフォルト）またはjoblibがインストールされていない場合、バッチ処理の並列化は行われません（ただし、OpenMP並列化は有効です）。

n_jobs=-1の場合は、すべてのCPUコアが使用されます(そして、joblibは自動的にOpenMP並列化を無効にするはずです)。

2つの並列化レイヤー（OpenMPとjoblib）は、それぞれ異なる使い方をします。メッシュは比較的小さいが、多くの海況を研究する場合は、joblib並列化を使用する必要があります。一方、メッシュが大きい場合や使用可能なRAMが少ない場合は、joblibの並列化をオフにしてOpenMPの並列化のみを使用することが有益な場合があります。

7. 出力

7.1 LinearPotentialFlowResultからデータセットを構築する

LinearPotentialFlowResultのリストがある場合、それらをxarrayデータセットにまとめると、より便利に後処理を行うことができます。次のような構文で行います。
from capytaine import assemble_dataset
dataset = assemble_dataset(list_of_results)

BEMSolver.fill_datasetメソッドにテスト行列を与えた場合、出力は直接xarrayデータセットになる。
assemble_datasetとfill_datasetは、データセットにもっと情報を格納するために、いくつかのオプションのキーワード引数を受け付けます。
wavenumber (default: False): 入力された波の波数をデータセットに追加する。
wavelength (default: False): 入力された波の波長をデータセットに追加する.
mesh (default: False): データセットのメッシュに関するいくつかの情報(面の数、求積法)を追加します。
hydrostatics (default: True): FloatingBodyに静水圧データがある場合、データセットに追加されます。

注
このコードは、DOFをユーザーから提供されたのと同じ順序に保つように最善を尽くしますが、常にそうであるという保証はありません。しかし、どのような場合でも、ラベルは常にデータと一致する。したがって、dofの名前で値を選択すると、常に正しいものが返されます。
data.sel(radiatating_dof="Heave",influated_dof="Heave")。

また、以下の構文で手動でdofを並べ替えることができます。
sorted_dofs = ["Surge", "Sway", "Heave", "Roll", "Pitch", "Yaw"]。
print(data.sel(radiating_dof=sorted_dofs,influcted_dof=sorted_dofs))する。

注
assemble_datasetで作成したデータセットは、自由表面を持つケースのデータのみを含みます。自由曲面を持たないケース(free_surface=infty)は無視される。

7.2 Bemioからデータセットを構築する

xarray データセットは、NEMOH、WAMIT、AQWA の流体力学出力データを読み込む Bemio パッケージを使って生成したデータ構造から作成することも可能です。これにより、他の BEM コードから生成された流体力学データを Capytaine で後処理できる。
BemioはCapytaineに同梱されていないため、単独でインストールする必要があります。Bemioのベースリポジトリはアーカイブされており、Python 2.7.xとしか互換性がないので、Python 3互換のフォークを使用することが推奨されます（ここで入手可能、または以下の方法でインストール可能）。
pip install git+https://github.com/michaelcdevin/bemio.git

Capytaineを使用してxarrayデータセットを構築するには、当該BEMプログラムからの出力ファイルをBemio data_structures.ben.HydrodynamicDataクラスに読み込み、assemble_datasetでこれを呼び出すことが必要です。例えば、WAMITの.outファイルからxarrayデータセットを作成する場合。

from bemio.io.wamit import read as read_wamit
import capytaine as cpt
bemio_data = read_wamit("myfile.out")
my_dataset = cpt.assemble_dataset(bemio_data, hydrostatics=False)

警告
作成されるデータセットには、元のデータセットで与えられた値から直接計算できる量しか含まれない。このため、Bemioで作成したxarrayデータセットとLinearPotentialFlowResultで作成したデータセットでは、形式は同じでも変数名に若干の違いがある場合がある。例えば、WAMITの.outファイルには慣性モーメントの計算に必要な回転半径が含まれておらず、剛体質量行列が計算できないため、上記の例ではmy_dataset['inertia_matrix']変数は含まれない。

7.3 データセットを NetCDF ファイルとして保存する

assemble_dataset (または fill_dataset) が生成する xarray データセットは NetCDF ファイルフォーマットに近い構造を持っており、簡単にこのフォーマットで保存する。
dataset.to_netcdf("path/to/dataset.nc")
詳細とオプションは xarray のドキュメントを参照。

しかし、いくつか注意しなければならない問題がある。


複素数

netCDF標準は複素数を扱わない。回避策として、capytaine.io.xarray モジュールの助けを借りて、複素数値配列をより大きな実数値配列として保存する。
from capytaine.io.xarray import separate_complex_values
separate_complex_values(dataset).to_netcdf("path/to/dataset.nc")。

これでデータセットの再読み込みができるようになる。
from capytaine.io.xarray import merge_complex_values
dataset = merge_complex_values(xr.open_dataset("path/to/dataset.nc"))

文字列の形式

xarrayでの文字列の扱いに問題があります。radiating_dofやinflucted_dofのような文字列をラベルに持つ座標に影響がある。これらは xarray に NetCDF 文字列オブジェクトとして格納され、NetCDF ファイルに書き込むことができるが、 Python 文字列として格納され、一般的な Python オブジェクトとして NetCDF ファイルには書けない。
問題は、xarrayライブラリが警告を出さずに一方から他方へ変更されることがあることである。これはデータセットをエクスポートしようとすると ValueError: unsupported dtype for netCDF4 variable: object というエラーにつながる。

これは、データセットをエクスポートする際に、文字列を正しいフォーマットに明示的に変換することで解決できる。

dataset.to_netcdf("dataset.nc",
  encoding={'radiating_dof': {'dtype': 'U'},
'influenced_dof': {'dtype': 'U'}})

こちらのGithub issueもご参照ください。


7.4 レガシーなTecplotファイルとしてデータを保存する

警告：この機能は実験的なものです。

以下のコードは、RadiationCoefficients.tecとExcitationForce.tecというファイルをNemoh 2.0と同じフォーマットで書き込みます。

from capytaine.io.legacy import write_dataset_as_tecplot_files
write_dataset_as_tecplot_files("path/to/directory", dataset)

8. 後処理

8.1 自由表面の高さ

自由表面の高さを計算するために、まずFreeSurfaceオブジェクトを初期化する。
from capytaine import FreeSurface
fs = FreeSurface(x_range=(-10, 10), nx=10, y_range=(-5, 5), ny=10)
上記のコードは、自由表面の規則的なメッシュ 10x10セルを生成する。このオブジェクトは、自由表面の高さを得るためにresultsオブジェクトと一緒にソルバで使用する。

fs_elevation = solver.get_free_surface_elevation(result, free_surface)
出力はnumpy配列である．周波数領域における各点での自由表面の高さを複素数として格納した（この例では 10x10の配列）。

結果オブジェクトはオプション keep_details=True で計算する必要がある.
ソルバは結果オブジェクトを計算したものである必要はありません.

乱れのない到来波(Airy波)は次のように計算されます。
incoming_waves = fs.incoming_waves(DiffractionProblem(omega=1.0, angle=pi/2))
3Dアニメーションの使用は Cookbookの例を参照してください。

8.2 インピーダンスとRAO

固有インピーダンスは、流体力学、静水圧、慣性特性に基づいて計算する。

import numpy as np
import xarray as xr
from capytaine import BEMSolver
from capytaine.bodies.predefined.spheres import Sphere
from capytaine.post_pro import impedance
f = np.linspace(0.1, 2.0)
omega = 2*np.pi*f
rho_water = 1e3
r = 1
sphere = Sphere(radius=r, ntheta=3, nphi=12, clip_free_surface=True)
sphere.center_of_mass = np.array([0, 0, 0])
sphere.add_all_rigid_body_dofs()
sphere.inertia_matrix = sphere.compute_rigid_body_inertia(rho=rho_water)
sphere.hydrostatic_stiffness = sphere.compute_hydrostatic_stiffness(rho=rho_water)
solver = BEMSolver()
test_matrix = xr.Dataset(coords={
    'rho': rho_water,
    'water_depth': [np.infty],
    'omega': omega,
    'wave_direction': 0,
    'radiating_dof': list(sphere.dofs.keys()),
    })
data = solver.fill_dataset(test_matrix, sphere_fb,
                           hydrostatics=True,
                           mesh=True,
                           wavelength=True,
                           wavenumber=True)
Zi = impedance(data)

慣性と剛性をinertia_matrixとhydrostatic_stiffnessというbodyの属性に割り当てることに注意してください。これらは、インピーダンス行列を計算するために fill_dataset 関数と impedance 関数が期待する名前です。

加振伝達関数response amplitude operator (RAO)を取り入れる
from capytaine.post_pro import rao
rao = rao(data)

9. 規約と他のコードとの相違点

備考

後日追記される。

他の多くのコードと異なり、Capytaineでは角度（入射波の方向など）はラジアン単位で示される。

9.1 WAMITとの関係

Capytaineでは，複素数値の振幅（位相）は，次の規則で定義される。
x=Re(X exp(-iwt))
WAMITはRe(X exp(+iwt))で定義している．

9.2 NemohとAquadynに関して

Capytaineは、ほとんどがNemohと同じ規約に従っており、Aquadynも同じ規約です。主な例外は、NemohとCapytaineでは加振力の位相角がAquadynと逆になっていることである。

10 Cookbook