3. Overview

3.1. Platforms

AGX Dynamics is built and tested under the following platforms and compilers:

  • Windows 11 & 10; both 32 and 64 bit. Compilers: Visual studio 2017, 2019, 2022

  • Linux, Ubuntu 18.04 (gcc 7.3.0), Ubuntu 20.04 (gcc 9.3.0), Ubuntu 22.04 (gcc 11.2.0)

  • Mac OS X (clang 11.0.0)

3.2. Libraries/namespaces

AGX consists of more than 500 classes which are separated into a number of namespaces and libraries listed in Table 1. Each namespace is accessible through the include directive:

#include <namespace/headerfile.h>
Table 3.1 Available namespaces and libraries

NAMESPACE

LIBRARY NAME

DESCRIPTION

agxSabre

agxSabre.lib

Sparse matrix solver.

agx, agxData

agxCore.lib

Core dynamics library, including basic data types: Vec3, Quat, Matrices ref_ptr, etc.

agxCollide

agxPhysics.lib

Geometric intersection system for calculating intersection points, normals, penetration depth. Classes for various primitives such as Sphere, Cylinder, Box, …

agxIO

agxPhysics.lib

Contain classes for reading / (writing) various file formats (images, models).

agxSDK

agxPhysics.lib

Simulation framework. Higher level abstraction for modeling a system, stepping collision and dynamics forward in discrete time steps. Event handling, materials.

agxStream

agxPhysics.lib

Serialization classes.

agxUtil

agxPhysics.lib

Various utilities for reading CPU information etc.

agxNet

agxPhysics.lib

Classes for communication, Sockets, Compression.

agxOSG

agxOSG.lib

Various classes for integrating AGX to OpenSceneGraph. For testing/debug purposes.

agxModel

agxModel.lib

High level modeling primitives, Tree, Terrain etc.

agxWire

agxPhysics.lib

Implementation of lumped element wire.

agxTerrain

agxTerrain.lib

Implementation of deformable Terrain.

An application using physics simulation typically requires the following functionality:

Table 3.2 Functionality in different namespaces

NAMESPACE

FUNCTIONALITY

CLASS NAME

agxCollide

Modeling API for geometry/shapes

Geometry, Shape, Sphere, Box, …

agxCollide

Intersection testing and extraction of contact data

Space, SweepAndPruneBroadPhase

agx

Modeling API for physical bodies

RigidBody, MassProperties

agx

Modeling API for constraints

Hinge, BallJoint, LockJoint, …

agx

Solving non-linear systems for constraints/forces. Integration of velocities and positions.

ZorroSolver

agxSDK

Integration of user code for event handlers

StepEventLister, ContactEventListener, …

agxIO

Reading Wave front files, Collada, images…

MeshReader, ImageReaderPNG, …

agxSensor

Reading gamepads etc.

Joystick, JoystickManager..

When building/running an application with dynamics, it is important that collision space (agxCollide::Space) and the dynamics are updated and stepped correctly. Therefore there is a class that integrates these two into a higher level functionality: agxSDK::Simulation. This class makes sure that both collision space and the dynamics system are updated when objects are added or removed. Even though it is possible to use agxCollide::Space individually, we recommend always using an agxSDK::Simulation instance when building simulations.

3.3. Debug vs. Release build mode

In the installer for Microsoft Windows, the debug libraries are available. Library/runtime files in debug build has the letter ‘d’ added to the library name:

BUILD MODE

SAMPLE NAMES

FLAGS

Release

agxCore.lib, agxPhysics.dll, agxModel.lib

C++: -D_NDEBUG Linker: /MD

Debug

agxCored.lib, agxPhysicsd.dll,agxModel.lib

C++: -D_DEBUG Linker: /MDd

When building applications with Microsoft Visual Studio it is very important to never mix libraries built in debug and release mode. The reason for this is how memory allocation/STL containers work in Visual Studio. Mixing libraries breaks the ODR (One Definition Rule) (see more on the web).

3.4. Solver/Dynamics

The solver in AGX Dynamics is a hybrid solver which is very flexible in terms of specifying solve type for various constraints. The solver contains both an iterative solver and a direct solver.

3.4.1. Iterative solver

