computation of the tangent in the volume and on the surface is used.

The resulting Newton scheme is thus given by $$\begin{aligned}
 R^{(i)} + \D_{\Delta \b{u}}R^{(i)}  + \h{R}^{}{(i)} + \D_{\Delta \h{\b{u}}}\h{R}{}^{(i)} & \equiv 0 \, , \notag \\
\implies \qquad \D_{\Delta \b{u}}R^{(i)}   + \D_{\Delta \h{\b{u}}}\h{R}{}^{(i)} & =  -R^{(i)} - \h{R}^{}{(i)} \, .
  \label{linear_combined}\end{aligned}$$

The fully-discrete problem {#sec_fully_discrete}
==========================

The material volume and the surface are partitioned into sets of non-overlapping cells (elements) individually denoted $\Omega_{e,0}$ and $\h{\Omega}_{e,0}$, respectively (see Fig. \[triangulation\]). The discrete form of the governing equations is obtained by approximating the test functions and trial solutions in the linearized weak form (\[linear\_combined\]) using vector-valued shape functions. The volume and surface shape functions, associated with an arbitrary degree of freedom in the volume $I$ and on the surface $\h{I}$, are denoted $\b{\Phi}^I$ and $\h{\b{\Phi}}{}^{\h{I}}$, respectively. The associated value of the degree of freedom (i.e. the increment in the motion) in the volume and on the surface are denoted $\mathsf{u}^I$ and $\h{\mathsf{u}}{}^{\h{I}}$, respectively. They form the entries of the global solution vectors $\mathsf{u}$ and $\h{\mathsf{u}}$. The trial solution in the volume and on the surface, and their gradients are approximated as follows[^1] $$\begin{aligned}
\b{u} \approx \sum_I \b{\Phi}^I \mathsf{u}^I    && \text{and}
    