Calculix 10本ノック: 1本目

Windows上で動かすLinux（WSL）で、Calculixを使えるようにした。
WSLでは基本GUIがないから、ひと手間必要だった。

Calculixとは

フリーの解析ソフト。
OSは基本的にLinux（Windows版もあるが、なぜかcgxの起動が遅い（僕だけ？））。
Abaqusと書式が似ているらしい。

• ccx ソルバー
• cgx　プリポスト（プリとしてあまり使われない。プリにはgmeshかFreeCADがスタンダード？）

〇本ノックとは

例題を繰り返すやつ。Qiitaとかでよく見る。

ノック1本目: ccx, cgxのインストール

windows上でLinuxが動く状態にする

windows上でUbuntu（Linuxの有名ディストリビューション）が動く（いわゆる Windows Subsystem for Linux）状態にする。

ubuntu

とすると、

yusuke@yusuke-PC-2:~$

と表示されて、ubuntuが起動した状態になればOK。

ccx、cgxをインストールする

この動画を参考にインストールできる。要するに、

sudo apt-get install calculix-ccx
sudo apt-get install calculix-cgx

をする。

インストールが終わったら、

cgx -b a.fbd

上記コマンドで実行できると思いきや、

yusuke@yusuke-PC-2:~$ cgx -b aaa
on a Linux machine, nodename yusuke-PC-2, release 4.4.0-18362-Microsoft, version #476-Microsoft Fri Nov 01 16:53:00 PST 2019, machine x86_64
parameters:3 arguments:2
freeglut (cgx): failed to open display ‘’

エラーで何も動かない。

エラーの解決策

このエラーの解決策は、このサイトのウィンドウが出ないの箇所 )に載っていて、

• Xming-6-9-0-31-setup.exe
• Xming-fonts-x-x-x-x-setup.exe(x-x-x-xはバージョン)

をダウンロードしてインストール 後に、下記のコマンドで環境変数の設定して、

export DISPLAY=localhost:0.0

再度cgx

cgx -b aaa.fbd

とすると無事に、cgxが立ち上がる。（ちなみに、Windowsを起動するたびに、Xmingを立ち上げる必要がある。Xmingなしだとやはり同じエラー。後述のbash.bashrcへの記入（Linuxのviとかで）も必要。）

yusuke@yusuke-PC-2:~$ cgx -b aaa.fbd
on a Linux machine, nodename yusuke-PC-2, release 4.4.0-18362-Microsoft, version #476-Microsoft Fri Nov 01 16:53:00 PST 2019, machine x86_64
parameters:3 arguments:2
GL_MAX_EVAL_ORDER:30
ERROR: The input file “aaa.fbd” could not be opened.

↑こんな表示が出る

（重要）exportは、永続性がない。

先ほどexportコマンドで、環境変数を設定した。
しかし、このままではexportは毎回入力しなければならない。
それを回避したいなら、

/etc/bash.bashrc

に、export文を追記する。 /etc/bash.bashrc は、毎回読み込まれるらしい。

サクラエディタで /etc/bash.bashrcを触ってはいけない

/etc/bash.bashrcの編集は、必ずLinuxのコマンドラインからviなどを使って行う。

上記フォルダを、windowsGUIからサクラエディタを使って編集してはいけない。

Windows上での場所が分かったからって、このLinuxファイルをWindowsのツールでいじってはいけない。

サクラエディタとかwindowsから編集すると、ubuntu起動時に、

-bash: /etc/bash.bashrc: Permission denied

と表示されて、bashが読み込まれなくなる。

ユーザー名の色も消える。

解決策は、

> ubuntu$ chmod 644 /etc/bash.bashrc

ここに載ってた。 https://tutorialmore.com/questions-114467.htm

 

Linuxのフォルダの場所

https://qiita.com/kalafinalice/items/70a76d35398ab11af778 にあったのだが、

/etc/bash.bashrcは、Windows上では

“C:\Users\yusuke\AppData\Local\Packages\CanonicalGroupLimited.UbuntuonWindows_79rhkp1fndgsc\LocalState\rootfs\etc\bash.bashrc”

にあった。yusukeの部分は各ユーザー名により違う。

 

ccx

ccxはGUI無いから、余計な事しないでも立ち上がる。

Windows上のファイルにアクセスするためには、

例えばデスクトップにファイルを置いて

ccx -i "/mnt/c/Users/yusuke/Desktop/eyebar"

c:￥にアクセスするために、

/mnt/c

と書き換えなくてはならないのが結構面倒くさい。

解析実行に関しての詳細は、ノック2本目で書く。

​

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と表示されていれば正しくインストールされている。

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

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なら、そのままアニメーションレンダリングが可能。

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

