in Computer Graphics Jan Bender Miles Macklin Matthias Mller - PowerPoint PPT Presentation

T6: Position-Based Simulation Methods in Computer Graphics Jan Bender Miles Macklin Matthias Müller

Jan Bender • Organizer • Professor at the Visual Computing Institute at Aachen University • Research topics – Rigid bodies, deformable solids, fluids – Collision detection, fracture, real-time visualization – Position based methods • Maintains open source PBD code base – github.com/InteractiveComputerGraphics/PositionBasedDynamics

Miles Macklin • Principal engineer at NVIDIA • Inventor and author of FLEX – Unified, particle based, position based solver, GPU accelerated – UE4 integration – developer.nvidia.com/flex • Research – Position based fluids – Inventor of XPBD, making PBD truly physical with a simple trick!

Matthias Müller • Leader of physics research group at NVIDIA • Co-initiator of PBD (with Thomas Jakobsen) • Co-founder of NovodeX which became physics group at NVIDIA • Research – Co-rotational FEM, SPH – Position based methods: cloth, soft bodies, shape matching, oriented particles, air meshes • www.matthiasmueller.info

Tutorial Outline • Matthias – Motivation, Basic Idea – The solver – Constraint examples for solids – Solver accelerations • Miles – Fluids – XPBD – Continuous materials – Rigid bodies

Motivation

Physical Simulations • Well studied problem in the computational sciences (since 1940s) • Complement / replace real experiments • Extreme conditions, spatial scale, time scale • Accuracy most important factor • Low accuracy – useless result

Computer Graphics • Early 1980s • Adopted methods: FEM, SPH, grid based fluids, .. • Applications – Special effects in movies and commercials – Computer games – VR • Requirements – Speed, stability, controllability – Only visual plausibility • New methods needed: e.g. PBD

Funhouse

Traditional Methods • Typically force based • Explicit integration – Simple and fast – Only conditionally stable (bad for real time apps) • Implicit integration – Expensive (multiple linearizations and solves per time step) – Numerical damping

Basic Idea

Force Based Update penetration forces velocities causes forces change velocities change positions • Reaction lag • Small spring stiffness → squashy system • Large spring stiffness → stiff system, overshooting

Position Based Update penetration move objects so that update velocities! detection only they do not penetrate • Controlled position change • Only as much as needed → no overshooting • Velocity update needed to get 2 nd order system!

Position Based Integration x 𝑜 , v 𝑜 , p, u ∈ ℝ 3𝑂 init x 0 , v 0 loop ← v 𝑜 + ∆𝑢 ∙ f 𝑓𝑦𝑢 ( x 𝑜 ) v 𝑜 velocity update ← x 𝑜 + ∆𝑢 ∙ v 𝑜 p prediction ← modify p x 𝑜+1 position correction ← (x 𝑜+1 − x 𝑜 )/∆𝑢 u velocity update ← modify u v 𝑜+1 velocity correction end loop

Position Correction • Example: Particle on circle prediction new velocity correction

Velocity Correction • External forces: v 𝑜+1 = u + ∆𝑢 g 𝑛 • Internal damping corrected • Friction velocity prediction • Restitution friction restitution collision correction

Distance Constraint 𝑥 1 x 1 − x 2 ∆x 2 ∆x 1 = − x 1 − x 2 − 𝑚 0 𝑥 1 + 𝑥 2 x 1 − x 2 ∆x 1 𝑚 0 𝑛 2 𝑥 2 x 1 − x 2 ∆x 2 = + x 1 − x 2 − 𝑚 0 𝑥 1 + 𝑥 2 x 1 − x 2 𝑛 1 𝑥 𝑗 = 1 • Conservation of momentum 𝑛 𝑗 • Stiffness: scale corrections by 𝑙 ∈ 0,1 – Easy to tune – Effect dependent on time step size and iteration count – Fixed! See XPBD

