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!


2 Answers 2


You've got it the wrong way around, though the principles are right.

In your write-up, you loop over all $i$ and $j$ and then ask which cells you need to consider to compute a particular matrix entry $M_{ij}$. This is expensive because most matrix entries are zero, and because you would have to visit the same cell multiple times.

In practice, you loop over all cells and for each cell ask "which are the degrees of freedom that live here" and then compute the corresponding contributions of this cell to these matrix entries. You might want to look at how the usual finite element libraries do this. Here is the assembly done in deal.II: https://dealii.org/developer/doxygen/deal.II/step_3.html#Step3assemble_system

  • $\begingroup$ I see that looping over each $i,j$ can be expensive and one should take advantage of the sparsity to reduce the cost. $\endgroup$
    – Tucker
    Jun 2, 2021 at 20:58
  • $\begingroup$ The final expression for the mass matrix seems to be what you are suggesting, that is there is a loop over elements/cells, the action of $G^{\left(e\right)}$ extracts those elements that live on the cell, etc. $\endgroup$
    – Tucker
    Jun 2, 2021 at 20:59

I didn't go through your steps, I am kind of short on time. Though, I will check the math if I am free this evening. Here are my book suggestions.

"The Finite Element Method: Theory, Implementation, and Applications" by Larsson and Bengzon is beginner-friendly. Also "Understanding and implementing the finite element method" by Gockenbach is pretty good.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.