Consider \( m+1 \) unknown functions: \( u^{(0)},\ldots, u^{(m)} \) governed by \( m+1 \) differential equations:
$$
\begin{align*}
\mathcal{L}_0(u^{(0)},\ldots,u^{(m)}) &= 0\\
&\vdots\\
\mathcal{L}_{m}(u^{(0)},\ldots,u^{(m)}) &= 0,
\end{align*}
$$
$$
\begin{align*}
\int_\Omega \mathcal{L}^{(0)}(u^{(0)},\ldots,u^{(m)}) v^{(0)}\dx &= 0\\
&\vdots\\
\int_\Omega \mathcal{L}^{(m)}(u^{(0)},\ldots,u^{(m)}) v^{(m)}\dx &= 0
\end{align*}
$$
Terms with second-order derivatives may be integrated by parts, with
Neumann conditions inserted in boundary integrals.
$$ V^{(i)} = \hbox{span}\{\basphi_0^{(i)},\ldots,\basphi_{N_i}^{(i)}\},$$
$$ u^{(i)} = B^{(i)}(\x) + \sum_{j=0}^{N_i} c_j^{(i)} \basphi_j^{(i)}(\x),
$$
Can derive \( m \) coupled linear systems for the unknowns \( c_j^{(i)} \), \( j=0,\ldots,N_i \), \( i=0,\ldots,m \).
The variational form is derived by taking the inner product of \( \boldsymbol{\mathcal{L}}(\u ) \) and \( \v \):
$$
\begin{equation*}
\int_\Omega \boldsymbol{\mathcal{L}}(\u )\cdot\v = 0\quad\forall\v\in\V
\end{equation*}
$$
$$
\begin{align*}
\mu \nabla^2 w &= -\beta\\
\kappa\nabla^2 T &= - \mu ||\nabla w||^2 \quad (= \mu \nabla w\cdot\nabla w)
\end{align*}
$$
Let \( w, (T-T_0) \in V \) with test functions \( v\in V \).
$$ V = \hbox{span}\{\basphi_0(x,y),\ldots,\basphi_N(x,y)\}, $$
$$
\begin{equation*}
w = \sum_{j=0}^N c^{(w)}_j \basphi_j,\quad T = T_0 +
\sum_{j=0}^N c^{(T)}_j\basphi_j
\end{equation*}
$$
Inserting \eqref{fem:sys:wT:ex:sum} in the PDEs, results in the residuals
$$
\begin{align*}
R_w &= \mu \nabla^2 w + \beta\\
R_T &= \kappa\nabla^2 T + \mu ||\nabla w||^2
\end{align*}
$$
Galerkin's method: make residual orthogonal to \( V \),
$$
\begin{align*}
\int_\Omega R_w v \dx &=0\quad\forall v\in V\\
\int_\Omega R_T v \dx &=0\quad\forall v\in V
\end{align*}
$$
Integrate by parts and use \( v=0 \) on \( \partial\Omega \) (Dirichlet conditions!):
$$
\begin{align*}
\int_\Omega \mu \nabla w\cdot\nabla v \dx &= \int_\Omega \beta v\dx
\quad\forall v\in V\\
\int_\Omega \kappa \nabla T\cdot\nabla v \dx &= \int_\Omega \mu
\nabla w\cdot\nabla w\, v\dx \quad\forall v\in V
\end{align*}
$$
$$ \int_{\Omega} (R_w, R_T)\cdot\v \dx = 0\quad\forall\v\in\V
$$
With \( \v = (v_0,v_1) \):
$$ \int_{\Omega} (R_w v_0 + R_T v_1) \dx = 0\quad\forall\v\in\V
$$
$$
\begin{equation*}
\int_\Omega (\mu\nabla w\cdot\nabla v_0 + \kappa\nabla T\cdot\nabla v_1)\dx
= \int_\Omega (\beta v_0 + \mu\nabla w\cdot\nabla w\, v_1)\dx,
\quad\forall \v\in\V
\end{equation*}
$$
Choosing \( v_0=v \) and \( v_1=0 \) gives the variational form \eqref{fem:sys:wT:ex:w:vf1}, while \( v_0=0 \) and \( v_1=v \) gives \eqref{fem:sys:wT:ex:T:vf1}.
$$
\begin{align*}
\mu (\nabla w,\nabla v) &= (\beta, v)
\quad\forall v\in V\\
\kappa(\nabla T,\nabla v) &= \mu(\nabla w\cdot\nabla w, v)\quad\forall v\in V
\end{align*}
$$
$$
\begin{align*}
\sum_{j=0}^N A^{(w)}_{i,j} c^{(w)}_j &= b_i^{(w)},\quad i=0,\ldots,N\\
\sum_{j=0}^N A^{(T)}_{i,j} c^{(T)}_j &= b_i^{(T)},\quad i=0,\ldots,N\\
A^{(w)}_{i,j} &= \mu(\nabla \basphi_j,\nabla\basphi_i)\\
b_i^{(w)} &= (\beta, \basphi_i)\\
A^{(T)}_{i,j} &= \kappa(\nabla \basphi_j,\nabla\basphi_i)\\
b_i^{(T)} &= (\mu\nabla w_{-}\cdot (\sum_k
c^{(w)}_k\nabla\basphi_k), \basphi_i)
\end{align*}
$$
Matrix-vector form (alternative notation):
$$
\begin{align*}
\mu K c^{(w)} &= b^{(w)}\\
\kappa K c^{(T)} &= b^{(T)}
\end{align*}
$$
where
$$
\begin{align*}
K_{i,j} &= (\nabla \basphi_j,\nabla \basphi_i)\\
b^{(w)} &= (b_0^{(w)},\ldots,b_{N}^{(w)})\\
b^{(T)} &= (b_0^{(T)},\ldots,b_{N}^{(T)})\\
c^{(w)} &= (c_0^{(w)},\ldots,c_{N}^{(w)})\\
c^{(T)} &= (c_0^{(T)},\ldots,c_{N}^{(T)})
\end{align*}
$$
First solve the system for \( c^{(w)} \), then solve the system for \( c^{(T)} \)
$$
\begin{align*}
\sum_{j=0}^N A^{(w,w)}_{i,j} c^{(w)}_j + \sum_{j=0}^N A^{(w,T)}_{i,j} c^{(T)}_j
&= b_i^{(w)},\quad i=0,\ldots,N,
\\
\sum_{j=0}^N A^{(T,w)}_{i,j} c^{(w)}_j + \sum_{j=0}^N A^{(T,T)}_{i,j} c^{(T)}_j &= b_i^{(T)},\quad i=0,\ldots,N,\\
A^{(w,w)}_{i,j} &= \mu(\nabla \basphi_j,\basphi_i)\\
A^{(w,T)}_{i,j} &= 0\\
b_i^{(w)} &= (\beta, \basphi_i)\\
A^{(w,T)}_{i,j} &= \mu(\nabla w_{-}\cdot\nabla\basphi_j), \basphi_i)\\
A^{(T,T)}_{i,j} &= \kappa(\nabla \basphi_j,\basphi_i)\\
b_i^{(T)} &= 0
\end{align*}
$$
$$
\begin{align*}
\mu K c^{(w)} &= b^{(w)}\\
L c^{(w)} + \kappa K c^{(T)} & =0
\end{align*}
$$
\( L \) is the matrix from the \( \nabla w_{-}\cdot\nabla \) operator:
\( L_{i,j} = A^{(w,T)}_{i,j} \).
Corresponding block form:
$$
\left(\begin{array}{cc}
\mu K & 0\\
L & \kappa K
\end{array}\right)
\left(\begin{array}{c}
c^{(w)}\\
c^{(T)}
\end{array}\right) =
\left(\begin{array}{c}
b^{(w)}\\
0
\end{array}\right)
$$
$$
\begin{align*}
V^{(w)} &= \hbox{span}\{\basphi_0^{(w)},\ldots,\basphi_{N_w}^{(w)}\}\\
V^{(T)} &= \hbox{span}\{\basphi_0^{(T)},\ldots,\basphi_{N_T}^{(T)}\}
\end{align*}
$$
$$
\begin{align*}
\int_\Omega \mu \nabla w\cdot\nabla v^{(w)} \dx &= \int_\Omega \beta v^{(w)}\dx
\quad\forall v^{(w)}\in V^{(w)}\\
\int_\Omega \kappa \nabla T\cdot\nabla v^{(T)} \dx &= \int_\Omega \mu
\nabla w\cdot\nabla w\, v^{(T)}\dx \quad\forall v^{(T)}\in V^{(T)}
\end{align*}
$$
Take the inner product with \( \v = (v^{(w)}, v^{(T)}) \) and integrate:
$$
\begin{equation*}
\int_\Omega (\mu\nabla w\cdot\nabla v^{(w)} +
\kappa\nabla T\cdot\nabla v^{(T)})\dx
= \int_\Omega (\beta v^{(w)} + \mu\nabla w\cdot\nabla w\, v^{(T)})\dx,
\end{equation*}
$$
valid \( \forall \v\in\V = V^{(w)}\times V^{(T)} \).