Spectral Element Method (SEM)を使って計算してみる(その3)

本当に正しいのかどうかはよくわからないものの,とりあえず,半無限弾性体を伝播する地震動の計算は最後まで走りました。SEMのコード(EFISPEC)のチュートリアルでは2層の成層構造の例題が示されています。先の計算ではチュートリアルの地盤モデルの表層と半無限の基盤の物性を同じに設定して実質的には単なる半無限弾性体として計算をして比較をしました。

今度は,チュートリアルと同じ2層地盤に挑戦してみます,って言うほどたいへんな話じゃないはずなのですが,よくわからんのでいろいろハマります。

1. モデル作成
何がたいへんかといって,メッシュ切りがたいへんです。チュートリアルでは表層を800mメッシュ,基盤層を2400mメッシュで切っています。しかし,Gmshを使ってボリュームごとに異なるサイズを切る方法がよくわかりません。もちろん,四面体要素でよければGmshはそれらしく切ってくれますが,EFISPECのほうがそれに対応していませんからいかんともしがたいものがあります。

難しいことを考えるのはやめて,まずは同じサイズの六面体要素を切ることにしました。先の計算では,地盤モデルである直方体を構成する全ての点の座標を入れて,それをもとに体積を構成する,というもっとも原始的な方法を使いました。しかし,このような簡単なモデルの場合,一つの平面を定義してそれを特定の方向へ引き延ばす(スキャン? Extrudeと呼ぶようです)ことで三次元モデルに拡張することがよく行われます。今度はそのような方法を試してみることにします。

1.1 ポイントの作成
まず,Extrudeするもとになる面を構成します。その面を構成する点を作成します。xz平面上に面を構成してその面をy方向にExtrudeしようという目論見です。そのため,xz平面上の点を定義します。

Point(1) = {-29000, -29000, 0, 800.0};
Point(2) = {29000, -29000, 0, 800.0};
Point(3) = {29000, -29000, -34000.0, 800.0};
Point(4) = {-29000, -29000, -34000.0, 800.0};
Point(5) = {-29000, -29000, -1000, 800.0};
Point(6) = {29000, -29000, -1000, 800.0};

前の計算との違いはPoint 5と6が増えていることです。この二つの点を繋いだラインが表層と基盤面の間の境界面になります。

1.2 ラインの作成
面を構成するときには,反時計回りが基本なので,それを意識してラインをつないでいきます。Line 2が境界面になる予定です。

Line(1) = {1, 5};
Line(2) = {5, 6};
Line(3) = {6, 2};
Line(4) = {2, 1};
Line(5) = {5, 4};
Line(6) = {4, 3};
Line(7) = {3, 6};

1.3 ラインループの作成
反時計回りにラインループを構成します。xz平面だけ考えているので,地盤モデルのxz断面だけ作ればOKです。Line Loop 21が表層,Line Loop 22が基盤層です。Line Loop 22を構成するLine 2がマイナスになっているのはLine Loop 22で反時計回りにループを構成するにはLine 2の向きが反対なので左向きのベクトルとなるようにするためです。

Line Loop(21) = {1, 2, 3, 4};
Line Loop(22) = {5, 6, 7, -2};

1.4 面の作成
作成したラインループを用いて面を作成します。これはラインループをそのまま面として定義すれば終わりです。

Plane Surface(31) = {21};
Plane Surface(32) = {22};

1.5 面から体積へ拡張
ここまででxz平面上の地盤モデルの断面が完成したので,これをy方向に引き伸ばして体積を作成します。このためにExtrudeを使います。

Extrude{0,58000,0}{
Surface{31};
Surface{32};
}

Extrudeのあとのカッコの中はどの方向に面を引き延ばすかを指定しています。この例では,y方向に58000だけ引き延ばす,ということを示しています。
このようにして二つの体積をもつ直方体が完成します。面を引き延ばすことで自動的に体積が作られるため面ループや体積を改めて定義する必要はありません。

これはこれで結構なのですが,Extrudeによって作成されたポイントやライン,面,体積の名前(番号)がよくわかりません。

なんらかの番号が自動的につけられているのですが,どういう規則で番号をつけているかわからないので,これは結構,ハマります。なぜなら,番号がわからないとこの後のステップで必要になるTransfiniteができないからです。

しょうがないのでGUI上で番号をコツコツ調べます。

