Appendix A - Computer Graphics and Multimedia

```R. Tomcin & D. Sibbing & L. Kobbelt / Efficient Enforcement of Hard Articulation Constraints
Appendix A: Velocity Verlet Time Integration
rections)
Let the position of a rigid body be x and the velocity v with
respect to its center of gravity, the orientation q (represented
by a unit quaternion) and the angular velocity ω . Let m be
its mass and I be its orientation dependent inertia tensor. The
sum of all external forces and torques acting on the body due
to gravity and springs are denoted f and τ, respectively.
Quantities at the beginning t = t (l) and end t (l+1) =
(l)
t + h of a time step of length h are denoted by ∗(l) or
∗(l+1) , respectively. To achieve more accurate results, we
slightly relaxed the restrictions mentioned in Section 3 and
show how to replace the forward Euler time integration of
[GBF03] with a more accurate Velocity Verlet method.
ω) represent an incremental rotation and ◦ be
Let q (4ω
the quaternion multiplication. The linear and angular accelerations are obtained from the equations of motion σ =
ω × Iω
ω)) and a = mf . The Verlet method is given
I−1 (ττ − (ω
by first updating the body’s location with
h2 (l)
a
2
(9)
h2 (l)
σ ) ◦ q (l)
2
(10)
x (l+1) = x (l) + h v (l) +
q (l+1) = q (h ω (l) +
and then the velocities with
h (l)
a + a (l+1)
2
h (l)
(l+1)
(l)
σ + σ (l+1) .
ω
=ω +
2
v (l+1) = v (l) +
(11)
(12)
We reformulated this update to meet the requirements of
the implicit contact- and collision distinction explained in
[GBF03] and Section 3. In particular, all external forces
should be applied in the linear and angular velocity update
and not appear in the position and orientation updates.
The restructured Verlet time integration for a single body
is obtained by reformulating equations (9) to (12) in terms
of
h
v− = v − a ,
2
h
v+ = v + a
2
and
(l)
x (l+1) = x (l) + h v +
(l+1)
q
(l)
h ω+ −
=q
(13)
3. Update velocities at the end of the time step
(l+1)
= v+
(l+1)
= ω+ +
v−
ω−
(l)
(l)
h
ω+ = ω + τ ,
2
which leads to the following 3 step procedure:
1. Update velocities by applying external impulses
(l)
(l)
v + = v − + h a (l)
(l)
(l)
ω + = ω − + h I(l)
−1 (l)
τ
2. Update body locations (done frequently during joint corc 2014 The Author(s)
c 2014 The Eurographics Association and John Wiley & Sons Ltd.
Computer Graphics Forum h ∗ (l)
(l+1)
σ
+ σ∗
2
The angular acceleration σ ∗ does not include torque directly,
meaning that
h
h
∗
σ = −ω
ω × Iω
ω = − ω+ − τ × I ω+ − τ
Iσ
2
2
This introduces a velocity update at the end of each time
step in addition to the application of external impulses (see
(l+1)
Figure 2). Note that σ ∗
in the velocity update is not
known and needs to be determined with Newton’s method
in a similar way as the original Verlet integration described
in [KE04].
We neglect torque-dependent terms in the orientation up(l)
date (we use ω+ instead of ω(l) in Equation (13)), because
of a more reliable convergence of joint corrections for cases
involving strong rotational motors. This way, up to 6 times
larger simulation step sizes could be chosen for the trike
and Strandbeest examples. At the same time, we still enjoy some of the benefits of the accuracy of Verlet integration such as long-term energy conservation and momentum
conservation. Although not as accurate as the Moser Veselov
integrator [MV91] that conserves both energy and momentum completely as used by [SKV∗ 12], our simulator is more
general in the sense that it features different coefficients of
restitution and sustained drift-free articulation and contact
constraints.
Appendix B: A heuristic for splitting the kinematic graph
To show, that our splitting strategy yields linear timecomplexity, we define the bandwidth β(A) as in [BBK05]:
β(A) := max {βi (A)} ,
h
ω− = ω − τ ,
2
!
h2 (l) −1 (l)
(l) (l)
ω ×I ω
◦ q (l)
I
2
1≤i≤n
βi (A) := i − min { j|Ai j 6= 0}
1≤ j≤i
Let the number of original constraints attached to each of
the subparts k be zk . The zk can be freely chosen, as long
as they add up to the total number of constraints of the
original body. The rows corresponding to extra joints have
the largest bandwidth (bottom row in Figure 4). Keeping in
mind that extra joints constrain 6 degrees of freedom this
means β(A) ≤ 11 + 2 · zmax with zmax = max1≤k≤nP (zk ).
Both the factorization and the solution has therefore linear
time-complexity. Note, that this result can easily be extended
to all multibody systems with a tree-like kinematic graph.
R. Tomcin & D. Sibbing & L. Kobbelt / Efficient Enforcement of Hard Articulation Constraints
In practice, it is desirable to find a splitting strategy that maximizes performance. We found, that in contact-heavy scenarios more computation time of our simulator is spent on
solving the system than factorizing it. Therefore the workload is roughly proportional to the number of off-diagonal
matrix entries nM . Considering again the different kinds of
matrix entries (sorted by color as in Figure 4 (bottom row)),
the number of off-diagonal matrix entries can be written for
nP ≥ 2 as
nP
}|
green
red
{
z }| { z }| {
nM = ∑ zk (zk − 1) + 30(nP − 1) + 72(nP − 2)
k=1
nP −1
+ 12(z0 + znP ) + 24
∑
zk
(14)
k=2
|
{z
}
blue
We allow breaking apart the constraints between a pair of
bodies. For example, joints 3 and 4 in Figure 4 (bottom row)
might be attached to the same pair of bodies in the original
configuration. For a joint with 6 constraints that is cut in half,
this causes a maximum of 3 · 3 + 3 · 3 = 18 additional matrix
entries. Note, that although this introduces a small closed
kinematic loop, there is no fill-in during the factorization.
For this reason we neglect this additional term and propose
a heuristic to minimize nM .
The number of matrix entries (14) can be written as
!
nP −1
nM =
∑
z2k + 23zk
+ z20 + z2nP + 11(z0 + znP )
k=2
+ 102 nP − 174 .
Minimizing this function is an integer programming problem. Except for the obvious difference between inner parts
and the two end parts, the zk should be chosen as equal as
possible for minimal nM . Therefore we constrain the number of inner (and outer) constraints per part to only differ by
1 at most. nM is reformulated in terms of the total number of
P
constraints n = ∑nk=1
zk , the number of constraints of inner
nP −1
parts m = ∑k=2 zk and the number of inner parts p = nP −2.
With the ceil and floor functions d·e and b·c, let a inner parts
contain zlower = b mp c, and the others zlower + 1 constraints.
We then have
nP −1
∑
z2k
1 1
+
m2 + 12m − mn + 102p +γ
p 2
{z
}
f (n,m,p)
p
with γ ∈ [0, d e + 1]
4
A splitting strategy can be obtained by minimizing f for a
given n with respect to m and p, yielding
We found, that this strategy was optimal for n ∈ [38, 218].
For n > 78, we noticed that 16 constraints for the 2 outer
parts and 10 for the inner parts seems to always be optimal. To get the maximum number of inner parts with 10
m
m
constraints, p needs to be b 10
c or b 10
c + 1.
Our heuristic to find p and m that minimize nM can be
given by the following algorithm that takes n as input:
if n < 27 // no split
p = -1
m = 0
else if n < 37 // one split
p = 0
m = 0
else if n < 219 // n < 79 should work, too
p = n==37 ? 1 : round((n-12)/sqrt{102}-2)
m = (n-12)*p/(p+2)
else
m = n-32
p = floor(m/10)
if p + 1 - m mod (p+1) < m mod p
p = p + 1
end if
We verified the optimality of this heuristic empirically for
n < 108 . A thorough proof, however, might be exaggerated
at this point, since nM is not a very accurate performance
measure anyway.
Appendix C: An alternative to matrix modification for a
single contact
There exist two closely related versions of the Cholesky decomposition of a symmetric positive definite matrix A ∈
RnC ×nC . The first is defined as
z0 + zn P = n − m
z20 + z2nP = d
1
nM = n2 + 11n +
2
|
n − 12
p∗ = √
−2
102
(n − 12)p∗
.
m∗ =
p∗ + 2
black
z
f (m, p) ≤ nM :
(n − m)2
e
2
A = GGT ,
= a z2lower + (p − a)(zlower + 1)2
k=2
m
m
= p(b c + 1)2 − (p − m mod p)(2b c + 1) .
|
{z
}
p
p
a
nM can be written in a form that yields a lower bound
(15)
where G is a lower triangular matrix and the second is given
by
A = LDLT .
(16)
L is also a lower triangular matrix, but has ones on the diagonal and D is a diagonal matrix.
c 2014 The Author(s)
c 2014 The Eurographics Association and John Wiley & Sons Ltd.
Computer Graphics Forum R. Tomcin & D. Sibbing & L. Kobbelt / Efficient Enforcement of Hard Articulation Constraints
We show how to add na asymmetric rows and columns to
the matrix A (as for the case of Coulomb’s model for dynamic friction) based on an existing Cholesky decomposition of A in the form (15) or (16). For our cases, na will not
be greater than three, because we only insert one (frictional)
contact at a time.
Since we add rows/columns to the bottom/right of A, the
decomposition does not need to be changed explicitly. Instead, this change can be emulated during the solution process. All this is done in a block-based manner, handling all
na constraints at once.
For our derivations, we denote modified parts of the system
matrix and the decomposition with ∗¯ . The additional na constraints yield the new system matrix
¯ 1,2 A
¯ = A
A
T
¯ 2,1 A
¯ 2,2 .
A
¯ 1,2 ∈
¯ 2,1 ∈ RnC ×na , A
The dimensions of the new parts are A
nC ×na
na ×na
¯
R
and A2,2 ∈ R
.
First version with A = GGT
T
¯ 1,2 G
¯ 1,2
A
A
0
U
G
=
T
T
¯ 2,2
¯ 2,1 A
¯ 2,2
¯ 2,1 1na
0
D
A
G
{z
}
|
{z
} |
{z
}|
¯
A
(17)
¯ T2,1 U
¯ 1,2 + D
¯ 2,2 = A
¯ 2,2
G
⇒
¯ 1,2
U
¯ 2,1
G
(19)
⇒
¯ 2,2
D
(20)
⇒
1. compute y 1by forward substitution
of Gyy1 = b 1
T
¯ −1
¯
2. x 2 ← D
b
−
G
y
2
2,1 1
2,2
{z
}
|
y2
¯ 1,2 x 2
3. compute x 1 by back substitution of GT x 1 = y 1 − U
Second version of decomposition with A = LDLT
In case the decomposition exists in the form A = LDLT ,
we can solve the modified system in a similar manner. The
modified decomposition is given by
T
¯ 1,2 L
¯ 1,2
0
A
A
D
0
L
U
=
T
T
¯ 2,2
¯ 2,1 1na
¯ 2,1 A
¯ 2,2
0 D
0
1 na
L
A
¯ 1,2 , L¯ 2,1 and D
¯ 2,2 can be obtained
This time, the matrices U
analogously from
¯ T2,1 DU
¯ 1,2 + D
¯ 2,2
L
¯ 2,2
=A
⇒
¯ 1,2
U
¯ 2,1
L
⇒
¯ 2,2
D
⇒
¯ 2,2 ∈ Rna ×na is not a diagonal matrix
Note, that the block D
and asymmetric in general, since this is actually a blockbased LU decomposition for asymmetric matrices, mixed
with a usual Cholesky decomposition. Because the addi¯ the
tional rows/columns are added to the bottom/right of A,
factor G does not change. For a purely symmetric modification method that is able to place constraints at arbitrary
positions in the matrix, see [DH05].
¯ 1,2 , A
¯ 2,1 and A
¯ 2,2 can be
Matrix equations for the blocks A
¯ 1,2 , G
¯ 2,1 and D
¯ 2,2 once for
isolated from (17) to compute U
every contact:
¯ 1,2 = A
¯ 1,2
GU
¯
¯
GG2,1 = A2,1
More precisely, the overall solution procedure to obtain x 1
and x 2 consists of the following steps:
¯ 1,2 = A
¯ 1,2
LDU
¯ 2,1 = A
¯ 2,1
LDL
¯
U
¯
G
and then a back substitution of
T
¯ 1,2
x1
y
G
U
= 1 .
¯ 2,2
x2
y2
0
D
(18)
¯ 1,2 and
Because G is a lower triangular matrix, computing U
¯ 2,1 from (18) and (19) amounts to na forward substitutions
G
¯ 2,2 in (20) involves a sparse muleach. The computation of D
T ¯
¯
tiplication G2,1 U1,2 and a subtraction.
The original, unmodified linear system Axx = b can be solved
with a forward substitution Gyy = b , followed by a back substitution GT x = y . To solve the modified system, we need to
first do a forward substitution of
G
0
y1
b1
=
¯ T2,1 1na
y2
b2
G
c 2014 The Author(s)
c 2014 The Eurographics Association and John Wiley & Sons Ltd.
Computer Graphics Forum 1. compute y 1by forward substitution
of Lyy1 = b 1
T
¯ −1
¯
2. x 2 ← D
b
−
L
y
2
2,1 1
2,2
|
{z
}
y2
¯ 1,2 x 2
3. compute x 1 by back subst. of LT x 1 = D−1 y 1 − U
Permutation
Due to fill-in minimization, the given factors G (or L and
D, respectively) are based on a reordered version A of the
original matrix A∗ , that can be denoted in terms of a permutation matrix P:
A = PA∗ PT
¯ 1,2 and A
¯ 2,1
For this reason, the additional matrix parts A
must be obtained by permuting their original counter-parts
¯ 1,2 , L
¯ 2,1 and D
¯ 2,2 :
before computing U
¯ 1,2 = PA
¯ ∗1,2
A
¯ 2,1 = PA
¯ ∗2,1
A
The given right-hand side b ∗1 is permuted accordingly with
b 1 = Pbb∗1 . To obtain the upper part x ∗1 of the final solution,
the inverse permutation is used:
x ∗1 = PT x 1
Of course, the reordering does not need to be performed explicitly. Instead, it is achieved by accessing the elements in
the right order, which amounts to a simple book-keeping.
```