Robust Treatment of Degenerate Elements in Interactive Corotational - - PowerPoint PPT Presentation

robust treatment of degenerate elements in interactive
SMART_READER_LITE
LIVE PREVIEW

Robust Treatment of Degenerate Elements in Interactive Corotational - - PowerPoint PPT Presentation

Robust Treatment of Degenerate Elements in Interactive Corotational FEM Simulations O. Civit-Flores and A. Sus n UPC-BarcelonaTech June 11, 2014 O. Civit-Flores & A. Susin (UPC) DAPD June 11, 2014 1 / 36 Outline Introduction 1


slide-1
SLIDE 1

Robust Treatment of Degenerate Elements in Interactive Corotational FEM Simulations

  • O. Civit-Flores and A. Sus´

ın

UPC-BarcelonaTech

June 11, 2014

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 1 / 36

slide-2
SLIDE 2

Outline

1

Introduction

2

Corotational FEM

3

Rotation extraction

4

Degeneration-Aware Polar Decomposition

5

Results

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 2 / 36

slide-3
SLIDE 3

Interactive FEM

Interactive simulation of deformable solids using FEM Applications: Virtual reality, surgery, training... Videogames Requirements: User interaction Efficiency Robustness Realism

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 3 / 36

slide-4
SLIDE 4

Contribution

Element degeneration threatens robustness and realism: We identify issues with existing degenerate element treatment schemes We propose a new method that avoids them

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 4 / 36

slide-5
SLIDE 5

Outline

1

Introduction

2

Corotational FEM

3

Rotation extraction

4

Degeneration-Aware Polar Decomposition

5

Results

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 5 / 36

slide-6
SLIDE 6

Finite Element Method

Partition the computational domain Ω into sub-domains Ωi with N shared nodes

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 6 / 36

slide-7
SLIDE 7

Linear FEM

Tetrahedral elements: re =    r1 . . . r4    , xe =    x1 . . . x4    , f e =    f 1 . . . f 4    The elastic forces on the nodes are: f e = −Keue, ue = xe − re Properties: Constant stiffness matrix Ke Invariant to translation, but not to rotation

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 7 / 36

slide-8
SLIDE 8

Linear FEM

linearization error

From [1]

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 8 / 36

slide-9
SLIDE 9

Corotational Linear FEM (I)

Idea: Apply LFEM in a reference system local to each element R r1 r2 r3 x2 x1 x3

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 9 / 36

slide-10
SLIDE 10

Corotational Linear FEM (II)

1

Compute deformation matrix F: F = D(xe)D(re)−1, D(ve) =

  • v2 − v1

v3 − v1 v4 − v1

  • 2

Factorize into rotation and scaling: F = RS

3

Apply linear elasticity in local coordinates: f e = −ReKe(RT

e xe − re)

Properties: Geometrically non-linear Invariant to translation and rotation

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 10 / 36

slide-11
SLIDE 11

Corotational Linear FEM (II)

1

Compute deformation matrix F: F = D(xe)D(re)−1, D(ve) =

  • v2 − v1

v3 − v1 v4 − v1

  • 2

Factorize into rotation and scaling: F = RS

3

Apply linear elasticity in local coordinates: f e = −ReKe(RT

e xe − re)

Properties: Geometrically non-linear Invariant to translation and rotation

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 10 / 36

slide-12
SLIDE 12

Corotational Linear FEM (II)

1

Compute deformation matrix F: F = D(xe)D(re)−1, D(ve) =

  • v2 − v1

v3 − v1 v4 − v1

  • 2

Factorize into rotation and scaling: F = RS

3

Apply linear elasticity in local coordinates: f e = −ReKe(RT

e xe − re)

Properties: Geometrically non-linear Invariant to translation and rotation

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 10 / 36

slide-13
SLIDE 13

Corotational Linear FEM (II)

1

Compute deformation matrix F: F = D(xe)D(re)−1, D(ve) =

  • v2 − v1

v3 − v1 v4 − v1

  • 2

Factorize into rotation and scaling: F = RS

3

Apply linear elasticity in local coordinates: f e = −ReKe(RT

e xe − re)

Properties: Geometrically non-linear Invariant to translation and rotation

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 10 / 36

slide-14
SLIDE 14

Corotational Linear FEM (II)

1

Compute deformation matrix F: F = D(xe)D(re)−1, D(ve) =

  • v2 − v1

v3 − v1 v4 − v1

  • 2

Factorize into rotation and scaling: F = RS

3

Apply linear elasticity in local coordinates: f e = −ReKe(RT

e xe − re)

