ゲームの物理

ήʔϜͷ෺ཧ NAOMASA MATSUBAYASHI ࣸਅࡱӨ OK SNS౤ߘ OK ಈը੾Γൈ͖ OK ∑
਺ࣜ஫ҙ https://github.com/Fadis/gct

ࠇͷλʔϯ

ࠇͷλʔϯ ࠇΛஔ͚ΔՕॴ

ࠇͷλʔϯ ࠇΛஔ͚ΔՕॴ ࠇ͕બͿ΂͖Օॴ

നͷλʔϯ

ࠇͷλʔϯ ന͕֯Λऔͬͯ͠·͏

നͷλʔϯ

നͷλʔϯ ࠇ͕֯ΛऔΕΔࣄ͕ ֬ఆ͢Δ

ͬͪ͜ͷํ͕ྑͦ͞͏ ঢ়ଶ ߦಈͷ݁Ռঢ়ଶ͕Ͳ͏มԽ͢Δ͔Λ ϧʔϧΛجʹ༧૝ͯ͠ ΑΓྑ͍ঢ়ଶʹભҠͦ͠͏ͳ ߦಈΛબͿ

ঢ়ଶ ߦಈͷ݁Ռঢ়ଶ͕Ͳ͏มԽ͢Δ͔Λ ϧʔϧΛجʹ༧૝ͯ͠ ΑΓྑ͍ঢ়ଶʹભҠͦ͠͏ͳ ߦಈΛબͿ ? ϓϨΠϠʔ͸ ϧʔϧΛཧղ͍ͯ͠ͳ͚Ε͹ͳΒͳ͍

ঢ়ଶ

ঢ়ଶ ༧૝

ΑΓྑ͍ঢ়ଶʹભҠͦ͠͏ͳߦಈ

ݹయྗֶʹै͏෺ମͷӡಈ͸ આ໌͕ແͯ͘΋୭΋͕݁ՌΛ༧૝Ͱ͖Δϧʔϧ

෺ମͷӡಈΛ ήʔϜͷϧʔϧͱ͢ΔࢼΈ͕ ݹ͔͘Βͳ͞Ε͖ͯͨ

m d2x dt2 = F ࣭఺ͷ࣭ྔ ࣭఺ͷҐஔ ࣭఺ʹՃΘΔ͍ΖΜͳྗ ӡಈํఔࣜ

ࣗવ௕ όωఆ਺ ͷόω n k p0 p1 p2 p3 ஄ੑମ

Fext p0 p1 p2 p3

m d2p0 dt2 = k p0 − p1 |p0 −
p1 | (|p0 − p1 | − n) + k p0 − p2 |p0 − p2 | (|p0 − p2 | − n) + Fext m d2p1 dt2 = k p1 − p0 |p1 − p0 | (|p1 − p0 | − n) + k p1 − p3 |p1 − p3 | (|p1 − p3 | − n) m d2p2 dt2 = k p2 − p0 |p2 − p0 | (|p2 − p0 | − n) + k p2 − p3 |p2 − p3 | (|p2 − p3 | − n) m d2p3 dt2 = k p3 − p1 |p3 − p1 | (|p3 − p1 | − n) + k p3 − p2 |p3 − p2 | (|p3 − p2 | − n)

m d2p0 dt2 = k p0 − p1 |p0 −
p1 | (|p0 − p1 | − n) + k p0 − p2 |p0 − p2 | (|p0 − p2 | − n) + Fext m d2p1 dt2 = k p1 − p0 |p1 − p0 | (|p1 − p0 | − n) + k p1 − p3 |p1 − p3 | (|p1 − p3 | − n) m d2p2 dt2 = k p2 − p0 |p2 − p0 | (|p2 − p0 | − n) + k p2 − p3 |p2 − p3 | (|p2 − p3 | − n) m d2p3 dt2 = k p3 − p1 |p3 − p1 | (|p3 − p1 | − n) + k p3 − p2 |p3 − p2 | (|p3 − p2 | − n) ͜ͷ෦෼͕ͪ͝Όͪ͝Ό͍ͯ͠ΔͷͰؔ਺ Ͱஔ͘ f

d2p0 dt2 = k m f0 (p0 , p1 ,
p2 , p3 , n) + Fext m d2p1 dt2 = k m f1 (p0 , p1 , p2 , p3 , n) d2p2 dt2 = k m f2 (p0 , p1 , p2 , p3 , n) d2p3 dt2 = k m f3 (p0 , p1 , p2 , p3 , n)

dv0 dt = k m f0 (p0 , p1 ,
p2 , p3 , n) + Fext m dv1 dt = k m f1 (p0 , p1 , p2 , p3 , n) dv2 dt = k m f2 (p0 , p1 , p2 , p3 , n) dv3 dt = k m f3 (p0 , p1 , p2 , p3 , n) dx0 dt = v0 dx1 dt = v1 dx2 dt = v2 dx3 dt = v3

v0 (t + Δt) = v0 (t) + k m
f0 (p0 (t), p1 (t), p2 (t), p3 (t), n)Δt + Fext m Δt v1 (t + Δt) = v1 (t) + k m f1 (p0 (t), p1 (t), p2 (t), p3 (t), n)Δt v2 (t + Δt) = v2 (t) + k m f2 (p0 (t), p1 (t), p2 (t), p3 (t), n)Δt v3 (t + Δt) = v3 (t) + k m f3 (p0 (t), p1 (t), p2 (t), p3 (t), n)Δt p0 (t + Δt) = p0 (t) + v0 (t)Δt p1 (t + Δt) = p1 (t) + v1 (t)Δt p2 (t + Δt) = p2 (t) + v2 (t)Δt p3 (t + Δt) = p3 (t) + v3 (t)Δt ΦΠϥʔ๏

v0 (t + Δt) = v0 (t) + k m
f0 (p0 (t), p1 (t), p2 (t), p3 (t), n)Δt + Fext m Δt v1 (t + Δt) = v1 (t) + k m f1 (p0 (t), p1 (t), p2 (t), p3 (t), n)Δt v2 (t + Δt) = v2 (t) + k m f2 (p0 (t), p1 (t), p2 (t), p3 (t), n)Δt v3 (t + Δt) = v3 (t) + k m f3 (p0 (t), p1 (t), p2 (t), p3 (t), n)Δt p0 (t + Δt) = p0 (t) + v0 (t)Δt p1 (t + Δt) = p1 (t) + v1 (t)Δt p2 (t + Δt) = p2 (t) + v2 (t)Δt p3 (t + Δt) = p3 (t) + v3 (t)Δt ط஌ͷ஋ ະ஌ͷ஋

v0 (t + Δt) = v0 (t) + k m
f0 (p0 (t), p1 (t), p2 (t), p3 (t), n)Δt + Fext m Δt v1 (t + Δt) = v1 (t) + k m f1 (p0 (t), p1 (t), p2 (t), p3 (t), n)Δt v2 (t + Δt) = v2 (t) + k m f2 (p0 (t), p1 (t), p2 (t), p3 (t), n)Δt v3 (t + Δt) = v3 (t) + k m f3 (p0 (t), p1 (t), p2 (t), p3 (t), n)Δt p0 (t + Δt) = p0 (t) + v0 (t)Δt p1 (t + Δt) = p1 (t) + v1 (t)Δt p2 (t + Δt) = p2 (t) + v2 (t)Δt p3 (t + Δt) = p3 (t) + v3 (t)Δt ྗ͔Β ଎౓Λ ٻΊΔ ଎౓͔Β ҐஔΛٻΊΔ

v0 (t + Δt) = v0 (t) + k m
f0 (p0 (t), p1 (t), p2 (t), p3 (t), n)Δt + Fext m Δt v1 (t + Δt) = v1 (t) + k m f1 (p0 (t), p1 (t), p2 (t), p3 (t), n)Δt v2 (t + Δt) = v2 (t) + k m f2 (p0 (t), p1 (t), p2 (t), p3 (t), n)Δt v3 (t + Δt) = v3 (t) + k m f3 (p0 (t), p1 (t), p2 (t), p3 (t), n)Δt p0 (t + Δt) = p0 (t) + v0 (t)Δt p1 (t + Δt) = p1 (t) + v1 (t)Δt p2 (t + Δt) = p2 (t) + v2 (t)Δt p3 (t + Δt) = p3 (t) + v3 (t)Δt ࣌ࠁ ͷ࣌఺ͷ ҐஔΛ࢖ͬͯ ࣌ࠁ ͷ ଎౓ΛٻΊ͍ͯΔ t t + Δt ཅղ๏ ͕େ͖͍ఔ ਖ਼͘͠ͳ͍ Ͱܭࢉ͢Δ͜ͱʹͳΔ ͦͷޡ͕ࠩ1αΠΫϧຖʹੵ΋͍ͬͯ͘ Δt p

v0 (t + Δt) = v0 (t) + k m
f0 (p0 (t + Δt), p1 (t + Δt), p2 (t + Δt), p3 (t + Δt), n)Δt + Fext m Δt v1 (t + Δt) = v1 (t) + k m f1 (p0 (t + Δt), p1 (t + Δt), p2 (t + Δt), p3 (t + Δt), n)Δt v2 (t + Δt) = v2 (t) + k m f2 (p0 (t + Δt), p1 (t + Δt), p2 (t + Δt), p3 (t + Δt), n)Δt v3 (t + Δt) = v3 (t) + k m f3 (p0 (t + Δt), p1 (t + Δt), p2 (t + Δt), p3 (t + Δt), n)Δt p0 (t + Δt) = p0 (t) + v0 (t + Δt)Δt p1 (t + Δt) = p1 (t) + v1 (t + Δt)Δt p2 (t + Δt) = p2 (t) + v2 (t + Δt)Δt p3 (t + Δt) = p3 (t) + v3 (t + Δt)Δt ࣌ࠁ ͷ࣌఺ͷҐஔΛ࢖ͬͯ ࣌ࠁ ͷ଎౓ΛٻΊ͍ͯΔ t + Δt t + Δt ӄղ๏ ͕େ͖ͯ͘΋ޡ͕ࠩੵ΋Γʹ͍͘ Δt

v0 (t + Δt) − k m f0 (p0 (t
+ Δt), p1 (t + Δt), p2 (t + Δt), p3 (t + Δt), n)Δt = v0 (t) + Fext m Δt v1 (t + Δt) − k m f1 (p0 (t + Δt), p1 (t + Δt), p2 (t + Δt), p3 (t + Δt), n)Δt = v1 (t) v2 (t + Δt) − k m f2 (p0 (t + Δt), p1 (t + Δt), p2 (t + Δt), p3 (t + Δt), n)Δt = v2 (t) v3 (t + Δt) − k m f3 (p0 (t + Δt), p1 (t + Δt), p2 (t + Δt), p3 (t + Δt), n)Δt = v3 (t) p0 (t + Δt) − v0 (t + Δt)Δt = p0 (t) p1 (t + Δt) − v1 (t + Δt)Δt = p1 (t) p2 (t + Δt) − v2 (t + Δt)Δt = p2 (t) p3 (t + Δt) − v3 (t + Δt)Δt = p3 (t) ӄղ๏ ط ஌ ͷ ஋ ະ஌ͷ஋ χϡʔτϯ๏΍ͦͷ೿ੜͷ൓෮ղ๏Ͱղ͘

v1 (t + Δt) − m f1 (p0 (t +
Δt), p1 (t + Δt), p2 (t + Δt), p3 (t + Δt), n)Δt = v1 (t) v2 (t + Δt) − k m f2 (p0 (t + Δt), p1 (t + Δt), p2 (t + Δt), p3 (t + Δt), n)Δt = v2 (t) v3 (t + Δt) − k m f3 (p0 (t + Δt), p1 (t + Δt), p2 (t + Δt), p3 (t + Δt), n)Δt = v3 (t) p0 (t + Δt) − v0 (t + Δt)Δt = p0 (t) p1 (t + Δt) − v1 (t + Δt)Δt = p1 (t) p2 (t + Δt) − v2 (t + Δt)Δt = p2 (t) p3 (t + Δt) − v3 (t + Δt)Δt = p3 (t) ط ஌ ͷ ஋ ະ஌ͷ஋ χϡʔτϯ๏΍ͦͷ೿ੜͷ൓෮ղ๏Ͱղ͘ ͜ͷܭࢉ͕ߴίετ ӄղ๏

ήʔϜͷϩδοΫ΍άϥϑΟοΫΛܭࢉͨ͠ ࢒ΓͷϦιʔεͰ ඵҎ಺ʹ ඵ෼Ҏ্ͷܭࢉΛ͍ͨ͠ 1 60 1 60

ӄղ๏Ͱղ͖͍ͨ Ͱ΋ඇઢܗͷ࿈ཱํఔࣜ͸ղ͖ͨ͘ͳ͍

Müller, M., Heidelberger, B., Hennix, M., & Ratcliff, J. (2006).
Position Based Dynamics. In Vriphys: 3rd Workshop in Virtual Realitiy, Interactions, and Physical Simulation. The Eurographics Association. https://doi.org/10.2312/PE/vriphys/vriphys06/071-080 https://matthias-research.github.io/pages/publications/posBasedDyn.pdf Position Based Dynamics

น Force Based Dynamics όωͰܨ͕͍ͬͯΔ ྡͷཻࢠ͔Βͷྗ Ϳ͔͍ͭͬͯΔน͔Β ड͚Δਨ௚߅ྗ ϓϨΠϠʔ͕͜ͷཻࢠΛ ௫ΜͰԡ͍ͯ͠Δ
·ཻͣࢠʹՃΘΔ ྗΛશͯٻΊΔ ྗͷ૯࿨Λ࢖ͬͯ଎౓Λߋ৽ ଎౓Λ࢖ͬͯ ཻࢠͷҐஔΛߋ৽ น

଎౓Λ࢖ͬͯ ཻࢠͷ ҐஔΛߋ৽ Position Based Dynamics

଎౓Λ࢖ͬͯ ཻࢠͷ ҐஔΛߋ৽ น Position Based Dynamics ෺ཧతʹ͓͔͍͠౓ ΛٻΊΔ

଎౓Λ࢖ͬͯ ཻࢠͷ ҐஔΛߋ৽ น Position Based Dynamics ෺ཧతʹ͓͔͍͠౓ ΛٻΊΔ ͓͔͍͠౓͕
খ͘͞ͳΔํ΁ ҐஔΛमਖ਼

