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



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

0 0


slides_fem_cedec2015

CEDEC2015 ショートセッション 「インタラクティブ弾性体シミュレーションのための有限要素法の実装とGPUによる高速化」 発表スライド

On Github sasekazu / slides_fem_cedec2015

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

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

佐瀬 一弥(北大),江間章斗(北大),福原洸(東北大),辻田哲平(防衛大),近野敦(北大)

About us

  • 北大 情報科学研究科 知能ロボットシステム研究室
    • ヒューマノイドロボット
    • 無人航空機
    • 脳外科手術シミュレータ

脳外科手術シミュレータ

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

できたこと その1 Polygon mesh embedding

できたこと その2 Stress visualization

できたこと その3 Fracture

Contents

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

有限要素法の概要

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

取り扱う内容

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

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

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

要素:メッシュの小領域, 節点:要素頂点の集合

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

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

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

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

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

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

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

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

変形の計算

剛性方程式

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

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

もっと詳しく…

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

動的解析

運動方程式

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

減衰項

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

参考:竹内,樫山,寺田,「計算力学」.森北出版.

時間積分

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

その場合,連立一次方程式に帰着する(導出は割愛).

参考

  • 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.

応力解析

ひずみ

応力

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による高速化

手法1

  • 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.

手法2

  • 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.

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

Georgii, J., and Westermann, R. (2008). Corotated Finite Elements Made Fast and Stable. Proceedings of the 5th Workshop on Virtual Reality Interaction and Physical Simulation, 11–19.

FEMのCG分野における商用化の動向

近年では商用ソフトへの導入も盛んに行われている.

  • 2008年 Pixelux EntertainmentがFEMベースの物理エンジンDMMを開発
  • 2009年 家庭用ゲーム機 (PS3, Xbox360) においてDMMを用いたリアルタイム破壊シミュレーションが実現 [1]
  • 2012年 MayaがDMMプラグインを搭載
  • 2013年 Houdiniが独自のFEMソルバを導入

今後はGPGPUによる高速化によりリアルタイムアプリケーションへのFEMの応用がさらに加速することが予想される.

[1] E.G. Parker and J.F. O'Brien, "Real-Time Deformation and Fracture in a Game Environment," SCA2009.

1
インタラクティブ弾性体シミュレーションのための 有限要素法の実装とGPUによる高速化 佐瀬 一弥(北大),江間章斗(北大),福原洸(東北大),辻田哲平(防衛大),近野敦(北大)