%File: ~/OOP/analysis/integrator/ArcLength.tex
%What: "@(#) ArcLength.tex, revA"

THE IMPLEMENTATION WILL HAVE TO CHANGE FOR DOMAIN-DECOMPOSITION
ANALYSIS .. AS DOES THE CONVERGENCE TEST STUFF .. THIS IS BECAUSE
USING DOT PRODUCTS OF VECTORS OBTAINED STRAIGHT FROM SYSTEM OF
EQUATION .. MAYBE MODIFY LinearSOE TO DO THE DOT PRODUCT .. WILL 
WORK IN DD IF ALL USE ONE SOE .. WHAT PetSC DOES, TALK WITH P. DEMMEL
ABOUT WHAT HE WILL PROVIDE. \\

\noindent {\bf Files}   \\
\indent \#include $<\tilde{ }$/analysis/integrator/ArcLength.h$>$  \\

\noindent {\bf Class Declaration}  \\
\indent class ArcLength: public StaticIntegrator  \\

\noindent {\bf Class Hierarchy} \\
\indent MovableObject \\
\indent\indent Integrator \\
\indent\indent\indent IncrementalIntegrator \\
\indent\indent\indent\indent StaticIntegrator \\
\indent\indent\indent\indent\indent {\bf ArcLength} \\

\noindent {\bf Description} \\ 
\indent ArcLength is a subclass of StaticIntegrator, it is
used to when performing a static analysis on the FE\_Model using an
arc length method. In the arc length method implemented by this class,
the following constraint equation is added to
equation~\ref{staticFormTaylor} of the StaticIntegrator class: 

\begin{equation}
\Delta \U_n^T \Delta \U_n  + \alpha^2 \Delta \lambda_n^2  = \Delta s^2
\end{equation}

where 

\[
\Delta \U_n = \sum_{j=1}^{i} \Delta \U_n^{(j)} = \Delta \U_n^{(i)} +
d\U^{(i)} 
\]

\[
\Delta \lambda_n = \sum_{j=1}^{i} \Delta \lambda_n^{(j)} = \Delta \lambda_n^{(i)} +
d\lambda^{(i)} 
\]

\noindent this equation cannot be added directly into
equation~\ref{staticFormTaylor} to produce a linear system of $N+1$
unknowns. To add this equation we make some assumptions ala Yang
(REF), which in so doing allows us to solve a system of $N$
unknowns using the method of ??(REF).  Rewriting
equation~\ref{staticFormTaylor} as  

\[
\K_n^{(i)} \Delta \U_n^{(i)} = \Delta \lambda_n^{(i)} \P +
\lambda_n^{(i)} \P - \F_R(\U_n^{(i)}) = \Delta \lambda_n^{(i)} \P + \R_n^{(i)}
\]

\noindent The idea of ?? is to separate this into two equations:

\def\Uh{\dot{\bf U}}
\def\Ub{\overline{\bf U}}

\[
\K_n^{(i)} \Delta \Uh_n^{(i)} = \P
\]

\[
\K_n^{(i)} \Delta \Ub_n^{(i)} = \R_n^{(i)}
\]

\noindent where now

\begin{equation}
 \Delta \U_n^{(i)} = \Delta \lambda_n^{(i)} \Delta \Uh_n^{(i)} + \Delta \Ub_n^{(i)}  
\label{splitForm}
\end{equation}

\noindent We now rewrite the constraint equation based on two conditions:

\begin{enumerate}
\item {\bf $i = 1$}: assuming $\R_n^{(1)} = \zero$, i.e. the system is
in equilibrium at the start of the iteration, the following is obtained

\[  \Delta \U_n^{(1)} = \Delta \lambda_n^{(1)} \Delta \Uh_n^{(1)} + \zero \]

\[ \Delta \lambda_n^{(1)} = \begin{array}{c} + \\ - \end{array}
\sqrt{\frac{\Delta s^2}{\Uh^T \Uh + \alpha^2}} \]

\noindent The question now is whether {\bf +} or {\bf -} should be
used. In this class, $d \lambda$ from the previous iteration $(n-1)$
is used, if it was positive {\bf +} is assumed, otherwise {\bf -}. This may
change. There are other ideas: ?(REF) number of negatives on diagonal
of decomposed matrix, ...

\item {\bf $i > 1$}

\[ \left( \Delta \U_n^{(i)} + d\U^{(i)} \right)^T \left( \Delta \U_n^{(i)} +
d\U^{(i)} \right) + \alpha^2 \left( \Delta \lambda_n^{(i)} + d\lambda^{(i)}
\right)^2 = \Delta s^2 \]

\[
\Delta {\U_n^{(i)}}^T\Delta \U_n^{(i)} + 2{d\U^{(i)}}^T\Delta \U_n^{(i)} + {d\U^{(i)}}^T d\U^{(i)}
+ \alpha^2 \Delta {\lambda_n^{(i)}}^2
+ 2 \alpha^2 d\lambda^{(i)} \Delta \lambda_n^{(i)} + \alpha^2 {d\lambda^{(i)}}^2
= \Delta s^2
\]

\noindent assuming the constraint equation was solved at $i-1$,
i.e. ${d\U^{(i)}}^T d\U^{(i)} + \alpha^2 {d\lambda^{(i)}}^2 = \Delta s^2$, this reduces to

