ArtiSynth Modeling Guide

John Lloyd and Antonio Sánchez

Last update: July, 2024

Contents

Preface

This guide describes how to create mechanical and biomechanical models in ArtiSynth using its Java API. Detailed information on how to use the ArtiSynth GUI for model visualization, navigation and simulation control is given in the ArtiSynth User Interface Guide. It is also possible to interface ArtiSynth with, or run it under, MATLAB. For information on this, see the guide Interfacing ArtiSynth to MATLAB.

Information on how to install and configure ArtiSynth is given in the installation guides for Windows, MacOS, and Linux.

It is assumed that the reader is familiar with basic Java programming, including variable assignment, control flow, exceptions, functions and methods, object construction, inheritance, and method overloading. Some familiarity with the basic I/O classes defined in java.io.*, including input and output streams and the specification of file paths using File, as well as the collection classes ArrayList and LinkedList defined in java.util.*, is also assumed.

How to read this guide

Section 1 offers a general overview of ArtiSynth’s software design, and briefly describes the algorithms used for physical simulation (Section 1.2). The latter section may be skipped on first reading. A more comprehensive overview paper is available online.

The remainder of the manual gives details instructions on how to build various types of mechanical and biomechanical models. Sections 3 and 4 give detailed information about building general mechanical models, involving particles, springs, rigid bodies, joints, constraints, and contact. Section 5 describes how to add control panels, controllers, and input and output data streams to a simulation. Section 6 describes how to incorporate finite element models. The required mathematics is reviewed in Section A.

If time permits, the reader will profit from a top-to-bottom read. However, this may not always be necessary. Many of the sections contain detailed examples, all of which are available in the package artisynth.demos.tutorial and which may be run from ArtiSynth using Models > All demos > tutorials. More experienced readers may wish to find an appropriate example and then work backwards into the text and preceding sections for any needed explanatory detail.

Chapter 1 ArtiSynth Overview

ArtiSynth is an open-source, Java-based system for creating and simulating mechanical and biomechanical models, with specific capabilities for the combined simulation of rigid and deformable bodies, together with contact and constraints. It is presently directed at application domains in biomechanics, medicine, physiology, and dentistry, but it can also be applied to other areas such as traditional mechanical simulation, ergonomic design, and graphical and visual effects.

1.1 System structure

An ArtiSynth model is composed of a hierarchy of models and model components which are implemented by various Java classes. These may include sub-models (including finite element models), particles, rigid bodies, springs, connectors, and constraints. The component hierarchy may be in turn connected to various agent components, such as control panels, controllers and monitors, and input and output data streams (i.e., probes), which have the ability to control and record the simulation as it advances in time. Agents are presented in more detail in Section 5.

The models and agents are collected together within a top-level component known as a root model. Simulation proceeds under the control of a scheduler, which advances the models through time using a physics simulator. A rich graphical user interface (GUI) allows users to view and edit the model hierarchy, modify component properties, and edit and temporally arrange the input and output probes using a timeline display.

1.1.1 Model components

Every ArtiSynth component is an instance of ModelComponent. When connected to the hierarchy, it is assigned a unique number relative to its parent; the parent and number can be obtained using the methods getParent() and getNumber(), respectively. Components may also be assigned a name (using setName()) which is then returned using getName().

A component’s number is not the same as its index. The index gives the component’s sequential list position within the parent, and is always in the range 0\ldots n-1, where n is the parent’s number of child components. While indices and numbers frequently are the same, they sometimes are not. For example, a component’s number is guaranteed to remain unchanged as long as it remains attached to its parent; this is different from its index, which will change if any preceding components are removed from the parent. For example, if we have a set of components with numbers

0 1 2 3 4 5

and components 2 and 4 are then removed, the remaining components will have numbers

0 1 3 5

whereas the indices will be 0 1 2 3. This consistency of numbers is why they are used to identify components.

A sub-interface of ModelComponent includes CompositeComponent, which contains child components. A ComponentList is a CompositeComponent which simply contains a list of other components (such as particles, rigid bodies, sub-models, etc.).

Components which contain state information (such as position and velocity) should extend HasState, which provides the methods getState() and setState() for saving and restoring state.

A Model is a sub-interface of CompositeComponent and HasState that contains the notion of advancing through time and which implements this with the methods initialize(t0) and advance(t0, t1, flags), as discussed further in Section 1.1.4. The most common instance of Model used in ArtiSynth is MechModel (Section 1.1.5), which is the top-level container for a mechanical or biomechanical model.

1.1.2 The RootModel

The top-level component in the hierarchy is the root model, which is a subclass of RootModel and which contains a list of models along with lists of agents used to control and interact with these models. The component lists in RootModel include:

models top-level models of the component hierarchy
inputProbes input data streams for controlling the simulation
controllers functions for controlling the simulation
monitors functions for observing the simulation
outputProbes output data streams for observing the simulation

Each agent may be associated with a specific top-level model.

1.1.3 Component path names

The names and/or numbers of a component and its ancestors can be used to form a component path name. This path has a construction analogous to Unix file path names, with the ’/’ character acting as a separator. Absolute paths start with ’/’, which indicates the root model. Relative paths omit the leading ’/’ and can begin lower down in the hierarchy. A typical path name might be

  /models/JawHyoidModel/axialSprings/lad

For nameless components in the path, their numbers can be used instead. Numbers can also be used for components that have names. Hence the path above could also be represented using only numbers, as in

  /0/0/1/5

although this would most likely appear only in machine-generated output.

1.1.4 Model advancement

ArtiSynth simulation proceeds by advancing all of the root model’s top-level models through a sequence of time steps. Every time step is achieved by calling each model’s advance() method:

  public StepAdjustment advance (double t0, double t1) {
     ... perform simulation ...
  }

This method advances the model from time t0 to time t1, performing whatever physical simulation is required (see Section 1.2). The method may optionally return a StepAdjustment indicating that the step size (t1 - t0) was too large and that the advance should be redone with a smaller step size.

The root model has it’s own advance(), which in turn calls the advance method for all of the top-level models, in sequence. The advance of each model is surrounded by the application of whatever agents are associated with that model. This is done by calling the agent’s apply() method:

   model.preadvance (t0, t1);
   for (each input probe p) {
      p.apply (t1);
   }
   for (each controller c) {
      c.apply (t0, t1);
   }
   model.advance (t0, t1);
   for (each monitor m) {
      m.apply (t0, t1);
   }
   for (each output probe p) {
      p.apply (t1);
   }

Agents not associated with a specific model are applied before (or after) the advance of all other models.

More precise details about model advancement are given in the ArtiSynth Reference Manual.

1.1.5 MechModel

Most ArtiSynth applications contain a single top-level model which is an instance of MechModel. This is aCompositeComponent that may (recursively) contain an arbitrary number of mechanical components, including finite element models, other MechModels, particles, rigid bodies, constraints, attachments, and various force effectors. The MechModel advance() method invokes a physics simulator that advances these components forward in time (Section 1.2).

For convenience each MechModel contains a number of predefined containers for different component types, including:

particles 3 DOF particles
points other 3 DOF points
rigidBodies 6 DOF rigid bodies
frames other 6 DOF frames
axialSprings point-to-point springs
connectors joint-type connectors between bodies
constrainers general constraints
forceEffectors general force-effectors
attachments attachments between dynamic components
renderables renderable components (for visualization only)

Each of these is a child component of MechModel and is implemented as a ComponentList. Special methods are provided for adding and removing items from them. However, applications are not required to use these containers, and may instead create any component containment structure that is appropriate. If not used, the containers will simply remain empty.

1.2 Physics simulation

Only a brief summary of ArtiSynth physics simulation is described here. Full details are given in [11] and in the related overview paper.

For purposes of physics simulation, the components of a MechModel are grouped as follows:

Dynamic components


Components, such as a particles and rigid bodies, that contain position and velocity state, as well as mass. All dynamic components are instances of the Java interface DynamicComponent.

Force effectors


Components, such as springs or finite elements, that exert forces between dynamic components. All force effectors are instances of the Java interface ForceEffector.

Constrainers


Components that enforce constraints between dynamic components. All constrainers are instances of the Java interface Constrainer.

Attachments


Attachments between dynamic components. While technically these are constraints, they are implemented using a different approach. All attachment components are instances of DynamicAttachment.

The positions, velocities, and forces associated with all the dynamic components are denoted by the composite vectors {\bf q}, {\bf u}, and {\bf f}. In addition, the composite mass matrix is 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} accounts for various “fictitious” forces.

Each integration step involves solving for the velocities {\bf u}^{k+1} at time step k+1 given the velocities and forces at step k. 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 h 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 (like 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) together comprise a simple symplectic Euler integrator.

In addition to forces, bilateral and unilateral constraints give rise to locally linear constraints on {\bf u} of the form

{\bf G}({\bf q}){\bf u}=0,\qquad{\bf N}({\bf q}){\bf u}\geq 0. (1.5)

Bilateral constraints may include rigid body joints, FEM incompressibility, and point-surface constraints, while unilateral constraints include contact and joint limits. 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}\hat{\bf M}^{k}&-{\bf G}^{T}&-{\bf N}^{T}\\
{\bf G}&0&0\\
{\bf N}&0&0\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\hat{\bf f}^{k}\\
-{\bf g}\\
-{\bf 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.6)

where {\bf w} is a slack variable, \tilde{\boldsymbol{\lambda}} and \tilde{\boldsymbol{\theta}} give the force constraint impulses over the time step, and {\bf g} and {\bf n} are derivative terms defined by

{\bf g}\equiv-h\dot{\bf G}{\bf u}^{k},\quad{\bf n}\equiv-h\dot{\bf N}{\bf u}^{%
k}, (1.7)

to account for time variations in {\bf G} and {\bf N}. In addition, \hat{\bf M} and \hat{\bf f} are {\bf M} and \bar{\bf f} augmented with stiffness and damping terms terms to accommodate implicit integration, which is often required for problems involving deformable bodies. The actual constraint forces \boldsymbol{\lambda} and \boldsymbol{\theta} can be determined by dividing the impulses by the time step h:

\boldsymbol{\lambda}=\tilde{\boldsymbol{\lambda}}/h,\quad\boldsymbol{\theta}=%
\tilde{\boldsymbol{\theta}}/h. (1.8)

We note here that ArtiSynth uses a full coordinate formulation, in which the position of each dynamic body is solved using full, or unconstrained, coordinates, with constraint relationships acting to restrict these coordinates. 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. 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 FEM models.

Attachments between components can be implemented by constraining the velocities of the attached components using special constraints of the form

{\bf u}_{j}=-{\bf G}_{j\alpha}{\bf u}_{\alpha} (1.9)

where {\bf u}_{j} and {\bf u}_{\alpha} denote the velocities of the attached and non-attached components. The constraint matrix {\bf G}_{j\alpha} is sparse, with a non-zero block entry for each master component to which the attached component is connected. The simplest case involves attaching a point j to another point k, with the simple velocity relationship

{\bf u}_{j}={\bf u}_{k} (1.10)

That means that {\bf G}_{j\alpha} has a single entry of -{\bf I} (where {\bf I} is the 3\times 3 identity matrix) in the k-th block column. Another common case involves connecting a point j to a rigid frame k. The velocity relationship for this is

{\bf u}_{j}={\bf u}_{k}-{\bf l}_{j}\times\boldsymbol{\omega}_{k} (1.11)

where {\bf u}_{k} and \boldsymbol{\omega}_{k} are the translational and rotational velocity of the frame and l_{j} is the location of the point relative to the frame’s origin (as seen in world coordinates). The corresponding {\bf G}_{j\alpha} contains a single 3\times 6 block entry of the form

\left(\begin{matrix}{\bf I}&[l_{j}]\end{matrix}\right) (1.12)

in the k-th block column, where

[l]\equiv\left(\begin{matrix}0&-l_{z}&l_{y}\\
l_{z}&0&-l_{x}\\
-l_{y}&l_{x}&0\end{matrix}\right) (1.13)

is a skew-symmetric cross product matrix. The attachment constraints {\bf G}_{j\alpha} could be added directly to (1.6), but their special form allows us to explicitly solve for {\bf u}_{j}, and hence reduce the size of (1.6), by factoring out the attached velocities before solution.

The MLCP (1.6) corresponds to a single step integrator. However, higher order integrators, such as Newmark methods, usually give rise to MLCPs with an equivalent form. Most ArtiSynth integrators use some variation of (1.6) to determine the system velocity at each time step.

To set up (1.6), 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 \hat{\bf f} and the stiffness/damping augmentation used to produce \hat{\bf M}; constrainers (Constrainer) supply {\bf G}, {\bf N}, {\bf g} and {\bf n}, and attachments (DynamicAttachment) provide the information needed to factor out attached velocities.

1.3 Basic packages

The core code of the ArtiSynth project is divided into three main packages, each with a number of sub-packages.

1.3.1 maspack

The packages under maspack contain general computational utilities that are independent of ArtiSynth and could be used in a variety of other contexts. The main packages are:

maspack.util               // general utilities
maspack.matrix             // matrix and linear algebra
maspack.graph              // graph algorithms
maspack.fileutil           // remote file access
maspack.properties         // property implementation
maspack.spatialmotion      // 3D spatial motion and dynamics
maspack.solvers            // LCP solvers and linear solver interfaces
maspack.render             // viewer and rendering classes
maspack.geometry           // 3D geometry and meshes
maspack.collision          // collision detection
maspack.widgets            // Java swing widgets for maspack data types
maspack.apps               // stand-alone programs based only on maspack

1.3.2 artisynth.core

The packages under artisynth.core contain the core code for ArtiSynth model components and its GUI infrastructure.

artisynth.core.util        // general ArtiSynth utilities
artisynth.core.modelbase   // base classes for model components
artisynth.core.materials   // materials for springs and finite elements
artisynth.core.mechmodels  // basic mechanical models
artisynth.core.femmodels   // finite element models
artisynth.core.probes      // input and output probes
artisynth.core.workspace   // RootModel and associated components
artisynth.core.driver      // start ArtiSynth and drive the simulation
artisynth.core.gui         // graphical interface
artisynth.core.inverse     // inverse tracking controller
artisynth.core.moviemaker  // used for making movies
artisynth.core.renderables // components that are strictly visual
artisynth.core.opensim     // OpenSim model parser (under development)
artisynth.core.mfreemodels // mesh free models (experimental, not supported)

1.3.3 artisynth.demos

These packages contain demonstration models that illustrate ArtiSynth’s modeling capabilities:

artisynth.demos.mech       // mechanical model demos
artisynth.demos.fem        // demos involving finite elements
artisynth.demos.inverse    // demos involving inverse control
artisynth.demos.tutorial   // demos in this manual

1.4 Properties

ArtiSynth components expose properties, which provide a uniform interface for accessing their internal parameters and state. Properties vary from component to component; those for RigidBody include position, orientation, mass, and density, while those for AxialSpring include restLength and material. Properties are particularly useful for automatically creating control panels and probes, as described in Section 5. They are also used for automating component serialization.

Properties are described only briefly in this section; more detailed descriptions are available in the Maspack Reference Manual and the overview paper.

The set of properties defined for a component is fixed for that component’s class; while property values may vary between component instances, their definitions are class-specific. Properties are exported by a class through code contained in the class definition, as described in Section 5.2.

1.4.1 Querying and setting property values

Each property has a unique name that can be used to access its value interactively in the GUI. This can be done either by using a custom control panel (Section 5.1) or by selecting the component and choosing Edit properties ... from the right-click context menu).

Properties can also be accessed in code using their set/get accessor methods. Unless otherwise specified, the names for these are formed by simply prepending set or get to the property’s name. More specifically, a property with the name foo and a value type of Bar will usually have accessor signatures of

  Bar getFoo()
  void setFoo (Bar value)

1.4.2 Property handles and paths

A property’s name can also be used to obtain a property handle through which its value may be queried or set generically. Property handles are implemented by the class Property and are returned by the component’s getProperty() method. getProperty() takes a property’s name and returns the corresponding handle. For example, components of type Muscle have a property excitation, for which a handle may be obtained using a code fragment such as

  Muscle muscle;
  ...
  Property prop = muscle.getProperty ("excitation");

Property handles can also be obtained for subcomponents, using a property path that consists of a path to the subcomponent followed by a colon ‘:’ and the property name. For example, to obtain the excitation property for a subcomponent located by axialSprings/lad relative to a MechModel, one could use a call of the form

  MechModel mech;
  ...
  Property prop = mech.getProperty ("axialSprings/lad:excitation");

1.4.3 Composite and inheritable properties

Figure 1.1: Inheritance of a property named stiffness among a component hierarchy. Explicit settings are in bold; inherited settings are in gray italic.

Composite properties are possible, in which a property value is a composite object that in turn has subproperties. A good example of this is the RenderProps class, which is associated with the property renderProps for renderable objects and which itself can have a number of subproperties such as visible, faceStyle, faceColor, lineStyle, lineColor, etc.

Properties can be declared to be inheritable, so that their values can be inherited from the same properties hosted by ancestor components further up the component hierarchy. Inheritable properties require a more elaborate declaration and are associated with a mode which may be either Explicit or Inherited. If a property’s mode is inherited, then its value is obtained from the closest ancestor exposing the same property whose mode is explicit. In Figure (1.1), the property stiffness is explicitly set in components A, C, and E, and inherited in B and D (which inherit from A) and F (which inherits from C).

1.5 Creating an application model

ArtiSynth applications are created by writing and compiling an application model that is a subclass of RootModel. This application-specific root model is then loaded and run by the ArtiSynth program.

The code for the application model should:

  • Declare a no-args constructor

  • Override the RootModel build() method to construct the application.

ArtiSynth can load a model either using the build method or by reading it from a file:

Build method

ArtiSynth creates an instance of the model using the no-args constructor, assigns it a name (which is either user-specified or the simple name of the class), and then calls the build() method to perform the actual construction.

Reading from a file

ArtiSynth creates an instance of the model using the no-args constructor, and then the model is named and constructed by reading the file.

The no-args constructor should perform whatever initialization is required in both cases, while the build() method takes the place of the file specification. Unless a model is originally created using a file specification (which is very tedious), the first time creation of a model will almost always entail using the build() method.

The general template for application model code looks like this:

