%----------------------------------------------------------
%

\documentclass[10pt,letterpaper,oneside,notitlepage]{article}
%\documentclass{report}%
\usepackage{algorithm}
\usepackage{algpseudocode}
\usepackage{enumitem}
\usepackage{nomencl}
\usepackage{amsmath}
\usepackage{amssymb}
%\usepackage{amsfonts}%
%\usepackage{graphicx}
%----------------------------------------------------------
\makenomenclature
%\theoremstyle{plain}
%\newtheorem{acknowledgement}{Acknowledgement}
%\newtheorem{definition}{Definition}
%\newtheorem{remark}{Remark}
%\numberwithin{equation}{section}
\renewcommand*{\thefootnote}{\fnsymbol{footnote}}
%-----------------------------------------------------------
\begin{document}
\title{Solve Algorithms in OpenFAST}
\author{Bonnie Jonkman}
%\begin{abstract}
%This document is used to describe the algorithms implemented in FAST v8.
%\end{abstract}
\maketitle

%\tableofcontents

\section{Definitions and Nomenclature}


\begin{table}[h]
   \centering
      \begin{tabular}{c|c|c}
      \textbf{Module} & \textbf{Abbreviation} & \textbf{Abbreviation}\\
      \textbf{Name}   & \textbf{in Module}    & \textbf{in this Document}\\
      \hline 
      ElastoDyn          & ED                    & ED        \\
      BeamDyn            & BD                    & BD        \\
      AeroDyn            & AD                    & AD        \\
      ServoDyn           & SrvD                  & SrvD      \\
      SubDyn             & SD                    & SD        \\
      ExtPtfm            & ExtPtfm               & ExtPtfm   \\
      HydroDyn           & HydroDyn              & HD        \\
      MAP++              & MAPp                  & MAP       \\
      FEAMooring         & FEAM                  & FEAM      \\
      MoorDyn            & MD                    & MD        \\
      OrcaFlexInterface  & Orca                  & Orca      \\
      InflowWind         & IfW                   & IfW       \\   
      IceFloe            & IceFloe               & IceF      \\   
      IceDyn             & IceD                  & IceD      \\   
      \end{tabular}
   \caption{Abbreviations for modules in FAST v8}
   \label{tab:Abbrev}
\end{table}


\nomenclature{$u\_ED$}{$ElastoDyn$ inputs}
\nomenclature{$u\_AD$}{$AeroDyn$ inputs}
\printnomenclature

\section{Initializations}


\pagebreak %break here for now so that it doesn't look so strange
\section{Input-Output Relationships}
\subsection {Input-Output Solves (Option 2 Before 1)}
This algorithm documents the procedure for the Input-Output solves in FAST, assuming
all modules are in use. If an individual module is not in use during a particular
simulation, the calls to that module's subroutines are omitted and the module's 
inputs and outputs are neither set nor used.

%\begin{algorithm}[ht]
%\caption{Input-Output Solves (Option 2 Before 1)}
%\label{IOSolves21}
\begin{algorithmic}[1]
\Procedure{CalcOutputs\_And\_SolveForInputs}{\null}

%%%%
% \start SolveOption2
\State
% SolveOption2a_Inp2BD
   \State $\mathit{y\_ED} \gets \Call{ED\_CalcOutput}{\mathit{p\_ED},\mathit{u\_ED},\mathit{x\_ED},\mathit{xd\_ED},\mathit{z\_ED}}$
   \State $\mathit{u\_BD}  \gets \Call{TransferOutputsToInputs}{\mathit{y\_ED,y\_SrvD}}$

\State
% SolveOption2b_Inp2IfW
   \State $\mathit{y\_BD} \gets \Call{BD\_CalcOutput}{\mathit{p\_BD},\mathit{u\_BD},\mathit{x\_BD},\mathit{xd\_BD},\mathit{z\_BD}}$
   \State $\mathit{u\_AD}($no IfW$)  \gets \Call{TransferOutputsToInputs}{\mathit{y\_ED,y\_BD}}$
   \State $\mathit{u\_IfW} \gets \Call{TransferOutputsToInputs}{\mathit{y\_ED} at \mathit{u\_AD} nodes}$

\State
% SolveOption2c_Inp2AD_SrvD
   \State $\mathit{y\_IfW} \gets \Call{IfW\_CalcOutput}{\mathit{u\_IfW} and other \mathit{IfW} data structures}$
   \State $\mathit{u\_AD}($InflowWind only$) \gets \Call{TransferOutputsToInputs}{\mathit{y\_IfW}}$
   \State $\mathit{u\_SrvD} \gets \Call{TransferOutputsToInputs}{\mathit{y\_ED},\mathit{y\_IfW},\mathit{y\_BD}}$

