Usage
-----

In this section, we will look at a few basic examples of how to read
different mesh properties.

Reading structured grid information
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

In this section, we will see how we can read multi-block structured grid
information from the library within the ``meshReaderLib`` folder. In all
cases, we need to first create a structured grid instance and then read
the mesh. This is achieved with the following function calls:

.. code:: cpp

  #include <filesystem>
  #include "meshReaderLib/meshReader.hpp"

  int main() {
    // provide a relative or absolute path to where the CGNS file is located that should be read
    auto structuredMeshFile = std::filesystem::path("path/to/structured/mesh.cgns");

    // create an instance of the structured mesh reading class
    ReadStructuredMesh structuredMesh(structuredMeshFile);

    // read all structured grid properties at once
    structuredMesh.readMesh();

    // perform any other function calls here to get grid properties, discussed below in detail
    // ...

    return 0;
  }

After the structured mesh has been set up in this way, we can proceed to
read specific mesh information.

Reading coordinate information
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

To read the coordinates of a structured grid, we call the dedicated
getter function. An example of how to loop over the mesh is provided as
well, as well as how to get the ``x`` and ``y`` coordinates.

.. code:: cpp

  // get the coordinates from the structured mesh
  auto coordinates = structuredMesh.getCoordinates();

  // get the number of grid blocks (could be a multi-block structured grid)
  auto numberOfBlocks = coordinates.size();

  // loop over all blocks
  for (unsigned block = 0; block < numberOfBlocks; ++block) {
    // get the number of x-coordinate
    auto numberOfIndicesInX = coordinates[block].size();

    // loop over all x-coordinates
    for (unsigned indexX = 0; indexX < numberOfIndicesInX; ++indexX) {
      // repeat the same for the y-coordinate
      auto numberOfIndicesInY =  coordinates[block][indexX].size();

      // loop over all y-coordinates
      for (unsigned indexY = 0; indexY < numberOfIndicesInY; ++indexY) {
        const auto &x = coordinates[block][indexX][indexY][COORDINATE::X];
        const auto &y = coordinates[block][indexX][indexY][COORDINATE::Y];

        // process results and use x and y as required
        // ...

      }
    }
  }

Reading interface connectivity information
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

If we have more than a single zone/block in our mesh, then we need to
know at which interface these blocks are connected, so that we can
request the correct neighboring vertices in our calculation later, where
we update values close the the interface and where the computational
stencil requires information from a neighboring block. The code below
shows how to read this information.

.. code:: cpp

  // get the interface connectivity from the structured mesh
  auto interfaceConnectivity = structuredMesh.getInterfaceConnectivity();

  // get the number of interfaces in the mesh
  auto numberOfInterfaces = interfaceConnectivity.size();

  // loop over all interfaces
  for (unsigned interface = 0; interface < numberOfInterfaces; ++interface) {
    auto ownerBlock = interfaceConnectivity[interface].zones[0];
    auto neighbourBlock = interfaceConnectivity[interface].zones[1];

    // get the number of vertices in interface
    auto numberOfVerticesInInterface = interfaceConnectivity[interface].ownerIndex.size();

    // loop over interface
    for (int vertex = 0; vertex < numberOfVerticesInInterface; ++vertex) {
      // get owner indices in x and y (i.e. indices i, j)
      const auto &owner_i = interfaceConnectivity[interface].ownerIndex[vertex][COORDINATE::X];
      const auto &owner_j = interfaceConnectivity[interface].ownerIndex[vertex][COORDINATE::Y];

      // get neighbour indices in x and y (i.e. indices i, j)
      const auto &neighbour_i = interfaceConnectivity[interface].neighbourIndex[vertex][COORDINATE::X];
      const auto &neighbour_j = interfaceConnectivity[interface].neighbourIndex[vertex][COORDINATE::Y];

      // use vertices here, for example, discretising a second order partial derivative d2u/dx2 as
      // u^new = (u[i           + 1] - 2 * u[i]       + u[i       - 1]) / (dx * dx), may be written as
      // u^new = (u[neighbour_i + 1] - 2 * u[owner_i] + u[owner_i - 1]) / (dx * dx)
      // ...
    }
  }

Reading boundary condition information
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Finally, we want to apply specific boundary conditions at domain
boundaries. For this, we first need to know which boundary condition is
prescribed at the boundary. To read boundary information, use the code
below as a starting point.

