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のとき，体積は変わらないというようす．

BCC Lattice and Cut point結果

球の関数に対して，Labelleさん達の論文Isosurface Stuffing: Fast Tetrahedral Meshes with Good Dihdral Anglesで記述されているcut pointを計算して四面体として出力する処理を施した結果です．
これがグッドアングルかどうかはあんまり気にもならんのですが，次に点のwarpingでもしてみましょうか．

BCC Lattice and Cut point失敗したけどなんとなく出た結果

まぁ失敗なんですけど，以前からシコシコプログラムしていたBCC Latticeのドロネー図を用いたFEM Mesh生成のcut pointを計算した結果です．
cut pointは，陰関数のzero-setの近似でFEM Meshのエッジを切った点です．
一つの四面体の中のどこにcut pointがあるかを調べて，cut pointと四面体の各頂点の関数値を調べて，四面体をzero-set surfaceの近似になるように分割します．

この結果はその分割でミスっているのですが，そこはプログラムの書き間違い程度だろうと予想できるので，直せば，ある程度きれいに陰関数のzero-set surfaceを近似したFEM Meshを作ることができます．

プログラムでややこしかった点は，データの構造です．
BCC Latticeをはじめに計算し，点群を保持して，そこからドロネー図を計算し，一様に四面体を計算するのですが，ここまでで点と四面体のデータをそろえてしまっているのに，さらにcut pointなるものを付加して，四面体を分割するといったことが必要になります．

四面体クラスは，点を参照しています．
この点はBCC Latticeの点でもあり，cutpointでもあります．

cut pointや四面体の分割によってオブジェクトの数が変化するので，リスト構造を用いることを考えましたが，ヘッドのポインタから尻尾に向かって潜っていく時，リストの数が増えればその分遅くなりますので，やめました．

こんな状況でも，たまにはリスト構造のコードも組む時あるのですが，今回は，配列を予想できる数だけ確保して，アクセスすることで，少しだけストレスを減らしています．
(リスト構造でストレスがかかるとか言っている時点で，javaで言うとVectorクラスしか知らない原始人みたいに思われますが，その通りなので，言い返しません．microsoftが独自に用意してくれたライブラリとかカッコよく使いこなしたいです．)

また，四面体の分割方法は各頂点を固定して8通りあり，4頂点の配置は24通りあるため，一つの四面体を192通り確かめて，分割方法を決定します．
そこで，4!のループの中に8通りのチェック関数を入れて，チェックしました．
for(int i=0;i<4;i++){
  for(int j=0;j<4;j++){
    if(i==j) continue;
    for(int k=0;k<4;k++){
      if(k==i) continue;
      if(k==j) continue;
      for(int l=0;l<4;l++){
        if(l==i) continue;
        if(l==j) continue;
        if(l==k) continue;
        四面体pointのindex入れ替え(i,j,k,l)
        8通りチェック関数(四面体);
      }
    }
  }
}
みたいな感じで．

BCC Lattice: without cut point and subdivision

BCC Latticeでcut pointもsubdivisionもしないで，陰関数(今回は球状)の内外だけで，形状を定義した例．

LabelleさんのIsosurface stuffingのcut point

