How to read a multi-block structured mesh from a CGNS file

In this article, we extend the work we have done in the previous article, where we have set up the basic library structure of our mesh reading library and started to provide an initial implementation of the base class. Specifically, we add a class that can read a multi-block structured mesh from a CGNS file. While we limit ourselves to 2D here, this can be easily extended to 3D (and we may get around one day to do so).

By the end of this article, you should have an excellent understanding of how to interact with the CGNS mid-level library and how we can use this to read the coordinates, interfaces, and boundary conditions from a CGNS file. You will start to build up an intuition for how the mid-level library calls are constructed and we even have a look at some advanced functionalities. I also make comments about how we could extend this class in case your mesh generator writes out CGNS files slightly differently to the one used in this series. Let’s get started.

Download Resources

All developed code and resources in this article are available for download. If you are encountering issues running any of the scripts, please refer to the instructions for running scripts downloaded from this website.

• Download: cgnsMeshReaderLib-Part-3.zip

In this series

• Part 1: What is the CGNS format and how to get started
• Part 2: How to inspect structured and unstructured grids using CGNS
• Part 3: How to set up a simple CGNS-based mesh reading library
• Part 4: How to read a multi-block structured mesh from a CGNS file
• Part 5: How to read a multi-block unstructured mesh from a CGNS file
• Part 6: How to test our CGNS-based mesh reading library

In this article

Overview

In the last article, we looked at setting up the basic library structure and proceeded with writing some code. We mainly concentrated on the ReadMeshBase class, which provides an interface for mesh reading classes to implement. It was structured in a way to read CGNS mesh files in particular and so provided already implementations to read the number of bases (and assertions that we only expect a single base), as well as code to read zone information (the number of vertices).

In this article, I want to continue the implementation of the mesh reading interface and write the code for the structured mesh reading class. This will require the addition of one header and one source file, as well as a few updates to existing files. The updated project tree is given below, and if you download the zip archive provided at the beginning of this article, you can inspect the project structure, as well as all of the code files as we go along.

meshReadingLibRoot
├── mesh
│   ├── structured2D.cgns
│   ├── structured2DNoFamily.cgns
│   ├── unstructured2D.cgns
│   └── unstructured2DNoFamily.cgns
├── meshReaderLib
│   ├── include
│   │   ├── readMeshBase.hpp
│   │   ├── readStructuredMesh.hpp    // new
│   │   └── types.hpp
│   ├── src
│   │   ├── readMeshBase.cpp
│   │   └── readStructuredMesh.cpp    // new
│   └── meshReader.hpp                // updated
├── buildAndRun.ps1                   // updated
├── buildAndRun.sh                    // updated
└── main.cpp                          // updated

We can see that we are going to add the meshReaderLib/include/readStructuredMesh.hpp header file, as well as the corresponding meshReaderLib/src/readStructuredMesh.cpp source file. These files will contain the definition and implementation of the ReadStructuredMesh class which will allow us, by the end of the article, to read a multi-block structured mesh.

Next, we have to provide a few updates to our build scripts, simply to include the added source file, we have to update our library header include file, so that we are including the new interface (class definition provided in the header file), and, finally, we want to update our main.cpp file to check that after we have read the mesh, we are not getting any error messages. Let’s look at the updated code first, and then we’ll jump into writing the structured mesh reading class.

Updating existing files

We have mentioned that a few files need updating. In this section, we will go through them quickly and see the changes we have to apply. Note that we have discussed each file in more depth in the previous article, so if some of the files do not make sense or you need additional explanations, have a look at the previous article which hopefully should clear things up.

Build scrips

We have to make minor tweaks to both the build script on Windows and UNIX to include the additional source file we are going to add. These are shown below

Windows

# clean up before building
Remove-Item .\build -Force -Recurse
New-Item -Name "build" -ItemType "directory"
Copy-Item "C:\libraries\bin\hdf5.dll" -Destination .\build
Copy-Item "C:\libraries\bin\zlib.dll" -Destination .\build

# compile source files into object files
cl.exe /nologo /EHsc /std:c++20 /I. /I "C:\libraries\include" /c /O2 /DCOMPILELIB .\meshReaderLib\src\readMeshBase.cpp /Fo".\build\readMeshBase.obj"
cl.exe /nologo /EHsc /std:c++20 /I. /I "C:\libraries\include" /c /O2 /DCOMPILELIB .\meshReaderLib\src\readStructuredMesh.cpp /Fo".\build\readStructuredMesh.obj"
cl.exe /nologo /EHsc /std:c++20 /I. /I "C:\libraries\include" /c /O2 /DCOMPILELIB .\main.cpp /Fo".\build\main.obj"

# link object files against CGNS library and its dependencies
cl.exe /nologo /EHsc /std:c++20 .\build\main.obj .\build\readMeshBase.obj .\build\readStructuredMesh.obj /Fe".\build\cgnsTest.exe" /link /MACHINE:x64 /LIBPATH:"C:\libraries\lib" cgns.lib hdf5.lib msvcrt.lib libcmt.lib
.\build\cgnsTest.exe

Here, we have added line 9 to include the meshReaderLib/src/readStructuredMesh.cpp source file during compilation. We then include the compiled object file, i.e. build/readStructuredMesh.obj, during the linking on line 13 as well, thus this line needs to be modified.

At the moment, we are compiling and linking both the library code and the main.cpp file at the same time. Once the library is completed, we are going to split both of these into separate compilation steps, i.e. we’ll compile the library code first into a static or dynamic library and we will then link our main.cpp file against that library, similar to how we set up our final library structure for our previous linear algebra solver library.

UNIX

#!/bin/bash

rm -rf build
mkdir -p build

# compile source files into object files
g++ -std=c++20 -I. -I ~/libs/include -c ./meshReaderLib/src/readMeshBase.cpp -o ./build/readMeshBase.o
g++ -std=c++20 -I. -I ~/libs/include -c ./meshReaderLib/src/readStructuredMesh.cpp -o ./build/readStructuredMesh.o
g++ -std=c++20 -I. -I ~/libs/include -c ./main.cpp -o ./build/main.o

# link object files against CGNS library and its dependencies
g++ -std=c++20 ./build/main.o ./build/readMeshBase.o ./build/readStructuredMesh.o -o ./build/cgnsTest -L~/libs/lib -lcgns

# run executable
./build/cgnsTest

Similar consideration to what we discussed for the Windows build script holds for UNIX (bash) as well. We are adding here line 8 to compile meshReaderLib/src/readStructuredMesh.cpp as part of the compilation step and then adding its object file during the linking on line 12. We’ll separate the compilation of the library and main.cpp file later as well for UNIX to have a separate library to link against once we compile the main.cpp file.

Header files

This change is a quick one, within the meshReaderLib/meshReader.hpp file, we simply add an include statement so that we are now including the class definition for the structured mesh reading class, i.e. line 2 #include "meshReaderLib/include/readStructuredMesh.hpp" is added here.

#include "meshReaderLib/include/readMeshBase.hpp"
#include "meshReaderLib/include/readStructuredMesh.hpp"
#include "meshReaderLib/include/types.hpp"

main.cpp for testing

The main.cpp file is the only modified file which is pretty much changed completely, at least in terms of its main content body. Let’s look at the file first and then discuss the changes.

#include <iostream>
#include <filesystem>
#include <cmath>
#include <string>

#include "meshReaderLib/meshReader.hpp"

