on the motion of the surface $\delta \h{\b{\varphi}} = \delta \b{\varphi}|_{\partial \c{B}_0}$, rendering the weak form as: $$\label{weak_combined}
\begin{split}
 0 &= \int_{\c{B}_0} \Grad  \delta\b{\varphi} : \b{P}(\b{F})  \, \d V 
     + \int_{\c{S}_0} \h{\Grad} \delta \h{\b{\varphi}} :\h{\b{P}}(\h{\b{F}})  \, \d A
    - \int_{\c{B}_0}  \delta \b{\varphi} \cdot \b{b}\prs_0 \, \d V 
    - \int_{\c{S}_0^\text{\tiny N}} \delta \h{\b{\varphi}} \cdot \h{\b{b}}\prs_0 \, \d A \\
    &\quad-\cancelto{0}{\biggl[\int_{\partial \c{B}_0} \delta\b{\varphi} \cdot [\b{P} \cdot \b{N}] \, \d A - \int_{\c{S}_0} \delta\h{\b{\varphi}} \cdot [\b{P} \cdot \b{N}] \, \d A \biggr]}
    \, , 
\end{split}$$ where the dependence of the Piola–Kirchhoff stress measures $\b{P}$ and $\h{\b{P}}$ on the solution is given via the constitutive relations in Summary \[constiutive\_relations\]. The Neumann part of the surface is denoted $\c{S}_0^\text{\tiny N}$. Note that from the assumption of a material surface $\delta\b{\varphi}\vert_{\c{S}_0} = \delta\h{\b{\varphi}}$, the final term of (\[weak\_combined\]) is thus zero due to traction continuity (as indicated). It will prove convenient, when constructing the finite element approximation, to decompose (\[weak\_combined\]) into volume and surface contributions. The volume and surface contributions will be recombined when the linear system is assembled.

The decomposed weak form of the governing equations is solved using a Newton–Raphson strategy. The time domain $[0,T]$ is decomposed into $N$ uniform intervals of duration $\Delta 