 \section{Mappings} \label{sec:mappings}
Fran\c{c}ois Faure


\newcommand{\mass}{\ensuremath{M}}             % matrice de masse
\newcommand{\vol}{\ensuremath{\mathcal V}} % volume d'intersection
\newcommand{\press}{\ensuremath{\rho}}

% In section~\ref{sec:pixelcontact}, we have shown how to efficiently compute repulsion forces between arbitrary polyhedra.
% In this section, we present a general framework to map the contact forces to the degrees of freedom of arbitrary physical models.
% We then demonstrate it on a variety of models.

\subsection{Motivation} \label{sec:geometryLayers}
Different geometrical models can be used to model objects in contact.
We organize them in a hierarchy of layers. An example is shown in figure~\ref{fig:hierarchy}, where a rigid object hits a shape embedded in deformable cells.

\begin{figure}
 \centering
 \includegraphics[width=\linewidth]{mappings.png}
 \caption{Mappings from the DOFs to a contact point. Top: two simulated objects in contact (red point). Bottom: hierarchy of geometrical layers. Positions and velocities are propagated top-down. The contact force $f_c$ is accumulated in the contact layers. Forces are then propagated bottom-up.
%Each object has its own state vectors and mappings.
}
 \label{fig:hierarchy}
\end{figure}