int main() {
  auto structuredMeshPathFamily = std::filesystem::path("mesh/structured2D.cgns");
  ReadStructuredMesh structuredMeshFamily(structuredMeshPathFamily);
  structuredMeshFamily.readMesh();

  std::cout << "structured mesh was read" << std::endl;

  return 0;
}

Since we have included the structured mesh reading class header file in the meshReaderLib/meshReader.hpp file, we do not need to make any changes to our library include header, e.g. line 6. The main body itself, however, has entirely changed, mainly as we were not able to instantiate any classes in the previous article, as the base class contains pure virtual functions, which does not allow the creation of objects from these classes.

We still provide a path to the structured mesh CGNS file, where boundary conditions are written under a family node. If you are wondering what families in the context of a CGNS file are and how they relate to boundary conditions, I have provided a detailed discussion of them in my article on inspecting structured and unstructured CGNS files.

On line 10, we path that mesh file path as a constructor argument to the ReadStructuredMesh class and then we use the generated class object structuredMeshFamily to read the mesh on line 11. If everything goes well and no error messages are encountered, we should see the message "structured mesh was read" printed to the console. This is our first indication that there were no issues during the mesh reading. Later, we will test the mesh reading more formally, where we also inspect the mesh (e.g. coordinates, boundary conditions, interface data) that the class is providing us.

Creating a derived class for structured grids

In this section, we will look then at the new code we provide to read a structured mesh from a CGNS file. We start, as usual, defining the interface (class structure defined in the header file) first, and then we proceed to look at each function implementation within the source file.

Creating the class interface

The class definition, or the interface, is defined in the new header file located at meshReaderLib/include/readStructuredMesh.hpp. The file is given below in its full length and we’ll go through it afterwards and highlight important code sections.

#pragma once

#if defined(WIN32) || defined(_WIN32) || defined(__WIN32__) || defined(__NT__)
  #if defined(COMPILEDLL)
    #define MESHREADERLIB_EXPORT __declspec(dllexport)
  #elif defined(COMPILELIB)
    #define MESHREADERLIB_EXPORT
  #else
    #define MESHREADERLIB_EXPORT __declspec(dllimport)
  #endif
#else
  #define MESHREADERLIB_EXPORT
#endif

#include <filesystem>
#include <iostream>
#include <vector>
#include <array>
#include <unordered_map>
#include <string>
#include <cmath>

#include "cgnslib.h"

#include "meshReaderLib/include/readMeshBase.hpp"
#include "meshReaderLib/include/types.hpp"

class MESHREADERLIB_EXPORT ReadStructuredMesh : public ReadMeshBase {
public:
  using ZoneType = typename std::vector<std::vector<cgsize_t>>;
  using CoordinateType = typename std::vector<std::vector<std::vector<std::vector<double>>>>;
  using IndexType = typename std::array<unsigned, 2>;
  using InterfaceConnectivityType = typename std::vector<InterfaceConnectivity<IndexType>>;
  using BoundaryConditionInformationType = typename std::vector<std::vector<BoundaryConditionInformation<IndexType>>>;

public:
  ReadStructuredMesh(std::filesystem::path cgnsFilePath);
  void readMesh() override final;

  const CoordinateType& getCoordinates() const { return _coordinates; }
  const InterfaceConnectivityType& getInterfaceConnectivity() const { return _interfaceConnectivity; }
  const BoundaryConditionInformationType& getBoundaryConditions() const { return _boundaryConditions; }

protected:
  void readCoorinates() override final;
  void readInterfaceConnectivity() override final;
  void readBoundaries() override final;

private:
  CoordinateType _coordinates;
  InterfaceConnectivityType _interfaceConnectivity;
  BoundaryConditionInformationType _boundaryConditions;
};

Preprocessor directives

Lines 1-13 are identical to the interface of the base class. We mentioned in the base class article, that this is only required for Windows, really, and only if we are targeting dynamic libraries. If all we wanted to do is build a static library, that boilerplate code can go. Now we get here into some programming issue, in that we have just repeated ourselves (i.e. we have copied and pasted code from the base class into this class). This isn’t great. In reality, you would want to put this into one central place, which is typically the header include file, so in our case, that would be the meshReaderLib/meshReader.hpp file.

Take a look at the header include file of the CGNS library, for example, i.e. cgnslib.h, which we include on line 23. If you open that file (e.g. C:\libraries\include\cgnslib.h if you kept the same library installation path as me), then you’ll find the following lines somewhere around line 70 (which may change with future releases):

#ifndef CGNSDLL
# ifdef _WIN32
#  if defined(BUILD_DLL)
#    define CGNSDLL __declspec(dllexport)
#  elif defined(USE_DLL)
#    define CGNSDLL __declspec(dllimport)
#  else
#    define CGNSDLL
#  endif
# else
#  define CGNSDLL
# endif
#endif

This looks very familiar, doesn’t it? The CGNS library handles this a bit differently than our current library (there is an additional if statement, i.e. lines 1 to 13), but the essence is the same. If you then go further down the cgnslib.h file, you’ll find the functions declaration for reading zone information (as an example), of which two selected functions look like this:

CGNSDLL int cg_nzones(int fn, int B, int *nzones);
CGNSDLL int cg_zone_read(int fn, int B, int Z, char *zonename, cgsize_t *size);

You’ll see that all functions now start with CGNSDLL and so will either have __declspec(dllexport) or __declspec(dllimport) in front of them (for dynamic libraries on Windows) or an empty string for any other case. In our case, we simply have to modify the class name, i.e. we have

class MESHREADERLIB_EXPORT ReadStructuredMesh : public ReadMeshBase

Now, all functions inside the class will have the correct __declspec directive set and we don’t have to prepend any of the function names inside the class with the MESHREADERLIB_EXPORT preprocessor directive. We discussed the __declspec directives in more depth in the article on compiling the linear algebra library in the previous series, so have a look if you need a refresher.

Specified types

In the header file, we define a few types on lines 30-34, that we use to define variables on lines 50-52. I just noticed as well that we still have the ZoneType on line 30, this one should have been deleted (it has moved to the base class). The CoordinateType is a 4-dimensional std::vector, the first index goes over all zones, the second and third over all points in the x and y direction, and the fourth index contains the x and y coordinate. If we had 3D mesh reading, we would need to make this a 5D vector. We could have also made the last argument a std::array<double, 2>, since in this case it is not changing, but it would not make much of a difference here.

We define an IndexType on line 32 of type std::array<double, 2>, i.e. this index type tells us that if we want to index any point in the 2D mesh, we need two indices, which in CFD terminology would be the i, j index. We use this index type as a template argument on the next two lines when we construct the InterfaceConnectivityType and the BoundaryConditionInformationType. These, in turn, make use of the InterfaceConnectivity and BoundaryConditionInformation struct which we defined in the base class header file. For example, the InterfaceConnectivity struct had this form:

template<typename IndexType>
struct MESHREADERLIB_EXPORT InterfaceConnectivity {
  std::array<int, 2> zones;
  std::vector<IndexType> ownerIndex;
  std::vector<IndexType> neighbourIndex;
};

It is templated on the IndexType. In the ReadStructuredMesh class, we define the IndexType on line 32, use it on line 33 to instantiate a version of the InterfaceConnectivity struct, which means that when we create the _interfaceConnectivity variable of type InterfaceConnectivityType on line 51, we can now store i, j index pairs for the ownerIndex and neighbourIndex as we go through all points in the interface (hence the IndexType being used within a std::vector).

