Filetype PDF | Posted on 19 Jan 2023 | 2 years ago
Authentication
digital resources
153x Filetype PDF File size 2.31 MB Source: animation.rwth-aachen.de
Position-Based Rigid Body Dynamics
Crispin Deul, Patrick Charrier and Jan Bender
Graduate School of Excellence Computational Engineering
¨
Technische Universitat Darmstadt, Germany
{deul | charrier | bender}@gsc.tu-darmstadt.de
Abstract
Weproposeaposition-basedapproachforlarge-
scale simulations of rigid bodies at interactive
frame-rates. Our method solves positional con-
straints between rigid bodies and therefore inte-
grates nicely with other position-based methods.
Interaction of particles and rigid bodies through
common constraints enables two-way cou-
pling with deformables. The method exhibits
exceptional performance and stability while
being user-controllable and easy to implement.
Various results demonstrate the practicabil-
ity of our method for the resolution of colli-
sions, contacts, stacking and joint constraints.
Keywords: real-time, rigid body dynamics,
two-way coupling, position-based dynamics
Figure 1: A large-scale mobile simulation in-
1 Introduction volving 127 constraints. A simulation
step requires 0.83ms on average.
Computergamesandvirtualreality applications
strive for ever greater realism in increasingly for methods that are both robust and fast at the
complex virtual environments. Such environ- sametimeis still ongoing.
ments typically include a vast number of dy-
namic deformable and rigid objects. Simulating The position-based dynamics method by
these objects according to the laws of motion re- Mueller et al. [3] has been applied to particle-
quires sophisticated methods and has a long his- based deformable bodies and fluids [4]. It
tory in computer graphics. In this context the allows for very stable animations at interactive
animation of rigid bodies and their interactions frame-rates. However, it does not yet support
under the influence of constraints such as col- rigid body dynamics under the influence of
lisions or joints is perhaps the most researched constraints and collisions.
concern. Numerous physics engines, like Bul- We present a novel approach that extends
let [1] or PhysX [2] have come to public atten- position-based dynamics to rigid bodies while
tion. While many methods achieve results of maintaining its significant properties. Our ap-
great fidelity, their computational effort can not proach excels in large-scale simulations at in-
be justified in real-time applications. The quest teractive frame-rates while being robust under a
vast number of constraints, as depicted in Fig- based on a prediction of the joint state [22], We-
ure1. Furthermore, the proposedmethodiseasy instein et al. [23] solved the nonlinear equation
to implement, controllable and supports two- by a Newton iteration method.
waycoupling of deformable and rigid bodies. Position-based methods became popular in
the last years since they are fast, robust and con-
2 Related Work trollable while no implicit time integration is re-
quired. A survey of position-based methods is
¨
presented by Bender et al. [24]. Muller et al. [3]
The field of rigid body dynamics has a long introduced position-based dynamics as a gener-
history in computer animation and may be de- alized framework capable of solving a large va-
composed into the four subfields integration riety of constraints. They demonstrated its ap-
schemes, collision detection, collision response plication on deformable solids [3, 25] and later
and constraint resolution. The survey of Bender also on fluids [4]. Diziol et al. [26] introduced
et al. [5] describes all subfields in more detail, a shape matching method with a volume conser-
but only the two latter ones are of primary inter- vation for deformable solids to achieve more re-
est in this paper. alistic results. In order to make shape matching
In previous work collision handling was per- ¨
more robust, Muller and Chentanez [25] added
formed by the application of forces [6], by solv- an orientation to the particles of their model.
inglinearcomplementaryproblems[7]orbyus- This allows them to simulate bodies with solid
¨
ing impulse-based methods [8, 9, 10]. Muller et components and basic joints. Although they
al. [3] even introduce a method that resolves col- treat each particle as a rigid body with posi-
lisions by directly modifying positions for the tion and orientation, they apply the standard
simplified case of individual particles. How- position-based constraints that depend solely on
ever, an extension to rigid bodies has not been particle positions. In constrast, our work incor-
shown yet. In contrast to that, Kaufman et porates position and orientation in the deriva-
al. [11] present a method which projects the ve- tion of constraint equations. This allows us
locities of the bodies in order to prevent pene- to form constraints between arbitrary points of
trations. Later, Kaufman et al. [12] introduced rigid bodies.
a staggering approach for frictional contacts.
Multi-impact problems were focused by Smith
et al. [13]. In order to increase the performance 3 Preliminaries
of large-scale simulations, shock propagation
methods were introduced [9, 14]. Furthermore, In this section we first start with a brief intro-
different GPU-based methods were presented to duction to the position-based dynamics frame-
simulate large systems in real-time [15, 16]. work (PBD). Then we continue with the basics
To simulate an articulated system with joints, of rigid body simulation.
equalityconstraintsaredefinedfortherigidbod- The original position-based dynamics ap-
ies. A classic approach to solve these constraints proach [3] is used to simulate a model defined
in real-time is to use reduced coordinate meth- byparticles and constraints between these parti-
ods. Featherstone demonstrated that a simula- cles. The constraints are solved at the position
tion with reduced coordinates can be performed level by applying correction displacements onto
in linear time for acyclic systems [17, 18]. Re- the particle positions. In order to measure the
don et al. [6] introduced an adaptive variant of constraint violation, preview positions of the
this approach which uses a reduced number of particles are computed in a first step by integrat-
degrees of freedom to improve the performance. ing the particle velocities under the influence
Baraff [19] introduced a Lagrange multiplier of external forces. Then the particle positions
method which also has a linear-time complex- are integrated with the new velocities yielding a
ity for acyclic models. Later, Bender [20, 21] symplectic Euler scheme. In the second step the
demonstrated that the idea of Baraff can also systemofconstraintsissolvedinaGauss-Seidel
be applied to impulse-based simulation. While typeiteration. The third step updates the particle
Bender proposed to solve a linearized equation velocities depending on the difference between
the particle positions at the start of the time step applying displacements to the according parti-
and the new particle positions multiplied by the cles. The displacements are computed by solv-
inverse time step size. In a final step the new ing constraints of the following form:
velocities are modified to handle friction and
damping. C(p+∆p)=0, (1)
As mentioned before, we want to solve con- where p = [p ,...,p ]T is the concatenated
straints on rigid bodies, which are defined by six 1 n
parameters. The translational motion parame- vector of particle preview positions and ∆p =
[∆p ,...,∆p ]T contains the corresponding
ters are the position x, the velocity v, and the 1 n
mass m, which rigid bodies have in common correction displacements.
with particles. In addition to that, the three ro- In order to solve a constraint function a first-
tational parameters are the orientation ϑ, the an- order Taylor approximation
gular velocity ω, and the inertia tensor I. C(p+∆p)≈C(p)+∇ C(p)·∆p=0, (2)
Theseparametersareemployedinasymplec- p
tic Euler scheme as follows. The equations for is used to linearize the constraint. However, this
velocity and position integration are defined by equation is underdetermined. By restricting the
F direction of the position correction to the gra-
v(t +h)=v(t )+ externalh dient direction this problem is solved [3]. The
0 0 m
final displacement vector is determined by
x(t +h)=x(t )+v(t +h)h
0 0 0
wC(p)
∆p =− i ∇ C(p), (3)
whiletheequationsfortherotationalparameters i p
P 2 i
w ∇ C(p)
j j p
are given by j
−1 where w is the inverse mass of particle p .
ω(t +h)=ω(t )+hI · i i
0 0 Asanextensiontotheparticleconstraint han-
(τexternal − (ω(t) × (I · ω(t)))) dling of PBD our approach handles constraints
q(t +h)=q(t )+ hω˜(t +h)·q(t ), between rigid bodies. In contrast to particles an
0 0 2 0 0
orientation is associated to rigid bodies. Points
p that are attached to a rigid body with index j
where ω˜ is the quaternion [0,ω ,ω ,ω ]. Af- i
x y z can be described by the following formula
ter the preview positions of the bodies are in-
tegrated with the equations above position con- p(x ,R(ϑ )) = x +R(ϑ )r , (4)
straints are solved in several iterations. How- i j j j j i
ever, in constrast to the original PBD approach, wherex isthecenterofmassandR(ϑ )thero-
that updates the velocities after the constraint j j
solver step, we update the velocities of con- tation matrix of the rigid body with index j. The
straints whenever a correction displacement is local position of the point i in the body frame
is encoded in r . Using the definition of body
applied to a rigid body. The required impulse to i
updatethevelocitiesiscomputedbemultiplying attached points (see Eq. 4) in the Taylor approx-
themassweighteddisplacementwiththeinverse imation of the constraint (see Eq. 2) yields:
time step size. By updating the velocities during C(p(x+∆x,R(ϑ+∆ϑ)))≈
the position iterations we can apply our friction C(p(x,R(ϑ)))+ (5)
resolution approach whose details are presented T T T T T
JC(x,ϑ)·[∆x ,∆ϑ ,...,∆x ,∆ϑ ] ,
in section 6. 1 1 n n
where x = [xT,··· .xT]T and ϑ =
1 n
T T T
4 Constrained Rigid Bodies [ϑ ,··· ,ϑ ] are the vectors containing all
1 n
positions and orientations of the n constrained
In this section we describe how to extend the bodies. Furthermore, the function p com-
PBD solver to solve constraints between rigid putes the concatenated vector of m positions
T T T
bodies. The standard solver works by itera- p = [p (x ,R(ϑ )) ,...,p (x ,R(ϑ )) ]
1 1 1 m n n
tively handling each constraint on its own by that are constrained by C(p). Due to the fact,
that the entries of p depend on a function (see 5 Joints
Eq. 4), the chain rule has to be applied to
the constraint function C(p) to compute its Defining translational joint constraints between
derivative. As a result, the gradient ∇ C(p) two rigid bodies is straight forward using the
p constraint definition of section 4. The transla-
in Eq. 2 is replaced by the k × 6n Jacobian
JC(x,ϑ) = ∂C ∂C of the constraint function tional constraints have in common that the dis-
∂x ∂ϑ tance of two joint points a and b, each attached
with respect to the rigid body positions and
orientations, where k is the dimension of the to one of the two bodies A and B, is constrained
codomain of C(p). to be zero. The only difference is the dimension
LetM bethe6×6massmatrixofrigidbody of the constraint space. A simple translational
j constraint is the ball joint
j with the mass on the first three entries of the
diagonal and the moment of inertia tensor Ij in C(a,b) = a−b=0
the lower right 3 × 3 submatrix. By rearrang- that removesalltranslationaldegreesoffreedom
ing Eq. 5 the same way that led from Eq. 2 to between the linked bodies. It follows that the
Eq. 3 and replacing the inverse particle mass wj joint points and with them the according bodies
with the inverse mass matrix M−1, the formula
j can not move away from each other. However,
to compute the vector of rigid body corrections the two bodies can freely rotate around the joint
for constraint C(p) is: position. In order to define a translational con-
straint that removes only two translational de-
[∆xT,∆ϑT,...,∆xT,∆ϑT]T = grees of freedom the constraint space has to be
i i n n (6)
−1 T −1 T−1 reduced to a plane. The plane is defined by the
−M JC JCM JC C(p), joint point of one of the two bodies and the plane
normal that is also attached to this body. It fol-
where the matrix M−1JT converts the mass lows that the joint points have to be projected
C
weighted displacement to a displacement in into the plane to measure the constraint viola-
−1 T−1
world space. The matrix JCM JC is the tion. Then, the correction is computed in the
mass matrix in constraint space. two-dimensionalconstraintspace. In a final step
The question now is how to parameterize the themassweightedcorrectiondisplacementmust
rotations of the rigid bodies to compute a Jaco- be projection back into the three-dimensional
bian JC that can be multiplied with the inverse space. Similarly, the constraint space of a joint
mass matrix M−1. Orientations can be param- constrainingonetranslationaldegreeoffreedom
j is defined by attaching a unit vector to one of the
eterized in different ways like Euler angles or
quaternions. Accordingly, all these parameteri- rigid bodies. Again, the constraint violation is
zations lead to a different Jacobian. Neverthe- measuredbyfirstprojectingthejointpointsonto
less, a solution can be found by starting at the the constraint space and then computing the dis-
velocity level. The relationship between angu- tance.
lar momentum and angular velocity is L = Iω. Constraints that remove only rotational de-
Taking the first order integration, with the as- grees of freedom can be defined analogously by
sumption that the axis of rotation stays constant constraining the orientation of the linked bodies.
during the time step, and rearranging the result
gives I−1Lh = ωh = ϑ, where h is the time 6 Collision Handling
step size. The vector ϑ represents a rotation of
|ϑ| about the axis ϑ/|ϑ| and is known as the ex- In this section we present our solution to handle
ponential map from R3 to S3 [27], where S3 is collisions, contact, and friction. In the follow-
the space of rotation quaternions . The paper of ing paragraphs we combine the terms collision
Grassia [27] also presents the derivation of the and contact under the term proximity and use
rotational part ∂R(ϑ)/∂ϑ of the Jacobian JC the special terms only where differences occur.
which is not repeated here.
Proximity Detection We perform only one
discrete proximity detection step to ensure a
no reviews yet
Please Login to review.
