slides_fem_cedec2015

インタラクティブ弾性体シミュレーションのための

有限要素法の実装とGPUによる高速化

佐瀬一弥（北大），江間章斗（北大），福原洸（東北大），辻田哲平（防衛大），近野敦（北大）

About us

北大情報科学研究科知能ロボットシステム研究室

ヒューマノイドロボット
無人航空機
脳外科手術シミュレータ

脳外科手術シミュレータ

臓器の力学挙動をできる限り正確に計算　→　有限要素法

できたことその１ Polygon mesh embedding

できたことその２ Stress visualization

できたことその３ Fracture

有限要素法の手続きの概要
Corotational Linear FEM
GPUによる高速化

有限要素法の概要

理論を学ぶには入門書
実装に必要と思われる最小限の内容をピックアップします．

取り扱う内容

解析領域の要素分割
剛性マトリクス
動的解析
応力解析

メッシュによる解析領域の設定

三次元では四面体，六面体などによる体積メッシュが必要．

要素：メッシュの小領域，　節点：要素頂点の集合

変位ベクトルによる変形の表現

FEMでは節点の物理量をひとまとめにして大きいベクトルを用いる．

変位ベクトル u，外力ベクトル f

線形FEMにおける要素レベルのつり合いの式

$\vec{f}_e = \mathbf{K}_e \vec{u}_e$

「変位ベクトル」に「剛性マトリクス」を作用させることで「外力ベクトル」が得られる．

Ke はヤング率，ポアソン比，要素頂点の初期位置によって作成する行列．

全体レベルのつり合いの式

全体剛性マトリクスはサイズを拡張した要素剛性マトリクスの総和

変形の計算

剛性方程式

$\vec{f} = \mathbf{K} \vec{u}$

外力から変位を求めるには連立方程式を解く．

$\vec{u} = \mathbf{K}^{-1} \vec{f}$

※実際には K は特異行列となるため幾何学的境界条件を付与する必要がある．
　（大抵の入門書に説明が書かれているので割愛）

もっと詳しく…

佐瀬，HTML5による有限要素法に基づいた2次元弾性体変形シミュレーション http://www.natural-science.or.jp/article/20130922210408.php
泉，酒井，実践有限要素法シミュレーション．森北出版．
竹内，樫山，寺田，計算力学．森北出版．

動的解析

運動方程式

$\mathbf{M} \vec{\ddot{u}} + \mathbf{K} \vec{u} = \vec{f}$

質量マトリクス M は，対角要素に等価質量を配置した集中質量行列を用いることが多い．

減衰項

$\mathbf{M} \vec{\ddot{u}} + \mathbf{C} \vec{\dot{u}} + \mathbf{K} \vec{u} = \vec{f}$

減衰マトリクス C は安定性向上を目的としたレイリー減衰が用いられることが多い．

$\mathbf{C} = \alpha \mathbf{M} + \beta \mathbf{K}$

参考：竹内，樫山，寺田，「計算力学」．森北出版．

時間積分

インタラクティブCGなら陰解法．20 ms など長い時間刻みに対しても安定に計算可能

その場合，連立一次方程式に帰着する（導出は割愛）．

$\mathbf{A} \vec{v}^{i+1} = \vec{b}$

参考

M. Müller, et al., "Stable real-time deformations," Proc. SCA 2002, pp. 49-54, 2002.
Gino van den Bergen(ed.), and Dirk Gregorius(ed.), Game Physics Pearls, CRC Press, 2010.

応力解析

ひずみ

$\vec{\epsilon}_e = \mathbf{B}_e \vec{u}_e$ $\vec{\epsilon}_e = [\epsilon_x, \epsilon_y, \epsilon_z, \epsilon_{xy}, \epsilon_{yz}, \epsilon_{zx}]^{T}$

応力

$\vec{\sigma}_e = \mathbf{D}_e \mathbf{\epsilon}_e = \mathbf{D}_e \mathbf{B}_e \vec{u}_e$ $\vec{\sigma}_e = [\sigma_x, \sigma_y, \sigma_z, \sigma_{xy}, \sigma_{yz}, \sigma_{zx}]^{T}$

De: ひずみ-応力マトリクス（ヤング率，ポアソン比で決まる）

Be: 変位-ひずみマトリクス（要素節点の初期位置で決まる）

主応力

応力テンソル