଎౓Λ࢖ͬͯ ཻࢠͷ ҐஔΛߋ৽ น Position Based Dynamics ෺ཧతʹ͓͔͍͠౓ ΛٻΊΔ ͓͔͍͠౓͕
খ͘͞ͳΔํ΁ ҐஔΛमਖ਼ मਖ਼͞ΕͨҐஔʹ Ҡಈ͢ΔΑ͏ʹ ଎౓Λमਖ਼

p0 p1 C(p0 , p1 ) = |p0 − p1
| − n ʹӨڹ͢Δ੍໿͕ ͱͷؒͷόωҎ֎ʹ Կ΋ແ͍৔߹όω͸ࣗવ௕ʹͳΔഺ p0 p1 ͓͔͍͠౓ ͸ ҎԼͷΑ͏ʹఆٛͰ͖Δ C ࣗવ௕n ڑ཭੍໿

p0 p1 Λ0ʹ͢ΔΑ͏ͳ ͷ मਖ਼ྔ ͸ҎԼͷΑ͏ʹͳΔ C p0 Δp Δp0
= − C(p0 , p1 ) ∑ j ∇pj C(p0 , p1 ) 2 ∇p0 C(p0 , p1 ) Δp0 ࣗવ௕n

p0 p1 ͱ ͷ࣭ྔ͕ҟͳΔ৔߹ ͱ ͸όω͔Βಉ͡େ͖͞ͷ ޲͖͕ٯͷྗΛड͚ ͍ܰํ͕ΑΓେ͖͘ಈ͘ഺ p0 p1
p0 p1 Δp0 = − 1 m0 C(p0 , p1 ) ∑ j 1 mj ∇pj C(p0 , p1 ) 2 ∇p0 C(p0 , p1 ) ͍ܰ ॏ͍ ͱ ͷ࣭ྔͷٯ਺ͷ࿨ p0 p1 ͷ࣭ྔͷٯ਺ p0 Ͱमਖ਼ྔΛεέʔϧ͢ΔࣄͰ ॏཻ͍ࢠ΄Ͳಈ͖ʹ͘͘͢Δ

C(p0 , p1 ) = p0 − p1 − n
∇p0 C(p0 , p1 ) = p0 − p1 p0 − p1 ∇p1 C(p0 , p1 ) = − p0 − p1 p0 − p1 Δp0 = − 1 m0 C(p0 , p1 ) ∑ j 1 mj ∇pj C(p0 , p1 ) 2 ∇p0 C(p0 , p1 ) = − 1 m0 ( p0 − p1 − n) 1 m0 p0 − p1 p0 − p1 2 + 1 m1 − p0 − p1 p0 − p1 2 p0 − p1 p0 − p1 Λల։͢Δ C

= − 1 m0 ( p0 − p1 − n)
1 m0 + 1 m1 p0 − p1 p0 − p1 ϕΫτϧ Λͦͷ௕͞ Ͱׂ͍ͬͯΔͷͰ ͸୯Ґ௕ ୯Ґ௕ͷϕΫτϧͷ௕͞ ͸ৗʹ Ͱ ΋ৗʹ p0 − p1 p0 − p1 p0 − p1 p0 − p1 p0 − p1 p0 − p1 1 p0 − p1 p0 − p1 2 1 = − 1 m0 ( p0 − p1 − n) 1 m0 p0 − p1 p0 − p1 2 + 1 m1 − p0 − p1 p0 − p1 2 p0 − p1 p0 − p1 Δp0

= − 1 m0 ( p0 − p1 − n)
1 m0 + 1 m1 p0 − p1 p0 − p1 Δp0 Δp1 = + 1 m1 ( p0 − p1 − n) 1 m0 + 1 m1 p0 − p1 p0 − p1 p0 p1 ৳ͼͨόω͕ݩͷ௕͞ʹ໭Δ࣌ ྆୺͸ٯํ޲ʹಈ͘ͷͰͦΕ͸ͦ͏ ΋ಉ༷ ͱ ͸ූ߸͕ٯͳͷͰ ͱ ͸ํ޲͕ٯʹͳΔ Δp1 ∇p0 C(p0 , p1 ) ∇p1 C(p0 , p1 ) Δp0 Δp1

p0 Δp0 p1 Δp1 Λमਖ਼ͨ͠ޙ Λमਖ਼͢Δͱ ͷҐஔ͸ ͕࠷খʹͳΔҐஔͰ͸ͳ͘ͳΔ p0 p1
p0 C

p0 = p0 − 1 m0 ( p0 − p1
− n) 1 m0 + 1 m1 p0 − p1 p0 − p1 p1 = p1 + 1 m1 ( p0 − p1 − n) 1 m0 + 1 m1 p0 − p1 p0 − p1 p0 = p0 − 1 m0 ( p0 − p1 − n) 1 m0 + 1 m1 p0 − p1 p0 − p1 p1 = p1 + 1 m1 ( p0 − p1 − n) 1 m0 + 1 m1 p0 − p1 p0 − p1 ⋮ ͷमਖ਼ͱ ͷमਖ਼Λ ަޓʹ܁Γฦ͢ͱ ͱ ͕ ڞʹ࠷খʹͳΔΑ͏ͳ ͱ ͕ಘΒΕΔ p0 p1 C(p0 , p1 ) C(p1 , p0 ) p0 p1 ݴ͍׵͑Δͱ Ψ΢εɾβΠσϧ๏Ͱղ͘

p0 p1 p2 p3 C(p0 , p1 ) = |p0
− p1 | − n C(p0 , p2 ) = |p0 − p2 | − n C(p1 , p3 ) = |p1 − p3 | − n C(p2 , p3 ) = |p2 − p3 | − n 1ͭͷཻࢠʹෳ਺ͷ੍໿͕͋ͬͨΒ

p0 = p0 − 1 m0 ( p0 − p1
− n) 1 m0 + 1 m1 p0 − p1 p0 − p1 p1 = p1 + 1 m1 ( p0 − p1 − n) 1 m0 + 1 m1 p0 − p1 p0 − p1 p0 = p0 − 1 m0 ( p0 − p2 − n) 1 m0 + 1 m2 p0 − p2 p0 − p2 p2 = p2 + 1 m2 ( p0 − p2 − n) 1 m0 + 1 m2 p0 − p2 p0 − p2 p2 = p2 − 1 m2 ( p2 − p3 − n) 1 m2 + 1 m3 p2 − p3 p2 − p3 p3 = p3 + 1 m3 ( p2 − p3 − n) 1 m2 + 1 m3 p2 − p3 p2 − p3 p1 = p1 − 1 m1 ( p1 − p3 − n) 1 m1 + 1 m3 p1 − p3 p1 − p3 p3 = p3 + 1 m3 ( p1 − p3 − n) 1 m1 + 1 m3 p1 − p3 p1 − p3 ⋮ p0 p1 p2 p3 ੍໿Λ1ͭͮͭฒ΂ͯΨ΢εɾβΠσϧ

p0 p1 ෺ମͷॊΒ͔͞(=߶ੑ஋=όωఆ਺)͕ҰఆͰͳ͍৔߹ p2 ݻ͍ ॊΒ͔͍ p0 p1 p2 όωఆ਺k
όωఆ਺l ݻ͍όωͱॊΒ͔͍όωʹҾͬுΒΕΔ఺ ͸ ݻ͍όωͷํʹدΔ p1

p0 = p0 − k′ 1 m0 ( p0
− p1 − n) 1 m0 + 1 m1 p0 − p1 p0 − p1 p1 = p1 + l′ 1 m1 ( p0 − p1 − n) 1 m0 + 1 m1 p0 − p1 p0 − p1 ⋮ 0.0Ҏ্1.0ҎԼͰද͞ΕΔ૬ରతͳ߶ੑ஋Λ ʹֻ͚Δ 2ͭͷόω͕ಉ͡௕͞৳ͼ͍ͯΔ࣌ ΑΓݻ͍όωͷํ͕ ʹେ͖ͳӨڹΛ༩͑ΔΑ͏ʹͳΔ Δp p

ϝογϡ͔ΒཻࢠΛੜ੒ ଎౓ͱ֎ྗΛ࢖ཻͬͯࢠͷҐஔΛߋ৽ ཻࢠͷҐஔΛमਖ਼ p0 = p0 − k′ 01
1 m0 ( p0 − p1 − n) 1 m0 + 1 m1 p0 − p1 p0 − p1 p1 = p1 + k′ 01 1 m1 ( p0 − p1 − n) 1 m0 + 1 m1 p0 − p1 p0 − p1 p0 = p0 − k′ 02 1 m0 ( p0 − p2 − n) 1 m0 + 1 m2 p0 − p2 p0 − p2 p2 = p2 + k′ 02 1 m2 ( p0 − p2 − n) 1 m0 + 1 m2 p0 − p2 p0 − p2 p2 = p2 − k′ 23 1 m2 ( p2 − p3 − n) 1 m2 + 1 m3 p2 − p3 p2 − p3 p3 = p3 + k′ 23 1 m3 ( p2 − p3 − n) 1 m2 + 1 m3 p2 − p3 p2 − p3 p1 = p1 − k′ 13 1 m1 ( p1 − p3 − n) 1 m1 + 1 m3 p1 − p3 p1 − p3 p3 = p3 + k′ 13 1 m3 ( p1 − p3 − n) 1 m1 + 1 m3 p1 − p3 p1 − p3 ཻࢠͷҐஔΛ࢖ͬͯ଎౓Λमਖ਼ ཻࢠͷҐஔΛ࢖ͬͯϝογϡΛߋ৽ ϨϯμϦϯά Ұఆճ਺܁Γฦͨ͠ GBMTF USVF ϝογϡ͔Β੍໿Λੜ੒

ݸʑͷܭࢉࣜͷؒʹ ෳࡶͳσʔλґଘ͕ੜ͡Δ p0 = p0 − k′ 01 1
m0 ( p0 − p1 − n) 1 m0 + 1 m1 p0 − p1 p0 − p1 p1 = p1 + k′ 01 1 m1 ( p0 − p1 − n) 1 m0 + 1 m1 p0 − p1 p0 − p1 p0 = p0 − k′ 02 1 m0 ( p0 − p2 − n) 1 m0 + 1 m2 p0 − p2 p0 − p2 p2 = p2 + k′ 02 1 m2 ( p0 − p2 − n) 1 m0 + 1 m2 p0 − p2 p0 − p2 p2 = p2 − k′ 23 1 m2 ( p2 − p3 − n) 1 m2 + 1 m3 p2 − p3 p2 − p3 p3 = p3 + k′ 23 1 m3 ( p2 − p3 − n) 1 m2 + 1 m3 p2 − p3 p2 − p3 p1 = p1 − k′ 13 1 m1 ( p1 − p3 − n) 1 m1 + 1 m3 p1 − p3 p1 − p3 p3 = p3 + k′ 13 1 m3 ( p1 − p3 − n) 1 m1 + 1 m3 p1 − p3 p1 − p3 Ψ΢εɾβΠσϧ๏Ͱ͸ ൪໨ͷܭࢉࣜ͸ ൪໨·Ͱͷܭࢉͷ݁ՌΛ ೖྗͱͯ͠༻͍Δ n n − 1 GPUͰฒྻԽͰ͖ͳͯ͘ࠔΔ

p0 = p0 − k′ 01 1 m0 (
p0 − p1 − n) 1 m0 + 1 m1 p0 − p1 p0 − p1 p1 = p1 + k′ 01 1 m1 ( p0 − p1 − n) 1 m0 + 1 m1 p0 − p1 p0 − p1 p0 = p0 − k′ 02 1 m0 ( p0 − p2 − n) 1 m0 + 1 m2 p0 − p2 p0 − p2 p2 = p2 + k′ 02 1 m2 ( p0 − p2 − n) 1 m0 + 1 m2 p0 − p2 p0 − p2 p2 = p2 − k′ 23 1 m2 ( p2 − p3 − n) 1 m2 + 1 m3 p2 − p3 p2 − p3 p3 = p3 + k′ 23 1 m3 ( p2 − p3 − n) 1 m2 + 1 m3 p2 − p3 p2 − p3 p1 = p1 − k′ 13 1 m1 ( p1 − p3 − n) 1 m1 + 1 m3 p1 − p3 p1 − p3 p3 = p3 + k′ 13 1 m3 ( p1 − p3 − n) 1 m1 + 1 m3 p1 − p3 p1 − p3 Λॻ͖׵͑Δ p0 Λॻ͖׵͑Δ p1 Λॻ͖׵͑Δ p2 Λॻ͖׵͑Δ p3 ಉཻ͡ࢠͷҐஔΛ ॻ͖׵͑Δ͕ࣜෳ਺͋Δ ϠίϏ๏Λ࢖ͬͯ ฒྻԽ͢Δͱ ෳ਺ͷࣜͷ͏ͪ 1ͭͷ݁Ռ͔͠ ൓ө͞Εͳ͍

p0 = p0 + 1 2 −k′ 01 1
m0 ( p0 − p1 − n) 1 m0 + 1 m1 p0 − p1 p0 − p1 − k′ 02 1 m0 ( p0 − p2 − n) 1 m0 + 1 m2 p0 − p2 p0 − p2 p1 = p1 + 1 2 k′ 01 1 m1 ( p0 − p1 − n) 1 m0 + 1 m1 p0 − p1 p0 − p1 − k′ 13 1 m1 ( p1 − p3 − n) 1 m1 + 1 m3 p1 − p3 p1 − p3 p2 = p2 + 1 2 k′ 02 1 m2 ( p0 − p2 − n) 1 m0 + 1 m2 p0 − p2 p0 − p2 − k′ 23 1 m2 ( p2 − p3 − n) 1 m2 + 1 m3 p2 − p3 p2 − p3 p3 = p3 + 1 2 k′ 23 1 m3 ( p2 − p3 − n) 1 m2 + 1 m3 p2 − p3 p2 − p3 + k′ 13 1 m3 ( p1 − p3 − n) 1 m1 + 1 m3 p1 − p3 p1 − p3 ಉཻ͡ࢠͷҐஔΛ ॻ͖׵͑Δࣜͷ ͷฏۉΛऔΔ Δp ͜ͷ࿈ཱํఔࣜΛ ϠίϏ๏Ͱղ͘ ࣜຖʹฒྻͰܭࢉ Symmetric Successive Over Relaxation(SSOR)

struct particle_type { vec3 position; vec3 previous_position; vec3 velocity; vec3
local_r; vec3 normal; float w; uint rigid; uint phase; uint attached; }; ཻࢠͷҐஔ 1αΠΫϧલͷཻࢠͷҐஔ ཻࢠͷ଎౓ ཻࢠͷ๏ઢ ཻࢠͷ࣭ྔͷٯ਺ 1 m Ͳͷ෺ମͷҰ෦͔Λද͢ID ཻࢠ͕ݻఆ͞Ε͍ͯΔ͔ ͜Εͷ഑ྻΛ༻ҙ͢Δ