ここまでで作成したgeoファイルをダブルクリックするとGmshが起動して,何か描かれているのでマウスで向きを適当に変えるとそれらしい直方体になっていることがわかります。この状態で,マウスをポイントやラインの上に載せるとポイントやラインの情報がポップアップで表示されます。数が多くなると覚えきれないので,ノートにメモしておかないと大混乱になります。

ところが,面や体積はマウスを載せても何も表示されません。本当はメニューのどこかを使って面や体積の番号を表示する方法があるのかもしれませんが,よくわからないので場当たり的に対応します。左ペインのメニューでModule -> Physical groups -> Add -> Surfaceと選んでいくと面を表す点線が表示されます。その点線の上にマウスを載せると面番号がポップアップで表示されます。同様にVolumeを選ぶと黄色い球が表示されて(この場合は2つ),この球の上にマウスを載せると体積番号が表示されます。体積は安直に1, 2という番号が振られていますが,面についてはパッと見ただけではよくわからない脈略のない番号がついています。

このようにして調べて見ると以下のようになっていることがわかりました。丸,四角,三角の数字はそれぞれ,ポイント,ライン,面の番号を表しています。ライン番号の横の数字はラインのベクトルの向きを示しています。こうやって改めて見ても番号の振り方にはなんだか脈略がありません。
numbering.png


1.6 面のtransfiniteの設定
Extrudeによって作成された面をメッシュに切るときに四角形で切られるようにtransfiniteを設定します。面番号は前のステップで調べているのでその番号を指定すればOKです。

Transfinite Surface {31};
Transfinite Surface {32};
Transfinite Surface {49};
Transfinite Surface {71};
Transfinite Surface {54};
Transfinite Surface {76};
Transfinite Surface {41};
Transfinite Surface {63};
Transfinite Surface {53};
Transfinite Surface {45};
Transfinite Surface {67};

さらに,メッシュが四角形となるように,recombineします。

Recombine Surface{31,32,49,71,54,76,41,63,53,45,67};

1.7 体積のtransfiniteの設定
最後に体積をtransfiniteして六面体要素が構成されるようにします。六面体メッシュを切るためには,面がそれぞれ4辺からなっていて,かつ体積が6面からなっていなくてはならない,という制約条件があります。そのうえで,transfiniteするに際に表面を構成するポイントの番号を反時計回りに指定してから裏面を構成するポイントの番号を反時計回りに指定しなくてはなりません。

Transfinite Volume{1} = {1,5,6,2,7,8,12,16};
Transfinite Volume{2} = {5,4,3,6,8,18,22,12};

2. メッシュ切り
geoファイルが正しく書けていればメッシュ切りは簡単です。geoファイルをダブルクリックしてGmshでモデルを表示させ,左ペインのメニューからModules -> Mesh -> 3Dと選べばメッシュを切ってくれます。これをAbaquasのinp形式でエクスポートします。

3. inpファイルの修正
基本的に先の計算の場合と同じです。ただ,表層と基盤層の間の境界面の取り扱いに注意が必要です。どうやらEFISPECは境界面としては,吸収境界と自由表面以外はまったく興味がないようなのでそれ以外の面を定義することはできません。

そのため,層境界であるSurface 45は丸ごと削除します。

また,自由表面はこの場合はSurface 53なのでこの面だけELSET=fsuを指定します。それ以外の面は吸収境界としてELSET=prxを一回だけ指定して重複する*ELEMENT行は削除します。

ELSET=Lineを指定しているELEMENTも全て丸ごと削除します。

ELSET=Volume1はELSET=l01に,ELSET=Volume2はELSET=l02に書き換えます。

typeを大文字のTYPEに書き換えて,面はS4R,体積はC3D8Rに修正して,体積(l01, l02),吸収境界(prx),自由表面(fsu)の順にELEMENTの定義を並べ替えます。

あとは,ヘッダを少し書きなおして,最後の行に*END PARTを追加すれば完成です。

4. 計算の実行
あとはチュートリアルと同じ条件で計算をするだけです。チュートリアルのt01に置かれているcfgファイル(設定全般),fsrファイル(地表面のレシーバーの座標),dcsファイル(ダブルカップルのソースの設定)をコピーしてきてprefixを適切に書きなおして(例えば,t03とだけ書いておく),さらに,node.txtに計算に用いる計算ノードのホスト名を以下のように列挙します。

Node01
Node03
Node05
Node06

