Finite element assembly

We now have many ingredients to assemble the finite element system at our disposal. We have the basis functions, the quadrature rule, the mesh, ... .

We will split the setup of the linear system in two parts:

  • the computation of local contributions on each element ("form integrals") and
  • putting together the local contributions to the global system (the "assembly").

Example: Mass matrix

We will start with the mass matrix, which is the simplest example. The mass matrix is given by $$ M_{ij} = \int_\Omega \phi_i \phi_j \, dx $$ The first step is to compute the local contributions. For the mass matrix, this is given by $$ M_{ij} = \int_{\Omega} \phi_i \phi_j \, dx = \sum_{T \in \mathcal{T}_h} \int_T \phi_i \phi_j \, dx $$

As most basis functions are zero outside of their element, we can restrict the integration on the element $T$ to basis functions that are non-zero on $T$.

For every $i$ so that $\phi_i|_T \neq 0$ we have that $\phi_i = \phi_{m(j)}$ for some $m(j)$ where $m$ is the local-to-global map (element_dofs of FESpace).

Instead of running of $i$ and $j$ for all dofs, the matrix $M$ is rather filled by running of all elements $T$ and collecting all contributions of the basis functions that are non-zero on $T$ first and then adding them to the matrix $M$.

$$ M = \sum_{T \in \mathcal{T}_h} E_T^T M_T E_T $$ where

  • $M_T \in \mathbb{R}^{n_{local} \times n_{local}}$ is the local mass matrix on $T$ ("form integrals") and
  • $E_T \in \mathbb{R}^{n_{local} \times n_{global}}$ is "connectivity" matrix (never used as a matrix though!) with $(E_T)_{i,j} = \delta_{m(i),j}$, i.e. the matrix that associates the global dof $i$ with the local dof $j$ according to the local-to-global map $m(j)$. (part of "assembly")

Computing the element matrix

We will now turn to the realization of $(M_T)$. To this end we note that with $\bar \phi_i = \phi_{m(i)}$ we have $$ (M_T)_{i,j} = \int_T \phi_{m(i)} \phi_{m(j)} \, dx = \int_T \bar \phi_{i} \bar \phi_{j} \, dx $$ Now, we want to translate this to the reference element and note that for most finite elements we have $\bar \phi_i(x) = \hat \phi_i(\hat x)$ with $\hat x = \Phi_T^{-1} (x)$ (or $x = \Phi_T(\hat x))$) where:

  • $\hat \phi_i$ is the basis function on the reference element $\hat T$ and
  • $\Phi_T: \hat T \to T$ is the (diffeomorphic) transformation from the reference element $\hat T$ to $T$.

This yields $$ (M_T)_{i,j} = \int_T \phi_{m(i)} \phi_{m(j)} \, dx = \int_T (\hat \phi_{i} \hat \phi_{j})(\Phi_T^{-1}(x)) \, dx = \int_{\hat T} |\operatorname{det}(D \Phi_T)| \hat \phi_{i} \hat \phi_{j} \, d \hat x. $$

Finally, we replace the integral with numerical quadrature (with abuse of notation with respect to $M^T$): $$ (M_T)_{i,j} = \sum_k^{n_{\text{quad}}} \hat \omega_k |\operatorname{det}(D \Phi_T)|(\hat x_k) \hat \phi_{i}(\hat x_k) \hat \phi_{j}(\hat x_k) $$ where $\{(\hat x_k, \hat \omega_k)\}_k$ is the quadrature rule for $\hat T$ (with sufficient accuracy).

Computing the (small and dense) matrix $M^T$ for a given element is the task of a FormIntegral, cf. formint.py. It needs to perform the following tasks:

  • determine requirements on the quadrature rule and pick one
  • evaluate the shape functions at the quadrature poins (for other form integrals derivatives or other differentials of the shape functions may be evaluated here)
  • evaluate transformation terms like $|\operatorname{det}(D \Phi_T)|$ at the quadrature points (for other form integrals the full Jacobian may also be needed)
  • sum all terms up

Depending on the concrete FormIntegral several matrix- or even higher order tensor terms may be set up as intermediate quantities. On the numpy level a convenient tool to apply tensor reduction (summing up on different axes) efficiently is the einsum function.

Next,

  • we will take a look at the (already mostly finished) implementation of the assemblies in forms.py (if possible together) and
  • afterwards turn to the implementation of the form integrals in the next tasks:

Task FormInt-1: The MassIntegral

Complete the implementation of the MassIntegral class in formint.py and test it with the test in test_formint_1d.py. You can also use the tests to see how the MassIntegral is called.

Note:

  • The implemenation of the MassIntegral is not very efficient can be computed dimension-independent, i.e. for the 1D and 2D (triangle) case at once. However, you can also implement it dimension-specific, if it helps.
  • Note that the FormIntegral implementation should also be (where possible) independent of the choice of finite element. However, the selection of the quadrature rule may be chosen in dependence of the order of the finite element.

Task FormInt-2: The SourceIntegral

Next, complete the implementation of the SourceIntegral class in formint.py and test it with the test in test_formint_1d.py and test_formint_2d.py. Note that now the source integral is a linear form integral corresponding to a given MeshFunction of the form $f_i = \int_T f(x) \phi_i(x) \, dx = \int_T |\operatorname{det}(D \Phi_T)| f(\Phi_T(\hat x)) \hat \phi_i (\hat x) \, d \hat x$.

Task FormInt-3: The LaplaceIntegral

To gather all components for the solution of Poisson(-like) problemc, we need to implement the Laplace term. Complete the implementation of the LaplaceIntegral class in formint.py and test it with the test in test_formint_1d.py. To this end, first write out the steps according to the steps above and be careful with transformations of derivatives of the shape functions (chain rule!).

Task Assembly-1:

Implement the solution of the following problem. Let $V_h$ be a proper finite element space (any of the implemented one is ok). Find $u_h \in V_h$ such that $$ \int_\Omega u_h v_h \, dx = \int_{\Omega} w v_h \, dx \quad \forall v_h \in V_h. $$

Write a demo (e.g. a jupyter notebook in demos) that solves the problem for a given $w$ and plots the solution $u_h$.

Write a test that checks that if $w \in V_h$ that $u_h = w$ exactly.