.. code:: cpp

  // get the boundary condition information
  auto boundaryConditions = structuredMesh.getBoundaryConditions();

  // get the number of blocks in the mesh
  auto numberOfBlocks = boundaryConditions.size();

  // loop over all blocks
  for (unsigned block = 0; block < numberOfBlocks; ++block) {
    // check how many boundaries are within each block
    auto numberOfBoundaries = boundaryConditions[block].size();

    // loop over all boundaries in current block
    for (unsigned boundary = 0; boundary < numberOfBoundaries; ++boundary) {
      // get the boundary type, can be of type BC::WALL, BC::SYMMETRY, BC::INLET, or BC::OUTLET
      auto bcType = boundaryConditions[block][boundary].boundaryType;

      // get the number of vertices in boundary
      auto numberOfVertices = boundaryConditions[block][boundary].indices.size();

      // loop over all vertices in boundary
      for (unsigned vertex = 0; vertex < numberOfVertices; ++vertex) {
        // get current vertex
        const auto bcVertex = boundaryConditions[block][boundary].indices[vertex];

        // apply boundary condition based on specified type
        if (bcType == BC::WALL) {
          // process wall ...
        } else if (bcType == BC::SYMMETRY) {
          // process symmetry ...
        } else if (bcType == BC::INLET) {
          // process inlet ...
        } else if (bcType == BC::OUTLET) {
          // process outlet ...
        }
      }
    }
  }

Reading unstructured grid information
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

As we saw for the structured grid, we need to first set up an instance
of the unstructured mesh reading class and read the mesh itself.
Afterwards, we have access to the mesh information via getters just as
with the structured grid. The information is now in an unstructured grid
format and thus different from the structured grid information.
Additional information on the vertex connectivity is available as well
and the examples shown below indicate how to read an unstructured grid
with all of its information. But first, we have to set up the
unstructured grid as mentioned above.

.. code:: cpp

  #include <filesystem>
  #include "meshReaderLib/meshReader.hpp"

  int main() {
    // provide a relative or absolute path to where the CGNS file is located that should be read
    auto unstructuredMeshFile = std::filesystem::path("path/to/unstructured/mesh.cgns");

    // create an instance of the unstructured mesh reading class
    ReadUnstructuredMesh unstructuredMesh(unstructuredMeshFile);

    // read all unstructured grid properties at once
    unstructuredMesh.readMesh();

    // perform any other function calls here to get grid properties, discussed below in detail
    // ...

    return 0;
  }

Reading coordinate (vertex) information
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The unstructured grid exposes the ``getCoordinates()`` function which
provides the ``x`` and ``y`` coordinates for each vertex in the
unstructured grid. Since the unstructured grid is allowed to be split
into several zones/blocks, we need to loop over the number of blocks
first, just as we did with our structured grid as well. A complete
example of accessing the ``x`` and ``y`` coordinates for an unstructured
grid is shown below.

.. code:: cpp

  // get the coordinates from the unstructured mesh
  auto coordinates = unstructuredMesh.getCoordinates();

  // get the number of grid blocks (could be a multi-block unstructured grid)
  auto numberOfBlocks = coordinates.size();

  // loop over all blocks
  for (unsigned block = 0; block < numberOfBlocks; ++block) {
    // get the number of vertices
    auto numberOfVertices = coordinates[block].size();

    // loop over all vertices
    for (unsigned vertex = 0; vertex < numberOfVertices; ++vertex) {
      const auto &x = coordinates[block][vertex][COORDINATE::X];
      const auto &y = coordinates[block][vertex][COORDINATE::Y];

      // process results and use x and y as required
      // ...

    }
  }

Reading element connectivity table
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

As alluded to above, unstructured grids require additional element
connectivity information, which specifies how vertices are connected
into a cell/element. The unstructured grid class has the
``getInternalCells()`` function available for this, which exposes, for
each cell/element the vertex ID that makes up the current cell. We can
then get the location for each vertex in the cell by looking up the
information in the coordinate array that we read in the previous code
example.

.. code:: cpp

  // get the internal elements from the unstructured mesh
  auto interalElements = unstructuredMesh.getInternalCells();

  // get the number of grid blocks (could be a multi-block unstructured grid)
  auto numberOfBlocks = interalElements.size();

  // loop over all blocks
  for (unsigned block = 0; block < numberOfBlocks; ++block) {
    // get the number of internal elements in current block
    auto numberOfVertices = interalElements[block].size();

    // loop over all elements
    for (unsigned element = 0; element < numberOfVertices; ++element) {
      // get number of vertices making up current element
      auto numberOfVertices = interalElements[block][element].size();

      if (numberOfVertices == 3) {
        auto v0 = interalElements[block][element][0];
        auto v1 = interalElements[block][element][1];
        auto v2 = interalElements[block][element][2];
      
        std::cout << "triangle with vertices: " << v0 << ", " << v1 << ", " << v2 << std::endl;

        std::cout << "triangle has vertices at: ";

        std::cout << "(" << coordinates[block][v0][COORDINATE::X] << ", ";
        std::cout << coordinates[block][v0][COORDINATE::Y] << "), ";

        std::cout << "(" << coordinates[block][v1][COORDINATE::X] << ", ";
        std::cout << coordinates[block][v1][COORDINATE::Y] << "), ";

        std::cout << "(" << coordinates[block][v2][COORDINATE::X] << ", ";
        std::cout << coordinates[block][v2][COORDINATE::Y] << "), ";  
      }

      // process results and use x and y as required
      // ...

    }
  }

