Quadrature over a triangle. Computing the dot products of the gradients

Question

From the book "Understanding the Finite Element Method" by Mark S. Gockenbach (Chapter 7), there is stated an integral:

$$ \int_{T_k} \kappa \nabla \phi_i \cdot \nabla\phi_j. $$

for a triangle denoted $T_k$. The gradient $\phi_i$ is denoted as:

$$\nabla \phi_i = \begin{bmatrix} b_i \\ c_i \end{bmatrix} {on} \hphantom a {T_k}$$

There is a matrix $M$ defined as: $$ \begin{bmatrix} 1 & x1 & y1; \\ 1 & x2 & y2; \\ 1 & x3 & y3; \\ \end{bmatrix} $$

Steps are given:

1. Estimate the integral $T_k$ over $\kappa$.

2. Invert the 3x3 matrix and extract the gradients of the basis functions from $M^{-1}$:

3. Compute the dot products of the gradients

$$\nabla \phi_i \cdot \nabla \phi_j, \quad i=1,2,3,\ j=i,\ldots,3$$

I don't really understand it yet. Given a triangle $(0,1),(-1,-2),(1,-3)$ how can $ \int_{T_k} \kappa \nabla \phi_i \cdot \nabla\phi_j. $ be computed?

How are the gradients of the basis functions extracted from $M^{-1}$?

Answer 1

The key piece here is that the basis functions on this triangle are linear, so $$\phi_i(x,y) = a_i + b_ix + c_iy.$$ Now notice the following:

1. The gradient of $\phi_i$ on the element is given explicitly as $(b_i,c_i)$, so if we can get those somehow for each basis element, we can just compute the appropriate dot products and we are done. Furthermore, these gradients (and their dot products!) are constant, so if we compute them, we can just move them out in front of an integral over the triangle. In particular, we will multiply them by an estimate of $\int_{T_k}\kappa~\mathrm{d}x$ computed via some quadrature scheme. This can be done by converting back to a reference triangle and using Gaussian quadrature or something super simple like a midpoint rule.

2. The 3 basis elements are defined by $\phi_i(x_i,y_i) = 1$ and $\phi_i(x_j,y_j) = 0$, $j\neq i$. Since our basis functions are linear, this leads to a $3\times 3$ linear system for the coefficients of a single basis function, $\beta_i = (a_i,b_i,c_i)^\top$. For $i=1$, this system is given by $M\beta_1 = (1,0,0)^\top$. More generally, we have $M\beta_i = e_i$, where $e_i$ is the $i$-th canonical basis vector consisting of all zeros and a 1 in the $i$-th coordinate. However, notice that the $e_i$'s form the columns of the $3\times 3$ identity matrix, so all 3 linear systems can be written as a single matrix equation $$M\boldsymbol{\beta} = I \implies\boldsymbol{\beta} = M^{-1},$$ where $\boldsymbol{\beta}$ is a matrix with $\beta_i$ as its columns. Here we can see that the coefficients of the basis function are given explicitly by the columns of $\boldsymbol{\beta}$, and hence we can now evaluate our basis functions anywhere and we and we have our gradients and we can compute each dot product $$\nabla\phi_i\cdot\nabla\phi_j = b_ib_j+c_ic_j.$$

To put it all together in terms of the original steps for your example problem:

1. Compute via quadrature $\tilde{\kappa} = \int_{T_k}\kappa~\mathrm{d}x$. For example, consider the midpoint rule and evaluate $\kappa$ at the centroid and scale by the area: $$\tilde{\kappa} = \kappa(\bar{x},\bar{y})\cdot A = \kappa(0,-4/3)\cdot \frac{7}{2}.$$
2. Compute $M^{-1}$ to get all basis coefficients and get the gradients $$M = \begin{pmatrix}1 & 0 & 1 \\ 1 & -1 & -2 \\ 1 & 1 & -3 \\ \end{pmatrix}\implies M^{-1} = \frac{1}{7}\begin{pmatrix} 5 & 1 & 1 \\ 1 & -4 & -2 \\ 2 & -1 & -1 \\ \end{pmatrix}.$$
3. The coefficients $b_i$ and $c_i$ are in the bottom 2 rows of this matrix and the inner products are given by the elements of the matrix $$L = \left(\frac{1}{7}\right)^2 \begin{pmatrix}1 & 2\\ -4 & -1\\ -2 & -1\\ \end{pmatrix}\begin{pmatrix}1 & -4 & -2 \\ 2 & -1 & -1\end{pmatrix} = \frac{1}{49}\begin{pmatrix}5 & -6 & -4 \\ -6 & 17 & 9 \\ -4 & 9 & 5\end{pmatrix}.$$ The full stiffness matrix of this bilinear form can then be written using the fact that the gradients are constant, so we have $$\int_{T_k}\kappa \nabla\phi_i\cdot\nabla\phi_j~\mathrm{d}x = \nabla\phi_i\cdot\nabla\phi_j\int_{T_k}\kappa~\mathrm{d}x\approx \nabla\phi_i\cdot\nabla\phi_j \tilde{\kappa}.$$ We can now compute our stiffness matrix as $K = \tilde{\kappa}L$: $$K = \frac{\kappa(0,-4/3)}{14}\begin{pmatrix}5 & -6 & -4 \\ -6 & 17 & 9 \\ -4 & 9 & 5\end{pmatrix}.$$