Chapter 3
Computational Exercises

3.1 One dimensional differentiation in a standard segment

In this first exercise we will demonstrate how to differentiate the function f(ξ) = ξ⁷ in the standard segment ξ ∈ [-1,1] using Gaussian quadrature. The Gaussian quadrature zeros and the differentiation matrix 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 zeros and differentiation matrix 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-differentiation/tutorial directory open the file named StdDifferentiation1D.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, a differentiation matrix and some arrays to hold the zeros and solution. Finally a PointsKey is defined which is then used to obtain the zeros in an array called quadZeros and the differentiation matrix in a pointer to a matrix called derivMatrix.

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 StdDifferentiation1D

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

./StdDifferentiation1D

You should now get some output similar to

====================================================== 
|      DIFFERENTIATION IN A 1D STANDARD REGION       | 
====================================================== 
Differentiate the function f(xi) = xi^7 in the 
standard segment xi=[-1,1] using quadrature points 
         Q = 7: Error = 1.49647

We can also use Gauss-Lobatto-Legendre type differentiation 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 differentiation 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. Differentiation in Q² = {-1 ≤ ξ₁,ξ₂ ≤ 1} in the ξ₁ direction is defined as

where h₁ and h₂ are the polynomials associated respectively with coordinates ξ₁ and ξ₂.

Because h₂ is not a function of ξ₁, we can rewrite the formula as

(3.1)

Typically, we only require the derivative at the nodal points (ξ_1i,ξ_2j). At these points, h₂ has trivial values, i.e. unit value at ξ_2j and null value at all other points. The summation over j therefore drops and we are left with the same formula as in our one-dimensional problem

(3.2)

where d_1ik is the differentiation matrix in the ξ₁ direction such that

Likewise, the derivative with respect to ξ₂ can be expressed by

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

=========================================================== 
|    DIFFERENTIATION IN 2D ELEMENT in Standard Region     | 
=========================================================== 
 
Differentiate the function f(x1,x2) = (x1)^7*(x2)^9 
on the standard quadrilateral element: 
         q1 = 7, q2 = 9: Error = 7.19196

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:

We denote an arbitrary quadrilateral region by Ω^e which is a function of the global Cartesian coordinate system (x₁,x₂) in two-dimensions. To differentiate in Ω^e we transform this region into the standard region Ω_st defined in terms of (ξ₁,ξ₂). We begin with the basic definition of differentiation:

Unlike differentiation in the standard region, both h₁ and h₂ are functions of both local coordinates. The chain rule can be applied to obtain a system similar to the previous exercise:

(3.4)

where h₁ and h₂ are functions of ξ₁ and ξ₂ respectively only and where and come from the inverse two-dimensional Jacobian matrix due to the transformation, defined as:

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 matrix. If the elemental region is straight-sided then we have seen that a mapping from (x₁,x₂) → (ξ₁,ξ₂) is given by equations (3.3).

Equation (3.4) turns out to be similar to equation (3.1) in the standard region and can be also rewritten in the style of equation (3.2):

Similarly to the differentiation in the standard region, derivatives with respect to the other local coordinate can also be found by substituting dx₂ for dx₁. We can also note that equation (3.2) is a particular case of equation (3.5) where = 1 and = 0. This effectively corresponds to the situation where the local region has the same shape and dimensions as the standard region (with the admission of a translation in the plane).

Task: 3.4 We now consider how to differentiate the function f(x₁,x₂) = x₁⁷ x₂⁹ in a local rectangular quadrilateral element. Consider the local quadrilateral element with vertices

This is clearly similar to the previous exercise. However, as we are calculating the derivatives of a function defined in a local element rather than in 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 function f(x₁,x₂) at the quadrature points.

2. Elements of the inverse Jacobian matrix of the transformation between local and reference coordinates should be taken into account when evaluating the derivatives.

In the file LocDifferentiation2D.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 matrix analytically and find its inverse. 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 LocDifferentiation2D

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

========================================================= 
|     DIFFERENTIATION IN 2D ELEMENT in Local Region     | 
========================================================= 
 
Differentiate the function f(x1,x2) = x1^7 * x2^9 
in a local quadrilateral element: 
         q1 = 8, q2 = 10: Average Error = 0.0346594

Advanced Task: 3.5

As it turns out in the previous task, the average error is not equal to zero. Why is that?

Try different values of q₁ and q₂ and plot the average error with respect to these two parameters, either one by one or simultaneously. Why are the values of q_1max and q_2max (at which the average error reaches computer precision) different from those expected in the standard region?³