Object orientation in C++
-----------------------------
In this section, a quick review on how object orientation works in C++ is given, afterwards the code architecture will be set up.

To get familiar with the concepts of object orientation, have a quick read of `this one page <https://www.tutorialspoint.com/cplusplus/cpp_object_oriented.htm>`_. 
It briefly describes *Objects*, *Classes*, *Abstraction*, *Encapsulation*, *Inheritance*, *Polymorphism* and *Overloading*. All of these concepts will be used in the first exercise.
We will exemplarily discuss object orientation with the following examples. 

Class Array2D - An example
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The first given class, `Array2D`, describes a two-dimensional field of values with given size. The values can be accessed by two indices :math:`(i,j)`. The file `array2d.h` has the following contents:

.. code-block:: c++
  :linenos:

  #pragma once

  #include <vector>
  #include <array>

  /** This class represents a 2D array of double values.
   *  Internally they are stored consecutively in memory.
   *  The entries can be accessed by two indices i,j.
   */
  class Array2D
  {
  public:
    //! constructor
    Array2D(std::array<int,2> size);

    //! get the size
    std::array<int,2> size() const;

    //! access the value at coordinate (i,j), declared not const, i.e. the value can be changed
    double &operator()(int i, int j);

    //! get the value at coordinate (i,j), declared const, i.e. it is not possible to change the value
    double operator()(int i, int j) const;

  protected:

    std::vector<double> data_;  //< storage array values, in row-major order
    const std::array<int,2> size_;    //< width, height of the domain
  };

Line 1 contains the `include guard`. Usually, header files will be included in multiple source and other header files, e.g. by `#include "array2d.h"`. The first stage of the compiler is to do `preprocessing`. It replaces all `#include` s by the content of the included file. The files that include other files form a graph and it is obvious that preprocessing can lead to the same header file getting copied multiple times to the same source file. This would mean that classes with the same name would be defined multiple times which is not allowed in C++. Here, `#pragma once` helps. It instructs the preprocessor to only include the respecitive file once, i.e. only the first time a corresponding `#include` is called.

The class definition contains a `public` section, lines 12-23 and a `protected` section, lines 25-28. Methods defined `public` can be called from anywhere where we have an object of the class. Protected methods and data members are only accessible from methods of this class or inherited classes. In this case, we have no protected methods, only the two members, `data_` and `size_`. 

Typically, a class does not have public data members, only public methods. This is in accordance with the `Encapsulation` principle. In order to access or change any data members, you need to call a method of the class. Methods that do not change any data members can (should) be declared `const`, as can be seen in lines 17 and 23.

Note the special methods in lines 20 and 23 with name `operator()`. This overloads the parantheses operator of the object. You can have code like this:

.. code-block:: c++
  :linenos:

  std::array<int,2> size{2,2}   
  Array2D array2D(size);        // object of class Array2D with size 2x2
  
  array2D(0,1) = 1.0;
  double value = array2D(0,1);

In lines 4 and 5, the `operator()` will be called with the parameters `i` and `j` set to 0 and 1. The non-const version of line 20 will be used in line 4, because the value 1.0 gets assigned to the internal member, the const version of line 23 will be used in line 5, because the value is only retrieved but not changed.

In the following, the implementation of the class is given, the file name is "array2d.cpp" (Both are stored in a subdirectory `storage`).

.. code-block:: c++
  :linenos:
  
  #include "storage/array2d.h"

  #include <cassert>

  Array2D::Array2D(std::array<int,2> size) :
    size_(size)
  {
    // allocate data, initialize to 0
    data_.resize(size_[0]*size_[1], 0.0);
  }

  //! get the size
  std::array<int,2> Array2D::size() const
  {
    return size_;
  }

  double &Array2D::operator()(int i, int j)
  {
    const int index = j*size_[0] + i;

    // assert that indices are in range
    assert(0 <= i && i < size_[0]);
    assert(0 <= j && j < size_[1]);
    assert(j*size_[0] + i < (int)data_.size());

    return data_[index];
  }

  double Array2D::operator()(int i, int j) const
  {
    const int index = j*size_[0] + i;

    // assert that indices are in range
    assert(0 <= i && i < size_[0]);
    assert(0 <= j && j < size_[1]);
    assert(j*size_[0] + i < (int)data_.size());

    return data_[index];
  }

If you read the code, you may notice the `assert` functions. These have a condition as argument. If the condition is false, the program will crash with an error in this line. In this example, it would be the case if the indices `i,j` were out of range, which should never happen, if the program has no bugs. The asserts are only in action in `debug` target, in `release` target, they are removed by the compiler.

A second note is on the memory layout of the stored data. `data_` is a `std::vector<double>`, i.e. a one-dimensional array of `double` values. You may have seen C code, where arrays were created like this:

.. code-block:: c++
  :linenos:
  
  double a[2];    // one-dimensional c-style array
  double b[2][2]; // two-dimensional c-style array
  double *c;      // only a pointer
  c = new double [2];   // allocate memory for 2 entries

You are explicitely discouraged to use these constructs. They deal with raw memory pointers or explicit allocation (line 4). If something fails here, it can easily result in undefined behaviour. This means that the program often works as it should and only fails sometimes and the crash occurs at a completely different location that where it was caused such that the bug is very hard to find.
Better use `std::array` and `std::vector` where the former is used when its size is fixed at compile-time while the latter is for when the size is only known at run-time. 

We need to transform the 1D array to a 2D array. There are two obvious ways to do this: store entries row-by-row (row-major) or column-by-column (column-major). None of them has advantages over the other. The presented code uses row-major storage. Care should be taken when accessing all data one after each other in a loop. When traversing all elements of the 2D array, it is more efficient to access the elements in their natural storage order than jumping around between memory locations. This should be considered when designing the two nested for loops for `i` and `j` that iterate over the array.

