## Chapter 3Computational Exercises

### 3.1 One dimensional differentiation in a standard segment

In this first exercise we will demonstrate how to differentiate the function f(ξ) = ξ7 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`.

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

Tip: To access individual elements in `MatrixSharedPtrType`, the pointer must first be dereferenced then accessed using the parantheses operator. This should look like `(*derivMatrix)(i, j)`. On the other hand, to access elements in `Array<OneD, NekDouble>`, the bracket operator can be directly used such as `quadZerosDir1[i]`. If the array is two-dimensional (`Array<TwoD, NekDouble>`), elements would be accessed through `quadDerivsDir1[i][j]`.

`make StdDifferentiation1D`

in the tutorial directory. When your code compiles successfully1 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

Task: 3.2 Evaluate the previous derivatives for a quadrature order of Q = Qmax where Qmax = 8 is the number of quadrature points required for an exact evaluation of the derivatives (calculate this value analytically). Verify that the error is zero (up to numerical precision).

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

LibUtilities::eGaussGaussLegendre;

with

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 Q2 = {-1 ≤ ξ12 ≤ 1} in the ξ1 direction is defined as

where h1 and h2 are the polynomials associated respectively with coordinates ξ1 and ξ2.

Because h2 is not a function of ξ1, we can rewrite the formula as

 (3.1)

Typically, we only require the derivative at the nodal points (ξ1i2j). At these points, h2 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 d1ik is the differentiation matrix in the ξ1 direction such that

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

Task: 3.3 Differentiate the function f(ξ12) = ξ17ξ29 in the standard quadrilateral element Q ∈ [-1,1] × [-1,1] using Gaussian quadrature points.

Using a series of one-dimensional Gaussian quadrature rules as outlined above evaluate the derivative in each direction and for each quadrature point by completing the first part of the code in the file `StdDifferentiation2D.cpp` in the directory

`\$NekDist/tutorial/fundamentals-differentiation/tutorial`.

The quadrature zeros and differentiation matrices in each of the coordinate directions have already been setup and are initially set to q1 = 7,q2 = 9 using a Gauss-Lobatto-Legendre quadrature rule. Complete the code by writing a structure of loops which implement the two-dimensional Gaussian quadrature rule2 . The expected output is given below). Also verify that the error is zero when q1 = 8,q2 = 10.

Recall that to compile the file you type

`make StdDifferentiation2D`

When executing the tutorial with the quadrature order q1 = 7,q2 = 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
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 (x1,x2) in two-dimensions. To differentiate in Ωe we transform this region into the standard region Ωst defined in terms of (ξ12). We begin with the basic definition of differentiation:

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

 (3.4)

where h1 and h2 are functions of ξ1 and ξ2 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., x1 = χ112), x2 = χ212)] 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 (x1,x2) → (ξ12) 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 dx2 for dx1. 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(x1,x2) = x17x29 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(x1,x2) 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