\State
% main
   \State $\mathit{y\_AD} \gets \Call{AD\_CalcOutput}{\mathit{p\_AD},\mathit{u\_AD},\mathit{x\_AD},\mathit{xd\_AD},\mathit{z\_AD}}$
   \State $\mathit{y\_SrvD} \gets \Call{SrvD\_CalcOutput}{}( \!
            \begin{aligned}[t]
                               & \mathit{p\_SrvD},\mathit{u\_SrvD}, \\
                               & \mathit{x\_SrvD},\mathit{xd\_SrvD},\mathit{z\_SrvD}) \\
            \end{aligned}$
   \State $\mathit{u\_ED} \gets \Call{TransferOutputsToInputs}{\mathit{y\_AD},\mathit{y\_SrvD}}$
   \State $\mathit{u\_BD} \gets \Call{TransferOutputsToInputs}{\mathit{y\_AD,y\_SrvD}}$
% \end SolveOption2
%%%%

%%%%
% \begin Transfer_ED_to_HD_SD_BD_Mooring
\State
%   \State $\mathit{u\_ED}($not platform reference point$) \gets \Call{TransferOutputsToInputs}{y\_SrvD,y\_AD}$ %\Comment{sets all but platform reference point inputs}

%   \State $\mathit{u\_BD}   \gets \Call{TransferMeshMotions}{\mathit{y\_ED}}$      % only if not BD_Solve_Option1
   \State $\mathit{u\_HD}   \gets \Call{TransferMeshMotions}{\mathit{y\_ED}}$
   \State $\mathit{u\_SD}   \gets \Call{TransferMeshMotions}{\mathit{y\_ED}}$
   \State $\mathit{u\_ExtPtfm}   \gets \Call{TransferMeshMotions}{\mathit{y\_ED}}$
   \State $\mathit{u\_MAP}  \gets \Call{TransferMeshMotions}{\mathit{y\_ED}}$
   \State $\mathit{u\_FEAM} \gets \Call{TransferMeshMotions}{\mathit{y\_ED}}$
   \State $\mathit{u\_MD}   \gets \Call{TransferMeshMotions}{\mathit{y\_ED}}$
   \State $\mathit{u\_Orca} \gets \Call{TransferMeshMotions}{\mathit{y\_ED}}$
   \State $\mathit{u\_SrvD\%PtfmStC} \gets \Call{TransferMeshMotions}{\mathit{y\_ED}}$\footnote{Only if using ServoDyn Structural control with platform TMD.}
% \end Transfer_ED_to_HD_SD_BD_Mooring
%%%%

\State
\State \Call{SolveOption1}{\null}
\State
   \State $\mathit{u\_IfW}  \gets \Call{TransferOutputsToInputs}{\mathit{u\_AD},\mathit{y\_ED}}$
   \State $\mathit{u\_AD}   \gets \Call{TransferOutputsToInputs}{\mathit{y\_ED}}$
   \State $\mathit{u\_SrvD} \gets \Call{TransferOutputsToInputs}{\mathit{y\_ED},\mathit{y\_AD},\mathit{y\_BD},\mathit{y\_SD}}$

\EndProcedure
\end{algorithmic}
%\end{algorithm}

Note that inputs to \emph{ElastoDyn} before calling CalcOutput() in the first step are not set in CalcOutputs\_And\_SolveForInputs(). 
Instead, the \emph{ElastoDyn} inputs are set depending on where CalcOutputs\_And\_SolveForInputs() is called:
\begin{itemize}[noitemsep] %i don't like the double spaces between bulleted items.
   \item At time 0, the inputs are the initial guess from \emph{ElastoDyn};
   \item On the prediction step, the inputs are extrapolated values from the time history of ElastoDyn inputs;
   \item On the first correction step, the inputs are the values calculated in the prediction step;
   \item On subsequent correction steps, the inputs are the values calculated in the previous correction step.
\end{itemize}


%\pagebreak %break here for now so that it doesn't look so strange
\subsection {Input-Output Solve for \textit{HydroDyn}, \textit{SubDyn}, \textit{OrcaFlexInterface}, \textit{BeamDyn}, \textit{ExtPtfm},  \textit{MAP}, \textit{FEAMooring}, \textit{MoorDyn},  
       \textit{FEAMooring}, \textit{IceFloe}, \textit{IceDyn}, and the Platform Reference Point Mesh in \textit{ElastoDyn}}

This procedure implements Solve Option 1 for the accelerations and loads in
\emph{HydroDyn},\emph{SubDyn},\emph{MAP},\emph{FEAMooring},\emph{OrcaFlexInterface},\emph{MoorDyn}, \emph{BeamDyn}, \emph{ExtPtfm}, \emph{IceFloe}, \emph{IceDyn}, and \emph{ElastoDyn} (at its platform reference point mesh). 
The other input-output relationships for these modules are solved using Solve Option 2.

%\begin{algorithm}[ht]
%\caption{Input-Output Solve for $HydroDyn$, $SubDyn$, $MAP$, $FEAMooring$, and the Platform Reference Point Mesh in $ElastoDyn$}
%\label{IOSolves_PlatformRef}
\begin{algorithmic}[1]