この場合,4 nodes x 2 ways x 8 coresで最大64コアが使えるので,設定ファイルがあるディレクトリで,

mpirun -np 64 --machinefile node.txt ../bin/efispec3d_1.0_avx.exe

とすれば64コアを使って計算をしてくれます(たぶん)。計算に関する情報は計算を行なっているディレクトリにlstという拡張子のファイルとしてログが残ります。このファイルを見ると,以下のようなことがわかります。

チュートリアルでは25,344要素でしたが,私の計算例ではメッシュを表層のサイズにあわせて全部同じサイズで切ったので,240,944要素もありました。かなりアホっぽいですが,練習なのでよしとします。

計算時間は同じように64コアを使ってチュートリアルはトータルで36.5秒,細かいメッシュを使った私の計算例では250.2秒でした。要素数が10倍弱で計算時間は7倍程度なので,計算時間が短い問題を解く場合,並列化のオーバーヘッドが馬鹿にならない,ということを示しているようです。ですのでむやみにコアをたくさん使ってもあまり賢いとは言えないかもしれません。

5. 計算結果
それでも24万要素の計算が4分あまりで終わってしまうなんていうのは遠い昔を知る者としては驚くばかりです。計算を専門にしている人からすればバカみたいなレベルの話ですが,ユーザーとして使うだけの立場からすれば計算機のありがたみを強く感じます。

スナップショットをparaviewで読み込んでjpegのパラパラアニメを出力して,ImageMagicで

convert -loop 0 2-layer_y.0* 2-layer_y_anim.gif

としてgifアニメを出力したものが下のアニメです。

2-layer_y_anim.gif

左がチュートリアルをそのまま計算したもの,右が細かいメッシュだけを使った私の計算例です。y方向の速度を表示しています。先の計算に比べてメッシュが細かいのでパッと見た限りでは両者はよくあっているように見えます。
posted by MOR at 16:23Comment(0)Computer

Spectral Element Method (SEM)を使って計算してみる(その2)

とりあえず,計算できることはわかったので次はモデルを作成してみることにします。

ところが,ここで,重大なことに気がつきました。というか最初から気がついておけ,という話なのですが,メッシュ切りに使うことが指定されているCUBITがアメリカの政府関係者しか使えないというライセンスになっていました。商用版も別途でているものの,高機能なかわりに高額,ということでそれほど複雑ではないメッシュ切りのためだけに高価なソフトを購入する財政的余裕はありません。

フリーソフトでなんとかならないものか,と考えたのですが,SEMのプログラムが六面体のメッシュを前提としているため,普通にやると四面体メッシュを切ってしまうフリーのソフトではなかなか難しそうです。また,メッシュデータはABAQUSのinpファイルでなくてはならないため,inpファイルをエクスポートできるソフトでなくてはなりません。

あれこれ探した結果,Gmshであればなんとかなりそうでしたのでこれで代替することにしました。とりあえずは,SEMのコードについてきたtutorialを再現できなくては話になりません。しかし,慣れない(というか,やったこともない)メッシュ切りはなかなかにたいへんな作業で,二層構造の地盤モデルをいきなり作れる気がしません。まずは,半無限弾性体をモデル化するところから手をつけることにします。

1. Gmshの入手
Gmshのページからバイナリをもらってきます。LinuxでもMacOSでも,もちろんWindowsでも動くようです。とりあえずMacOS版をもらってきました。2019年7月1日にv.4.4.0がリリースされたようですが,私がダウンロードした時にはなかったので,v.4.3.0を使っています。

インストールは特にやることもないので,落としてきたファイルを展開して中身を適当なディレクトリにコピーすれば使えます(MacOSの場合)。チュートリアルが付いているのでそれを見ると参考になりそうです。

2. Gmshの使い方
最初からマニュアルを読む気力もなく,安直にネットに落ちている情報を拾い読みします。概念を知るには,ここの文書が役にたちました。また,実践的な情報はここを参考にさせてもらいました。六面体メッシュの切り方について具体的に書かれているのは,後者のページくらいしか見つけられませんでした。英語のページは最初から読む気力がなかったので...

GUIでモデルを作成してメッシュを切っていくこともできますが,結局,思った通りにできないので,スクリプトを直接書いてしまうのが一番手堅い,ということがわかってきました(いっぱい失敗した,ということです)。