\[ \mathbf{S} = \left( \begin{matrix} \sigma_x & \sigma_{xy} & \sigma_{xz} \\ \sigma_{xy} & \sigma_{y} & \sigma_{yz} \\ \sigma_{xz} & \sigma_{yz} & \sigma_{z} \end{matrix} \right) \]

主応力は応力テンソルの固有値 σ1, σ2, σ3．

四面体の応力をスカラーで表したときは最大主応力などを用いることができる．

Demo 2d

http://scc.ist.hokudai.ac.jp/~sase/workshop/cedec2015/app/fem/fem.html

Corotational FEM

Linear FEM 回転による不自然な変形

(左) Linear FEM, (右) Corotational Linear FEM

微小変形理論における線形ひずみは回転が圧縮として評価される．

幾何学的非線形性：回転に起因する非線形性

非線形FEMの枠組みで解くべきだが，Corotational FEM と呼ばれる定式化によって線形FEMを拡張することである程度対応できる．

Linear FEM

\[ \vec{f}_e = \mathbf{K}_e \left( \vec{x}_e - \vec{x}_{0e} \right) \]

Corotational Linear FEM

\[ \begin{align} \vec{f}_e & = \mathbf{R}_e \mathbf{K}_e \left( \mathbf{R}_e^T \vec{x}_e - \vec{x}_{0e} \right) = \mathbf{R}_e \mathbf{K}_e \mathbf{R}_e^T \vec{x}_e - \mathbf{R}_e \mathbf{K}_e \vec{x}_{0e} = \mathbf{K}'_e \vec{x}_e - \vec{f}_{0e} \end{align} \]

M. Müller, et al., "Stable real-time deformations," SCA 2002.
M. Müller and M. Gross, "Interactive Virtual Materials," GI 2004.

回転成分の抽出

変形勾配テンソルの右極分解（SVDを用いる）により，剛体回転成分を抽出する．

G. Irving, et al., "Invertible finite elements for robust simulation of large deformation," SIGGRAPH 2004.

ヤング率

減衰

メッシュ埋め込み

Gino van den Bergen(ed.), and Dirk Gregorius(ed.), Game Physics Pearls, CRC Press, 2010.

Limitation: 大ひずみにより振動が生じる（発散はしない）

Aleka McAdams, et al., "Efficient elasticity for character skinning with contact and collisions," SIGGRAPH2011

GPUによる高速化

手法１

Reduction array を用いたマトリクスアセンブリ and 共役勾配法による連立一次方程式求解
C. Cecka, et al., "Application of assembly of finite element methods on graphics processors for real-time elastodynamics," GPU Computing Gems, 2011.

手法２

Element-by-element 法による全体剛性マトリクスを作成せずに共役勾配法を実行する方法
J. Allard, et al., "Implicit FEM Solver on GPU for Interactive Deformation Simulation," GPU Computing Gems Jade Edition, 2012.
土木学会編，計算力学の常識，丸善．

手法1 GPUによる剛性マトリクス組み立て

疎行列格納形式を用いる CSR (Complessed Sparse Row format) など
非零要素ごとに総和元インデックスリストを作成する
非零要素ごとに並列に総和を実行

手法1 連立一次方程式の解法

共役勾配法を用いる
疎行列-ベクトル積などの行列演算を cuspase, cublass で実装
反復数は 40 ほど必要
対角スケーリング前処理が有効

金田(編)，並列数値処理，コロナ社．

結果

GPU実装したのは Matrix assembly と Solve (対角スケーリング付き共役勾配法反復数40 )
対象モデル：半脳，節点数 8647, 四面体要素数 32639
Matrix assembly 3.8 ms のうち 2.4 ms が CPU-GPU 間データ転送
Solve 14.7 ms のうち 8.2 ms が CPU-GPU 間データ転送

まとめ

Corotaional Linear FEM が実用的．
GPU実装は，CPU-GPU間データ転送を最小化すれば，5-10 倍のスピードアップが見込める．（対象規模，マシンスペックによる）

デモのご案内

さらなる高速化

マルチグリッド法

階層メッシュを用いて反復法の収束を加速させる方法．

ボクセルメッシュ，六面体要素に対して適用するアプローチが多い．

C. Dick, et al., "A real-time multigrid finite hexahedra method for elasticity simulation using CUDA," Simulation Modelling Practice and Theory, 19(2), 801–816. 2011.