Lab Manual

User Manual:

Open the PDF directly: View PDF .
Page Count: 23

Download
Open PDF In Browser	View PDF

MEAN-DG
Mixed Elements Applications with Nodal Discontinuous Galerkin
methods

Lab Manual

by

Subodh Joshi, Tapan Mankodi, Prof. S. Gopalakrishnan

Indian Institute of Technology, Bombay
Mumbai
2018

Contents
1 Introduction
1.1 What . . . . . .
1.2 Why . . . . . .
1.3 Who . . . . . .
1.4 How to refer this

. . . . . . .
. . . . . . .
. . . . . . .
manual . . .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

2 Mathematical Background
2.1 Introduction . . . . . . . . . . . . . . .
2.2 Mapping on the Reference Element . .
2.2.1 Hexahedral cells . . . . . . . . .
2.3 Mapping onto the Reference Face . . .
2.3.1 Quadrilateral Faces . . . . . . .
2.4 Modal and Nodal Basis Functions . . .
2.4.1 Modal Basis Functions . . . . .
2.4.2 Nodal Basis Functions . . . . .
2.4.3 Multidimensional Interpolation
3 Code Structure
3.1 Code structure . . . . . . . .
3.1.1 Tensor . . . . . . . . .
3.1.2 Point . . . . . . . . . .
3.1.3 Cell . . . . . . . . . .
3.1.4 Face . . . . . . . . . .
3.1.5 GeometryIO . . . . . .
3.1.6 Functional Space . . .
3.1.7 Math functions . . . .
3.1.8 Gasdynamics, Riemann
3.1.9 DG . . . . . . . . . . .

.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.

. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
solvers etc.
. . . . . . .

4 Applications

.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.

2
2
2
3
3

.
.
.
.
.
.
.
.
.

4
4
5
5
8
8
12
12
13
14

.
.
.
.
.
.
.
.
.
.

15
16
16
16
16
17
17
17
17
18
18
19

5 License
20
5.1 LICENSE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

i

List of Figures
2.1
2.2
2.3

a. Reference cell and faces relative numbering b.
vertices numbering . . . . . . . . . . . . . . . . . . .
First 9 Legendre polynomials . . . . . . . . . . . . .
First 9 Chebyshev polynomials . . . . . . . . . . . . .

3.1

Class dependency (crudely) . . . . . . . . . . . . . . . . . . . . . . . . . . 15

ii

Reference
. . . . . .
. . . . . .
. . . . . .

cell and
. . . . . . 8
. . . . . . 12
. . . . . . 13

List of Tables

1

Chapter 1
Introduction
1.1

What

MEAN-DG (Mixed Elements Applications with Nodal Discontinuous Galerkin method)
aims to be a versatile, industry-applications oriented efficient Discontinuous Galerkin
code. However, all efforts have been taken to keep it as a flexible framework such that
in addition to the Galerkin family of schemes, other schemes such as the FV scheme,
spectral methods etc. can be written as required. The code in its finished form is expected to cater four types of 3D elements: hexahedral, tetrahedral, prism and pyramid
cells. It can be customized for a wide array of equations ranging from as simple as the
linear advection equation, Laplace equation to as complicated as the Navier Stokes (N-S)
equations and Magnetohydrodynamics (MHD) equations. The code is mostly a c++ code
with some utilities written in Python language. The data-structure and the file reading
format is inspired from the OpenFOAM file format. Thus, test-cases and applications for
OpenFOAM can be readily solved using MEAN-DG (depending on the solvers coded).
After the current phase the code will have following features:
• Solve advection dominated flows efficiently
• Achieve arbitrary higher-order accuracy in space
• Open-MP parallelization
• Capability to handle flexible geometries and different mesh element types
• Capability to code other systems of equations

1.2

Why