An iterative solver gradually improves the solution with the number of iterations, and after a finite number of iterations it thus delivers a solution with finite precision. The most common type of iterative solver is the Gauss-Seidel relaxational method and an extension called Projected Gauss-Seidel (PGS) which is used to solve constrained linear systems of equations. It solves the equations one-by-one and updates all the dependent variables of the system before proceeding with the next equation. One iteration corresponds to a full pass through all the equations. Since the constraint equations are treated one-by-one to solve for the constraint force and corresponding updated velocity, and the information about this local solution propagates to the coupled equations, it is often referred to as an iterative impulse propagation method. However, the term propagation is somewhat misleading since the ordering of the equations is not well defined in PGS, so information propagates randomly, unless specifically ordered. Relaxational solvers are moderately efficient and conceptually simple, at least when the local constraint violations are small, which is often true for systems with small mass ratios and weak couplings. Therefore they can be fast and efficient for large scale stacking problems with homogeneous types of bodies, while they are less efficient for e.g. lines and wires with large mass and force ratios, where they result in rubber band behavior unless the number of iterations is set to an extremely high value. They are also not so good at treating loop closures in physical systems, unless these problems are given special treatment.

3.4.2. Direct solver

A direct solver finds the exact solution (down to machine precision) in a limited number of operations (provided there exist a well-defined solution). For sparsely connected systems, a sparse type of solver is much preferred over a general dense solver, and can give enormous performance benefits. The chances of finding a proper physical solution are improved through regularization, which introduces an extra inverse mass/inertia term in all equations.

Our direct solver computes solutions of quadratic programming (QP) problems using a block pivot method. This is similar to Newton-Raphson iterations without damping or line search and is generally very fast and almost always returns a good solution in just a few steps. The linear algebra operations rely on our own sparse symmetric indefinite factorization code which is tailored for multi-body dynamics, and can also handle ill-conditioned systems of equations arising from joint singularities, also known as gimbal lock. The overall algorithm is generally faster than publicly available QP solvers based on active set or interior point methods.

The direct solver is very robust and can handle large multi-body systems with very high mass ratios, and yet maintain constraints accurately, to second order with respect to the time step. This is slower than the iterative solver but the latter makes larger errors and cannot process systems that have large mass ratios, or many joints. The simulation of a heavy vehicle for instance requires the direct solver for stability.

3.4.3. Split solver

This is the default solver in AGX Dynamics. It is a hybrid solver utilizing both the DIRECT and the ITERATIVE solver. By default, all ordinary constraints (hinges, prismatics, etc.) and all contacts normal forces are in a first stage solved with the direct solver. The normal forces are used as input to the iterative solver which solves for friction in the second stage. In the third and last pass, the constraints and normal forces are solved again in the direct solver given the frictional forces in the system. This results in a very fast solution, even for moderately large contact scenes with friction. For even larger systems with stacking etc., it is more efficient to select a pure iterative solution (described in Section 11.16).

Selected constraints can also explicitly solver solved in the iterative solver, (see Section 12.17). The DIRECT solver utilizes an in-house implementation of a sparse matrix solver: SABRE. SABRE is an optimized solver for rigid body simulations. It makes heavy use of BLAS-level3 which results in a fully utilized CPU pipeline.

The split solver is slower than the purely iterative one, but it produces much better results with residual errors that are at least two orders of magnitude smaller. This generally produces more stable stacking with much better dry friction in general. Because of a fundamental mathematical limitation however, the direct solver becomes slower and slower for piling or volume filling configurations. The slowdown - the rate of change of the execution time divided by the number of objects depends on the topology of the problem. It is worst for volume packing, moderate for planar constructions, i.e. walls, and nil for linear ones, i.e., stacks.

The execution time for the purely iterative solver grows only linearly with the number of bodies if the number of iteration is kept constant.

3.5. Multi-core support

Modern processors usually have two or more execution units, cores, which can be used to parallelize systems.

AGX is built upon a task model, which allows for a whole tree of tasks, each depending on other tasks. This task-tree can be split and run in parallel based on each tasks predecessors. For example, update bounding volumes can operate in parallel on each geometry, which makes it parallelizable. On the other hand, the broad phase tasks, require the bounding volume to be updated, so it need to wait until update bounding volume is done. The directory Components in the binary directory of AGX consists of many task files (both binary and XML) which together build the complete structure of the AGX Dynamics engine. More on parallelization in Section 31.5.

3.6. Time integration

