% ----------------------------------------------------------------
% Report Class (This is a LaTeX2e document)  *********************
% ----------------------------------------------------------------
\documentclass[11pt]{report}
\usepackage[english]{babel}
\usepackage{amsmath,amsthm}
\usepackage{amsfonts}
\usepackage{graphicx}
% THEOREMS -------------------------------------------------------
\newtheorem{thm}{Theorem}[chapter]
\newtheorem{cor}[thm]{Corollary}
\newtheorem{lem}[thm]{Lemma}
\newtheorem{prop}[thm]{Proposition}
\theoremstyle{definition}
\newtheorem{defn}[thm]{Definition}
\theoremstyle{remark}
\newtheorem{rem}[thm]{Remark}
% ----------------------------------------------------------------
\begin{document}

\section{Physics Based Animation}

\section{The Particle System}
\par
In the monolithic view, we lump all the particles together and form a large monolithic system.  This state of the system is a phase space vector and is represented as an array of each particles' position in 3 dimensions long with each particles velocity in 3 dimensions.  We then treat this larger vector as the mechanical system and integrate its motion.

\subsection{The Integrator}
\par
We used a symplectic Euler integration scheme for this simulation.

\section{Solving for Constraint forces}

\par
During constraint-based animation, we are computing forces which, if applied, exactly cancel those forces which might cause a particle to go where it is not supposed to go.  For instance, imagine a rigid rod which is free only to pivot about one end.  A particle at the free end of this rigid rod is constrained to be no farther from the pivot than the length of the rod.  Thus, if the rod is rotating about the pivot, the particle at the end of the rod wants to fly off in the direction of its instantaneous velocity.  A force, usually understood as a tension along the rod itself, must be applied to keep that particle moving along the expected trajectory.
\par
It is precisely this force we want to compute during an animation if we wanted to simulate this kind of situation.

\subsection{A Constraint Class}

\par
A major ingredient in creating a general constrained mechanics system is just how one understands a constraint.  Mathematically, a constraint is an equation such as
\begin{equation}\label{stretchConstr}
    \| \vec{x_1} - \vec{x_2} \| - l =0
\end{equation}
which is a stretch constraint in Stam's lingo \cite{nucleus}.  This states that the distance between two points $\vec{x_1}$ and $\vec{x_2}$ should be exactly of length $l$.  The two "points" should refer to the positions of a particle.  One may have other constraints, such as pinning constraints.  A pinning constraint will constrain a particle in a single dimension and would have the form
\begin{equation}\label{xPin}
   \vec{x_1}_x=0
\end{equation}
This would fix the $x$ coordinate of a particle to be exactly zero, and if forces acted on the particle, it would move, but only about the $y$-$z$plane.


\section{Constrained Mechanics}

One can find the exact form of this system in the Pixar notes \cite{PixarNotes}.


\begin{equation}\label{linSystem}
GM^{-1}G^{T} \lambda = -\dot{G} \dot{x}- GMF-k_{sc} g - k_{sd}\dot{g}
\end{equation}

\section{Solving Linear Systems}

\par
To find the constraint forces, we need to solve the linear system \ref{linSystem}.  Thus we are trying to find $\lambda$ such that

%how do we know that GM^{-1}G^{T} is symmetric and positive definite?  Should we include a proof

\begin{equation}\label{linSystemInv}
\lambda = (GM^{-1}G^{T} )^{-1}(-\dot{G} \dot{x}- GMF-k_{sc} g - k_{sd}\dot{g})
\end{equation}

which we recognize as being of the form $x=A^{-1}b$.  To solve this linear system, we take advantage of the fact that the matrix $GM^{-1}G^{T}$ is sparse as well as symmetric and positive definite.

\section{Conjugate Gradients}
\par
A system that is symmetric and positive definite can be solved using the method of conjugate gradients.  An algorithm for conjugate gradients can be found in \cite{MatrixComp}.  Luckily, the matrix which we want to invert in constraint based animations is both symmetric and positive definite.
\par
One of the goals of this project was to do a nice, fast linear solve.  To achieve this, we avoid building the matrices and certainly avoid actually inverting any large dense matrix.  This came back to the constraints class.  This had to harmonize with the method we were going to use to solve the linear system and that method was conjugate gradients.
\par
If we look closely at the algorithm, the matrix $A$ in the linear system $Ax=b$, appears in two places.

