Numerical Integration

To form the finite element integrals we need to be able to evaluate integrals of the form $\int_{\Omega} \phi_i \phi_j dx$, $\int_{\Omega} \nabla \phi_i \cdot \nabla \phi_j dx$ and $\int_{\Omega} f ~ \phi_i~ dx$.

As we will see later, we decompose the domain into a set of elements $\{T\} =: \mathcal{T}_h$ and then decompose the integration problem into the element problems.

For now, we will assume that the domain of integration is a single element $T$.

This element furthermore can be mapped to a reference element of same type (segment, triangle, quadrilateral, tetrahedron, hexahedron) so that the fundamental integration problem to solve is the integration over the reference element.

On the reference interval the previously mentioned integrals could be solved analytically. However, for coefficient functions that are not polynomials, e.g. in the r.h.s. source terms or the diffusivity, this is soon not possible. Therefore, most often numerical integration is used in finite element codes.

1D Numerical integration

As Numerical integration has been considered in the numerics lecture (Numerik II), we will not do a detailed lecture on this topic here. Instead, we will just briefly recap the main ideas and then go to the implement of some standard 1D numerical integration schemes.

The integration problem:

$$ I(f) := \int_{0}^{1} f(x) dx \approx Q(f) := \sum_{i=1}^{n} w_i f(x_i) \tag{Q} $$

Let's do a brief recap of numerical integration based on a few repition questions (if needed we will do a more detailed lecture/session on this):

Interpolation quadrature:

• What is the idea of interpolation quadrature?

Interpolation Exactness:

• What is the exactness degree of a quadrature scheme?
• How can it be used constructively to come up with a quadrature scheme of a certain exactness degree?

Gaussian quadrature:

• What is the idea of Gaussian quadrature?
• Which accuracy can be achieved with Gaussian quadrature?
• How are Gaussian quadrature points and weights computed? How are they related to orthogonal polynomials?

Implementation: $\leadsto$ intrule.py, intrule_1d.py

Tasks:

Task Int1D-1:

Complete the implementation of the Newton-Cotes Rules of arbitrary order $n$ (see intrule_1d.py) and make sure that the unit test test_newtoncotes in tests/test_intrule_1d.py passes.

Task Int1D-2:

Implement the Gaussian-Legendre quadrature rule of order $n$ and make sure that the corresponding unit tests tests/test_intrule_1d.py passes. Add your own unit tests, e.g. to check for positivity of the weights.

To compute the nodes and weights you can choose different approaches:

• formulate a nonlinear problem based on the exactness property
• compute the zeros of the corresponding orthogonal polynomials, e.g. by a combination of a bisection and Newton's method

To test your implementation you can compare to the numpy-implementation of the Gauss-Legendre quadrature rule which is wrapped into the NP_GaussLegendreRule in intrule_1d.py.

Finally, the integration rule should be implemented w.r.t. the interval $(0,1)$.

Task Int1D-3:

For integration tasks of the form $\int_{-1}^{1} (1-x)^\alpha (1+x)^\beta f(x) dx$, i.e. the quadrature problem with the positive weight function $w(x) = (1-x)^\alpha (1+x)^\beta$, the Gauss-Rules are called Gauss-Jacobi rules and the corresponding orthogonal polynomials are the Jacobi polynomials that we saw before in fe.ipynb.

Hint: Note that with $x = 2\tilde x-1$ we have

• $dx = 2 d\tilde x$ and
• $\tilde x = \frac{1}{2} (x + 1)$ and
• $1+x = 2 \tilde x$ and $1-x = 2 (1-\tilde x)$ and hence
• $(1-x)^\alpha (1+x)^\beta = 2^{\alpha+\beta} (1-\tilde x)^\alpha \tilde x^\beta$ which gives
• $\int_{-1}^{1} (1-x)^\alpha (1+x)^\beta f(x) dx = 2^{\alpha+\beta+1} \int_{0}^{1} (1-\tilde x)^\alpha \tilde x^\beta f(2 \tilde x -1 ) d\tilde x$

Now, implement the Gauss-Jacobi quadrature rule of order $n$ and make sure that the corresponding unit test (for $\alpha=\beta=0$) in tests/test_intrule_1d.py passes. Add your own unit tests (for $(\alpha,\beta)\neq (0,0)$). Again, you may compare to existing implementations from scipy which is wrapped into the SP_GaussJacobiRule in intrule_1d.py.

For the implementation you can use the same approaches as for the Gauss-Legendre rule above. Finally, the integration rule should be implemented w.r.t. the interval $(0,1)$.

These integration rules will be important for arbitrary order integration rules on triangles later.

2D Numerical integration

For low and moderate orders of the integration problem we can proceed as for the 1D case while for higher orders we typically want to resort to tensor product rules at some point. These are easily set up for quadrilaterals and hexahedrons but need some tricks to be applied on triangles and tetrahedrons.

In the context of this course we will concencrate on triangles in 2D and use the following triangle as the reference triangle:

