1 ArtiSynth Overview 1.1 System structure 1.3 Basic packages

1.2 Simulation overview

An overview of the ArtiSynth simulation mechanism is given here. More details are given in [11] and in the related overview paper.

1.2.1 Model advancement

ArtiSynth simulation works by advancing models forward in time using a series of time steps, with a physics integrator updating positions, velocities and forces for all of the dynamic components over each step. The step size is usually fixed, but can vary between models and also over different time ranges, as described in more detail in the “Advancing Models in Time” section of the ArtiSynth Reference Manual.

Simulation happens under the control of a scheduler, which determines the next time to advance to and then advances the root model using its advance() method. In pseudo-code, this looks like:

⬇

double t0 = 0;

while (simulating) {

double t1 = getNextAdvanceTime (t0);

rootModel.advance (t0, t1);

t0 = t1;

}

Most often, the next advance time is determined by adding the value of the root model’s maxStepSize property to t0.

Figure 1.1: Typical model advancement sequence, from time $t_{0}$ to time $t_{1}$ .

Within its advance() method, the root model in turn calls the advance() methods for its collection of models, while also calling the apply() methods for the model agents. For each model, the advancement sequence typically looks like that shown in Figure 1.1, with any input probes and controllers called before the advance() method and any monitors and output probes called after. The preadvance() method is an additional method that a model can use to perform pre-step initialization before input probes and controllers are called.

By default, the update interval in Figure 1.1 equals the maxStepSize of the root model. However, will be less if required by the model’s own maxStepSize value, or a smaller (explicit) update interval for any of the output probes. Figure 1.1 also assumes that all of the model agents are explicitly associated with the model. Input probes and controllers that are not associated with the model will be applied before the preadvance() method, and monitors and output probes not associated with the model will be called after those that are. Again, more precise details on this are given in the “Advancing Models in Time” section of the ArtiSynth Reference Manual.

1.2.2 Full coordinate formulation

It is important to understand that ArtiSynth uses a full coordinate formulation, in which the position and velocity of each dynamic body is both represented and solved using full, or unconstrained, coordinates, with extra equations added to the system solve to account for constraints such as joints and contact. In particular, point positions and velocities are represented using 3-vectors, while the “positions” of RigidBody and Frame components are represented as the translation and orientation of their coordinate frames. Internally, these frames are stored as a 3-vector plus a quaternion, but they are most often presented to the application as the body’s pose, which is the $4\times 4$ homogeneous transformation ${\bf T}_{BW}$ that maps points from body to world coordinates:

${\bf T}_{BW}=\left(\begin{matrix}{\bf R}_{BW}&{\bf p}_{BW}\\ 0&1\end{matrix}\right).$

Here ${\bf R}_{BW}$ is a $3\times 3$ rotation matrix and ${\bf p}_{BW}$ is a 3-vector giving the location of the frame origin in world coordinates. The “velocity” of RigidBody and Frame components is the spatial velocity

$\hat{\bf v}_{BW}\equiv\left(\begin{matrix}{\bf v}_{BW}\\ \boldsymbol{\omega}_{BW}\end{matrix}\right)$

where ${\bf v}_{BW}$ and $\boldsymbol{\omega}_{BW}$ give the coordinate frame’s translational and angular velocity in world coordinates. See Sections A.2 and A.5 for more details.

In contrast, some other simulation systems, including OpenSim [7], use reduced coordinates, in which the system dynamics are formulated using a smaller set of coordinates (such as joint angles) that implicitly take the system’s constraints into account. For example, a two degree-of-freedom (DOF) double pendulum in ArtiSynth is described using 12 degrees of freedom (six for each of the two bodies), with 10 constraint equations added to reduce the net DOF to two, while in OpenSim, a double pendulum would be described using only its two joint angles.

Each methodology has its own advantages. Reduced formulations yield systems with fewer degrees of freedom and no constraint errors. On the other hand, full coordinates make it easier to combine and connect a wide range of components, including rigid bodies and FE models, and they also allow modeling of situations where constraints are compliant (i.e., soft) such that the bodies connected to them can assume arbitrary poses with respect to each other. The speed advantages of reduced coordinates are much less important in ArtiSynth, where the main performance bottleneck is usually the computational requirements of FE models.

1.2.3 Time stepping and integration physics

This section describes the mechanical simulation that takes place inside the advance() method of MechModel. An integrator is used to advance the model from time $t_{0}$ to time $t_{1}$ , corresponding to a time step given by $h\equiv t_{1}-t_{0}$ . Several integrators are available and may be selected via the model’s integrator property. The ones discussed here are the first order explicit integrator ConstrainedForwardEuler and the first-order implicit integrator ConstrainedBackwardEuler (which is the system default).