Properties: Geometrically non-linear Invariant to translation and rotation

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 10 / 36

slide-15
SLIDE 15

Dynamics and Quasistatics

Node positions x ∈ R3N are the DOF: Dynamics: M¨ x = f s(x) + f d(x, ˙ x) + f ext Quasistatics: f s(x) = −f ext

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 11 / 36

slide-16
SLIDE 16

Element Degeneration

Collapse with | det(F)| < ǫ or inversion with det(F) < 0 Unphysical Unavoidable with (finite) linear forces Unavoidable due to discretization Unavoidable due to user interaction det(F) < ǫ affects F = RS factorization

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 12 / 36

slide-17
SLIDE 17

Element Degeneration

Collapse with | det(F)| < ǫ or inversion with det(F) < 0 Unphysical Unavoidable with (finite) linear forces Unavoidable due to discretization Unavoidable due to user interaction det(F) < ǫ affects F = RS factorization

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 12 / 36

slide-18
SLIDE 18

Element Degeneration

Collapse with | det(F)| < ǫ or inversion with det(F) < 0 Unphysical Unavoidable with (finite) linear forces Unavoidable due to discretization Unavoidable due to user interaction det(F) < ǫ affects F = RS factorization

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 12 / 36

slide-19
SLIDE 19

Element Degeneration

Collapse with | det(F)| < ǫ or inversion with det(F) < 0 Unphysical Unavoidable with (finite) linear forces Unavoidable due to discretization Unavoidable due to user interaction det(F) < ǫ affects F = RS factorization From [2]

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 12 / 36

slide-20
SLIDE 20

Element Degeneration

Collapse with | det(F)| < ǫ or inversion with det(F) < 0 Unphysical Unavoidable with (finite) linear forces Unavoidable due to discretization Unavoidable due to user interaction det(F) < ǫ affects F = RS factorization

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 12 / 36

slide-21
SLIDE 21

Element Degeneration

Collapse with | det(F)| < ǫ or inversion with det(F) < 0 Unphysical Unavoidable with (finite) linear forces Unavoidable due to discretization Unavoidable due to user interaction det(F) < ǫ affects F = RS factorization

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 12 / 36

slide-22
SLIDE 22

Outline

1

Introduction

2

Corotational FEM

3

Rotation extraction

4

Degeneration-Aware Polar Decomposition

5

Results

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 13 / 36

slide-23
SLIDE 23

Methods

Several methods to extract R from F are possible: Polar Decomposition [1] QR Factorization [3] Hybrid PD-QR [5] Modified Singular Value Decomposition (SVD1) [2] Coherent Singular Value Decomposition (SVD2) [4] Project/Reflect [6] Degeneration-Aware Polar Decomposition [?]

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 14 / 36

slide-24
SLIDE 24

Polar Decomposition

Factorizes F = RS, where R is orthonormal and S is symmetric Best matching, minimizes F − R2

F

Fails if | det(F)| ≤ ǫ (collapsed) Reflected R with det(R) = −1 if det(F) < 0 (inverted)

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 15 / 36

slide-25
SLIDE 25

Polar Decomposition

Factorizes F = RS, where R is orthonormal and S is symmetric Best matching, minimizes F − R2

F

Fails if | det(F)| ≤ ǫ (collapsed) Reflected R with det(R) = −1 if det(F) < 0 (inverted)

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 15 / 36

slide-26
SLIDE 26

Polar Decomposition

Factorizes F = RS, where R is orthonormal and S is symmetric Best matching, minimizes F − R2

F

Fails if | det(F)| ≤ ǫ (collapsed) Reflected R with det(R) = −1 if det(F) < 0 (inverted)

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 15 / 36

slide-27
SLIDE 27

QR Factorization

Factorizes F = RE using Gram-Schmidt orthonormalization, where R is orthonormal and E is upper-triangular Fast and Robust Handles collapsed and inverted elements seamlessly Induces Anisotropy Critical point on collapse plane

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 16 / 36

slide-28
SLIDE 28

QR Factorization

Factorizes F = RE using Gram-Schmidt orthonormalization, where R is orthonormal and E is upper-triangular Fast and Robust Handles collapsed and inverted elements seamlessly Induces Anisotropy Critical point on collapse plane

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 16 / 36

slide-29
SLIDE 29

QR Factorization

Factorizes F = RE using Gram-Schmidt orthonormalization, where R is orthonormal and E is upper-triangular Fast and Robust Handles collapsed and inverted elements seamlessly Induces Anisotropy Critical point on collapse plane

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 16 / 36

slide-30
SLIDE 30

Hybrid PD-QR

