Chapter 3
Computational Exercises

3.1 One dimensional integration in a standard segment

In this first exercise we will demonstrate how to integrate the function f(ξ) = ξ¹² on the standard segment ξ ∈ [-1,1] using Gaussian quadrature. The Gaussian quadrature weights and zeros are coded in the LIbUtilities library and for future reference this can be found under the directory $NEKDIST/library/LibUtilities/Foundations/. For the following exercises we will access the zero and points from the PointsManager. The PointsManager is a type of map (or manager) which requires a key defining known Gaussian quadrature types called PointsKey.

In the $NEKDIST/tutorial/fundamentals-integration/tutorial directory open the file named StdIntegration1D.cpp. Look over the comments supplied in the file which outline how to define the number of quadrature points to apply, the type of Gaussian quadrature and some arrays to hold the zeros, weights and solution. Finally a PointsKey is defined which is then used to obtain the zeros and weights in two arrays called quadZeros and quadWeights.

Task: 3.1

Implement a short block of code where you see the comments “Write your code here” which evaluates the loop

To compile your code type

make StdIntegation1D

in the tutorial directory. When your code compiles successfully¹ then type

./StdIntegration1D

You should now get some output similar to

====================================================== 
|        INTEGRATION ON A 1D STANDARD REGION         | 
====================================================== 
Integrate the function f(xi) = xi^12 on the standard 
segment xi=[-1,1] with Gaussian quadrature 
         Q = 4: Error = 0.179594

We can also use Gauss-Lobatto-Legendre type integration rather than Gauss-Legendre type in the previous exercises. To do this we replace

LibUtilities::PointsType quadPointsType = 
                        LibUtilities::eGaussGaussLegendre;

with

LibUtilities::PointsType quadPointsType = 
                        LibUtilities::eGaussLobattoLegendre;

3.2 Two-dimensional integration in a standard and local region

3.2.1 Quadrilateral element in a standard region

A straightforward extension of the one-dimensional Gaussian rule is to the two-dimensional standard quadrilateral region and similarly to the three-dimensional hexahedral region. Integration over Q² = {-1 ≤ ξ₁,ξ₂ ≤ 1} is mathematically defined as two one-dimensional integrals of the form

So if we replace the right-hand-side integrals with our one-dimensional Gaussian integration rules we obtain

where q₁ and q₂ are the number of quadrature points in the ξ₁ and ξ₂ directions. This expression will be exact if u(ξ₁,ξ₂) is a polynomial and q₁,q₂ are chosen appropriately. To numerically evaluate this expression the summation over ‘i’ must be performed q₁ times at every ξ_2i point, that is,

When executing the tutorial with the quadrature order q₁ = 6,q₂ = 7 you should get an output of the form:

============================================================= 
|        INTEGRATION ON 2 DIMENSIONAL ELEMENTS              | 
============================================================= 
 
Integrate the function f(x1,x2) = (x1)^12*(x2)^14 
on the standard quadrilateral element: 
         q1 = 6, q2 = 7: Error = 0.00178972

3.2.2 General straight-sided quadrilateral element

For elemental shapes with straight sides a simple mapping may be constructed using a linear mapping similar to the vertex modes of a hierarchical/modal expansion. For the straight-sided quadrilateral with vertices labeled as shown in figure 3.1(b) the mapping can be defined as:

If we denote an arbitrary quadrilateral region by Ω^e which is a function of the global Cartesian coordinate system (x₁,x₂) in two-dimensions. To integrate over Ω^e we transform this region into the standard region Ω_st defined in terms of (ξ₁,ξ₂) and we have

where J_2D is the two-dimensional Jacobian due to the transformation, defined as:

(3.2)

As we have assumed that we know the form of the mapping [i.e., x₁ = χ₁(ξ₁,ξ₂), x₂ = χ₂(ξ₁,ξ₂)] we can evaluate all the partial derivatives required to determine the Jacobian. If the elemental region is straight-sided then we have seen that a mapping from (x₁,x₂) → (ξ₁,ξ₂) is given by equations (3.1).

Task: 3.5 We now consider how to integrate the function f(x₁,x₂) = x₁¹² x₂¹⁴ on a local rectangular quadrilateral element using Gaussian quadrature. Consider the local quadrilateral element with vertices

This is clearly similar to the previous exercise. However, as we are calculating the integral of a function defined on a local element rather than on a reference element, we have to take into account the geometry of the element. Therefore, the implementation is altered in two ways:

1. The quadrature zeros should be transformed to local coordinates to evaluate the integrand f(x₁,x₂) at the quadrature points.

2. The Jacobian of the transformation between local and reference coordinates should be taken into account when evaluating the integral.

In the file LocIntegration2D.cpp you are provided with the same set up as the previous task but now with a definition of the coordinate mapping included. Evaluate the expression for the Jacobian analytically. Then write a line of code in the loop for the Jacobian as indicated by the comments “Write your code here”. When you have written your expression you can compile the code with the command

make LocIntegration2D

Verify that the error is not equal to zero when q₁ = 8,q₂ = 9. Why might this be the case?³ .

Using the quadrature order specified in the file your output should look like:

=========================================================== 
|      INTEGRATION ON 2D ELEMENT in Local Region          | 
=========================================================== 
 
Integrate the function f(x1,x2) = x1^12 * x2^14 
on a local quadrilateral element: 
         Error = 0.000424657