Higher-order accuracy in space is extremely important for many applications such as
boundary layer flows, vortex dominated flows, linear propagating waves (such as in electromagnetics and aeroacoustics) etc. Finite volume (FV) methods are well established
and mature enough to cater to the industry demands. However, FV methods suffer from
certain limitations, such as:
• Poor hardware scalability and parallelization capacity (because of wider stencils
used and the inter-cell communications it demands)
2

• Limitations on order of accuracy because of non-compact nature of the reconstruction
• Difficult boundary treatments
On the other hand the spectral methods and global collocation methods are extremely
accurate but lack flexibility and simplicity of the finite volume methods. The DG methods
fall right in the sweet spot between the element based Galerkin methods and the finite
volume methods. As of now, not many large-scale codes are available based on the DG
framework. A few codes are available such as,
1. Deal.II (http://www.dealii.org/)
2. Dune (https://www.dune-project.org/about/dune/)
3. Fenix (https://fenicsproject.org/documentation/)
The idea behind MEAN-DG is to have a DG code similar to the OpenFOAM FV code
for industry level applications. It also serves the purpose of capacity building, indigenous
technology building and have an in-house computation-platform.

1.3

Who

The PI of the project is Prof. S. Gopalakrishnan. Primary contributions to the project
are made by Tapan Mankodi and Subodh Joshi. The project is funded by Siemens Technologies, Bangalore.

1.4

How

to refer this manual

Chapter 2 covers the necessary mathematical background to get started with the code.
Chapter 3 explains the structure of the code. Chapter 4 explains the work-flow for a few
applications solved with MEAN-DG. Chapter 5 consists of the license. User is strongly
recommended to go through it before using the code.

3

Chapter 2
Mathematical Background
2.1

Introduction

Discontinuous Galerkin finite element method is an amalgamation of the element based
Galerkin schemes and the finite volume method. It is a compact, higher-order accurate
scheme where the solution can be discontinuous between two elements. The DG scheme
requires the physical domain (Ω) discretized with N e finite elements (written as Ωk ). The
finite element cells are assumed to be one of the following types.
1. Hexahedral (8 edges, 8 nodes, 6 faces)
2. Tetrahedral (6 edges, 4 nodes, 4 faces)
3. Prism (9 edges, 6 nodes, 5 faces)
4. Pyramid (8 edges, 5 nodes, 5 faces)
The solution is assumed to be continuous (along with continuous derivatives depending on
the governing equations) inside the element. However, the solution can be discontinuous
at the edges and faces. The value of the conserved variable (and the numerical fluxes) at
the faces and edges are found out by solving the Riemann problem which arises because
of the discontinuity.
The DG method essentially starts with the ‘variational’ formulation of the problem.
Consider a differential equation given by,
L(u) = f (u), u ∈ H 1 (Ω)

(2.1)

d
and
(and suitable boundary conditions, not mentioned here). For example, L = dx
f (u) = u results in,
du
=u
(2.2)
dx
We will be using the general notation given by equation 2.1 henceforth. Then, if uk is the
approximate solution inside of the element, then instead of the solving equation 2.1 with
uk , we solve the following variational formulation of equation 2.1.
Z
Z
k
φL(u )dΩ =
φf (uk )dΩ
(2.3)
Ωk

Ωk

4

φ is the test function. The approximate solution uk is obtained by polynomial approximation of u inside Ωk . A modal interpolation (basis functions consisting of modes, i.e.
frequencies) or a nodal interpolation (basis functions consist of the cardinal polynomials
such as the Lagrange polynomials) can be used for this. In the DG methods, the basis
functions serve as the test functions. The integrations in equation 2.2 are solved numerically using quadrature formulas. The resulting algebraic system is solved iteratively (or
ODE is solved numerically if time dependant). This process is repeated till convergence
(for steady state problems) or till the final time is reached (for transient problems). It
is easier to solve the problem on a reference cell (canonical element) and then map the
solution on to the physical cell. This requires a mapping to be defined between the reference cell and the physical cell. Based on this general process, the following mathematical
concepts are revised in this chapter.
1. Mapping functions (Jacobian (J) of transformation, inverse Jacobian (J−1 ) and
(|J|))
2. Modal and nodal basis functions
3. Polynomial interpolation and differentiation
4. Numerical quadrature (integration)
5. Cell matrices (Vandermonde matrix, Mass matrix, differentiation matrix and flux
matrix)
6. Other mathematical functions
The discussion is kept brief and concise and only those concepts concerned with MEANDG are discussed. The user is encouraged to go through the references for in-depth theory
of DG schemes.
Note 2.1. We use the c++ style indexing. i.e. indices start at 0 (not 1).

2.2

Mapping on the Reference Element

Corresponding to each physical cell, there exists a canonical (reference) element. The
physical cell is expressed in terms of the physical coordinate system (x, y, z). The reference
coordinate system is written in terms of (r, s, t). Thus, each of the physical coordinates
can be written in terms of the (r, s, t) coordinate system. In the code, we describe the
physical coordinate system as XY Z, and the reference system as RST . The mapping from
the physical coordinate system to the reference coordinate system is characterized by the
Jacobian of transformation J. In this section, we describe the Jacobian of transformation
for the four types of the cells.

2.2.1

Hexahedral cells

Hexahedral, or Hex in short, cells are mapped to a cube in 3D. The Hex cell has 8 nodes
given by the array, ((x0 , y0 , z0 ), (x1 , y1 , z1 ), ..., (x7 , y7 , z7 )). We define the linear shape
functions on the reference element as follows:
1
N0 = (1 − r)(1 − s)(1 − t),
8

{r, s, t ∈ V |∀v ∈ V ⊂ R, v ∈ [−1, 1]}
5

(2.4)

1
N1 = (1 + r)(1 − s)(1 − t)
8
1
N2 = (1 + r)(1 + s)(1 − t)
8
1
N3 = (1 − r)(1 + s)(1 − t)
8
1
N4 = (1 − r)(1 − s)(1 + t)
8
1
N5 = (1 + r)(1 − s)(1 + t)
8
1
N6 = (1 + r)(1 + s)(1 + t)
8
1
N7 = (1 − r)(1 + s)(1 + t)
8
The mapping then is simply given as,
x(r, s, t) =

7
X

(2.5)
(2.6)
(2.7)
(2.8)
(2.9)
(2.10)
(2.11)

Ni xi

(2.12)

Ni yi

(2.13)

Ni zi

(2.14)

i=0

y(r, s, t) =

7
X
0

z(r, s, t) =

7
X
0

Definition 1. The Jacobian of transformation (J) is defined as :

 




J=




∂x
∂r

∂x
∂s

∂x
∂t

∂y
∂r

∂y
∂s

∂y
∂t

∂z
∂r

∂z
∂s

∂z
∂t

  J00 J01 J02
 
 
= J
  10 J11 J12
 
 
J20 J21 J22









(2.15)

The partial derivatives in equation 2.15 can be easily computed from equations 2.4-2.11
and 2.12, 2.13, 2.14. In the code, the above implementation is found in Cell.cpp.
Cell::calculateJacobian3DTensor(double r,double s,double t)
Note 2.2. The first name indicates the name of the class (Cell), the second symbol (::)
indicates the scope (in this case, indicating ’belongs to Cell’ the third name is the name of
the class method (calculateJacobian3DTensor) and terms in the bracket being arguments
passed to this function.
It is to be noted that the name ‘Tensor’ in the above function indicates that the functional
space for a Hex cell is obtained by a tensor product (outer product) of 1D functional
spaces.
The inverse of the Jacobian is found out by the standard process of inversion of a
matrix.
6


J−1




=





∂r
∂x

∂r
∂y

∂r
∂z

∂s
∂x

∂s
∂y

∂s
∂z

∂t
∂x

∂t
∂y

∂t
∂z









(2.16)

These are given as:
∂r
1
=
(J22 J11 − J21 J12 )
∂x
|J|

(2.17)

∂r
−1
=
(J22 J01 − J21 J02 )
∂y
|J|

(2.18)

1
∂r
=
(J12 J01 − J11 J02 )
∂z
|J|

(2.19)

∂s
−1
=
(J22 J10 − J20 J12 )
∂x
|J|

(2.20)

∂s
1
=
(J22 J00 − J20 J02 )
∂y
|J|

(2.21)

−1
∂s
=
(J12 J00 − J10 J02 )
∂z
|J|

(2.22)

1
∂t
=
(J21 J10 − J20 J11 )
∂x
|J|

(2.23)

∂t
−1
=
(J21 J00 − J20 J01 )
∂y
|J|

(2.24)

1
∂t
=
(J11 J00 − J10 J01 )
∂z
|J|

(2.25)

where, the determinant of the Jacobian |J| is calculated as:
|J| = J00 (J22 J11 − J21 J12 ) − J01 (J22 J10 − J20 J12 ) + J02 (J21 J10 − J20 J11 )

(2.26)

The inverse of Jacobian is implemented in
Cell::calculateInverseJacobianMatrix3DTensor(arguments)
Note 2.3. The inverse of Jacobian is called as dRST_by_dXYZ in the code. Each
cell stores these matrices for all the quadrature points within that cell. The metric terms
given by 2.17–2.25 are used in computations involving derivatives as seen in the subsequent
sections and chapters.
The inverse Jacobian is stored for all the quadrature points for all the cells. It is an
attribute of the cell (i.e. Cell::dRST_by_dXYZ) and has dimensions nQuad × 3 × 3.
Use of the Jacobian Consider the following integration:
Z
I=
φ(x, y, z)dΩk
(2.27)
Ωk

7

The integration can be written as,
Z
I=

φ(r, s, t)|J|dΩs

(2.28)

Ωs

where, Ωs is a reference element and J is the Jacobian matrix. The main benefit of
performing the above conversion is that, the function φ(r, s, t) needs to be computed only
once for the standard element, each cell will have its own Jacobian matrix, so number of
computations reduce drastically.
In the current code, the vertices and faces are numbered in a standardized fashion as
shown in Fig.2.1.

5

1 Top
2

3 LHS

F
ro

4
4

n

t

6

7

Ba

5 RHS

ck

1

2

0 Base
0
Figure 2.1: a. Reference cell and faces relative numbering
numbering

3

b. Reference cell and vertices

The corresponding implementation can be found in:
Cell::orderDefiningPoints()

2.3

Mapping onto the Reference Face

DG method involves evaluation of surface integrations in order to find the numerical flux
at the cell-interface. The surface integrations in turn require mapping of a given arbitrary
surface on the reference surface. The parameterization of a plane surface requires two
parameters r and s.

2.3.1

Quadrilateral Faces

For a quadrilateral faces, the parameters assume the following values: (r, s) ∈ [−1, 1].
However, the given surface can be arbitrarily oriented in the space, and thus requires
(x, y, z) position to be completely defined. We thus need to orient the surface such as one
of the variables becomes redundant. This means rotating the surface such as to make it
parallel to xy plane, for example. Thus, each surface needs to be re-oriented in space such
that the face normal n̂ becomes parallel to the z axis. After the we can map the rotated
face onto the standard element. Let θz be the angle made by the face normal with the z
8

axis and θx be the angle made by the face normal with the x axis. The rotation matrix
rotates the face through θx first (counterclockwise) and the through θz counterclockwise.
The resulting rotation matrix R is given as:
R =< Rz , Rx >

(2.29)

where,


cos(θx ) sin(θx ) 0


Rx =  −sin(θx ) cos(θx ) 0

0
0
1







;




Rz = 


cos(θz ) 0 −sin(θz )
0

1

sin(θz ) 0

0
cos(θz )




;


(2.30)

and < a, b > indicates the dot product between matrices a and b. Matrix R−1 indicates
the reverse transformation. The rotation matrix is calculated in,
Face::calculateRotationMatrixParallelToXY()
and the inverse matrix is computed in,
Face::calculateInverseRotationMatrix()
In order to describe the procedure to compute the Jacobian of transformation (|Jf |) for
face, we use the following nomencleture:
• (x, y, z): physical coordinate system (actual coordinates)
• (x0 , y 0 , z 0 ): Rotated coordinates such that the face is || to the xy plane
• (r, s): Reference coordinates such that r, s ∈ [−1, 1]
There are two ways the Jacobian can be computed.
Method 1
First, rotate each vertex of the face using R. i.e.




0
x
x
 i 
 i 
 0 


 yi  = R  yi  ,




0
zi
zi

i ∈ {0, 1, 2, 3}

(2.31)

Each of the rotated points is then mapped to the standard element. After rotation of the
0
0
= 0, ∂z
= 0 since z 0 remains constant. The mapping from the rotated face f 0 to
face, ∂z
∂r
∂s
the standard face f s can be defined through the shape functions.
1
N0 = (1 − r)(1 − s),
4

{r, s ∈ V |∀v ∈ V ⊂ R, v ∈ [−1, 1]}

1
N1 = (1 + r)(1 − s)
4
9

(2.32)
(2.33)

1
N2 = (1 + r)(1 + s)
4
1
N3 = (1 − r)(1 + s)
4

(2.34)
(2.35)

Then,
0

x =

3
X

Ni x0i ,

0

y =

3
X

i=0

Ni yi0 ,

(2.36)

i=0

Equation2.36 can be used to find partial derivatives of x0 and y 0 with respect to r and s.
The required face Jacobian is:




f
J =


∂x
∂r

∂x
∂s

∂y
∂r

∂y
∂s

∂z
∂r

∂z
∂s






(2.37)

which can be written as,




x
  h


Jf =  y  .
 
z



∂
∂r

i

∂
∂s

0



x

 h


−1
0
= R  y .


0
z

i.e.
Jf = R−1 Jf

∂
∂r

∂
∂s

i

0

(2.38)

(2.39)

where,


∂x0
∂r

∂x0
∂s




0

Jf = 


∂y 0
∂r

∂y 0
∂s

0

0






(2.40)

The procedure is summarized as,
1. Rotate the face so as to make it parallel to the xy plane.
0

2. Find the Jacobian of transformation from the rotated face to the standard face (Jf ).
3. Find the Jacobian of the face as given by equation 2.39.
q
f
4. |J | is given as | (Jf )T (Jf ) |
Method 2
In this method, we parametarize the surface x, y, z with parameters r, s. Let
X = xî + y ĵ + z k̂
We use the standard shape functions to map the physical element to parameters r, s. i.e.
x=

3
X
i=0

Ni xi ,

y=

3
X
i=0

10

Ni yi ,

z=

3
X
i=0

Ni zi ,

(2.41)

where same shape functions are used as earlier.
Note 2.4. In order to find the vertex mapping of a randomly oriented face with the
standard element, it would be advisable to first rotate the face to make it parallel to the xy
plane and then find the vertex mapping. Then find the shape function and proceed with
equation 2.41. However, we do not use the coordinates of the rotated face here.
After that, we find
are implemented in

∂X
∂r

and

∂X
.
∂s

Then |Jf | is simply

∂X
∂r

×

∂X
∂s

Both of these methods

Face::calculateJacobianQuadFace(double r, double s)
Any given quadrature point (r, s) can be mapped into (x, y, z) space using the two operations:
1. First find mapping from (r, s) to (x0 , y 0 , z 0 ) using equations 2.32-2.35 and 2.36.
2. Then map from (x0 , y 0 , z 0 ) to physical coordinates using R−1 .
The global positioning of all the quadrature points on each face is stored in array
Face::quadPointsGlobalLocation
Index i runs over the surface quadrature points and index j runs over the x, y, z axes.
Note 2.5. Some of the volume quadrature points coincide with the surface quadrature
points. Function
DG::mapFaceQuadPointsToCellQuadPoints()
performs the mapping by direct comparison of global positioning of surface and volume
quadrature points.
The mapping of surface to volume quadrature points is stored in arrays
Face::mapOwnerQuadPoints[ ] and Face::mapNeighbourQuadPoints[ ].
Thus, if we are finding a 2nd -order polynomial for a face, then it should have 9 quadrature points (QP). The corresponding cell will meanwhile have 27 quadrature points, 9 of
which will coincide with the face QP. Then if
Face::mapOwnerQuadPoints = [3, 2, 21, 14, ... , 7]
Then, it means the 0th quad point on face is same as the 3rd quad point of owner cell and
so on.
Note 2.6. In the improved version of the code, the face-vertex numbering is standardized,
such that the vertex numbering follows the right-hand rule. i.e. if the curled fingers
indicate the sequence of vertices, then the stretched thumb points in the direction of hte
face-normal. In that case, the rotation matrices are not required. i.e. both method 1 and
method 2 are obsolete. Refer:
Face::orderDefiningPoints()

11

2.4

Modal and Nodal Basis Functions

Basis functions are the independant dimensions on which any function can be projected.
In this chapter, we present a few modal and nodal basis functions.

2.4.1

Modal Basis Functions

These are the set of basis functions, where each dimension corresponds to a ‘mode’ or
frequency. A simple example is a set of monomials defined on the interval [a, b]
Note 2.7. This is just an example, in reality x could be valid over (−∞, ∞)
K
= {1, x, x2 , x3 , ..., xK−1 },
Vmono

x ∈ C 0 [a, b]

(2.42)

A polynomial up to order K − 1 can be represented exactly using K modal functions.
Thus, any general polynomial can be written as
K

P (x) =

i=K−1
X

αi xi , ,

x ∈ C 0 [a, b]

(2.43)

i=0

The monomials is not a ‘good’ choice of the basis functions as we will see later, thus
we need a set of basis functions where each of the functions is ‘orthogonal’ (or better
‘orthonormal’) to every other functions. One such important set of orthogonal functions
is called as ‘Jacobi’ polynomials. A special case of ‘Jacobi’ polynomials is ‘Legendre’
polynomials. In our code, we have extensively used Legendre polynomials. A general ith
modal basis function is indicated by symbol ψi . Other popular example of modal basis
functions is a set of Chebyshev polynomials.
legendre polynomials
ψ0
ψ1
ψ2
ψ3
ψ4
ψ5
ψ6
ψ7
ψ8

1.00
0.75
0.50

ψ

0.25
0.00
−0.25
−0.50
−0.75
−1.00
−1.00 −0.75 −0.50 −0.25

0.00

x

0.25

0.50

0.75

Figure 2.2: First 9 Legendre polynomials

12

1.00

Fig.2.2 shows first 9 Legendre polynomials. Observe that the ‘frequency’ of the polynomials increases (for increasing count) indicating higher modes. Fig.2.3 shows Chebyshev
polynomials.
chebyshev polynomials
ψ0
ψ1
ψ2
ψ3
ψ4
ψ5
ψ6
ψ7
ψ8

1.00
0.75
0.50

ψ

0.25
0.00
−0.25
−0.50
−0.75
−1.00
−1.00 −0.75 −0.50 −0.25

0.00

x

0.25

0.50

0.75

1.00

Figure 2.3: First 9 Chebyshev polynomials

The 1D modal functions are implemented in:
FunctionalSpace::monomials1D(args),
FunctionalSpace::legendrePolynomial1D(args),
FunctionalSpace::chebyshevPolynomial1D(args) and
FunctionalSpace::jacobiPolynomial1D(args)
It is indeed difficult to comprehend how the polynomials shown in Fig.2.2 and Fig.2.3 be
‘orthogonal’. The main reason for this is that we try to stick to our understanding of the
3D space and the notion of orthogonality that comes with it, however, the polynomial
space is essentially ∞ dimensional space. For this space, the definition of orthogonality
depends on the definition of the Inner-Product and ‘metric’ of that space. For L2 space,
the orthogonality is defined as,
Z
ψi ψj dΩ = αδij , ∀ψi , ψj ∈ L2 {Ω} & α ∈ R
(2.44)
Ω

where, δij is a ‘Kronecker Delta’ which means δij = 1 if i = j and 0 otherwise. The
functions ψi and ψj are orthonormal if α = 1.

2.4.2

Nodal Basis Functions

The ‘nodal’ basis functions are the set of functions which are ‘Cardinal’ in nature. What
that means is, at the interpolation point, the basis function will have value 1 and for all
13

other collocation points, the value of the basis function will be zero. Thus, the interpolated
value of the function agrees perfectly at the collocation points but may have some error
at other places. The general symbol used for nodal basis functions is φ. The interpolation
itself is given as,
N
−1
X
N
f (x) =
f (xi )φi (x), x ∈ C 0 [a, b]
(2.45)
i=0

and, φi (xj ) = δij . We use Lagrange polynomials as the nodal basis functions. The
polynomials in 1D are give as,
φi (x) =

N
−1
Y

(x − xj )
(xi − xj )
j=0,j6=i

(2.46)

it is clear that at x = xi , φi (xi ) = 1. One crucial point that needs to be taken care
of is selection of the collocation points [x0 , x1 , ..., xN −1 ]. It is seen that the equi-spaced
arrangement of collocation points gives poor interpolation for very high orders of accuracies. Therefore we take these points as roots of orthonormal modal polynomials. Popular choices for collocation points are Gauss-Legendre (LG) collocation points, LegendreGauss-Lobatto (LGL) points or Chebyshev points. An important distinction between LG
and LGL points is that, LGL points include the ‘end-points’ whereas LG points include
only the internal points. For DG methods, this doesn’t have a significant difference since
solving a Riemann problem at the face remains unchanged.

2.4.3

Multidimensional Interpolation

For number of dimensions greater than 1, the interpolating polynomials are obtained
depending on the type of the element we wish to interpolate on. For hex elements, these
are obtained using tensor product of 1D polynomials. i.e.
φ(x, y, z) = φ1D (x)φ1D (y)φ1D (z)

(2.47)

Whereas, for other types of elements, special multidimensional polynomials are used.
However, the concept of modal and nodal basis functions remains essentially the same
irrespective of the number of dimensions.

14

Chapter 3
Code Structure
Roughly the code structure is displayed in Fig.3.1. For more accurate and detailed description, run Doxygen utility and refer the generated documentation.
OF Input

GeometryIO

Tensor

Point

Cell

Functional Space

Face

Math Functions

DG

Figure 3.1: Class dependency (crudely)

Each of the the constituent elements are briefly described next.

15

3.1
3.1.1

Code structure
Tensor

MEANDG has a central datastructure called as ‘Tensor’, which forms all the important
arrays, matrices, multidimensional storage lists etc. Refer ‘include/Tensor’ for details.
Tensor class is made of 4 sub-classes.
TensorO1, TensorO2, TensorO3, Matrix
From these, O(1), O(2), O(3) tensors can be created of arbitrary dimensions. The tensor
objects interact with one another and have several in-build methods such as L1 , L2 , L∞
norms etc. All the maths functions accept objects of Tensor library.
Following example code shows usage of tensor class:
TensorO1 A(5); // creates a O(1) tensor of size 5 with
// data-members of the type double
TensorO1 C(5);
for (int i=0; i<5; i++){
// generate a random number between 0 and 100
double valueA = getRandom(0.0,100.0);
// generate another random number between 0 and 100
double valueC = getRandom(0.0,100.0);
// set value at ith index
A.setValue(i, valueA);
C.setValue(i, valueC);
};
// now A and C have random values
double dotAnswer = Math::dot(A,C); // perform math product between two arrays
cout << dotAnswer << endl; // print the answer

For further details and the available functions, refer the code documentation.

3.1.2

Point

Point is a physical point in (x, y, z) coordiante system. It has a unique id, and description
of its location.

3.1.3

Cell

If the domain is indicated by Ω, then the descritization of the domain (Ωh ) is described
as,
n
[
Ωh =
Si
(3.1)
l=1

16

where Si are the finite volume cells. These are non-overlapping pieces of the domain.
Each cell has the following important attributes (among others):
• Vertices (of the type ‘Point’)
• Faces (i.e. boundary of the cell)
• Degrees of Freedom (i.e. locations where unkowns are evaluated)
• Data (such as the conserved variables)
• Volume
• Mathematical constructs (Mass matrix, Stiffness matrix etc.)
• Mapping onto the standart cell (Jacobian etc.)
Thus, each of these parameters (among others) can be accessed using the accessor
(get, set) functions of the class ‘Cell’. For further details, refer Cell.h and Cell.cpp.

3.1.4

Face

A face is the boundary of a Cell. The face has following important attributes:
• Geometrical information (area, normal, vertices etc.)
• Neighbour and Owner cells
• Data (numerical flux and variables)
• Mapping (Jacobian etc.)

3.1.5

GeometryIO

This module reads the input file in the OpenFOAM format and fills the data arrays.
Refer:
GeometryIO::fillDataArrays(args)

3.1.6

Functional Space

This module is responsible for everything related to L2 , H 1 functional spaces. It contains
nodal, modal basis functions and derivatives. It also computes various matrices used in
DG formulation.

3.1.7

Math functions

These are various mathematical functions.

17

3.1.8

Gasdynamics, Riemann solvers etc.

Depending on the system of equations which one wants to solve, include the physicsspecific functions and files. For example, for Euler’s equations of Gasdynamics, we have
Gasdyanmics.cpp module and corresponding Riemann solvers. Apart from these, there
are test functions, display functions etc. Refer the code files in ‘include/’ and ‘src/’ folders.

3.1.9

DG

This is the core (central) module of the code. It performs following operations:
• Memory management
• Cells, faces creation
• Calling appropriate functions for computation of cell matrices
• Time integration
• PDE solving
Currently it is configured to solve the hyperbolic system of PDEs of the type:
∂Q
+ ∇.F(Q) = 0
∂t

18

(3.2)

Chapter 4
Applications

19

Chapter 5
License
MEANDG is licensed under BSD-3 clause license. Refer:
https://en.wikipedia.org/wiki/BSD_licenses
The library Eigen which is used in this software (see ./include/Eigen/) is licensed under
MPL2 (see ./include/Eigen/COPYING.README).
Use preprocessor flag -DEIGEN_MPL2_ONLY to ensure that only MPL2 (and possibly
BSD) licensed code from Eigen is compiled. Refer
http://eigen.tuxfamily.org/index.php?title=Main\_Page#License
for details. No files of the Eigen library are modified for use with MEANDG.

5.1

LICENSE

Copyright (c) 2018, Subodh Joshi, Tapan Mankodi, Shivasubramanian Gopalakrishnan
All rights reserved.
Redistribution and use in source and binary forms, with or without modification, are
permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice, this list of
conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice, this list
of conditions and the following disclaimer in the documentation and/or other materials
provided with the distribution.
3. Neither the name of the copyright holder nor the names of its contributors may be used
to endorse or promote products derived from this software without specific prior written
permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

20



---

Source Exif Data:

File Type                       : PDF
File Type Extension             : pdf
MIME Type                       : application/pdf
PDF Version                     : 1.5
Linearized                      : No
Page Count                      : 23
Producer                        : pdfTeX-1.40.16
Creator                         : TeX
Create Date                     : 2019:03:07 16:04:54+05:30
Modify Date                     : 2019:03:07 16:04:54+05:30
Trapped                         : False
PTEX Fullbanner                 : This is pdfTeX, Version 3.14159265-2.6-1.40.16 (TeX Live 2015/Debian) kpathsea version 6.2.1

EXIF Metadata provided by EXIF.tools