package artisynth.models.experimental; // package where the model resides
import artisynth.core.workspace.RootModel;
... other imports ...
public class MyModel extends RootModel {
   // no-args constructor
   public MyModel() {
      ... basic initialization ...
   }
   // build method to do model construction
   public void build (String[] args) {
      ... code to build the model ....
   }
}

Here, the model itself is called MyModel, and is defined in the (hypothetical) package artisynth.models.experimental (placing models in the super package artisynth.models is common practice but not necessary).

Note: The build() method was only introduced in ArtiSynth 3.1. Prior to that, application models were constructed using a constructor taking a String argument supplying the name of the model. This method of model construction still works but is deprecated.

1.5.1 Implementing the build() method

As mentioned above, the build() method is responsible for actual model construction. Many applications are built using a single top-level MechModel. Build methods for these may look like the following:

   public void build (String[] args) {
      MechModel mech = new MechModel("mech");
      addModel (mech);
      ... create and add components to the mech model ...
      ... create and add any needed agents to the root model ...
   }

First, a MechModel is created (with the name "mech" in this example, although any name, or no name, may be given) and added to the list of models in the root model using the addModel() method. Subsequent code then creates and adds the components required by the MechModel, as described in Sections 3, 4 and 6. The build() method also creates and adds to the root model any agents required by the application (controllers, probes, etc.), as described in Section 5.

When constructing a model, there is no fixed order in which components need to be added. For instance, in the above example, addModel(mech) could be called near the end of the build() method rather than at the beginning. The only restriction is that when a component is added to the hierarchy, all other components that it refers to should already have been added to the hierarchy. For instance, an axial spring (Section 3.1) refers to two points. When it is added to the hierarchy, those two points should already be present in the hierarchy.

The build() method supplies a String array as an argument, which can be used to transmit application arguments in a manner analogous to the args argument passed to static main() methods. Build arguments can be specified when a model is loaded directly from a class using Models > Load from class ..., or when the startup model is set to automatically load a model when ArtiSynth is first started (Settings > Startup model). Details are given in the “Loading, Simulating and Saving Models” section of the User Interface Guide.

Build arguments can also be listed directly on the ArtiSynth command line when specifying a model to load using the -model <classname> option. This is done by enclosing the desired arguments within square brackets [ ] immediately following the -model option. So, for example,

> artisynth -model projects.MyModel [ -size 50 ]

would cause the strings "-size" and "50" to be passed to the build() method of MyModel.

1.5.2 Making models visible to ArtiSynth

In order to load an application model into ArtiSynth, the classes associated with its implementation must be made visible to ArtiSynth. This usually involves adding the top-level class folder associated with the application code to the classpath used by ArtiSynth.

The demonstration models referred to in this guide belong to the package artisynth.demos.tutorial and are already visible to ArtiSynth.

In most current ArtiSynth projects, classes are stored in a folder tree separate from the source code, with the top-level class folder named classes, located one level below the project root folder. A typical top-level class folder might be stored in a location like this:

  /home/joeuser/artisynthProjects/classes

In the example shown in Section 1.5, the model was created in the package artisynth.models.experimental. Since Java classes are arranged in a folder structure that mirrors package names, with respect to the sample project folder shown above, the model class would be located in

  /home/joeuser/artisynthProjects/classes/artisynth/models/experimental

At present there are three ways to make top-level class folders known to ArtiSynth:

Add projects to your Eclipse launch configuration

If you are using the Eclipse IDE, then you can add the project in which are developing your model code to the launch configuration that you use to run ArtiSynth. Other IDEs will presumably provide similar functionality.

Add the folders to the external classpath

You can explicitly add the class folders to ArtiSynth’s external classpath. The easiest way to do this is to select “Settings > External classpath ...” from the Settings menu, which will open an external classpath editor which lists all the classpath entries in a large panel on the left. (When ArtiSynth is first installed, the external classpath has no entries, and so this panel will be blank.) Class folders can then by added via the “Add class folder” button, and the classpath is saved using the Save button.

Add the folders to your CLASSPATH environment variable

If you are running ArtiSynth from the command line, using the artisynth command (or artisynth.bat on Windows), then you can define a CLASSPATH environment variable in your environment and add the needed folders to this.

All of these methods are described in more detail in the “Installing External Models and Packages” section of the ArtiSynth Installation Guide (available for Linux, Windows, and MacOS).

1.5.3 Loading and running a model

If a model’s classes are visible to ArtiSynth, then it may be loaded into ArtiSynth in several ways:

Loading from the Model menu

If the root model is contained in a package located under artisynth.demos or artisynth.models, then it will appear in the default model menu (Models in the main menu bar) under the submenu All demos or All models.

Loading by class path

A model may also be loaded by choosing “Load from class ...” from the Models menu and specifying its package name and then choosing its root model class. It is also possible to use the -model <classname> command line argument to have a model loaded directly into ArtiSynth when it starts up.

Loading from a file

If a model has been saved to a .art file, it may be loaded from that file by choosing File > Load model ....

These methods are described in detail in the section “Loading and Simulating Models” of the ArtiSynth User Interface Guide.

The demonstration models referred to in this guide should already be present in the model menu and may be loaded from the submenu Models > All demos > tutorial.

Once a model is loaded, it can be simulated, or run. Simulation of the model can then be started, paused, single-stepped, or reset using the play controls (Figure 1.2) located at the upper right of the ArtiSynth window frame. Starting and stopping a simulation is done by clicking play/pause, while reset resets the simulation to time 0. The single-step button advances the simulation by one time step. The stop-all button will also stop the simulation, along with any Jython commands or scripts that are running.

Figure 1.2: The ArtiSynth play controls. From left to right: step size control, current simulation time, and the reset, skip-back, play/pause, single-step, skip-forward and stop-all buttons.

Comprehensive information on exploring and interacting with models is given in the ArtiSynth User Interface Guide.

Chapter 2 Supporting classes

ArtiSynth uses a large number of supporting classes, mostly defined in the super package maspack, for handling mathematical and geometric quantities. Those that are referred to in this manual are summarized in this section.

2.1 Vectors and matrices

Among the most basic classes are those used to implement vectors and matrices, defined in maspack.matrix. All vector classes implement the interface Vector and all matrix classes implement Matrix, which provide a number of standard methods for setting and accessing values and reading and writing from I/O streams.

General sized vectors and matrices are implemented by VectorNd and MatrixNd. These provide all the usual methods for linear algebra operations such as addition, scaling, and multiplication:

  VectorNd v1 = new VectorNd (5);        // create a 5 element vector
  VectorNd v2 = new VectorNd (5);
  VectorNd vr = new VectorNd (5);
  MatrixNd M = new MatrixNd (5, 5);      // create a 5 x 5 matrix
  M.setIdentity();                       // M = I
  M.scale (4);                           // M = 4*M
  v1.set (new double[] {1, 2, 3, 4, 5}); // set values
  v2.set (new double[] {0, 1, 0, 2, 0});
  v1.add (v2);                           // v1 += v2
  M.mul (vr, v1);                        // vr = M*v1
  System.out.println ("result=" + vr.toString ("%8.3f"));

As illustrated in the above example, vectors and matrices both provide a toString() method that allows their elements to be formatted using a C-printf style format string. This is useful for providing concise and uniformly formatted output, particularly for diagnostics. The output from the above example is

  result=   4.000   12.000   12.000   24.000   20.000

Detailed specifications for the format string are provided in the documentation for NumberFormat.set(String). If either no format string, or the string "%g", is specified, toString() formats all numbers using the full-precision output provided by Double.toString(value).

For computational efficiency, a number of fixed-size vectors and matrices are also provided. The most commonly used are those defined for three dimensions, including Vector3d and Matrix3d:

  Vector3d v1 = new Vector3d (1, 2, 3);
  Vector3d v2 = new Vector3d (3, 4, 5);
  Vector3d vr = new Vector3d ();
  Matrix3d M = new Matrix3d();
  M.set (1, 2, 3,  4, 5, 6,  7, 8, 9);
  M.mul (vr, v1);        // vr = M * v1
  vr.scaledAdd (2, v2);  // vr += 2*v2;
  vr.normalize();        // normalize vr
  System.out.println ("result=" + vr.toString ("%8.3f"));

2.2 Rotations and transformations

maspack.matrix contains a number classes that implement rotation matrices, rigid transforms, and affine transforms.

Rotations (Section A.1) are commonly described using a RotationMatrix3d, which implements a rotation matrix and contains numerous methods for setting rotation values and transforming other quantities. Some of the more commonly used methods are:

   RotationMatrix3d();         // create and set to the identity
   RotationMatrix3d(u, angle); // create and set using an axis-angle
   setAxisAngle (u, ang);      // set using an axis-angle
   setRpy (roll, pitch, yaw);  // set using roll-pitch-yaw angles
   setEuler (phi, theta, psi); // set using Euler angles
   invert ();                  // invert this rotation
   mul (R)                     // post multiply this rotation by R
   mul (R1, R2);               // set this rotation to R1*R2
   mul (vr, v1);               // vr = R*v1, where R is this rotation

Rotations can also be described by AxisAngle, which characterizes a rotation as a single rotation about a specific axis.

Rigid transforms (Section A.2) are used by ArtiSynth to describe a rigid body’s pose, as well as its relative position and orientation with respect to other bodies and coordinate frames. They are implemented by RigidTransform3d, which exposes its rotational and translational components directly through the fields R (a RotationMatrix3d) and p (a Vector3d). Rotational and translational values can be set and accessed directly through these fields. In addition, RigidTransform3d provides numerous methods, some of the more commonly used of which include:

   RigidTransform3d();         // create and set to the identity
   RigidTransform3d(x, y, z);  // create and set translation to x, y, z
   // create and set translation to x, y, z and rotation to roll-pitch-yaw
   RigidTransform3d(x, y, z, roll, pitch, yaw);
   invert ();                  // invert this transform
   mul (T)                     // post multiply this transform by T
   mul (T1, T2);               // set this transform to T1*T2
   mulLeftInverse (T1, T2);    // set this transform to inv(T1)*T2

Affine transforms (Section A.3) are used by ArtiSynth to effect scaling and shearing transformations on components. They are implemented by AffineTransform3d.

Rigid transformations are actually a specialized form of affine transformation in which the basic transform matrix equals a rotation. RigidTransform3d and AffineTransform3d hence both derive from the same base class AffineTransform3dBase.

2.3 Points and Vectors

The rotations and transforms described above can be used to transform both vectors and points in space.

Vectors are most commonly implemented using Vector3d, while points can be implemented using the subclass Point3d. The only difference between Vector3d and Point3d is that the former ignores the translational component of rigid and affine transforms; i.e., as described in Sections A.2 and A.3, a vector v has an implied homogeneous representation of

{\bf v}^{*}\equiv\left(\begin{matrix}{\bf v}\\
0\end{matrix}\right), (2.1)

while the representation for a point p is

{\bf p}^{*}\equiv\left(\begin{matrix}{\bf p}\\
1\end{matrix}\right). (2.2)

Both classes provide a number of methods for applying rotational and affine transforms. Those used for rotations are

  void transform (R);             // this = R * this
  void transform (R, v1);         // this = R * v1
  void inverseTransform (R);      // this = inverse(R) * this
  void inverseTransform (R, v1);  // this = inverse(R) * v1

where R is a rotation matrix and v1 is a vector (or a point in the case of Point3d).

The methods for applying rigid or affine transforms include:

  void transform (X);             // transforms this by X
  void transform (X, v1);         // sets this to v1 transformed by X
  void inverseTransform (X);      // transforms this by the inverse of X
  void inverseTransform (X, v1);  // sets this to v1 transformed by inverse of X

where X is a rigid or affine transform. As described above, in the case of Vector3d, these methods ignore the translational part of the transform and apply only the matrix component (R for a RigidTransform3d and A for an AffineTransform3d). In particular, that means that for a RigidTransform3d given by X and a Vector3d given by v, the method calls

  v.transform (X.R)
  v.transform (X)

produce the same result.

2.4 Spatial vectors and inertias

The velocities, forces and inertias associated with 3D coordinate frames and rigid bodies are represented using the 6 DOF spatial quantities described in Sections A.5 and A.6. These are implemented by classes in the package maspack.spatialmotion.

Spatial velocities (or twists) are implemented by Twist, which exposes its translational and angular velocity components through the publicly accessible fields v and w, while spatial forces (or wrenches) are implemented by Wrench, which exposes its translational force and moment components through the publicly accessible fields f and m.

Both Twist and Wrench contain methods for algebraic operations such as addition and scaling. They also contain transform() methods for applying rotational and rigid transforms. The rotation methods simply transform each component by the supplied rotation matrix. The rigid transform methods, on the other hand, assume that the supplied argument represents a transform between two frames fixed within a rigid body, and transform the twist or wrench accordingly, using either (A.27) or (A.29).

The spatial inertia for a rigid body is implemented by SpatialInertia, which contains a number of methods for setting its value given various mass, center of mass, and inertia values, and querying the values of its components. It also contains methods for scaling and adding, transforming between coordinate systems, inversion, and multiplying by spatial vectors.

2.5 Meshes

ArtiSynth makes extensive use of 3D meshes, which are defined in maspack.geometry. They are used for a variety of purposes, including visualization, collision detection, and computing physical properties (such as inertia or stiffness variation within a finite element model).

A mesh is essentially a collection of vertices (i.e., points) that are topologically connected in some way. All meshes extend the abstract base class MeshBase, which supports the vertex definitions, while subclasses provide the topology.

Through MeshBase, all meshes provide methods for adding and accessing vertices. Some of these include:

  int numVertices();                 // return the number of vertices
  Vertex3d getVertex (int idx);      // return the idx-th vertex
  void addVertex (Vertex3d vtx);     // add vertex vtx to the mesh
  Vertex3d addVertex (Point3d p);    // create and return a vertex at position p
  void removeVertex (Vertex3d vtx);  // remove vertex vtx for the mesh
  ArrayList<Vertex3d> getVertices(); // return the list of vertices

Vertices are implemented by Vertex3d, which defines the position of the vertex (returned by the method getPosition()), and also contains support for topological connections. In addition, each vertex maintains an index, obtainable via getIndex(), that equals the index of its location within the mesh’s vertex list. This makes it easy to set up parallel array structures for augmenting mesh vertex properties.

Mesh subclasses currently include:

PolygonalMesh

Implements a 2D surface mesh containing faces implemented using half-edges.

PolylineMesh

Implements a mesh consisting of connected line-segments (polylines).

PointMesh

Implements a point cloud with no topological connectivity.

PolygonalMesh is used quite extensively and provides a number of methods for implementing faces, including:

  int numFaces();                 // return the number of faces
  Face getFace (int idx);         // return the idx-th face
  Face addFace (int[] vidxs);     // create and add a face using vertex indices
  void removeFace (Face f);       // remove the face f
  ArrayList<Face> getFaces();     // return the list of faces

The class Face implements a face as a counter-clockwise arrangement of vertices linked together by half-edges (class HalfEdge). Face also supplies a face’s (outward facing) normal via getNormal().

Some mesh uses within ArtiSynth, such as collision detection, require a triangular mesh; i.e., one where all faces have three vertices. The method isTriangular() can be used to check for this. Meshes that are not triangular can be made triangular using triangulate().

2.5.1 Mesh creation

Meshes are most commonly created using either one of the factory methods supplied by MeshFactory, or by reading a definition from a file (Section 2.5.5). However, it is possible to create a mesh by direct construction. For example, the following code fragment creates a simple closed tetrahedral surface:

   // a simple four-faced tetrahedral mesh
   PolygonalMesh mesh = new PolygonalMesh();
   mesh.addVertex (0, 0, 0);
   mesh.addVertex (1, 0, 0);
   mesh.addVertex (0, 1, 0);
   mesh.addVertex (0, 0, 1);
   mesh.addFace (new int[] { 0, 2, 1 });
   mesh.addFace (new int[] { 0, 3, 2 });
   mesh.addFace (new int[] { 0, 1, 3 });
   mesh.addFace (new int[] { 1, 2, 3 });

Some of the more commonly used factory methods for creating polyhedral meshes include:

  MeshFactory.createSphere (radius, nslices, nlevels);
  MeshFactory.createIcosahedralSphere (radius, divisons);
  MeshFactory.createBox (widthx, widthy, widthz);
  MeshFactory.createCylinder (radius, height, nslices);
  MeshFactory.createPrism (double[] xycoords, height);
  MeshFactory.createTorus (rmajor, rminor, nmajor, nminor);

Each factory method creates a mesh in some standard coordinate frame. After creation, the mesh can be transformed using the transform(X) method, where X is either a rigid transform ( RigidTransform3d) or a more general affine transform (AffineTransform3d). For example, to create a rotated box centered on (5,6,7), one could do:

  // create a box centered at the origin with widths 10, 20, 30:
  PolygonalMesh box = MeshFactory.createBox (10, 20, 20);
  // move the origin to 5, 6, 7 and rotate using roll-pitch-yaw
  // angles 0, 0, 45 degrees:
  box.transform (
     new RigidTransform3d (5, 6, 7,  0, 0, Math.toRadians(45)));

One can also scale a mesh using scale(s), where s is a single scale factor, or scale(sx,sy,sz), where sx, sy, and sz are separate scale factors for the x, y and z axes. This provides a useful way to create an ellipsoid:

   // start with a unit sphere with 12 slices and 6 levels ...
  PolygonalMesh ellipsoid = MeshFactory.createSphere (1.0, 12, 6);
  // and then turn it into an ellipsoid by scaling about the axes:
  ellipsoid.scale (1.0, 2.0, 3.0);

MeshFactory can also be used to create new meshes by performing Boolean operations on existing ones:

  MeshFactory.getIntersection (mesh1, mesh2);
  MeshFactory.getUnion (mesh1, mesh2);
  MeshFactory.getSubtraction (mesh1, mesh2);

2.5.2 Setting normals, colors, and textures

Meshes provide support for adding normal, color, and texture information, with the exact interpretation of these quantities depending upon the particular mesh subclass. Most commonly this information is used simply for rendering, but in some cases normal information might also be used for physical simulation.

For polygonal meshes, the normal information described here is used only for smooth shading. When flat shading is requested, normals are determined directly from the faces themselves.

Normal information can be set and queried using the following methods:

  setNormals (
     List<Vector3d> nrmls, int[] indices);  // set all normals and indices
  ArrayList<Vector3d> getNormals();         // get all normals
  int[] getNormalIndices();                 // get all normal indices
  int numNormals();                         // return the number of normals
  Vector3d getNormal (int idx);             // get the normal at index idx
  setNormal (int idx, Vector3d nrml);       // set the normal at index idx
  clearNormals();                           // clear all normals and indices

The method setNormals() takes two arguments: a set of normal vectors (nrmls), along with a set of index values (indices) that map these normals onto the vertices of each of the mesh’s geometric features. Often, there will be one unique normal per vertex, in which case nrmls will have a size equal to the number of vertices, but this is not always the case, as described below. Features for the different mesh subclasses are: faces for PolygonalMesh, polylines for PolylineMesh, and vertices for PointMesh. If indices is specified as null, then normals is assumed to have a size equal to the number of vertices, and an appropriate index set is created automatically using createVertexIndices() (described below). Otherwise, indices should have a size of equal to the number of features times the number of vertices per feature. For example, consider a PolygonalMesh consisting of two triangles formed from vertex indices (0, 1, 2) and (2, 1, 3), respectively. If normals are specified and there is one unique normal per vertex, then the normal indices are likely to be

   [ 0 1 2  2 1 3 ]

As mentioned above, sometimes there may be more than one normal per vertex. This happens in cases when the same vertex uses different normals for different faces. In such situations, the size of the nrmls argument will exceed the number of vertices.

The method setNormals() makes internal copies of the specified normal and index information, and this information can be later read back using getNormals() and getNormalIndices(). The number of normals can be queried using numNormals(), and individual normals can be queried or set using getNormal(idx) and setNormal(idx,nrml). All normals and indices can be explicitly cleared using clearNormals().

Color and texture information can be set using analogous methods. For colors, we have

  setColors (
     List<float[]> colors, int[] indices);  // set all colors and indices
  ArrayList<float[]> getColors();           // get all colors
  int[] getColorIndices();                  // get all color indices
  int numColors();                          // return the number of colors
  float[] getColor (int idx);               // get the color at index idx
  setColor (int idx, float[] color);        // set the color at index idx
  setColor (int idx, Color color);          // set the color at index idx
  setColor (
     int idx, float r, float g, float b, float a); // set the color at index idx
  clearColors();                            // clear all colors and indices

When specified as float[], colors are given as RGB or RGBA values, in the range [0,1], with array lengths of 3 and 4, respectively. The colors returned by getColors() are always RGBA values.

With colors, there may often be fewer colors than the number of vertices. For instance, we may have only two colors, indexed by 0 and 1, and want to use these to alternately color the mesh faces. Using the two-triangle example above, the color indices might then look like this:

   [ 0 0 0 1 1 1 ]

Finally, for texture coordinates, we have

  setTextureCoords (
     List<Vector3d> coords, int[] indices); // set all texture coords and indices
  ArrayList<Vector3d> getTextureCoords();   // get all texture coords
  int[] getTextureIndices();                // get all texture indices
  int numTextureCoords();                   // return the number of texture coords
  Vector3d getTextureCoords (int idx);      // get texture coords at index idx
  setTextureCoords (int idx, Vector3d coords);// set texture coords at index idx
  clearTextureCoords();                     // clear all texture coords and indices

When specifying indices using setNormals, setColors, or setTextureCoords, it is common to use the same index set as that which associates vertices with features. For convenience, this index set can be created automatically using

   int[] createVertexIndices();

Alternatively, we may sometimes want to create a index set that assigns the same attribute to each feature vertex. If there is one attribute per feature, the resulting index set is called a feature index set, and can be created using

   int[] createFeatureIndices();

If we have a mesh with three triangles and one color per triangle, the resulting feature index set would be

   [ 0 0 0 1 1 1 2 2 2 ]

Note: when a mesh is modified by the addition of new features (such as faces for PolygonalMesh), all normal, color and texture information is cleared by default (with normal information being automatically recomputed on demand if automatic normal creation is enabled; see Section 2.5.3). When a mesh is modified by the removal of features, the index sets for normals, colors and textures are adjusted to account for the removal.
For colors, it is possible to request that a mesh explicitly maintain colors for either its vertices or features (Section 2.5.4). When this is done, colors will persist when vertices or features are added or removed, with default colors being automatically created as necessary.

Once normals, colors, or textures have been set, one may want to know which of these attributes are associated with the vertices of a specific feature. To know this, it is necessary to find that feature’s offset into the attribute’s index set. This offset information can be found using the array returned by

  int[] getFeatureIndexOffsets()

For example, the three normals associated with a triangle at index ti can be obtained using

   int[] indexOffs = mesh.getFeatureIndexOffsets();
   ArrayList<Vector3d> nrmls = mesh.getNormals();
   // get the three normals associated with the triangle at index ti:
   Vector3d n0 = nrmls.get (indexOffs[ti]);
   Vector3d n1 = nrmls.get (indexOffs[ti]+1);
   Vector3d n2 = nrmls.get (indexOffs[ti]+2);

Alternatively, one may use the convenience methods

   Vector3d getFeatureNormal (int fidx, int k);
   float[] getFeatureColor (int fidx, int k);
   Vector3d getFeatureTextureCoords (int fidx, int k);

which return the attribute values for the k-th vertex of the feature indexed by fidx.

In general, the various get methods return references to internal storage information and so should not be modified. However, specific values within the lists returned by getNormals(), getColors(), or getTextureCoords() may be modified by the application. This may be necessary when attribute information changes as the simulation proceeds. Alternatively, one may use methods such as setNormal(idx,nrml) setColor(idx,color), or setTextureCoords(idx,coords).

Also, in some situations, particularly with colors and textures, it may be desirable to not have color or texture information defined for certain features. In such cases, the corresponding index information can be specified as -1, and the getNormal(), getColor() and getTexture() methods will return null for the features in question.

2.5.3 Automatic creation of normals and hard edges

For some mesh subclasses, if normals are not explicitly set, they are computed automatically whenever getNormals() or getNormalIndices() is called. Whether or not this is true for a particular mesh can be queried by the method

   boolean hasAutoNormalCreation();

Setting normals explicitly, using a call to setNormals(nrmls,indices), will overwrite any existing normal information, automatically computed or otherwise. The method

   boolean hasExplicitNormals();

will return true if normals have been explicitly set, and false if they have been automatically computed or if there is currently no normal information. To explicitly remove normals from a mesh which has automatic normal generation, one may call setNormals() with the nrmls argument set to null.

More detailed control over how normals are automatically created may be available for specific mesh subclasses. For example, PolygonalMesh allows normals to be created with multiple normals per vertex, for vertices that are associated with either open or hard edges. This ability can be controlled using the methods

   boolean getMultipleAutoNormals();
   setMultipleAutoNormals (boolean enable);

Having multiple normals means that even with smooth shading, open or hard edges will still appear sharp. To make an edge hard within a PolygonalMesh, one may use the methods

   boolean setHardEdge (Vertex3d v0, Vertex3d v1);
   boolean setHardEdge (int vidx0, int vidx1);
   boolean hasHardEdge (Vertex3d v0, Vertex3d v1);
   boolean hasHardEdge (int vidx0, int vidx1);
   int numHardEdges();
   int clearHardEdges();

which control the hardness of edges between individual vertices, specified either directly or using their indices.

2.5.4 Vertex and feature coloring

The method setColors() makes it possible to assign any desired coloring scheme to a mesh. However, it does require that the user explicitly reset the color information whenever new features are added.

For convenience, an application can also request that a mesh explicitly maintain colors for either its vertices or features. These colors will then be maintained when vertices or features are added or removed, with default colors being automatically created as necessary.

Vertex-based coloring can be requested with the method

   setVertexColoringEnabled();

This will create a separate (default) color for each of the mesh’s vertices, and set the color indices to be equal to the vertex indices, which is equivalent to the call

   setColors (colors, createVertexIndices());

where colors contains a default color for each vertex. However, once vertex coloring is enabled, the color and index sets will be updated whenever vertices or features are added or removed. Meanwhile, applications can query or set the colors for any vertex using getColor(idx), or any of the various setColor methods. Whether or not vertex coloring is enabled can be queried using

   getVertexColoringEnabled();

Once vertex coloring is established, the application will typically want to set the colors for all vertices, perhaps using a code fragment like this:

   mesh.setVertexColoringEnabled();
   for (int i=0; i<mesh.numVertices(); i++) {
      ... compute color for the vertex ...
      mesh.setColor (i, color);
   }

Similarly, feature-based coloring can be requested using the method

   setFeatureColoringEnabled();

This will create a separate (default) color for each of the mesh’s features (faces for PolygonalMesh, polylines for PolylineMesh, etc.), and set the color indices to equal the feature index set, which is equivalent to the call

   setColors (colors, createFeatureIndices());

where colors contains a default color for each feature. Applications can query or set the colors for any vertex using getColor(idx), or any of the various setColor methods. Whether or not feature coloring is enabled can be queried using

   getFeatureColoringEnabled();

2.5.5 Reading and writing mesh files

PolygonalMesh, PolylineMesh, and PointMesh all provide constructors that allow them to be created from a definition file, with the file format being inferred from the file name suffix:

   PolygonalMesh (String fileName) throws IOException
   PolygonalMesh (File file) throws IOException
   PolylineMesh (String fileName) throws IOException
   PolylineMesh (File file) throws IOException
   PointMesh (String fileName) throws IOException
   PointMesh (File file) throws IOException
Suffix Format PolygonalMesh PolylineMesh PointMesh
.obj Alias Wavefront X X X
.ply Polygon file format X X
.stl STereoLithography X
.gts GNU triangulated surface X
.off Object file format X
.vtk VTK ascii format X
.vtp VTK XML format X X
Table 2.1: Mesh file formats which are supported for different mesh types

The currently supported file formats, and their applicability to the different mesh types, are given in Table 2.1. For example, a PolygonalMesh can be read from either an Alias Wavefront .obj file or an .stl file, as show in the following example:

   PolygonalMesh mesh0 = null;
   PolygonalMesh mesh1 = null;
   try {
      mesh0 = new PolygonalMesh ("meshes/torus.obj");
   }
   catch (IOException e) {
      System.err.println ("Can’t read mesh:");
      e.printStackTrace();
   }
   try {
      mesh1 = new PolygonalMesh ("meshes/cylinder.stl");
   }
   catch (IOException e) {
      System.err.println ("Can’t read mesh:");
      e.printStackTrace();
   }

The file-based mesh constructors may throw an I/O exception if an I/O error occurs or if the indicated format does not support the mesh type. This exception must either be caught, as in the example above, or thrown out of the calling routine.

In addition to file-based constructors, all mesh types implement read and write methods that allow a mesh to be read from or written to a file, with the file format again inferred from the file name suffix:

   read (File file) throws IOException
   write (File file) throws IOException
   read (File file, boolean zeroIndexed) throws IOException
   write (File file, String fmtStr, boolean zeroIndexed) throws IOException

For the latter methods, the argument zeroIndexed specifies zero-based vertex indexing in the case of Alias Wavefront .obj files, while fmtStr is a C-style format string specifying the precision and style with which the vertex coordinates should be written. (In the former methods, zero-based indexing is false and vertices are written using full precision.)

As an example, the following code fragment writes a mesh as an .stl file:

   PolygonalMesh mesh;
   ... initialize ...
   try {
      mesh.write (new File ("data/mymesh.obj"));
   }
   catch (IOException e) {
      System.err.println ("Can’t write mesh:");
      e.printStackTrace();
   }

Sometimes, more explicit control is needed when reading or writing a mesh from/to a given file format. The constructors and read/write methods described above make use of a specific set of reader and writer classes located in the package maspack.geometry.io. These can be used directly to provide more explicit read/write control. The readers and writers (if implemented) associated with the different formats are given in Table 2.2.

Suffix Format Reader class Writer class
.obj Alias Wavefront WavefrontReader WavefrontWriter
.ply Polygon file format PlyReader PlyWriter
.stl STereoLithography StlReader StlWriter
.gts GNU triangulated surface GtsReader GtsWriter
.off Object file format OffReader OffWriter
.vtk VTK ascii format VtkAsciiReader
.vtp VTK XML format VtkXmlReader
Table 2.2: Reader and writer classes associated with the different mesh file formats

The general usage pattern for these classes is to construct the desired reader or writer with a path to the desired file, and then call readMesh() or writeMesh() as appropriate:

   // read a mesh from a .obj file:
   WavefrontReader reader = new WavefrontReader ("meshes/torus.obj");
   PolygonalMesh mesh = null;
   try {
      mesh = reader.readMesh();
   }
   catch (IOException e) {
      System.err.println ("Can’t read mesh:");
      e.printStackTrace();
   }

Both readMesh() and writeMesh() may throw I/O exceptions, which must be either caught, as in the example above, or thrown out of the calling routine.

For convenience, one can also use the classes GenericMeshReader or GenericMeshWriter, which internally create an appropriate reader or writer based on the file extension. This enables the writing of code that does not depend on the file format:

   String fileName;
   ...
   PolygonalMesh mesh = null;
   try {
      mesh = (PolygonalMesh)GenericMeshReader.readMesh(fileName);
   }
   catch (IOException e) {
      System.err.println ("Can’t read mesh:");
      e.printStackTrace();
   }

Here, fileName can refer to a mesh of any format supported by GenericMeshReader. Note that the mesh returned by readMesh() is explicitly cast to PolygonalMesh. This is because readMesh() returns the superclass MeshBase, since the default mesh created for some file formats may be different from PolygonalMesh.

2.5.6 Reading and writing normal and texture information

When writing a mesh out to a file, normal and texture information are also written if they have been explicitly set and the file format supports it. In addition, by default, automatically generated normal information will also be written if it relies on information (such as hard edges) that can’t be reconstructed from the stored file information.

Whether or not normal information will be written is returned by the method

   boolean getWriteNormals();

This will always return true if any of the conditions described above have been met. So for example, if a PolygonalMesh contains hard edges, and multiple automatic normals are enabled (i.e., getMultipleAutoNormals() returns true), then getWriteNormals() will return true.

Default normal writing behavior can be overridden within the MeshWriter classes using the following methods:

   int getWriteNormals()
   setWriteNormals (enable)

where enable should be one of the following values:

0

normals will never be written;

1

normals will always be written;

-1

normals will written according to the default behavior described above.

When reading a PolygonalMesh from a file, if the file contains normal information with multiple normals per vertex that suggests the existence of hard edges, then the corresponding edges are set to be hard within the mesh.

2.5.7 Constructive solid geometry

ArtiSynth contains primitives for performing constructive solid geometry (CSG) operations on volumes bounded by triangular meshes. The class that performs these operations is maspack.collision.SurfaceMeshIntersector, and it works by robustly determining the intersection contour(s) between a pair of meshes, and then using these to compute the triangles that need to be added or removed to produce the necessary CSG surface.

Figure 2.1: Dumbbell shaped mesh produced from the CSG union of two balls and a cylinder.

The CSG operations include union, intersection, and difference, and are implemented by the following methods of SurfaceMeshIntersector:

  findUnion (mesh0, mesh1);         // volume0 U volume1
  findIntersection (mesh0, mesh1);  // volume0 ^ volume1
  findDifference01 (mesh0, mesh1);  // volume0 - volume1
  findDifference10 (mesh0, mesh1);  // volume1 - volume

Each takes two PolyhedralMesh objects, mesh0 and mesh1, and creates and returns another PolyhedralMesh which represents the boundary surface of the requested operation. If the result of the operation is null, the returned mesh will be empty.

The example below uses findUnion to create a dumbbell shaped mesh from two balls and a cylinder:

  // first create two ball meshes and a bar mesh
  double radius = 1.0;
  int division = 1; // number of divisons for icosahedral sphere
  PolygonalMesh ball0 = MeshFactory.createIcosahedralSphere (radius, division);
  ball0.transform (new RigidTransform3d (0, -2*radius, 0));
  PolygonalMesh ball1 = MeshFactory.createIcosahedralSphere (radius, division);
  ball1.transform (new RigidTransform3d (0, 2*radius, 0));
  PolygonalMesh bar = MeshFactory.createCylinder (
     radius/2, radius*4, /*ns=*/32, /*nr=*/1, /*nh*/10);
  bar.transform (new RigidTransform3d (0, 0, 0, 0, 0, Math.PI/2));
  // use a SurfaceMeshIntersector to create a CSG union of these meshes
  SurfaceMeshIntersector smi = new SurfaceMeshIntersector();
  PolygonalMesh balls = smi.findUnion (ball0, ball1);
  PolygonalMesh mesh = smi.findUnion (balls, bar);

The balls and cylinder are created using the MeshFactory methods createIcosahedralSphere() and createCylinder(), where the latter takes arguments ns, nr, and nh giving the number of slices along the circumference, end-cap radius, and length. The final resulting mesh is shown in Figure 2.1.

2.6 Reading source relative files

ArtiSynth applications frequently need to read in various kinds of data files, including mesh files (as discussed in Section 2.5.5), FEM mesh geometry (Section 6.2.2), probe data (Section 5.4.4), and custom application data.

Often these data files do not reside in an absolute location but instead in a location relative to the application’s class or source files. For example, it is common for applications to store geometric data in a subdirectory "geometry" located beneath the source directory. In order to access such files in a robust way, and ensure that the code does not break when the source tree is moved, it is useful to determine the application’s source (or class) directory at run time. ArtiSynth supplies several ways to conveniently handle this situation. First, the RootModel itself supplies the following methods:

  // find path to the root model’s source directory
  String findSourceDir ();
  // get path to a file specified relative to the root model’s source directory
  String getSourceRelativePath (String relpath);

The first method returns the path to the source directory of the root model, while the second returns the path to a file specified relative to the root model source directory. If the root model source directory cannot be found (see discussion at the end of this section) both methods return null. As a specific usage example, assume that we have an application model whose build() method needs to load in a mesh torus.obj from a subdirectory meshes located beneath the source directory. This could be done as follows:

  String pathToMesh = getSourceRelativePath ("meshes/torus.obj");
  // read the mesh from a .obj file :
  WavefrontReader reader = new WavefrontReader (pathToMesh);
  PolygonalMesh mesh = null;
  try {
     mesh = reader.readMesh () ;
  }
  catch (IOException e) {
     System.err.println ("Can’t read mesh:");
     e.printStackTrace () ;
  }

