20200201

DualSPHysicsとBlenderで動画を作った

DesignSPHysics(FreeCADのアドオンで、中にDualSPHysicsが入っている)とBlender(+VisualSPHysics)を使って動画を作成した。

試した

粒子で豆腐っぽい形の液体を作って、それを下に落とす。boundaryは描いてない。妥当な結果か分からない。


以下、やり方等のメモ。

FreeCAD+DesignSPHysicsを使ってSPH計算

インストール

  1. FreeCADをインストールする。

  2. DesignSPHysicsをInstallation instructionを参考にインストールする。

    a. GitHubからダウンロードする。

    b. フォルダの名前をDesignSPHysics に変更して、 Windowsなら、%appdata%/FreeCAD/Mod Linux なら ~/.FreeCAD/Mod に入れる。

    c. 同フォルダ内のDesignSPHysics.FCMacro を、macroフォルダに入れる (Windowsなら%appdata%/FreeCAD/Macro Linuxなら; ~/.FreeCAD/Macro 。)

    d. FreeCADを立ち上げて、マクロ(M)>マクロ…で、DesignSPHysics.FCMacroと表示されていれば正しくインストールされている。

09.png

DesignSPHysics.FCMacroの使用方法

公式の動画がYouTubeに上がっていた。

始める前に見たほうがいいかもしれない。

バグ?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をロードして結果が確認できる。

10.png

vtkの準備

まずはvtkファイルの準備。これがないと結果が見れない。

  • 粒子で出力したいなら

    • PartVTKOut4でvtkを出力→点群が出てくる。
      • 連番のvtkが出てくる。
  • 液体っぽい表現で出力したいなら

    • IsoSurface4でvtkを出力→流体の表面を模したメッシュが出てくる。
      • Blenderのpythonスクリプトでmarching cube(後述)する必要はない。

ParaView

  1. ParaViewを起動

  2. VTKをすべて選択して、ParaViewにドラッグドロップする。

  3. 左側のメニューの目のマークを押して表示。
  4. 緑色の三角ボタン(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%/Surface

5行目のnameの部分は、解析ファイルごとに書き換えないといけない。

Blender+VisualSPHysicsを使った流体のレンダリング

VisualSPHysicsを使うと、連番vtkファイルを、アニメーションのシーケンスにしてくれる。

isosurfaceなら、そのままアニメーションレンダリングが可能。

点群のままならば、レンダリングされないので注意。


11.png12.gif

必要なもの

情報は随時更新されるはずなので、最新のものに従っていただくとして、このサイトを参考に必要なものをインストールする。

  • 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月だが。)

posted by yuchan at 07:00 | Comment(0) | SPH
この記事へのコメント
コメントを書く
お名前:

メールアドレス:

ホームページアドレス:

コメント: