Deformationテスト

デフォメーションのプログラムを組んでみた．
詳しいことは後でまとめるとして，まずは絵だけ．


初期状態の四面体です．ここからてっぺんの頂点を右にずらします．
以下に有限要素法を適用していない結果と適用した結果を並べます．


有限要素法を適用していない結果．てっぺんの頂点が移動しただけです．


有限要素法を適用した結果．てっぺんの頂点が移動するにつれて，他の頂点にも力が加わり移動しています．

今回の古典力学のシミュレーションはオイラー法で計算しているので，精度はすげー悪いはずです．論文に出ているようなタイムステップを実装しているわけではないので，そのうち誤差がたまって，四面体が爆発してしまいます．次はまともなタイムステップの計算をしてみます．

追記(2009/9/6):このデフォメーションはオイラー法でしか，古典力学の計算をできないかもしれない．ルンゲクッタ法では導関数を求めなければならない．しかしデフォメーションでは次のタイムステップの加速度は次の計算を行わない限り求まらない．つまり，ルンゲクッタ法に用いるための導関数は次の計算まで求まらないが，次の計算をするためには，次の導関数が必要となり，計算ができないということになる．だから，オイラー法で逐次計算していくしかない．ということは，みんな計算誤差なんて無視しているのか．

メモ
ちなみに，2007年のシグラフで公開されたBargteilさん達の論文に書いてあるFEMをそのまま実装しても全然うまくいきません．理由は，若干の記述を間違えている箇所があるからです．まじめに読んで実装しても，頂点は変な回転をしながら移動します．この論文が参照している2003年に公開されたTeranさん達のFinite Volume Methods for the Simulation of Skeletal Muscleか，2004年に公開されたIrvingさん達のInvertible Finite Elements For Robust Simulation of Large Deformationに記述されているFEMを実装するべきです．

F=XBはXが先に決まる

先日のF=XBで，deformation gradient Fと変形後の座標を意味するXのどちらが先に与えられるのかという疑問があったが，変形後の座標Xだということに気がついた．
Xは，力が加わり，加速度，速度が頂点を移動させた結果であるため，タイムステップがひとつ前の状態で決まる．変形後，Fが求まることで，さらに四面体の各頂点に力が伝搬するような仕組みになっている．
有限要素法の意味を考えてみるとすぐにわかるはずだった．

A Finite Element Method for Animating Large Viscoplastic Flowを再度読んでみる

Irvingさんたちの論文Invertible Finite Elements For Robust Simulation of Large Deformationを読んで，Bargteilさんたちの論文A Finite Element Method for Animating Large Viscoplastic Flowで記述されている有限要素法の触りの部分がわかった．

まず四面体を3x3行列で表現する．x₀とx₁，x₂，x₃を頂点とすると行列X_ij = x _ij - x_0jのように定義する．iは1,2,3でjは次元．
X₀を変形前の形状として，B=X₀^-1を基底行列とする．
変形後をXとして，deformation gradient FはF=XBとする．
SVDでF=UF^V^Tのように分解し，応力P^をP^=2u(F^-I)+λ Tr(F^-I)Iで計算する．uとλは素材の係数とのこと．
i番目の頂点に加わる力は，g_i=UP^V^T∑_jA_jN_jで計算する．

ここでトーシロー（素人）の私がわからないのが，変形後のXがあっての話なのか，変形させるためのFがあっての話なのか．このあと読み進めてみる．

Deformation Gradientに関するメモ

Deformation Gradientという言葉に遭遇したが，調べずに放置していたので，このあたりでメモ程度残しておく．

とりあえずググって出てきたのが，wikipedia．
http://en.wikipedia.org/wiki/Finite_strain_theory

これによると，デフォメーション前後の位置を見たとき，この近隣のサンプル点との相対的な位置の変化量を示した行列のことを指しているらしい．
具体的にはどうすればいいのかな．

剛体のアニメーションの計算では，物体内部のある任意の点は，近隣の点との相対的な距離は変化しない．しかし，デフォメーションでは，相対的な距離が変化する．準備として近隣ノードとの隣接となる点を見つけておきたい．

詳しくは，リンク先に書いてある．


Deformation Gradientが1のとき，体積は変わらないというようす．

BargteilさんらのViscoplasticモデル

BargteilさんたちのA Finite Element method for Animating Large Viscoplastic Flowの3章をさっと読んでみた．
前回読んで意味がわからず途中で断念したIrvingさんたちの論文を参考に，粘塑性の物体の変形と動きをモデル化している．今のところ，この論文に関しても，ほとんど理解できていない．
Irvingらの考えも理解できていないので，この分野の知見が足りないのだろう．

まぁ，よくわからなかったが，私の知見と照らし合わせて，何が書いてあったかまとめてみる．
3.1 Finite Element Methodでは，オブジェクトを四面体にして，四面体を変形するということが記述されていた．
四面体を行列X（ラージエックス）で表記している．この行列Xは3x3行列で，3次元上のあるノードから他の3つのノードの相対的な位置を表している．このとき，行をノードのインデックスで，列を次元としている．
シミュレーションをはじめるために，basis matrix（基底行列，と訳すのかな）を計算する．この基底行列βは初期のX0の逆行列とする．
このとき，与えられた変形状態がXのとき，Deformation gradient(変形勾配，と訳すかな)Fは，F=Xβで計算できる．これは，初期の四面体X0にDeformation gradient Fという変形行列をかけることで，変形後の座標Xが求まるという見慣れた計算．論文ではF=∂x/∂u=Xβと記述されている．uの意味が不明だが，きっと変形ベクトルを意味するんだろう．
あとの計算で使うため，FはF=UF^Vのように対角化する．
ここでは，Deformation gradient Fを計算するために，1st Piola-Kirchhoff stress(何これ?と思ってググったらwikipediaに書いてあった)を計算する．
http://en.wikipedia.org/wiki/Stress_(physics)
どうやら，素材と応力を関係させるとのことらしい．その応力P^は，P^ = 2μ(F^-I)+λTr(F^-I)I，で表されている．ここでは，λとμは材質のパラメーター．Trは行列のトレース．
それで力の計算は，g=UP^V^T∑ANで計算できるとのこと．ANは面法線の重み付き面積とする．
ここまでで，材質が関係する応力を計算し，座標変換の行列Fを求めたいということがわかった．たぶん1st Piola-Kirchhoff stressを知っていれば，ここまで読む必要なかったように思えた．

次に，3.2 Plasticity Modelについて読んだが，これまたよくわからない．
かかっている力eを可塑性e_pと弾性e_eに分けられると仮定する．体積の変化によるひずみによって加速が起こるとしている．（で，何よって感じがする）
可塑性のF_Pと弾性F_eの掛け算モデルは適切であり，F=F_P・F_eで表す．このとき，det(F_P)=1とすることで，体積の変化が起こらないことを保障するとのこと．
F^* =(det(F^))^-1/3 F^を計算するために対角化されたDeformation gradientのF^を使う．このとき，det(F^*)=1とする．そして，F^_P = (F^*)^γとする．γはねばりの係数．
つまり，(det(F^))^-1/3 F^に適当な数をかけた行列を座標変換の行列としてbasis matrix βにかけてやれば，よいということだろうか．

Deformation gradient Fの行列式が1になるように，かつ，弾性と可塑性に分けて計算すれば，それなりに見えるのではないかと，思ってきた．

Plasticモデル

今回読む論文は，Bargteilさん達の
A Finite Element Method for Animating Large Viscoplastic Flow
です．
PlasticをYahoo!辞書で引いてみますと，「プラスチックの」「ビニールの」等の素材を意味する形容詞と「可塑性の」という変形の性質を意味する形容詞が出てきます．まだ論文はちゃんと読んでいないので，はっきり言えませんが，変形について論じられているため，ここでのPlasticは「可塑性の」なのだと思います．
この論文を読み切った後は，可塑性の変形を実装しようと思います．