The InterfaceConnectivityType is stored in a 1D std::vector, while the BoundaryConditionInformationType is stored in a 2D std::vector. This has to do with how interfaces and boundaries are stored. Interfaces are stored in a global sense for structured grids (at least this is how Pointwise is generating the CGNS file, other mesh readers may store interfaces per zone, which is possible, and for unstructured grids, this is what Pointwise is doing, nothing we can do about). So, we read all interfaces and store them in the 1D std::vector we have defined.

Boundary conditions are available per zone, even if the boundary conditions are stored as a family under the base node. We read boundary conditions per zone and only if the name and type are stored under a family node, do we read that and then look it up while we do the boundary condition reading. The base class has already read any family-type boundary conditions and made them available for us, so if that is the case, we can simply ask for the information from the base class. But since boundary information is available per zone, we store the information in a 2D std::vector (i.e. the first index is over zones, and the second index is over all boundary conditions in that zone).

Function declarations

Let’s look at the function declarations next, there are only a few. I have copied them below so we don’t have to keep scrolling.

public:
  ReadStructuredMesh(std::filesystem::path cgnsFilePath);
  void readMesh() override final;

  const CoordinateType& getCoordinates() const { return _coordinates; }
  const InterfaceConnectivityType& getInterfaceConnectivity() const { return _interfaceConnectivity; }
  const BoundaryConditionInformationType& getBoundaryConditions() const { return _boundaryConditions; }

protected:
  void readCoorinates() override final;
  void readInterfaceConnectivity() override final;
  void readBoundaries() override final;

We define the constructor with the filename as an argument and we start to override all virtual functions from the base class. As a result, once all pure virtual functions are overridden, we can instantiate an object from this class, which we saw already when we modified the main.cpp tester function. The readMesh() function, in particular, is defined in the public area of the class, whereas all the other pure virtual functions, e.g. readCoorinates(), readInterfaceConnectivity(), and readBoundaries(), are defined in the protected area and thus only internal class methods can call these functions.

Then we have some controversial functions. The primary aim of a class is to implement some logic (hidden from the user, e.g. here read a mesh from a CGNS file) and then expose the end product to the user (e.g. here we want users to be able to access things like the coordinate, interface, and boundary condition arrays). We typically use so-called getter functions that allow users to get information from the class, and these functions typically start with the word get. We have 3 such functions on lines 5-7 above, but these were not defined in the base class.

Typically, you want to define all getters in the base class and then make sure that each derived class implements them. In this case, you can make sure that if a user requests some data, each class will implement that getter function and will return some data. In our case, we want to return data for either a structured or an unstructured grid. The data we will return will have a different type (mainly the IndexType is going to change, unstructured grids store indices in 1D arrays, i.e. they only store the index of a vertex). We could address this issue with templates at the expense of readability.

There is one more complication, though. The structured and unstructured grids return different information. The structured grid only exposes information for the coordinate, interface and boundary conditions. The unstructured grid provides the same information and additionally some element connectivity data. So we couldn’t define all getter functions in the base class, as there is no connectivity to return (we’ll look into that in more detail in the next article). So even if we introduced templates to store getters in the base class, we would still have to define additional getter functions in the derived class.

I have decided that for readability purposes, it is better to just define all getters in each derived class. If I didn’t care about documenting this class or I could expect that anyone reading my classes are all C++ experts, then I would probably find a more elegant solution with templates. Be aware that this is not a clean design but one which will help us understand our class better (hopefully).

Implementing the interface

Let us now turn to the source file and look at how the above-mentioned classes are implemented on a structured grid.

Constructor

The constructor is rather straightforward, we first include the header file of the class definition for the ReadStructuredMesh on line 1, and then we simply call the constructor and pass the filename, which we receive as a constructor argument, directly to the constructor of the base class, which here is the ReadMeshBase class. We pass the filename to the base class as it is the base class which is storing the filename.

#include "meshReaderLib/include/readStructuredMesh.hpp"

ReadStructuredMesh::ReadStructuredMesh(std::filesystem::path cgnsFilePath) : ReadMeshBase(cgnsFilePath) { }

Implementing the readMesh() function

Next, let’s look at the readMesh() function. The implementation is shown below.

void ReadStructuredMesh::readMesh() {
  readZones(CGNS_ENUMV(Structured));
  readCoorinates();
  readInterfaceConnectivity();
  readBoundaries();
}

The first thing we do is to read the zones. The argument we pass to the function is an enum, which tells the function that we expect the zones to be of type Structured and not Unstructured. We discussed the readZones() function in the previous article, which is implemented in the base class. We then simply call all the virtual functions that we are overwriting and implementing in this class. These definitions are given in the following.

Reading coordinates from a CGNS file

We start with the read coordinate function. The code definition is given below, and we will, as usual, discuss the implementation afterwards.

void ReadStructuredMesh::readCoorinates() {
  DataType_t type; char coordname[128];
  int numberOfDimensions = _zoneSizesMax[0].size();

  // resize the coordinates array, we store one per zone for both x and y
  _coordinates.resize(_numberOfZones);
  
  for (int zone = 0; zone < _numberOfZones; ++zone) {
    // set the number of max vertices for x and y. Coordinate reading will require this data structure later.
    cgsize_t maxVertices[2] = { _zoneSizesMax[zone][COORDINATE::X], _zoneSizesMax[zone][COORDINATE::Y] };
    cgsize_t minVertices[2] = { _zoneSizesMin[zone][COORDINATE::X], _zoneSizesMin[zone][COORDINATE::Y] };

    // based on the max vertices, resize the coordinates array for both x and y
    _coordinates[zone].resize(maxVertices[COORDINATE::X]);
    for (int i = 0; i < maxVertices[0]; ++i) {
      _coordinates[zone][i].resize(maxVertices[COORDINATE::Y]);
      for (int j = 0; j < maxVertices[1]; ++j) {
        _coordinates[zone][i][j].resize(numberOfDimensions);
      }
    }

    // loop over each direction, e.g. x and y (numberOfDimensions should be 2)
    for (int dimension = 0; dimension < numberOfDimensions; ++dimension) {
      // get the coordinates information, mainly the type and name. The type is either single or double preciosion and
      // the coordname is a SIDS-compliant identification, e.g. 'CoordinateX' and 'CoordinateY'      
      if (cg_coord_info(_fileIndex, 1, zone + 1, dimension + 1, &type, coordname)) cg_error_exit();

      // temporary (1D) coordinate array to store coordinates in. 2D std::vectors do no guarantee contiguous memory 
      // layout, thus we need to first read in all coordinates into a 1D array and then map that to a 2D array 
      std::vector<double> temp(maxVertices[COORDINATE::X] * maxVertices[COORDINATE::Y]);      
      if (cg_coord_read(_fileIndex, 1, zone + 1, coordname, type, minVertices, maxVertices, &temp[0])) cg_error_exit();

      // now map temporary (1D) coordinate array to 2D coordinate array
      for (int vertex = 0; vertex < maxVertices[COORDINATE::X] * maxVertices[COORDINATE::Y]; ++vertex) {
        unsigned i = vertex % maxVertices[0];
        unsigned j = vertex / maxVertices[1];
        _coordinates[zone][i][j][dimension] = temp[vertex];
      }
    }
  }
}

After defining a few variables on line 2, we seek the number of dimensions on line 3. We could have hard-coded this information here to 2, as it is a 2D mesh reading class, but we can also look into the _zoneSizesMax or _zoneSizesMin array, as they store the maximum and minimum index in each coordinate direction for a structured grid. For 2D, we have 2 indices to store the maximum number of points in the x and y direction, and so we get the number of dimensions in this way on line 3.