Use Polar Decomposition for undegenerate elements and QR for degenerate ones below a threshold det(F) < α Inherits good properties of PD and QR... ...but also the drawbacks of QR... ...and adds discontinuity across transition

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 17 / 36

slide-31
SLIDE 31

Hybrid PD-QR

Use Polar Decomposition for undegenerate elements and QR for degenerate ones below a threshold det(F) < α Inherits good properties of PD and QR... ...but also the drawbacks of QR... ...and adds discontinuity across transition

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 17 / 36

slide-32
SLIDE 32

Hybrid PD-QR

Use Polar Decomposition for undegenerate elements and QR for degenerate ones below a threshold det(F) < α Inherits good properties of PD and QR... ...but also the drawbacks of QR... ...and adds discontinuity across transition

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 17 / 36

slide-33
SLIDE 33

Modified Singular Value Decomposition (SVD1)

SVD factorizes F = U ˆ FVT, where U and VT are orthonormal and ˆ F is diagonal, with singular values σi ≥ σi+1 ≥ 0 Handles collapsed elements Handles inverted elements negating the smallest σi ¯ R = ¯ U ¯ VT equivalent to “invertible” PD Computationally expensive Critical point at repeated singular values

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 18 / 36

slide-34
SLIDE 34

Modified Singular Value Decomposition (SVD1)

SVD factorizes F = U ˆ FVT, where U and VT are orthonormal and ˆ F is diagonal, with singular values σi ≥ σi+1 ≥ 0 Handles collapsed elements Handles inverted elements negating the smallest σi ¯ R = ¯ U ¯ VT equivalent to “invertible” PD Computationally expensive Critical point at repeated singular values

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 18 / 36

slide-35
SLIDE 35

Modified Singular Value Decomposition (SVD1)

SVD factorizes F = U ˆ FVT, where U and VT are orthonormal and ˆ F is diagonal, with singular values σi ≥ σi+1 ≥ 0 Handles collapsed elements Handles inverted elements negating the smallest σi ¯ R = ¯ U ¯ VT equivalent to “invertible” PD Computationally expensive Critical point at repeated singular values

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 18 / 36

slide-36
SLIDE 36

Coherent Singular Value Decomposition (SVD2)

Tries to improve SVD1 by enforcing temporal coherence in the inferred inversion direction (smallest σi) Better aligned forces Computationally expensive Discontinuous

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 19 / 36

slide-37
SLIDE 37

Coherent Singular Value Decomposition (SVD2)

Tries to improve SVD1 by enforcing temporal coherence in the inferred inversion direction (smallest σi) Better aligned forces Computationally expensive Discontinuous

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 19 / 36

slide-38
SLIDE 38

Coherent Singular Value Decomposition (SVD2)

Tries to improve SVD1 by enforcing temporal coherence in the inferred inversion direction (smallest σi) Better aligned forces Computationally expensive Discontinuous

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 19 / 36

slide-39
SLIDE 39

Comparison: Force Fields

Identity and PD QR SVD1 and SVD2 Project and Reflect Forcefields experienced by x1 (red) by in a [−3, 3]2 region when x2, x3 (green, blue) are fixed

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 20 / 36

slide-40
SLIDE 40

Comparison: Trajectories

DAPD QR MSVD

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 21 / 36

slide-41
SLIDE 41

Discussion

The main problems are: Discrete-time heuristic ignores trajectories: QR, SVD1, SVD2 Critical points induce counter-intuitive rotations: QR, SVD1 Discontinuity induces jitter: Hybrid PD-QR, SVD2 Our solution: Degeneration-Aware Polar Decomposition

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 22 / 36

slide-42
SLIDE 42

Discussion

The main problems are: Discrete-time heuristic ignores trajectories: QR, SVD1, SVD2 Critical points induce counter-intuitive rotations: QR, SVD1 Discontinuity induces jitter: Hybrid PD-QR, SVD2 Our solution: Degeneration-Aware Polar Decomposition

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 22 / 36

slide-43
SLIDE 43

Outline

1

Introduction

2

Corotational FEM

3

Rotation extraction

4

Degeneration-Aware Polar Decomposition

5

Results

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 23 / 36

slide-44
SLIDE 44

Idea

Use continuous time to detect collapse Compute R avoiding rotation-past-collapse

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 24 / 36

slide-45
SLIDE 45

Detecting collapse (I)

Given a timestep [t0, t1]:

1

New degeneration if det(F(t0)) > α and det(F(t1)) ≤ α

2

Compute tc ∈ [t0, t1] solving det(F(tc)) = α

3