void main() { const resource_pair_type id = resource_pair[ push_constants.instance ];
const instance_resource_index_type inst = instance_resource_index[ id.inst ]; const primitive_resource_index_type prim = primitive_resource_index[ id.prim ]; const mesh_type mesh = mesh_pool[ prim.mesh ]; const mat4 l2w = matrix_pool[ inst.world_matrix ]; const uint vertex_id = gl_GlobalInvocationID.x; const bool in_range = mesh.unique_vertex_count > vertex_id; if( !in_range ) return; if( mesh.particle_offset == 0xFFFFFFFF ) return; const vec3 center_of_mass = ( mesh.rigid == 0xFFFFFFFF ) ? vec3( 0, 0, 0 ) : ( l2w * vec4( rigid_pool[ mesh.rigid ].center_of_mass.xyz, 1.0 ) ).xyz; const uint accessor_id = mesh.accessor; const vec4 lp0 = read_vertex( accessor_pool[ accessor_id + 1 ], vertex_id, vec4( 0.0, 0.0, 0.0, 1.0 ) ); const vec4 wp0 = l2w * lp0; const vec3 wn0 = ( mesh.topology == GCT_SHADER_TOPOLOGY_ID_POINT ) ? vec3( 0.0, 0.0, 0.0 ) : mat3(l2w) * read_vertex( accessor_pool[ accessor_id + 2 ], vertex_id, vec4( 0.0, 0.0, 0.0, 1.0 ) ).xyz; particle_pool[ mesh.particle_offset + vertex_id ].position = wp0.xyz / wp0.w; particle_pool[ mesh.particle_offset + vertex_id ].previous_position = wp0.xyz / wp0.w; particle_pool[ mesh.particle_offset + vertex_id ].velocity = vec3( 0, 0, 0 ); particle_pool[ mesh.particle_offset + vertex_id ].local_r = lp0.xyz - center_of_mass; particle_pool[ mesh.particle_offset + vertex_id ].normal = wn0; particle_pool[ mesh.particle_offset + vertex_id ].w = 1.0/0.01; particle_pool[ mesh.particle_offset + vertex_id ].rigid = mesh.rigid; particle_pool[ mesh.particle_offset + vertex_id ].phase = prim.mesh; } ϝογϡ͔ΒཻࢠΛੜ੒ 1ϓϦϛςΟϒ=1εϨουͰ࣮ߦ͢Δ

const instance_resource_index_type inst = instance_resource_index[ id.inst ]; const primitive_resource_index_type prim = primitive_resource_index[ id.prim ]; const mesh_type mesh = mesh_pool[ prim.mesh ]; const mat4 l2w = matrix_pool[ inst.world_matrix ]; const uint vertex_id = gl_GlobalInvocationID.x; const bool in_range = mesh.unique_vertex_count > vertex_id; if( !in_range ) return; if( mesh.particle_offset == 0xFFFFFFFF ) return; const vec3 center_of_mass = ( mesh.rigid == 0xFFFFFFFF ) ? vec3( 0, 0, 0 ) : ( l2w * vec4( rigid_pool[ mesh.rigid ].center_of_mass.xyz, 1.0 ) ).xyz; const uint accessor_id = mesh.accessor; const vec4 lp0 = read_vertex( accessor_pool[ accessor_id + 1 ], vertex_id, vec4( 0.0, 0.0, 0.0, 1.0 ) ); const vec4 wp0 = l2w * lp0; const vec3 wn0 = ( mesh.topology == GCT_SHADER_TOPOLOGY_ID_POINT ) ? vec3( 0.0, 0.0, 0.0 ) : mat3(l2w) * read_vertex( accessor_pool[ accessor_id + 2 ], vertex_id, vec4( 0.0, 0.0, 0.0, 1.0 ) ).xyz; particle_pool[ mesh.particle_offset + vertex_id ].position = wp0.xyz / wp0.w; particle_pool[ mesh.particle_offset + vertex_id ].previous_position = wp0.xyz / wp0.w; particle_pool[ mesh.particle_offset + vertex_id ].velocity = vec3( 0, 0, 0 ); particle_pool[ mesh.particle_offset + vertex_id ].local_r = lp0.xyz - center_of_mass; particle_pool[ mesh.particle_offset + vertex_id ].normal = wn0; particle_pool[ mesh.particle_offset + vertex_id ].w = 1.0/0.01; particle_pool[ mesh.particle_offset + vertex_id ].rigid = mesh.rigid; particle_pool[ mesh.particle_offset + vertex_id ].phase = prim.mesh; } ෺ཧγϛϡϨʔγϣϯͰ͸ িಥ͢ΔՄೳੑ͕͋Δෳ਺ͷ෺ମ͕ ಉ͡࠲ඪܥʹஔ͔Ε͍ͯΔඞཁ͕͋Δ ௖఺ͷ࠲ඪΛϫʔϧυ࠲ඪܥʹม׵͢Δ

const instance_resource_index_type inst = instance_resource_index[ id.inst ]; const primitive_resource_index_type prim = primitive_resource_index[ id.prim ]; const mesh_type mesh = mesh_pool[ prim.mesh ]; const mat4 l2w = matrix_pool[ inst.world_matrix ]; const uint vertex_id = gl_GlobalInvocationID.x; const bool in_range = mesh.unique_vertex_count > vertex_id; if( !in_range ) return; if( mesh.particle_offset == 0xFFFFFFFF ) return; const vec3 center_of_mass = ( mesh.rigid == 0xFFFFFFFF ) ? vec3( 0, 0, 0 ) : ( l2w * vec4( rigid_pool[ mesh.rigid ].center_of_mass.xyz, 1.0 ) ).xyz; const uint accessor_id = mesh.accessor; const vec4 lp0 = read_vertex( accessor_pool[ accessor_id + 1 ], vertex_id, vec4( 0.0, 0.0, 0.0, 1.0 ) ); const vec4 wp0 = l2w * lp0; const vec3 wn0 = ( mesh.topology == GCT_SHADER_TOPOLOGY_ID_POINT ) ? vec3( 0.0, 0.0, 0.0 ) : mat3(l2w) * read_vertex( accessor_pool[ accessor_id + 2 ], vertex_id, vec4( 0.0, 0.0, 0.0, 1.0 ) ).xyz; particle_pool[ mesh.particle_offset + vertex_id ].position = wp0.xyz / wp0.w; particle_pool[ mesh.particle_offset + vertex_id ].previous_position = wp0.xyz / wp0.w; particle_pool[ mesh.particle_offset + vertex_id ].velocity = vec3( 0, 0, 0 ); particle_pool[ mesh.particle_offset + vertex_id ].local_r = lp0.xyz - center_of_mass; particle_pool[ mesh.particle_offset + vertex_id ].normal = wn0; particle_pool[ mesh.particle_offset + vertex_id ].w = 1.0/0.01; particle_pool[ mesh.particle_offset + vertex_id ].rigid = mesh.rigid; particle_pool[ mesh.particle_offset + vertex_id ].phase = prim.mesh; } ௖఺ʹ๏ઢ͕͋ͬͨΒ ͦΕ΋ϫʔϧυ࠲ඪܥʹม׵͓ͯ͘͠

const instance_resource_index_type inst = instance_resource_index[ id.inst ]; const primitive_resource_index_type prim = primitive_resource_index[ id.prim ]; const mesh_type mesh = mesh_pool[ prim.mesh ]; const mat4 l2w = matrix_pool[ inst.world_matrix ]; const uint vertex_id = gl_GlobalInvocationID.x; const bool in_range = mesh.unique_vertex_count > vertex_id; if( !in_range ) return; if( mesh.particle_offset == 0xFFFFFFFF ) return; const vec3 center_of_mass = ( mesh.rigid == 0xFFFFFFFF ) ? vec3( 0, 0, 0 ) : ( l2w * vec4( rigid_pool[ mesh.rigid ].center_of_mass.xyz, 1.0 ) ).xyz; const uint accessor_id = mesh.accessor; const vec4 lp0 = read_vertex( accessor_pool[ accessor_id + 1 ], vertex_id, vec4( 0.0, 0.0, 0.0, 1.0 ) ); const vec4 wp0 = l2w * lp0; const vec3 wn0 = ( mesh.topology == GCT_SHADER_TOPOLOGY_ID_POINT ) ? vec3( 0.0, 0.0, 0.0 ) : mat3(l2w) * read_vertex( accessor_pool[ accessor_id + 2 ], vertex_id, vec4( 0.0, 0.0, 0.0, 1.0 ) ).xyz; particle_pool[ mesh.particle_offset + vertex_id ].position = wp0.xyz / wp0.w; particle_pool[ mesh.particle_offset + vertex_id ].previous_position = wp0.xyz / wp0.w; particle_pool[ mesh.particle_offset + vertex_id ].velocity = vec3( 0, 0, 0 ); particle_pool[ mesh.particle_offset + vertex_id ].local_r = lp0.xyz - center_of_mass; particle_pool[ mesh.particle_offset + vertex_id ].normal = wn0; particle_pool[ mesh.particle_offset + vertex_id ].w = 1.0/0.01; particle_pool[ mesh.particle_offset + vertex_id ].rigid = mesh.rigid; particle_pool[ mesh.particle_offset + vertex_id ].phase = prim.mesh; } ॳ଎0m/sɺ࣭ྔ0.01kgͰཻࢠͷ഑ྻʹه࿥͢Δ

struct distance_constraint_type { uint to_id; float stiffness; float natural_distance; float
reserved; }; ઀ଓઌͷཻࢠͷID όωͷ͔ͨ͞ όωͷࣗવ௕ ͜͏͍͏ߏ଄ମͷ഑ྻΛཻࢠͷ࠷େ਺ͷ32ഒͷ਺༻ҙ͢Δ ൪໨ͷཻࢠ͔Βܨ͕͍ͬͯΔόωΛ ͔Β ൪໨ͷߏ଄ମʹه࿥͢Δ i 32i 32 (i + 1) − 1

ϝογϡ͔Β੍໿Λੜ੒ uint i2, vec3 p0, vec3 p1, vec3 p2, uint
id_offset, float stiffness ) { distance_constraint_insert( i0, i1, id_offset, stiffness, distance( p0, p1 ) ); distance_constraint_insert( i1, i2, id_offset, stiffness, distance( p1, p2 ) ); distance_constraint_insert( i2, i0, id_offset, stiffness, distance( p2, p0 ) ); } void distance_constraint_insert( mat4 l2w, uint accessor_id, uint input_primitive_id, uint id_offset, float stiffness ) { const uint i0 = read_index( accessor_pool[ accessor_id + 0 ], input_primitive_id * 3u + 0u ); const uint i1 = read_index( accessor_pool[ accessor_id + 0 ], input_primitive_id * 3u + 1u ); const uint i2 = read_index( accessor_pool[ accessor_id + 0 ], input_primitive_id * 3u + 2u ); const vec4 p0 = l2w * read_vertex( accessor_pool[ accessor_id + 1 ], i0, vec4( 0.0, 0.0, 0.0, 1.0 ) ); const vec4 p1 = l2w * read_vertex( accessor_pool[ accessor_id + 1 ], i1, vec4( 0.0, 0.0, 0.0, 1.0 ) ); const vec4 p2 = l2w * read_vertex( accessor_pool[ accessor_id + 1 ], i2, vec4( 0.0, 0.0, 0.0, 1.0 ) ); distance_constraint_insert( i0, i1, i2, p0.xyz/p0.w, p1.xyz/p1.w, p2.xyz/p2.w, id_offset, stiffness ); } 1ϓϦϛςΟϒ=1εϨουͰ࣮ߦ͢Δ

uint i2, vec3 p0, vec3 p1, vec3 p2, uint id_offset,
float stiffness ) { distance_constraint_insert( i0, i1, id_offset, stiffness, distance( p0, p1 ) ); distance_constraint_insert( i1, i2, id_offset, stiffness, distance( p1, p2 ) ); distance_constraint_insert( i2, i0, id_offset, stiffness, distance( p2, p0 ) ); } void distance_constraint_insert( mat4 l2w, uint accessor_id, uint input_primitive_id, uint id_offset, float stiffness ) { const uint i0 = read_index( accessor_pool[ accessor_id + 0 ], input_primitive_id * 3u + 0u ); const uint i1 = read_index( accessor_pool[ accessor_id + 0 ], input_primitive_id * 3u + 1u ); const uint i2 = read_index( accessor_pool[ accessor_id + 0 ], input_primitive_id * 3u + 2u ); const vec4 p0 = l2w * read_vertex( accessor_pool[ accessor_id + 1 ], i0, vec4( 0.0, 0.0, 0.0, 1.0 ) ); const vec4 p1 = l2w * read_vertex( accessor_pool[ accessor_id + 1 ], i1, vec4( 0.0, 0.0, 0.0, 1.0 ) ); const vec4 p2 = l2w * read_vertex( accessor_pool[ accessor_id + 1 ], i2, vec4( 0.0, 0.0, 0.0, 1.0 ) ); distance_constraint_insert( i0, i1, i2, p0.xyz/p0.w, p1.xyz/p1.w, p2.xyz/p2.w, id_offset, stiffness ); } ௖఺഑ྻ͔Βࡾ֯ܗͷ3ͭͷ௖఺ΛऔΓग़͢