適当なエディタで拡張子がgeoというファイルを作成します。ここにモデルの設定をどんどん書き込んでいきます。

2.1 やりたいこと
SEMのチュートリアルに載っているモデルを作成することが目的です。x, y方向に長さ58km, z (深さ)方向に34km,表層が1kmでその下は半無限弾性体です。2層構造を作れる気がしないのでまずは,深さ34kmまでの半無限弾性体を作って,これを1kmの立方体でメッシュ分割します。

基本的にはここで例が示されている。box.geoを書き換えながらモデルを作ることにします。

2.2 ポイントの作成
まず,モデルの外形を決めるポイントを決定します。直方体なので8個のポイントの座標を定義します。

Point(1) = {-29000, -29000, 0, 1000.0};
Point(2) = {29000, -29000, 0, 1000.0};
Point(3) = {29000, 29000, 0, 1000.0};
Point(4) = {-29000, 29000, 0, 1000.0};
Point(5) = {-29000, -29000, -34000.0, 1000.0};
Point(6) = {29000, -29000, -34000.0, 1000.0};
Point(7) = {29000, 29000, -34000.0, 1000.0};
Point(8) = {-29000, 29000, -34000.0, 1000.0};

左辺のカッコ内はポイントの番号,右辺のカッコ内は順にx, y, z座標で,最後の数字はそのポイント周辺で切りたいおよそのメッシュサイズです。

2.3 ラインの作成
ポイントができたらラインを構成します。直方体なので全部で12個の辺がありますから,ポイントをつないで辺を定義します。

Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 5};
Line(9) = {1, 5};
Line(10) = {2, 6};
Line(11) = {3, 7};
Line(12) = {4, 8};

左辺のカッコ内はラインの番号,右辺の番号はラインを構成するポイントの始点と終点です。ラインはベクトルなので始点から終点に向かう方向があります。これを意識しておかないとどんどん混乱してきます。

2.4 ラインループの作成
ラインをつないでループを構成します。このループを使って次のステップで面を定義します。ラインループだけでは表面は構成されません。また,ラインループは左回りに構成しなくてはなりません。ラインがベクトルだったのでその向きに注意が必要です。

Line Loop(21) = {4, 1, 2, 3};
Line Loop(22) = {8, 5, 6, 7};
Line Loop(23) = {-1, 9, 5, -10};
Line Loop(24) = {-2, 10, 6, -11};
Line Loop(25) = {-3, 11, 7, -12};
Line Loop(26) = {-4, 12, 8, -9};

左辺のカッコ内はラインループの番号,右辺の番号はラインループを構成するラインの番号です。このとき,左回りになるようにラインを並べなくてはなりません。ラインの始点と終点が数珠つなぎになるように並べます。ラインの番号の前にマイナスをつけるとラインの方向を逆転させることができます。ラインの方向を意識して適切にマイナスをつけないとラインループを構成できなくてエラーになります。

2.5 面の作成
ラインループから面を構成します。今は直方体を作るだけなのでラインループの一つ一つが面に対応します。

Plane Surface(51) = {21};
Plane Surface(52) = {22};
Plane Surface(53) = {23};
Plane Surface(54) = {24};
Plane Surface(55) = {25};
Plane Surface(56) = {26};

左辺のカッコ内は面の番号,右辺の数字はラインループの番号です。

2.6 面ループの作成
次に面をつなぎ合わせた面ループを構成します。面ループによって囲まれた内側が次のステップで体積となります。面によって囲まれた部分が閉じていれば面ループを並べる順番はなんでもOKです。

Surface Loop(101) = {51, 52, 53, 54, 55, 56};

左辺のカッコ内は面ループの番号,右辺の数字は面ループを構成する面の番号です。

2.7 体積の作成
最後に体積を作成します。これは面ループから構成します。体積から一部をくり抜いたりすることもできるようですが,とりあえず,今はそういう難しいことは考えません。

Volume(201) = {101};

左辺のカッコ内は体積の番号,右辺の番号は体積を構成する面ループの番号です。

2.8 面のtransfiniteの設定
ここでなぜtransfiniteという語を使っているのかよくわかりませんが,transfiniteの翻訳は数学でいうところの超限数です。無限の濃度(アレフゼロとか)と云うことだろう,と理解していますが,無限の理論は得意ではないのであまり突っ込まれると困ります。兎に角,ここでやりたいことはtransfiniteという処理をして,面を四角形で,体積を六面体でメッシュに切ることです。

