

% \begin{abstract}
% This is the documentation of the SOFA library. Chapter~\ref{chapter:pba} gives theoretical background on physically-based animation. Chapter~\ref{chapter:as} shows how to implement simulations using \sofa. Chapter~\ref{chapter:es} describes how to integrate new components in \sofa.
% \end{abstract}

%\tableofcontents





\section{Commented example} \label{sec:commentedExample}
Figure~\ref{fig:mixedPendulum} shows a simple scene composed of two different objects, one rigid body and one particle system, and linked by a spring.
\begin{figure}
 \centering
 \includegraphics[width=0.9\linewidth]{mixedPendulum.png}
 \caption{A pendulum composed of a rigid body (reference frame and yellow point) attached to an elastic string (green) fixed at one end (pink point). 
 The corresponding scene graph is displayed on the left.}
 \label{fig:mixedPendulum}
\end{figure}
This scene is modeled and simulated in C++ as shown in section~\ref{cpp:hybrid}. 
The corresponding scene graph is shown in figure~\ref{fig:mixedPendulum-graph}. 
Note that the graph in the left of figure~\ref{fig:mixedPendulum} only displays a hierarchical view, while the whole graph includes additional pointers displayed as dashed arrows in figure~\ref{fig:mixedPendulum-graph}.
\begin{figure}
 \centering
 \includegraphics[width=\linewidth]{mixedPendulum-graph}
 \caption{The scene graph of the mixed pendulum. 
 The nodes are displayed as grey hexagons, while the components are displayed as rectangles with colors associated with their types or roles. 
 The bold plain arrows denote node hierarchy, while the thin plain arrows point to the components attached to the nodes, and the dotted arrows denote pointers between components.}
 \label{fig:mixedPendulum-graph}
\end{figure}

The scene is modeled as a tree structure with four nodes:
\begin{itemize}
 \item \texttt{root}
 \item \texttt{deformableBody} corresponds to the elastic string
 \item \texttt{rigidBody} corresponds to the rigid object
 \item \texttt{rigidParticles} corresponds to a set of particles (only one in this case) attached to the rigid body
\end{itemize}
Each node can have children nodes and \textit{components}. 
Each component implements a reduced set of functionalities.


One of the most important type of component is the \texttt{MechanicalObject}, which contains a list of \textit{degrees of freedom} (DOF), i.e. coordinates, velocities, and associated auxiliary vectors such as forces and accelerations.
All the coordinates in a \texttt{MechanicalObject} have the same type, e.g. 3D vectors for particles, or (translation, rotation) pairs for rigid bodies. 
\texttt{MechanicalObject}, like many other \sofa{} classes, is a generic (C++ template) class instantiated on the types of DOF it stores.
The particle DOFs are drawn as white points, whereas the rigid body DOFs are drawn as red, green, blue reference frame axes.
There can be at most one \texttt{MechanicalObject} attached to a given node. 
This guarantees that all the components attached to the same node process the same types of DOF. 
Consequently, the particles and the rigid body necessarily belong to different nodes.

In this example, the masses are stored in \texttt{UniformMass} components.
The types of their values are related to the types of their associated DOF.
\texttt{UniformMass} is derived from the abstract \texttt{Mass} class, and stores only one value, for the case where all the associated objects have the same mass. 
If necessary, it can be replaced by a \texttt{DiagonalMass} instantiated on the same DOF types, for the case where the associated objects have different masses. 
This is an important feature of \sofa: each component can be replaced by another one deriving from the same abstract class and instantiated on the same DOF types. 
This results in a high flexibility.

The \texttt{FixedConstraint} component attached a particle to a fixed point in world space, drawn in pink. 
The constraints act as filters which cancel the forces and displacements applied to their associated particle(s). 
They do not model more complex constraints such as maintaining three points aligned.

The \texttt{StiffSpringForceField} stores a list of springs, each of them modeled by a pair of indices, as well as the standard physical parameters, stiffness, damping and rest length.

The rigid body is connected to the deformable string by a spring.
Since this spring is shared by the two bodies, it is modeled in the \texttt{StiffSpringForceField} attached to a common ancestor, the graph root in this example.
Our springs can only connect particles. 
We thus need to attach a particle to the rigid body. 
Since the particle DOFs types are different from the rigid body DOF types, they have to be stored in another \texttt{MechanicalObject}, called \texttt{rigidParticleDOF} in this example, and attached to a different node.
However, \texttt{rigidParticleDOF} is not a set of independent DOF, since they are fixed in the reference frame of the rigid body. We thus attach it to a child node of the rigid body, and connect it to \texttt{rigidDOF} using a \texttt{RigidMapping}. 
This component stores the coordinates of the particle in the reference frame of the rigid body. 
Its task is to propagate the position, velocity and displacement of the rigid body down to the yellow particle, and conversely, to propagate the forces applied to the particle up to the rigid body.