uint id_offset, float stiffness, float natural_distance ) { distance_constraint_insert_unidirectional( from_id,
to_id, id_offset, stiffness, natural_distance ); distance_constraint_insert_unidirectional( to_id, from_id, id_offset, stiffness, natural_distance ); } void distance_constraint_insert( uint i0, uint i1, uint i2, vec3 p0, vec3 p1, vec3 p2, uint id_offset, float stiffness ) { distance_constraint_insert( i0, i1, id_offset, stiffness, distance( p0, p1 ) ); distance_constraint_insert( i1, i2, id_offset, stiffness, distance( p1, p2 ) ); distance_constraint_insert( i2, i0, id_offset, stiffness, distance( p2, p0 ) ); } void distance_constraint_insert( mat4 l2w, uint accessor_id, uint input_primitive_id, uint id_offset, float stiffness 3ͭͷ௖఺Λ݁Ϳ3ͭͷลͰ

if( orig == 0 ) { distance_constraint_pool[ gcid ].stiffness =
stiffness; distance_constraint_pool[ gcid ].natural_distance = natural_distance; break; } gcid = distance_constraint_next( gcid ); } } void distance_constraint_insert( uint from_id, uint to_id, uint id_offset, float stiffness, float natural_distance ) { distance_constraint_insert_unidirectional( from_id, to_id, id_offset, stiffness, natural_distance ); distance_constraint_insert_unidirectional( to_id, from_id, id_offset, stiffness, natural_distance ); } void distance_constraint_insert( uint i0, uint i1, uint i2, vec3 p0, vec3 p1, vec3 p2, uint id_offset, float stiffness from_idͱto_idΛ݁Ϳล1ͭʹରͯ͠ from_id͔Βto_id΁ͷ੍໿ͱ to_id͔Βfrom_id΁ͷ੍໿Λผʑʹه࿥͢Δ

void distance_constraint_insert_unidirectional( uint from_id, uint to_id, uint id_offset, float stiffness,
float natural_distance ) { uint gcid = distance_constraint_begin( from_id, id_offset ); for( uint cid = 0; cid < 32; cid++ ) { const uint orig = atomicCompSwap( distance_constraint_pool[ gcid ].to_id, 0, to_id + 1 ); if( orig == to_id + 1 ) break; if( orig == 0 ) { distance_constraint_pool[ gcid ].stiffness = stiffness; distance_constraint_pool[ gcid ].natural_distance = natural_distance; break; } gcid = distance_constraint_next( gcid ); } } void distance_constraint_insert( uint from_id, 32εϩοτͷ઀ଓ৘ใͷ͏ͪ ·੍ͩ໿͕ه࿥͞Ε͍ͯͳ͍εϩοτʹ atomicԋࢉͰॻ͖ࠐΉ ॻ͖ࠐ΋͏ͱ͍ͯ͠Δͷͱ ಉ੍͡໿͕طʹه࿥͞Ε͍ͯͨΒ Կ΋ͤͣʹॻ͖ࠐΊͨࣄʹ͢Δ

const instance_resource_index_type inst = instance_resource_index[ id.inst ]; const primitive_resource_index_type prim = primitive_resource_index[ id.prim ]; const mesh_type mesh = mesh_pool[ prim.mesh ]; const mat4 l2w = matrix_pool[ inst.world_matrix ]; const uint vertex_id = gl_GlobalInvocationID.x; const bool in_range = mesh.unique_vertex_count > vertex_id; if( !in_range ) return; if( mesh.particle_offset == 0xFFFFFFFF ) return; if( particle_pool[ mesh.particle_offset + vertex_id ].attached != 0 ) return; const float delta_t = 1.0/60.0; const float g = 9.80665; const uint accessor_id = mesh.accessor; const vec3 n0 = ( mesh.topology == GCT_SHADER_TOPOLOGY_ID_POINT ) ? vec3( 0.0, 0.0, 0.0 ) : normalize( mat3(l2w) * read_vertex( accessor_pool[ accessor_id + 2 ], vertex_id, vec4( 0.0, 0.0, 0.0, 1.0 ) ).xyz ); vec3 v = particle_pool[ mesh.particle_offset + vertex_id ].velocity; v.y += g * delta_t; vec3 p = particle_pool[ mesh.particle_offset + vertex_id ].position; p += v * delta_t; particle_pool[ mesh.particle_offset + vertex_id ].position = p; particle_pool[ mesh.particle_offset + vertex_id ].normal = n0; } ଎౓ͱ֎ྗΛ࢖ཻͬͯࢠͷҐஔΛߋ৽ 1ཻࢠ=1εϨουͰ࣮ߦ͢Δ

const instance_resource_index_type inst = instance_resource_index[ id.inst ]; const primitive_resource_index_type prim = primitive_resource_index[ id.prim ]; const mesh_type mesh = mesh_pool[ prim.mesh ]; const mat4 l2w = matrix_pool[ inst.world_matrix ]; const uint vertex_id = gl_GlobalInvocationID.x; const bool in_range = mesh.unique_vertex_count > vertex_id; if( !in_range ) return; if( mesh.particle_offset == 0xFFFFFFFF ) return; if( particle_pool[ mesh.particle_offset + vertex_id ].attached != 0 ) return; const float delta_t = 1.0/60.0; const float g = 9.80665; const uint accessor_id = mesh.accessor; const vec3 n0 = ( mesh.topology == GCT_SHADER_TOPOLOGY_ID_POINT ) ? vec3( 0.0, 0.0, 0.0 ) : normalize( mat3(l2w) * read_vertex( accessor_pool[ accessor_id + 2 ], vertex_id, vec4( 0.0, 0.0, 0.0, 1.0 ) ).xyz ); vec3 v = particle_pool[ mesh.particle_offset + vertex_id ].velocity; v.y += g * delta_t; vec3 p = particle_pool[ mesh.particle_offset + vertex_id ].position; p += v * delta_t; particle_pool[ mesh.particle_offset + vertex_id ].position = p; particle_pool[ mesh.particle_offset + vertex_id ].normal = n0; } Δt = 1 60 ඵ ॏྗՃ଎౓ g = 9.80665

const instance_resource_index_type inst = instance_resource_index[ id.inst ]; const primitive_resource_index_type prim = primitive_resource_index[ id.prim ]; const mesh_type mesh = mesh_pool[ prim.mesh ]; const mat4 l2w = matrix_pool[ inst.world_matrix ]; const uint vertex_id = gl_GlobalInvocationID.x; const bool in_range = mesh.unique_vertex_count > vertex_id; if( !in_range ) return; if( mesh.particle_offset == 0xFFFFFFFF ) return; if( particle_pool[ mesh.particle_offset + vertex_id ].attached != 0 ) return; const float delta_t = 1.0/60.0; const float g = 9.80665; const uint accessor_id = mesh.accessor; const vec3 n0 = ( mesh.topology == GCT_SHADER_TOPOLOGY_ID_POINT ) ? vec3( 0.0, 0.0, 0.0 ) : normalize( mat3(l2w) * read_vertex( accessor_pool[ accessor_id + 2 ], vertex_id, vec4( 0.0, 0.0, 0.0, 1.0 ) ).xyz ); vec3 v = particle_pool[ mesh.particle_offset + vertex_id ].velocity; v.y += g * delta_t; vec3 p = particle_pool[ mesh.particle_offset + vertex_id ].position; p += v * delta_t; particle_pool[ mesh.particle_offset + vertex_id ].position = p; particle_pool[ mesh.particle_offset + vertex_id ].normal = n0; } ֎ྗ(ॏྗ)Ͱཻࢠͷ଎౓Λߋ৽

const instance_resource_index_type inst = instance_resource_index[ id.inst ]; const primitive_resource_index_type prim = primitive_resource_index[ id.prim ]; const mesh_type mesh = mesh_pool[ prim.mesh ]; const mat4 l2w = matrix_pool[ inst.world_matrix ]; const uint vertex_id = gl_GlobalInvocationID.x; const bool in_range = mesh.unique_vertex_count > vertex_id; if( !in_range ) return; if( mesh.particle_offset == 0xFFFFFFFF ) return; if( particle_pool[ mesh.particle_offset + vertex_id ].attached != 0 ) return; const float delta_t = 1.0/60.0; const float g = 9.80665; const uint accessor_id = mesh.accessor; const vec3 n0 = ( mesh.topology == GCT_SHADER_TOPOLOGY_ID_POINT ) ? vec3( 0.0, 0.0, 0.0 ) : normalize( mat3(l2w) * read_vertex( accessor_pool[ accessor_id + 2 ], vertex_id, vec4( 0.0, 0.0, 0.0, 1.0 ) ).xyz ); vec3 v = particle_pool[ mesh.particle_offset + vertex_id ].velocity; v.y += g * delta_t; vec3 p = particle_pool[ mesh.particle_offset + vertex_id ].position; p += v * delta_t; particle_pool[ mesh.particle_offset + vertex_id ].position = p; particle_pool[ mesh.particle_offset + vertex_id ].normal = n0; } ଎౓ͰཻࢠͷҐஔΛߋ৽ ͜ͷཻ࣌ࢠͷҠಈΛ๦͛Δ෺͸ߟྀ͠ͳ͍

vec4 pbd_distance_constraint_dx( uint particle_id, uint particle_offset, uint distance_constraint_offset, uint constraint_index
) { const uint iter = distance_constraint_begin( particle_id, distance_constraint_offset ); const bool is_end = distance_constraint_is_end( iter + constraint_index ); const distance_constraint_type dc = is_end ? distance_constraint_type( 0, 0.0, 0.0, 0.0 ) : distance_constraint_get( iter + constraint_index ); const vec3 p0 = particle_pool[ particle_offset + particle_id ].position.xyz; const vec3 p1 = particle_pool[ particle_offset + ( is_end ? particle_id : dc.to_id ) ].position.xyz; const float w0 = particle_pool[ particle_offset + particle_id ].w; const float w1 = particle_pool[ particle_offset + ( is_end ? particle_id : dc.to_id ) ].w; const float d = distance( p0, p1 ); const float c = d - dc.natural_distance; const vec3 grad_c = ( p0 - p1 ) / max( d, 0.01 ); const float dl = 1.0 / ( w0 + w1 ); const vec3 dx = -dc.stiffness * dl * w0 * c * grad_c; return subgroupAdd( vec4( dx, is_end ? 0.0 : 1.0 ) ); } ཻࢠͷҐஔΛमਖ਼ 1੍໿=1εϨουͰ࣮ߦ͢Δ

) { const uint iter = distance_constraint_begin( particle_id, distance_constraint_offset ); const bool is_end = distance_constraint_is_end( iter + constraint_index ); const distance_constraint_type dc = is_end ? distance_constraint_type( 0, 0.0, 0.0, 0.0 ) : distance_constraint_get( iter + constraint_index ); const vec3 p0 = particle_pool[ particle_offset + particle_id ].position.xyz; const vec3 p1 = particle_pool[ particle_offset + ( is_end ? particle_id : dc.to_id ) ].position.xyz; const float w0 = particle_pool[ particle_offset + particle_id ].w; const float w1 = particle_pool[ particle_offset + ( is_end ? particle_id : dc.to_id ) ].w; const float d = distance( p0, p1 ); const float c = d - dc.natural_distance; const vec3 grad_c = ( p0 - p1 ) / max( d, 0.01 ); const float dl = 1.0 / ( w0 + w1 ); const vec3 dx = -dc.stiffness * dl * w0 * c * grad_c; return subgroupAdd( vec4( dx, is_end ? 0.0 : 1.0 ) ); } ੍໿Λ1ͭऔΓग़͢ 1੍໿=1εϨουͰ࣮ߦ͢Δ ಉҰཻࢠΛ࢝఺ͱ͢Δ32εϨουΛಉҰSubgroupʹ͢Δ

) { const uint iter = distance_constraint_begin( particle_id, distance_constraint_offset ); const bool is_end = distance_constraint_is_end( iter + constraint_index ); const distance_constraint_type dc = is_end ? distance_constraint_type( 0, 0.0, 0.0, 0.0 ) : distance_constraint_get( iter + constraint_index ); const vec3 p0 = particle_pool[ particle_offset + particle_id ].position.xyz; const vec3 p1 = particle_pool[ particle_offset + ( is_end ? particle_id : dc.to_id ) ].position.xyz; const float w0 = particle_pool[ particle_offset + particle_id ].w; const float w1 = particle_pool[ particle_offset + ( is_end ? particle_id : dc.to_id ) ].w; const float d = distance( p0, p1 ); const float c = d - dc.natural_distance; const vec3 grad_c = ( p0 - p1 ) / max( d, 0.01 ); const float dl = 1.0 / ( w0 + w1 ); const vec3 dx = -dc.stiffness * dl * w0 * c * grad_c; return subgroupAdd( vec4( dx, is_end ? 0.0 : 1.0 ) ); } ΛٻΊΔ Δp

) { const uint iter = distance_constraint_begin( particle_id, distance_constraint_offset ); const bool is_end = distance_constraint_is_end( iter + constraint_index ); const distance_constraint_type dc = is_end ? distance_constraint_type( 0, 0.0, 0.0, 0.0 ) : distance_constraint_get( iter + constraint_index ); const vec3 p0 = particle_pool[ particle_offset + particle_id ].position.xyz; const vec3 p1 = particle_pool[ particle_offset + ( is_end ? particle_id : dc.to_id ) ].position.xyz; const float w0 = particle_pool[ particle_offset + particle_id ].w; const float w1 = particle_pool[ particle_offset + ( is_end ? particle_id : dc.to_id ) ].w; const float d = distance( p0, p1 ); const float c = d - dc.natural_distance; const vec3 grad_c = ( p0 - p1 ) / max( d, 0.01 ); const float dl = 1.0 / ( w0 + w1 ); const vec3 dx = -dc.stiffness * dl * w0 * c * grad_c; return subgroupAdd( vec4( dx, is_end ? 0.0 : 1.0 ) ); } ಉҰSubgroupͰܭࢉͨ͠࠷େ32ݸͷ ͷ૯࿨ΛٻΊΔ Δp Ұॹʹ੍໿ͷݸ਺΋ٻΊ͓͍ͯͯޙͰฏۉʹͰ͖ΔΑ͏ʹ͢Δ

const instance_resource_index_type inst = instance_resource_index[ id.inst ]; const primitive_resource_index_type prim = primitive_resource_index[ id.prim ]; const mesh_type mesh = mesh_pool[ prim.mesh ]; const uint vertex_id = gl_GlobalInvocationID.x; const bool in_range = mesh.unique_vertex_count > vertex_id; if( !in_range ) return; if( mesh.particle_offset == 0xFFFFFFFF ) return; if( particle_pool[ mesh.particle_offset + vertex_id ].attached != 0 ) return; const float delta_t = 1.0/60.0; vec3 p = particle_pool[ mesh.particle_offset + vertex_id ].position; vec3 prev_p = particle_pool[ mesh.particle_offset + vertex_id ].previous_position; vec3 v = ( p - prev_p ) / delta_t; particle_pool[ mesh.particle_offset + vertex_id ].velocity = v; particle_pool[ mesh.particle_offset + vertex_id ].previous_position = p; } ཻࢠͷҐஔΛ࢖ͬͯ଎౓Λमਖ਼ 1ཻࢠ=1εϨουͰ࣮ߦ͢Δ