On line 6, we resize the coordinate array for the number of zones as we want to store all of them in the same array. We then loop over all zones on line 8 and start resizing the second and third dimensions of the coordinate array based on the maximum number of vertices in the x and y directions. Remember, this information is available for each zone, and we read that information in the base class when we were reading the zones.

We then loop over each dimension on line 23 and start reading the coordinates for each direction (here x and y). We first read coordinate information on line 26 which mainly provides us with the precision of the stored coordinates, i.e. either single- or double-precision. In C++ this maps to either a float or double type, respectively. Notice that we have to provide here the file index, base index, zone index and now also the coordinate index (or dimension, e.g. 1 for x and 2 for y) to the cg_coord_info() function call.

The coordname variable should return a name for the coordinate which is standardised, i.e. it should be either "CoordinateX" or "CoordinateY" depending on which coordinate we read. If the name is anything else, it is not a SIDS-compliant CGNS file. If we already know that the coordinates are stored as either single or double precision, than we can skip the function call on line 26 and go directly to lines 30-31.

We attempt to read the coordinates now on line 31. We have one problem. The CGNS library is written in C, so it can only work with C data types. We can still use C++ std::vector and then get the underlying raw C data structure (in the end, all the C++ data types store data with primitive C style types and arrays) but our coordinate vector is defined as a 4D std::vector. While a 1D std::vector is contiguous in space, a 4D vector isn’t, and so we can’t use this array here directly and have to create a new, temporary 1D std::vector to read the coordinates into, which we do on line 30.

Using cg_coord_read() on line 31, we can read the coordinate information, which expects a file index and zone index as usual, and then further the data type (single- or double-precision), coordinate name (CoordinateX or CoordinateY), the minimum and maximum vertex to read, and, finally, the 1D coordinate array itself.

Let’s assume we have a zone where we have 5 vertices in the x direction and 3 in the y direction. Then, the temporary coordinate array that we defined on line 30, which has a size of number of vertices in x * number of vertices in y, i.e. 15 in this example, will have all the x or y coordinates stored in a 1D array (depending on which coordinate we read on line 31). We now need to map this 1D array back into a 2D array and store it in our 4D global coordinate array.

We do this on lines 34-38. We loop over all the vertices in the 1D array and then calculate what the corresponding indices are for the 2D array. The i index is calculated as unsigned i = vertex % maxVertices[0];. maxVertices[0] is the maximum number of vertices in the x direction, in this example 5. The modulo operator % is a division but will return the remainder of the division. A good explanation of the modulo operator can be found on StackOverflow. So, in our example, i would be calculated as a series of 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4.

For the j direction, we have unsigned j = vertex / maxVertices[1];. Since maxVertices[1], i.e. the maximum number of vertices in the y direction, in our example 3, is an integer, integer division rules apply here. Integer division works by dividing the number as if we are dividing two floating point numbers and then rounding the result down. For example, 1/2 is 0.5, rounding it down is 0. 9/10 is 0.9, but this is still 0 after rounding it down. So, for the j index, we will get values of 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2 in our example.

Thus, the 2D mapping on lines 34-38 allows us to store the temporary 1D array with all the coordinates in our 4D array by finding the correct i and j indices first.

Reading interface connectivity for multi-block structured grids from a CGNS file

We have to break this part up into two pieces, as I want to discuss first why we have to read the interface and what information is available, as well as how we can read it with CGNS, as there are several possibilities. Afterwards, we’ll look at the mesh reading itself.

What interface connectivity information do we need and what does CGNS offer?

et’s remind ourselves what we are trying to achieve in this case, for which we need to look at our structured grid again. The figure below shows the structured grid itself, with all vertex pairs i,j for both zones and also the corresponding boundary conditions for the east, west, north, and south boundary (which we will use in the next section):

Between the two zones, i.e. the black and the red grid, we have an interface. If we loop over our mesh, then we will need information from our neighbouring cells to construct derivatives. For example, a simple finite difference discretisation for \partial^2 u/\partial x^2 may be written as \partial u/\partial x\approx (u_{i+1,j}-2u_{i,j}+u_{i-1,j})/(\Delta x)^2. If we are on the back grid at index 5,3, then we u_{i-1,j}=u_{4,3}^{zone1}, u_{i,j}=u_{5,3}^{zone1}, and u_{i+1,j}=u_{2,3}^{zone2}, i.e. the first two values can be read from zone 1, but the last value for u needs to come from zone 2. Once we know which index on zone 2 the index pair of 5,3 on zone 1 is connected to, (in this case, 1,3 in zone 2), we can find i+1,j, which in this case is the index pair 2,3.

This is what the interface will provide us. We see that on the interface, we have vertices that are part of both the black and red grid. Thus, we will have a set of i,j indices both for the black and red grid that index the same vertex. This connection between the indices is what the interface will provide us, so that once we calculate gradients and request information from neighbouring cells, we have a way to obtain this.

The CGNS mid-level library defines 4 general types of connectivity: 1-to-1, general, periodic, and overset holes. Yes, we’ve found it, the periodic boundary condition is hidden here, under connectivity, not boundary conditions. In our case, we have a 1-to-1 connectivity (each vertex in the black grid as a corresponding vertex in the red grid, or, we have a conformal connection). For the 1-to-1 connectivity type, there are 5 different functions defined:

• cg_n1to1_global(): Get the total number of 1-to-1 interfaces in a database
• cg_1to1_read_global(): Read data for all 1-to-1 interfaces in a database
• cg_1to1_write(): Write 1-to-1 connectivity data for a zone
• cg_n1to1(): Get the number of 1-to-1 interfaces in a zone
• cg_1to1_read(): Read 1-to-1 connectivity data for a zone

You see that there are two options here, either we have the interface written for each zone, in which case we loop over all zones and then read interface information with cg_1to1_read(), or, all information is stored globally, i.e. under the base node and not under each zone node, in which case we have to read all connectivity information with cg_1to1_read_global(). Which connectivity information is available to us depends on the grid generation software.

What I have done here is to call cg_n1to1_global() first and then loop over all zones and call cg_n1to1() afterwards. Then I looked at which of these calls was giving me a total interface of 1. In this case, Pointwise (my mesh generation software) stored interfaces as global interfaces, so I have to use cg_1to1_read_global() to get all connectivity. This is somewhat annoying, as I would have preferred to store it under each zone, but this is what we have to do.

before we get into the code, I noticed now that I did not implement a fail-fast strategy in this case, if you want to make the following code more robust, you should also loop over all cells and assert that there are no interface connectivity stored under a zone (because we are not reading it). I would implement it as well, but can’t generate a mesh file that produces an interface under each zone, so I wouldn’t be able to test the implementation.

CGNS implementation of reading global interface connectivity information.

With this out of the way, let’s have a look at the code, which is a bit longer, and then discuss it afterwards.

