We want to find matrices $X,Z$ such that
\begin{equation*}
\mathbf{y}=X\mathbf{\beta}+Z\mathbf{u}+\mathbf{\varepsilon}
\end{equation*}
where $\mathbf{\beta}$ contains the fixed effects/parameters and $\mathbf{u}$ contains the random effects.
What we get is the following
\begin{eqnarray*}
\left ( \begin{array}{c}
y_{111}\\y_{112}\\y_{221}\\y_{222}\\y_{331}\\y_{332}\\y_{441}\\y_{442}\\
y_{511}\\y_{512}\\y_{621}\\y_{622}\\y_{731}\\y_{732}\\y_{841}\\y_{842}\\

\end{array} \right )
&=&
\left ( \begin{array}{ccccccccccccccc}
1&1&0&0&0&1&0&1&0&0&0&0&0&0&0\\
1&1&0&0&0&0&1&0&1&0&0&0&0&0&0\\
1&0&1&0&0&1&0&0&0&1&0&0&0&0&0\\
1&0&1&0&0&0&1&0&0&0&1&0&0&0&0\\
1&0&0&1&0&1&0&0&0&0&0&1&0&0&0\\
1&0&0&1&0&0&1&0&0&0&0&0&1&0&0\\
1&0&0&0&1&1&0&0&0&0&0&0&0&1&0\\
1&0&0&0&1&0&1&0&0&0&0&0&0&0&1\\

1&1&0&0&0&1&0&1&0&0&0&0&0&0&0\\
1&1&0&0&0&0&1&0&1&0&0&0&0&0&0\\
1&0&1&0&0&1&0&0&0&1&0&0&0&0&0\\
1&0&1&0&0&0&1&0&0&0&1&0&0&0&0\\
1&0&0&1&0&1&0&0&0&0&0&1&0&0&0\\
1&0&0&1&0&0&1&0&0&0&0&0&1&0&0\\
1&0&0&0&1&1&0&0&0&0&0&0&0&1&0\\
1&0&0&0&1&0&1&0&0&0&0&0&0&0&1\\


\end{array} \right )
\left ( \begin{array}{c}
\mu\\\iota_1\\\iota_2\\\iota_3\\\iota_4\\\nu_1\\\nu_2\\(\iota\nu)_{11}\\(\iota\nu)_{12}\\(\iota\nu)_{21}\\(\iota\nu)_{22}\\(\iota\nu)_{31}\\(\iota\nu)_{32}\\(\iota\nu)_{41}\\(\iota\nu)_{42}\\
\end{array} \right )
\\ &+&
\left ( \begin{array}{cccccccc}
1&0&0&0&0&0&0&0\\
1&0&0&0&0&0&0&0\\
0&1&0&0&0&0&0&0\\
0&1&0&0&0&0&0&0\\
\multicolumn{8}{c}{\ddots}\\
0&0&0&0&0&0&0&1\\
0&0&0&0&0&0&0&1\\
\end{array} \right )
\left ( \begin{array}{c}
\phi_1\\\phi_2\\\phi_3\\\phi_4\\\phi_5\\\phi_6\\\phi_7\\\phi_8\\
\end{array} \right )
+
\mathbf{\varepsilon}.
\end{eqnarray*}
Because of the restrictions
\[
\sum_i \iota_i=\sum_v \nu_v = \sum_{i}(\iota\nu)_{iv}=\sum_v (\iota\nu)_{iv}=0
\]
we can reduce the $X$ matrix and write
\begin{eqnarray*}
\left ( \begin{array}{c}
y_{111}\\y_{112}\\y_{221}\\y_{222}\\y_{331}\\y_{332}\\y_{441}\\y_{442}\\
y_{511}\\y_{512}\\y_{621}\\y_{622}\\y_{731}\\y_{732}\\y_{841}\\y_{842}\\

\end{array} \right )
&=&
\left ( \begin{array}{cccccccc}
1&0&0&0&0&0&0&0\\
1&0&0&0&1&0&0&0\\
1&1&0&0&0&0&0&0\\
1&1&0&0&1&1&0&0\\
1&0&1&0&0&0&0&0\\
1&0&1&0&1&0&1&0\\
1&0&0&1&0&0&0&0\\
1&0&0&1&1&0&0&1\\

1&0&0&0&0&0&0&0\\
1&0&0&0&1&0&0&0\\
1&1&0&0&0&0&0&0\\
1&1&0&0&1&1&0&0\\
1&0&1&0&0&0&0&0\\
1&0&1&0&1&0&1&0\\
1&0&0&1&0&0&0&0\\
1&0&0&1&1&0&0&1\\


\end{array} \right )
\left ( \begin{array}{c}
\mu \\\iota_2\\\iota_3\\\iota_4\\\nu_2\\(\iota\nu)_{22}\\(\iota\nu)_{32}\\(\iota\nu)_{42}\\
\end{array} \right )
\\ &+&
\left ( \begin{array}{cccccccc}
1&0&0&0&0&0&0&0\\
1&0&0&0&0&0&0&0\\
0&1&0&0&0&0&0&0\\
0&1&0&0&0&0&0&0\\
\multicolumn{8}{c}{\ddots}\\
0&0&0&0&0&0&0&1\\
0&0&0&0&0&0&0&1\\
\end{array} \right )
\left ( \begin{array}{c}
\phi_1\\\phi_2\\\phi_3\\\phi_4\\\phi_5\\\phi_6\\\phi_7\\\phi_8\\
\end{array} \right )
+
\mathbf{\varepsilon}.
\end{eqnarray*}
In further calculations $X,Z,\mathbf{\beta},\mathbf{u}$ will be in the reduced form.