Let the positions, velocities, and forces associated with all the dynamic components be denoted by the composite vectors ${\bf q}$ , ${\bf u}$ , and ${\bf f}$ . In addition, let the composite mass matrix be given by ${\bf M}$ . Newton’s second law then gives

${\bf f}=\frac{d{\bf M}{\bf u}}{dt}={\bf M}\dot{\bf u}+\dot{\bf M}{\bf u},$

(1.1)

where the $\dot{\bf M}{\bf u}$ term accounts for various “fictitious” forces.

Each integration step involves solving for the velocities ${\bf u}^{k+1}$ at time step given the velocities and forces at step . One way to do this is to solve the expression

${\bf M}\,{\bf u}^{k+1}={\bf M}{\bf u}^{k}+h\bar{\bf f}$

(1.2)

for ${\bf u}^{k+1}$ , where is the step size and $\bar{\bf f}\equiv{\bf f}-\dot{\bf M}{\bf u}$ . Given the updated velocities ${\bf u}^{k+1}$ , one can determine $\dot{\bf q}^{k+1}$ from

$\dot{\bf q}^{k+1}={\bf Q}{\bf u}^{k+1},$

(1.3)

where ${\bf Q}$ accounts for situations (such as rigid bodies) where $\dot{\bf q}\neq{\bf u}$ , and then solve for the updated positions using

${\bf q}^{k+1}={\bf q}^{k}+h\dot{\bf q}^{k+1}.$

(1.4)

(1.2) and (1.4) by themselves comprise a simple symplectic Euler integrator.

In addition to forces, bilateral and unilateral constraints give rise to (generally) nonlinear constraint equations of the form

${\bf g}({\bf q})=0,\qquad{\bf n}({\bf q})\geq 0.$

(1.5)

Bilateral constraints are equality constraints that may include rigid body joints, FEM incompressibility, and point-surface constraints, while unilateral constraints are inequality constraints that include contact forces, friction and joint limits. Linearizing (1.5) leads to locally linear constraints on the velocity ${\bf u}$ of the form

${\bf G}({\bf q}){\bf u}={\bf d}_{g},\quad{\bf N}({\bf q}){\bf u}\geq{\bf d}_{n},$

(1.6)

where ${\bf G}$ and ${\bf N}$ are the derivative matrices of ${\bf g}({\bf q})$ . Because ${\bf G}$ and ${\bf N}$ are evaluated at step but generally vary in time, we can use their time derivatives to improve constraint satisfaction for ${\bf u}^{k+1}$ by modifying (1.6) to

${\bf G}{\bf u}^{k+1}={\bf d}_{g},\quad{\bf N}{\bf u}^{k+1}\geq{\bf d}_{n},$

(1.7)

where

${\bf d}_{g}\equiv-h\dot{\bf G}{\bf u}^{k},\quad{\bf d}_{n}\equiv-h\dot{\bf N}{% \bf u}^{k}.$

(1.8)

Constraints give rise to constraint forces (in the directions ${\bf G}({\bf q})^{T}$ and ${\bf N}({\bf q})^{T}$ ) which supplement the forces of (1.1) in order to enforce the constraint conditions. In addition, for unilateral constraints, we have a complementarity condition in which ${\bf N}{\bf u}>0$ implies no constraint force, and a constraint force implies ${\bf N}{\bf u}=0$ . Any given constraint usually involves only a few dynamic components and so ${\bf G}$ and ${\bf N}$ are generally sparse.

Adding constraints to the velocity solve (1.2) leads to a mixed linear complementarity problem (MLCP) of the form

	$\displaystyle\left(\begin{matrix}{\bf M}&-{\bf G}^{T}&-{\bf N}^{T}\\ {\bf G}&{\bf R}_{g}&0\\ {\bf N}&0&{\bf R}_{n}\end{matrix}\right)\left(\begin{matrix}{\bf u}^{k+1}\\ \tilde{\boldsymbol{\lambda}}\\ \tilde{\boldsymbol{\theta}}\end{matrix}\right)+\left(\begin{matrix}-{\bf M}{% \bf u}^{k}-h{\bf f}^{k}\\ -{\bf b}_{g}\\ -{\bf b}_{n}\end{matrix}\right)=\left(\begin{matrix}0\\ 0\\ {\bf w}\end{matrix}\right),$
	$\displaystyle 0\leq\boldsymbol{\theta}\perp{\bf w}\geq 0,$		(1.9)