Mappings are one of the major features of \sofa. 
They allow us to use different geometric models for a given body, e.g. a coarse tetrahedral mesh for viscoelastic internal forces, a set of spheres for collision detection and modeling, and a fine triangular mesh for rendering.

The gravity applied to the scene is modeled in the \texttt{Gravity} component near the root. 
It applies to all the scene, unless locally overloaded by another gravity component inside a branch of the tree.

The abstract component classes are defined in namespace \texttt{core::componentmodel}.

So far, we have discussed the physical model of the scene.
To animate it, we need to solve an \textit{Ordinary Differential Equation} (ODE) in time.
There are plenty of ODE solvers, and \sofa{} allows the design and the re-use of a wide variety of them.
Here we use a simple explicit Euler method, modeled using an \texttt{EulerSolver} component.
It triggers computations such as force accumulation, acceleration computation and linear operations on state vectors.
More sophisticated solvers are available in \sofa, and can be used by simply replacing the  \texttt{EulerSolver} component by another one, e.g \texttt{RungeKutta4} or \texttt{CGImplicit}.

Other capabilities of \sofa, such as collision detection and response, will be discussed in subsequent sections.

\section{Multi-model objects} \label{sec:multimodel}
An important feature of Sofa is the possibility of using different models of a single physical object. Figure~\ref{fig:liver} shows a scene graph representing a liver, and three different images of it.
The liver exhibits three different geometries for mechanics, rendering and collision.
\begin{figure}
 \centering
 \includegraphics[width=0.95\linewidth]{liver_graph}\\
 \includegraphics[width=0.3\linewidth]{liver_smooth_visu}
 \includegraphics[width=0.3\linewidth]{liver_behavior}
 \includegraphics[width=0.3\linewidth]{liver_collision}
 \caption{A liver. Top: scene graph. Bottom: visual model, mechanical model, collision model, respectively.}
 \label{fig:liver}
\end{figure}
The corresponding xml code is given in section~\ref{xml:liver}.
\label{bla:liver}

On top of the scene, collision-related components allow a user to interact with the collision models using rays cast from the mouse pointer and hitting collision models. 
%Collision is discussed in section~\ref{sec:collision}.

The liver is modeled using three nodes, in two levels. 
The parent level contains the mechanical DOFs (particle positions and velocities) in a \texttt{MechanicalObject} component. 
These DOFs are the mechanically independent degrees of freedom of the object, in Lagrange's formalism. 
The node also contains components related to the dynamics of the particles, such as mass and internal forces. 
We call it the \textit{behavior model}.

The two other nodes are in the lower level because during the simulation, their coordinates are totally defined by the coordinates of their parent node. 
Thus, they do not belong to the set of mechanically independent DOFs. 
\emph{Mappings} are used to compute their positions and velocities based on their parent's, using the pointers represented as dashed arrows. 
Mappings are not symmetric. 
The motion of the parent DOFs is mapped to the children DOFs, whereas the motion of the children DOFs is not mapped to their parent. 
This ensures consistency.

The \texttt{VisualModel} has vertices which are used for rendering, along with other rendering data such as a list of polygons, normals, etc. 
The mapping is one-way and the mapped DOFs have no mechanical influence.

The \texttt{SphereModel} class derives from \texttt{MechanicalObject}, with an additional radius value. 
It also derives from \texttt{CollisionModel}, which allows it to be processed by the collision detection and modeling pipeline. 
When contact or mouse interaction forces are applied to the spheres, the forces are propagated bottom-up to their parent DOFs by the mapping (see section~\ref{sec:mappings}). 
This allows the contact forces to be taken into account in the dynamics equations. 
The mapping is thus two-ways and derives from \texttt{MechanicalMapping} instead of \texttt{Mapping}. 
This is why it has a different color in the image of the scene graph.
Again, the mechanical mappings are not symmetric: the forces are propagated from the children to the parents, not the other way round.

Mappings only propagate positions top-down, whereas MechanicalMappings additionally propagate velocities top-down and forces bottom-up.

Mapped models can be designed independently of their parent models, provided that the adequate (mechanical) mapping is available. 
This results in a high flexibility. 
For example, collision spheres can be replaced by collision triangles without changing anything in the behavior model or in the visual model. 
Similarly, other visual models can be used without modifying the behavior and collision models, and different behavior models can be used with the same collision and visual models, as illustrated in figure~\ref{fig:behaviormodels}.

\begin{figure}
 \centering
 \includegraphics[width=0.4\linewidth]{demoLiverFall1.png}
 \includegraphics[width=0.4\linewidth]{demoLiverFall2.png}
 % demo.: 1179666x1179666 pixel, 0dpi, infxinf cm, bb=
 \caption{Left: four behavior models (from left to right: deformable grid, springs, rigid, tetrahedral FEM) combined with the same collision model (right).}
 \label{fig:behaviormodels}
\end{figure}


\section{Recursive data processing}
A typical simulation program, controlled by an application such as the Graphics User Interface (GUI), looks like the one given in figure~\ref{pc:animationloop}.
\begin{figure}
\begin{code_cpp}
init();
repeat {
    animate();
    draw();
}
\end{code_cpp}
\caption{Pseudocode for a standard simulation program.}
\label{pc:animationloop}
\end{figure}
In \sofa, each of the simulation methods is implemented as a recursive graph traversal, \texttt{InitVisitor}, \texttt{AnimateVisitor} and \texttt{VisualDrawVisitor}, respectively. 
Visitors are explained in the next section.

\subsection{Visitors}
The data structure is processed using objects called \emph{visitors}.
They recursively traverse the tree structure and call appropriate virtual methods to a subset of components during the \textit{Top-Down Traversal} (TDT), using virtual method \texttt{Visitor::processNodeTopDown}, then during the \textit{Bottom-Up Traversal} (BUT), using virtual method \texttt{Visitor::processNodeBottomUp}.

For example, the \texttt{VisualDrawVisitor} draws the \texttt{VisualModel} components during the TDT, and does nothing during the BUT.
The \texttt{MechanicalComputeForceVisitor} accumulates the forces in the appropriate DOF vectors during the TDT, then propagates the forces to the parent DOFs using the mechanical mappings during the BUT.

When processed by a visitor $a$, a component can fire another visitor $b$ through its associated sub-tree. 
Visitor $a$ can continue once visitor $b$ is finished.
During the TDT, each traversed component decides whether the calling visitor continues, or prunes the sub-tree associated with the component, or terminates.

The components directly access their sibling components only, except for the mappings.
A component traversed by a visitor can indirectly access the data in its associated sub-tree in read-write mode using visitors, whereas data in its parent graph is read-only and only partially accessible using method \texttt{getContext}.
Sibling nodes of the same type can be traversed by visitors in arbitrary order.

The visitors belong to namespace \texttt{sofa::simulation}.

\subsection{ODE Solvers}
When an \texttt{AnimateVisitor} traverses a node with an \texttt{OdeSolver} component,
the solver takes the control of its associated subtree and prunes the \texttt{AnimateVisitor}. 
The solver triggers visitors in its associated subtree to perform the standard mechanical computations and integrate time.

The simplest solver is the explicit Euler method, implemented in \texttt{EulerSolver}. 
The algorithm is shown as pseudocode in figure~\ref{pc:expliciteuler}.
\begin{figure}
\begin{code_cpp}
f = 0
accumulateForces(f,x,v);
a = f/M;
a = filter(a);
x += v * dt;
v += a * dt;
\end{code_cpp}
\caption{Pseudocode for explicit Euler integration.}
\label{pc:expliciteuler}
\end{figure}
Net force is computed in the first line.
In the second line, the acceleration is deduced by dividing the force by the mass.
Then the accelerations of the fixed points are canceled.
Finally, position and velocity are updated.

This algorithm can not be directly implemented in \sofa{} because there are no state vectors x,v,f,a which gather the state values of all the objects in the scene.
The solver processes an arbitrary number of objects, of possibly different types, such as particles and rigid bodies. Each physical object carries its state values and auxiliary vectors in its own \texttt{MechanicalObject} component, which is not directly accessible to the solver.

The solvers represent state vectors as \texttt{MultiVector} objects using symbolic identificators implemented in class \texttt{VecId}.
There are four staticly predefined identificators: \texttt{VecId::position()}, \texttt{VecId::velocity()}, \texttt{VecId::force()} and \texttt{VecId::dx()}.
A \texttt{Multivector} declared by a solver with a given VecId implicitly refers to all the state vectors in the different \texttt{MechanicalObject} components with the same \texttt{VecId} in the solver's subtree.