void ReadStructuredMesh::readInterfaceConnectivity() {
  int numberOfInterfaces;

  // read the number of interfaces in the mesh
  if (cg_n1to1_global(_fileIndex, 1, &numberOfInterfaces)) cg_error_exit();
  
  // allocate temporary memory. This is required to ensure memory layout is contiguous
  char **interfaceName = new char*[numberOfInterfaces];
  char **zoneName = new char*[numberOfInterfaces];
  char **donorName = new char*[numberOfInterfaces];
  cgsize_t **owner = new cgsize_t*[numberOfInterfaces];
  cgsize_t **neighbour = new cgsize_t*[numberOfInterfaces];
  int **transform = new int*[numberOfInterfaces];

  int numberOfDimensions = _zoneSizesMax[0].size();
  for (int interface = 0; interface < numberOfInterfaces; ++interface) {
    interfaceName[interface] = new char[128];
    zoneName[interface] = new char[128];
    donorName[interface] = new char[128];
    owner[interface] = new cgsize_t[2 * numberOfDimensions];
    neighbour[interface] = new cgsize_t[2 * numberOfDimensions];
    transform[interface] = new int[numberOfDimensions];
  }
  
  // read all connectivity (interface) information all at once
  if (cg_1to1_read_global(_fileIndex, 1, interfaceName, zoneName, donorName, owner, neighbour, transform))
    cg_error_exit();

  _interfaceConnectivity.resize(numberOfInterfaces);
  for (int interface = 0; interface < numberOfInterfaces; ++interface) {
    // store zone index of the interface owner and neighbour 
    _interfaceConnectivity[interface].zones[INTERFACE::OWNER] = _zoneNamesMap[std::string(zoneName[0])];
    _interfaceConnectivity[interface].zones[INTERFACE::NEIGHBOUR] = _zoneNamesMap[std::string(donorName[0])];

    // construct (temporary) 2x2 transformation matrix based on information read from the interface. Details about the
    // transformation matrix can be found here: https://cgns.github.io/CGNS_docs_current/sids/cnct.html#Transform
    std::vector<std::vector<int>> tMatrix(numberOfDimensions, std::vector<int>(numberOfDimensions));
    for (int dimension = 0; dimension < numberOfDimensions; ++dimension) {
      auto value = transform[interface][dimension];
      tMatrix[std::abs(value) - 1][dimension] = (value > 0) ? 1 : ((value < 0) ? -1 : 0);
    }
    
    // number of vertices in the interface alogn the x and y direction
    unsigned numberOfIndicesX = std::abs(owner[interface][2] - owner[interface][0]);
    unsigned numberOfIndicesY = std::abs(owner[interface][3] - owner[interface][1]);
    unsigned numberOfIndices  = std::max(numberOfIndicesX + 1, numberOfIndicesY + 1);

    // a map of indices, mapping i,j indices from the owning zone to the i,j indices of the neighbouring zone
    _interfaceConnectivity[interface].ownerIndex.resize(numberOfIndices);
    _interfaceConnectivity[interface].neighbourIndex.resize(numberOfIndices);

    for (unsigned index = 0; index < numberOfIndices; ++index) {
      const unsigned ownerStartX = owner[interface][COORDINATE::X] - 1;
      const unsigned ownerStartY = owner[interface][COORDINATE::Y] - 1;
      const unsigned neighbourStartX = neighbour[interface][COORDINATE::X] - 1;
      const unsigned neighbourStartY = neighbour[interface][COORDINATE::Y] - 1;

      const unsigned indexX = numberOfIndicesX == 0 ? 0 : index;
      const unsigned indexY = numberOfIndicesY == 0 ? 0 : index;

      unsigned ownerIndexX = ownerStartX + indexX;
      unsigned ownerIndexY = ownerStartY + indexY;

      // calculate the i, j indices of the neighbour zone with the transformation matrix and owner indices. See
      // https://cgns.github.io/CGNS_docs_current/sids/cnct.html#Transform for more information about this equation
      unsigned neighbourIndexX = tMatrix[0][0] * indexX + tMatrix[0][1] * indexY + neighbourStartX;
      unsigned neighbourIndexY = tMatrix[1][0] * indexX + tMatrix[1][1] * indexY + neighbourStartY;

      // write connectivty information to interface connectivity structure
      _interfaceConnectivity[interface].ownerIndex[index][COORDINATE::X] = ownerIndexX;
      _interfaceConnectivity[interface].ownerIndex[index][COORDINATE::Y] = ownerIndexY;
      _interfaceConnectivity[interface].neighbourIndex[index][COORDINATE::X] = neighbourIndexX;
      _interfaceConnectivity[interface].neighbourIndex[index][COORDINATE::Y] = neighbourIndexY;
    }
  }

  // deallocate memory to avoid memory leakage
  for (int interface = 0; interface < numberOfInterfaces; ++interface) {
    delete [] interfaceName[interface];
    delete [] zoneName[interface];
    delete [] donorName[interface];
    delete [] owner[interface];
    delete [] neighbour[interface];
    delete [] transform[interface];
  }
  delete [] interfaceName;
  delete [] zoneName;
  delete [] donorName;
  delete [] owner;
  delete [] neighbour;
  delete [] transform;
}

C-style array memory allocation for CGNS functions

There are two places in the code that we can from the beginning glance over, lines 7-23 and lines 78-91. We spend 30 lines of code here just to allocate and deallocate memory. I hate that it takes that much space and I hate that I have to do it. If you are using the C++ STL containers like a std::vector, and I have outlined my reasons for why you should do so as well, then this waste of space is just endlessly annoying. But CGNS is written in C, for a good reason (it can be used with a variety of programming languages), in which case primitive data types are important.

So, this memory allocation and deallocation are variables that we need to read the global interface connectivity information. We see here on lines 7-23 that we have to allocate memory for raw 2D pointers using C-style arrays. We have two dimensions here so that we can loop over each interface (first index/dimension) and then for each interface we have to store several entries in an array (second index/dimension), which we will see in the remaining section of the code.

Reading interface connectivity information in the interface

Before that mess of a memory allocation, we have to read the number of interfaces that are globally stored under the base node. We do that on line 5. Since they are under the base node, we only have to provide the file index and base index as an argument (and then receive the number of interfaces as an output). Remember that we asserted before that we only have 1 base, so using here a hard-coded base index of 1, as we did in other places, will not break the code.

On line 26, we can read all interface connectivity information at once, as we have made space (memory allocation) before for all of these variables. Let’s look at the arguments in a bit more detail that we get out for each interface:

cg_1to1_read_global(_fileIndex, 1, interfaceName, zoneName, donorName, owner, neighbour, transform)

The inputs are again the file index, base index (1), and then we have a few outputs. We get the interface name, zone and donor name. In CGNS’s terminology, whenever we have an owner/neighbour situation (like for an interface, where there will be two interfaces), the neighbour zone/grid/etc. will always be identified as the donor. So when we say zone and donor name, we mean the names of the zones from the owning and neighbouring (or left and right, top and bottom, front and back, etc.) grid.

We only get the names of the zones, but not the actual zone indices, which would be far easier, hence we set up a zone name to zone index map during the zone reading in the base class. If you need a reminder, look at the zone reading in the base class implementation.

The owner and neighbour arrays store the actual i,j indices for both zones which is the interface information we are after. Though, we will only use one of them and then calculate the other, for reasons we will see shortly. We’ll need the transform for this as well, which provides a compact transformation matrix, that tells us how the two zones are orientated to each other. If they both use locally the same coordinate system as the global coordinate system, then we technically don’t need this transformation, but to make sure we can support any type of global interface connection, we have to implement it.

We don’t actually care about the interface name, but we do want to store the owner and neighbour zone index. This is what we do on lines 32-33 and we make use of the aforementioned zone name to zone index map which we set up previously to get the zone index from a given zone name.

Constructing the transformation matrix for owner-neighbour indices

Next we look at the transformation, for which we are constructing a temporary matrix. There is a bit of background reading required in the SIDS which you should read before continuing here. Otherwise you may get lost as it is not that straightforward.

The transform variable that we received from the interface reading (last argument). It stores a compact transformation matrix. In our case, we only have two dimensions, but let’s assume that we had a 3D code and a 3D interface, we may have received a transform array with the following three entries transform = [-2, +3, +1]. This is the same example as in the SIDS linked to above. Then, we can construct the full 3×3 transformation matrix as:

