&& \Grad \b{u} \approx \sum_I \Grad \b{\Phi}^I \mathsf{u}^I \, , \\
\h{\b{u}} \approx \sum_{\h{I}} \h{\b{\Phi}}{}^{\h{I}} \h{\mathsf{u}}^{\h{I}}    && \text{and}
    && \h{\Grad} \h{\b{u}} \approx \sum_{\h{I}} \h{\Grad} \h{\b{\Phi}}{}^{\h{I}} \h{\mathsf{u}}^{\h{I}} \, .\end{aligned}$$ The same approximations are used for the test functions (i.e. a Bubnov–Galerkin spatial discretization is employed).

The resulting linearised incremental problem to be solved is given in matrix form as $$\begin{aligned}
\underbrace{\begin{bmatrix}
\mathsf{A}^{(i)} + \h{\mathsf{A}}^{(i)}
\end{bmatrix}}_{\overset{\square}{\mathsf{A}}}
\begin{bmatrix}
\Delta \mathsf{u}
\end{bmatrix}
=
-
\underbrace{\begin{bmatrix}
 \mathsf{R}^{(i)} + \h{\mathsf{R}}^{(i)}
\end{bmatrix}}_{\overset{\square}{\mathsf{R}}} \, , \label{newton_step}\end{aligned}$$ where the global residual vector $\overset{\square}{\mathsf{R}}$ and system matrix $\overset{\square}{\mathsf{A}}$ are assembled from the element contributions as follows: $$\begin{aligned}
\mathsf{R}^{(i)} = {\mathop{\mbox{\LARGE {\textsf {\textbf A}}}}}_e \sum_{I} \mathsf{R}^I_e
&& \text{and} &&
\mathsf{A}^{(i)} = {\mathop{\mbox{\LARGE {\textsf {\textbf A}}}}}_e \sum_{I}\sum_{J} \mathsf{A}^{IJ}_e \, , \\
\h{\mathsf{R}}^{(i)} = {\mathop{\mbox{\LARGE {\textsf {\textbf A}}}}}_e \sum_{\h{I}} \h{\mathsf{R}}^{\h{I}}_e
&& \text{and} &&
\h{\mathsf{A}}^{(i)} = {\mathop{\mbox{\LARGE {\textsf {\textbf A}}}}}_e \sum_{\h{I}}\sum_{\h{J}} \h{\mathsf{A}}^{\h{I}\h{J}}_e \, ,\end{aligned}$$ where $$\begin{aligned}
 \mathsf{R}^I_e &:= \int_{\Omega_{e,0}} \Grad  \b{\Phi}^I : \b{P}  \, \d V 
    - \int_{\Omega_{e,0}}   \b{\Phi}^I \cdot \b{b}\prs_0 \, \d V \, , \\
 \h{\mathsf{R}}^{\h{I}}_e &:=   \int_{\h{\Omega}_{e,0}} \h{\Grad} \h{\b{\Phi}}{}^{\h{I}} :\h{\b{P}}  \, \d A 
    - \int_{\Omega_{e,0}^\text{\tiny N}}  \h{\b{\Phi}}{}^{\h{I}} \cdot \h{\b{b}}\prs_0 \, \d A \, , \\
 \mathsf{A}^{IJ}_e &:= \int_{\Omega_{e,0}} \Grad  {\b{\Phi}}^I : \v{A} : \Grad {\b{\Phi}}^J \, \d V \, , \\
  