Vector operations can be remotely triggered by a solver using a visitor of a given type, which defines the operator, and given \texttt{VecId}s, which define the operands.
During the subtree traversal, the operator is applied to the given vectors of the traversed \texttt{MechanicalObject} components.

For example, let us comment the visitors performed by the \texttt{EulerSolver} shown in figure~\ref{fig:mixedPendulum}. Its implementation is in method \texttt{component::odesolver::EulerSolver::solve(double)}.
First, multivectors are declared.

Then method \texttt{core::componentmodel::behavior::OdeSolver::computeForce(VecId)} is called. 
It first fires a \texttt{MechanicalResetForceVisitor} to reset the force vectors of all the \texttt{MechanicalObject} components. 
It then fires a  \texttt{MechanicalComputeForceVisitor}. 
During the TDT, each component derived from \texttt{core::componentmodel::behavior::BaseForceField} computes and accumulates its force in its sibling \texttt{MechanicalObject}. 
In the example shown in figure~\ref{fig:mixedPendulum}, \texttt{F13} adds its contribution to \texttt{Dof1} and \texttt{Dof3}, then \texttt{F1} and \texttt{M1} add their contributions to \texttt{Dof1}, then \texttt{M2} to \texttt{Dof2}. 
Then during the BUT, the mechanical mappings sum up the forces of their child DOF to their parent DOF, \textit{i.e.}, the force in \texttt{Dof3} to \texttt{Dof2} through \texttt{M23} in the same example.
Note that branches \texttt{deformableBody} and \texttt{rigidBody} can be processed in parallel.
At the end, the force vector in \texttt{Dof1} contains the net force applied to the particles, and the force vector in \texttt{Dof2} contains the net (six-dimensional) force applied to the rigid body.

Then method \texttt{OdeSolver::accFromF(VecId,VecId)} fires a \texttt{MechanicalAccFromFVisitor}. 
Each component derived from \texttt{core::componentmodel::behavior::BaseMass} computes the accelerations corresponding to the forces in its sibling \texttt{MechanicalObject}.

Then method \texttt{OdeSolver::projectResponse} fires a \texttt{MechanicalApplyConstraintsVisitor}. 
All the \texttt{core::componentmodel::behavior::BaseConstraint} components (component \texttt{C} in the example) filter the acceleration vector to maintain some points fixed.

Once the acceleration is computed, multivector methods are used to update the positions and velocities. 
Here again, visitors are used to perform the desired operation in each traversed \texttt{MechanicalObject}.

MultiVector operations are pruned at the first level for efficiency, because the solvers deal with the mechanically independent state variables rather than the mapped variables.
Moreover, the mapped coordinates can not be assumed to vary linearly along with their parent variables.
Applying a \\ \texttt{MechanicalPropagatePositionAndVelocityVisitor} is thus necessary to update the mapped DOFs based on the mechanically independent DOFs.
This visitor is automatically performed after time integration, as one can see in the code of method \texttt{MechanicalIntegrationVisitor::fwdOdeSolver}.
It is also used by some solvers when auxiliary states are needed, as discussed in section~\ref{sec:statevectors}, in order to update the mapped DOFs.

\todo{Call tree of AnimateVisitor, restricted to sofa::core and sofa::simulation}

\todo{Discuss independent and shared solvers}

\section{State vectors} \label{sec:statevectors}
The state vectors contain the coordinates, velocites, and other DOF-related values such as force and acceleration.
They are stored in \texttt{MechanicalObject} components.
This template class can be instantiated on a variety of types to model particles, rigid bodies or other types of bodies.
The template parameter is a \texttt{DataTypes} class which describes data and data containers, such as the the type of coordinates and coordinate derivatives used.
These two types are the same in the case of particles, but they are different in the case of rigid bodies.

Each \texttt{MechanicalObject} can represent a set of physical objects of the same type, such as particles.
The coordinate state vectors are defined by the \texttt{VecCoord} type, while the derivatives (velocity, acceleration, force, small displacement) are defined by the \texttt{VecDeriv} type.
Each \texttt{MechanicalObject} stores two arrays of state vectors, one for coordinates and the other for derivatives, as illustrated in figure~\ref{fig:mechanicalobject}.
\begin{figure}
 \centering
 \includegraphics[width=0.45\linewidth]{MechanicalObject1}