例にもよってLabelleさんのIsosurface stuffingの続きです．
BCC latticeによって構成した四面体にcut-pointの位置を推定します．
cut pointはzeroset surfaceを表現するために四面体のエッジ上にある点であり，このzeroset surfaceを近似するように定義します．このとき，zerosetを定義するためには，陰関数表現による関数が必要です．
陰関数モデリングの表面のポリゴン化で有名なBloomenthal[http://www.unchainedgeometry.com/jbloom/]の要領で，cut pointを推定します．
この論文では，四面体の頂点からcut pointの推定位置を事前に定義しておくとか言ってますけど，逐次的に計算してもいいと思います．
速度の問題ですから，そんなのどっちでもいいと思います．

時間作って，次はこのあたりを実装してみたいと思います．

とりあえずBCC Lattice

まぁ取り合えず，実装してみた．
NxNのグリッド状の格子を作り，各格子に中心点をつけ，その点群からドロネー図を計算した図．
構造はred green strategyの論文と同じメッシュ構造をしています．
隣り合う格子の中心点が接続されているところもあるし，興味深い構造していますね．

Isosurface stuffingのさわりだけ読んでみた

LabelleさんのIsosurface stuffing: fast tetrahedral Meshes with Good Dihedral Anglesを少し読んだ．
手法についての概略を読んだだけだが，実装は簡単そう．
例にもよってBCC格子なのだが，すっきりとした説明でわかりやすかった．
BCCを立方体（8頂点）にその中心となる点(1点)をグリッド状に配置してドロネー図をとったような格子だという．
そんで，陰関数を用いたモデリングみたいに内外を定義して，それに格子をフィッティングさせるという手法．
ここまでは，前回読んだA crystalline, red green strategy for meshing highly deformable objects with tetrahedraとあまり変わらないようだ．これ以降は読んでいないので，はっきりとした違いはつかめないままかもしれないが，八分木を使わないってことが違いなのかな．それってそんなに大きな違いなのかよく知らないけど，今回のは説明が十分そうなので実装はできそう．
明日か明後日，BCC格子くらいは実装してみたい．

どうやら予想通り

先週から読んでいた論文をだいたい読み終わりました．A crystalline, red green strategy for meshing highly deformable objects with tetrahedra
陰関数で定義したモデルで内外判定を行って，モデルの内側に面心立方みたいな接続の仕方をしている四面体を敷き詰めるということをやっています．面心立方みたいな，というのは，近隣の格子の中心にあるノード同士をつなげて新しい格子としている構造のこと．
もっとも内側はでかい格子で，表面になるほど細かい格子で構成し，形状を表現します．
形状の変形は，格子を構成しているエッジにバネみたいなモデルの動きをさせることで表現します．
モデルの表面に行くほど細かい格子を構成するという手法については他の手法を参照しているので，また孫引きして読まなきゃ実装できんのですが，ここらで一旦こいつからは手を引いて，次の論文に行こうと思います．
次はこれ．
Isosurface stuffing; fast tetrahedral meshes with good dihedral angles
Labelleさんたちの手法について読んでみようと思います．
テーマは引き続きBCC Latticeです．

昨日の続き，なんとなくわかってきた，かな？

今日は仕事を早く上がれたのでよかった．
帰ってからは昨日に引き続き，Neil MolinoさんたちのA crystalline, red green strategy for meshing highly deformable objects with tetrahedraを一気に読み込もうと意気込んだが，またもや少ししか進まなかった．
といっても，昨日より元気だったので，ずいぶん理解がマシになったと思う．
グリッド上にBBCを構成し，構成された四面体ってのを各辺の中心点で分割して，入力の面を近似する．このとき，内外の判定が必要なので，陰関数を用いたモデリングによって内外を判定する．そして，そのエッジをバネとして扱うことで，中身があるような変形ができる．ということか？な？
まだよくわかってないな．
ということで，また明日．

よーわからんものにいきなり当たったな

昨日続きでこの論文を読んでいたんだけど，仕事が遅かったので，今日はあまり進まず．
Neil MolinoさんたちのA crystalline, red green strategy for meshing highly deformable objects with tetrahedra
BCCってのは，Body-centered cubicの頭文字をとったもので，立方体の中心にノードがあるような構造の格子を言うみたい．
で？てって感じだが，有限要素法とかに使えるらしい．
しかしそれ以上がまだ理解できない．
つまりは，この論文の価値ある部分が全く理解できない．
なんだというんだこれは．
昨日は楽勝だぜ，とか思ってたのに．

どうやってかしらんが，BBCの四面体を分割していくようで，解像度に合わせて，分割を変えていけばいいじゃない，というような感じのことは読み取れるのだが，なんとも曖昧な状態．
これは読んだうちにはいらんな．
また明日，読むわ．
今週中に実装と実験はしたいな．とその前に，実装の価値があるかを見極めたい．

手始めに

有限要素法とやらが気になってきたので，最新のしぐらふ論文を漁ってみることにした．
で見つけたのが，これ．
Chris Wojtanさん達のFast Viscoelastic Behavior with Thin Features
まぁ，ざっというと，モデルを四面体で再メッシュ化して，有限要素法で粘弾性の変形をシミュレーションする，ということ．
実装も簡単そうだから，やってみようと思う．
手法のほとんどは他の論文に書かれた手法を実装することになりそうなので，それらをまずは実装しなきゃならない．
四面体でのメッシュ化の方法は
A crystalline, red green strategy for meshing highly deformable objects with tetrahedra
Isosurface stuffing; fast tetrahedral meshes with good dihedral angles
自己衝突は，
Robust Treatment of collisions, contact and friction for cloth animation
Elasticモデルは，
Volume conserving finite element simulations of deformable models
Plasticモデルは，
A finite element method for animating large viscoplastic flow
を読めば十分そうな感じ．
なんとも簡単な．

論文読んでいる時点でいくつかできない操作があることに気付いたので，実装して本当にできないのか確認してみようと思う．