Computational Exercises

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

`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]`

.
To compile your code type

`make StdDifferentiation1D`

in the tutorial directory. When your code compiles
successfully^{1}
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

| 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

with

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^{2} = {-1 ≤ ξ_{1},ξ_{2} ≤ 1} in the ξ_{1} direction is defined as

where h_{1} and h_{2} are the polynomials associated respectively with coordinates ξ_{1} and
ξ_{2}.

Because h_{2} is not a function of ξ_{1}, we can rewrite the formula as

| (3.1) |

Typically, we only require the derivative at the nodal points (ξ_{1}_{i},ξ_{2}_{j}). At these points, h_{2} has
trivial values, i.e. unit value at ξ_{2}_{j} 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_{1}_{i}k is the differentiation matrix in the ξ_{1} direction such that

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

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 q_{1} = 7,q_{2} = 9 using
a Gauss-Lobatto-Legendre quadrature rule. Complete the code by writing a
structure of loops which implement the two-dimensional Gaussian quadrature
rule^{2} .
The expected output is given below). Also verify that the error is zero when q_{1} = 8,q_{2} = 10.

Recall that to compile the file you type

`make StdDifferentiation2D`

When executing the tutorial with the quadrature order q_{1} = 7,q_{2} = 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

| 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

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_{1},x_{2}) in two-dimensions. To differentiate in Ω^{e} we transform this region
into the standard region Ω_{st} defined in terms of (ξ_{1},ξ_{2}). We begin with the basic definition of
differentiation:

Unlike differentiation in the standard region, both h_{1} and h_{2} 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_{1} and h_{2} 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., x_{1} = χ_{1}(ξ_{1},ξ_{2}),
x_{2} = χ_{2}(ξ_{1},ξ_{2})] 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_{1},x_{2}) → (ξ_{1},ξ_{2}) 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_{2} for dx_{1}. 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).

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:

- The quadrature zeros should be transformed to local coordinates to evaluate the
function f(x
_{1},x_{2}) at the quadrature points. - 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

| 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_{1} and q_{2} and plot the average error with respect to these two parameters,
either one by one or simultaneously. Why are the values of q_{1}_{max} and q_{2}_{max} (at which the
average error reaches computer precision) different from those expected in the standard
region?^{3}