\includegraphics[width=0.45\linewidth]{MechanicalObject2}
  \caption{A \texttt{MechanicalObject} and a component addressing it. Left: using the default state vectors. Right: using auxiliary state vectors.}
 \label{fig:mechanicalobject}
\end{figure}

Auxiliary vectors are necessary for complex solvers, such as \texttt{RungeKutta2Solver}. 
This solver first performs a half-length Euler step, then evaluates the derivative of this new state (called the \emph{midpoint}), and finally uses this derivative to update the initial state over a whole time step.

To compute the forces at the midpoint while keeping the initial state for further use, we use the auxiliary vectors \texttt{newX} and \texttt{newV}.
However, components such as forces and constraints use state vectors, and we have to make sure that they use the right ones.
To ensure consistency and make the use of auxiliary states transparent, the other components get access to the state vectors using methods \texttt{MechanicalObject::getX()}, \texttt{getV()}, \texttt{getF()} and  \texttt{getDx()}.
These methods return pointers to the appropriate vectors, as illustrated in figure~\ref{fig:mechanicalobject}.

Internal  \texttt{MechanicalObject} switches are performed by methods \texttt{MechanicalObject::setX()}, \texttt{setV()}, \texttt{setF()} and  \texttt{setDx()}.
These methods are applied by the visitors which take multivectors as parameters, before they use other components. 
See, for example, method\\ \texttt{MechanicalPropagatePositionAndVelocityVisitor::fwdMechanicalState}.

Note that some constraint-based animation methods require large state vectors and matrices encompassing all the mechanical objects of the scene.
Such methods are currently under development in \sofa, and they are not yet documented.
They use visitors to count the total number of scalar DOFs and to gather them in large state vectors, as well as to build mechanical matrices  such as mass, stiffness, damping and compliance etc.

\subsection{Mechanical groups}
During the simulation, each solver prunes the \texttt{AnimateVisitor} which traverses it and manages its associated subtree by itself using other visitors.
The objects animated by a given solver are called a \emph{mechanical group}.
Each mechanical group corresponds to a subtree in the scene graph.
In the example discussed in section~\ref{sec:commentedExample}, there is one mechanical group because a single solver located near the root manages the whole scene.
However, using separate solvers for different objects can sometimes increase efficiency.
In the example shown in figure~\ref{fig:twoSolvers}, the same deformable body is animated using a \texttt{RungeKutta2Solver} while the rigid body is animated using an \texttt{EulerSolver}.
\begin{figure}
 \centering
 \includegraphics[width=0.95\linewidth]{twoSolvers}
  \caption{A scene graph with objects animated using different ODE solvers.}
 \label{fig:twoSolvers}
\end{figure}


A mechanical group can include interaction forces between elements of the group, and such interaction forces are handled by the solver as expected.
Interaction forces can also occur between objects which do not belong to the same group.
In this case, the interaction force is located at a higher hierarchical level than the objects it applies to, as shown in figure~\ref{fig:twoSolvers}.
It can not be traversed by visitors fired by the solvers.
Its evaluation is performed by the \texttt{AnimateVisitor}, and accumulated as external forces in the associated \texttt{MechanicalObject} components.
Consequently, it acts as a constant force during each whole animation step.
In a \texttt{RungeKutta2Solver}, during the force computation at midpoint, its value is the same as at the starting point.
In a \texttt{CGImplicitSolver}, its stiffness is not taken into account, which may introduce instabilities if its actual stiffness is high.

The default collision manager of Sofa circumvents this problem by dynamically gathering the objects in contact in a common mechanical group.

% \section{Mostly used components and methods}
% \subsection{ForceFields}
% \subsection{Constraints}
% \subsection{Masses}
% \subsection{Mappings}
% \subsection{ContextObjects}
%
% \section{Collision detection} \label{sec:collision}
%
%
% \section{Limitations}

% \pagebreak
% \appendix
\section{Code of the examples}
\subsection{The hybrid pendulum}\label{cpp:hybrid}
This is the code of the example commented in section~\ref{sec:commentedExample} :
\includecode{C++}{../../applications/tutorials/mixedPendulum/Main.cpp}

\subsection{A liver}\label{xml:liver}
The XML code of the liver discussed in section~\ref{sec:multimodel} page~\pageref{bla:liver} is in \texttt{../examples/Demos/liver.scn}

% \include{sofaMechanical}
%
% \include{extendingSofa}
%
% \include{physicallyBasedAnimation}
%
% \include{appendices}