A more general path finding utility is provided by maspack.util.PathFinder, which provides several static methods for locating source and class directories:

  // find path to the source directory associated with classObj
  String findSourceDir (Object classObj);
  // get path to a file specified relative to classObj source directory
  String getSourceRelativePath (Object classObj, String relpath);
  // find path to the class directory associated with classObj
  String findClassDir (Object classObj);
  // get path to a file specified relative to classObj class directory
  String getClassRelativePath (Object classObj, String relpath);

The “find” methods return a string path to the indicated class or source directory, while the “relative path” methods locate the class or source directory and append the additional path relpath. For all of these, the class is determined from classObj, either directly (if it is an instance of Class), by name (if it is a String), or otherwise by calling classObj.getClass(). When identifying a package by name, the name should be either a fully qualified class name, or a simple name that can be located with respect to the packages obtained via Package.getPackages(). For example, if we have a class whose fully qualified name is artisynth.models.test.Foo, then the following calls should all return the same result:

   Foo foo = new Foo();
   PathFinder.findSourceDir (foo);
   PathFinder.findSourceDir (Foo.class);
   PathFinder.findSourceDir ("artisynth.models.test.Foo");
   PathFinder.findSourceDir ("Foo");

If the source directory for Foo happens to be /home/projects/src/artisynth/models/test, then

   PathFinder.getSourceRelativePath (foo, "geometry/mesh.obj");

will return /home/projects/src/artisynth/models/test/geometry/mesh.obj.

When calling PathFinder methods from within the relevant class, one can specify this as the classObj argument.

With respect to the above example locating the file "meshes/torus.obj", the call to the root model method getSourceRelativePath() could be replaced with

  String pathToMesh = PathFinder.getSourceRelativePath (
     this, "meshes/torus.obj");

Since this is assumed to be called from the root model’s build method, the “class” can be indicated by simply passing this to getSourceRelativePath().

As an alternative to placing data files in the source directory, one could place them in the class directory, and then use findClassDir() and getClassRelativePath(). If the data files were originally defined in the source directory, it will be necessary to copy them to the class directory. Some Java IDEs will perform this automatically.

The PathFinder methods work by climbing the class’s resource hierarchy. Source directories are assumed to be located relative to the parent of the root class directory, via one of the paths specified by getSourceRootPaths(). By default, this list includes "src", "source", and "bin". Additional paths can be added using addSourceRootPath(path), or the entire list can be set using setSourceRootPaths(paths).

At preset, source directories will not be found if the reference class is contained in a jar file.

2.7 Reading and caching remote files

ArtiSynth applications often require the use of large data files to specify items such as FEM mesh geometry, surface mesh geometry, or medical imaging data. The size of these files may make it inconvenient to store them in any version control system that is used to store the application source code. As an alternative, ArtiSynth provides a file manager utility that allows such files to be stored on a separate server, and then downloaded on-demand and cached locally. To use this, one starts by creating an instance of a FileManager, using the constructor

  FileManager (String downloadPath, String remoteSourceName)

where downloadPath is a path to the local directory where the downloaded file should be placed, and remoteSourceName is a URI indicating the remote server location of the files. After the file manager has been created, it can be used to fetch remote files and cache them locally, using various get methods:

  File get (String destName);
  File get (String destName, String sourceName);

Both of these look for the file destName specified relative to the local directory, and return a File handle for it if it is present. Otherwise, they attempt to download the file from the remote source location, place it in the local directory, and return a File handle for it. The location of the remote file is given relative to the remote source URI by destName for the first method and sourceName for the second.

A simple example of using a file manager within a RootModel build() method is given by the following fragment:

  // create the file manager ...
  FileManager fm = new FileManager (
    getSourceRelativePath("geometry"),
    "http://myserver.org/artisynth/data/geometry");
  // ... and use it to get a bone mesh file
  File meshFile = fm.get ("tibia.obj");

Here, a file manager is created that uses a local directory "geometry", located relative to the RootModel source directory (see Section 2.6), and looks for missing files relative to the URI

  http://myserver.org/artisynth/data/geometry

The get() method is then used to obtain the file "tibia.obj" from the local directory. If it is not already present, it is downloaded from the remote location.

The FileManager contains other features and functionality, and one should consult its API documentation for more information.

Chapter 3 Mechanical Models I

This section details how to build basic multibody-type mechanical models consisting of particles, springs, rigid bodies, joints, and other constraints.

3.1 Springs and particles

The most basic type of mechanical model consists simply of particles connected together by axial springs. Particles are implemented by the class Particle, which is a dynamic component containing a three-dimensional position state, a corresponding velocity state, and a mass. It is an instance of the more general base class Point, which is used to also implement spatial points such as markers which do not have a mass.

3.1.1 Axial springs and materials

An axial spring is a simple spring that connects two points and is implemented by the class AxialSpring. This is a force effector component that exerts equal and opposite forces on the two points, along the line separating them, with a magnitude f that is a function f(l,\dot{l}) of the distance l between the points, and the distance derivative \dot{l}.

Each axial spring is associated with an axial material, implemented by a subclass of AxialMaterial, that specifies the function f(l,\dot{l}). The most basic type of axial material is a LinearAxialMaterial, which determines f according to the linear relationship

f(l,\dot{l})=k(l-l_{0})+d\dot{l} (3.1)

where l_{0} is the rest length and k and d are the stiffness and damping terms. Both k and d are properties of the material, while l_{0} is a property of the spring.

Axial springs are assigned a linear axial material by default. More complex, nonlinear axial materials may be defined in the package artisynth.core.materials. Setting or querying a spring’s material may be done with the methods setMaterial() and getMaterial().

3.1.2 Example: a simple particle-spring model

Figure 3.1: ParticleSpring model loaded into ArtiSynth.

An complete application model that implements a simple particle-spring model is given below.

1 package artisynth.demos.tutorial;
2
3 import java.awt.Color;
4 import maspack.matrix.*;
5 import maspack.render.*;
6 import artisynth.core.mechmodels.*;
7 import artisynth.core.materials.*;
8 import artisynth.core.workspace.RootModel;
9
10 /**
11  * Demo of two particles connected by a spring
12  */
13 public class ParticleSpring extends RootModel {
14
15    public void build (String[] args) {
16
17       // create MechModel and add to RootModel
18       MechModel mech = new MechModel ("mech");
19       addModel (mech);
20
21       // create the components
22       Particle p1 = new Particle ("p1", /*mass=*/2, /*x,y,z=*/0, 0, 0);
23       Particle p2 = new Particle ("p2", /*mass=*/2, /*x,y,z=*/1, 0, 0);
24       AxialSpring spring = new AxialSpring ("spr", /*restLength=*/0);
25       spring.setPoints (p1, p2);
26       spring.setMaterial (
27          new LinearAxialMaterial (/*stiffness=*/20, /*damping=*/10));
28
29       // add components to the mech model
30       mech.addParticle (p1);
31       mech.addParticle (p2);
32       mech.addAxialSpring (spring);
33
34       p1.setDynamic (false);                // first particle set to be fixed
35
36       // increase model bounding box for the viewer
37       mech.setBounds (/*min=*/-1, 0, -1, /*max=*/1, 0, 0);
38       // set render properties for the components
39       RenderProps.setSphericalPoints (p1, 0.06, Color.RED);
40       RenderProps.setSphericalPoints (p2, 0.06, Color.RED);
41       RenderProps.setCylindricalLines (spring, 0.02, Color.BLUE);
42    }
43 }

Line 1 of the source defines the package in which the model class will reside, in this case artisynth.demos.tutorial. Lines 3-8 import definitions for other classes that will be used.

The model application class is named ParticleSpring and declared to extend RootModel (line 13), and the build() method definition begins at line 15. (A no-args constructor is also needed, but because no other constructors are defined, the compiler creates one automatically.)

To begin, the build() method creates a MechModel named "mech", and then adds it to the models list of the root model using the addModel() method (lines 18-19). Next, two particles, p1 and p2, are created, with masses equal to 2 and initial positions at 0, 0, 0, and 1, 0, 0, respectively (lines 22-23). Then an axial spring is created, with end points set to p1 and p2, and assigned a linear material with a stiffness and damping of 20 and 10 (lines 24-27). Finally, after the particles and the spring are created, they are added to the particles and axialSprings lists of the MechModel using the methods addParticle() and addAxialSpring() (lines 30-32).

At this point in the code, both particles are defined to be dynamically controlled, so that running the simulation would cause both to fall under the MechModel’s default gravity acceleration of (0,0,-9.8). However, for this example, we want the first particle to remain fixed in place, so we set it to be non-dynamic (line 34), meaning that the physical simulation will not update its position in response to forces (Section 3.1.3).

The remaining calls control aspects of how the model is graphically rendered. setBounds() (line 37) increases the model’s “bounding box” so that by default it will occupy a larger part of the viewer frustum. The convenience method RenderProps.setSphericalPoints() is used to set points p1 and p2 to render as solid red spheres with a radius of 0.06, while RenderProps.setCylindricalLines() is used to set spring to render as a solid blue cylinder with a radius of 0.02. More details about setting render properties are given in Section 4.3.

To run this example in ArtiSynth, select All demos > tutorial > ParticleSpring from the Models menu. The model should load and initially appear as in Figure 3.1. Running the model (Section 1.5.3) will cause the second particle to fall and swing about under gravity.

3.1.3 Dynamic, parametric, and attached 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’s dynamic property to false, so that it does not respond to force inputs. As shown in the example above, this can be done using the method setDynamic():

  comp.setDynamic (false);

The method isDynamic() can be used to query the dynamic property.

Dynamic components can also be attached to other dynamic components (as mentioned in Section 1.2) so that their positions and velocities are controlled by the master components that they are attached to. To attach a dynamic component, one creates an AttachmentComponent specifying the attachment connection and adds it to the MechModel, as described in Section 3.7. The method isAttached() can be used to determine if a component is attached, and if it is, getAttachment() can be used to find the corresponding AttachmentComponent.

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. The component will move so as to follow the other master component(s) to which it is attached.

3.1.4 Custom axial materials

Application authors may create their own axial materials by subclassing AxialMaterial and overriding the functions

  double computeF (l, ldot, l0, excitation);
  double computeDFdl (l, ldot, l0, excitation);
  double computeDFdldot (l, ldot, l0, excitation);
  boolean isDFdldotZero ();

where excitation is an additional excitation signal a, which is used to implement active springs and which in particular is used to implement axial muscles (Section 4.5), for which a is usually in the range [0,1].

The first three methods should return the values of

f(l,\dot{l},a),\quad\frac{\partial f(l,\dot{l},a)}{\partial l},\quad\text{and}%
\quad\frac{\partial f(l,\dot{l},a)}{\partial\dot{l}}, (3.2)

respectively, while the last method should return true if \partial f(l,\dot{l},a)/\partial\dot{l}\equiv 0; i.e., if it is always equals to 0.

3.1.5 Damping parameters

Mechanical models usually contain damping forces in addition to spring-type restorative forces. Damping generates forces that reduce dynamic component velocities, and is usually the major source of energy dissipation in the model. Damping forces can be generated by the spring components themselves, as described above.

A general damping can be set for all particles by setting the MechModel’s pointDamping property. This causes a force

{\bf f}_{i}=-d_{p}{\bf v}_{i} (3.3)

to be applied to all particles, where d_{p} is the value of the pointDamping and {\bf v}_{i} is the particle’s velocity.

pointDamping can be set and queried using the MechModel methods

  setPointDamping (double d);
  double getPointDamping();

In general, whenever a component has a property propX, that property can be set and queried in code using methods of the form

  setPropX (T d);
  T getPropX();

where T is the type associated with the property.

pointDamping can also be set for particles individually. This property is inherited (Section 1.4.3), so that if not set explicitly, it inherits the nearest explicitly set value in an ancestor component.

3.2 Rigid bodies

Rigid bodies are implemented in ArtiSynth by the class RigidBody, which is a dynamic component containing a six-dimensional position and orientation state, a corresponding velocity state, an inertia, and an optional surface mesh.

A rigid body is associated with its own 3D spatial coordinate frame, and is a subclass of the more general Frame component. The combined position and orientation of this frame with respect to world coordinates defines the body’s pose, and the associated 6 degrees of freedom describe its “position” state.

3.2.1 Frame markers

Figure 3.2: A force {\bf f} applied to a frame marker attached to a rigid body. The marker is located at the point {\bf r} with respect to the body coordinate frame B.

ArtiSynth makes extensive use of markers, which are (massless) points attached to dynamic components in the model. Markers are used for graphical display, implementing attachments, and transmitting forces back onto the underlying dynamic components.

A frame marker is a marker that can be attached to a Frame, and most commonly to a RigidBody (Figure 3.2). They are frequently used to provide the anchor points for attaching springs and, more generally, applying forces to the body.

Frame markers are implemented by the class FrameMarker, which is a subclass of Point. The methods

  Point3d getLocation();
  void setLocation (Point3d r);

get and set the marker’s location {\bf r} with respect to the frame’s coordinate system. When a 3D force {\bf f} is applied to the marker, it generates a spatial force \hat{\bf f} (Section A.5) on the frame given by

\hat{\bf f}=\left(\begin{matrix}{\bf f}\\
{\bf r}\times{\bf f}\end{matrix}\right). (3.4)

Frame markers can be created using a variety of constructors, including

  FrameMarker ();
  FrameMarker (String name);
  FrameMarker (Frame frame, Point3d loc);

where FrameMarker() creates an empty marker, FrameMarker(name) creates an empty marker with a name, and FrameMarker(frame,loc) creates an unnamed marker attached to frame at the location loc with respect to the frame’s coordinates. Once created, a marker’s frame can be set and queried with

  void setFrame (Frame frame);
  Frame getFrame ();

A frame marker can be added to a MechModel with the MechModel methods

  void addFrameMarker (FrameMarker mkr);
  void addFrameMarker (FrameMarker mkr, Frame frame, Point3d loc);

where addFrameMarker(mkr,frame,loc) also sets the frame and the marker’s location with respect to it.

MechModel also supplies convenience methods to create a marker, attach it to a frame, and add it to the model:

  FrameMarker addFrameMarker (Frame frame, Point3d loc);
  FrameMarker addFrameMarkerWorld (Frame frame, Point3d locw);

Both methods return the created marker. The first, addFrameMarker(frame,loc), places it at the location loc with respect to the frame, while addFrameMarkerWorld(frame,pos) places it at pos with respect to world coordinates.

3.2.2 Example: a simple rigid body-spring model

Figure 3.3: RigidBodySpring model loaded into ArtiSynth.

A simple rigid body-spring model is defined in

  artisynth.demos.tutorial.RigidBodySpring

This differs from ParticleSpring only in the build() method, which is listed below:

1    public void build (String[] args) {
2
3       // create MechModel and add to RootModel
4       MechModel mech = new MechModel ("mech");
5       addModel (mech);
6
7       // create the components
8       Particle p1 = new Particle ("p1", /*mass=*/2, /*x,y,z=*/0, 0, 0);
9       // create box and set its pose (position/orientation):
10       RigidBody box =
11          RigidBody.createBox ("box", /*wx,wy,wz=*/0.5, 0.3, 0.3, /*density=*/20);
12       box.setPose (new RigidTransform3d (/*x,y,z=*/0.75, 0, 0));
13       // create marker point and connect it to the box:
14       FrameMarker mkr = new FrameMarker (/*x,y,z=*/-0.25, 0, 0);
15       mkr.setFrame (box);
16
17       AxialSpring spring = new AxialSpring ("spr", /*restLength=*/0);
18       spring.setPoints (p1, mkr);
19       spring.setMaterial (
20          new LinearAxialMaterial (/*stiffness=*/20, /*damping=*/10));
21
22       // add components to the mech model
23       mech.addParticle (p1);
24       mech.addRigidBody (box);
25       mech.addFrameMarker (mkr);
26       mech.addAxialSpring (spring);
27
28       p1.setDynamic (false);               // first particle set to be fixed
29
30       // increase model bounding box for the viewer
31       mech.setBounds (/*min=*/-1, 0, -1, /*max=*/1, 0, 0);
32       // set render properties for the components
33       RenderProps.setSphericalPoints (p1, 0.06, Color.RED);
34       RenderProps.setSphericalPoints (mkr, 0.06, Color.RED);
35       RenderProps.setCylindricalLines (mkr, 0.02, Color.BLUE);
36    }

The differences from ParticleSpring begin at line 9. Instead of creating a second particle, a rigid body is created using the factory method RigidBody.createBox(), which takes x, y, z widths and a (uniform) density and creates a box-shaped rigid body complete with surface mesh and appropriate mass and inertia. As the box is initially centered at the origin, moving it elsewhere requires setting the body’s pose, which is done using setPose(). The RigidTransform3d passed to setPose() is created using a three-argument constructor that generates a translation-only transform. Next, starting at line 14, a FrameMarker is created for a location (-0.25,0,0)^{T} relative to the rigid body, and attached to the body using its setFrame() method.

The remainder of build() is the same as for ParticleSpring, except that the spring is attached to the frame marker instead of a second particle.

To run this example in ArtiSynth, select All demos > tutorial > RigidBodySpring from the Models menu. The model should load and initially appear as in Figure 3.3. Running the model (Section 1.5.3) will cause the rigid body to fall and swing about under gravity.

3.2.3 Creating rigid bodies

As illustrated above, rigid bodies can be created using factory methods supplied by RigidBody. Some of these include:

  createBox (name, widthx, widthy, widthz, density);
  createCylinder (name, radius, height, density, nsides);
  createSphere (name, radius, density, nslices);
  createEllipsoid (name, radx, rady, radz, density, nslices);

The bodies do not need to be named; if no name is desired, then name and can be specified as null.

In addition, there are also factory methods for creating a rigid body directly from a mesh:

  createFromMesh (name, mesh, density, scale);
  createFromMesh (name, meshFileName, density, scale);

These take either a polygonal mesh (Section 2.5), or a file name from which a mesh is read, and use it as the body’s surface mesh and then compute the mass and inertia properties from the specified (uniform) density.

