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 ...
         }
       }
     }
   }