Compute the collapse feature pair (Ac, Bc) at tc Details: In 2D only V-E case, in 3D V-F and E-E cases Degeneration threshold α > 0 improves numerical robustness

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 25 / 36

slide-46
SLIDE 46

Detecting collapse (I)

Given a timestep [t0, t1]:

1

New degeneration if det(F(t0)) > α and det(F(t1)) ≤ α

2

Compute tc ∈ [t0, t1] solving det(F(tc)) = α

3

Compute the collapse feature pair (Ac, Bc) at tc Details: In 2D only V-E case, in 3D V-F and E-E cases Degeneration threshold α > 0 improves numerical robustness

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 25 / 36

slide-47
SLIDE 47

Detecting collapse (II)

x1 x2 x3

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 26 / 36

slide-48
SLIDE 48

Computing R (I)

We compute ¯ R for degenerate elements as follows:

1

Compute degeneration direction ˆ dc from (Ac, Bc)

2

Displace nodes of feature Ac as ¯ xi = xi + λ(α, β)ˆ dc

3

Compute ¯ F from displaced nodes ¯ xi

4

Compute ¯ R using Polar Decomposition on ¯ F Details: λ(α, β) computes a displacement length that reverts element degeneration and guarantees det( ¯ F) > ǫ Parameter β controls force alignment

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 27 / 36

slide-49
SLIDE 49

Computing R (I)

We compute ¯ R for degenerate elements as follows:

1

Compute degeneration direction ˆ dc from (Ac, Bc)

2

Displace nodes of feature Ac as ¯ xi = xi + λ(α, β)ˆ dc

3

Compute ¯ F from displaced nodes ¯ xi

4

Compute ¯ R using Polar Decomposition on ¯ F Details: λ(α, β) computes a displacement length that reverts element degeneration and guarantees det( ¯ F) > ǫ Parameter β controls force alignment

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 27 / 36

slide-50
SLIDE 50

Computing R (II)

∆h ˆ dc Reflect ˆ dc Project ˆ dc xc Project (β = 1) and Reflect (β = 2). Height ∆h depends on α

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 28 / 36

slide-51
SLIDE 51

Computing R (III)

The effect of parameter β for values 1, 2 and 10

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 29 / 36

slide-52
SLIDE 52

Computing R (VI)

ˆ dc Reflect ˆ dc x(b)

c (t)

det(F) = 0 x(a)

c (t)

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 30 / 36

slide-53
SLIDE 53

Discussion

DAPD: Cached pair (Ac, Bc) ensures temporal coherency Coherently aligned recovery forces Continuous and differentiable everywhere... ...except when Ac or Bc collapse during inversion

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 31 / 36

slide-54
SLIDE 54

Outline

1

Introduction

2

Corotational FEM

3

Rotation extraction

4

Degeneration-Aware Polar Decomposition

5

Results

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 32 / 36

slide-55
SLIDE 55

Results

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 33 / 36

slide-56
SLIDE 56

Results

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 33 / 36

slide-57
SLIDE 57

Results

(a) QR, (b) SVD1, (c) DAPD β = 1, (d) DAPD β = 5, (e) Nonlinear

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 33 / 36

slide-58
SLIDE 58

Results

http://www.lsi.upc.edu/˜ocivit/videos/DCFEM_Full.ogg http://www.lsi.upc.edu/˜ocivit/videos/DCNLFEM-Draft.avi

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 33 / 36

slide-59
SLIDE 59

Conclusions

Advantages of DAPD: Increased realism Shorter recovery time Faster than SVD1 Limitations: Cannot handle initially degenerate elements Limited to triangular/tetrahedral elements Slower than QR

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 34 / 36

slide-60
SLIDE 60

References

M.M¨ uller, M.Gross Interactive virtual materials Proceedings of Graphics Interface 2004 G.Irving, J.Teran, R.Fedkiw Invertible finite elements for robust simulation of large deformation, Symposium on Computer Animation 2004 M.Nesme, Y.Payan, F.Faure Efficient, Physically Plausible Finite Elements Eurographics 2005 R.Schmedding, M.Teschner Inversion handling for stable deformable modeling The Visual Computer 2008 E.Parker, J.O’Brien Real-time deformation and fracture in a game environment Symposium on Computer Animation 2009

  • O. Civit-Flores, A. Sus´

ın Robust treatment of degenerate elements in interactive corotational FEM simulations (to appear) Computer Graphics Forum 2014

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 35 / 36

slide-61
SLIDE 61

Q&A

?

  • O. Civit-Flores & A. Susin (UPC)

DAPD June 11, 2014 36 / 36