必要なもの

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

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

Blenderのpythonでクオータニオン（Quaternion）を勉強した

Blenderのpythonでクオータニオン（Quaternion）を勉強した。しかし、BlenderのQuaternionの機能は一切使わない。Quaternionを使った計算方法、モデルの操作方法が知りたいからだ。

まとめ

先に、まとめを書く。

クォータニオンで

• オブジェクトの回転ができる。

• 回転が、任意のベクトルでできる。

注意点

• アフィンと違い、行列の積ではない。Quaternion同士の積は、独自ルールのかけ算が要る。
• 原点を中心として回転になるので、原点を通らないベクトルで回転させたいときは、座標系の変換がいる。

やりたい事、やりたくない事

• やりたくない事
  • Quaternionの数学的な理解。
• やりたいこと
  • Quaternionの数学的操作。

ググったら、

blenderのpythonを使って、オブジェクトを回転させている記事はあった。)

こことかも結構わかりやすくて面白かった。

しかし、オブジェクトの操作のAPIが、クォータニオンを直接受け取れるから実行できるコマンドの記事だった。クォータニオンを使って、計算をしているわけではなかった。
クォータニオンが、回転のどこら辺に効いているのかを知りたい。

そこで、自分で用意した3×3×3＝27個の点の回転をクォータニオンで計算する。

クォータニオンは…

その前に、クォータニオンについて、参考書やネットなどで情報収集した。

• クォータニオンは複素数の拡張版。実数一つに、虚数三つ、計四つ。これは暗記した知識として持っている。「先頭が角度で、その後にx,y,zに関わる何か？」な状態。何も具体的にイメージできない。
• プログラムとかだと、独自の型をもったクォータニオン型の変数として表されることが多い。行列みたく、単なる長さ４の配列じゃなめなのか。
• 回転を表すクォータニオンはサイズが１らしい。だったらサイズ１以外だったらどうなっちゃうのか。

そもそも、回転を表すクォータニオンを準備できたとして、どうやって物を（例えばxyz座標上の点を）回転操作するのか。掛け算？

「どうせアフィンみたく行列の掛け算でできるんだろ？もったいぶるなよ。」とか思ってたが違った。

Blenderのモデル

下記コードでできたモデルを使った。

coord = [[[[i,j,k] for i in range(3)] for j in range(3)] for k in range(3)]

for i in range(3):
    for j in range(3):
        for k in range(3):
            bpy.ops.mesh.primitive_uv_sphere_add(size=1, view_align=False, enter_editmode=False, location=(coord[i][j][k][0]*3, coord[i][j][k][1]*3, coord[i][j][k][2]*3), layers=(True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False))
            mat = bpy.data.materials.new("Material_%s_%s_%s"%(coord[i][j][k][0]+1, coord[i][j][k][1]+1, coord[i][j][k][2]+1))
            bpy.context.object.data.materials.append(mat)
            bpy.context.object.active_material_index = 0
            bpy.context.object.active_material.diffuse_color = (coord[i][j][k][0]/3.0, coord[i][j][k][1]/3.0, coord[i][j][k][2]/3.0)
            bpy.context.object.data.name = "ball_%s_%s_%s"%(coord[i][j][k][0]+1, coord[i][j][k][1]+1, coord[i][j][k][2]+1)
            bpy.context.object.name = "ball_%s_%s_%s"%(coord[i][j][k][0]+1, coord[i][j][k][1]+1, coord[i][j][k][2]+1)

拾ってきた式

クォータニオンの操作は、この本を参考にした。。面白かった。

クォータニオン

クォータニオンは、こんな感じで4つ1セット書かれている。

座標

位置ベクトルをクォータニオンで表す。

座標が[x,y,z]だとしたら、

クォータニオンは[0,x,y,z]と最初に0をつけるだけ。（本当は0+xi+yj+zkと書くべき？）

回転を表すルール

上記が成立するクォータニオンは、回転を表現できる。

回転を表しつつ、上記ルールを守るための、三角関数を使ったクォータニオン

「二つしかないからクォータニオンじゃないじゃん！」→実部と虚部を分けて書いてあるだけ。

• 右側（虚部）の数値は、虚部で分け合って使う。
• 分け合い方で、回転軸のベクトルを作る。
• さらに上記「二乗して合計して1」というルールも守って作らなければならない。
• は、逆四元数。虚部の正負がひっくり返った四元数。

cosとsinがどんな塩梅で合計１になりそうかは、ちょっと分からん。

↑こんな感じで配分してみた。

1/3は、三等分だから。

√2は、わからん。

とにかく、√2/3しとくと、「二乗して合計して1」ルールを満たす。

回転操作

クォータニオンのかけ算

回転にはクォータニオンのかけ算を使う。

