# Time integration¶

Solfec-1.0 implements three variants of time integration of rigid body kinematics:

• RIG_POS – a semi–explicit scheme with a positive energy drift, cf. NEW1 in 
• RIG_NEG – a semi–explicit scheme with a negative energy drift, cf. NEW2 in 
• RIG_IMP – a semi–implicit without energy drift (but rather, oscillation), cf. NEW3 in 

And two variants of time integration for pseudo–rigid and finite–element kinematics:

• DEF_EXP – an explicit leap–frog scheme, cf. Section 5.1 in 
• DEF_LIM – a linearly implicit leap–frog scheme,  and TR1

A time integration scheme is selected by modifying the scheme parameter of the BODY object.

## Rigid integration¶

Linear motion is integrated like deformable motion. Rigid rotations are integrated as follows

(31)$\mathbf{\Lambda}^{t+\frac{h}{2}}=\mathbf{\Lambda}^{t}\exp\left[\frac{h}{2}\mathbf{\Omega}^{t}\right]$
(32)$\mathbf{T}^{t+\frac{h}{2}}=\left(\mathbf{\Lambda}^{t+\frac{h}{2}}\right)^{T}\mathbf{t}^{t+\frac{h}{2}}$
(33)$\mathbf{\Omega}^{t+\frac{h}{2}}=\mathbf{J}^{-1}\left[\exp\left[-\frac{h}{2}\mathbf{\Omega}^{t}\right]\mathbf{J}\mathbf{\Omega}^{t}+\frac{h}{2}\mathbf{T}^{t+\frac{h}{2}}\right]$
(34)$\mathbf{\Omega}_{1}^{t+h}=\mathbf{\Omega}^{t}+\mathbf{J}^{-1}h\left[\mathbf{T}^{t+\frac{h}{2}}-\mathbf{\Omega}^{t+\frac{h}{2}}\times\mathbf{J}\mathbf{\Omega}^{t+\frac{h}{2}}\right]$

If explicit

(35)$\mathbf{\Lambda}^{t+h}=\mathbf{\Lambda}^{t+\frac{h}{2}}\exp\left[\frac{h}{2}\mathbf{\Omega}_{1}^{t+h}\right]$
(36)$\mathbf{\Omega}_{2}^{t+h}=\mathbf{J}^{-1}\exp\left[-\frac{h}{2}\mathbf{\Omega}_{1}^{t+h}\right]\left[\exp\left [-\frac{h}{2}\mathbf{\Omega}^{t}\right]\mathbf{J}\mathbf{\Omega}^{t}+h\mathbf{T}^{t+\frac{h}{2}}\right]$

otherwise

(37)$\textbf{solve} \left(\exp\left[\frac{h}{2}\mathbf{\Omega}_{3}^{t+h}\right]\mathbf{J}\mathbf{\Omega}_{3}^{t+h}= \exp\left[-\frac{h}{2}\mathbf{\Omega}^{t}\right]\mathbf{J}\mathbf{\Omega}^{t}+h\mathbf{T}^{t+\frac{h}{2}}\right)$
(38)$\mathbf{\Lambda}^{t+h}=\mathbf{\Lambda}^{t+\frac{h}{2}}\exp\left[\frac{h}{2}\mathbf{\Omega}_{3}^{t+h}\right]$

The scheme ending at (35) is DEF_POS, ending at (36) is DEF_NEG, and using instead (37) and (38) is DEF_IMP. Above, $$\exp\left[\cdot\right]$$ is the exponential map defined by the Rodrigues formula

$\exp\left[\mathbf{\Psi}\right]=\mathbf{I}+\frac{\sin\left\Vert \mathbf{\Psi}\right\Vert }{\left\Vert \mathbf{\Psi}\right\Vert } \hat{\mathbf{\Psi}}+\frac{1-\cos\left\Vert \mathbf{\Psi}\right\Vert }{\left\Vert \mathbf{\Psi}\right\Vert ^{2}}\hat{\mathbf{\Psi}}^{2}$

where $$\mathbf{I}$$ is the $$3\times3$$ identity operator, $$\hat{\mathbf{\Psi}}$$ creates the skew symmetric matrix out of a 3-vector $$\mathbf{\Psi}$$, and $$\left\Vert \cdot\right\Vert$$ stands for the Euclidean norm. The time step is denoted as $$h$$.

## Deformable integration¶

$\mathbf{q}^{t+\frac{h}{2}}=\mathbf{q}^{t}+\frac{h}{2}\mathbf{u}^{t}$
$\mathbf{u}^{t+h}=\mathbf{u}^{t}+\mathbf{A}^{-1}h\mathbf{f}\left(\mathbf{q}^{t+\frac{h}{2}},\mathbf{u}^{t}\right)$
$\mathbf{q}^{t+h}=\mathbf{q}^{t+\frac{h}{2}}+\frac{h}{2}\mathbf{u}^{t+h}$

where in the explicit case

(39)$\mathbf{A}=\mathbf{M}\text{ for DEF_EXP}$

and in the linearly implicit case

(40)$\mathbf{A}=\mathbf{M}+\left(\frac{\eta h}{2}+\frac{h^{2}}{4}\right)\mathbf{K}\left(\mathbf{q}^{t+h/2}\right)\text{ for DEF_LIM}$

The time step is denoted as $$h$$. See TR1 for technical details.

## Implementation¶

Time integration is implement in bod.c (rigid, pseudo–rigid) and fem.c (finite–element) files. Inverse generalized inertia matrix $$\mathbf{A}^{-1}$$ is declared in bod.h as follows:

struct general_body
{
/* ... */

MX *inverse;      /* generalized inverse inertia oprator */

/* ... */
}


Rigid integration formulae (31)-(34) are in bod.c:BODY_Dynamic_Step_Begin.
Rigid integration formulae (35)-(38) are in bod.c:BODY_Dynamic_Step_End.
Pseudo–rigid integration is included in the same routines: first half–step and second half–step.
Finite–element, total Lagrangian formulation based, integration formulae (5) and (6) are in fem.c:TL_dynamic_step_begin.
Finite–element, total Lagrangian formulation based, integration formula (7) is in fem.c:TL_dynamic_step_end.

  (1, 2, 3) IJNME, 81(9):1073–1092, 2010.