\begin{bmatrix}0 & 0 & 1 \\ -1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}

So, the number in the transform array tells us which row to index (ignoring the sign), i.e. the first column has a value of -2, which means that we index the second row. The second column has a value of 3, meaning we index the third row, etc. Each of these indexed position get a one, all other entries get a zero. The sign within the transform array is then used to determine if it is plus or minus one in the matrix. If all coordinate systems have the same orientation, then we have the identity matrix, i.e. transform = [1, 2, 3].

The above logic is implemented in the code on lines 37-41, which is copied below to avoid excessive scrolling.

std::vector<std::vector<int>> tMatrix(numberOfDimensions, std::vector<int>(numberOfDimensions));
for (int dimension = 0; dimension < numberOfDimensions; ++dimension) {
  auto value = transform[interface][dimension];
  tMatrix[std::abs(value) - 1][dimension] = (value > 0) ? 1 : ((value < 0) ? -1 : 0);
}

First, we create a 2×2 matrix called tMatrix. We loop over all dimensions (which we could have also called columns, i.e. we are looping over all columns here) and then we store a temporary value which is just the value from the transform array (to make the next line more readable). On line 4 above, we index the matrix with the same logic described above, i.e. the row is given by the value in the transform array, and we take the absolute value. Since C++ is 0-based, we have to subtract 1. Then we simply check if the value is larger or smaller than 0, and write +1 or -1 accordingly.

If the above shorthand if/else syntax looks confusing, have a look at learncpp and their write-up, as always, they have a rather good explanation of the topic.

Calculate the number of indices in the interface

Lines 44-46 determine the number of vertices that we have in the interface. This is repeated below for convenience:

unsigned numberOfIndicesX = std::abs(owner[interface][2] - owner[interface][0]);
unsigned numberOfIndicesY = std::abs(owner[interface][3] - owner[interface][1]);
unsigned numberOfIndices  = std::max(numberOfIndicesX + 1, numberOfIndicesY + 1);

Both the variables owner and neighbour, which we read on line 26 in the cg_1to1_read_global() mid-level library call, have 4 entries, as seen on lines 20-21, where we allocated twice the number of dimensions (the number of dimensions is 2 due to a 2D mesh). The values stored in these arrays are the start and end indices in the interfaces on both the owning and neighbouring zones. The structure of the array is [start_i, start_j, end_i, end_j]. Thus, if we want to know how many vertices there are in the x (i) direction, we have to calculate end_i-start_i, and, end_j-start_j in the y (j) direction, respectively. This is done on lines 1-2 above.

Since we are on a structured grid, at the interface, one of these indices will be zero, as either the i or j coordinate will be constant. Looking back at our structured grid example above, we have a constant value for the x coordinate and thus i=constant. For the index pairs shown above in the structured grid, we have start_i=5, start_j=1, end_i=5, end_j=5. We start our indexing again with 1, as CGNS files use Fortran-based indexing rather than zero-based C-style indexing.

Using these numbers, we have numberOfIndicesX = 0 and numberOfIndicesY = 4 on lines 1-2 above. But, since we know there are 5 vertices in the interface, we have to add 1 to each of these values, which we do on line 3. Line 3 is then also checking which value is the largest, i.e. one of them will always be 1 and the other will always be the number of vertices/indices in the interface, which we store in the numberOfIndices variable on line 3. Knowing the number of vertices, allows us to allocate memory in our _interfaceConnectivity variable on lines 49-50 that will store all i,j index pairs in the interface.

Constructing i,j index pairs of owner/neighbour vertices using the transformation matrix

We have now the number of indices for which we need to generate pairs of i,j indices for both the owning and neighbouring zones, and this is what is done on lines 52-74 above. I have again copied it below so we avoid scrolling back and forth.

for (unsigned index = 0; index < numberOfIndices; ++index) {
  const unsigned ownerStartX = owner[interface][COORDINATE::X] - 1;
  const unsigned ownerStartY = owner[interface][COORDINATE::Y] - 1;
  const unsigned neighbourStartX = neighbour[interface][COORDINATE::X] - 1;
  const unsigned neighbourStartY = neighbour[interface][COORDINATE::Y] - 1;

  const unsigned indexX = numberOfIndicesX == 0 ? 0 : index;
  const unsigned indexY = numberOfIndicesY == 0 ? 0 : index;

  unsigned ownerIndexX = ownerStartX + indexX;
  unsigned ownerIndexY = ownerStartY + indexY;

  // calculate the i, j indices of the neighbour zone with the transformation matrix and owner indices. See
  // https://cgns.github.io/CGNS_docs_current/sids/cnct.html#Transform for more information about this equation
  unsigned neighbourIndexX = tMatrix[0][0] * indexX + tMatrix[0][1] * indexY + neighbourStartX;
  unsigned neighbourIndexY = tMatrix[1][0] * indexX + tMatrix[1][1] * indexY + neighbourStartY;

  // write connectivty information to interface connectivity structure
  _interfaceConnectivity[interface].ownerIndex[index][COORDINATE::X] = ownerIndexX;
  _interfaceConnectivity[interface].ownerIndex[index][COORDINATE::Y] = ownerIndexY;
  _interfaceConnectivity[interface].neighbourIndex[index][COORDINATE::X] = neighbourIndexX;
  _interfaceConnectivity[interface].neighbourIndex[index][COORDINATE::Y] = neighbourIndexY;
}

We first identify the starting index in x and y (i.e. i, and j) on lines 2-5. We determine if the index is increasing in the x or y direction on lines 7-8. If the number of indices in the x direction is zero, then we leave the indexX variable set to zero, otherwise, if it is not zero, then we set it equal to the index, which is the loop variable specified on line 1 in the for loop. This essentially just checks if we are looping along the i or j direction in the interface.

Next, we have to set and compute the index pairs of i,j values that allow us to index vertices on the interface for both the owning and neighbouring zones. We set the owner x and y indices on lines 10-11 and then compute the corresponding i,j index pair for the neighbour zone on lines 15-16. This step may seem confusing and I would encourage you to read the documentation on the transformation in the SIDS again. But let’s consider a special case, which should make this calculation easier to understand.

I mentioned previously when we constructed the transformation matrix that, if all zones have the same local coordinate system, then we have a transformation matrix that is the identity matrix. If we have the identity matrix, then all the off-diagonal components are zero and lines 15-16 above would reduce to

unsigned neighbourIndexX = tMatrix[0][0] * indexX + neighbourStartX;
unsigned neighbourIndexY = tMatrix[1][1] * indexY + neighbourStartY;

Keeping further in mind that the values inside the transformation matrix are either positive or negative one (and, in this case, they are positive one since we have the identity matrix, i.e. not transformation), then we can further simplify the above to:

unsigned neighbourIndexX = indexX + neighbourStartX;
unsigned neighbourIndexY = indexY + neighbourStartY;

Comparing this with the way we set the owner zone indices (lines 10-11):

unsigned ownerIndexX = ownerStartX + indexX;
unsigned ownerIndexY = ownerStartY + indexY;

The order is different, but we see that if there is no transformation between the two different zones, then we essentially just take the start index in i and j, and then we add the loop index indexX and indexY to it in both cases. If you have a look at the SIDS linked above, you will find an example where you have two non-matching coordinate systems for two different zones, in which case the transformation is required.

Finally, on lines 19-22 we store now all the i,j index pairs for both the owning and neighbouring zones in the interface, at which point we have gathered all the information in the interface that we want to store.