When a body is created directly from a surface mesh, its center of mass will typically not be coincident with the origin of its coordinate frame. Section 3.2.6 discusses the implications of this and how to correct it.

Alternatively, one can create a rigid body directly from a constructor, and then set the mesh and inertia properties explicitly:

  PolygonalMesh femurMesh;
  SpatialInertia inertia;
  ... initialize mesh and inertia appropriately ...
  RigidBody body = new RigidBody ("femur");
  body.setMesh (femurMesh);
  body.setInertia (inertia);

3.2.4 Pose and velocity

A body’s pose can be set and queried using the methods

  setPose (RigidTransform3d T);   // sets the pose to T
  getPose (RigidTransform3d T);   // gets the current pose in T
  RigidTransform3d getPose();     // returns the current pose (read-only)

These use a RigidTransform3d (Section 2.2) to describe the pose. Body poses are described in world coordinates and specify the transform from body to world coordinates. In particular, the pose for a body A specifies the rigid transform {\bf T}_{AW}.

Rigid bodies also expose the translational and rotational components of their pose via the properties position and orientation, which can be queried and set independently using the methods

  setPosition (Point3d p);       // sets the position to p
  getPosition (Point3d p);       // gets the current position in p
  Point3d getPosition();         // returns the current position (read-only)
  setOrientation (AxisAngle a);  // sets the orientation to a
  getOrientation (AxisAngle a);  // gets the current orientation in a
  AxisAngle getOrientation();    // returns the current orientation (read-only)

The velocity of a rigid body is described using a Twist (Section 2.4), which contains both the translational and rotational velocities. The following methods set and query the spatial velocity as described with respect to world coordinates:

  setVelocity (Twist v);         // sets the spatial velocity to v
  getVelocity (Twist v);         // gets the current spatial velocity in v
  Twist getVelocity();           // returns current spatial velocity (read-only)

During simulation, unless a rigid body has been set to be parametric (Section 3.1.3), its pose and velocity are updated in response to forces, so setting the pose or velocity generally makes sense only for setting initial conditions. On the other hand, if a rigid body is parametric, then it is possible to control its pose during the simulation, but in that case it is better to set its target pose and/or target velocity, as described in Section 5.3.1.

3.2.5 Inertia and the surface mesh

The “mass” of a rigid body is described by its spatial inertia, which is a 6\times 6 matrix relating its spatial velocity to its spatial momentum (Section A.6). Within ArtiSynth, spatial inertia is described by a SpatialInertia object, which specifies its mass, center of mass (with respect to body coordinates), and rotational inertia (with respect to the center of mass).

Most rigid bodies are also associated with a polygonal surface mesh, which can be set and queried using the methods

  setSurfaceMesh (PolygonalMesh mesh);
  setSurfaceMesh (PolygonalMesh mesh, String meshFileName);
  PolygonalMesh getSurfaceMesh();

The second method takes an optional fileName argument that can be set to the name of a file from which the mesh was read. Then if the model itself is saved to a file, the model file will specify the mesh using the file name instead of explicit vertex and face information, which can reduce the model file size considerably.

Rigid bodies can also have more than one mesh, as described in Section 3.2.9.

The inertia of a rigid body can be explicitly set using a variety of methods including

  setInertia (M)                    // set using SpatialInertia M
  setInertia (mass, Jxx, Jyy, Jzz); // mass and diagonal rotational inertia
  setInertia (mass, J);             // mass and full rotational inertia
  setInertia (mass, J, com);        // mass, rotational inertia, center-of-mass

and can be queried using

  getInertia (M);                   // get SpatialInertia in M
  getInertia ();                    // return read-only SpatialInertia

In practice, it is often more convenient to simply specify a mass or a density, and then use the geometry of the surface mesh (and possibly other meshes, Section 3.2.9) to compute the remaining inertial values. How a rigid body’s inertia is computed is determined by its inertiaMethod property, which can be one

EXPLICIT

Inertia is set explicitly.

MASS

Inertia is determined implicitly from the mesh geometry and the body’s mass.

DENSITY

Inertia is determined implicitly from the mesh geometry and the body’s density (which is multiplied by the mesh volume(s) to determine a mass).

When using DENSITY to determine the inertia, it is generally assumed that the contributing meshes are both polygonal and closed. Meshes which are either open or non-polygonal generally do not have a well-defined volume which can be multiplied by the density to determine the mass.

The inertiaMethod property can be set and queried using

  setInertiaMethod (InertiaMethod method);
  InertiaMethod getInertiaMethod();

and its default value is DENSITY. Explicitly setting the inertia using one of setInertia() methods described above will set inertiaMethod to EXPLICIT. The method

  setInertiaFromDensity (density);

will (re)compute the inertia using the mesh geometry and a density value and set inertiaMethod to DENSITY, and the method

  setInertiaFromMass (mass);

will (re)compute the inertia using the mesh geometry and a mass value and set inertiaMethod to MASS.

Finally, the (assumed uniform) density of the body can be queried using

   getDensity();

There are some subtleties involved in determining the inertia using either the DENSITY or MASS methods when the rigid body contains more than one mesh. Details are given in Section 3.2.9.

3.2.6 Coordinate frames and the center of mass

Figure 3.4: Left: rigid body whose coordinate frame B is not coincident with the center of mass (COM). Right: same body, with its coordinate frame translated to be coincident with the COM.

It is important to note that the origin of a body’s coordinate frame will not necessarily coincide with its center of mass (COM), and in fact the frame origin does not even have to lie inside the body’s surface (Figure 3.4). This typically occurs when a body’s inertia is computed directly from its surface mesh (or meshes), as described in Section 3.2.5.

Having the COM differ from the frame origin may lead to some undesired effects. For instance, since the body’s spatial velocity is defined with respect to the frame origin and not the COM, if the two are not coincident, then a purely angular body velocity will cause the COM to translate. The body’s spatial inertia also becomes more complicated, with non-zero 3 x 3 blocks in the lower left and upper right (Section A.6), which can have a small effect on computational accuracy. Finally, manipulating a body’s pose in the ArtiSynth UI (as described in the section “Model Manipulation” in the ArtiSynth User Interface Guide) can also be more cumbersome if the origin is located far from the COM.

There are several ways to ensure that the COM and frame origin are coincident. The most direct is to call the method centerPoseOnCenterOfMass() after the body has been created:

   String meshFilePath = "/project/geometry/bodyMesh.obj";
   double density = 1000;
   PolygonalMesh mesh = new PolygonalMesh (meshFilePath); // read in a mesh
   RigidBody bodyA = RigidBody.createFromMesh (
      "bodyA", mesh, density, /*scale=*/1); // create body from the mesh
   bodyA.centerPoseOnCenterOfMass();        // center body on the COM

This will shift the body’s frame to be coincident with the COM, while at the same time translating its mesh vertices in the opposite direction so that its mesh (or meshes) don’t move with respect to world coordinates. The spatial inertia is updated as well.

Alternatively, if the body is being created from a single mesh, one may transform that mesh to be centered on its COM before it is used to define the body. This can be done using the PolygonalMesh method translateToCenterOfVolume(), which centers a mesh’s vertices on its COM (assuming a uniform density):

   PolygonalMesh mesh = new PolygonalMesh (meshFilePath); // read in a mesh
   mesh.translateToCenterOfVolume();        // center mesh on its COM
   RigidBody bodyA = RigidBody.createFromMesh (
      "bodyA", mesh, density, /*scale=*/1); // create body from the mesh

3.2.7 Damping parameters

As with particles, it is possible to set damping parameters for rigid bodies. Damping can be specified in two different ways:

  1. 1.

    Translational/rotational damping which is proportional to a body’s translational and rotational velocity;

  2. 2.

    Inertial damping, which is proportional to a body’s spatial inertia multiplied by its spatial velocity.

Translational/rotational damping is controlled by the MechModel properties frameDamping and rotaryDamping, and generates a spatial force centered on each rigid body’s coordinate frame given by

\hat{\bf f}=-\left(\begin{matrix}d_{f}{\bf v}\\
d_{r}\boldsymbol{\omega}\end{matrix}\right), (3.5)

where d_{f} and d_{r} are the frameDamping and rotaryDamping values, and {\bf v} and \boldsymbol{\omega} are the translational and angular velocity of the body’s coordinate frame. The damping parameters can be set and queried using the MechModel methods

  setFrameDamping (double df)
  setRotaryDamping (double dr)
  double getFrameDamping()
  double getRotaryDamping()

These damping parameters can also be set for individual bodies using their own (inherited) frameDamping and rotaryDamping properties.

For models involving rigid bodies, it is often necessary to set rotaryDamping to a non-zero value because frameDamping will provide no damping at all when a rigid body is simply rotating about its coordinate frame origin.

Inertial damping is controlled by the MechModel property inertialDamping, and generates a spatial force centered on a rigid body’s coordinate frame given by

\hat{\bf f}=-d_{I}\,{\bf M}\,\hat{\bf v},\quad\hat{\bf v}\equiv\left(\begin{%
matrix}{\bf v}\\
\boldsymbol{\omega}\end{matrix}\right), (3.6)

where d_{I} is the inertialDamping, {\bf M} is the body’s 6\times 6 spatial inertia matrix (Section A.6), and \hat{\bf v} is the body’s spatial velocity. The inertial damping property can be set and queried using the MechModel methods

  setInertialDamping (double di)
  double getInertialDamping()

This parameter can also be set for individual bodies using their own (inherited) inertialDamping property.

Inertial damping offers two advantages over translational/rotational damping:

  1. 1.

    It is independent of the location of the body’s coordinate frame with respect to its center of mass;

  2. 2.

    There is no need to adjust two different translational and rotational parameters or to consider their relative sizes, as these considerations are contained within the spatial inertia itself.

3.2.8 Rendering rigid bodies

A rigid body is rendered in ArtiSynth by drawing its mesh (or meshes, Section 3.2.9) and/or coordinate frame.

Meshes are drawn using the face rendering properties described in more detail in Section 4.3. The most commonly used of these are:

  • faceColor: A value of type java.awt.Color giving the color of mesh faces. The default value is GRAY.

  • shading: A value of type Renderer.Shading indicating how the mesh should be shaded, with the options being FLAT, SMOOTH, METAL, and NONE. The default value is FLAT.

  • alpha: A double value between 0 and 1 indicating transparency, with transparency increasing as value decreases from 1. The default value is 1.

  • faceStyle: A value of type Renderer.FaceStyle indicating which face sides should be drawn, with the options being FRONT, BACK, FRONT_AND_BACK, and NONE. The default value is FRONT.

  • drawEdges: A boolean indicating whether the mesh edges should also be drawn, using either the edgeColor rendering property, or the lineColor property if edgeColor is not set. The default value is false.

  • edgeWidth: An integer giving the width of the mesh edges in pixels.

These properties, and others, can be set either interactively in the GUI, or in code. To set the render properties in the GUI, select the rigid body or its mesh component, and then right click the mouse and choose Edit render props .... More details are given in the section “Render properties” in the ArtiSynth User Interface Guide.

Figure 3.5: Different rendering settings for a rigid body hip mesh showing the default (left), smooth rendering with a lighter color (center), and wireframe (right).

Properties can also be set in code, usually during the build() method. Typically this is done using a static method of the RenderProps class that has the form

  RenderProps.setXXX (comp, value)

where XXX is the property name, comp is the component for which the property should be set, and value is the desired value. Some examples are shown in Figure 3.5 for a rigid body hip representation with a fairly coarse mesh. The left image shows the default rendering, using a gray color and flat shading. The center image shows a lighter color and smooth shading, which could be set by the following code fragment:

import maspack.render.*;
import maspack.render.Renderer.*;
  ...
  RigidBody hipBody;
  ...
  RenderProps.setFaceColor (hipBody, new Color (255, 255, 204));
  RenderProps.setShading (hipBody, Shading.SMOOTH);

Finally, the right image shows the body rendered as a wire frame, which can by done by setting faceStyle to NONE and drawEdges to true:

  RenderProps.setFaceStyle (hip, FaceStyle.NONE);
  RenderProps.setDrawEdges (hip, true);
  RenderProps.setEdgeWidth (hip, 2);
  RenderProps.setEdgeColor (hip, Color.CYAN);

Render properties can also be set in higher level model components, from which their values will be inherited by lower level components that have not explicitly set their own values. For example, setting the faceColor render property in the MechModel will automatically set the face color for all subcomponents which have not explicitly set faceColor. More details on render properties are given in Section 4.3.

Figure 3.6: Rigid body axes rendered with axisDrawStyle set to LINE (left) and ARROW (right).

In addition to mesh rendering, it is often useful to draw a rigid body’s coordinate frame, which can be done using its axisLength and axisDrawStyle properties. Setting axisLength to a positive value will cause the body’s three coordinate axes to be drawn, with the indicated length, with the x, y and z axes colored red, green, and blue, respectively. The axisDrawStyle property controls how the axes are rendered (Figure 3.6). It has the type Renderer.AxisDrawStyle, and can be set to the following values:

OFF

Axes are not rendered.

LINE

Axes are rendered as simple red-green-blue lines, with a width given by the joint’s lineWidth rendering property.

ARROW

Axes are rendered as solid red-green-blue arrows.

As with the rendering proprieties, the axisLength and axisDrawStyle properties can be managed either interactively in the GUI (by selecting the body, right clicking and choosing Edit properties ...), or in code, using the following methods:

  double getAxisLength()
  void setAxisLength (double len)
  AxisDrawStyle getAxisDrawStyle()
  void setAxisDrawStyle (AxisDrawStyle style)

3.2.9 Multiple meshes

A RigidBody may contain multiple meshes, which can be useful for various reasons:

  • It may be desirable to use different meshes for collision detection, inertia computation, and visual presentation;

  • Different render properties can be set for different mesh components, allowing the body to be rendered in a more versatile way;

  • Different mesh components can be selected individually.

Each rigid body mesh is encapsulated inside a RigidMeshComp component, which is in turn stored in a subcomponent list called meshes. Meshes do not need to be instances of PolygonalMesh; instead, they can be any instance of MeshBase, including PointMesh and PolylineMesh.

The default surface mesh, returned by getSurfaceMesh(), is also stored inside a RigidMeshComp in the meshes list. By default, the surface mesh is the first mesh in the list, but is otherwise defined to be the first mesh in meshes which is also an instance of PolygonalMesh. The RigidMeshComp containing the surface mesh can be obtained using the method getSurfaceMeshComp().

A RigidMeshComp contains a number of properties that control how the mesh is displayed and interacts with its rigid body:

renderProps

Render properties controlling how the mesh is rendered (see Section 4.3).

hasMass

A boolean, which if true means that the mesh will contribute to the body’s inertia when the inertiaMethod is either MASS or DENSITY. The default value is true.

massDistribution

An enumerated type defined by MassDistribution which specifies how the mesh’s inertia contribution is determined for a given mass. VOLUME, AREA, LENGTH, and POINT indicate, respectively, that the mass is distributed evenly over the mesh’s volume, area (faces), length (edges), or points. The default value is determined by the mesh type: VOLUME for a closed PolygonalMesh, AREA for an open PolygonalMesh, LENGTH for a PolylineMesh, and POINT for a PointMesh. Applications can specify an alternate value providing the mesh has the features to support it. Specifying DEFAULT will restore the default value.

isCollidable

A boolean, which if true, and if the mesh is a PolygonalMesh, means that the mesh will take part in collision and wrapping interactions (Chapter 8 and Section 9.3). The default value is true, and the get/set accessors have the names isCollidable() and setIsCollidable().

volume

A double whose value is the volume of the mesh. If the mesh is a PolygonalMesh, this is the value returned by its computeVolume() method. Otherwise, the volume is 0, unless setVolume(vol) is used to explicitly set a non-zero volume value.

mass

A double whose default value is the product of the density and volume properties. Otherwise, if mass has been explicitly set using setMass(mass), the value is the explicit mass.

density

A double whose default value is the rigid body’s density. Otherwise, if density has been explicitly set using setDensity(density), the value is the explicit density, or if mass has been explicitly set using setMass(mass), the value is the explicit mass divided by volume.

Note that by default, the density of a RigidMeshComp is simply the density setting for the rigid body, and the mass is this times the volume. However, it is possible to set either an explicit mass or a density value that will override this. (Also, explicitly setting a mass will unset any explicit density, and explicitly setting the density will unset any explicit mass.)

When the inertiaMethod of the rigid body is either MASS or DENSITY, then its inertia is computed from the sum of all the inertias {\bf M}_{k} of the component meshes k for which hasMass is true. Each {\bf M}_{k} is computed by the mesh’s createInertia(mass,massDistribution) method, using the mass and massDistribution properties of its RigidMeshComp.

When forming the body inertia from the inertia components of individual meshes, no attempt is made to account for mesh overlap. If this is important, the meshes themselves should be modified in advance so that they do not overlap, perhaps by using the CSG primitives described in Section 2.5.7.

Instances of RigidMeshComp can be created directly, using constructions such as

  PolygonalMesh mesh;
  ... initialize mesh ...
  RigidMeshComp mcomp = new RigidMeshComp (mesh);

or

  RigidMeshComp mcomp = new RigidMeshComp ("meshName");
  mcomp.setMesh (mesh);

after which they can be added or removed from the meshes list using the methods

  void addMeshComp (RigidMeshComp mcomp)
  void addMeshComp (RigidMeshComp mcomp, int idx)
  int numMeshComps()
  boolean removeMeshComp (RigidMeshComp mcomp)
  boolean removeMeshComp (String name)
  void clearMeshComps()

It is also possible to add meshes directly to the meshes list, using the methods

  RigidMeshComp addMesh (MeshBase mesh)
  RigidMeshComp addMesh (MeshBase mesh, boolean hasMass, boolean collidable)

each of which creates a RigidMeshComp, adds it to the mesh list, and returns it. The second method also specifies the values of the hasMass and collidable properties (both of which are true by default).

3.2.10 Example: a composite rigid body

Figure 3.7: RigidCompositeBody loaded into ArtiSynth and run for 0.75 seconds. The ball on the right falls less because it has a lower density than the rest of the body.

An example of constructing a rigid body from multiple meshes is defined in

  artisynth.demos.tutorial.RigidCompositeBody

This uses three meshes to construct a rigid body whose shape resembles a dumbbell. The code, with the include files omitted, is listed below:

1 public class RigidCompositeBody extends RootModel {
2
3    public void build (String[] args) {
4
5       // create MechModel and add to RootModel
6       MechModel mech = new MechModel ("mech");
7       addModel (mech);
8
9       // create the component meshes
10       PolygonalMesh ball1 = MeshFactory.createIcosahedralSphere (0.8, 1);
11       ball1.transform (new RigidTransform3d (1.5, 0, 0));
12       PolygonalMesh ball2 = MeshFactory.createIcosahedralSphere (0.8, 1);
13       ball2.transform (new RigidTransform3d (-1.5, 0, 0));
14       PolygonalMesh axis = MeshFactory.createCylinder (0.2, 2.0, 12);
15       axis.transform (new RigidTransform3d (0, 0, 0, 0, Math.PI/2, 0));
16
17       // create the body and add the component meshes
18       RigidBody body = new RigidBody ("body");
19       body.setDensity (10);
20       body.setFrameDamping (10); // add damping to the body
21       body.addMesh (axis);
22       RigidMeshComp bcomp1 = body.addMesh (ball1);
23       RigidMeshComp bcomp2 = body.addMesh (ball2);
24       mech.addRigidBody (body);
25
26       // connect the body to a spring attached to a fixed particle
27       Particle p1 = new Particle ("p1", /*mass=*/0, /*x,y,z=*/0, 0, 2);
28       p1.setDynamic (false);
29       mech.addParticle (p1);
30       FrameMarker mkr = mech.addFrameMarkerWorld (body, new Point3d (0, 0, 0.2));
31       AxialSpring spring =
32          new AxialSpring ("spr", /*k=*/150, /*d=*/0, /*restLength=*/0);
33       spring.setPoints (p1, mkr);
34       mech.addAxialSpring (spring);
35
36       // set the density for ball1 to be less than the body density
37       bcomp1.setDensity (8);
38
39       // set render properties for the component, with the ball
40       // meshes having different colors
41       RenderProps.setFaceColor (body, new Color (250, 200, 200));
42       RenderProps.setFaceColor (bcomp1, new Color (200, 200, 250));
43       RenderProps.setFaceColor (bcomp2, new Color (200, 250, 200));
44       RenderProps.setSphericalPoints (mech, 0.06, Color.WHITE);
45       RenderProps.setCylindricalLines (spring, 0.02, Color.BLUE);
46    }
47 }

As in the previous examples, the build() method starts by creating a MechModel (lines 6-7). Three different meshes (two balls and an axis) are then constructed at lines 10-15, using MeshFactory methods (Section 2.5) and transforming each result to an appropriate position/orientation with respect to the body’s coordinate frame.

The body itself is constructed at lines 18-24. Its default density is set to 10, and its frame damping (Section 3.2.7) is also set to 10 (the previous rigid body example in Section 3.2.2 relied on spring damping to dissipate energy). The meshes are added using addMesh(), which allocates and returns a RigidMeshComp. For the ball meshes, these are saved in bcomp1 and bcomp2 and used later to adjust density and/or render properties.

Lines 27-34 create a simple linear spring, connected to a fixed point p0 and a marker mkr. The marker is created and attached to the body by the MechModel method addFrameMarkerWorld(), which places the marker at a known position in world coordinates. The spring is created using an AxialSpring constructor that accepts a name, along with stiffness, damping, and rest length parameters to specify a LinearAxialMaterial.

At line 37, bcomp1 is used to set the density of ball1 to 8. Since this is less than the default body density, the inertia component of ball1 will be lighter than that of ball2. Finally, render properties are set at lines 41-45. This includes setting the default face colors for the body and for each ball.

To run this example in ArtiSynth, select All demos > tutorial > RigidCompositeBody from the Models menu. The model should load and initially appear as in Figure 3.7. Running the model (Section 1.5.3) will cause the rigid body to fall and swing about under gravity, with the right ball (ball1) not falling as far because it has less density.

3.3 Joints and connectors

In a typical mechanical model, many of the rigid bodies are interconnected, either using spring-type components that exert binding forces on the bodies, or through joints and connectors that enforce the connection using hard constraints. This section describes the latter. While the discussion focuses on rigid bodies, joints and connectors can be used more generally with any body that implements the ConnectableBody interface. In particular, this allows joints to also interconnect finite element models, as described in Section 6.6.2.

3.3.1 Joints and coordinate frames

Consider two rigid bodies A and B. The pose of body B with respect to body A can be described by the 6 DOF rigid transform {\bf T}_{BA}. If A and B are unconnected, {\bf T}_{BA} may assume any possible value and has a full six degrees of freedom. A joint between A and B constrains the set of poses that are possible between the two bodies and reduces the degrees of freedom available to {\bf T}_{BA}. For ease of use, the constraining action of a joint is described with respect to a pair of local coordinate frames C and D that are connected to frames A and B, respectively, by auxiliary transformations. This allows joints to be placed at locations that do not correspond directly to frames A or B.

Figure 3.8: Coordinate frames D and C for a hinge joint.

The joint frames C and D move with respect to each other as the joint moves. The allowed joint motions therefore correspond to the allowed values of the joint transform {\bf T}_{CD}. Although both frames typically move with their attached bodies, D is considered the base frame and C the motion frame (this is because when a joint is used to connect a single body to ground, body B is set to null and the world frame takes its place). As an example of a joint’s constraining effect, consider a hinge joint (Figure 3.8), which allows C to move with respect to D only by rotating about the z axis while the origins of C and D remain coincident. Other motions are prohibited. If we let \theta describe the counter-clockwise rotation angle of C about the z axis, then {\bf T}_{CD} should always have the form

{\bf T}_{CD}=\left(\begin{matrix}\cos(\theta)&-\sin(\theta)&0&0\\
\sin(\theta)&\cos(\theta)&0&0\\
0&0&1&0\\
0&0&0&1\end{matrix}\right). (3.7)
Figure 3.9: Transforms connecting joint coordinate frames C and D with rigid bodies A and B.

When a joint is attached to bodies A and B, frame C is fixed to body A and frame D is fixed to body B. Except in special cases, the joint frames C and D are not coincident with the body frames A and B. Instead, they are located relative to A and B by the transforms {\bf T}_{CA} and {\bf T}_{DB}, respectively (Figure 3.9). Since {\bf T}_{CA} and {\bf T}_{DB} are both fixed, the joint constraints on {\bf T}_{CD} constrain the relative poses of A and B, with {\bf T}_{AB} determined from

{\bf T}_{AB}={\bf T}_{DB}\,{\bf T}_{CD}\,{\bf T}_{CA}^{-1}. (3.8)

(See Section A.2 for a discussion of determining transforms between related coordinate frames).

3.3.2 Joint coordinates, constraints, and errors

Each different joint and connector type restricts the motion between two bodies to M degrees of freedom, for some M<6. Sometimes, the joint also defines a set of M coordinates that parameterize these M DOFs. For example, the hinge joint described above is parameterized by \theta. Other examples are given in Section 3.4: a 2 DOF cylindrical has coordinates z and \theta, a 3 DOF gimbal joint is parameterized by the roll-pitch-yaw angles \theta, \phi, and \psi, etc. When {\bf T}_{CD}={\bf I} (where {\bf I} is the identity transform), the coordinates are usually all equal to zero, and the joint is said to be in the zero state.

As explained in Section 1.2, ArtiSynth uses a full coordinate formulation for dynamic simulation. That means that instead of using joint coordinates to describe system state, it uses the combined full coordinates {\bf q} of all dynamic components. For example, a model consisting of a single rigid body connected to ground by a hinge joint will have 6 DOF (corresponding to the 6 DOF of the body), rather than the 1 DOF implied by the hinge joint. The DOF restrictions imposed by the joints are then enforced by a set of linearized constraint relationships

{\bf G}({\bf q}){\bf u}={\bf g},\qquad{\bf N}({\bf q}){\bf u}\geq{\bf n} (3.9)

that restrict the body velocities {\bf u} computed at each simulation step, usually by solving an MLCP like (1.6). As explained in Section 1.2, the right side vectors {\bf g} and {\bf n} in (3.9) contain time derivative terms, which for simplicity much of the following presentation will assume to be 0.

Each joint contributes its own set of constraint equations to (3.9). Typically these take the form of bilateral, or equality, constraints

{\bf G}_{\text{J}}({\bf q}){\bf u}=0 (3.10)

which are added to the system’s global bilateral constraint matrix {\bf G}. {\bf G}_{\text{J}} contains 6-M rows providing 6-M individual constraints {\bf G}_{k}. During simulation, these give rise to 6-M constraint forces (corresponding to \boldsymbol{\lambda} in (1.8)) which enforce the constraints.

In some cases, the joint also maintains unilateral, or inequality constraints, to keep {\bf T}_{CD} out of inadmissible regions. These take the form

{\bf N}_{\text{J}}({\bf q}){\bf u}\geq 0 (3.11)

and are added to the system’s global unilateral constraint matrix {\bf N}. They give rise to constraint forces corresponding to \boldsymbol{\theta} in (1.8). A common use of unilateral constraints is to enforce range limits of the joint coordinates (Section 3.3.5), such as

\theta_{\text{min}}\leq\theta\leq\theta_{\text{max}}. (3.12)

A specific unilateral constraint is added to {\bf N}_{\text{J}} only when {\bf T}_{CD} is on or within the boundary of the inadmissible region associated with that constraint. The constraint is then said to be engaged. The combined number of bilateral and engaged unilateral constraints for a particular joint should not exceed 6; otherwise, the joint would be overconstrained.

Joint coordinates, when supported for a particular joint, can be both read and set. Setting a coordinate causes the joint transform {\bf T}_{CD} to change. To accommodate this, the system adjusts the poses of one or both bodies connected to the joint, along with adjacent bodies connected to them, with preference given to bodies that are not attached to “ground”. However, if this is done during simulation, and particularly if one or both of the bodies connected to the joint are moving dynamically, the results will be unpredictable and will likely conflict with the simulation.

Joint coordinates are also often exported as properties. For example, the HingeJoint class (Section 3.4) exports its \theta coordinate as the property theta, which can be accessed in the GUI, or via the accessor methods

  double getTheta()          // get theta in degrees
  void setTheta (double deg) // set theta in degrees

Since joint constraints are generally nonlinear, their linearized enforcement at the velocity level by (3.9) will usually produce small errors as the simulation proceeds. These errors are reduced using a position correction step described in Section 4.9.1 and [11]. Errors can also be caused by joint compliance (Section 3.3.8). Both effects mean that the joint transform {\bf T}_{CD} may deviate from the allowed values dictated by the joint type. In ArtiSynth, this is accounted for by introducing an additional constraint frame G between D and C (Figure 3.10). G is computed to be the nearest frame to C that lies exactly in the joint constraint space. {\bf T}_{GD} is therefore a valid joint transform, {\bf T}_{GC} accommodates the error, and the whole joint transform is given by the composition

{\bf T}_{CD}={\bf T}_{GD}\;{\bf T}_{CG}. (3.13)

If there is no compliance or joint error, then frames G and C are identical, {\bf T}_{CG}={\bf I}, and {\bf T}_{CD}={\bf T}_{GD}. Because {\bf T}_{CG} describes the joint error, we sometimes refer to it as {\bf T}_{err}={\bf T}_{CG}.

Figure 3.10: 2D schematic showing the joint frames D and C, along with the intermediate frame G that accounts for numeric error and complaint motion.

3.3.3 Creating joints

Joint and connector components in ArtiSynth are both derived from the superclass BodyConnector, with joints being further derived from JointBase, which provides support for coordinates. Some of the commonly used joints and connectors are described in Section 3.4.

An application creates a joint by constructing it and adding it to a MechModel. Many joints have constructors of the form

  XXXJoint (bodyA, bodyB, TDW)

which specifies the bodies A and B which the joint connects, along with the transform {\bf T}_{DW} giving the pose of the joint base frame D in world coordinates. The constructor then assumes that the joint is in the zero state, so that C and D are the same and {\bf T}_{CD}={\bf I} and {\bf T}_{CW}={\bf T}_{DW}, and then computes {\bf T}_{CA} and {\bf T}_{DB} from

\displaystyle{\bf T}_{CA} \displaystyle={\bf T}_{AW}^{-1}\;{\bf T}_{CW} (3.14)
\displaystyle{\bf T}_{DB} \displaystyle={\bf T}_{BW}^{-1}\;{\bf T}_{DW}, (3.15)

where {\bf T}_{AW} and {\bf T}_{BW} are the current poses of A and B.

After the joint is created, it should be added to the system’s MechModel using addBodyConnector(), as shown in the following code fragment:

   MechModel mech;
   RigidBody bodyA, bodyB;
   RigidTransform3d TDW;
   ... initialize mech, bodyA, bodyB, and TDW ...
   HingeJoint joint = new HingeJoint (bodyA, bodyB, TDW);
   mech.addBodyConnector (joint);

It is also possible to create a joint using its default constructor and attach it to the bodies afterward, using the method setBodies(bodyA,bodyB,TDW), as in the following:

   HingeJoint joint = new HingeJoint ();
   joint.setBodies (bodyA, bodyB, TDW);
   mech.addBodyConnector (joint);

One reason for doing this is that it allows the joint transform {\bf T}_{CD} to be modified (by setting coordinate values) before setBodies() is called; this is discussed further in Section 3.3.4.

Joints usually offer a number of other constructors that let its world location and body relationships to be specified in different ways. These may include:

  XXXJoint (bodyA, TCA, bodyB, TDB)
  XXXJoint (bodyA, bodyB, TCW, TDW)

The first, which is restricted to rigid bodies, allows the application to explicitly specify transforms {\bf T}_{CA} and {\bf T}_{DB} connecting frames C and D to the body frames A and B, and is useful when {\bf T}_{CA} and {\bf T}_{DB} are explicitly known, or the initial value of {\bf T}_{CD} is not the identity. Likewise, the second constructor allows {\bf T}_{CW} and {\bf T}_{DW} to be explicitly specified, with {\bf T}_{CD}\neq{\bf I} if {\bf T}_{CW}\neq{\bf T}_{DW}. For instance, suppose {\bf T}_{CD} and {\bf T}_{DW} are both known. Then we can use the relationship

{\bf T}_{CW}={\bf T}_{DW}\,{\bf T}_{CD} (3.16)

to create the joint as in the following code fragment:

   MechModel mech;
   RigidBody bodyA, bodyB;
   RigidTransform3d TDW, TCD;
   ... initialize mech, bodyA, bodyB, TDW, and TCD ...
   // compute TCW:
   RigidTransform3d TCW = new RigidTransform3d();
   TCW.mul (TDW, TCD);
   HingeJoint joint = new HingeJoint (bodyA, bodyB, TCW, TDW);
   mech.addBodyConnector (joint);

As an alternative to specifying {\bf T}_{DW} or its equivalents, some joint types provide constructors that let the application locate specific joint features. These may be easier to use in some cases. For instance, HingeJoint provides a constructor

   HingeJoint (bodyA, bodyB, originD, zaxis)

that specifies origin of D and its z axis (which is the rotation axis), with the remaining orientation of D aligned as closely as possible with the world. SphericalJoint provides a constructor

   SphericalJoint (bodyA, bodyB, originD)

that specifies origin of D and aligns its orientation with the world. Users should consult the source code or API documentation for specific joints to see what special constructors may be available.

Finally, it is possible to use joints to connect a single body to ground (by convention, this is the A body). Most joints provide a constructor of the form

  XXXJoint (bodyA, TDW)

which allows this to be done explicitly. Alternatively, most joint constructors which supply body B will allow this to be specified as null, so that body A will be connected to ground by default.

3.3.4 Working with coordinates

As mentioned in Section 3.3.2, some joints support coordinates that parameterize the valid motions within the joint transform {\bf T}_{CD}. All such joints are subclasses of JointBase, which provides some generic methods for querying and setting coordinate values (JointBase is in turn a subclass of BodyConnector).

The number of coordinates is returned by the method numCoordinates(); if this returns 0, then coordinates are not supported. Each coordinate has an index in the range 0\ldots M-1, where M is the number of coordinates. Coordinate values can be queried or set using the following methods:

  getCoordinate (int idx)               // get coordinate value with index idx
  getCoordinates (VectorNd coords)      // get all coordinates values
  setCoordinate (int idx, double value) // set coordinate value with index idx
  setCoordinates (VectorNd values)      // set all coordinates values

Specific joint types usually also provide names for their joint coordinates, along with integer constants describing their indices and methods for accessing their values. For example, CylindricalJoint supports two coordinates, z and \theta, along with the following:

  // coordinate indices
  static final int Z_IDX = 0;
  static final int THETA_IDX = 1;
  // set/get z value and range
  double getZ()
  void setZ (double z)
  // set/get theta value and range in degrees
  double getTheta()
  void setTheta (double theta)

The coordinate values are also exported as the properties z and theta, allowing them to be set in the GUI. For convenience, particularly in GUI applications, the properties and methods for controlling specific angular coordinates generally use degrees instead of radians.

As discussed in Section 3.3.2, unlike in some multibody simulation systems (such as OpenSim), joint coordinates are not fundamental quantities that describe system state. As such, then, coordinates can usually only be set in specific circumstances that avoid simulation conflicts. In general, when joint coordinates are set, the system adjusts the poses of one or both bodies connected to this joint, along with adjacent bodies connected to them, with preference given to bodies that are not attached to “ground”. However, if this is done during simulation, and particularly if one or both of the bodies connected to the joint are moving dynamically, the results will be unpredictable and will likely conflict with the simulation.

If a joint has been created with its default constructor and not yet attached to any bodies, then setting joint values will simply set the joint transform {\bf T}_{CD}. This can be useful in situations where one needs to initialize a joint’s {\bf T}_{CD} to a non-identity value corresponding to a particular set of joint coordinates:

  RigidTransform3d TDW; // known location for D frame
  double z, theta; // desired initial coordinate values
  ...
  CylindricalJoint joint = new CylindricalJoint();
  joint.setZ (z);
  joint.setTheta (thetaDeg);
  joint.setBodies (bodyA, bodyB, TDW);