where ${\bf w}$ is a slack variable, $\tilde{\boldsymbol{\lambda}}$ and $\tilde{\boldsymbol{\theta}}$ give the constraint force impulses over the time step, and the constraint offsets ${\bf b}_{g}$ and ${\bf b}_{n}$ are typically set to ${\bf d}_{g}$ and ${\bf d}_{n}$ as per (1.8). ${\bf R}_{g}$ and ${\bf R}_{n}$ are optional diagonal regularization terms, which if non-zero serve to soften the constraints, either to handle redundancy or to simulate physical compliance, as discussed later in Sections 3.4.9 and 8.6.1. If ${\bf R}_{g}\neq 0$ or ${\bf R}_{n}\neq 0$ , then ${\bf d}_{g}$ or ${\bf d}_{n}$ will typically also include a term related to the constraint error.

The actual constraint forces $\boldsymbol{\lambda}$ and $\boldsymbol{\theta}$ can be determined by dividing the constraint impulses by the time step :

$\boldsymbol{\lambda}=\tilde{\boldsymbol{\lambda}}/h,\quad\boldsymbol{\theta}=% \tilde{\boldsymbol{\theta}}/h.$

(1.10)

By itself, the MLCP (1.9) implements a constrained forward Euler integrator corresponding to setting the MechModel’s integrator property to ConstrainedForwardEuler.

For stiff systems, which often arise with systems containing FE models, explicit integrators can be unstable unless the step size is set to a very small value. This can typically be corrected by using an implicit integrator. A first order example of this is a backward Euler integrator in which ${\bf f}^{k}$ is replaced with an approximation of ${\bf f}^{k+1}$ given by the Taylor series expansion

${\bf f}^{k+1}\approx{\bf f}^{k}+\frac{\partial{\bf f}}{\partial{\bf q}}\delta{% \bf q}+\frac{\partial{\bf f}}{\partial{\bf u}}\delta{\bf u},$

(1.11)

where ${\partial{\bf f}}/{\partial{\bf q}}$ and ${\partial{\bf f}}/{\partial{\bf u}}$ are the negatives of the system stiffness and damping matrices, respectively. Since

$\delta{\bf u}={\bf u}^{k+1}-{\bf u}^{k}\quad\text{and}\quad\delta{\bf q}% \approx h{\bf u}^{k+1},$

(1.11) can be reformulated as

${\bf f}^{k+1}\approx{\bf f}^{k}+h\frac{\partial{\bf f}}{\partial{\bf q}}{\bf u% }^{k+1}+\frac{\partial{\bf f}}{\partial{\bf u}}({\bf u}^{k+1}-{\bf u}^{k}),$

(1.12)

which leads to the modified MLCP

	$\displaystyle\left(\begin{matrix}\hat{\bf M}&-{\bf G}^{T}&-{\bf N}^{T}\\ {\bf G}&{\bf R}_{g}&0\\ {\bf N}&0&{\bf R}_{n}\end{matrix}\right)\left(\begin{matrix}{\bf u}^{k+1}\\ \tilde{\boldsymbol{\lambda}}\\ \tilde{\boldsymbol{\theta}}\end{matrix}\right)+\left(\begin{matrix}-{\bf b}_{f% }\\ -{\bf b}_{g}\\ -{\bf b}_{n}\end{matrix}\right)=\left(\begin{matrix}0\\ 0\\ {\bf w}\end{matrix}\right),$
	$\displaystyle 0\leq\boldsymbol{\theta}\perp{\bf w}\geq 0,$		(1.13)

where

$\hat{\bf M}\equiv{\bf M}-h\frac{\partial{\bf f}}{\partial{\bf u}}-h^{2}\frac{% \partial{\bf f}}{\partial{\bf q}}\quad\text{and}\quad{\bf b}_{f}\equiv\left({% \bf M}-h\frac{\partial{\bf f}}{\partial{\bf u}}\right){\bf u}^{k}+h{\bf f}^{k}.$

The MLCP (1.13) corresponds to the default ArtiSynth integrator ConstrainedBackwardEuler. Higher order integrators, such as Newmark methods, often give rise to MLCPs with an equivalent form. Most ArtiSynth integrators use some variation of (1.13) to determine the system velocity at each time step; a description of the second order Trapezoidal integrator is given in [11].

To assemble (1.13), the MechModel component hierarchy is traversed and the methods of the different component types are queried for the required values. Dynamic components (type DynamicComponent) provide ${\bf q}$ , ${\bf u}$ , and ${\bf M}$ ; force effectors (ForceEffector) determine ${\bf b}_{f}$ and the stiffness/damping augmentation used to produce $\hat{\bf M}$ ; and constrainers (Constrainer) supply ${\bf G}$ , ${\bf N}$ , ${\bf b}_{g}$ and ${\bf b}_{n}$ .

1.2.4 Parametric and attached dynamic components

By default, a dynamic component is advanced through time in response to the forces applied to it. However, it is also possible to set a dynamic component to be parametric, such that it does not respond to forces and its position and/or velocity is either fixed or specified directly by the simulation. This is done by setting the component’s dynamic property to false, either through the GUI (by selecting the component(s) and choosing Edit properties ... from the context menu), or in code using the methods