Inheritance
^^^^^^^^^^^^^^^^
Inherited classes have all members and methods of their parent (=base) class. They add can some more data members and methods. C++ code of an inherited class can look like the following:

.. code-block:: c++
  :linenos:
  
  ...
  /** A field variable is the discretization of a scalar function f(x) with x in the computational domain.
   *  More specifically, a scalar value is stored at discrete nodes/points. The nodes are arranged in an equidistant mesh
   *  with specified mesh width. 
   */
  
  class FieldVariable :
   public Array2D
  {
    ...
  
Polymorphism
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Using polymorphism it is possible to define functional "interfaces" of what a class can do. There might be different implementations to this interface. Which one to use can be chosen at run-time, for example depending on some parameter value in the "settings" file.

One example within the NumSim code is the different discretizations. There is the normal `Central Differences` discretization and the `Donor Cell` scheme with upwind differences.

The implementation of an interface is usually an *abstract base class* with virtual methods. The following is an exemplary `Discretization` interface that defines virtual methods to compute various derivatives. "= 0" means that the method has no implementation in this class.

.. code-block:: c++
  :linenos:
  
  #pragma once

  #include "discretization/0_staggered_grid.h"

  class Discretization : 
    public StaggeredGrid
  {
  public:

    //! construct the object with given number of cells in x and y direction
    Discretization(std::array<int,2> nCells, std::array<double,2> meshWidth);

    //! compute the 1st derivative ∂ u^2 / ∂x
    virtual double computeDu2Dx(int i, int j) const = 0;

    //! compute the 1st derivative ∂ v^2 / ∂x
    virtual double computeDv2Dy(int i, int j) const = 0;

    //! compute the 1st derivative ∂ (uv) / ∂x
    virtual double computeDuvDx(int i, int j) const = 0;
    
  ...
  
Then there exist two implementations that both derive/inherit from `Discretization`: `CentralDifferences` and `DonorCell`. E.g. the header file for `DonorCell` begins like this:

.. code-block:: c++
  :linenos:
  
  #pragma once

  #include "discretization/1_discretization.h"

  class DonorCell : public Discretization
  {
  public:

    //! use the constructor of the base class
    DonorCell(std::array<int,2> nCells, std::array<double,2> meshWidth, double alpha);

    //! compute the 1st derivative ∂ u^2 / ∂x
    virtual double computeDu2Dx(int i, int j) const;

    //! compute the 1st derivative ∂ v^2 / ∂x
    virtual double computeDv2Dy(int i, int j) const;

    //! compute the 1st derivative ∂ (uv) / ∂x
    virtual double computeDuvDx(int i, int j) const;
    
  ...
  
There is an actual implementation of the methods in the corresponding source file. The same holds for `CentralDifferences`, only the implementation of the methods differ.

Now, we have the interface class `Discretization` and the two derived classes `CentralDifferences` and `DonorCell`. The following can be done: 

* Create an object of type `CentralDifferences`
* Get a pointer to it
* Cast the pointer to a pointer of the interface class, `Discretization`
* Call any method of the pointer of the interface class. The call to the method will be *dispatched* to the original `CentralDifferences` object.

Pointers to variables/objects of a given `Type` are represented by `std::shared_ptr<Type>` or `std::unique_ptr<Type>`. (`unique_ptr` cannot be copied, so better just always use `shared_ptr`.)
There exist also the old C-style pointers, e.g. `double *a`. With them you had to do memory allocation and deallocation yourself. This often lead to errors that were hard to find or to memory leaks. In C++, these old pointers are considered bad practice, so please don't use them.

The following example demonstrates the use of `std::shared_ptr`.

.. code-block:: c++
  :linenos:
  
  std::shared_ptr<CentralDifferences> centralDifferences1 = std::make_shared<CentralDifferences>(settings_.nCells, meshWidth_, settings_.alpha);
  double value1 = centralDifferences1->computeDu2Dx(0,1);
  
  std::shared_ptr<Discretization> centralDifferences2 = std::make_shared<CentralDifferences>(settings_.nCells, meshWidth_, settings_.alpha);
  double value2 = centralDifferences2->computeDu2Dx(0,1);
  
In line 1 a pointer to an object of class `CentralDifferences` is defined. `std::make_shared` creates a new object of this class and return the pointer to it which is then assigned to the pointer `centralDifferences1`. The methods of this object can now be access using the "`->`" operator. In line 2 the method `computeDu2Dx` is called with parameters `i=0` and `j=1`. The result is stored in the double `value1`.

Lines 4-5 demonstrate the use of polymorphism. This time, a pointer of the base class, `Discretization`, is defined in line 4. Again, `std::make_shared` creates a new `CentralDifferences` object. The pointer to it is downcasted to the pointer of type `Discretization` and assigned to `centralDifferences2`. But the object it points to is still of type `CentralDifferences`. Therefore, in line 5, we can access the method `computeDu2Dx`, too.
  
This makes it possible to have code like this:
  
.. code-block:: c++
  :linenos:
  
  std::shared_ptr<Discretization> discretization;
  
  // create discretization
  if (settings_.useDonorCell)   // depending on a settings value
  {
    // create donor cell discretization
    discretization = std::make_shared<DonorCell>(settings_.nCells, meshWidth_, settings_.alpha);
  }
  else 
  {
    // create central differences
    discretization = std::make_shared<CentralDifferences>(settings_.nCells, meshWidth_);
  }