$\hat{T} = \operatorname{conv}\{(0,0),(1,0),(0,1)\} = \{ (x,y) \in \mathbb{R}^2 ~|~ x\geq 0, y\geq 0, x+y\leq 1 \}$

Implementation: intrule_2d.py

Low order integration rules

Based on exactness demands and counting arguments (dimension of the space of polynomials) and a selection of points, weights can be computed.

Task Int2D-1:

Implement the a vertex integration rules for triangles, i.e. find the weights $w_i$ such that the integration rule $Q(f) = \sum_{i=1}^{3} w_i f(x_i,y_i)$ is exact for all polynomials of degree $\leq 1$ where $(x_0,y_0) = (0,0)$, $(x_1,y_1) = (1,0)$ and, $(x_2,y_2) = (0,1)$.

Here, you may want to compute the weights on paper and simply hack them into the code or use a generic approach to compute the weights (which is more involved, but may help for the proceeding tasks).

Check the exactness degree and implement a test for the integration rule in a new test file tests/test_intrule_2d.py.

Task Int2D-2:

Combining the edge midpoints and the vertices, we have $6$ points in total. Find the weights $w_i$ such that the integration rule is exact for all polynomials of degree $\leq 2$. Implement the rule and test it.

Task Int2D-3:

For a Gauss-Legendre rule with 3 nodes we have $6$ degrees of freedom (3 nodes and 3 weights). Therefore, we can expect to be able to integrate polynomials of degree $\leq 2$ exactly.

Finding the nodes and weights is a nonlinear problem and will most certainly require a numerical approach. Write a proper nonlinear solver for the problem and implement the integration rule. Afterwards test it.

High order integration rules and the Duffy transformation

To generate generic higher order integration rules on triangles, one tool is the Duffy transformation which maps between the reference triangle and the reference square:

The idea is then to transform the integral on the triangle to an integral on the square and then use a tensor product rule on the square (i.e. iterated integrals).

With $\Phi : \hat Q \to \hat T, \Phi(x_1,x_2) = ((1-x_2)x_1 ,x_2)$ (and hence $\Phi^{-1}(x_1,x_2) = (\frac{x_1}{1-x_2},x_2)$) being the Duffy transformation, there holds $$ I_{\hat{T}}(f) := \int_{\hat{T}} f(x_1,x_2) d(x_1,x_2) = \int_{\hat{Q}} |\operatorname{det}(D \Phi)| (f\circ \Phi)(x_1,x_2) d(x_1,x_2) $$ with $|\operatorname{det}(D \Phi)| = (1-x_2)$. Let $\tilde f = f \circ \Phi$, then the new integration problem is $$ \int_{\hat{Q}} (1-x_2) \tilde f(x_1,x_2) d(x_1,x_2) = \int_0^1 \underbrace{\int_0^1 (1-x_2) \tilde f(x_1,x_2) ~ d x_2}_{I^{(1,0)}(\tilde f(x_1,\cdot))} ~ d x_1 = I^{(0,0)}(\tilde I(x)) \text{ with } \tilde I(x) := I^{(1,0)}(\tilde f(x,\cdot)) $$ Now, we can replace the 1D-Integrals with quadrature rules where we replace the inner integral $I^{(1,0)}$ with a 1D quadrature rule $Q^{(1,0)}$ and the outer integral $I^{(0,0)}$ with a 1D quadrature rule $Q^{(0,0)}$: $$ \int_{\hat{T}} f(x_1,x_2) d(x_1,x_2) \approx Q^{(0,0)}( \tilde Q(x)) \text{ with } \tilde Q(x) := Q^{(1,0)}(\tilde f(x,\cdot)) $$ This yields $$ \int_{\hat{T}} f(x_1,x_2) d(x_1,x_2) \approx \sum_i w_i^{(0,0)} \sum_j w_j^{(1,0)} \tilde f((x_1)_i, (x_2)_j) = \sum_i w_i^{(0,0)} \sum_j w_j^{(1,0)} f((1-(x_2)_j) (x_1)_i, (x_2)_j) \tag{D} $$

Task Int1D-4:

Assume that $Q^{(0,0)}$ is a Gauss-Legendre rule with exactness degree $n$, i.e. $Q^{(0,0)}(f) = I^{(0,0)}(f)$ holds for all polynomials $f$ of degree $\leq n$ and $Q^{(1,0)}$ is a Gauss-Jacobi rule with exactness degree $n$, i.e. $Q^{(1,0)}(f) = I^{(1,0)}(f)$ holds for all polynomials $f$ of degree $\leq n$. Here, $I^{(\alpha,\beta)}(f) = \int_0^1 (1-x)^\alpha x^\beta f(x) dx$. Show that the integration rule (D) is exact for all polynomials of degree $\leq n$.

Task Int1D-5:

Based on the previous observations implement an arbitrary order accurate integration rule for the triangle and test it. Use 1D Gauss-Legendre and Gauss-Jacobi rules and corresponding transformation formulas to compute and make use of the 1D integral rules.