The problem is first solve for scalar field \(\theta \) with the interpolating polynomial \(a_{1}+a_{2}x+a_{3}y\). Writing\begin{equation} \theta =\begin{bmatrix} 1 & x & y \end{bmatrix}\begin{bmatrix} a_{1}\\ a_{2}\\ a_{3}\end{bmatrix} \tag{1} \end{equation} Evaluating the field \(\theta \) at each node gives\[\begin{bmatrix} \theta _{1}\\ \theta _{2}\\ \theta _{3}\end{bmatrix} =\begin{bmatrix} 1 & x_{1} & y_{1}\\ 1 & x_{2} & y_{2}\\ 1 & x_{3} & y_{3}\end{bmatrix}\begin{bmatrix} a_{1}\\ a_{2}\\ a_{3}\end{bmatrix} \] Hence\begin{align} \begin{bmatrix} a_{1}\\ a_{2}\\ a_{3}\end{bmatrix} & =\begin{bmatrix} 1 & x_{1} & y_{1}\\ 1 & x_{2} & y_{2}\\ 1 & x_{3} & y_{3}\end{bmatrix} ^{-1}\begin{bmatrix} \theta _{1}\\ \theta _{2}\\ \theta _{3}\end{bmatrix} \nonumber \\ & =\frac{1}{\Delta }\begin{bmatrix} x_{2}y_{3}-x_{3}y_{2} & x_{3}y_{1}-x_{1}y_{3} & x_{1}y_{2}-x_{2}y_{1}\\ y_{2}-y_{3} & y_{3}-y_{1} & y_{1}-y_{2}\\ x_{3}-x_{2} & x_{1}-x_{3} & x_{2}-x_{1}\end{bmatrix}\begin{bmatrix} \theta _{1}\\ \theta _{2}\\ \theta _{3}\end{bmatrix} \tag{2} \end{align}
Where \(\Delta \) is the determinant \(x_{1}y_{2}-x_{2}y_{1}-x_{1}y_{3}+x_{3}y_{1}+x_{2}y_{3}-x_{3}y_{2}\). Substituting (2) into (1) gives\begin{align} \theta & =\overset{\mathbf{N}}{\overbrace{\frac{1}{\Delta }\begin{bmatrix} 1 & x & y \end{bmatrix}\begin{bmatrix} x_{2}y_{3}-x_{3}y_{2} & x_{3}y_{1}-x_{1}y_{3} & x_{1}y_{2}-x_{2}y_{1}\\ y_{2}-y_{3} & y_{3}-y_{1} & y_{1}-y_{2}\\ x_{3}-x_{2} & x_{1}-x_{3} & x_{2}-x_{1}\end{bmatrix} }}\begin{bmatrix} \theta _{1}\\ \theta _{2}\\ \theta _{3}\end{bmatrix} \nonumber \\ & =\begin{bmatrix} N_{1} & N_{2} & N_{3}\end{bmatrix}\begin{bmatrix} \theta _{1}\\ \theta _{2}\\ \theta _{3}\end{bmatrix} \tag{3} \end{align}
Where \begin{align} N_{1} & =\frac{1}{\Delta }\left [ x_{2}y_{3}-x_{3}y_{2}+x\left ( y_{2}-y_{3}\right ) +y\left ( x_{3}-x_{2}\right ) \right ] \tag{4}\\ N_{2} & =\frac{1}{\Delta }\left [ x_{3}y_{1}-x_{1}y_{3}+x\left ( y_{3}-y_{1}\right ) +y\left ( x_{1}-x_{3}\right ) \right ] \nonumber \\ N_{3} & =\frac{1}{\Delta }\left [ x_{1}y_{2}-x_{2}y_{1}+x\left ( y_{1}-y_{2}\right ) +y\left ( x_{2}-x_{1}\right ) \right ] \nonumber \end{align}
For constant stress triangle, the field is a vector field. Hence replacing \(\theta \) with \(\begin{bmatrix} u\\ v \end{bmatrix} \) equation (3) becomes\[\begin{bmatrix} u\\ v \end{bmatrix} =\begin{bmatrix} N_{1} & 0 & N_{2} & 0 & N_{3} & 0\\ 0 & N_{1} & 0 & N_{2} & 0 & N_{3}\end{bmatrix}\begin{bmatrix} u_{1}\\ v_{1}\\ u_{2}\\ v_{2}\\ u_{3}\\ v_{3}\end{bmatrix} \] From the above \begin{align*} \frac{\partial u}{\partial x} & =\frac{\partial N_{1}}{\partial x}u_{1}+\frac{\partial N_{2}}{\partial x}u_{2}+\frac{\partial N_{3}}{\partial x}u_{3}\\ \frac{\partial v}{\partial y} & =\frac{\partial N_{1}}{\partial y}v_{1}+\frac{\partial N_{2}}{\partial y}v_{2}+\frac{\partial N_{3}}{\partial y}v_{3}\\ \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} & =\frac{\partial N_{1}}{\partial y}u_{1}+\frac{\partial N_{2}}{\partial y}u_{2}+\frac{\partial N_{3}}{\partial y}u_{3}+\frac{\partial N_{1}}{\partial x}v_{1}+\frac{\partial N_{2}}{\partial x}v_{2}+\frac{\partial N_{3}}{\partial x}v_{3} \end{align*}
Hence\begin{align} \begin{bmatrix} \varepsilon _{x}\\ \varepsilon _{y}\\ \gamma _{xy}\end{bmatrix} & =\begin{bmatrix} \frac{\partial }{\partial x} & 0\\ 0 & \frac{\partial }{\partial y}\\ \frac{\partial }{\partial y} & \frac{\partial }{\partial x}\end{bmatrix}\begin{bmatrix} u\\ v \end{bmatrix} \nonumber \\ & =\begin{bmatrix} \frac{\partial }{\partial x}u\\ \frac{\partial }{\partial y}v\\ \frac{\partial }{\partial y}u+\frac{\partial }{\partial x}v \end{bmatrix} \nonumber \\ & =\begin{bmatrix} \frac{\partial N_{1}}{\partial x}u_{1}+\frac{\partial N_{2}}{\partial x}u_{2}+\frac{\partial N_{3}}{\partial x}u_{3}\\ \frac{\partial N_{1}}{\partial y}v_{1}+\frac{\partial N_{2}}{\partial y}v_{2}+\frac{\partial N_{3}}{\partial y}v_{3}\\ \frac{\partial N_{1}}{\partial y}u_{1}+\frac{\partial N_{2}}{\partial y}u_{2}+\frac{\partial N_{3}}{\partial y}u_{3}+\frac{\partial N_{1}}{\partial x}v_{1}+\frac{\partial N_{2}}{\partial x}v_{2}+\frac{\partial N_{3}}{\partial x}v_{3}\end{bmatrix} \nonumber \\ & =\overset{B}{\overbrace{\begin{bmatrix} \frac{\partial N_{1}}{\partial x} & 0 & \frac{\partial N_{2}}{\partial x} & 0 & \frac{\partial N_{3}}{\partial x} & 0\\ 0 & \frac{\partial N_{1}}{\partial y} & 0 & \frac{\partial N_{2}}{\partial y} & 0 & \frac{\partial N_{3}}{\partial y}\\ \frac{\partial N_{1}}{\partial y} & \frac{\partial N_{1}}{\partial x} & \frac{\partial N_{2}}{\partial y} & \frac{\partial N_{2}}{\partial x} & \frac{\partial N_{3}}{\partial y} & \frac{\partial N_{3}}{\partial x}\end{bmatrix} }}\begin{bmatrix} u_{1}\\ v_{1}\\ u_{2}\\ v_{2}\\ u_{3}\\ v_{3}\end{bmatrix} \tag{5} \end{align}
From (4) all of the \(\frac{\partial N_{i}}{\partial x},\frac{\partial N_{j}}{\partial y}\)terms are evaluated. Substituting the result into (5) gives the \(B\) matrix\begin{align*} \frac{\partial N_{1}}{\partial x} & =\frac{1}{\Delta }\left ( y_{2}-y_{3}\right ) \\ \frac{\partial N_{2}}{\partial x} & =\frac{1}{\Delta }\left ( y_{3}-y_{1}\right ) \\ \frac{\partial N_{3}}{\partial x} & =\frac{1}{\Delta }\left ( y_{1}-y_{2}\right ) \end{align*}
And\begin{align*} \frac{\partial N_{1}}{\partial y} & =\frac{1}{\Delta }\left ( x_{3}-x_{2}\right ) \\ \frac{\partial N_{2}}{\partial y} & =\frac{1}{\Delta }\left ( x_{1}-x_{3}\right ) \\ \frac{\partial N_{3}}{\partial y} & =\frac{1}{\Delta }\left ( x_{2}-x_{1}\right ) \end{align*}
Hence \(B\) becomes\begin{equation} B=\frac{1}{\Delta }\begin{bmatrix} y_{2}-y_{3} & 0 & y_{3}-y_{1} & 0 & y_{1}-y_{2} & 0\\ 0 & x_{3}-x_{2} & 0 & x_{1}-x_{3} & 0 & x_{2}-x_{1}\\ x_{3}-x_{2} & y_{2}-y_{3} & x_{1}-x_{3} & y_{3}-y_{1} & x_{2}-x_{1} & y_{1}-y_{2}\end{bmatrix} \tag{6} \end{equation} Letting \(y_{i}-y_{j}=y_{ij}\) and \(x_{i}-x_{j}=x_{ij}\), the above becomes\[ B=\frac{1}{\Delta }\begin{bmatrix} y_{23} & 0 & y_{31} & 0 & y_{12} & 0\\ 0 & x_{32} & 0 & x_{13} & 0 & x_{21}\\ x_{32} & y_{23} & x_{13} & y_{31} & x_{21} & y_{12}\end{bmatrix} \] But the area of triangle is given by\begin{align*} A & =\frac{1}{2}\begin{vmatrix} \mathbf{i} & \mathbf{j} & \mathbf{k}\\ x_{2}-x_{1} & y_{2}-y_{1} & 0\\ x_{3}-x_{1} & y_{3}-y_{1} & 0 \end{vmatrix} \\ 2A & =\left ( x_{2}-x_{1}\right ) \left ( y_{3}-y_{1}\right ) -\left ( y_{2}-y_{1}\right ) \left ( x_{3}-x_{1}\right ) \\ & =x_{1}y_{2}-x_{2}y_{1}-x_{1}y_{3}+x_{3}y_{1}+x_{2}y_{3}-x_{3}y_{2} \end{align*}
And the determinant \(\Delta \) was found above to be \(x_{1}y_{2}-x_{2}y_{1}-x_{1}y_{3}+x_{3}y_{1}+x_{2}y_{3}-x_{3}y_{2}\), hence\[ 2A=\Delta \] Substituting the above into \(B\) found above in equation (6) gives\[ B=\frac{1}{2A}\begin{bmatrix} y_{2}-y_{3} & 0 & y_{3}-y_{1} & 0 & y_{1}-y_{2} & 0\\ 0 & x_{3}-x_{2} & 0 & x_{1}-x_{3} & 0 & x_{2}-x_{1}\\ x_{3}-x_{2} & y_{2}-y_{3} & x_{1}-x_{3} & y_{3}-y_{1} & x_{2}-x_{1} & y_{1}-y_{2}\end{bmatrix} \]