This can also be done in vector form:

  RigidTransform3d TDW; // known location for D frame
  VectorNd coordValues; // desired initial coordinate values
  ...
  CylindricalJoint joint = new CylindricalJoint();
  joint.setCoordinates (coordValues);
  joint.setBodies (bodyA, bodyB, TDW);

In either of these cases, setBodies() will not use {\bf T}_{CD}={\bf I} but instead use the value determined by the initial coordinate values.

To determine the {\bf T}_{CD} corresponding to a particular set of coordinates, one may use the method

  void coordinatesToTCD (RigidTransform3d TCD, VectorNd coords)

In some cases, within a model’s build() method, one may wish to set initial coordinates after a joint has been attached to its bodies, in order to move those bodies (along with the bodies attached to them) into an initial configuration without having to explicitly calculate the poses from the joint coordinates. As mentioned above, the system will make a decision about which attached bodies are most “free” and adjust their poses accordingly. This is done in the example of the next section.

3.3.5 Coordinate limits and locking

It is possible to set limits on a joint coordinate’s range, and also to lock a coordinate in place at its current value.

When a joint coordinate hits either an upper or lower range limit, a unilateral constraint is invoked to prevent it from violating the limit, and remains engaged until the joint moves away from the limit. Each range constraint that is engaged reduces the number of joint DOFs by one.

By default, joint range limits are usually disabled (i.e., they are set to (-\inf,\inf)). They can be queried and set, for a given joint with index idx, using the methods:

  DoubleInterval getCoordinateRange (int idx)
  double getMinCoordinate (int idx)
  double getMaxCoordinate (int idx)
  void setCoordinateRange (idx, DoubleInterval rng)

where range limits for angular coordinates are specified in radians. For convenience, the following methods are also provided which use degrees instead of radians for angular coordinates:

  DoubleInterval getCoordinateRangeDeg (int idx)
  double getMinCoordinateDeg (int idx)
  double getMaxCoordinateDeg (int idx)
  void setCoordinateRangeDeg (idx, DoubleInterval rng)

Range checking can be disabled by setting the range to (-\inf,\inf), or by specifying rng as null, which implicitly does the same thing.

Ranges for angular coordinates are not limited to \pm 180^{\circ} but can instead be set to larger values; the joint will continue to wrap until the limit is reached.

Joint coordinates can also be locked, so that they hold their current value and don’t move. A joint is locked using a bilateral constraint that prevents motion in either direction and reduces the joint’s DOF count by one. The following methods are available for querying or setting a coordinate’s locked status:

  boolean isLocked (int idx)
  void setLocked (int idx, boolean locked)

As with coordinate values, specific joint types usually provide methods for controlling the ranges and locking status of individual coordinates, with ranges for angular coordinates specified in degrees instead of radians. For example, CylindricalJoint supplies the methods

  // set/get z range
  DoubleInterval getZRange()
  void setZRange (double min, double max)
  // control z locking
  boolean isZLocked()
  void setZLocked (boolean locked)
  // set/get theta range in degrees
  DoubleInterval getThetaRange()
  void setThetaRange (double min, double max)
  void setThetaRange (DoubleInterval rng)
  // control theta locking
  boolean isThetaLocked()
  void setThetaLocked (boolean locked)

The range and locking information is also exported as the properties zRange, thetaRange, zLocked, and thetaLocked, allowing them to be set in the GUI.

3.3.6 Example: a simple hinge joint

Figure 3.11: RigidBodyJoint model loaded into ArtiSynth.

A simple model showing two rigid bodies connected by a joint is defined in

  artisynth.demos.tutorial.RigidBodyJoint

The build method for this model is given below:

1    public void build (String[] args) {
2
3       // create MechModel and add to RootModel
4       mech = new MechModel ("mech");
5       mech.setGravity (0, 0, -98);
6       mech.setFrameDamping (1.0);
7       mech.setRotaryDamping (4.0);
8       addModel (mech);
9
10       PolygonalMesh mesh;  // bodies will be defined using a mesh
11
12       // create first body and set its pose
13       mesh = MeshFactory.createRoundedBox (lenx1, leny1, lenz1, /*nslices=*/8);
14       RigidTransform3d TMB =
15          new RigidTransform3d (0, 0, 0, /*axisAng=*/1, 1, 1, 2*Math.PI/3);
16       mesh.transform (TMB);
17       bodyB = RigidBody.createFromMesh ("bodyB", mesh, /*density=*/0.2, 1.0);
18       bodyB.setPose (new RigidTransform3d (0, 0, 1.5*lenx1, 1, 0, 0, Math.PI/2));
19       bodyB.setDynamic (false);
20
21       // create second body and set its pose
22       mesh = MeshFactory.createRoundedCylinder (
23          leny2/2, lenx2, /*nslices=*/16, /*nsegs=*/1, /*flatBottom=*/false);
24       mesh.transform (TMB);
25       bodyA = RigidBody.createFromMesh ("bodyA", mesh, 0.2, 1.0);
26       bodyA.setPose (new RigidTransform3d (
27                         (lenx1+lenx2)/2, 0, 1.5*lenx1, 1, 0, 0, Math.PI/2));
28
29       // create the joint
30       RigidTransform3d TDW =
31          new RigidTransform3d (lenx1/2, 0, 1.5*lenx1, 1, 0, 0, Math.PI/2);
32       HingeJoint joint = new HingeJoint (bodyA, bodyB, TDW);
33
34       // add components to the mech model
35       mech.addRigidBody (bodyB);
36       mech.addRigidBody (bodyA);
37       mech.addBodyConnector (joint);
38
39       joint.setTheta (35);  // set joint position
40
41       // set render properties for components
42       RenderProps.setFaceColor (joint, Color.BLUE);
43       joint.setShaftLength (4);
44       joint.setShaftRadius (0.2);
45    }

A MechModel is created as usual at line 4. However, in this example, we also set some parameters for it: setGravity() is used to set the gravity acceleration vector to (0,0,-98)^{T} instead of the default value of (0,0,-9.8)^{T}, and the frameDamping and rotaryDamping properties (Section 3.2.7) are set to provide appropriate damping.

Each of the two rigid bodies are created from a mesh and a density. The meshes themselves are created using the factory methods MeshFactory.createRoundedBox() and MeshFactory.createRoundedCylinder() (lines 13 and 22), and then RigidBody.createFromMesh() is used to turn these into rigid bodies with a density of 0.2 (lines 17 and 25). The pose of the two bodies is set using RigidTransform3d objects created with x, y, z translation and axis-angle orientation values (lines 18 and 26).

The hinge joint is implemented using HingeJoint, which is constructed at line 32 with the joint coordinate frame D being located in world coordinates by TDW as described in Section 3.3.3.

Once the joint is created and added to the MechModel, the method setTheta() is used to explicitly set the joint parameter to 35 degrees. The joint transform {\bf T}_{CD} is then set appropriately and bodyA is moved to accommodate this (bodyA being chosen since it is the most free to move).

Finally, joint rendering properties are set starting at line 42. We render the joint as a cylindrical shaft about the rotation axis, using its shaftLength and shaftRadius properties. Joint rendering is discussed in more detail in Section 3.3.10).

To run this example in ArtiSynth, select All demos > tutorial > RigidBodyJoint from the Models menu. The model should load and initially appear as in Figure 3.11. Running the model (Section 1.5.3) will cause bodyA to fall and swing under gravity.

3.3.7 Constraint forces

During each simulation solve step, the joint velocity constraints described by (3.10) and (3.11) are enforced by bilateral and unilateral constraint forces {\bf f}_{g} and {\bf f}_{n}:

{\bf f}_{g}={\bf G}_{\text{J}}^{T}\boldsymbol{\lambda}_{J},\quad{\bf f}_{n}={%
\bf N}_{\text{J}}^{T}\boldsymbol{\theta}_{J}. (3.17)