void setDynamic (boolean enable)	Sets the dynamic property.
boolean isDynamic()	Queries the dynamic property.

Usage of setDynamic() is usually restricted to dynamic components that have mass, which include Particle and RigidBody but not Point or Frame.

Because the positions and velocities of parametric components are given, the MLCP (1.13) needs to be modified to reflect this. This involves partitioning the system velocity and force vectors into

${\bf u}=\left(\begin{matrix}{\bf u}_{\alpha}\\ {\bf u}_{\rho}\end{matrix}\right),\qquad{\bf f}=\left(\begin{matrix}{\bf f}_{% \alpha}\\ {\bf f}_{\rho}\end{matrix}\right),$

where $\alpha$ and $\rho$ denote the active (non-parametric) and parametric components, respectively. The MLCP is then likewise partitioned. Considering for simplicity the case with only bilateral constraints, this takes the form

$\left(\begin{matrix}\hat{\bf M}_{\alpha\alpha}&\hat{\bf M}_{\alpha\rho}&{\bf G% }_{\alpha}^{T}\\ \hat{\bf M}_{\rho\alpha}&\hat{\bf M}_{\rho\rho}&{\bf G}_{\rho}^{T}\\ {\bf G}_{\alpha}&{\bf G}_{\rho}&0\\ \end{matrix}\right)\left(\begin{matrix}{\bf u}_{\alpha}^{k+1}\\ {\bf u}_{\rho}^{k+1}\\ \tilde{\boldsymbol{\lambda}}\end{matrix}\right)=\left(\begin{matrix}{\bf b}_{f% \alpha}\\ {\bf b}_{f\rho}\\ {\bf b}_{g}\end{matrix}\right).$

Since ${\bf u}_{\rho}^{k+1}$ is given, the system can be reduced to

$\left(\begin{matrix}\hat{\bf M}_{\alpha\alpha}&{\bf G}_{\alpha}^{T}\\ {\bf G}_{\alpha}&{\bf R}_{g}\\ \end{matrix}\right)\left(\begin{matrix}{\bf u}_{\alpha}^{k+1}\\ \tilde{\boldsymbol{\lambda}}\end{matrix}\right)=\left(\begin{matrix}{\bf b}_{f% \alpha}-\hat{\bf M}_{\alpha\rho}{\bf u}_{\rho}^{k+1}\\ {\bf b}_{g}-{\bf G}_{\rho}{\bf u}_{\rho}^{k+1}\end{matrix}\right).$

ArtiSynth also provides an attachment mechanism that allows dynamic components to be directly attached to other dynamic components. This is one of the primary means by which different parts of a model are connected together: points and particles (including FE nodes) can be attached to rigid bodies or other particles, and Frame components can be attached to other rigid bodies or connected to the nodes of an FE model.

When a dynamic component is attached, its position becomes a direct function of the positions of the one or more master components to which it is attached. For example, if a component is attached to a master component , then its position ${\bf q}_{a}$ is a function of the position ${\bf q}_{m}$ of :

${\bf q}_{a}=g({\bf q}_{m}).$

Differentiating this yields a linear relationship between the velocities ${\bf u}_{a}$ and ${\bf u}_{m}$ of and ,

${\bf u}_{a}={\bf G}{\bf u}_{m},$

(1.14)

where ${\bf G}$ is the derivative of .

Attachments can be considered as constraints in which the positions/velocities of the attached components are constrained to be functions of the positions/velocities of the master components. While enforcement could be effected by adding velocity constraints of the form (1.14) to the system MLCP (1.13), in practice ArtiSynth actually uses the velocity constraints to preprocess the MLCP and remove the attached components from the system solve, which generally improves solve times. The attached positions/velocities are then updated after the solve based on the master positions/velocities.

As mentioned in Section 1.1.1, attachments are implemented by components which implement DynamicAttachment. These may be added to the model hierarchy separately, or, in the case of Marker and AttachedFrame components, embedded within specialized Point and Frame subclasses which are always attached.

Overall, a dynamic component can be in one of three states:

active: Component is dynamic and unattached. The method isActive() returns true. The component will move in response to forces.
parametric: Component is not dynamic, and is unattached. The method isParametric() returns true. The component will either remain fixed, or will move around in response to external inputs specifying the component’s position and/or velocity. One way to supply such inputs is to use controllers or input probes, as described in Section 5.
attached: Component is attached. The method isAttached() returns true, and the method getAttachment() will return the actual attachment component. The component will move so as to follow the other master component(s) to which it is attached.

1.1 System structure Bibliography 1.3 Basic packages

Generated on Sun Jun 28 14:20:46 2026 by LaTeXML [LOGO]