The time stepper used in AGX is derived from a discrete variational principle applied to the discrete Lagrangian representing the constrained dynamic system. This results in a stepping equation that also is valid for non-smooth and dry frictional forces, which are represented through Rayleigh dissipation functions. In the presence of smooth forces only this stepping equation corresponds to the classical Leapfrog integration scheme. For constraint impacts, a second order term is used to correct the velocities and the corresponding impulses. The rotational degrees of freedom require special treatment since the integration of the angular momentum and the quaternions contain a non-linear coupling in the gyroscopic term. In many cases the gyroscopic term is neglected (always in game physics engines, and usually the approximation is not even known, and thought to be “exact”). Neglecting gyroscopic forces corresponds to conservation of angular velocity rather than angular momentum, which in turn can be said to correspond to an approximation where the inertia tensor is homogeneous, which is true for symmetric bodies (i.e. a sphere with homogeneous mass distribution).

The integration of rigid body motion in AGX is based on a Cartesian representation of body coordinates in world frame \(x\) and their velocities \(v = \dot{x}\) Their orientations are represented with quaternions \(q\) which represent the transformation from the body frame to world coordinates, as well as the angular velocity \(\omega\) in world frame. If we label \(f\) as the sum of external and constraint forces, and \(\tau\) as the sum of external and constraint torques, the time integration is

(3.1)\[ \begin{align}\begin{aligned}v_{k+1}=v_k + \frac{h}{m} f_k\\\omega_{k+1}= \omega_k + hJ_k^{-1} \tau_k\\x_{k+1}= x_k + hv_{k+1}\\q_{k+1}= q \ \exp(\frac{h}{2} \omega_{k+1})\end{aligned}\end{align} \]
where the quantities are defined according to
\(m\) the scalar mass
\(J\) the inertia tensor in world frame
\(h\) the timestep
\(k\) the discrete time
in addition, the exponential of the angular velocity is defined as the quaternion

(3.2)\[\begin{split}\exp(\cfrac{h}{2}\omega) = \begin{bmatrix} \cos(\cfrac{h}{2})\parallel\omega\parallel \\ \cfrac{1}{\parallel\omega\parallel} \sin\left(\cfrac{h}{2}\right)\omega \end{bmatrix}\end{split}\]

where quaternions are defined so that the first component is the real part, and the other three are the imaginary part. For a quaternion \(q\), we denote the real and imaginary parts as \(q_0\) and \(q_v\) respectively. We also use the right handed quaternion product so that for two quaternions \(p\), \(q\) we have

(3.3)\[\begin{split}pq = \begin{bmatrix} p_0 \\ p_v \end{bmatrix} \begin{bmatrix} q_0 \\ q_v \end{bmatrix} = \begin{bmatrix} p_0q_0 - p_v \cdot q_v \\ p_0q_v + q_0p_v + p_v × q_v \end{bmatrix}\end{split}\]

and (·) and (×) are the ordinary scalar and vector products, respectively. The latter is also right handed so that for two vectors \(x,y \ \epsilon \ \mathbb{R}^3\) then

(3.4)\[\begin{split}x×y= \begin{bmatrix} -x_3y_2 + x_2y_3 \\ x_3y_1 - x_1y_3 \\ -x_2y_1 + x_1y_2 \end{bmatrix}\end{split}\]

The time integration described in Eqn. (3.1) corresponds to the velocity version of the Störmer-Verlet stepper [1]. The equations used to solve for the constraint forces are described elsewhere[see 2, Eqn.4.31 p. 100], and this is a symmetrized variant of the Shake algorithm [1] which includes robust constraint stabilization as well as constraint relaxation to prevent ill-conditioning in the linear algebra. Contact constraints are also described in [2], and the contact forces are computed by solving a linear complementarity problem.

3.7. Contact mechanics (Elastic Foundation Model)

Contact mechanics in AGX Dynamics is governed by a dynamic set of local contact constraints. These constraints exist where contacts are detected and are constantly updated as a simulation progresses. Mechanically they represent a linear-elastic material constitution which is specified individually for every ContactMaterial. Furthermore, the constraints dynamically assess and utilize the local contact surface, penetration and material volume in order to collectively achieve highly stable and powerful contact mechanics.

Note however that this contact mechanics model is not intended nor generally suitable for systems that are heavily influenced by detailed or global structural deformations. Such demanding contact interactions may require complementing high-resolution nonlinear FEA.