クォータニオンのかけ算のルールは、以下の式。独自ルール。甘くはなかったか。

これは、関数を使って実現しておく。スカラー,i,j,kは別々に計算して返せばいいだろう。自作なので、わざわざクォータニオン型とか作らんで、配列でやっとけばいっか。自分がわかってればいいもんね。

def quakakezan(q=[1,2,3,4],r=[5,6,7,8]):
    a=q[0]*r[0]-q[1]*r[1]-q[2]*r[2]-q[3]*r[3]
    b=q[0]*r[1]+q[1]*r[0]+q[2]*r[3]-q[3]*r[2]
    c=q[0]*r[2]+q[2]*r[0]+q[3]*r[1]-q[1]*r[3]
    d=q[0]*r[3]+q[3]*r[0]+q[1]*r[2]-q[2]*r[1]
    e =[a,b,c,d]
    return e

想定外の変数が入ってきたら、エラーで止まるかも。
最初から値が入っているのは、自分がどんな値を入れたらいいか忘れそうだし、関数のすぐ下にコメントアウトで説明書き入れたって、どうせ読まないから、あえて入れちゃう。

回転操作

1. 座標の位置クォータニオンを作る。

1. 前述、回転のためのクォータニオンを作る。
2. 2.を、前と後ろから、かける

もちろん、四元数のかけ算のルールで計算する。

テクスチャを貼っていないから気づかないだろうけど、ボールは回転していない。同じ向きを画面に向けたまま。
BlenderのQuaternionの機能を使って回転させると、ボールの向きも変わるが、それはしなかった。

コード全体は、YouTubeの説明のところに書いてあるので、ぜひ見てください。

振り返ってみると

テクスチャを貼って、オブジェクトのvertexを回転させるのでもよかった気がする。

が、あえてオブジェクトをすべて動かすことで、何か理解した気になった。

Shade3Dで画面スクロールするスクリプト

やったこと

​ Shade3D ver.16 professionalを使い続けています。毎バージョンの最終アップデートで、画面のスクロールができなくなって困っていたので、画面移動のスクリプトを書きました。

Shade3Dは次のバージョンが出ると画面のスクロールができなくなる。

​ 初めて触った3Dソフトがshade（今はshade3D）なので、自分はshadeが好きです。

​ 只、このソフトは、新しいバージョンが出るたびに画面のスクロールができなくなるupdateを配布して、ソフトの利便性を壊してしまいます。多分、意図して？かしないでか？わざと？指定座標位置に移動するようなスクリプトの紹介記事だけ、スクリプト配布か排除されているくらいの徹底ぶり。

​ 現在はサブスクリプションに変わったので、こういうことをしないと思います。amazonで前のバージョンの購入はできないようになってしまいました。フォトショとか3DCGのソフトって、大きな機能のアップデートがないと、新しく買ってもらえなくなる。だから、サブスクリプション形式の売り方になるんだと思っている。もう、あっと驚く機能の追加はないのかぁ。

導入方法、使用方法、欠点など

​ そこで、やや利便性は劣りますが、最後にクリックした位置に画面を移動させるpythonスクリプト（一行）を書いてみました。

• 導入手順は下の通りになります。

1. 下記スクリプトをテキストファイルとして、適当な場所に「scroll.py」などという名前で保存する。

xshade.scene().scroll(xshade.scene().cursor_position)

1. Shade3Dを立ち上げて、表示(V)>ショートカット(R)と選択し、ショートカット画面を開きます。

2. ショートカットを設定したいボタンをクリックして、スクリプト(S)>参照…と選択する→ファイルダイアログ開く。

3. 先ほど書いた「scroll.py」を選択する。

• 使用方法は、

クリック→設定したショートカットキーで実行されるはず。

• 欠点は、

すべての図形ウィンドウが動いてしまうことです。つまり四面図の内、レンダリングのカメラ画面以外は全て移動されてしまいます。

画像は、著作権とか怖そうなので貼らないことにしました。そういう怖さがあるから、商用ソフトの使用方法ってネットに情報が載ってないんだと思う。この記事もアウトなんだろうか。

今月は引っ越しで忙しかった

引っ越しをした。

元居た場所から200mほどしか離れていなかったため、台車を買って自分で運んだ。

さすがに大きな白物家電は、業者に頼んだ。

結論から言うと、自分で荷物を運ぶのは時間の無駄。

体力的にきつかったとかは問題ではない。

自分ですべて運ばなければならないので、時間的な拘束がある。

それがエラいプレッシャーになるし。面倒くさい。面倒くさくて手が止まる。

あの時間があったとしても、大したことはしてなかったと思うが、引っ越しは業者に頼むのが賢明。

学んだ。