まず,面を四角形でメッシュ切りできるように

Transfinite Surface {51};
Transfinite Surface {52};
Transfinite Surface {53};
Transfinite Surface {54};
Transfinite Surface {55};
Transfinite Surface {56};

として各面をtranfiniteします。そのうえで,recombineします。recombineというのはどうやら三角形を二つくっつけて四角形にする,ということのようです。transfiniteしなくても四角形のメッシュに切ることはできますが,長方形や正方形のような整った形にはなりません。整った形にするためにはtransfiniteをしておかねばならないようです。transfiniteしたうえで,

Recombine Surface {51, 52, 53, 54, 55, 56};

とやるとカッコ内で指定した面が四角形でメッシュ切りできるようになります。

2.9 体積のtransfiniteの設定
六面体でメッシュ切りするためには体積についてもtransfiniteをしなくてはなりません。ただ,体積の場合はrecombineを指定しないようです。ここらへんのポリシーがよくわかりません。六面体メッシュを切るためには,面がそれぞれ4辺からなっていて,かつ体積が6面からなっていなくてはならない,という制約条件があります。そのうえで,transfiniteするに際に表面を構成するポイントの番号を反時計回りに指定してから裏面を構成するポイントの番号を反時計回りに指定しなくてはなりません。表面と裏面はどこを基準にするか,ということはどうでもよいようで,向かい合う2面を決めてそれぞれの面を構成するポイントを反時計回りに指定していく,ということのようです。

Transfinite Volume{201} = {5, 6, 2, 1, 8, 7, 3, 4};

3. メッシュ切り
3.1 GmshのGUIの起動
ここまでのコマンドを書いて拡張子をgeoとして適当な名前で保存します。MacOSの場合は保存したgeoファイルをダブルクリックするとGmshが起動してモデルを表示します。すると,以下のように作成したモデルが表示されます。下の絵はマウスでモデルを少し回した状態でスクリーンショットを撮っています。
Gmsh.png

3.2 メッシュ分割
左側のメニューのMeshの左側の三角をクリックしてメニューを広げて,3Dというメニューをクリックします。すると,メッシュ切りをしてくれます。この程度のメッシュであれば一瞬で終わります。メッシュに切れると以下のようになります。
Gmsh_2.png

3.3 メッシュデータのエクスポート
メッシュデータをABAQUSのinp形式で出力します。メニューのFile -> Export...と選ぶと出力ファイル名や出力形式を選ぶ画面が表示されます。
Gmsh_3.png
ファイル名のところには拡張子として.inpを書いたうえで一番下のファイル形式を選ぶところで「Mesh - Abaqus INP (*.inp)」を選んでSaveをクリックすればinp形式のファイルが保存されます。

4. メッシュデータの修正
以上で一応,Abaquas形式のメッシュデータが用意できたことになりますが,SEMのプログラムはGmshが吐き出したinpファイルを処理することができません。SEMのプログラムはCUBITが吐き出すinpファイルを前提としているのでそれにあうように少し手作業によって修正をしなくてはなりません。

修正すべき点は以下の通りです。

4.1 *PART行の挿入
*NODEという行の直上に*PARTという文字列を書いた行を挿入します。

4.2 *NODE行の修正
*NODEとだけ書かれている行を
*NODE, NSET=ALLNODES
と修正します。

4.3 *ELEMENT行の修正
*ELEMENTで始まる行の「type=」という部分のtypeをTYPEに書き直します。

*ELEMENTで始まる行のうち,ELSET=Line1などと書かれているラインの要素を定義されている部分は要素ごと全て削除します。この例ではLine12まであるのでそれらを全て(要素の定義ごと)削除します。

*ELEMENTで始まる行のうち,ELSET=Surface51などと書かれている部分については,この例ではSurface51が地表の自由表面,Surface52からSurface56までが吸収境界にならなくてはなりません。そのため,Surface51の行は

*ELEMENT, TYPE=S4R, ELSET=fsu

に修正します。TYPEはもともとはCPS4が指定されていますが,これをS4Rに修正しておきます。また,Surface52からSurface56については,Surface52の行のみ(最初に出てくる吸収境界面)を以下のように修正し,Surface53からSuraface56についてはその行を削除します(*ELEMENT行のみを削除して要素の定義は残します)。