const instance_resource_index_type inst = instance_resource_index[ id.inst ]; const primitive_resource_index_type prim = primitive_resource_index[ id.prim ]; const mesh_type mesh = mesh_pool[ prim.mesh ]; const uint vertex_id = gl_GlobalInvocationID.x; const bool in_range = mesh.unique_vertex_count > vertex_id; if( !in_range ) return; if( mesh.particle_offset == 0xFFFFFFFF ) return; if( particle_pool[ mesh.particle_offset + vertex_id ].attached != 0 ) return; const float delta_t = 1.0/60.0; vec3 p = particle_pool[ mesh.particle_offset + vertex_id ].position; vec3 prev_p = particle_pool[ mesh.particle_offset + vertex_id ].previous_position; vec3 v = ( p - prev_p ) / delta_t; particle_pool[ mesh.particle_offset + vertex_id ].velocity = v; particle_pool[ mesh.particle_offset + vertex_id ].previous_position = p; } ཻࢠͷҐஔΛ࢖ͬͯ଎౓Λमਖ਼ 1ཻࢠ=1εϨουͰ࣮ߦ͢Δ मਖ਼ޙͷҐஔ - લͷλΠϜεςοϓͰͷҐஔ Δt ଎౓ = લͷλΠϜεςοϓͰͷҐஔʹमਖ਼ޙͷҐஔΛ୅ೖ

͜ͷ2఺͚ͩΛγϛϡϨʔγϣϯ͔Βআ֎ͯ͠ಈ͔ͳ͘͢Δ

͸ ͷಈ͖ʹ੍໿ΛՃ͑Δ ΋ ͷಈ͖ʹ੍໿ΛՃ͑Δ িಥ

͸γϛϡϨʔγϣϯ։࢝લʹ ݟ͚ͭΔࣄ͕Ͱ͖Δ ͸1εςοϓຖʹ ݟ͚ͭͳ͚Ε͹ͳΒͳ͍ িಥ

୔ࢁ͋Δཻࢠͷத͔Β͋Δཻࢠͱিಥ͍ͯ͠Δ͔΋͠Εͳ͍ ཻࢠΛૉૣ͘ݟ͚ͭग़͢ʹ͸ ϘΫηϧʹ ཻࢠͷIDΛ ه࿥͢Δ

͋Δཻࢠ͕ه࿥͞ΕͨϘΫηϧͱͦͷपғͷϘΫηϧΛಡΉ ͦ͜ʹه࿥͞Εཻͨࢠ͸িಥ͍ͯ͠Δ͔΋͠Εͳ͍

໰୊఺ ຆͲͷϘΫηϧ͕ۭʹͳΔͨΊϝϞϦͷར༻ޮ཰͕ѱ͍ શͯͷཻࢠ͕ऩ·Δ࠷খͷάϦουΛࣄલʹݟੵ΋Δͷ͕ࠔ೉

#define spatial_hash_iterator uint spatial_hash_iterator spatial_hash_function( hash_table_type table, ivec3 voxel )
{ const spatial_hash_iterator h = uint( voxel.x * 92837111 ) ^ uint( voxel.y * 689287499 ) ^ uint( voxel.z * 283923481 ); return table.offset + ( h % table.size ); } ۭؒతϋογϡςʔϒϧ ϘΫηϧͷΠϯσοΫεΛద౰ͳϋογϡؔ਺ʹ͔͚ͯ εΧϥͷ੔਺஋ʹ͢Δ

ۭؒతϋογϡςʔϒϧ ϋογϡؔ਺ 6100 6100 6110 ⋯ ⋯ ͕ ʹ͋Δ p1
(763,12,5) p1

(763,12,5) p0 ಉ͡ϋογϡ஋ͷҐஔʹ طʹଞͷཻࢠ͕ه࿥͞Ε͍ͯͨΒ

(763,12,5) p0 1ͭྡʹه࿥ΛࢼΈΔ p1

ۭؒతϋογϡςʔϒϧ ϋογϡؔ਺ 6100 6100 6110 ⋯ ⋯ ʹ͋Δཻࢠ͕ཉ͍͠ (763,12,5) p0
ϘΫηϧͷ࠲ඪΛϋογϡؔ਺ʹ͔͚Δ ϋογϡ஋ʹରԠ͢ΔҐஔ͔Βॱʹ ه࿥͞Ε͍ͯΔ஋ΛಡΉ p1 ه࿥͞Εཻ͍ͯͨࢠ͕࣮ࡍʹͦͷϘΫηϧͷதʹ͋ͬͨΒ ͦΕΛ݁Ռͱͯ͠औΓग़͢ Կ΋ه࿥͞Ε͍ͯͳ͍ॴʹୡͨ͠Β୳ࡧऴྃ