\Procedure{SolveOption1}{\null}
   \State
   \State $\mathit{y\_MAP}     \gets \Call{CalcOutput}{\mathit{p\_MAP},\mathit{u\_MAP},\mathit{x\_MAP},\mathit{xd\_MAP},\mathit{z\_MAP}}$ 
   \State $\mathit{y\_MD}      \gets \Call{CalcOutput}{\mathit{p\_MD},\mathit{u\_MD},\mathit{x\_MD},\mathit{xd\_MD},\mathit{z\_MD}}$
   \State $\mathit{y\_FEAM}    \gets \Call{CalcOutput}{\mathit{p\_FEAM},\mathit{u\_FEAM},\mathit{x\_FEAM},\mathit{xd\_FEAM},\mathit{z\_FEAM}}$
   \State $\mathit{y\_IceF}    \gets \Call{CalcOutput}{\mathit{p\_IceF},\mathit{u\_IceF},\mathit{x\_IceF},\mathit{xd\_IceF},\mathit{z\_IceF}}$
   \State $\mathit{y\_IceD(:)} \gets \Call{CalcOutput}{\mathit{p\_IceD(:)},\mathit{u\_IceD(:)},\mathit{x\_IceD(:)},\mathit{xd\_IceD(:)},\mathit{z\_IceD(:)}}$
   \State $\mathit{y\_SrvD}    \gets \Call{CalcOutput}{\mathit{p\_SrvD},\mathit{u\_SrvD},\mathit{x\_SrvD},\mathit{xd\_SrvD},\mathit{z\_SrvD}}$\footnote{Only if using ServoDyn Structural control with platform TMD.}
   \State
   \State\Comment{Form $u$ vector using loads and accelerations from $\mathit{u\_HD}$, $\mathit{u\_BD}$,  $\mathit{u\_SD}$, $\mathit{u\_Orca}$, $\mathit{u\_ExtPtfm}$, $\mathit{u\_SrvD}$\footnote{Only if using ServoDyn Structural control with platform TMD and SubDyn.} and platform reference input from $\mathit{u\_ED}$}
   \State
   \State $u \gets \Call{u\_vec}{\mathit{u\_HD},\mathit{u\_SD},\mathit{u\_ED},\mathit{u\_BD},\mathit{u\_Orca},\mathit{u\_ExtPtfm}}$
   \State $k \gets 0$
   \Loop\Comment{Solve for loads and accelerations (direct feed-through terms)}
      \State $y\_ED   \gets \Call{ED\_CalcOutput}{\mathit{p\_ED},\mathit{u\_ED},\mathit{x\_ED},\mathit{xd\_ED},\mathit{z\_ED}}$
      \State $y\_SD   \gets \Call{SD\_CalcOutput}{\mathit{p\_SD},\mathit{u\_SD},\mathit{x\_SD},\mathit{xd\_SD},\mathit{z\_SD}}$
      \State $y\_HD   \gets \Call{HD\_CalcOutput}{\mathit{p\_HD},\mathit{u\_HD},\mathit{x\_HD},\mathit{xd\_HD},\mathit{z\_HD}}$
      \State $y\_BD   \gets \Call{BD\_CalcOutput}{\mathit{p\_BD},\mathit{u\_BD},\mathit{x\_BD},\mathit{xd\_BD},\mathit{z\_BD}}$
      \State $y\_Orca \gets \Call{Orca\_CalcOutput}{\mathit{p\_Orca},\mathit{u\_Orca},\mathit{x\_Orca},\mathit{xd\_Orca},\mathit{z\_Orca}}$
      \State $\mathit{y\_ExtPtfm}     \gets \Call{CalcOutput}{\mathit{p\_ExtPtfm},\mathit{u\_ExtPtfm},\mathit{x\_ExtPtfm},\mathit{xd\_ExtPtfm},\mathit{z\_ExtPtfm}}$ 
   
      \If{ $k \geq k\_max$}
         \State exit loop
      \EndIf
      
      \State$\mathit{u\_BD\_tmp}      \gets \Call{TransferMeshMotions}{y\_ED}$
      \State$\mathit{u\_MAP\_tmp}     \gets \Call{TransferMeshMotions}{y\_ED}$
      \State$\mathit{u\_FEAM\_tmp}    \gets \Call{TransferMeshMotions}{y\_ED}$
      \State$\mathit{u\_Orca\_tmp}    \gets \Call{TransferMeshMotions}{y\_ED}$
      \State$\mathit{u\_MD\_tmp}      \gets \Call{TransferMeshMotions}{y\_ED}$
      \State$\mathit{u\_IceF\_tmp}    \gets \Call{TransferMeshMotions}{y\_SD}$
      \State$\mathit{u\_IceD\_tmp(:)} \gets \Call{TransferMeshMotions}{y\_SD}$
      \State$\mathit{u\_HD\_tmp}      \gets \Call{TransferMeshMotions}{y\_ED,y\_SD}$
      \State$\mathit{u\_SrvD\_tmp}    \gets \Call{TransferMeshMotions}{y\_BD,y\_ED,y\_SD}$\footnote{Only if using ServoDyn Structural control.}
      \State$\mathit{u\_SD\_tmp}      \gets \!
            \begin{aligned}[t]
           & \Call{TransferMeshMotions}{\mathit{y\_ED}}  \\
                & \cup \Call{TransferMeshLoads}{}(\!
                   \begin{aligned}[t] 
                        & \mathit{y\_SD},                                  \\
                        & \mathit{y\_HD},       \mathit{u\_HD\_tmp},       \\
                        & \mathit{y\_IceF},     \mathit{u\_IceF\_tmp},     \\
                        & \mathit{y\_IceD(:)},  \mathit{u\_IceD\_tmp(:)},  \\
                     \end{aligned}
         \end{aligned}$
      \State$\mathit{u\_ED\_tmp} \gets \Call{TransferMeshLoads}{}( \!
                 \begin{aligned}[t]   & \mathit{y\_ED}, \\
                                      & \mathit{y\_HD},  \mathit{u\_HD\_tmp},  \\
                                      & \mathit{y\_SD},  \mathit{u\_SD\_tmp},  \\  
                                      & \mathit{y\_MAP}, \mathit{u\_MAP\_tmp}, \\ 
                                      & \mathit{y\_FEAM},\mathit{u\_FEAM\_tmp},\\
                                      & \mathit{y\_AD},  \mathit{u\_AD\_tmp}\footnote{Only if AD buoyancy at of hub enabled.},  \\  
                                      & \mathit{y\_SrvD},\mathit{u\_SrvD\_tmp}\footnote{Only if using ServoDyn Structural control.} )    % SrvD%PtfmStC only
                     \end{aligned}$

      \State
      \State$\mathit{U\_Residual} \gets u - \Call{u\_vec}{}( \!
                 \begin{aligned}[t]   & \mathit{u\_HD\_tmp}, \\
                                      & \mathit{u\_SD\_tmp}, \\
                                      & \mathit{u\_ED\_tmp}, \\
                                      & \mathit{u\_BD\_tmp}, \\
                                      & \mathit{u\_Orca\_tmp},\\
                                      & \mathit{u\_ExtPtfm\_tmp})
                     \end{aligned}$
      \State
      
      \If{ last Jacobian was calculated at least $\mathit{DT\_UJac}$ seconds ago }
         \State Calculate $\frac{\partial U}{\partial u}$
      \EndIf
      
      %\textit{
      %\State Perturb each input in $u$
      %\State Call $\Call{CalcOutput}{p,u,x,xd,z}$ for each module
      %\State Transfer perturbed outputs to inputs
      %\State Form new $u$
      %\State Compare new $u$ with $U\_Residual$
      %\State}

      \State Solve  $\frac{\partial U}{\partial u} \Delta u = - \mathit{U\_Residual}$ for $\Delta u$

      \State    
      \If{$\lVert \Delta u \rVert_2 < $ tolerance } \Comment{To be implemented later}
         \State exit loop
      \EndIf
      \State
      \State $u \gets u + \Delta u$
      \State Transfer $u$ to $\mathit{u\_HD}$, $\mathit{u\_SD}$, $\mathit{u\_BD}$, $\mathit{u\_Orca}$, $\mathit{u\_ExtPtfm}$, and $\mathit{u\_ED}$\Comment{loads and accelerations only}
      \State $k=k+1$
      
   \EndLoop   
   
   \State\Comment{Transfer non-acceleration fields to motion input meshes}
   \State 
   
   \State$\mathit{u\_BD}($not accelerations$) \gets \Call{TransferMeshMotions}{\mathit{y\_ED}}$
   \State$\mathit{u\_HD}($not accelerations$) \gets \Call{TransferMeshMotions}{\mathit{y\_ED},\mathit{y\_SD}}$
   \State$\mathit{u\_SD}($not accelerations$) \gets \Call{TransferMeshMotions}{\mathit{y\_ED}}$
   \State$\mathit{u\_Orca}($not accelerations$) \gets \Call{TransferMeshMotions}{\mathit{y\_ED}}$   \State$\mathit{u\_ExtPtfm}($not accelerations$) \gets \Call{TransferMeshMotions}{\mathit{y\_ED}}$
  \State 
   \State $\mathit{u\_MAP}     \gets \Call{TransferMeshMotions}{\mathit{y\_ED}}$
   \State $\mathit{u\_MD}      \gets \Call{TransferMeshMotions}{\mathit{y\_ED}}$
   \State $\mathit{u\_FEAM}    \gets \Call{TransferMeshMotions}{\mathit{y\_ED}}$
   \State $\mathit{u\_IceF}    \gets \Call{TransferMeshMotions}{\mathit{y\_SD}}$
   \State $\mathit{u\_IceD(:)} \gets \Call{TransferMeshMotions}{\mathit{y\_SD}}$
   \State $\mathit{u\_SrvD}    \gets \Call{TransferMeshMotions}{\mathit{y\_BD,y\_ED,y\_SD}}$\footnote{Only if using ServoDyn Structural control.}    % For SrvD%PtfmStC
         
\EndProcedure
\end{algorithmic}


\subsection {Implementation of line2-to-line2 loads mapping}
The inverse-lumping of loads is computed by a block matrix solve for the distributed forces and moments, 
using the following equation:

\begin{equation}
\label{EqLump}
   \begin{bmatrix}
   F^{DL} \\
   M^{DL} \\
   \end{bmatrix}
=
   \begin{bmatrix}
      A & 0 \\
      B & A \\   
   \end{bmatrix}
   \begin{bmatrix}
   F^{D} \\
   M^{D} \\
   \end{bmatrix}
\end{equation}

Because the forces do not depend on the moments, we first solve for the distributed forces, $F^D$:
\begin{equation}
\label{EqLumpF}
   \begin{bmatrix}   F^{DL} \\   \end{bmatrix}
=
   \left[      A    \right]
   \left[   F^{D}  \right]
\end{equation}

We then use the known values to solve for the distributed moments, $M^D$:
\begin{equation}
\label{EqLumpM1}
   \left[   M^{DL} \right]
=
   \begin{bmatrix} B & A \\   \end{bmatrix}
   \begin{bmatrix}
   F^{D} \\
   M^{D} \\
   \end{bmatrix}
= \left[   B \right] \left[   F^D \right] + \left[   A \right] \left[   M^D \right] 
\end{equation}
or
\begin{equation}
\label{EqLumpM2}
\left[   M^{DL} \right] - \left[   B \right] \left[   F^D \right] = \left[   A \right] \left[   M^D \right]
\end{equation}
Rather than store the matrix $B$, we directly perform the cross products that the matrix $B$ represents.
This makes the left-hand side of Equation \ref{EqLumpM2} known, leaving us with one matrix solve. This 
solve uses the same matrix $A$ used to obtain the distributed forces in Equation \ref{EqLumpF}; $A$ depends 
only on element reference positions and connectivity. We use 
the $LU$ factorization of matrix $A$ so that the second solve does not introduce much additional overhead.



\pagebreak %break here for now so that it doesn't look so strange
\section{Solve Option 2 Improvements}
\subsection {Input-Output Solves inside AdvanceStates}
This algorithm documents the procedure for advancing states with option 2 
Input-Output solves in FAST, assuming
all modules are in use. If an individual module is not in use during a particular
simulation, the calls to that module's subroutines are omitted and the module's 
inputs and outputs are neither set nor used.

\begin{algorithmic}[1]
\Procedure{FAST\_AdvanceStates}{\null}
\State $\Call{ED\_UpdateStates}{\mathit{p\_ED},\mathit{u\_ED},\mathit{x\_ED},\mathit{xd\_ED},\mathit{z\_ED}}$
\State $\mathit{y\_ED} \gets \Call{ED\_CalcOutput}{\mathit{p\_ED},\mathit{u\_ED},\mathit{x\_ED},\mathit{xd\_ED},\mathit{z\_ED}}$
\State
\State $\mathit{u\_BD}($hub and root motions$) \gets \Call{TransferOutputsToInputs}{\mathit{y\_ED}}$
\State $\Call{BD\_UpdateStates}{\mathit{p\_BD},\mathit{u\_BD},\mathit{x\_BD},\mathit{xd\_BD},\mathit{z\_BD}}$
\State $\mathit{y\_BD} \gets \Call{BD\_CalcOutput}{\mathit{p\_BD},\mathit{u\_BD},\mathit{x\_BD},\mathit{xd\_BD},\mathit{z\_BD}}$
\State
\State $\mathit{u\_AD}($not InflowWind$) \gets \Call{TransferOutputsToInputs}{\mathit{y\_ED},\mathit{y\_BD}}$
\State $\mathit{u\_IfW} \gets \Call{TransferOutputsToInputs}{\mathit{y\_ED},\mathit{y\_BD}$ at $\mathit{u\_AD}$ nodes$}$
\State $\Call{IfW\_UpdateStates}{\mathit{p\_IfW},\mathit{u\_IfW},\mathit{x\_IfW},\mathit{xd\_IfW},\mathit{z\_IfW}}$
\State $\mathit{y\_IfW} \gets \Call{IfW\_CalcOutput}{\mathit{u\_IfW}$ and other $\mathit{IfW}$ data structures$}$
\State
\State $\mathit{u\_AD}($InflowWind only$) \gets \Call{TransferOutputsToInputs}{\mathit{y\_IfW}}$
\State $\mathit{u\_SrvD} \gets \Call{TransferOutputsToInputs}{\mathit{y\_BD},\mathit{y\_ED},\mathit{y\_IfW},\mathit{y\_SD}}$
\State $\Call{AD\_UpdateStates}{\mathit{p\_AD},\mathit{u\_AD},\mathit{x\_AD},\mathit{xd\_AD},\mathit{z\_AD}}$
\State $\Call{SrvD\_UpdateStates}{\mathit{p\_SrvD},\mathit{u\_SrvD},\mathit{x\_SrvD},\mathit{xd\_SrvD},\mathit{z\_SrvD}}$
\State
\State All other modules (used in Solve Option 1) advance their states
\EndProcedure
\end{algorithmic}

Note that AeroDyn and ServoDyn outputs get calculated inside the ${CalcOutputs\_And\_SolveForInputs}$ routine. ElastoDyn, BeamDyn, and
InflowWind outputs do not get recalculated in ${CalcOutputs\_And\_SolveForInputs}$ except for the first time the routine is called
(because CalcOutput is called before UpdateStates at time 0).


\section {Linearization}
\subsection{Loads Transfer}
The loads transfer can be broken down into four components, all of which are used in the Line2-to-Line2 loads transfer:
\begin{enumerate}
  \item Augment the source mesh with additional nodes.
  \item Lump the distributed loads on the augmented Line2 source mesh to a Point mesh.
  \item Perform Point-to-Point loads transfer.
  \item Distribute (or "unlump") the point loads.   
\end{enumerate}
The other loads transfers are just subsets of the Line2-to-Line2 transfer:
\begin{itemize}
  \item Line2-to-Line2: Perform steps 1, 2, 3, and 4.
  \item Line2-to-Point: Perform steps 1, 2, and 3.
  \item Point-to-Line2: Perform steps 3 and 4.
  \item Point-to-Point: Perform step 3.
\end{itemize}


Each of the four steps can be represented with a linear equation. The linearization of the loads transfers is just multiplying the 
appropriate matrices generated in each of the steps.


\subsubsection{Step 1: Augment the source mesh}
The equation that linearizes mesh augmentation is
\begin{equation}
\label{Augment}
   \left\{   \begin{matrix} \vec{u}^D \\ \vec{u}^{SA} \\ \vec{f}^{SA} \\ \vec{m}^{SA} \end{matrix} \right\} 
=
   \begin{bmatrix}
   I_{\mathit{N_D}} & 0   & 0   & 0   \\
   0                      & M^A & 0   & 0   \\
   0                      & 0   & M^A & 0   \\
   0                      & 0   & 0   & M^A \\
   \end{bmatrix}
   \left\{   \begin{matrix} \vec{u}^D \\ \vec{u}^S \\ \vec{f}^S \\ \vec{m}^S \end{matrix} \right\} 
\end{equation}
where $M^A \in \mathbb{R}^{ \mathit{N_{SA}},\mathit{N_{S}}}$ indicates the mapping of nodes from the source mesh (with $N_S$ nodes) to the augmented source mesh
(with $N_{SA}$ nodes). The destination mesh (with $N_D$ nodes) is unchanged, as is indicated by matrix $I_{\mathit{N_D}}$.



\subsubsection{Step 2: Lump loads on a Line2 mesh to a Point mesh}
The equation that linearizes the lumping of loads is
\begin{equation}
\label{Lump}
   \left\{   \begin{matrix} \vec{u}^{SA} \\ \vec{F}^{SAL} \\ \vec{M}^{SAL} \end{matrix} \right\} 
=
   \begin{bmatrix}
    I_{\mathit{N_{SA}}}   & 0 \         & 0   \\
    0                     & M_{li}^{SL} & 0   \\
    M_{uS}^{SL}           & M_{f}^{SL}  & M_{li}^{SL} \\
   \end{bmatrix}
   \left\{   \begin{matrix} \vec{u}^{SA} \\ \vec{f}^{SA} \\ \vec{m}^{SA} \end{matrix} \right\} 
\end{equation}
where $M_{li}^{SL}, M_{uS}^{SL}, M_{f}^{SL} \in \mathbb{R}^{ \mathit{N_{SA}},\mathit{N_{SA}}}$ are block matrices that indicate the mapping of the lumped values to distributed values. $M_{li}^{SL}$ is matrix $A$ in Equation \ref{EqLumpF}, which depends only on element reference positions and connectivity. Matrices $M_{uS}^{SL}$and $M_{f}^{SL}$ also depend on values at their operating point.


\subsubsection{Step 3: Perform Point-to-Point loads transfer}
The equation that performs Point-to-Point load transfer can be written as
\begin{equation}
\label{P2P}
   \left\{   \begin{matrix} \vec{u}^D \\ \vec{u}^{S} \\ \vec{F}^{D} \\ \vec{M}^{D} \end{matrix} \right\} 
=
   \begin{bmatrix}
   I_{\mathit{N_D}} & 0                      & 0          & 0          \\
   0                & I_{\mathit{N_{S}}}     & 0          & 0          \\
   0                & 0                      & M_{li}^{D} & 0          \\
   M_{uD}^{D}       & M_{uS}^{D}             & M_{f}^{D}  & M_{li}^{D} \\
   \end{bmatrix}
   \left\{   \begin{matrix} \vec{u}^D \\ \vec{u}^{S}  \\ \vec{F}^{S}  \\ \vec{D}^{S} \end{matrix} \right\} 
\end{equation}
where $M_{li}^{D}, M_{uS}^{D}, M_{f}^{D} \in \mathbb{R}^{ \mathit{N_{D}},\mathit{N_{S}}}$ are block matrices that indicate the transfer of loads from one source
node to a node on the destination mesh. $M_{uD}^{D} \in \mathbb{R}^{ \mathit{N_{D}},\mathit{N_{D}}}$ is a diagonal matrix that indicates how the destination mesh's displaced position effects the transfer.


\subsubsection{Step 4: Distribute Point loads to a Line2 mesh}
Distributing loads from a Point mesh to a Line2 mesh is the inverse of step 2.

From Equation \ref{Lump} the equation that linearizes the lumping of loads on a destination mesh is
\begin{equation}
\label{LumpD}
   \left\{   \begin{matrix} \vec{u}^D \\ \vec{F}^{D} \\ \vec{M}^{D} \end{matrix} \right\} 
=
   \begin{bmatrix}
   I_{\mathit{N_D}} & 0           & 0   \\
   0                & M_{li}^{DL} & 0   \\
   M_{uD}^{DL}      & M_{f}^{DL}  & M_{li}^{DL} \\
   \end{bmatrix}
   \left\{   \begin{matrix} \vec{u}^D \\ \vec{f}^{D} \\ \vec{m}^{D} \end{matrix} \right\} 
\end{equation}
where $M_{li}^{DL}, M_{uD}^{DL}, M_{f}^{DL} \in \mathbb{R}^{ \mathit{N_{D}},\mathit{N_{D}}}$ are block matrices that indicate the mapping of the lumped values to distributed values. It follows that the inverse of this equation is
\begin{equation}
\label{InvLumpD}
   \left\{   \begin{matrix} \vec{u}^D \\ \vec{f}^{D} \\ \vec{m}^{D} \end{matrix} \right\} 
=
   \begin{bmatrix}
   I_{\mathit{N_D}} & 0                & 0   \\
   0                       & \left[ M_{li}^{DL} \right]^{-1} & 0   \\
   -\left[ M_{li}^{DL} \right]^{-1} M_{uD}^{DL} & 
   -\left[ M_{li}^{DL} \right]^{-1} M_{f}^{DL} \left[ M_{li}^{DL} \right]^{-1}  &
    \left[ M_{li}^{DL} \right]^{-1} \\
   \end{bmatrix}
   \left\{   \begin{matrix} \vec{u}^D \\ \vec{F}^{D} \\ \vec{M}^{D} \end{matrix} \right\} 
\end{equation}
The only inverse we need is already formed (stored as an LU decomposition) from the loads transfer, so we need not form it again.


\subsubsection{Putting it together}
To form the matrices for loads transfers for the various mappings available, we now need to multiply a few matrices to return the linearization
matrix that converts loads from the source mesh to loads on the line mesh:
\begin{equation}
\label{LinearEqn}
   \left\{   \begin{matrix} \vec{f}^{D} \\ \vec{m}^{D} \end{matrix} \right\} 
=
   \begin{bmatrix}
   0        & 0         & M_{li} & 0   \\
   M_{uD}   & M_{uS}    & M_f    & M_{li}
   \end{bmatrix}
   \left\{   \begin{matrix} \vec{u}^D \\ \vec{u}^S \\ \vec{f}^{D} \\ \vec{m}^{D} \end{matrix} \right\} 
\end{equation}




\begin{itemize}
  \item Line2-to-Line2: Perform steps 1, 2, 3, and 4.
         \begin{multline}
            \left\{   \begin{matrix} \vec{f}^{D} \\ \vec{m}^{D} \end{matrix} \right\} 
           =
            \begin{bmatrix}
            0                       & \left[ M_{li}^{DL} \right]^{-1} & 0   \\
            -\left[ M_{li}^{DL} \right]^{-1} M_{uD}^{DL} & 
            -\left[ M_{li}^{DL} \right]^{-1} M_{f}^{DL} \left[ M_{li}^{DL} \right]^{-1}  &
             \left[ M_{li}^{DL} \right]^{-1} \\
            \end{bmatrix}                  
            \\                                 
           \begin{bmatrix}
            I_{N_D}              & 0                      & 0          & 0          \\
            0                    & 0                      & M_{li}^{D} & 0          \\
            M_{uD}^{D}           & M_{uS}^{D}             & M_{f}^{D}  & M_{li}^{D} \\
            \end{bmatrix}            
            \begin{bmatrix}
            I_{N_D}              & 0                     & 0           & 0                \\
            0                    & I_{\mathit{N_{SA}}}   & 0           & 0                \\
            0                    & 0                     & M_{li}^{SL} & 0                \\
            0                    & M_{uS}^{SL}           & M_{f}^{SL}  & M_{li}^{SL}      \\
            \end{bmatrix}
            \\            
            \begin{bmatrix}
            I_{\mathit{N_D}}     & 0   & 0   & 0   \\
            0                    & M^A & 0   & 0   \\
            0                    & 0   & M^A & 0   \\
            0                    & 0   & 0   & M^A \\
            \end{bmatrix}
            \left\{   \begin{matrix} \vec{u}^D \\ \vec{u}^{S}  \\ \vec{f}^{S}  \\ \vec{m}^{S} \end{matrix} \right\}             
         \end{multline}            
      \begin{align}
       M_{li} &= \left(M_{li}^{DL}\right)^{-1}M_{li}^D M_{li}^{SL} M_A \\ 
        M_{uD} &= \left(M_{li}^{DL}\right)^{-1}\left[M_{uD}^D - M_{uD}^{DL}\right] \\
        M_{uS} &= \left(M_{li}^{DL}\right)^{-1} \left[ M_{uS}^D + M_{li}^D M_{uS}^{SL}\right] M_A \\
        M_{f}  &= \left(M_{li}^{DL}\right)^{-1}\left( \left[M_{f}^D - M_{f}^{DL}\left(M_{li}^{DL}\right)^{-1}M_{li}^D\right] M_{li}^{SL} +
                         M_{li}^D M_{f}^{SL}  \right)M_A \end{align} 
   
   
  \item Line2-to-Point: Perform steps 1, 2, and 3.
         \begin{multline}
            \left\{   \begin{matrix} \vec{F}^{D} \\ \vec{M}^{D} \end{matrix} \right\} 
           =
           \begin{bmatrix}
            0                    & 0                      & M_{li}^{D}  & 0            \\
            M_{uD}^{D}           & M_{uS}^{D}             & M_{f}^{D}   & M_{li}^{D}   \\
            \end{bmatrix}            
            \begin{bmatrix}
            I_{N_D}              & 0                     & 0            & 0            \\
            0                    & I_{\mathit{N_{SA}}}   & 0            & 0            \\
            0                    & 0                     & M_{li}^{SL}  & 0            \\
            0                    & M_{uS}^{SL}           & M_{f}^{SL}   & M_{li}^{SL}  \\
            \end{bmatrix}
            \\            
            \begin{bmatrix}
            I_{\mathit{N_D}}     & 0      & 0      & 0   \\
            0                    & M^A    & 0      & 0   \\
            0                    & 0      & M^A    & 0   \\
            0                    & 0      & 0      & M^A \\
            \end{bmatrix}
            \left\{   \begin{matrix} \vec{u}^D \\ \vec{u}^{S}  \\ \vec{f}^{S}  \\ \vec{m}^{S} \end{matrix} \right\}             
         \end{multline}               
      The linearization routine returns these four matrices:
   \begin{align}
    M_{li} &= M_{li}^D M_{li}^{SL} M_{A} \\ 
    M_{uD} &= M_{uD}^D \\ 
    M_{uS} &= \left[ M_{uS}^D + M_{li}^{D} M_{uS}^{SL}\right]M_{A} \\ 
    M_{f}  &= \left[ M_{f}^D M_{li}^{SL} + M_{li}^{D} M_{f}^{SL} \right]M_{A}  \end{align} 
   
   
  \item Point-to-Line2: Perform steps 3 and 4.
         %\begin{equation}
         \begin{multline}
            \left\{   \begin{matrix} \vec{f}^{D} \\ \vec{m}^{D} \end{matrix} \right\} 
           =
            \begin{bmatrix}
            0                       & \left[ M_{li}^{DL} \right]^{-1} & 0   \\
            -\left[ M_{li}^{DL} \right]^{-1} M_{uD}^{DL} & 
            -\left[ M_{li}^{DL} \right]^{-1} M_{f}^{DL} \left[ M_{li}^{DL} \right]^{-1}  &
             \left[ M_{li}^{DL} \right]^{-1} \\
            \end{bmatrix}                  \\                                 
          \begin{bmatrix}
            I_{N_D}             & 0                      & 0          & 0          \\
            0                   & 0                      & M_{li}^{D} & 0          \\
            M_{uD}^{D}          & M_{uS}^{D}             & M_{f}^{D}  & M_{li}^{D} \\
            \end{bmatrix}            
            \left\{   \begin{matrix} \vec{u}^D \\ \vec{u}^{S}  \\ \vec{F}^{S}  \\ \vec{M}^{S} \end{matrix} \right\}             
         \end{multline}            
         %\end{equation}   
      The linearization routine returns these four matrices:
   \begin{align} M_{li} &= \left(M_{li}^{DL}\right)^{-1}M_{li}^D \\
    M_{uD} &= \left(M_{li}^{DL}\right)^{-1}\left[M_{uD}^D - M_{uD}^{DL}\right] \\ 
    M_{uS} &= \left(M_{li}^{DL}\right)^{-1}M_{uS}^D \\ 
    M_{f}  &= \left(M_{li}^{DL}\right)^{-1}\left[M_{f}^D - M_{f}^{DL} M_{li} \right] \end{align} 

   
  \item Point-to-Point: Perform step 3.
            \begin{equation}
            \left\{   \begin{matrix} \vec{F}^{D} \\ \vec{M}^{D} \end{matrix} \right\} 
         =
            \begin{bmatrix}
            0                    & 0                     & M_{li}^{D} & 0          \\
            M_{uD}^{D}           & M_{uS}^{D}            & M_{f}^{D}  & M_{li}^{D} \\
            \end{bmatrix}
            \left\{   \begin{matrix} \vec{u}^D \\ \vec{u}^{S}  \\ \vec{F}^{S}  \\ \vec{M}^{S} \end{matrix} \right\} 
         \end{equation}      
   The linearization routine returns these four matrices:
   \begin{align} M_{li} &= M_{li}^D \\ 
    M_{uD} &= M_{uD}^D \\ 
    M_{uS} &= M_{uS}^D \\ 
    M_{f}  &= M_{f}^D  \end{align} 
   
\end{itemize}

\end{document}