% bases de la dynamique
The state of a simulated system can be described by the values and time derivatives of its independent degrees of freedom (DOF) gathered in two vectors $x_0$ and $v_0$.
The dynamics equation (Newton's law) relates the second time derivative $a_0$ of the DOF to the forces $f_0$ acting on them: $f_0 = \mass a_0$, where \mass~ is a matrix modeling the mass of the system.
%We call $f_0$ the \textit{net force}.

% attacher de la g�om�trie
A geometrical model can be attached to the DOF for visualization or contact computation. 
Its DOF positions, velocities and associated forces are stored in vectors $x_1$, $v_1$ and $f_1$, respectively.
They are not independent variables, since the positions and velocities are bound to the independent DOF.
We say that the child geometrical model $1$ is mapped from the parent model $0$,
 using a kinematic operator which we call \textit{mapping}:
\begin{eqnarray*} %\label{eq:mapV}
x_1 &=&\JNL_1(x_0)\\ 
v_1 &=& J_1 v_0
\end{eqnarray*}
Typical mappings include polygonal shapes attached to rigid bodies using local coordinates, or embedded in deformable cells using barycentric coordinates, as well as skin surrounding articulated bodies using vertex blending techniques.
Matrix $J_1 = \frac{\partial x_1}{\partial x_0}$ encodes the linear relation between the DOF velocities and the shape velocities. Due to linearity, the same relation holds on elementary displacements $dx$.
It also holds on accelerations, with an additional offset due to velocities when the position mapping \JNL is nonlinear.
In most cases, operators \JNL~ and \J~ are the same, but in the case of rigid bodies, \JNL~ is nonlinear with respect to $x_0$ and it can not be written as a matrix.
For surfaces embedded in deformable cells, matrix \J~contains the barycentric coordinates. 
For surfaces attached to rigid bodies, each row of the matrix encodes the usual relation $v = \dot o + \omega \times (x-o)$ for each vertex. 
Similarly, skins around articulated bodies involve, at each vertex, the weighted  contributions of the rigid bodies. 


% g�om�trie suppl�mentaire d�e aux contacts
When shapes collide, additional geometry can be necessary to model the contact.
For instance, when an edge intersects another one, a contact force is applied to the intersection points.
These points are defined by their barycentric coordinates with respect to their edge vertices. 
Other relations can be used, depending on the kind of geometrical primitives in contact.
This additional geometry requires another geometrical layer connected to the shape by a mapping, as illustrated in figure~\ref{fig:hierarchy}.
% This layer is also connected to the shape using mappings:
% \begin{eqnarray*} %\label{eq:mapV}
% x_2 &=&\JNL_2(x_1)\\ 
% v_2 &=&J_2 v_1 
% \end{eqnarray*}
We extend this approach to tree-like hierarchies of geometries, with the independent DOFs at the root. 
For instance, the independent DOF may have two children, one for collision using a coarse mesh, and the other for rendering using a finer mesh. The synchronization between these siblings is automatically guaranteed by their attachment to their common ancestor, the DOF layer.

Positions and velocities are propagated top-down in the hierarchy. Conversely,
in order to take all the forces into account in the dynamics equation, the forces are propagated bottom-up, up to the independent DOFs, where Newton's law $f=ma$ is applied. 
Given forces $f_c$ applied to a child model, the mapping computes and accumulates the equivalent forces $f_p$ applied to its parent. 
Since equivalent forces must have the same power, the following relation holds:
$$
v_{p}^T f_p = v_c^T f_c
$$
The kinematic relation $v_{c} = J v_{p}$ allows us to rewrite the previous equation as
$$
v_{p}^T f_{p} = v_{p}^T J^T f_c
$$
Since this relation holds for any velocity $v_p$, the principle of virtual work allows us to simplify the previous equation to obtain:
\begin{equation} \label{eq:mapF}
f_{p} = J^T f_c
\end{equation}

\subsection{Mapping functions}
As seen in section~\ref{sec:geometryLayers}, the mappings propagate positions, velocities, displacements and accelerations top-down, and they propagate forces bottom-up.
The top-down propagation methods are:
\begin{itemize}
 \item \texttt{apply (const MechanicalParams*, MultiVecCoordId outPos, ConstMultiVecCoordId inPos )} for positions,
 \item \texttt{applyJ(const MechanicalParams*, MultiVecDerivId outVel, ConstMultiVecDerivId inVel )} for velocities and small displacements,
\item \texttt{computeAccFromMapping(const MechanicalParams*, MultiVecDerivId outAcc, ConstMultiVecDerivId inVel, ConstMultiVecDerivId inAcc )} for accelerations, taking into account velocity-dependent accelerations in nonlinear mappings.
\end{itemize}
The bottom-up propagation methods are:
\begin{itemize}
 \item \texttt{applyJT(const MechanicalParams*, MultiVecDerivId inForce, ConstMultiVecDerivId outForce )} for child forces or changes of child forces,
\item \texttt{applyDJT(const MechanicalParams*, MultiVecDerivId parentForce, ConstMultiVecDerivId  childForce )} for changes of parent force due to a change of mapping with constant child force,
\item \texttt{applyJT(const ConstraintParams*, MultiMatrixDerivId inConst, ConstMultiMatrixDerivId outConst )} for constraint Jacobians,
\end{itemize}


The name of the methods used to propagate velocities or small displacements top-down contain $J$, which denotes the kinematic matrix, while the names of the methods used to propagate forces or constraint Jacobians bottom-up contain $JT$, which denotes the transpose of the same.
Method \texttt{applyJT(const MechanicalParams*, MultiVecDerivId inForce, ConstMultiVecDerivId outForce )} is used to accumulate forces from a child model to its parent. It performs a cumulative write (+=) since a model may have several children:
\begin{equation}
 \label{eq.applyjt}
 f_{p} += J^T f_c
\end{equation}

Some differential equation solvers need compute the change of force df, given a change of position dx.
The displacement dx, considered small, is propagated top-down using the linear operator \texttt{applyJ(const MechanicalParams*, MultiVecDerivId outVel, ConstMultiVecDerivId inVel )}, then the force changes are accumulated bottom-up.
Differentiating eq.\ref{eq.applyjt}, we get:
\begin{equation}
 \label{eq.applydjt}
 \delta f_{p} += J^T \delta f_c + \delta J^T f_c
\end{equation}
Once the change of child force $\delta f_c$ is computed (see the section on force fields), method \texttt{applyJT(const MechanicalParams*, MultiVecDerivId inForce, ConstMultiVecDerivId outForce )} is used to accumulate it in the parent, corresponding to the first term in the right of eq.\ref{eq.applydjt}.
Method \texttt{applyDJT(const MechanicalParams*, MultiVecDerivId parentForce, ConstMultiVecDerivId  childForce )} is used to accumulate the second term, which is due to the change of matrix J due to a displacement.
It is null in linear mappings, such as BarycentricMapping.
This method queries the last displacement propagated and the child force using the MechanicalParams.


Constraints enforced using Lagrange multipliers are represented using linear equations.
If a linear constraint on the child DOFS in expressed as $L_c v_c = a$, where $L_c$ is the Jacobian of the constraint, then the equivalent constraint at the parent level is: $L_p v_p = a$, where $L_p = J^T L_c$.
Method applyJT(const ConstraintParams*, MultiMatrixDerivId inConst, ConstMultiMatrixDerivId outConst ) is used to compute $L_p$.
Since the Jacobians are generally sparse, they are encoded in sparse matrices instead of the dense vectors used for forces.



%This simple relation allows us to map forces bottom-up through the hierarchy. 
%The corresponding acceleration can then be computed in the DOF layer.

% We emphasize the simplicity and the generality of equation~\ref{eq:mapF}. 
% It makes no assumption about the constitutive laws of the simulated objects.
% It allows us to connect arbitrary bodies using any geometry, provided that a kinematic operator linking the velocity of the geometry to the body DOFs is available.
% %It is not limited to polygonal geometries.

% It encompasses and generalizes \cite{Sifakis07}'s hard bindings which were introduced for Finite Elements and empirically extended to rigid bodies.
% It enables us to efficiently and straightforwardly connect or make collide a wide variety of bodies.

% Let $x_n$, $y_n$, $f_n$ be, respectively, the positions, velocities and net forces in the geometrical model $n$, and $x_{n+1}$, $y_{n+1}$, $f_{n+1}$ the positions, velocities and net forces in model $n+1$ mapped from model $n$.
% The mappings have the following propagation methods:
% \begin{itemize}
%  \item apply: propa
% \end{itemize}


