I will now explain what I understand to be the process of a finite element mass matrix assembly. I would like a reference which does something similar, or if I am mistaken about the process please let me know.
Let $\Omega$ be a bounded domain and let $\lbrace \Omega^{\left(e\right)}\rbrace_{e=1}^{E}$ be disjoint subsets of $\Omega$ which comprise a triangulation of $\Omega$. Let $\mathcal{V}$ be a finite element space over $\Omega$, with basis elements $\lbrace \phi_{i}\rbrace_{i=1}^{n}$. Let $M\in\mathbb{R}^{n\times n}$ be the mass-matrix
\begin{align} M_{i,j}=\sum_{e=1}^{E}\int_{\Omega^{e}}\phi_{i}\left(x\right)\phi_{j}\left(x\right)\mathrm{d}x. \end{align}
We then map the triangle $\Omega^{\left(e\right)}$ onto a fixed reference element $\Omega^{r}$ by a change of independent variable, achieved by a linear transformation.
\begin{align} M_{i,j}=\sum_{e=1}^{E}\int_{\Omega^{r}}\tilde{\phi}_{i}\left(x\right)\tilde{\phi}_{j}\left(x\right)|J^{\left(e\right)}|\mathrm{d}x. \end{align}
Because the transformation is linear the determinant of the element Jacobian $|J^{\left(e\right)}|$ is constant.
\begin{align} M_{i,j}=\sum_{e=1}^{E}|J^{\left(e\right)}|\int_{\Omega^{r}}\tilde{\phi}_{i}\left(x\right)\tilde{\phi}_{j}\left(x\right)\mathrm{d}x. \end{align}
Each element of the basis for $\mathcal{V}$ is assumed to be supported on a few elements $\Omega^{\left(e\right)}$, let $m$ be the number of basis functions with support on $\Omega^{\left(e\right)}$. $m$ is assumed to be independent of the element $\Omega^{\left(e\right)}$.
\begin{align} M_{i,j}=\sum_{e=1}^{E}|J^{\left(e\right)}|\int_{\Omega^{r}} \sum_{\ell,p=1}^{m} \delta_{i,e_{\ell}}\,\tilde{\phi}_{i}\left(x\right)\tilde{\phi}_{j}\left(x\right) \, \delta_{e_{p},j}\, \mathrm{d}x. \end{align}
where $\delta$ is the delta-Kronecker symbol and $\lbrace e_{p}\rbrace_{p=1}^{m}$ are the global indicies associated to a finite element function with support on $\Omega^{\left(e\right)}$. Let $\lbrace \psi_{i}\rbrace_{i=1}^{m}\subset \lbrace \tilde{\phi}_{i}\rbrace_{i=1}^{n}$ be the transformed basis functions that are supported on the reference element, those functions that are supported are element independent.
\begin{align} M_{i,j}=\sum_{e=1}^{E}|J^{\left(e\right)}|\int_{\Omega^{r}} \sum_{\ell,p=1}^{m} \delta_{i,e_{\ell}}\,\psi_{\ell}\left(x\right)\psi_{p}\left(x\right) \, \delta_{e_{p},j}\, \mathrm{d}x. \end{align}
\begin{align} M_{i,j}=\sum_{e=1}^{E}|J^{\left(e\right)}| \sum_{\ell,p=1}^{m} \delta_{i,e_{\ell}} \int_{\Omega^{r}} \psi_{\ell}\left(x\right)\psi_{p}\left(x\right) \mathrm{d}x\,\delta_{e_{p},j}. \end{align}
The only thing left to do is quadrature, so for that purpose let $\lbrace x_{k}\rbrace_{k=1}^{q}\subset\Omega^{r}$ be quadrature points and $\lbrace w_{k}\rbrace_{k=1}^{q}\subset\mathbb{R}$ be the quadrature weights.
\begin{align} M_{i,j}\approx\sum_{e=1}^{E}|J^{\left(e\right)}| \sum_{\ell,p=1}^{m} \delta_{i,e_{\ell}} \sum_{k=1}^{q} \psi_{\ell}\left(x_{k}\right)\,w_{k}\,\psi_{p}\left(x_{k}\right) \,\delta_{e_{p},j}. \end{align}
Let $G^{\left(e\right)}\in\mathbb{R}^{m\times n}$ be the global to local element matrix defined by $\left(G^{\left(e\right)}\right)_{p,j}=\delta_{e_{p},j}$, $B\in\mathbb{R}^{q\times m}$ be the local dof to basis functoin evaluation at the quadrature points $B_{k,p}=\psi_{p}\left(x_{k}\right)$ and $D^{\left(e\right)}\in\mathbb{R}^{q\times q}$ be a diagonal matrix which holds the quadrature weights $D^{\left(e\right)}_{i,j}=|J^{\left(e\right)}|\delta_{i,j}w_{i}$. With this we have
\begin{align} M\approx\sum_{e=1}^{E}\left(G^{\left(e\right)}\right)^{\top}\,B^{\top}\,D^{\left(e\right)}\,B\,G^{\left(e\right)} \end{align}
Again I am looking for textbooks which also describe this procedure or any feedback if I have any errors thanks!