Here, {\bf f}_{g} and {\bf f}_{n} are spatial forces (or wrenches, Section A.5) acting in the joint coordinate frame C, and \boldsymbol{\lambda}_{J} and \boldsymbol{\theta}_{J} are the Lagrange multipliers computed as part of the mechanical system solve (see (1.6) and (1.8)). The sizes of \boldsymbol{\lambda}_{J} and \boldsymbol{\theta}_{J} equal the number of bilateral and engaged unilateral constraints in the joint; these numbers can be queried for a particular joint using the methods numBilateralConstraints() and numEngagedUnilateralConstraints(). (The number of engaged unilateral constraints may be less than the total number of unilateral constraints; the latter may be queried with numUnilateralConstraints(), while the total number of constraints is returned by numConstraints().

Applications may sometimes need to query the current constraint force values, typically from within a controller or monitor (Section 5.3). The Lagrange multipliers themselves may be obtained with

  void getBilateralForces (VectorNd lam)
  void getUnilateralForces (VectorNd the)

which load the multipliers into lam or the and set their sizes to the number of bilateral or engaged unilateral constraints. Alternatively, one can retrieve the individual multiplier for the constraint indexed by idx using

  double getConstraintForce (int idx);

Typically, it is more useful to find the spatial constraint forces {\bf f}_{g} and {\bf f}_{n}, which can be obtained with respect to frame C:

  // place the forces in the wrench f
  void getBilateralForcesInC (Wrench f)
  void getUnilateralForcesInC (Wrench f)
  // convenience methods that allocate the wrench and return it
  Wrench getBilateralForcesInC();
  Wrench getUnilateralForcesInC();

If the attached bodies A and B are rigid bodies, it is also possible to obtain the constraint wrenches experienced by those bodies:

  // place the forces in the wrench f
  void getBilateralForcesInA (Wrench f)
  void getUnilateralForcesInA (Wrench f)
  void getBilateralForcesInB (Wrench f)
  void getUnilateralForcesInB (Wrench f)
  // convenience methods that allocate the wrench and return it
  Wrench getBilateralForcesInA();
  Wrench getUnilateralForcesInA();
  Wrench getBilateralForcesInB();
  Wrench getUnilateralForcesInB();

Constraint wrenches obtained for bodies A or B are given in world coordinates, which is consistent with the forces reported by rigid bodies via their getForce() method. To orient the forces into body coordinates, one may use the inverse of the rotation matrix {\bf R} of the body’s pose. For example:

   RigidBody bodyA;
   // ... body A initialized, etc. ...
   Wrench force = joint.getBilateralForceInA();
   force.inverseTransform (bodyA.getPose().R);

3.3.8 Compliance and regularization

By default, the constraints used to implement joints and couplings are treated as hard, so that the system tries to respect the constraint conditions (3.9) as exactly as possible as the simulation proceeds. Sometimes, however, it is desirable to introduce some “softness” into the constraints, whereby constraint forces are determined as a linear function of their distance from the constraint. Adding compliance also allows an application to regularize a system of joint constraints that would otherwise be overconstrained, as illustrated in Section 3.3.9.

To describe compliance precisely, consider the bilateral constraint portion of the MLCP in (1.6), which solves for the updated system velocities {\bf u}^{k+1} at each time step:

\left(\begin{matrix}\hat{\bf M}^{k}&-{\bf G}^{T}\\
{\bf G}&0\end{matrix}\right)\left(\begin{matrix}{\bf u}^{k+1}\\
\tilde{\boldsymbol{\lambda}}\end{matrix}\right)=\left(\begin{matrix}{\bf M}{%
\bf u}^{k}-h\hat{\bf f}^{k}\\
0\end{matrix}\right). (3.18)

Here {\bf G} is the system’s bilateral constraint matrix, \tilde{\boldsymbol{\lambda}} denotes the constraint impulses (from which the constraint forces \boldsymbol{\lambda} can be determined by \boldsymbol{\lambda}=\tilde{\boldsymbol{\lambda}}/h), and for simplicity we have assumed that {\bf G} is constant and so the {\bf g} term on the lower right side is 0.

Solving (3.18) results in constraint forces that satisfy {\bf G}{\bf u}^{k+1}=0 precisely, corresponding to hard constraints. To implement soft constraints, start by defining a function \boldsymbol{\phi}({\bf q}) that defines the distances from each constraint, where {\bf q} is the vector of system positions; these distances are the local translational and rotational deviations from each constraint’s correct position and are discussed in more detail in Section 4.9.1. Then assume that the constraint forces are a linear function of these distances:

\boldsymbol{\lambda}=-{\bf C}^{-1}\boldsymbol{\phi}({\bf q}), (3.19)

where {\bf C} is a diagonal compliance matrix that is equivalent to an inverse stiffness matrix. We also note that \boldsymbol{\phi} will be time varying, and that we can approximate its change between time steps as

\boldsymbol{\phi}^{k+1}\approx\boldsymbol{\phi}^{k}+h\dot{\boldsymbol{\phi}}^{%
k+1},\quad\text{with}\quad\dot{\boldsymbol{\phi}}^{k+1}={\bf G}{\bf u}^{k+1}. (3.20)

Next, assume that in using (3.19) to determine \boldsymbol{\lambda} for a particular time step, we use the average value of \boldsymbol{\phi} over the step, represented by \bar{\boldsymbol{\phi}}=(\boldsymbol{\phi}^{k+1}+\boldsymbol{\phi}^{k})/2. Substituting this and (3.20) into (3.19), multiplying by {\bf C}, and rearranging yields:

{\bf G}{\bf u}^{k+1}+\frac{2{\bf C}}{h}\boldsymbol{\lambda}=-\frac{2}{h}%
\boldsymbol{\phi}^{k}. (3.21)

Then noting that \tilde{\boldsymbol{\lambda}}=h\boldsymbol{\lambda}, we obtain a revised form of (3.18),

\left(\begin{matrix}\hat{\bf M}^{k}&-{\bf G}^{T}\\
{\bf G}&2{\bf C}/h^{2}\end{matrix}\right)\left(\begin{matrix}{\bf u}^{k+1}\\
\tilde{\boldsymbol{\lambda}}\end{matrix}\right)=\left(\begin{matrix}{\bf M}{%
\bf u}^{k}-h\hat{\bf f}^{k}\\
-2\boldsymbol{\phi}^{k}/h\end{matrix}\right), (3.22)

in the which the zeros in the matrix and right hand side have been replaced by compliance terms. The resulting constraint behavior is different from that of (3.18) in two important ways:

  1. 1.

    The joint now allows 6 DOF, with motion along the constrained directions limited by restoring spring constants given by the reciprocals of the diagonal entries of {\bf C}.

  2. 2.

    If {\bf C} has no zero diagonal entries, then the system (3.22) is regularized by the 2{\bf C}/h^{2} term in the lower right matrix block. This means that the matrix is always non-singular, even if {\bf G} is rank deficient, and so compliance offers a way to handle overconstrained models, as discussed further in Section 3.3.9.

Unilateral constraints can be regularized using the same approach, with a distance function defined such that \boldsymbol{\phi}({\bf q})\leq 0.

The reason for specifying soft constraints using compliance instead of stiffness is that by setting {\bf C}=0 we can easily handle the case of infinite stiffness where the constraints are strictly enforced. The ArtiSynth compliance implementation uses a slightly more complex version of (3.22) that accounts for non-constant {\bf G} and also allows for a damping term -{\bf D}\dot{\boldsymbol{\phi}}, where {\bf D} is again a diagonal matrix. For more details, see [9] and [21].

When using compliance, damping is often needed for stability, and, in the case of unilateral constraints, to prevent “bouncing”. A good choice for damping is usually critical damping, which is discussed further below.

Any joint which is a subclass of BodyConnector allows individual compliance values C_{i} and damping values D_{i} to be set for each of the joint’s i constraints. These values comprise the diagonal entries in the compliance and damping matrices {\bf C} and {\bf D}, and can be queried and set using the methods

  VectorNd getCompliance()
  void setCompliance (VectorNd compliance)
  VectorNd getDamping()
  void setCompliance (VectorNd damping)

The vectors supplied to the above set methods contain the requested compliance or damping values. If their size n is less than numConstraints(), then compliance or damping will be set for the first n constraints. Damping for a specific constraint only has an effect if the compliance for that constraint is nonzero.

What compliance and damping values should be specified? Compliance is usually relatively easy to figure out. Each of the joint’s individual constraints i corresponds to a row in its bilateral constraint matrix {\bf G}_{\text{J}} or unilateral constraint matrix {\bf N}_{\text{J}}, and represents a specific 6 DOF direction along which the spatial velocity \hat{\bf v}_{CD} (of frame C with respect to D) is restricted (more details on this are given in Section 4.9.1). Each of these constraint directions is usually predominantly linear or rotational; specific descriptions for the constraints of different joints are provided in Section 3.4. To determine compliance for a constraint i, estimate the typical force f likely to act along its direction, decide how much displacement \delta q (translational or rotational) along that constraint is desirable, and then set the compliance C_{i} to the associated inverse stiffness:

C_{k}=\delta q/f. (3.23)

Once C_{k} is determined, the damping D_{k} can be estimated based on the desired damping ratio \zeta, using the formula

D_{k}=2\zeta\sqrt{M/C_{k}} (3.24)

where M is total mass of the bodies attached to the joint. Typically, the desired damping will be close to critical damping, for which \zeta=1.

Constraints associated with linear motion will typically require different compliance values from those associated with rotation. To make this process easier, joint components allow the setting of collective compliance values for their linear and rotary constraints, using the methods

  void setLinearCompliance (double c)
  double getLinearCompliance()
  void setRotaryCompliance (double c)
  double getRotaryCompliance()

The set() methods will set a uniform compliance for all linear or rotary constraints, except for unilateral constraints associated with coordinate limits. At the same time, they will also set an automatically computed critical damping value. Likewise, the get() methods query these linear or rotary constraints for uniform compliance values (with the corresponding critical damping), and return either that value, or -1 if it does not exist.

Most of the demonstration models for the joints described in Section 3.4 allow these linear and rotary compliance settings to be adjusted interactively using a control panel, enabling users to experimentally gain a feel for their behavior.

To determine programmatically whether a particular constraint is linear or rotary, one can use the joint method

  VectorNi getConstraintFlags()

which returns a vector of information flags for all its constraints. Linear and rotary constraints are indicated by the flags LINEAR and ROTARY, defined in RigidBodyConstraint.

3.3.9 Example: an overconstrained linkage

Situations may occasionally arise in which a model is overconstrained, which means that the rows of the bilateral constraint matrix {\bf G} in (3.9) are not all linearly dependent, or in other words, {\bf G} does not have full row rank. At present, the ArtiSynth solver has difficultly handling overconstrained models, but these situations can often be handled by adding a small amount of compliance to the constraints. (Overconstraining is not a problem with unilateral constraints {\bf N}, because of the way they are handled by the solver.)

One possible symptom of an overconstrained system is a error message in the application’s terminal output, such as

Pardiso: num perturbed pivots=12

Overconstraining frequently occurs in closed-chain linkages, involving loops in which a jointed sequence of links is connected back on itself. Depending on how the constraints are configured and how redundant they are, the system may still be able to move. A classical example is the four-bar linkage, a common version of which consists of four links, or “bars”, arranged as a parallelogram and connected by hinge joints at the corners. One link is usually connected to ground, and so the remaining three links together have 18 DOF, while the four hinge joints together remove 20 DOF, overconstraining the system. However, the constraints are redundant in such as way that the linkage still actually has 1 DOF.

Figure 3.12: FourBarLinkage model, several steps into the simulation.

To model a four-bar in ArtiSynth presently requires adding compliance to the hinge joints. An example of this is defined by the demo program

  artisynth.demos.tutorial.FourBarLinkage

shown in Figure 3.12. The code for the build() method and a couple of supporting methods is given below:

1    /**
2     * Create a link with a length of 1.0, width of 0.25, and specified depth
3     * and add it to the mech model. The parameters x, z, and deg specify the
4     * link’s position and orientation (in degrees) in the x-z plane.
5     */
6    protected RigidBody createLink (
7       MechModel mech, String name,
8       double depth, double x, double z, double deg) {
9       int nslices = 20; // num slices on the rounded mesh ends
10       PolygonalMesh mesh =
11          MeshFactory.createRoundedBox (1.0, 0.25, depth, nslices);
12       RigidBody body = RigidBody.createFromMesh (
13          name, mesh, /*density=*/1000.0, /*scale=*/1.0);
14       body.setPose (new RigidTransform3d (x, 0, z, 0, Math.toRadians(deg), 0));
15       mech.addRigidBody (body);
16       return body;
17    }
18
19    /**
20     * Create a hinge joint connecting one end of link0 with the other end of
21     * link1, and add it to the mech model.
22     */
23    protected HingeJoint createJoint (
24       MechModel mech, String name, RigidBody link0, RigidBody link1) {
25       // easier to locate the link using TCA and TDB since we know where frames
26       // C and D are with respect the link0 and link1
27       RigidTransform3d TCA = new RigidTransform3d (0, 0,  0.5, 0, 0, Math.PI/2);
28       RigidTransform3d TDB = new RigidTransform3d (0, 0, -0.5, 0, 0, Math.PI/2);
29       HingeJoint joint = new HingeJoint (link0, TCA, link1, TDB);
30       joint.setName (name);
31       mech.addBodyConnector (joint);
32       // set joint render properties
33       joint.setAxisLength (0.4);
34       RenderProps.setLineRadius (joint, 0.03);
35       return joint;
36    }
37
38    public void build (String[] args) {
39       // create a mech model and set rigid body damping parameters
40       MechModel mech = new MechModel ("mech");
41       addModel (mech);
42       mech.setFrameDamping (1.0);
43       mech.setRotaryDamping (4.0);
44
45       // create four ’bars’ from which to construct the linkage
46       RigidBody[] bars = new RigidBody[4];
47       bars[0] = createLink (mech, "link0", 0.2, -0.5,  0.0, 0);
48       bars[1] = createLink (mech, "link1", 0.3,  0.0,  0.5, 90);
49       bars[2] = createLink (mech, "link2", 0.2,  0.5,  0.0, 180);
50       bars[3] = createLink (mech, "link3", 0.3,  0.0, -0.5, 270);
51       // ground the left bar
52       bars[0].setDynamic (false);
53
54       // connect the bars using four hinge joints
55       HingeJoint[] joints = new HingeJoint[4];
56       joints[0] = createJoint (mech, "joint0", bars[0], bars[1]);
57       joints[1] = createJoint (mech, "joint1", bars[1], bars[2]);
58       joints[2] = createJoint (mech, "joint2", bars[2], bars[3]);
59       joints[3] = createJoint (mech, "joint3", bars[3], bars[0]);
60
61       // Set uniform compliance and damping for all bilateral constraints,
62       // which are the first 5 constraints of each joint
63       VectorNd compliance = new VectorNd(5);
64       VectorNd damping = new VectorNd(5);
65       for (int i=0; i<5; i++) {
66          compliance.set (i, 0.000001);
67          damping.set (i, 25000);
68       }
69       for (int i=0; i<joints.length; i++) {
70          joints[i].setCompliance (compliance);
71          joints[i].setDamping (damping);
72       }
73    }

Two helper methods are used to construct the model: createLink() (lines 6-17), and createJoint() (lines 23-36). createLink() makes the individual rigid bodies used to build the linkage: a mesh is produced defining the body’s shape (a box with rounded ends), and then passed to the RigidBody createFromMesh() method which creates the body and sets its inertia according to a specified density. The body’s pose is then set so as to center it at (x,0,z) while rotating it about the y axis by the angle deg (in degrees). The completed body is then added to the MechModel mech and returned.

The second helper method, createJoint(), connects two rigid bodies (link0 and link1) together using a HingeJoint. Because we know the location of the joint in body-relative coordinates, it is easier to create the joint using the transforms {\bf T}_{CA} and {\bf T}_{DB} instead of {\bf T}_{DW}: {\bf T}_{CA} locates the joint at the top end of link0, at (0,0,0.5), with the z axis parallel to the body’s y axis, while {\bf T}_{DB} similarly locates the joint at the bottom of link1. After the joint is created and added to the MechModel, its render properties are set so that its axis drawn as a blue cylinder.

The build() method itself begins by creating a MechModel and setting damping parameters for the rigid bodies (lines 40-43). Next, createLink() is used to create and store the four links (lines 46-50), and the left bar is attached to ground by making it non-dynamic (line 52). The links are then connected together using joints created by createJoint() (lines 55-59). Finally, uniform compliance and damping values are set for each of the joint’s bilateral constraints, using the setCompliance() and setDamping() methods (lines 63-72). Values are set for the first five constraints, since for a HingeJoint these are the bilateral constraints. The compliance value of C=10^{-6} was found experimentally to be low enough so as to not cause noticeable deflections in the joints. Given C and an average mass of around M=150 for each link pair, (3.24) suggests the damping factor of D=25000. Note that for this example, very similar settings could be achieved by simply calling

  for (int i=0; i<joints.length; i++) {
     joints[i].setLinearCompliance (0.000001);
     joints[i].setRotaryCompliance (0.000001);
  }

In principle, we only need to set compliance for the constraints that are redundant, but it can sometimes be difficult to determine exactly which these are. Also, different values are often needed for linear and rotary constraints; that is not necessary here because the links have unit length and so the linear and rotary units have similar scales.

3.3.10 Rendering joints

Most joints provide a means to render themselves in order to provide a graphical representation of their position and configuration. Control over this is achieved by setting various properties in the joint component, including both specialized properties and the standard render properties (Section 4.3) used by all renderable components.

All joints which are subclasses of JointBase support rendering of both their C and D coordinate frames, through the properties drawFrameC, drawFrameD, and axisLength. The first two properties are of the type Renderer.AxisDrawStyle (described in detail in Section 3.2.8), and can be set to LINE or ARROW to enable the coordinate axes to be drawn either as lines or solid arrows. The axisLength property has type double and specifies the length with which the axes are drawn. As with all properties, these properties can be set either in the GUI, or in code using accessor methods supplied by the joint:

  void setAxisLength (double l)
  double getAxisLength()
  void setDrawFrameC (AxisDrawStyle style)
  (AxisDrawStyle getDrawFrameC()
  void setDrawFrameD (AxisDrawStyle style)
  (AxisDrawStyle getDrawFrameD()

Another pair of properties used by several joints is shaftLength and shaftRadius, which specify the length and radius used to draw shaft or axis structures associated with the joint. These are rendered as solid cylinders, using the color indicated by the faceColor rendering property. The default value of both properties is 0; if shaftLength is 0, then the structures are not drawn, while if shaftRadius is 0, a default value proportional to shaftLength is used. For example, to enable rendering of a blue shaft along the rotation axis of a hinge joint, one may use the code fragment

  HingeJoint joint;
  ...
  joint.setShaftLength (0.5); // set shaft dimensions
  joint.setShaftRadius (0.05);
  RenderProps.setFaceColor (joint, Color.BLUE); // set the color

As another example, to enable rendering of a green ball about the center of a spherical joint, one may use the fragment

  SphericalJoint joint;
  ...
  joint.setJointRadius (0.02); // set the ball size
  RenderProps.setFaceColor (joint, Color.GREEN); // set the color

Specific joints may define additional properties to control how they are rendered.

3.4 Joint components

ArtiSynth supplies a number of basic joints and connectors in the package artisynth.core.mechmodels, the most common of which are described here.

Many of the descriptions are associated with a demonstration model, named XXXJointDemo, where XXX is the joint type. These demos are located in the package artisynth.demos.mech, and can be loaded by selecting All demos > mech > XXXJointDemo from the Models menu. When run, they can be interactively controlled, using either the pull tool (see the section “Pull Manipulation” in the ArtiSynth User Interface Guide), or the interactive control panel. The control panel allows the adjustment of coordinate values and ranges (if supported), some of the render properties, and the different compliance and damping properties (Section 3.3.8). One can inspect the source code for each demo in its .java file located in the folder <ARTISYNTH_HOME>/src/artisynth/demos/mech.

3.4.1 Hinge joint

Figure 3.13: Coordinate frames (left) and demo model (right) for the hinge joint.

The HingeJoint (Figure 3.13) is a 1 DOF joint that constrains motion between frames C and D to a simple rotation about the z axis of D. It implements six constraints and one coordinate \theta (Table 3.1), to which the joint transform {\bf T}_{CD} is related by

{\bf T}_{CD}=\left(\begin{matrix}\cos(\theta)&-\sin(\theta)&0&0\\
\sin(\theta)&\cos(\theta)&0&0\\
0&0&1&0\\
0&0&0&1\end{matrix}\right).

The value and ranges for \theta are exported by the properties theta and thetaRange, and the \theta coordinate index is defined by the constant THETA_IDX. For rendering, the properties shaftLength and shaftRadius control the size of a shaft drawn about the rotation axis, using the faceColor rendering property. A demo is provided by
artisynth.demos.mech.HingeJointDemo.

In addition to the standard constructors described in Section 3.3.3,

  HingeJoint (bodyA, bodyB, originD, zaxis)

creates a hinge joint with a specified origin and z axis direction for frame D (in world coordinates), and frames C and D coincident.

Index type/name description
0 bilateral restricts translation along x
1 bilateral restricts translation along y
2 bilateral restricts translation along z
3 bilateral restricts rotation about x
4 bilateral restricts rotation about y
5 unilataral enforces limits on \theta
0 \theta counter-clockwise rotation of C about the z axis
Table 3.1: Constraints (top) and coordinates (bottom) for the hinge joint.

3.4.2 Slider joint

Figure 3.14: Coordinate frames (left) and demo model (right) for the slider joint.

The SliderJoint (Figure 3.14) is a 1 DOF joint that constrains motion between frames C and D to a simple translation along the z axis of D. It implements six constraints and one coordinate z (Table 3.2), to which the joint transform {\bf T}_{CD} is related by

{\bf T}_{CD}=\left(\begin{matrix}1&0&0&0\\
0&1&0&0\\
0&0&1&z\\
0&0&0&1\end{matrix}\right).

The value and ranges for z are exported by the properties z and zRange, and the z coordinate index is defined by the constant Z_IDX. For rendering, the properties shaftLength and shaftRadius control the size of a shaft drawn about the sliding axis, using the faceColor rendering property. A demo is provided by artisynth.demos.mech.SliderJointDemo.

In addition to the standard constructors described in Section 3.3.3,

  SliderJoint (bodyA, bodyB, originD, zaxis)

creates a slider joint with a specified origin and z axis direction for frame D (in world coordinates), and frames C and D coincident.

Index type/name description
0 bilateral restricts translation along x
1 bilateral restricts translation along y
2 bilateral restricts rotation about x
3 bilateral restricts rotation about y
4 bilateral restricts rotation about z
5 unilataral enforces limits on the z coordinate
0 z translation of C along the z axis
Table 3.2: Constraints (top) and coordinates (bottom) for the slider joint.

3.4.3 Cylindrical joint

Figure 3.15: Coordinate frames (left) and demo model (right) for the cylindrical joint.

The CylindricalJoint (Figure 3.15) is a 2 DOF joint that constrains motion between frames C and D to translation and rotation along and about the z axis of D. It implements six constraints and two coordinates z and \theta (Table 3.3), to which the joint transform {\bf T}_{CD} is related by

{\bf T}_{CD}=\left(\begin{matrix}\cos(\theta)&-\sin(\theta)&0&0\\
\sin(\theta)&\cos(\theta)&0&0\\
0&0&1&z\\
0&0&0&1\end{matrix}\right).

The value and ranges for z and \theta are exported by the properties z, theta, zRange and thetaRange, and the coordinate indices are defined by the constants Z_IDX and THETA_IDX. For rendering, the properties shaftLength and shaftRadius control the size of a shaft drawn about the sliding/rotation axis, using the faceColor rendering property. A demo is provided by artisynth.demos.mech.CylindricalJointDemo.

In addition to the standard constructors described in Section 3.3.3,

  CylindricalJoint (bodyA, bodyB, originD, zaxis)

creates a cylindrical joint with a specified origin and z axis direction for frame D (in world coordinates), and frames C and D coincident.

Index type/name description
0 bilateral restricts translation along x
1 bilateral restricts translation along y
2 bilateral restricts rotation about x
3 bilateral restricts rotation about y
4 unilataral enforces limits on the z coordinate
5 unilataral enforces limits on the \theta coordinate
0 z translation of C along the z axis
1 \theta rotation of C about the z axis
Table 3.3: Constraints (top) and coordinates (bottom) for the cylindrical joint.

3.4.4 Slotted hinge joint

Figure 3.16: Coordinate frames (left) and demo model (right) for the slotted hinge joint.

The SlottedHingeJoint (Figure 3.16) is a 2 DOF joint that constrains motion between frames C and D to translation along the x axis and rotation about the z axis of D. It implements six constraints and two coordinates x and \theta (Table 3.4), to which the joint transform {\bf T}_{CD} is related by

{\bf T}_{CD}=\left(\begin{matrix}\cos(\theta)&-\sin(\theta)&0&x\\
\sin(\theta)&\cos(\theta)&0&0\\
0&0&1&0\\
0&0&0&1\end{matrix}\right). (3.25)

The value and ranges for x and \theta are exported by the properties x, theta, xRange and thetaRange, and the coordinate indices are defined by the constants X_IDX and THETA_IDX. For rendering, the properties shaftLength and shaftRadius control the size of a shaft drawn about the rotation axis, while slotWidth and slotDepth control the width and depth of a slot drawn along the sliding (x) axis; both are drawn using the faceColor rendering property. When rendering the slot, its bounds along the x axis are set to xRange by default. However, this may be too large, particularly if xRange is unbounded. As an alternate, the property slotRange will be used instead if its range (i.e., the upper bound minus the lower bound) exceeds 0. A demo of SlottedHingeJoint is provided by artisynth.demos.mech.SlottedHingeJointDemo.

In addition to the standard constructors described in Section 3.3.3,

  SlottedHingeJoint (bodyA, bodyB, originD, zaxis)

creates a slotted hinge joint with a specified origin and z axis direction for frame D (in world coordinates), and frames C and D coincident.

Index type/name description
0 bilateral restricts translation along y
1 bilateral restricts translation along z
2 bilateral restricts rotation about x
3 bilateral restricts rotation about y
4 unilataral enforces limits on the x coordinate
5 unilataral enforces limits on the \theta coordinate
0 x translation of C along the x axis
1 \theta rotation of C about the z axis
Table 3.4: Constraints (top) and coordinates (bottom) for the slotted hinge joint.

3.4.5 Universal joint

Figure 3.17: Coordinate frames (left) and demo model (right) for the universal joint.
Index type/name description
0 bilateral restricts translation along x
1 bilateral restricts translation along y
2 bilateral restricts translation along z
3 bilateral restricts rotation about the final x axis of C
4 unilataral enforces limits on the roll coordinate
5 unilataral enforces limits on the pitch coordinate
0 \theta (roll) first rotation of C about the z axis of D
1 \phi (pitch) second rotation of C about the rotated y^{\prime} axis
Table 3.5: Constraints (top) and coordinates (bottom) for the universal joint.

The UniversalJoint (Figure 3.17) is a 2 DOF joint that allows C two rotational degrees of freedom with respect to D: a roll rotation \theta about D’s z axis, followed by a pitch rotation \phi about the rotated y^{\prime} axis. It implements six constraints and the two coordinates \theta and \phi (Table 3.5), to which the joint transform {\bf T}_{CD} is related by

{\bf T}_{CD}=\left(\begin{matrix}c_{\theta}c_{\phi}&-s_{\theta}&c_{\theta}s_{%
\phi}&0\\
s_{\theta}c_{\phi}&c_{\theta}&s_{\theta}s_{\phi}&0\\
-s_{\phi}&0&c_{\phi}&0\\
0&0&0&1\end{matrix}\right),

where

c_{\theta}\equiv\cos(\theta),\;s_{\theta}\equiv\sin(\theta),\;c_{\phi}\equiv%
\cos(\phi),\;s_{\phi}\equiv\sin(\phi).

The value and ranges for \theta and \phi are exported by the properties roll, pitch, rollRange and pitchRange, and the coordinate indices are defined by the constants ROLL_IDX and PITCH_IDX. For rendering, the properties shaftLength and shaftRadius control the size of shafts drawn about the roll and pitch axes, while jointRadius specifies the radius of a ball drawn around the origin of D; both are drawn using the faceColor rendering property. A demo is provided by artisynth.demos.mech.UniversalJointDemo.

3.4.6 Skewed universal joint

The SkewedUniversalJoint (Figure 3.18) is a version of the universal joint in which the pitch axis is skewed relative to its nominal direction by an angle \alpha. More precisely, let x^{\prime} and y^{\prime} be the x and