The approach used in AGX Dynamics will lead to slight interpenetration between the bodies, which models local linear deformation. The amount of interpenetration will depend on material parameters such as Young’s Modulus, as well as external forces influencing the contact. For parameters influencing the contact mechanics, (see Material).

There are two different approaches for calculating the resulting stiffness between two interacting geometries/shapes: contact based and area based. Which version that is being used, can be controlled per ContactMaterial. See: ContactMaterial.

3.7.1. Area Based approach

Collision detection will provide points, which are a discretization of the contact patch. For each point, a separating normal, as well as a contact depth (in case of overlap) and a contact area (a part of the complete contact patch) will be computed. Each contact point can be considered as a spring. Assuming the ContactMaterial’s Young’s modulus \(E\) and a material rest length \(h\), the spring’s contact stiffness per unit area \(k\) can be defied as:

\[k = E / h\]

Note that the rest length \(h\) measures the total length of the bodies in contact direction, and it varies with the position of the contact point and the direction of the contact normal. The total length of the involved shapes belonging to each of the two bodies is measured along the point and normal, with an upper cut-off representing the maximum length of the locally elastic domain, as well as a lower cut-off for numerical reasons.

Given a depth \(d\) computed from collision detection, the local pressure \(P\) can then be computed as:

\[P = k \cdot d\]

Given the contact area A for the contact point, the force F for the contact point is computed as:

\[F = P \cdot A\]

This contact model is a simplified version of the Elastic Foundation Model, which is defined in the publications by Johnson [3] and Perez [4]. Note that the version presented here differs from the ones presented in the literature in the following ways:

  • The effect of Poisson’s ratio is ignored, thus corresponding to a Poisson’s ratio of 0 (a Poisson’s ratio of 0.3 would increase k with a factor of 1.35).

  • The Young’s Modulus is defined for each pair of materials, and the rest lengths of the two bodies are added, thus giving a common spring stiffness per unit area. This is a simplification of the original model, in which the stiffness per unit area is computed for each body given its Young’s modulus, rest length (andPoisson’s ratio), and a common spring stiffness per unit area is computed by \(k = \frac{k_1 \cdot k_2}{k_1 + _k2}\)

3.7.2. Contact point based approach (default)

The default contact model is an approximation of the area-based model above, where \(h\) is set to 1m and the contact area is approximated.

This model has the advantage that it is computationally faster and gives good-enough stability and fidelity in scenes with objects of sizes 1dm and larger as well as heavy masses, but it is less exact, especially when dealing with smaller and lighter objects.

3.8. Unicode

AGX uses utf-8 internally. One of the reasons is that it is backwards compatible with ascii. It will also save storage. std::string/agx::String works just as it is as a storage container for utf8 strings.

In <agx/Encoding.h> there are several utility methods for converting between Wide and utf8.

So for example, reading a string encoded in utf8 from the windows registry requires you to convert from utf8 to Wide:

const char *keyInUtf8; // Initialized with some utf8 string
std::wstring wKey = agx::utf8ToWide( key ); // Convert to Wide character string

agx::String& value; // Here is where we store the result

dwType=0;
DWORD dwLen=0;
status = RegQueryValueExW(hkey, wKey.c_str(), 0,&dwType,nullptr, &dwLen);

if (status== ERROR_SUCCESS &&  dwType == REG_SZ)
{
agx::Vector<wchar_t> data;

// Since we use wchar, this will overallocate, but that's ok.
data.resize(dwLen+1);

// Now get the value
status = RegQueryValueExW(hkey, wKey.c_str(), nullptr,nullptr,(PBYTE)(data.ptr()), &dwLen);
if (status== ERROR_SUCCESS)
{
  // If we use RegQueryValueEx and UNICODE is not set, then the API call
  // will be RegQueryValueExA that will convert from wstring to ANSI if
  // the data in the registry key is written as unicode.
  //
  // We don't want that conversion since it might not be safe and having to mess
  // around with code pages is not fun.

  value = agx::wideToUtf8( std::wstring(data.ptr() ) );
  found = true;
}
RegCloseKey(hkey);
}

3.8.1. Limitations

Applications such as agxViewer will not be able to read files in a path containing utf8 characters.

C++ streams do not handle unicode paths in windows, and there are no open methods that handle wstring. Hence basic_ios and classes such as ifstream must not be used for accessing files in paths with non-ascii characters.