*ELEMENT, TYPE=S4R, ELSET=prx

ここでもTYPEはS4Rに修正しています。

*ELEMENTで始まる行のうち,ELSET=Volume201と書かれている部分については,以下のように修正します。

*ELEMENT, TYPE=C3D8R, ELSET=l01

l01は媒質の番号で*.cfgで設定しているmaterial = 1という番号と対応する番号でなくてはなりません(l01の最初の文字は小文字のエルです)。またTYPEはもともとはC3D8ですが,これをC3D8Rに修正しています。

要素の順番を入れ替えます。SEMのプログラムは体積,吸収境界,自由表面という順番に要素が並んでいるものと想定しています。Gmshはライン要素,面要素,体積要素の順に出力するので,エディタで無理やり順番を入れ替えます。要素の番号が一番左のカラムに連番でつけられていますが,番号をつけなおす必要はありません。たんに,プログラムが読み込む順番にあわせて要素の定義の順番を入れ替えればOKです。

最後に,一番最後の行に

*END PART

を追記します。この行がないとメッシュファイルを読み終わらずに,読み込み待ちの状態でプログラムが止まってしまいます。

以上で,プログラムが読み込むことを想定しているメッシュデータのファイルになります。このようにしてできたinpファイルをSEMのプログラムの設定ファイルを置いているディレクトリにコピーして,prefixで指定したファイル名.inpとリネームした上でefispec3d_1.0_avx.exeをmpirunを使って実行すればなんだかわからないですが計算をしてくれます。

ただし,チュートリアルにinpファイルがある問題と近い問題の場合はこれでうまくいくと思いますが,複雑な境界面の場合,ちゃんと動くのかどうかまったく想像もつきません。

5. 結果

5.1 スナップショット
こうやって計算した結果をparaviewでアニメーションにしてjpegのパラパラアニメを出力してImageMagicのconvertでアニメーションgifに変換しました。
左がチュートリアル,右が上のような条件で計算した結果でy成分の速度です。チュートリアルは2層構造ですが,媒質の特性を2つの層で同じものに書き換えています。
計算条件の違いはチュートリアルは1km以深が24kmメッシュ,1km以浅が800mメッシュであるのに対して,上の計算条件では全てのメッシュが1kmメッシュです。
1-layer_y_anim.gif

計算結果を見ると,似たような波動の伝播の様子ですが,右のほうは時々パルス的なおかしなノイズが中心あたりのメッシュに入っているように見えます。メッシュの切り方とソースの位置関係によって変なノイズがはいってしまうのかもしれませんが,ちょっと簡単には原因が解明できる気がしません。

5.2 時刻歴波形
fsrファイルで設定した座標においたレシーバーの記録はgplファイルに出力されています。バイナリファイルですが,どうやら単精度の浮動小数点としてただ数字が並んでいるだけのようです。以下のようなperlスクリプトで安直にテキストに変換してしまいます。

#!/usr/bin/perl
#
$len_time_step = 40; # 40 bytes for one time step
#
$infile = $ARGV[0];
$size = -s "$infile";
$time_step = $size / $len_time_step;
$time_id = 0;
open(IN, $infile);
binmode IN;
for($ii=0;$ii<$time_step;$ii++){
read(IN, $buf, $len_time_step);
($time[$ii],$dx[$ii],$dy[$ii],$dz[$ii],$vx[$ii],$vy[$ii],$vz[$ii],$ax[$ii],$ay[$ii],$az[$ii]) = unpack("f"x10, $buf);
#
print "$time[$ii]\t$dx[$ii]\t$dy[$ii]\t$dz[$ii]\t$vx[$ii]\t$vy[$ii]\t$vz[$ii]\t$ax[$ii]\t$ay[$ii]\t$az[$ii]\n";
}
close(IN);
#

適当なファイル名(たとえば,read_gpl.pl)で保存して実行可能フラグをつけて変換したいファイル名を引数にして実行すればOKです。出力は標準出力に出てくるので適当にリダイレクトしてしまいます。

chmod +x read_gpl.pl
./read_gpl.pl t02.fsr.000001.gpl > hoge.dat

これでテキストに変換されます。左カラムから順に,時刻,変位のx, y, z成分,速度のx, y, z成分,加速度のx, y, z成分です。なんだかチュートリアルと全然違う計算結果になってるような気がします。
posted by MOR at 19:07Comment(0)Computer