I'm following the derivation from Finite Element Method using Matlab 2nd Edition, pg 311-315, which derives of the local stiffness matrix for planar isotropic linear elasticity as follows:
Force Balance Equations
$\frac{\partial\sigma_x}{\partial x}+\frac{\tau_{xy}}{\partial y} + f_x=0$
$\frac{\partial\tau_{xy}}{\partial x}+\frac{\sigma_y}{\partial y} + f_y=0$
Using galerkin method, we multiply the first and second equation by test functions $w_1$ and $w_2$, respectively and integrate over the domain $\Omega$. Integrating by parts, I see that we obtain the weak formulation:
$$\int_\Omega \begin{bmatrix} \frac{\partial w_1}{\partial x} & 0 & \frac{\partial w_1}{\partial y}\\ 0 & \frac{\partial w_2}{\partial y} & \frac{\partial w_2}{\partial x} \end{bmatrix} \begin{bmatrix} \sigma_x \\ \sigma_y \\ \tau_{xy} \end{bmatrix} = \int_\Omega \begin{bmatrix} w_1 f_x \\ w_2 f_y \end{bmatrix} + \int_{\partial\Gamma} \begin{bmatrix} w_1 \Phi_x \\ w_2 \Phi_y \end{bmatrix}$$ where $\Gamma$ is the portion of the boundary with the neumann (traction) boundary condition in the x and y directions $\Phi_x$ and $\Phi_y$.
Let's just consider the integral on the left hand side of this equation. Using the linear isotropic stress strain relationship in two dimensions we can rewrite this equation as
$$\int_\Omega M D \epsilon$$
where
$M=\begin{bmatrix} \frac{\partial w_1}{\partial x} & 0 & \frac{\partial w_1}{\partial y}\\ 0 & \frac{\partial w_2}{\partial y} & \frac{\partial w_2}{\partial x}\end{bmatrix}$, $D=\frac{E}{1-\nu^2} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2} \end{bmatrix}$, and $\epsilon = \begin{bmatrix} \frac{\partial u}{\partial x} \\ \frac{\partial v}{\partial y} \\ \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \end{bmatrix}$.
Suppose the domain is tesselated into triangular elements. Consider a single element $e$ and the three basis functions (unit hat functions) on this element as $H_i(x)$ for i=1,2,3. We can charaterize the displacement functions on this element as $u(x,y)=\sum_{i=1}^3 u_iH_i$ and $v(x,y)=\sum_{i=1}^3 v_iH_i$, then we can rewrite
$$\epsilon=Bd$$
where $B=\begin{bmatrix} \frac{\partial H_1}{\partial x} & 0 & \frac{\partial H_2}{\partial x} & 0 & \frac{\partial H_3}{\partial x} & 0 \\ 0 & \frac{\partial H_1}{\partial y} & 0 & \frac{\partial H_2}{\partial y} & 0 & \frac{\partial H_3}{\partial y} \\ \frac{\partial H_1}{\partial y} & \frac{\partial H_1}{\partial x} & \frac{\partial H_2}{\partial y} & \frac{\partial H_2}{\partial x} & \frac{\partial H_3}{\partial y} & \frac{\partial H_3}{\partial x}\end{bmatrix}$ and $d=\begin{bmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ u_3 \\ v_3 \end{bmatrix}$.
Thus, we can write the integral over each element as
$$\int_e M D \epsilon= \int_e MDBd$$.
The author claims that the matrix $M$ becomes $B^T$ when we only consider the test functions equivalent to basis functions with support on element $e$. That is, when $w_1,w_2 = H_i$ for $i=1,2,3$ we obtain
$$\int_e B^TDBd$$
It's not immediately obvious to me why the matrix $M$ becomes $B^T$ over the element $e$. How can I arrive at this conclusion just by letting the test functions be the basis functions? Any help with this would be greatly appreciated! :)