.. _reading-interface-connectivity-information-1:

Reading interface connectivity information
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Similar to the structured grid, we may have several zones or blocks in
our grid. At the interface between two zones/blocks, we need to know
which faces of the unstructured cells are connected to each other so
that we can request information from neighboring cells on different
zones/blocks for the calculation of flow properties. An example to
obtain the matching vertex pairs that make up a face is shown below.

.. code:: cpp

  // get the interface connectivity from the unstructured mesh
  auto interfaceConnectivity = unstructuredMesh.getInterfaceConnectivity();

  // get the number of blocks in the mesh
  auto numberOfBlocks = interfaceConnectivity.size();

  // loop over all interfaces
  for (unsigned block = 0; block < numberOfBlocks; ++block) {
    // get the number of interfaces in current block
    auto numberOfInterfaces = interfaceConnectivity[block].size();

    // loop over all interfaces in current block
    for (unsigned interface = 0; interface < numberOfInterfaces; ++interface) {
      auto ownerBlock = interfaceConnectivity[block][interface].zones[0];
      auto neighbourBlock = interfaceConnectivity[block][interface].zones[1];  

      // get the number of vertices in interface
      auto numberOfFacesInInterface = interfaceConnectivity[block][interface].ownerIndex.size();

      // loop over faces in interface
      for (int face = 0; face < numberOfFacesInInterface; ++face) {
        // get the starting vertex of the owning face
        const auto &ownerFaceStart = interfaceConnectivity[block][interface].ownerIndex[face][0];
        
        // get the ending vertex of the owning face
        const auto &ownerFaceEnd = interfaceConnectivity[block][interface].ownerIndex[face][1];

        // get the starting vertex of the neighbour face
        const auto &neighbourFaceStart = interfaceConnectivity[block][interface].neighbourIndex[face][0];
        
        // get the ending vertex of the neighbour face
        const auto &neighbourFaceEnd = interfaceConnectivity[block][interface].neighbourIndex[face][1];
      }
    }
  }

.. _reading-boundary-condition-information-1:

Reading boundary condition information
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Getting boundary conditions is also no different from reading boundary
conditions on a structured grid. The only difference here is that
instead of reading vertices that are on the boundary, we read vertex
pairs that make up the face of the boundary. The example below shows how
to extract the type of boundary, as well as the faces that are contained
within that boundary patch.

.. code:: cpp

  // get the boundary condition information
  auto boundaryConditions = unstructuredMesh.getBoundaryConditions();

  // get the number of blocks in the mesh
  auto numberOfBlocks = boundaryConditions.size();

  // loop over all blocks
  for (unsigned block = 0; block < numberOfBlocks; ++block) {
    // check how many boundaries are within each block
    auto numberOfBoundaries = boundaryConditions[block].size();

    // loop over all boundaries in current block
    for (unsigned boundary = 0; boundary < numberOfBoundaries; ++boundary) {
      // get the boundary type, can be of type BC::WALL, BC::SYMMETRY, BC::INLET, or BC::OUTLET
      auto bcType = boundaryConditions[block][boundary].boundaryType;

      // get the number of faces in boundary
      auto numberOfFaces = boundaryConditions[block][boundary].indices.size();

      // loop over all faces in boundary
      for (unsigned face = 0; face < numberOfFaces; ++face) {
        // get current starting location of boundary face
        const auto bcFaceStart = boundaryConditions[block][boundary].indices[face][0];

        // get current ending location of boundary face
        const auto bcFaceEnd = boundaryConditions[block][boundary].indices[face][1];

        // apply boundary condition based on specified type
        if (bcType == BC::WALL) {
          // process wall ...
        } else if (bcType == BC::SYMMETRY) {
          // process symmetry ...
        } else if (bcType == BC::INLET) {
          // process inlet ...
        } else if (bcType == BC::OUTLET) {
          // process outlet ...
        }
      }
    }
  }