\[
\Delta {\U_n^{(i)}}^T\Delta \U_n^{(i)} + 2{d\U^{(i)}}^T\Delta \U_n^{(i)} + 
\alpha^2 \Delta {\lambda_n^{(i)}}^2
+ 2 \alpha^2 d\lambda^{(i)} \Delta \lambda_n^{(i)} 
= 0
\]

\noindent substituting for $\Delta {\U_n^{(i)}} $ using
equation~\ref{splitForm} this can be expressed as:

\[
\Delta \lambda_n^{(i)^2} \left( \Delta \Uh_n^{(i)} \Delta \Uh_n^{(i)} +
\alpha^2 \right) +
2* \Delta \lambda_n^{(i)} \left( \Delta \Uh_n^{(i)} \Delta \Ub_n^{(i)}
+ d\U^{(i)} \Delta \Uh_n^{(i)} 
+ \alpha^2d \lambda^{(i)} \right)
\]
\[
+ \left (\Delta \Ub_n^{(i)} \Delta \Ub_n^{(i)} + d\U^{(i)} \Delta
\Ub_n^{(i)}
\right) =0 
\]

which is a quadratic in $\Delta \lambda_n^{(i)}$, which can be solved for two roots.
The root chosen is the one which will keep a positive angle between
the incremental displacement before and after this step.


\end{enumerate}
\noindent {\bf Class Interface} \\
\indent // Constructors \\
\indent {\em ArcLength(double arc, double $\alpha$ = 1.0);}\\ \\
\indent // Destructor \\
\indent {\em $\tilde{ }$ArcLength();}\\  \\
\indent // Public Methods \\
\indent {\em int newStep(void);} \\
\indent {\em int update(const Vector \&$\Delta U$);} \\
\indent {\em int domainChanged(void); }\\ \\
\indent // Public Methods for Output\\
\indent {\em int sendSelf(int commitTag, Channel \&theChannel);}\\ 
\indent {\em int recvSelf(int commitTag, Channel \&theChannel,
FEM\_ObjectBroker \&theBroker);}\\ 
\indent {\em int Print(OPS_Stream \&s, int flag = 0);}\\

\noindent {\bf Constructors} \\
\indent {\em ArcLength(double dS, double alpha = 1.0);}\\ \\
The integer INTEGRATOR\_TAGS\_ArcLength (defined in
$<$classTags.h$>$) is passed to the StaticIntegrator classes
constructor. The value of $\alpha$ is set to {\em alpha} and 
$\Delta s$ to {\em dS}. \\

\noindent {\bf Destructor} \\
\indent {\em $\tilde{ }$ArcLength();}\\ 
Invokes the destructor on the Vector objects created in {\em
domainChanged()}. \\

\noindent {\bf Public Methods}\\

{\em int newStep(void);} \\
{\em newStep()} performs the first iteration, that is it solves for 
$\lambda_n^{(1)}$ and $\Delta \U_n^{(1)}$ and updates the model with
$\Delta \U_n^{(1)}$ and increments the load factor by
$\lambda_n^{(1)}$. To do this it must set the rhs of the LinearSOE to
$\P$, invoke {\em formTangent()} on itself and solve the LinearSOE to
get $\Delta \Uh_n^{(1)}$. \\

{\em int update(const Vector \&$\Delta U$);} \\
Note the argument $\Delta U$ should be equal to $\Delta \Ub_n^{(i)}$.
The object then determines $\Delta \Uh_n^{(i)}$ by setting the rhs of
the linear system of equations to be $\P$ and then solving the
linearSOE. It then solves for
$\Delta \lambda_n^{(i)}$ and $\Delta \U_n^{(i)}$ and updates the model with
$\Delta \U_n^{(i)}$ and increments the load factor by $\Delta
\lambda_n^{(i)}$. Sets the vector $x$ in the LinearSOE object to be
equal to $\Delta \U_n^{(i)}$ before returning (this is for the
convergence test stuff. \\


\indent {\em int domainChanged(void); }\\ 
The object creates the Vector objects it needs. Vectors are created to
stor $\P$, $\Delta \Ub_n^{(i)}$, $\Delta \Uh_n^{(i)}$, $\Delta
\Ub_n^{(i)}$, $dU^{(i)}$. To form $\P$, the current load factor is
obtained from the model, it is incremented by $1.0$, {\em
formUnbalance()} is invoked on the object, and the $b$ vector is
obtained from the linearSOE. This is $\P$, the load factor on the
model is then decremented by $1.0$. \\

{\em int sendSelf(int commitTag, Channel \&theChannel); } \\ 
Places the values of $\Delta s$ and $\alpha$ in a
vector of size $2$ and invokes {\em sendVector()} on {\em theChannel}.
Returns $0$ if successful, a warning message is printed and
a $-1$ is returned if {\em theChannel} fails to send the Vector. \\

{\em int recvSelf(int commitTag, Channel \&theChannel, 
FEM\_ObjectBroker \&theBroker); } \\ 
Receives in a Vector of size 2 the values of $\Delta s$ and $\alpha$.
Returns $0$ if successful, a warning message is printed, $\delta
\lambda$ is set to $0$, and a $-1$ is returned if {\em theChannel}
fails to receive the Vector.\\

{\em int Print(OPS_Stream \&s, int flag = 0);}\\
The object sends to $s$ its type, the current value of $\lambda$, and
$\delta \lambda$. 