General Internal Constraint • Define constraint via scalar function: 𝐷 𝑡𝑢𝑠𝑓𝑢𝑑ℎ x 1 , x 2 = x 1 − x 2 − 𝑚 0 𝐷 𝑤𝑝𝑚𝑣𝑛𝑓 x 1 , x 2 , x 3 , x 4 = x 2 − x 1 × x 3 − x 1 ∙ x 4 − x 1 − 6𝑤 0 • Find configuration for which 𝐷 = 0 • Search along 𝛼𝐷 𝐷 = 0 𝛼𝐷 rigid body modes

Constraint Projection 𝐷 x + ∆x = 0 • Linearization (equal for distance constraint) 𝐷 x + ∆x ≈ 𝐷 x + 𝛼𝐷 x 𝑈 ∆x = 0 • Correction vectors ∆x = 𝜇 M −1 𝛼𝐷 x ∆x = 𝜇 𝛼𝐷 x 𝐷 x 𝐷 x λ = − λ = − 𝛼𝐷 x 𝑈 M −1 𝛼𝐷 x 𝛼𝐷 x 𝑈 𝛼𝐷 x M = 𝑒𝑗𝑏𝑕 𝑛 1 , 𝑛 2 , . . , 𝑛 𝑜

The Solver

Constraint Solver • Gauss-Seidel – Iterate through all constraints and apply projection – Perform multiple iterations – Simple to implement • Modified Jacobi – Process all constraints in parallel – Accumulate corrections – After each iteration, average corrections [Bridson et al., 2002] • Both known for slow convergence

Global Solver [Goldenthal et al., 2007] • Constraint vector 𝛼𝐷 1 x 𝑈 𝐷 1 x 𝜇 1 C x = ⋯ 𝛼C x = ⋯ ⋯ λ = 𝛼𝐷 𝑁 x 𝑈 𝐷 𝑁 x 𝜇 𝑁 𝐷 x ∆x = M −1 𝛼𝐷 x 𝜇 λ = − 𝛼𝐷 x 𝑈 M −1 𝛼𝐷 x 𝛼𝐷 x M −1 𝛼C x 𝑈 λ = −C x ∆x = M −1 𝛼C x 𝑈 λ

Global vs. Gauss-Seidel • Gradients fixed 𝛼𝐷 1 𝛼𝐷 2 • Linear solution ≠ true solution 𝑚 1 𝑚 2 • Multiple Newton steps necessary • Current gradients at each constraint projection 𝑚 1 𝑚 2 • Solver converges to the true solution

Other Speedup Tricks • Use as smoother in a multi-grid method • Long range distance constraints (LRA) • Hierarchy of meshes • Shape matching → more details later

Powerful Gauss-Seidel • Can handle inequality constraints trivially (LCPs, QPs)! – Fluids: separating boundary conditions [Chentanez at al., 2012] – Rigid bodies: LCP solver [Tonge et al., 2012] – Deformable objects: Long range attachments [Kim et al., 2012] • Works on non-linear problem directly • Handles under and over-constrained problems • GS + PBD: garbage in, simulation out (almost  ) • Fine grained interleaved solver trivial • Easy to implement and parallelize

Constraint Examples

Bending 𝐨 2 𝐲 2 𝐨 1 𝐲 1 𝐲 4 𝐲 3 𝐲 2 − 𝐲 1 × 𝐲 3 − 𝐲 1 𝐲 2 − 𝐲 1 × 𝐲 4 − 𝐲 1 𝐷 𝑐𝑓𝑜𝑒𝑗𝑜𝑕 (𝐲 1 , 𝐲 2 , 𝐲 𝟒 , 𝐲 𝟓 ) = 𝑏𝑑𝑝𝑡 ∙ − 𝜒 0 𝐲 2 − 𝐲 1 × 𝐲 3 − 𝐲 1 𝐲 2 − 𝐲 1 × 𝐲 4 − 𝐲 1 • More expensive than constraint 𝐷 𝑡𝑢𝑠𝑓𝑢𝑑ℎ (𝐲 3 , 𝐲 4 ) • But: Orthogonal to stretching

Stretching – Bending Independence bending resistence stretching resistance