bool spatial_hash_write_table( hash_table_type table, spatial_hash_iterator current, ivec3 key, uint value
) { const uint orig = atomicCompSwap( spatial_hash_pool[ current ].index, 0u, value + 1u ); if( orig == 0 ) { spatial_hash_pool[ current ].voxel = key; } return orig == 0; } bool spatial_hash_not_match( hash_table_type table, spatial_hash_iterator current, ivec3 voxel ) { const spatial_hash_type key_value = spatial_hash_read_table( table, current ); return key_value.index != 0 && voxel != key_value.voxel; } uint spatial_hash_next( hash_table_type table, spatial_hash_iterator current, ivec3 voxel ) { uint candidate = spatial_hash_next( table, current ); for( uint i = 0u; i != 32u; i++ ) { ۭؒతϋογϡςʔϒϧʹॻ͘ 1ཻࢠ=1εϨουͰ࣮ߦ͢Δ

return candidate; } bool spatial_hash_is_end( hash_table_type table, spatial_hash_iterator current )
{ return spatial_hash_pool[ current ].index == 0; } spatial_hash_iterator spatial_hash_find( hash_table_type table, ivec3 voxel ) { uint current = spatial_hash_function( table, voxel ); for( uint i = 0u; i != 32u; i++ ) { if( !spatial_hash_not_match( table, current, voxel ) ) break; current = spatial_hash_next( table, current, voxel ); } return current; } void spatial_hash_insert( hash_table_type table, ivec3 voxel, uint value ) { uint current = spatial_hash_find( table, voxel ); while( spatial_hash_write_table( table, current, voxel, value ) == false ) { current = spatial_hash_next( table, current ); } } ϘΫηϧʹཻࢠ͕͋Δͱͨ͠Β ͦΕ͕ه࿥͞Ε͍ͯΔഺͷҐஔΛٻΊΔ ͦͷҐஔʹ৽ཻ͍͠ࢠͷه࿥ΛࢼΈΔ طʹԿ͔ॻ͔Ε͍ͯͨΒྡΛࢼ͢

const instance_resource_index_type inst = instance_resource_index[ id.inst ]; const primitive_resource_index_type prim = primitive_resource_index[ id.prim ]; const mesh_type mesh = mesh_pool[ prim.mesh ]; const uint particle_id = gl_GlobalInvocationID.x / 27u; const uint offset_id = gl_GlobalInvocationID.x % 27u; const bool in_range = mesh.unique_vertex_count > particle_id; if( !in_range ) return; if( mesh.particle_offset == 0xFFFFFFFF ) return; if( mesh.constraint_offset == 0xFFFFFFFF ) return; if( particle_pool[ mesh.particle_offset + particle_id ].attached != 0 ) return; const vec3 p0 = particle_pool[ mesh.particle_offset + particle_id ].position; const uint t0 = particle_pool[ mesh.particle_offset + particle_id ].phase; const ivec3 center = spatial_hash_position_to_voxel( spatial_hash_config.config, p0 ); uint dc[32]; for( uint i = 0u; i < 32u; i++ ) { dc[ i ] = mesh.particle_offset + distance_constraint_pool[ mesh.distance_constraint_offset + particle_id * 32u + i ].to_id; } const ivec3 voxel = center + voxel_offset[ offset_id ]; spatial_hash_iterator iter = spatial_hash_find( spatial_hash_config.config, voxel ); for( uint i = 0u; i < 32u; ++i ) { if( spatial_hash_is_end( spatial_hash_config.config, iter ) ) break; if( spatial_hash_not_match( spatial_hash_config.config, iter, voxel ) ) break; const uint to_id = spatial_hash_get( spatial_hash_config.config, iter ); if( mesh.particle_offset + particle_id != to_id ) { bool is_dc = false; for( uint j = 0u; j != 32u; j++ ) { ۭؒతϋογϡςʔϒϧΛಡΉ 1ཻࢠ=27εϨουͰ࣮ߦ͢Δ

if( particle_pool[ mesh.particle_offset + particle_id ].attached != 0 ) return;
const vec3 p0 = particle_pool[ mesh.particle_offset + particle_id ].position; const uint t0 = particle_pool[ mesh.particle_offset + particle_id ].phase; const ivec3 center = spatial_hash_position_to_voxel( spatial_hash_config.config, p0 ); uint dc[32]; for( uint i = 0u; i < 32u; i++ ) { dc[ i ] = mesh.particle_offset + distance_constraint_pool[ mesh.distance_constraint_offset + particle_id * 32u + i ].to_id; } const ivec3 voxel = center + voxel_offset[ offset_id ]; spatial_hash_iterator iter = spatial_hash_find( spatial_hash_config.config, voxel ); for( uint i = 0u; i < 32u; ++i ) { if( spatial_hash_is_end( spatial_hash_config.config, iter ) ) break; if( spatial_hash_not_match( spatial_hash_config.config, iter, voxel ) ) break; const uint to_id = spatial_hash_get( spatial_hash_config.config, iter ); if( mesh.particle_offset + particle_id != to_id ) { bool is_dc = false; for( uint j = 0u; j != 32u; j++ ) { if( dc[ j ] == 0u ) break; if( dc[ j ] == to_id ) { is_dc = true; break; } } if( !is_dc ) { const vec3 p1 = particle_pool[ to_id ].position; const uint t1 = particle_pool[ to_id ].phase; if( t0 == t1 ) { if( distance( p0, p1 ) < spatial_hash_config.config.scale / 2 ) { constraint_insert_unidirectional( িಥ͍ͯ͠ΔՄೳੑ͕͋Δཻࢠ͕ ڑ཭੍໿ͷؔ܎ʹ͋ͬͨΒແࢹ͢Δ ͦͷཻࢠ͸ڑ཭੍໿ͷࣜʹΑͬͯద੾ʹॲཧ͞ΕΔ

break; } } if( !is_dc ) { const vec3 p1
= particle_pool[ to_id ].position; const uint t1 = particle_pool[ to_id ].phase; if( t0 == t1 ) { if( distance( p0, p1 ) < spatial_hash_config.config.scale / 2 ) { constraint_insert_unidirectional( particle_id, to_id, mesh.constraint_offset ); } } else if( distance( p0, p1 ) < spatial_hash_config.config.scale ) { constraint_insert_unidirectional( particle_id, to_id, mesh.constraint_offset ); } } } iter = spatial_hash_next(spatial_hash_config.config, iter, voxel ); } } ࣗݾিಥͩͬͨ৔߹2ͭͷཻࢠ͕10cmະຬͷڑ཭ʹ͋ͬͨΒিಥΛٙ͏ ࣗݾিಥͰͳ͔ͬͨΒ20cmະຬͷڑ཭ʹ͋ͬͨΒিಥΛٙ͏ ٙΘཻ͍͠ࢠ͸িಥ੍໿Ϧετʹ௥Ճ͢Δ

෍ͷް͞ িಥΛٙ͏൒ܘ িಥ͔ͨ͠΋͠Εͳ͍ ཻࢠͷҎલͷҐஔ িಥ͔ͨ͠΋͠Εͳ͍ ཻࢠͷݱࡏͷҐஔ ͱিಥ͔ͨ͠΋͠Εͳཻ͍ࢠ ͕ࢵͷྖҬΛ௨͍ͬͯͨΒ 2ͭͷཻࢠ͸িಥ͍ͯ͠ΔͷͰ ͸ࢵͷྖҬͷ֎·Ͱԡ͠໭͞ΕΔ
ͷ๏ઢ িಥ੍໿

vec4 pbd_collision_constraint_dx( uint particle_id, uint particle_offset, uint constraint_offset, uint constraint_index
) { const uint iter = constraint_begin( particle_id, constraint_offset ); const bool is_end = constraint_is_end( iter + constraint_index ); const uint dc = is_end ? 0u : constraint_get( iter + constraint_index ); const vec3 p0 = particle_pool[ particle_offset + particle_id ].position; const vec3 pp0 = particle_pool[ particle_offset + particle_id ].previous_position; const vec3 p1 = particle_pool[ ( is_end ? ( particle_offset + particle_id ) : dc ) ].position; vec3 n1 = particle_pool[ ( is_end ? ( particle_offset + particle_id ) : dc ) ].normal; const bool has_normal = n1 != vec3( 0.0, 0.0, 0.0 ); n1 = has_normal ? n1 : normalize( p0 - p1 ); const bool front = dot( ( pp0 - p1 ), n1 ) >= 0.0; n1 = front ? n1 : -n1; const float thickness = 0.1; const float w0 = particle_pool[ particle_offset + particle_id ].w; const float w1 = particle_pool[ ( is_end ? ( particle_offset + particle_id ) : dc ) ].w; const float d = distance( p0, p1 ); const float dx_length = -( w0*( dot( n1, p0 - p1 ) - thickness ) )/( w0 + w1 ); const vec3 dx = ( is_end ? vec3( 0, 0, 0 ) : ( dx_length > 0.0 ? n1 * dx_length : vec3( 0, 0, 0 ) ) ); return subgroupAdd( vec4( dx, ( dx == vec3( 0, 0, 0 ) ) ? 0.0 : 1.0 ) ); ཻࢠͷҐஔΛमਖ਼ 4੍໿=1εϨουͰ࣮ߦ͢Δ

vec4 pbd_collision_constraint_dx( uint particle_id, uint particle_offset, uint constraint_offset, uint constraint_index
) { const uint iter = constraint_begin( particle_id, constraint_offset ); const bool is_end = constraint_is_end( iter + constraint_index ); const uint dc = is_end ? 0u : constraint_get( iter + constraint_index ); const vec3 p0 = particle_pool[ particle_offset + particle_id ].position; const vec3 pp0 = particle_pool[ particle_offset + particle_id ].previous_position; const vec3 p1 = particle_pool[ ( is_end ? ( particle_offset + particle_id ) : dc ) ].position; vec3 n1 = particle_pool[ ( is_end ? ( particle_offset + particle_id ) : dc ) ].normal; const bool has_normal = n1 != vec3( 0.0, 0.0, 0.0 ); n1 = has_normal ? n1 : normalize( p0 - p1 ); const bool front = dot( ( pp0 - p1 ), n1 ) >= 0.0; n1 = front ? n1 : -n1; const float thickness = 0.1; const float w0 = particle_pool[ particle_offset + particle_id ].w; const float w1 = particle_pool[ ( is_end ? ( particle_offset + particle_id ) : dc ) ].w; const float d = distance( p0, p1 ); const float dx_length = -( w0*( dot( n1, p0 - p1 ) - thickness ) )/( w0 + w1 ); const vec3 dx = ( is_end ? vec3( 0, 0, 0 ) : ( dx_length > 0.0 ? n1 * dx_length : vec3( 0, 0, 0 ) ) ); return subgroupAdd( vec4( dx, ( dx == vec3( 0, 0, 0 ) ) ? 0.0 : 1.0 ) ); িಥ੍໿ͷϦετ͔Βিಥ૬खͷཻࢠͷIDΛऔΓग़͢

const uint iter = constraint_begin( particle_id, constraint_offset ); const bool
is_end = constraint_is_end( iter + constraint_index ); const uint dc = is_end ? 0u : constraint_get( iter + constraint_index ); const vec3 p0 = particle_pool[ particle_offset + particle_id ].position; const vec3 pp0 = particle_pool[ particle_offset + particle_id ].previous_position; const vec3 p1 = particle_pool[ ( is_end ? ( particle_offset + particle_id ) : dc ) ].position; vec3 n1 = particle_pool[ ( is_end ? ( particle_offset + particle_id ) : dc ) ].normal; const bool has_normal = n1 != vec3( 0.0, 0.0, 0.0 ); n1 = has_normal ? n1 : normalize( p0 - p1 ); const bool front = dot( ( pp0 - p1 ), n1 ) >= 0.0; n1 = front ? n1 : -n1; const float thickness = 0.1; const float w0 = particle_pool[ particle_offset + particle_id ].w; const float w1 = particle_pool[ ( is_end ? ( particle_offset + particle_id ) : dc ) ].w; const float d = distance( p0, p1 ); const float dx_length = -( w0*( dot( n1, p0 - p1 ) - thickness ) )/( w0 + w1 ); const vec3 dx = ( is_end ? vec3( 0, 0, 0 ) : ( dx_length > 0.0 ? n1 * dx_length : vec3( 0, 0, 0 ) ) ); return subgroupAdd( vec4( dx, ( dx == vec3( 0, 0, 0 ) ) ? 0.0 : 1.0 ) ); } িಥཻͨ͠ࢠΛিಥ͞Εཻͨࢠͷ๏ઢํ޲ʹҠಈͤͯ͞ ؏௨͠ͳ͍Α͏ʹ͢Δ িಥ͞Εཻͨࢠଆ΋ಉ͡ܭࢉΛ͢ΔͷͰ૒ํͷཻࢠ͕ ཭ΕΔํ޲ʹҠಈ͢Δ

ϝογϡ͔ΒཻࢠΛੜ੒ ଎౓ͱ֎ྗΛ࢖ཻͬͯࢠͷҐஔΛߋ৽ ϨϯμϦϯά ϝογϡ͔Β੍໿Λੜ੒ ϝογϡ͔ΒٯҾ͖Λ࡞Δ ͜͜Ͱڑ཭੍໿ͱ Ұॹʹमਖ਼ ۭؒతϋογϡςʔϒϧʹॻ͘ ۭؒతϋογϡςʔϒϧΛಡΉ ͜͜Ͱিಥ੍໿ͷ
ϦετΛੜ੒ ཻࢠͷҐஔΛमਖ਼ ཻࢠͷҐஔΛ࢖ͬͯ଎౓Λमਖ਼ ཻࢠͷҐஔΛ࢖ͬͯϝογϡΛߋ৽ Ұఆճ਺܁Γฦͨ͠ GBMTF USVF

p0 = p0 − k′ 1 m0 ( p0
− p1 − n) 1 m0 + 1 m1 p0 − p1 p0 − p1 p1 = p1 + l′ 1 m1 ( p0 − p1 − n) 1 m0 + 1 m1 p0 − p1 p0 − p1 ⋮ 0.0Ҏ্1.0ҎԼͰද͞ΕΔ૬ରతͳ߶ੑ஋Λ ʹֻ͚Δ 2ͭͷόω͕ಉ͡௕͞৳ͼ͍ͯΔ࣌ ΑΓݻ͍όωͷํ͕ ʹେ͖ͳӨڹΛ༩͑ΔΑ͏ʹͳΔ Δp p ૬ରతͳ߶ੑ஋ͬͯԿͩ ܭଌͰ͖Δ߶ੑ஋Λ࢖͍͍ͨ

Miles Macklin, Matthias Müller, and Nuttapong Chentanez. 2016. XPBD: position-based
simulation of compliant constrained dynamics. In Proceedings of the 9th International Conference on Motion in Games (MIG '16). Association for Computing Machinery, New York, NY, USA, 49–54. https://doi.org/10.1145/2994258.2994272 XPBD- Position-Based Simulation of Compliant Constrained Dynamics

M d2p dt2 = F = − ∇p U (p)
ཻࢠʹՃΘΔྗ=ϙςϯγϟϧͷޯ഑ M ( pt+Δt − 2pt + pt−Δt Δt2 ) = − ∇p U (pt+Δt) U = 1 2 kd2 C01 (p0 , p1 ) = p0 − p1 − n U01 (p0 , p1) = 1 2 k ( p0 − p1 − n) 2 = 1 2 kC01 (p0 , p1) C01 (p0 , p1) όωͷϙςϯγϟϧ ߶ੑ஋

−∇p U (p) = − ∇p C (p) T α−1C
(p) U01 (p0 , p1) = 1 2 k ( p0 − p1 − n) 2 = 1 2 kC01 (p0 , p1) C01 (p0 , p1) C = [C01 , C02 , . . . ] α = 1 k01 0 0 ⋯ 0 1 k02 0 0 0 1 k03 ⋮ ⋱ શͯͷڑ཭੍໿ΛϕΫτϧʹ͢Δͱ ҎԼͷΑ͏ʹॻ͚Δ U (p) = 1 2 C (p) T α−1C (p)

= ∇C (pt+Δt) T λt+Δt 1 Δt2 ( M (
pt+Δt − 2pt + pt−Δt Δt2 ) + ∇p U (pt+Δt) ) Δt2 = 0 M (pt+Δt − ˜ p) − ∇p C (pt+Δt) T λt+Δt = 0 C (pt+Δt) + ˜ αλt+Δt = 0 ͱ ͱ Λ͜ͷΑ͏ʹఆٛ͢Δͱ ˜ p ˜ α λ ͜Ε͕ಘΒΕΔ M ( pt+Δt − 2pt + pt−Δt Δt2 ) = − ∇p U (pt+Δt) ˜ p = 2pt − pt−Δt = pt + Δtvt ˜ α = α Δt2 λ = − ˜ α−1C (p) −∇p U (p) = − ∇p C (p) T α−1C (p)

M (pt+Δt − ˜ p) − ∇p C (pt+Δt) T
λt+Δt = 0 C (pt+Δt) + ˜ αλt+Δt = 0 g (pt+Δt , λt+Δt) = M (pt+Δt − ˜ p) − ∇p C (pt+Δt) T λt+Δt h (pt+Δt , λt+Δt) = C (pt+Δt) + ˜ αλt+Δt ࠨลͷؔ਺ʹ໊લΛ෇͚Δ g (pt+Δt , λt+Δt ) + ∂g (pt+Δt , λt+Δt ) ∂p Δp + ∂g (pt+Δt , λt+Δt ) ∂λ Δλ = 0 h (pt+Δt , λt+Δt ) + ∂h (pt+Δt , λt+Δt ) ∂p Δp + ∂h (pt+Δt , λt+Δt ) ∂λ Δλ = 0 ͱ Λ2ม਺1֊ͷςΠϥʔల։ʹ͔͚Δ g h g (pt+Δt , λt+Δt ) + ∂g (pt+Δt , λt+Δt ) ∂p Δp − ∇p C (xt+Δt) T Δλ = 0 h (pt+Δt , λt+Δt ) + ∇p C (pt+Δt) Δp + ˜ αΔλ = 0

∂p ∂λ g (pt+Δt , λt+Δt ) + ∂g (pt+Δt
, λt+Δt ) ∂p Δp − ∇p C (xt+Δt) T Δλ = 0 h (pt+Δt , λt+Δt ) + ∇p C (pt+Δt) Δp + ˜ αΔλ = 0 ∂g(pt+Δt , λt+Δt) ∂p −∇p C (xt+Δt) T ∇p C (pt+Δt) ˜ α [ Δp Δλ] = − g (pt+Δt , λt+Δt ) h (pt+Δt , λt+Δt ) g (pt+Δt , λt+Δt) = M (pt+Δt − ˜ p) − ∇p C (pt+Δt) T λt+Δt ∂g (pt+Δt , λt+Δt) ∂p = M − ∂∇p C (pt+Δt) T λt+Δt ∂p ͜ΕΛແࢹ

M −∇p C (pt+Δt) T ∇p C (pt+Δt) ˜ α
[ Δp Δλ] = − g (pt+Δt , λt+Δt ) h (pt+Δt , λt+Δt ) ཻࢠ͕ಈ͍͍ͯͳ͍( )͔ͭ ੍໿͕ຬͨ͞Ε͍ͯΔ( )ͷ࣌ ʹͳΔ Δp = 0 C(pt+Δt ) = 0 g (pt+Δt , λt+Δt) = 0 1αΠΫϧຖͷཻࢠͷมҠ ͕খ͘͞ 1αΠΫϧຖʹཻࢠͷҐஔ੍͕໿Λຬͨ͢ঢ়ଶʹमਖ਼͞Ε͍ͯΔͳΒ ͸ৗʹখ͞ͳ஋ʹͳΔͷͰ ͰۙࣅͰ͖Δ Δp g (pt+Δt , λt+Δt) 0

γϡʔΞิߦྻ ਖ਼ଇͳϒϩοΫߦྻ ͕͋Δ࣌ ͳΒ [ A B C D] [
A B C D] [ x y] = [ a b] [D − CA−1B] [y] = [b] M −∇p C (pt+Δt) T ∇p C (pt+Δt) ˜ α [ Δp Δλ] = − [ 0 h (pt+Δt , λt+Δt )]

[ ˜ α + ∇p C (pt+Δt) M−1 ∇p C
(xt+Δt) T ] Δλ = − h (pt+Δt , λt+Δt ) h (pt+Δt , λt+Δt) = C (pt+Δt) + ˜ αλt+Δt [ ˜ α + ∇p C (pt+Δt) M−1 ∇p C (xt+Δt) T ] Δλ = − C (pt+Δt) − ˜ αλt+Δt Δλij = −Cij (pi,t+Δt) − 1 kij Δt2 λij,t+Δt ∇p C (pi,t+Δt) M−1 ∇p C (xi,t+Δt) T + 1 kij Δt2 M −∇p C (pt+Δt) T ∇p C (pt+Δt) ˜ α [ Δp Δλ] = − [ 0 h (pt+Δt , λt+Δt )]

Δλij = −Cij (pi,t+Δt) − 1 kij Δt2 λij,t+Δt ∇p
C (pi,t+Δt) M−1 ∇p C (xi,t+Δt) T + 1 kij Δt2 M −∇p C (pt+Δt) T ∇p C (pt+Δt) ˜ α [ Δp Δλ] = − [ 0 h (pt+Δt , λt+Δt )] MΔp − ∇p C (pt+Δt) T Δλ = 0 Δp − M−1 ∇p C (pt+Δt) T Δλ = 0 Δp = M−1 ∇p C (pt+Δt) T Δλ Δpi = 1 mi ∇p C (pi,t+Δt) T Δλij

C (pi,t+Δt) M−1 ∇p C (xi,t+Δt) T + 1 kij Δt2 Δpi = 1 mi ∇p C (pi,t+Δt) T Δλij ߶ੑ஋ ͜Ε͕ਖ਼͍͠߶ੑ஋ͷֻ͚ํ ͸ະ஌ͷ஋ͳͷͰ Λߋ৽͢Δࣜͱ Λߋ৽͢ΔࣜΛ Ұॹʹฒ΂ͯΨ΢εɾβΠσϧ๏ʹ͔͚Δ λij,t+Δt λij,t+Δt p

C (pi,t+Δt) M−1 ∇p C (xi,t+Δt) T + 1 kij Δt2 Δpi = 1 mi ∇p C (pi,t+Δt) T Δλij ݁Ռతʹղ͖ํ͕ ϥάϥϯδϡͷະఆ৐਺๏ͷܗʹͳΔҝ ͸͠͹͠͹ Lagrange Multiplierͱݺ͹ΕΔ Δλij

C (pi,t+Δt) M−1 ∇p C (xi,t+Δt) T + 1 kij Δt2 ੍໿͕ ͱ ͷؒͷڑ཭੍໿ͷ৔߹ ͸ ͕ ·ͨ͸ ͷ৔߹͚ͩ୯Ґ௕ͷϕΫτϧΛฦ͢ pi pj ∇p C p pi pj Δλij = −( pi − pj − n) − 1 kij Δt2 λij,t+Δt 1 mi + 1 mj + 1 kij Δt2 Αͬͯ͜͏ͳΔ

p0 = p0 − 1 m0 (( p0 − p1
− n) − 1 k01 Δt2 λ01,t+Δt) 1 m0 + 1 m1 + 1 k01 Δt2 p0 − p1 p0 − p1 p0 = p0 − k′ 01 1 m0 ( p0 − p1 − n) 1 m0 + 1 m1 p0 − p1 p0 − p1 XPBD PBD

vec4 xpbd_distance_constraint_dx( uint particle_id, uint particle_offset, uint distance_constraint_offset, uint lambda_offset,
uint constraint_index, float dt ) { const uint iter = distance_constraint_begin( particle_id, distance_constraint_offset ); const uint lambda_iter = distance_lambda_begin( particle_id, lambda_offset ); const bool is_end = distance_constraint_is_end( iter + constraint_index ); const distance_constraint_type dc = is_end ? distance_constraint_type( 0, 0.0, 0.0, 0.0 ) : distance_constraint_get( iter + constraint_index ); float lambda = ( is_end || lambda_offset == 0xFFFFFFFF ) ? 0.0 : distance_lambda_get( lambda_iter + constraint_index ); const vec3 p0 = particle_pool[ particle_offset + particle_id ].position.xyz; const vec3 p1 = particle_pool[ particle_offset + ( is_end ? particle_id : dc.to_id ) ].position.xyz; const float w0 = particle_pool[ particle_offset + particle_id ].w; const float w1 = particle_pool[ particle_offset + ( is_end ? particle_id : dc.to_id ) ].w; const float d = distance( p0, p1 ); const float alpha_dt2 = 1.0 / ( dc.stiffness * dt * dt ); const float c = d - dc.natural_distance; const vec3 grad_c = ( p0 - p1 ) / max( d, 0.01 ); const float dl = ( -c - alpha_dt2 * lambda ) / ( w0 + w1 + alpha_dt2 ); const vec3 dx = w0 * grad_c * dl; ཻࢠͷҐஔΛमਖ਼ 1੍໿=1εϨουͰ࣮ߦ͢Δ

vec4 xpbd_distance_constraint_dx( uint particle_id, uint particle_offset, uint distance_constraint_offset, uint lambda_offset,
uint constraint_index, float dt ) { const uint iter = distance_constraint_begin( particle_id, distance_constraint_offset ); const uint lambda_iter = distance_lambda_begin( particle_id, lambda_offset ); const bool is_end = distance_constraint_is_end( iter + constraint_index ); const distance_constraint_type dc = is_end ? distance_constraint_type( 0, 0.0, 0.0, 0.0 ) : distance_constraint_get( iter + constraint_index ); float lambda = ( is_end || lambda_offset == 0xFFFFFFFF ) ? 0.0 : distance_lambda_get( lambda_iter + constraint_index ); const vec3 p0 = particle_pool[ particle_offset + particle_id ].position.xyz; const vec3 p1 = particle_pool[ particle_offset + ( is_end ? particle_id : dc.to_id ) ].position.xyz; const float w0 = particle_pool[ particle_offset + particle_id ].w; const float w1 = particle_pool[ particle_offset + ( is_end ? particle_id : dc.to_id ) ].w; const float d = distance( p0, p1 ); const float alpha_dt2 = 1.0 / ( dc.stiffness * dt * dt ); const float c = d - dc.natural_distance; const vec3 grad_c = ( p0 - p1 ) / max( d, 0.01 ); const float dl = ( -c - alpha_dt2 * lambda ) / ( w0 + w1 + alpha_dt2 ); const vec3 dx = w0 * grad_c * dl; ΛಡΉ λij,t+Δt ͸1λΠϜεςοϓຖʹθϩΫϦΞ͢Δ λ

const distance_constraint_type dc = is_end ? distance_constraint_type( 0, 0.0, 0.0,
0.0 ) : distance_constraint_get( iter + constraint_index ); float lambda = ( is_end || lambda_offset == 0xFFFFFFFF ) ? 0.0 : distance_lambda_get( lambda_iter + constraint_index ); const vec3 p0 = particle_pool[ particle_offset + particle_id ].position.xyz; const vec3 p1 = particle_pool[ particle_offset + ( is_end ? particle_id : dc.to_id ) ].position.xyz; const float w0 = particle_pool[ particle_offset + particle_id ].w; const float w1 = particle_pool[ particle_offset + ( is_end ? particle_id : dc.to_id ) ].w; const float d = distance( p0, p1 ); const float alpha_dt2 = 1.0 / ( dc.stiffness * dt * dt ); const float c = d - dc.natural_distance; const vec3 grad_c = ( p0 - p1 ) / max( d, 0.01 ); const float dl = ( -c - alpha_dt2 * lambda ) / ( w0 + w1 + alpha_dt2 ); const vec3 dx = w0 * grad_c * dl; if( !( is_end || lambda_offset == 0xFFFFFFFF ) ) { distance_lambda_set( lambda_iter + constraint_index, is_end ? 0.0 : lambda + dl ); } return subgroupAdd( vec4( dx, is_end ? 0.0 : 1.0 ) ); } Δλij = −( pi − pj − n) − 1 kij Δt2 λij,t+Δt 1 mi + 1 mj + 1 kij Δt2

const distance_constraint_type dc = is_end ? distance_constraint_type( 0, 0.0, 0.0,
0.0 ) : distance_constraint_get( iter + constraint_index ); float lambda = ( is_end || lambda_offset == 0xFFFFFFFF ) ? 0.0 : distance_lambda_get( lambda_iter + constraint_index ); const vec3 p0 = particle_pool[ particle_offset + particle_id ].position.xyz; const vec3 p1 = particle_pool[ particle_offset + ( is_end ? particle_id : dc.to_id ) ].position.xyz; const float w0 = particle_pool[ particle_offset + particle_id ].w; const float w1 = particle_pool[ particle_offset + ( is_end ? particle_id : dc.to_id ) ].w; const float d = distance( p0, p1 ); const float alpha_dt2 = 1.0 / ( dc.stiffness * dt * dt ); const float c = d - dc.natural_distance; const vec3 grad_c = ( p0 - p1 ) / max( d, 0.01 ); const float dl = ( -c - alpha_dt2 * lambda ) / ( w0 + w1 + alpha_dt2 ); const vec3 dx = w0 * grad_c * dl; if( !( is_end || lambda_offset == 0xFFFFFFFF ) ) { distance_lambda_set( lambda_iter + constraint_index, is_end ? 0.0 : lambda + dl ); } return subgroupAdd( vec4( dx, is_end ? 0.0 : 1.0 ) ); } Δpi = 1 mi ∇p C (pi,t+Δt) T Δλij

( is_end || lambda_offset == 0xFFFFFFFF ) ? 0.0 :
distance_lambda_get( lambda_iter + constraint_index ); const vec3 p0 = particle_pool[ particle_offset + particle_id ].position.xyz; const vec3 p1 = particle_pool[ particle_offset + ( is_end ? particle_id : dc.to_id ) ].position.xyz; const float w0 = particle_pool[ particle_offset + particle_id ].w; const float w1 = particle_pool[ particle_offset + ( is_end ? particle_id : dc.to_id ) ].w; const float d = distance( p0, p1 ); const float alpha_dt2 = 1.0 / ( dc.stiffness * dt * dt ); const float c = d - dc.natural_distance; const vec3 grad_c = ( p0 - p1 ) / max( d, 0.01 ); const float dl = ( -c - alpha_dt2 * lambda ) / ( w0 + w1 + alpha_dt2 ); const vec3 dx = w0 * grad_c * dl; if( !( is_end || lambda_offset == 0xFFFFFFFF ) ) { distance_lambda_set( lambda_iter + constraint_index, is_end ? 0.0 : lambda + dl ); } return subgroupAdd( vec4( dx, is_end ? 0.0 : 1.0 ) ); } ৽͍͠ ͷ஋Λอଘ λij,t+Δt ಉҰSubgroupͰܭࢉͨ͠࠷େ32ݸͷ ͷ૯࿨ΛٻΊΔ Δp Ұॹʹ੍໿ͷݸ਺΋ٻΊ͓͍ͯͯޙͰฏۉʹͰ͖ΔΑ͏ʹ͢Δ

ϝογϡ͔ΒཻࢠΛੜ੒ ଎౓ͱ֎ྗΛ࢖ཻͬͯࢠͷҐஔΛߋ৽ ϨϯμϦϯά ϝογϡ͔Β੍໿Λੜ੒ ϝογϡ͔ΒٯҾ͖Λ࡞Δ ۭؒతϋογϡςʔϒϧʹॻ͘ ۭؒతϋογϡςʔϒϧΛಡΉ ΛθϩΫϦΞ͢Δ λ ͜͜ͷ಺༰͕
มΘΔ ཻࢠͷҐஔΛमਖ਼ ཻࢠͷҐஔΛ࢖ͬͯ଎౓Λमਖ਼ ཻࢠͷҐஔΛ࢖ͬͯϝογϡΛߋ৽ Ұఆճ਺܁Γฦͨ͠ GBMTF USVF

Miles Macklin and Matthias Müller. 2013. Position based fluids. ACM
Trans. Graph. 32, 4, Article 104 (July 2013), 12 pages. https://doi.org/10.1145/2461912.2461984 Position Based Fluids

Dv Dt = − 1 ρ ∇p + μ∇2v +
Fext ྲྀମͷӡಈํఔࣜ φϏΤɾετʔΫεํఔࣜ ߲࣌ؒ ѹྗ߲ ೪ੑ߲ ֎ྗ߲ ྲྀମ͸ѹྗ ͕ ߴ͍ํ͔Β௿͍ํʹྲྀΕ͕ͨΔ p ྲྀΕ͍ͯΔྲྀମͷपғͷྲྀମ͸ ҰॹʹྲྀΕ͕ͨΔ ॏྗ౳֎͔Β ՃΘΔྗ

Dv Dt = − 1 ρ ∇p + μ∇2v +
Fext ྲྀମΛཻࢠͷू·Γͱݟ၏͢ͱ ྲྀମ ཻࢠ ཻࢠಉ࢜ͷিಥ ཻࢠಉ࢜ͷຎࡲ

Dv Dt = − 1 ρ ∇p + μ∇2v +
Fext ਫ౳ͷྲྀମͷ෼ࢠ͸͕͍ͨʹিಥͯ͠֎ʹඈΜͰߦ͜͏ͱ͢Δ͕ େؾѹͰԡ͚͑ͭ͞ΒΕ͍ͯΔҝ·ͱ·ΓΛอ͍ͬͯΔ େؾѹ

Dv Dt = − 1 ρ ∇p + μ∇2v +
Fext ۭؾͷཻࢠ͸༻ҙ͍ͯ͠ͳ͍ͷͰ େؾѹ͸ࣗવʹ͸ੜ͡ͳ͍ ୯७ʹিಥ੍໿ͰཻࢠͷҐஔΛमਖ਼͢Δͱ ྲྀମ͕രൃ࢛ࢄ͢Δ

Dv Dt = − 1 ρ ∇p + μ∇2v +
Fext ྲྀମͷཻࢠ͸Ұఆͷີ౓Λอ͕ͪͨΔੑ࣭͕͋Δͱߟ͑Δ Ci (p0 , p1 , ⋯) = ρi ρinitial − 1 ൪໨ͷཻࢠͷपғͷྲྀମͷີ౓ͱ i ྲྀମͷॳظঢ়ଶͰͷີ౓͕ ͍ۙ΄Ͳྑ͍ ີ౓੍໿

Dv Dt = − 1 ρ ∇p + μ∇2v +
Fext ྲྀମͷཻࢠ͸Ұఆͷີ౓Λอ͕ͪͨΔੑ࣭͕͋Δͱߟ͑Δ Ci (p0 , p1 , ⋯) = ρi ρinitial − 1 ൪໨ͷཻࢠͷपғͷྲྀମͷີ౓ͱ i ྲྀମͷॳظঢ়ଶͰͷີ౓͕ ͍ۙ΄Ͳྑ͍ ີ౓੍໿ ͜ΕͲ͏΍ͬͯٻΊΔͷ?

MONAGHAN, J. J. 1992. Smoothed particle hydrodynamics. Annual Review of
Astronomy and Astrophysics 30, 1, 543– 574. SMOOTHED PARTICLE HYDRODYNAMICS https://doi.org/10.1146/annurev.aa.30.090192.002551

ρi = ∑ j mj W (pi − pj ,
h) ൪໨ͷཻࢠͷपғͷྲྀମͷີ౓͸ i ۙ͘ʹ͋Δશͯͷཻࢠʹ͍ͭͯ ڑ཭ʹԠͨ͡Өڹ౓ͱ ཻࢠͷ࣭ྔͷੵΛٻΊͯ ૯࿨ΛऔΔͱٻ·ΔΑ શͯͷྲྀମཻࢠ͕ಉ࣭͡ྔͳΒ͔ࣜΒ ΛফͤΔ mj

ͱͳΔΑ͏ͳཻࢠͷҠಈ ΛٻΊ͍ͨ C (p + Δp) = 0 Δp ΛςΠϥʔల։
C C (p + Δp) ≃ C (p) + ∇CTΔp = 0

C (p + Δp) ≃ C (p) + ∇CT ∇Cλ
ͳͷͰ Δpi = − Ci (p) ∑ k ∇pk Ci (p) 2 ∇p Ci (p) C (p + Δp) ≃ C (p) + ∇CTΔp = 0 λ = − C0 (p) ∑ k ∇pk C0 (p) 2 − C1 (p) ∑ k ∇pk C1 (p) 2 ⋮

Ci (p) = 1 ρinitial ∑ j W (pi −
pj , h) − 1 ∇pk Ci (p) = 1 ρinitial ∑ j ∇pk W (pi − pj , h) = 1 ρinitial ∑ j ∇pk W (pi − pj , h) (k = i) −∇pk W (pi − pj , h) (k = j)

= 1 ρinitial ∑ j ∇pk W (pi − pj
, h) (k = i) −∇pk W (pi − pj , h) (k = j) ∇pk Ci (p) C (p + Δp) ≃ C (p) + ∇CT ∇Cλ = 0 Δp = − C(p) ∑ k ∇pk C(p) 2 ∇p C(p) Δpi = 1 ρinitial ∑ j (λi + λj)∇W (pi − pj , h) λ = − C0 (p) ∑ k ∇pk C0 (p) 2 − C1 (p) ∑ k ∇pk C1 (p) 2 ⋮ ͜ͷ ͰཻࢠͷҐஔΛमਖ਼͍͚ͯ͠͹ྑ͍ Δpi

λi = − Ci (p) ∑ k ∇pk Ci (p)
2 ཻࢠؒͷڑ཭ W ∇W ∇pk Ci (p) = 1 ρinitial ∑ j ∇pk W (pi − pj , h) ͸Χʔωϧؔ਺ͷڥք෇ۙͰ0ʹͳΔ ∇W

λi = − Ci (p) ∑ k ∇pk Ci (p)
2 ཻࢠؒͷڑ཭ W ∇W ∇pk Ci (p) = 1 ρinitial ∑ j ∇pk W (pi − pj , h) ۙ๣ͷཻࢠ͕શͯ Χʔωϧؔ਺ͷڥք෇ۙʹ͋Δͱ ͜ͷ෦෼͕௕͞0ʹͳΔ

λi = − Ci (p) ∑ k ∇pk Ci (p)
2 ཻࢠؒͷڑ཭ W ∇W ∇pk Ci (p) = 1 ρinitial ∑ j ∇pk W (pi − pj , h) ͦͷ݁Ռ͜ͷ෼਺ͷ෼฼͕0ʹͳΓ NaN͕ᷓΕग़͢

λi = − Ci (p) ∑ k ∇pk Ci (p)
2 + ϵ ཻࢠؒͷڑ཭ W ∇W খ͍͞ఆ਺ ΛՃ͑ͯ ෼฼͕ઈରʹ0ʹͳΒͳ͍Α͏ʹ͢Δ ϵ

float poly6_kernel( vec3 r, float h ) { const float
rd = length( r ); const float alpha = 315.0 / ( 64.0 * pi * pow( h, 9 ) ); return ( rd < h ) ? ( alpha * pow( h * h - rd * rd, 3 ) ) : 0.0; } vec3 poly6_kernel_gradient( vec3 r, float h ) { const float rd = length( r ); const float alpha = 315.0 / ( 64.0 * pi * pow( h, 9 ) ); return ( rd < h ) ? ( float( -6 * alpha * pow( h * h - rd * rd, 2 ) ) * r ) : vec3( 0.0, 0.0, 0.0 ); } W ∇W SPHͰΑ͘࢖ΘΕΔPoly6Χʔωϧؔ਺Λ࣮૷

float pbd_fluid_get_rho( uint particle_id, uint particle_offset, uint constraint_offset, float radius
) { const uint iter = constraint_begin( particle_id, constraint_offset ); const vec3 center = particle_pool[ particle_offset + particle_id ].position; float rho = 0.0; for( uint constraint_block_id = 0u; constraint_block_id != 4u; constraint_block_id++ ) { uint gcid = iter + gl_SubgroupSize * constraint_block_id + gl_SubgroupInvocationID; const bool is_end = constraint_is_end( gcid ); const uint dc = is_end ? 0u : constraint_get( gcid ); const vec3 neighbor = particle_pool[ ( is_end ? ( particle_offset + particle_id ) : dc ) ].position; rho += subgroupAdd( is_end ? 0.0 : poly6_kernel( center - neighbor, radius ) ); } return rho; } float pbd_fluid_get_lambda( float rho0, uint particle_id, uint particle_offset, uint constraint_offset, float radius ) { const uint iter = constraint_begin( particle_id, constraint_offset ); const vec3 center = particle_pool[ particle_offset + particle_id ].position; float dpici = 0.0; ཻࢠͷҐஔΛमਖ਼ 4੍໿=1εϨουͰ࣮ߦ͢Δ ۙ๣ཻࢠ͸িಥ੍໿Ͱ࢖ͬͨ΋ͷͱಉ͡ ۭؒతϋογϡςʔϒϧΛ࢖ͬͯಘΔ

float pbd_fluid_get_rho( uint particle_id, uint particle_offset, uint constraint_offset, float radius
) { const uint iter = constraint_begin( particle_id, constraint_offset ); const vec3 center = particle_pool[ particle_offset + particle_id ].position; float rho = 0.0; for( uint constraint_block_id = 0u; constraint_block_id != 4u; constraint_block_id++ ) { uint gcid = iter + gl_SubgroupSize * constraint_block_id + gl_SubgroupInvocationID; const bool is_end = constraint_is_end( gcid ); const uint dc = is_end ? 0u : constraint_get( gcid ); const vec3 neighbor = particle_pool[ ( is_end ? ( particle_offset + particle_id ) : dc ) ].position; rho += subgroupAdd( is_end ? 0.0 : poly6_kernel( center - neighbor, radius ) ); } return rho; } float pbd_fluid_get_lambda( float rho0, uint particle_id, uint particle_offset, uint constraint_offset, float radius ) { const uint iter = constraint_begin( particle_id, constraint_offset ); const vec3 center = particle_pool[ particle_offset + particle_id ].position; float dpici = 0.0; ཻࢠ ʹ͓͚Δྲྀମͷີ౓ ΛٻΊΔ i ρi ∑ j W (pi − pj , h) pj

uint constraint_offset, float radius ) { const uint iter =
constraint_begin( particle_id, constraint_offset ); const vec3 center = particle_pool[ particle_offset + particle_id ].position; float dpici = 0.0; float dpjci = 0.0; { vec3 sum = vec3( 0.0, 0.0, 0.0 ); for( uint constraint_block_id = 0u; constraint_block_id != 4u; constraint_block_id++ ) { uint gcid = iter + gl_SubgroupSize * constraint_block_id + gl_SubgroupInvocationID; const bool is_end = constraint_is_end( gcid ); const uint dc = is_end ? 0u : constraint_get( gcid ); const vec3 neighbor = particle_pool[ ( is_end ? ( particle_offset + particle_id ) : dc ) ].position; vec3 grad = is_end ? vec3( 0.0, 0.0, 0.0 ) : poly6_kernel_gradient( center - neighbor, radius ) / rho0; sum += subgroupAdd( grad ); dpjci += subgroupAdd( dot( grad, grad ) ); } dpici = dot( sum, sum ); } const float eps = 1000.0; return -( pbd_fluid_get_rho( particle_id, particle_offset, constraint_offset, radius )/rho0 - 1.0f ) / ( dpici + dpjci + eps ); } vec4 pbd_fluid_constraint( float rho0, uint particle_id, λi = − Ci (p) ∑ k ∇pk Ci (p) 2 + ϵ −Ci (p) ∑ k ∇pk Ci (p) 2 + ϵ ∑ j − 1 ρinitial ∇pk W (pi − pj , h) 2 1 ρinitial ∑ j ∇pk W (pi − pj , h) 2

uint constraint_offset, float radius ) { const uint iter =
constraint_begin( particle_id, constraint_offset ); const vec3 center = particle_pool[ particle_offset + particle_id ].position; float lambda_i = pbd_fluid_get_lambda( rho0, particle_id, particle_offset, constraint_offset, radius ); vec3 sum = vec3( 0.0f, 0.0f, 0.0f ); uint count = 0; for( uint constraint_block_id = 0u; constraint_block_id != 4u; constraint_block_id++ ) { uint gcid = iter + gl_SubgroupSize * constraint_block_id + gl_SubgroupInvocationID; const bool is_end = constraint_is_end( gcid ); const uint dc = is_end ? 0u : constraint_get( gcid ); const vec3 neighbor = particle_pool[ ( is_end ? ( particle_offset + particle_id ) : dc ) ].position; float lambda_j = is_end ? 0.0 : pbd_fluid_get_lambda( rho0, dc - particle_offset, particle_offset, constraint_offset, radius ); vec3 dxdw = poly6_kernel_gradient( center - neighbor, radius ); dxdw *= lambda_i + lambda_j; sum += subgroupAdd( is_end ? vec3( 0.0, 0.0, 0.0 ) : dxdw ); count += subgroupAdd( is_end ? 0 : 1 ); } sum /= rho0; return vec4( sum, float( count ) ); } Δpi = 1 ρinitial ∑ j (λi + λj)∇W (pi − pj , h) ∇W (pi − pj , h) ∑ j (λi + λj)∇W (pi − pj , h) 1 ρinitial ∑ j (λi + λj)∇W (pi − pj , h) λi λj

Dv Dt = − 1 ρ ∇p + μ∇2v +
Fext ྲྀମͷཻࢠ͸ۙ͘ͷཻࢠʹҾͬுΒΕΔ ཻࢠؒͷ૬ର଎౓Λখ͘͢͞Δํ޲ʹྗ͕ಇ͘

Dv Dt = − 1 ρ ∇p + μ∇2v +
Fext ྲྀମͷཻࢠ͸ۙ͘ͷཻࢠʹҾͬுΒΕΔ ཻࢠؒͷ૬ର଎౓Λখ͘͢͞Δํ޲ʹྗ͕ಇ͘ ଎౓ʹґଘ͢Δݱ৅ Position BasedͰѻ͍ʹ͍͘

ϝογϡ͔ΒཻࢠΛੜ੒ ଎౓ͱ֎ྗΛ࢖ཻͬͯࢠͷҐஔΛߋ৽ ϨϯμϦϯά ϝογϡ͔Β੍໿Λੜ੒ ϝογϡ͔ΒٯҾ͖Λ࡞Δ ۭؒతϋογϡςʔϒϧʹॻ͘ ۭؒతϋογϡςʔϒϧΛಡΉ ΛθϩΫϦΞ͢Δ λ ཻࢠͷҐஔΛमਖ਼
ཻࢠͷҐஔΛ࢖ͬͯ଎౓Λमਖ਼ ཻࢠͷҐஔΛ࢖ͬͯϝογϡΛߋ৽ Ұఆճ਺܁Γฦͨ͠ GBMTF USVF XSPHͷਓ޻೪ੑͰ଎౓Λमਖ਼ ଎౓ΛٻΊͨޙͰܭࢉ

vnew i = vi + c∑ j (vj − vi)
W (pi − pj , h) ೪ੑͷӨڹͷ େ͖͞ΛܾΊΔఆ਺ ཻࢠͷҐஔ͔ΒٻΊͨ଎౓ ೪ੑΛߟྀͨ͠଎౓ ۙ๣ͷཻࢠͱͷ଎౓ࠩ XSPH ਓ޻೪ੑ

vec3 pbd_fluid_xsph_viscosity( uint particle_id, uint particle_offset, uint constraint_offset, float radius
) { const uint iter = constraint_begin( particle_id, constraint_offset ); const vec3 center = particle_pool[ particle_offset + particle_id ].position; const vec3 center_v = particle_pool[ particle_offset + particle_id ].velocity; vec3 sum = vec3( 0.0f, 0.0f, 0.0f ); uint count = 0; for( uint constraint_block_id = 0u; constraint_block_id != 4u; constraint_block_id++ ) { uint gcid = iter + gl_SubgroupSize * constraint_block_id + gl_SubgroupInvocationID; const bool is_end = constraint_is_end( gcid ); const uint dc = is_end ? 0u : constraint_get( gcid ); const vec3 neighbor = particle_pool[ ( is_end ? ( particle_offset + particle_id ) : dc ) ].position; const vec3 neighbor_v = particle_pool[ ( is_end ? ( particle_offset + particle_id ) : dc ) ].velocity; const float w = poly6_kernel( center - neighbor, radius ); const vec3 vij = neighbor_v - center_v; sum += subgroupAdd( is_end ? vec3( 0.0, 0.0, 0.0 ) : vij * w ); count += subgroupAdd( is_end ? 0 : 1 ); } sum *= 0.001; return sum; } 941)ͷਓ޻೪ੑͰ଎౓Λमਖ਼ 4੍໿=1εϨουͰ࣮ߦ͢Δ W (pi − pj , h) vj vi ∑ j (vj − vi) W (pi − pj , h) c∑ j (vj − vi) W (pi − pj , h)

Miles Macklin, Kier Storey, Michelle Lu, Pierre Terdiman, Nuttapong Chentanez,
Stefan Jeschke, and Matthias Müller. 2019. Small steps in physics simulation. In Proceedings of the 18th annual ACM SIGGRAPH/ Eurographics Symposium on Computer Animation (SCA '19). Association for Computing Machinery, New York, NY, USA, Article 2, 1–7. https://doi.org/10.1145/3309486.3340247 Small Steps in Physics Simulation

ϝογϡ͔ΒཻࢠΛੜ੒ ଎౓ͱ֎ྗΛ࢖ཻͬͯࢠͷҐஔΛߋ৽ ϨϯμϦϯά ϝογϡ͔Β੍໿Λੜ੒ ϝογϡ͔ΒٯҾ͖Λ࡞Δ ۭؒతϋογϡςʔϒϧʹॻ͘ ۭؒతϋογϡςʔϒϧΛಡΉ ϨϯμϦϯά ΛθϩΫϦΞ͢Δ λ
ཻࢠͷҐஔΛमਖ਼ ཻࢠͷҐஔΛ࢖ͬͯ଎౓Λमਖ਼ ཻࢠͷҐஔΛ࢖ͬͯϝογϡΛߋ৽ Ұఆճ਺܁Γฦͨ͠ GBMTF USVF

ϝογϡ͔ΒཻࢠΛੜ੒ ଎౓ͱ֎ྗΛ࢖ཻͬͯࢠͷҐஔΛߋ৽ ཻࢠͷҐஔΛ࢖ͬͯ଎౓Λमਖ਼ ϨϯμϦϯά Ұఆճ਺܁Γฦͨ͠ GBMTF ϝογϡ͔Β੍໿Λੜ੒ ϝογϡ͔ΒٯҾ͖Λ࡞Δ ۭؒతϋογϡςʔϒϧʹॻ͘ ۭؒతϋογϡςʔϒϧΛಡΉ
ཻࢠͷҐஔΛ࢖ͬͯϝογϡΛߋ৽ ΛθϩΫϦΞ͢Δ λ ཻࢠͷҐஔΛमਖ਼ USVF

U01 (p0 , p1) = 1 2 k ( p0
− p1 − n) 2 ·ͨ͸ ͕֎ྗʹΑͬͯҠಈ͠ ͜ͷ෦෼ͷઈର஋͕૿͑Δ p0 p1 ϙςϯγϟϧ͸ͦͷ2৐Ͱ૿͑Δ ͕େ͖͍΄Ͳ1λΠϜεςοϓͰͷ ͷมԽ͸େ͖͘ͳΔ Δt p ͦͷ2৐Ͱ૿͑ΔϙςϯγϟϧΛ ࠷খʹ͢ΔΑ͏ʹ Λमਖ਼͢Δͷ͕൓෮ܭࢉͷ໨త p ൓෮ճ਺Λ૿΍͢ΑΓ Λࡉ͔͘ࠁΉํ͕ શମͱͯ͠গͳ͍ܭࢉճ਺Ͱ Λमਖ਼Ͱ͖Δ Δt p

Matthias Müller, Miles Macklin, Nuttapong Chentanez, Stefan Jeschke, and Tae-Yong
Kim. 2020. Detailed rigid body simulation with extended position based dynamics. In Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation (SCA '20). Eurographics Association, Goslar, DEU, Article 10, 1–12. Detailed Rigid Body Simulation with Extended Position Based Dynamics https://doi.org/10.1111/cgf.14105

TO BE CONTINUED

Position Based Dynamics ·ͱΊ ඵʹؒʹ߹Θͳ͍ͱࠔΔਓͷͨΊͷ ӄղ๏෺ཧγϛϡϨʔγϣϯ 1 60 ஄ੑମ΋ྲྀମ΋߶ମ΋OK ࣍ճʹଓ͘!

https://fadis.booth.pm/ ٕज़ॻయͰຊग़ͯ͠·͢ ϝογϡ γΣʔμʔ FADIS PRESS Ver. 1.0 /&8

ゲームの物理

ゲームの物理

More Decks by Fadis

Other Decks in Programming

Featured

Transcript