Reading boundary conditions from a CGNS file

Let’s look again at the detailed drawing of our structured grid and try to establish what we want to achieve here.

We see that we have 4 different boundary conditions for the east (outlet), west (inlet), north (symmetry), and south (wall) boundary on this grid. Both the north and south boundary is split among the two zones. We want to read now for each zone what the boundary condition is. We expect, for a 2D mesh, that we have 4 boundary conditions per zone. This number is reduced if we have interfaces as well. In this case, both structured zones have 1 interface, so we expect to read 3 boundary conditions per interface.

The full boundary condition reading is given below. We will break this code up into separate parts and then discuss the various aspects of the code in the following discussion.

void ReadStructuredMesh::readBoundaries() {
  _boundaryConditions.resize(_numberOfZones);
  for (int zone = 0; zone < _numberOfZones; ++zone) {
    // get the number of boundary conditions for the current zone
    int numberOfBoundaries;
    if (cg_nbocos(_fileIndex, 1, zone + 1, &numberOfBoundaries)) cg_error_exit();
    _boundaryConditions[zone].resize(numberOfBoundaries);

    for (int boundary = 0; boundary < numberOfBoundaries; ++boundary) {
      // temporary boundary condition information;
      char boundaryName[128]; BCType_t boundaryType; PointSetType_t pointSetType; cgsize_t numberOfPointsInBoundary;
      int normalIndex, numberOfDataset, normalList; cgsize_t normalListSize; DataType_t normalDataType;

      // read boundary condition information
      if (cg_boco_info(_fileIndex, 1, zone + 1, boundary + 1, boundaryName, &boundaryType, &pointSetType,
        &numberOfPointsInBoundary, &normalIndex, &normalListSize, &normalDataType, &numberOfDataset)) cg_error_exit();

      assert(pointSetType == PointRange && "Expected to read a start and end location for boundary conditions only");

      // start and end location for boundary conditions in i, j index
      cgsize_t boundaryPoints[2][2]; 

      // read boundary condition start and end location
      if (cg_boco_read(_fileIndex, 1, zone + 1, boundary + 1, boundaryPoints[0], &normalList)) cg_error_exit();

      // number of vertices in the interface along the x and y direction
      unsigned numberOfIndicesX = std::abs(boundaryPoints[1][0] - boundaryPoints[0][0]);
      unsigned numberOfIndicesY = std::abs(boundaryPoints[1][1] - boundaryPoints[0][1]);
      unsigned numberOfIndices  = std::max(numberOfIndicesX + 1, numberOfIndicesY + 1);
      _boundaryConditions[zone][boundary].indices.resize(numberOfIndices);

      // get all indices along current boundary condition and store them in boundary structure. Reconstructed based on
      // the start and end location of the current boundary.
      for (unsigned index = 0; index < numberOfIndices; ++index) {
        const unsigned indexX = numberOfIndicesX == 0 ? 0 : index;
        const unsigned indexY = numberOfIndicesY == 0 ? 0 : index;
        _boundaryConditions[zone][boundary].indices[index][COORDINATE::X] = boundaryPoints[0][COORDINATE::X] + indexX - 1;
        _boundaryConditions[zone][boundary].indices[index][COORDINATE::Y] = boundaryPoints[0][COORDINATE::Y] + indexY - 1;
      }
      
      // Writing the boundary name and type to the boundary structure. We need to check if we need to get this
      // information from a family or if we can read this directly.
      if (boundaryType == CGNS_ENUMV(FamilySpecified)) {
        // go to the current boundary node within the ZoneBC_t node, which must be present for each zone
        if (cg_goto(_fileIndex, 1, "Zone_t", zone + 1, "ZoneBC_t", 1, "BC_t", boundary + 1,"end")) cg_error_exit();
      
        // once we have navigated to the boundary node, read the family name, this will be checked against the family
        // names we read previously. From the family name, we can get the boundary condition type through the map we set
        // up earlier.
        char familyName[128];
        if (cg_famname_read(familyName)) cg_error_exit();

        // store boundary condition name and type in the boundary condition structure
        _boundaryConditions[zone][boundary].boundaryName = familyName;
        _boundaryConditions[zone][boundary].boundaryType = getBoundaryConditionID(_boundaryConditionMap[familyName]);
      } else {
        // store boundary condition name and type in the boundary condition structure
        _boundaryConditions[zone][boundary].boundaryName = boundaryName;
        _boundaryConditions[zone][boundary].boundaryType = getBoundaryConditionID(boundaryType);
      }
    }
  }
}

Reading the boundary condition information

By the time we read the boundary conditions, we know how many zones we have, and we want to store the boundary conditions for each zone separately. So we first make space for as many zones as we have on line 2 and then loop over all zones on line 2.

Lines 5-7 reads the number of boundary conditions per zone. Remember, that number should be between 0-4, depending on how many interfaces we have. If we have only a single zone, then we would expect to have exactly 4 boundary conditions. We could add these as additional assertions into the code, but I guess I have thought of this too late. But feel free to modify your code.

On line 9, we loop over all boundaries in this zone and then define some variables on lines 11-12 that we need to read boundary conditions from the corresponding CGNS mid-level library functions.

Line 15 reads some general information about the current boundary. Notice that we index the first base, zone+1 and now also boundary+1, which we defined on line 9 as the for loop variable. In particular, we are interested here boundaryName, boundaryType, and pointSetType. The name is what we gave during mesh generation and lets us identify the boundary later if we want to separate between boundaries. The type is one of the 21 different CGNS boundary conditions, and the point set type is either of type point range or point list (other possibilities exist but these are the most common choices).

A point list will give us all the i,j indices that are attached to a certain boundary, and a point range will give us the start and end point in the boundary, similar to what we saw for the interface (where we also only received the start and end indices and had to construct the indices in-between ourselves).

The remaining variables are of no interest to us. You may surprised to see that we don’t care about the numberOfPointsInBoundary, this is because for a point range, this value will always be 2 (i.e. we have a start and end location (2 entries)). If we receive a point list instead, this variable would return the number of indices in the boundary. In this case, Pointwise returns the boundaries as a point list.

Constructing i,j index pairs for each boundary condition

On lines 18-39, we need to construct now the i,j index pairs for each boundary condition and store that in our boundary condition struct. These lines are repeated below for convenience.

assert(pointSetType == PointRange && "Expected to read a start and end location for boundary conditions only");

// start and end location for boundary conditions in i, j index
cgsize_t boundaryPoints[2][2]; 

// read boundary condition start and end location
if (cg_boco_read(_fileIndex, 1, zone + 1, boundary + 1, boundaryPoints[0], &normalList)) cg_error_exit();

// number of vertices in the interface along the x and y direction
unsigned numberOfIndicesX = std::abs(boundaryPoints[1][0] - boundaryPoints[0][0]);
unsigned numberOfIndicesY = std::abs(boundaryPoints[1][1] - boundaryPoints[0][1]);
unsigned numberOfIndices  = std::max(numberOfIndicesX + 1, numberOfIndicesY + 1);
_boundaryConditions[zone][boundary].indices.resize(numberOfIndices);

// get all indices along current boundary condition and store them in boundary structure. Reconstructed based on
// the start and end location of the current boundary.
for (unsigned index = 0; index < numberOfIndices; ++index) {
  const unsigned indexX = numberOfIndicesX == 0 ? 0 : index;
  const unsigned indexY = numberOfIndicesY == 0 ? 0 : index;
  _boundaryConditions[zone][boundary].indices[index][COORDINATE::X] = boundaryPoints[0][COORDINATE::X] + indexX - 1;
  _boundaryConditions[zone][boundary].indices[index][COORDINATE::Y] = boundaryPoints[0][COORDINATE::Y] + indexY - 1;
}