Triangle Collision 𝐨 𝐲 2 𝐨 𝐫 𝐲 1 𝐲 2 𝐫 𝐲 3 𝐲 1 𝐲 3 𝐲 2 − 𝐲 1 × 𝐲 3 − 𝐲 1 𝐷 𝑑𝑝𝑚𝑚 (𝐲 1 , 𝐲 2 , 𝐲 𝟒 , 𝐲 𝟓 ) = 𝐫 − 𝐲 1 ∙ − ℎ 𝐲 2 − 𝐲 1 × 𝐲 3 − 𝐲 1

Cloth Example King of Wushu

Tetra Volume 𝐲 4 𝐲 2 𝐲 3 𝐲 1 𝐷 𝑏𝑗𝑠 x 1 , x 2 , x 3 , x 4 = 𝑒𝑓𝑢 x 2 − x 1 , x 3 − x 1 , x 4 − x 1 − 6𝑊 0

Soft Body Example

Global Volume - Balloons 𝐲 2 𝐲 3 𝐲 1 𝐲 4 𝐷 𝑐𝑏𝑚𝑚𝑝𝑝𝑜 (𝐲 1 , … , 𝐲 𝑶 ) = origin 𝑜 𝑢𝑠𝑗𝑏𝑜𝑕𝑚𝑓𝑡 𝐲 7 1 𝐲 𝑢 1 𝑗 × 𝐲 𝑢 2 ∙ 𝐲 𝑢 3 − 𝑙 𝑞𝑠𝑓𝑡𝑡𝑣𝑠𝑓 𝑊 𝐲 5 𝑗 𝑗 0 6 𝐲 6 𝑗=1

Air Meshes • Triangulate air • Prevent volume from inverting • Add one unilateral constraint per cell: 𝐷 𝑏𝑗𝑠 x 1 , x 2 , x 3 = x 2 − x 1 × x 3 − x 1 ≥ 0

Locking • Elements can invert without collisions • Solution: Mesh optimization (edge flips)

2D Boxes

Boxes Recovery

3D Air Meshes • Per tetra unilateral constraint: 𝐷 𝑏𝑗𝑠 x 1 , x 2 , x 3 = 𝑒𝑓𝑢 x 2 − x 1 , x 3 − x 1 , x 4 − x 1 ≥ 0 • Mesh optimization more expensive!

3D Air Meshes • Two cases that work well without optimization • Multi-layered clothing • Tissue collision • No large relative translations / rotations

Multi-Layered Clothing

Untangling

High Resolution Air Mesh

Tissue Collision Handling

Position Based Fluids [Macklin et al. 2013] • Particle based • Pair-wise lower distance constraints  granular behavior • Move particles in local neighborhood such that density = rest density • Density constraint 𝐷 x 1 , . . , x 𝑜 = 𝜍 𝑇𝑄𝐼 x 1 , . . , x 𝑜 − 𝜍 0

Position Based Fluids

Shape Matching • Optimally match rest with deformed shape • Only allow translation and rotation ∆x 𝑗 p 𝑗 • Global correction, no propagation needed • No mesh needed!

2d Demo

Optimal Translation • Given rest positions 𝐲 𝑗 , current positions 𝐲 𝑗 and masses 𝑛 𝑗 • Compute 𝐝 = 1 𝐝 𝑁 𝑛 𝑗 𝐲 𝑗 𝑁 = 𝑛 𝑗 𝑗 𝑗 𝐝 = 1 𝐮 = 𝐝 − 𝐝 𝑁 𝑛 𝑗 𝐲 𝑗 𝑗 𝐝 𝐮 = 𝐝 − 𝐝

Optimal Transformation • The optimal linear transformation is: −1 𝑈 𝑈 𝐁 = 𝑛 𝑗 𝐬 𝑗 𝐬 𝑗 𝑛 𝑗 𝐬 𝑗 𝐬 𝑗 r 𝑗 𝐝 𝑗 𝑗 = 𝐁 𝑠 𝐁 𝑡 𝐬 𝑗 = 𝐲 𝑗 − 𝐝 𝐬 𝑗 𝐝 𝐬 𝑗 = 𝐲 𝑗 − 𝐝