\begin{equation}
\alpha_k = r_{k-1}^{T}r_{k-1}/p_{k}^T A p_k
\end{equation}
and
\begin{equation}
    r_k   = r_{k-1}-\alpha_k A p_k
\end{equation}
Each time it appears uninverted and as a part of a matrix vector multiplication.  Recall that in our case $A=GM^{-1}G^{T}$.  Thus, the conjugate gradients method could be programmed as long as we had some methods which could take a bunch of constraints, and a vector and compute the matrix vector computations.  We use associativity so that we do the first matrix multiplication, then the second
\begin{equation}
   ( GM^{-1}G^{T} )x=    (G(M^{-1}(G^{T} x)))
\end{equation}
Thus, we only needed to create methods like \emph{multiplyByG(aVector V)} and \emph{multiplyByGTranspose(aVector V)} which spit out the resulting vector.


\section{Cloth Simulation}
\par
The main goal in this project was to simulate cloth.  In particular, we wanted to simulate inextensible cloth as was done in \cite{inextCloth}.  Hard constraints are a great way to make an inextensible object like cloth.
\par
The setup is quite simple.  The cloth consists of an array of particles, each understood to have a mass, a velocity and forces that act upon it which would then go towards their integration.  The particles are constrained to be a fixed distance from their four neighbors (making a plus sign with the particle in the middle).

\section{Stabilization}
\par
When simulating a constrained system, we find that the state vector in position space will drift off of the constraint manifold.  The prototypical example of this is a simple pendulum which oscillates but also slowly falls in the direction of gravity, eventually slowing down do to the lengthening of the pendulum.  The constraint manifold for this system, a one particle system, is the circle along which we would normally expect the pendulum tip to move.
\par
There are two methods used to stabilize a constrained system, namely Baumgarte and post-step stabilization.  The fast projection algorithm in \cite{inextCloth} is a form of post-step stabilization.  It is an iterative algorithm that steps an unconstrained system back onto the constraint manifold.
\par

\subsection{Error and Efficiency}

\par
We want to look at the errors due to drift.  To do this we need a distance metric.  We follow the method in \cite{KlineNPai} which is to take the error as the largest absolute value over all the constraints in the constraint vector.  We also look at an error metric that uses the constraint that has the largest as the norm of the error.

\section{Minimum Spanning Tree}
\par
To compute a minimum spanning tree, there are two main algorithms Prim and Kruskal (cite the algorithms book).  These algorithms have well known complexities which are functions of the number of nodes and edges in the graph.  The complexity function for Kruskal is $\mathcal{O}(ElogV)$ and Prim is$\mathcal{O}(E+VlogV)$.  We can compute their exact behaviour as a function of square sheet size since there is a function that takes the number of nodes in the sheet and produces the number of edges.
\begin{equation}
E=
\end{equation}
One can compute this function, by looking at the Euler characteristic of a graph.  This relates the vertices, $V$, edges, $E$, and faces,$F$, of a graph as follows 

\begin{equation}
E=V+F-2
\end{equation}
Since we can relate faces to vertices as
\begin{equation}
F=(\sqrt{V}-1)^2+1
\end{equation}
We get the
making the complexity functions
\begin{eqnarray}
C_{Kruskal}=\mathcal{O}() \\
C_{Prim}=\mathcal{O}() 
\end{eqnarray}
At this point, we can compute which one increases faster as sheet size increases.  To do this, we find the ratio of the two functions and then take the limit of this expression as the number of nodes goes to infinity.
\begin{equation}
\lim_{V \rightarrow \infty } \frac{ }{}=\lim_{V \rightarrow \infty } \log(V) =
\end{equation}
This tells us that in this case, where we know that edges are a function of vertices, Kruskal performs worse than prim.  In fact, if we look at the graph of this ratio, we can see that near zero, the rate of increase of the ratio is extremely steep, indicating that even near zero Kruskal is working much harder per node than Prim.

%\begin{figure}
  % Requires \usepackage{graphicx}
%  \includegraphics[width=3]{fig:RatioOfKruskalToPrim.eps}\\
%  \caption{look at this!}\label{kpRatio}
%\end{figure}



Of course, we are interested in only a range of values for the number of nodes which comprise our sheet.  Specifically, we are interested in only that number of nodes that gives an animation which \emph{looks} realistic and no more.  That said, the graph of this ratio when looked at in the range of 0 to 40, which is the largest range we might consider will give hints as to which algorithm to    use.



\bibliographystyle{plain}
\bibliography{clothSimBiblio}


\end{document}
% ----------------------------------------------------------------