Since we expect a point range (start and end indices) rather than a point list (all indices), we first assert that we will indeed receive the boundary information as a point range and not a point list. This is done on line 1. We then set up a boundaryPoints array on line 4 where we store the start and end locations for the i and j direction. We read this information on line 7 (we can again ignore here the remaining arguments, in this case just the normalList).

Similar to how we found the number of vertices along an interface, we are now using the same technique to find the number of indices/vertices for each boundary condition on lines 10-12. We find the total number of indices (vertices) on line 12 using the same comparison we used earlier for the interface. On line 13, we allocate space for all index pairs i,j for the current boundary and then on lines 17-22 first construct the i,j index (lines 18-19) and then store it inside the boundary condition variable on lines 20-21.

The logic on lines 18-19 follows exactly the same logic that we used during the interface i,j index pair creation; we check if the indices should increase in the i or j direction and then add only increase the index along either the i or j direction.

Extending the boundary condition reading for point lists

For completeness, if your mesh generator returns a point list instead, then your code simplifies significantly. Without testing the code below (as pointwise does not return a point list), you would simplify the code to

assert(pointSetType == PointList && "Expected to read all i,j index pairs");

std::vector<cgsize_t> boundaryPoints(2 * numberOfPointsInBoundary);

// read boundary condition with point list information
if (cg_boco_read(_fileIndex, 1, zone + 1, boundary + 1, &boundaryPoints[0], &normalList)) cg_error_exit();

// map 1D array of boundary points to 2D i,j index pair
for (unsigned index = 0; index < 2 * numberOfPointsInBoundary; ++++index) {
  _boundaryConditions[zone][boundary].indices[index][COORDINATE::X] = boundaryPoints[index + COORDINATE::X];
  _boundaryConditions[zone][boundary].indices[index][COORDINATE::Y] = boundaryPoints[index + COORDINATE::Y];
}

First, we assert that we have a point list instead of a point range, and then we set up a 1D std::vector that will read all i,j pairs into a 1D array and then we loop over all vertices in the boundary. I am assuming here that the numberOfPointsInBoundary variable will return the number of vertices in a boundary (in our case 5, for example, for each boundary). If that is the case, then we need to loop over twice as many entries, as we have both i and j stored. This is done on line 9 and we also have to increase the loop counter by 2, instead of 1, hence ++++index and not ++index.

Remember that COORDINATE::X=0 and COORDINATE::Y=1, so when we index into the boundaryPoints on lines 10-11, where index will always increase by (i.e. we have index=[0, 2, 4, 6, ...]), we can index the i and j coordinates in the 1D boundaryPoints array using this notation.

Storing boundary condition name and type

Finally, we want to store the boundary name and type in our boundary condition variable. This is done on lines 41-60 in the above-shown implementation and repeated again below for convenience:

// Writing the boundary name and type to the boundary structure. We need to check if we need to get this
// information from a family or if we can read this directly.
if (boundaryType == CGNS_ENUMV(FamilySpecified)) {
  // go to the current boundary node within the ZoneBC_t node, which must be present for each zone
  if (cg_goto(_fileIndex, 1, "Zone_t", zone + 1, "ZoneBC_t", 1, "BC_t", boundary + 1,"end")) cg_error_exit();

  // once we have navigated to the boundary node, read the family name, this will be checked against the family
  // names we read previously. From the family name, we can get the boundary condition type through the map we set
  // up earlier.
  char familyName[128];
  if (cg_famname_read(familyName)) cg_error_exit();

  // store boundary condition name and type in the boundary condition structure
  _boundaryConditions[zone][boundary].boundaryName = familyName;
  _boundaryConditions[zone][boundary].boundaryType = getBoundaryConditionID(_boundaryConditionMap[familyName]);
} else {
  // store boundary condition name and type in the boundary condition structure
  _boundaryConditions[zone][boundary].boundaryName = boundaryName;
  _boundaryConditions[zone][boundary].boundaryType = getBoundaryConditionID(boundaryType);
}

We first check if the boundary condition is stored under a family node or the boundary node within the zone. If the boundary type, read before with the cg_boco_info() function, is of type FamilySpecified, then we know that we have to go to the family node and read the boundary information from there. This is what we check on line 3.

If it is stored under a family, then we have to use the cg_goto() function, which allows us to index a specific node directly. It is similar like changing a directory in your file system, we go to a specific folder and then afterwards we can do some specific operations just in that folder. So in this case, we navigate to the first base, the current zone (zone + 1) and for the current zone we go to the first boundary node (there should just be one, i.e. all boundary information should be stored under a single ZoneBC_t node per zone. In this boundary node, we go to the current boundary condition, which is indexed with boundary + 1.

Each call to cg_goto() has to end with the string "end" to signal that we have reached our desired location. And, if you are wondering how you can find out the type of each node, these are given in the mid-level library documentation, the SIDS, and, if you prefer the lazy approach (I do), then you can also ask the CGNS file to tell you the type of each node by running cgnslist as we have seen in our article on inspecting structured and unstructured grids. You will need to run cgnslist with the -a flag, i.e. cgnslist -a meshFile.cgns. Then, you will not just get the structure of the file but for each node you also get the type of the node. Use that to construct your cg_goto() path.

Once we have found the correct node in the CGNS file, we can start reading the information that is available in this particular node. We do that on line 11, where we read the family name that is storing the corresponding boundary conditions. Using the family boundary condition lookup table that we have created in the base class (see the previous article), we can query what boundary condition should be used for the current family name.

On line 14, we store the family name, which is the boundary condition name we specified in the mesh generator, and then we store the boundary condition type, which we first have to get from our boundary condition lookup table (mentioned above, i.e we map a family name to a boundary condition as discussed in the previous article) and then we further map the corresponding CGNS boundary condition (of which there are 21) to one of the 4 boundary conditions that we currently support with our mesh reading.

Should the boundary condition be available directly under the boundary condition node, instead of under the family node, then we can use the read boundary name and boundary condition directly (which we obtained earlier when we called the cg_boco_info() function). This information is then directly stored in our boundary condition variable, as shown on lines 18-19. We have to, again, map one of the 21 CGNS-defined boundary conditions to one of 4 that we have defined in our mesh reading classes. The getBoundaryConditionID() is defined in the base class and was discussed in the previous article.

Summary

And there you have it, if you arrived at this point and feel comfortable with the write-up above, then you should now be in a good position to read a structured grid from a CGNS file. There are a few things that you could add here, for example, adding ghost cells if your structured solver requires them, but Pointwise does not offer support for writing these (as far as I am aware) so we skipped this part.

I hope you have obtained a good understanding and some intuition for how to interact with the mid-level library of the CGNS library. We always make calls to the library using a set of indices, where we specify which node in the CGNS file we want to read, and receive some information back, that we can then either directly store or further process into data that is useful to us, as we saw in the case of reconstructing i,j index pairs in the interface and on boundaries.

In the next article, we will tackle the topic of reading an unstructured grid. There are many similarities so if you understood this article, understanding unstructured mesh reading will not get much more complicated. The complication only arises when we are using the unstructured grid in a CFD solver, where the discretisation becomes a bit more involved, but this is a topic best saved for another time.

---

Tom-Robin Teschner is a senior lecturer in computational fluid dynamics and course director for the MSc in computational fluid dynamics and the MSc in aerospace computational engineering at Cranfield University.