試した
粒子で豆腐っぽい形の液体を作って、それを下に落とす。boundaryは描いてない。妥当な結果か分からない。
以下、やり方等のメモ。
FreeCAD+DesignSPHysicsを使ってSPH計算
インストール
FreeCADをインストールする。
DesignSPHysicsをInstallation instructionを参考にインストールする。
a. GitHubからダウンロードする。
b. フォルダの名前を
DesignSPHysicsに変更して、 Windowsなら、%appdata%/FreeCAD/ModLinux なら~/.FreeCAD/Modに入れる。c. 同フォルダ内の
DesignSPHysics.FCMacroを、macroフォルダに入れる （Windowsなら%appdata%/FreeCAD/MacroLinuxなら;~/.FreeCAD/Macro。）d. FreeCADを立ち上げて、マクロ(M)>マクロ…で、DesignSPHysics.FCMacroと表示されていれば正しくインストールされている。
DesignSPHysics.FCMacroの使用方法
始める前に見たほうがいいかもしれない。
バグ？FreeCAD0.18、DesignSPHysics v0.6.0.1912-12-3
円柱をx500mm,y500mmに置いて、Run GenCaseした形状をParaViewで確認したところ、おかしな形状になった。
xmlを見ると、下のポイントのx,y座標が、(500,500)ではなく（0,0）のままだった。
xmlを手動で修正して、手動でGenCaseする必要がある。この場合、FreeCAD上のRun GenCaseではだめ。
cmdで手動で”~~\DualSPHysics_v4.4\bin\windows\GenCase4_win64.exe”を実行。
やり方は、RUN_DIRECTORYフォルダにあるexampleの.batファイルを参考に、自分でcmdで操作。
出力ファイルを作ったら、FreeCADのRUNでも解析できる。
（FreeCADに頼らない解析実行の勉強にはなる。）
より難解な形状を作ったり解析したりするには、DesignSPHysics使うしかない？バグあったけど、XMLの手入力はさすがに厳しい。形状が複雑になるにつれ、手作業でxmlを書くのは難しくなる。
ParaViewで解析結果の確認
上の動画はBlenderで作ったけど、結果の確認は基本ParaViewで行う。
ParaViewにvtkをロードして結果が確認できる。
vtkの準備
まずはvtkファイルの準備。これがないと結果が見れない。
粒子で出力したいなら
- PartVTKOut4でvtkを出力→点群が出てくる。
- 連番のvtkが出てくる。
- PartVTKOut4でvtkを出力→点群が出てくる。
液体っぽい表現で出力したいなら
- IsoSurface4でvtkを出力→流体の表面を模したメッシュが出てくる。
- Blenderのpythonスクリプトでmarching cube（後述）する必要はない。
- IsoSurface4でvtkを出力→流体の表面を模したメッシュが出てくる。
ParaView
ParaViewを起動
VTKをすべて選択して、ParaViewにドラッグドロップする。
- 左側のメニューの目のマークを押して表示。
- 緑色の三角ボタン（Play）を押すとアニメーションされる。
paraviewは、 科学技術可視化 のためのソフトなので、レンダリングは充実していない。
paraview作ってるロスアラモス国立研究所って、アメリカのニューメキシコ州にあるマンハッタン計画で有名な研究所らしい。知らんかった。
FreeCADを使わずに計算を実行する例
DesignSPHysicsのフォルダ内に、examplesが入っている。←FreeCADを使わずに計算を実行している。
DualSPHysicsは、コマンドラインから実行するソフト。
examplesでは、コマンドライン作業をバッチファイルを使って実行させている。
examplesのフォルダに入っている.batファイルがバッチファイル。バッチファイルは、cmdの操作を記述できる。
バッチファイルで作業を簡略化
バッチファイルはダブルクリックして使う。テキストエディタで編集して、処理の順序を変えたりもできる。面倒なコマンドプロンプトの作業もダブルクリックだけで終わる。バッチファイルとは。 )
examplesの.batファイルは、プログラムファイルの場所が正しく書かれていないため、自分で書き直す必要がある。
.batファイルで、isosurfaceを作るには
例えば、isosurfaceの作成を、.batファイルで書くと
@echo off
rem "name" and "dirout" are named according to the testcase
set name=1
set dirout=%name%_out
set isosurface="X:\app\SPH\DualSPHysics_v4.4\bin\windows\IsoSurface4_win64.exe"
%isosurface% -dirin %dirout% -saveiso %dirout%/Surface5行目のnameの部分は、解析ファイルごとに書き換えないといけない。
Blender＋VisualSPHysicsを使った流体のレンダリング
VisualSPHysicsを使うと、連番vtkファイルを、アニメーションのシーケンスにしてくれる。
isosurfaceなら、そのままアニメーションレンダリングが可能。
点群のままならば、レンダリングされないので注意。
必要なもの
情報は随時更新されるはずなので、最新のものに従っていただくとして、このサイトを参考に必要なものをインストールする。
- Blender(2.8以上 2.7なら導入に相当苦労する)
- Python3.?(cmdでpythonと打った時に2.7が出る状態だと、導入できなかった。Blenderの中のPythonにPATHを通したらできた。)
Visual Studio 2015 redistributable package
visualSPHysics (エラーが起きて導入難しい時は、Forkされた古いほうを使う。)
- Blenderが2.7 なら、最新版は無理。（Pythonが2.7でも無理かも？）
- OpenGLのバージョン問題でBlender2.8が導入できなかった。
visualSPHysicsの使い方
Blenderのアドオンを無事導入できたら、
3D Viewの中にマウスカーソルを入れて、Shift+A>Mesh>DualSPHysics objectを押すとファイルダイアログが出現。
- 左下に、操作できるGUIある。
- 連番ファイルは、どれか一つを選ぶだけで全部選択してくれる。
アプライすると、アニメーションのシーケンスが下に出ているかと思います。
繰り返すけど、inportされたのが点群データだった場合、Blenderでレンダリングされない。（isosurfaceを出力しよう！）
Blenderのvertexとして結果をinportし、marching cubeでメッシュしてレンダリング
！苦行！なんでそんな事やるの！？isosurface4でmarching cube法で処理された結果出力できますけど？（マニュアル（DualSPHysics_v4.0_GUIDE.pdf）の11.6 Surface representationに書いてありました。）
それでもやりたい人向けに、書き残しておく。
marching cube法を試すには
マーチングキューブ法は、YouTubeの解説が充実している。
githubからmarching cubeのコードをダウンロードして使う。使い方ものってた。
しかし、これはmain関数を書き直す必要があった。理由は、点群データをメッシュする用途で書かれていなかったから。
書き直した部分だけ抜粋。
import mathutils
def main(a):
print("start calculation of isosurface")
bpy.context.scene.objects.active = bpy.data.objects[a]
obj = context.object
mesh = obj.data
size = len(mesh.vertices)
kd = mathutils.kdtree.KDTree(size)
for i, v in enumerate(mesh.vertices):
kd.insert(v.co, i)
kd.balance()
def scalarfield(pos):
co_find = (pos[0],pos[1],pos[2])
co, index, dist = kd.find(co_find)
return dist
p0=-0.1,-0.1,-0.1
p1=1.1,1.1,1
res=77#200
resolution=(res,res,res)
isolevel=0.02
start = time.time()
isosurface(p0,p1,resolution,isolevel,scalarfield)
elapsed = time.time()-start
print("end test %r"%elapsed)
del obj, kd
###########間は省略###############
if __name__=="__main__":
# __
# ||
# ||
#_||_
#\ /
# \/ENTER OBJECT NAME.
a="ENTER_OBJECT_NAME"
bpy.ops.object.select_all(action='DESELECT')
bpy.data.objects[a].select = True
main(a)かなり遅い。
結論
Blenderで液体っぽく表現したいなら、IsoSurface4を使いましょう。そうすれば、isosurfaceのvtkを出力してくれます。BlenderのPythonでやる必要はありません。これに気付かず、正月休みを無駄にしました。
（正月に記す。出すのは2月だが。）
