$\newcommand{\Dev}{\operatorname{Dev}}$ $\newcommand{\dev}{\operatorname{dev}}$ $\newcommand{\grad}{\operatorname{grad}}$ $\newcommand{\Grad}{\operatorname{Grad}}$ $\newcommand{\divg}{\operatorname{div}}$ $\newcommand{\Ei}{\operatorname{Ei}}$ $\newcommand{\tr}{\operatorname{tr}}$ $\newcommand{\Divg}{\operatorname{Div}}$ $\newcommand{\cay}{\operatorname{cay}}$ $\newcommand{\rot}{\operatorname{rot}}$

Consider a spring between points $\mathbf{x}^{\left(i\right)}$ on rigid bodies $\left(1\right)$ and $\left(2\right)$ . The spring force is oriented along $\mathbf{n}$ , where $$\begin{equation} \mathbf{n}=\frac{\mathbf{x}^{\left(2\right)}-\mathbf{x}^{\left(1\right)}}{\left|\mathbf{x}^{\left(2\right)}-\mathbf{x}^{\left(1\right)}\right|}\,.\label{eq:rbs-1} \end{equation}$$ We may write $$\begin{equation} L=\left|\mathbf{x}^{\left(2\right)}-\mathbf{x}^{\left(1\right)}\right|\,,\quad L_{0}=\left|\mathbf{X}^{\left(2\right)}-\mathbf{X}^{\left(2\right)}\right|\,.\label{eq:rbs-2} \end{equation}$$ Then, the spring force acting on body $\left(1\right)$ is $$\begin{equation} \mathbf{f}^{\left(1\right)}=k\left(L-L_{0}\right)\mathbf{n}=k\left(\mathbf{x}^{\left(2\right)}-\mathbf{x}^{\left(1\right)}-L_{0}\mathbf{n}\right)\,,\label{eq:rbs-3} \end{equation}$$ where $k$ is the spring constant. It is understood that the special case $L_{0}=0$ produces $\mathbf{f}^{\left(1\right)}=k\left(\mathbf{x}^{\left(2\right)}-\mathbf{x}^{\left(1\right)}\right)$ . We do not need to use the augmented Lagrangian method here, therefore $\mathbf{f}^{\left(1\right)}=k\mathbf{c}$ , where $$\begin{equation} \mathbf{c}=\mathbf{x}_{n+\alpha}^{\left(2\right)}-\mathbf{x}_{n+\alpha}^{\left(1\right)}-L_{0}\frac{\mathbf{x}_{n+\alpha}^{\left(2\right)}-\mathbf{x}_{n+\alpha}^{\left(1\right)}}{\left|\mathbf{x}_{n+\alpha}^{\left(2\right)}-\mathbf{x}_{n+\alpha}^{\left(1\right)}\right|}\,.\label{eq:rbs-4} \end{equation}$$ Since the spring may freely pivot about its insertion point, we let $\mathbf{m}^{\left(1\right)}=\mathbf{0}$ . Thus, $$\begin{equation} \boxed{\begin{aligned}\mathbf{f}^{\left(1\right)} & =k\left(\mathbf{x}_{n+\alpha}^{\left(2\right)}-\mathbf{x}_{n+\alpha}^{\left(1\right)}-L_{0}\frac{\mathbf{x}_{n+\alpha}^{\left(2\right)}-\mathbf{x}_{n+\alpha}^{\left(1\right)}}{\left|\mathbf{x}_{n+\alpha}^{\left(2\right)}-\mathbf{x}_{n+\alpha}^{\left(1\right)}\right|}\right)\\ \mathbf{m}^{\left(1\right)} & =\mathbf{0} \end{aligned} }\label{eq:rbs-5} \end{equation}$$ and the virtual work associated with this spring is $$\begin{equation} \delta G=\left[\begin{array}{cccc} \delta\mathbf{r}^{\left(1\right)} & \delta\boldsymbol{\omega}^{\left(1\right)} & \delta\mathbf{r}^{\left(2\right)} & \delta\boldsymbol{\omega}^{\left(2\right)}\end{array}\right]\left[\begin{array}{c} \mathbf{f}^{\left(1\right)}\\ \hat{\mathbf{z}}_{n+\alpha}^{\left(1\right)}\cdot\mathbf{f}^{\left(1\right)}\\ -\mathbf{f}^{\left(1\right)}\\ -\hat{\mathbf{z}}_{n+\alpha}^{\left(2\right)}\cdot\mathbf{f}^{\left(1\right)} \end{array}\right]\,.\label{eq:rbs-6} \end{equation}$$ The linearization of $\delta G$ may be expressed as $$\begin{equation} -D\delta G=\left[\begin{array}{cccc} \delta\mathbf{r}^{\left(1\right)} & \delta\boldsymbol{\omega}^{\left(1\right)} & \delta\mathbf{r}^{\left(2\right)} & \delta\boldsymbol{\omega}^{\left(2\right)}\end{array}\right]\left[\mathbf{K}\right]\left[\begin{array}{c} \Delta\mathbf{r}^{\left(1\right)}\\ \Delta\boldsymbol{\theta}^{\left(1\right)}\\ \Delta\mathbf{r}^{\left(2\right)}\\ \Delta\boldsymbol{\theta}^{\left(2\right)} \end{array}\right]\,,\label{eq:rbs-7} \end{equation}$$ where the stiffness matrix is $$\begin{equation} \left[\mathbf{K}\right]=k\alpha\left[\begin{array}{cccc} \mathbf{P} & -\mathbf{P}\cdot\hat{\mathbf{z}}_{n+1}^{\left(1\right)} & -\mathbf{P} & \mathbf{P}\cdot\hat{\mathbf{z}}_{n+1}^{\left(2\right)}\\ \hat{\mathbf{z}}_{n+\alpha}^{\left(1\right)}\cdot\mathbf{P} & -\hat{\mathbf{z}}_{n+\alpha}^{\left(1\right)}\cdot\mathbf{P}\cdot\hat{\mathbf{z}}_{n+1}^{\left(1\right)} & -\hat{\mathbf{z}}_{n+\alpha}^{\left(1\right)}\cdot\mathbf{P} & \hat{\mathbf{z}}_{n+\alpha}^{\left(1\right)}\cdot\mathbf{P}\cdot\hat{\mathbf{z}}_{n+1}^{\left(2\right)}\\ -\mathbf{P} & \mathbf{P}\cdot\hat{\mathbf{z}}_{n+1}^{\left(1\right)} & \mathbf{P} & -\mathbf{P}\cdot\hat{\mathbf{z}}_{n+1}^{\left(2\right)}\\ -\hat{\mathbf{z}}_{n+\alpha}^{\left(2\right)}\cdot\mathbf{P} & \hat{\mathbf{z}}_{n+\alpha}^{\left(2\right)}\cdot\mathbf{P}\cdot\hat{\mathbf{z}}_{n+1}^{\left(1\right)} & \hat{\mathbf{z}}_{n+\alpha}^{\left(2\right)}\cdot\mathbf{P} & -\hat{\mathbf{z}}_{n+\alpha}^{\left(2\right)}\cdot\mathbf{P}\cdot\hat{\mathbf{z}}_{n+1}^{\left(2\right)} \end{array}\right]\,,\label{eq:rbs-8} \end{equation}$$ with $$\begin{equation} \mathbf{P}=\left(1-\frac{L_{0}}{L_{n+\alpha}}\right)\mathbf{I}+\frac{L_{0}}{L_{n+\alpha}}\mathbf{n}_{n+\alpha}\otimes\mathbf{n}_{n+\alpha}\,.\label{eq:rbs-9} \end{equation}$$