Biomedical Engineering Reference
In-Depth Information
The presence of FE models means that the system is often stiff, requiring the
use of an implicit integration step where an approximation of
f
k
+
1
is used in place
of
f
k
. This can be achieved by replacing
M
and
f
k
M
and
f
k
, which
in (
11.1
) with
contain additional terms derived from the force Jacobians
u
[
34
].
Combining all this into a matrix form with the constraint conditions leads to a mixed
linear complementarity problem (MLCP), which we solve at each time step:
⊡
∂
f
/∂
q
and
∂
f
/∂
⊤
⊡
⊤
⊡
⊤
⊡
⊤
M
h
f
k
N
T
G
00
N
00
G
T
u
k
+
1
z
Mu
k
0
0
w
−
−
+
⊣
⊦
=
⊣
⊦
+
⊣
⊦
=
⊣
⊦
g
n
0
∈
z
√
w
ₒ
0
.
(11.2)
Nu
k
+
1
, and the
∗
Here
g
and
n
arise from the time derivatives of
G
and
N
,
w
∈
√
ₒ
complementarity condition 0
z
w
0 arises from the fact that for unilateral
0 and
Nu
k
+
1
constraints,
z
0 must be mutually exclusive. The system (
11.2
)
is also applicable to higher order integrators such as the trapezoidal rule [
34
], and
is also used to compute position corrections
>
>
q
that remove errors due to constraint
δ
drift and contact interpenetration.
11.3.2 Coupling FE Models and Rigid Bodies
In models such as our orofacial model, it is necessary to connect FE models to other
FE models and rigid bodies. In ArtiSynth, connecting FE and rigid-body is done
using point-based attachments, whereby an FE node is attached either to another FE
node, an FE element, or a rigid body. In all cases, this results in the state (position and
velocity) of the attached node becoming an explicit function of the states of several
other nodes or bodies. If we let
α
denote all
unattached (or
master
) nodes and bodies, and denote these sets' respective velocities
by
u
β
and
u
α
, then at any time
u
β
can be determined by the velocity constraint
β
denote the set of all attached nodes, and
u
β
+
G
β α
u
α
=
0
where
G
β α
is time varying and sparse. In other words, attachments can be imple-
mented as a special kind of bilateral constraint. If we partition system (
11.2
) into the
sets
h
f
, and ignore unilateral constraints for simplicity, we
Mu
k
β
and
α
,let
b
∗
+
obtain
⊡
⊤
⊡
⊤
⊡
⊤
M
αα
M
α β
−
G
T
G
T
β α
u
k
+
1
u
k
+
1
β
λ
α
λ
β
αα
−
b
b
⊣
⊦
⊣
⊦
=
⊣
⊦
.
M
β α
M
ββ
−
G
T
α β
−
I
β
(11.3)
g
G
αα
G
α β
0
0
g
G
β α
I
0
0
β
