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

In this article, we continue developing our CGNS mesh reader library and implement the unstructured mesh reading class to complement our structured mesh reading class developed in the previous article. There are a few similarities between reading a structured and unstructured grid, and where that is the case, I will go through the code a bit faster and refer back to the structured mesh reading. This allows me to concentrate more on the differences, which are mainly in the element connectivity array that we need to read and process.

By the end of this article, you will have a good understanding of how to read an unstructured mesh from a CGNS file. Unstructured grids do not have to be written into several blocks/zones, as we saw for example with structured grids. However, to be as general as possible, we do allow for multi-block unstructured grids here and do read interface information as well if present. Therefore, the way that mesh reading is described here is as universal as possible and you should be able to read any unstructured grid with this implementation.

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

At this point, we have created the library structure and the base class, which includes some shared functionality, and in the previous article, we implemented a class that is capable of reading multi-block structured grids. In this article, I want to further extend the capability of our library by reading unstructured grids as well. For this, we will generate another derived class from our base mesh reading class, and then look at the implementation.

There is quite a bit of overlap between this article and the previous one on reading a structured grid. Even if your intentions are to only read unstructured grids, the previous article contains a lot of explanations that will not be repeated here to keep things to the point (and, if you have gone through the previous two articles already, you’ll appreciate this). There are a few new things we have to discuss in this article, specifically, reading cells which require element connectivity data, and we will focus on this area.

Let’s go through the files we need to update first and then concentrate on the new unstructured mesh reading class.

Updating existing files

As we did in the previous example, we need to update a few files. Specifically, we need to make modifications to the build scripts, the header include file, and the main.cpp file, which will now include the unstructured mesh reading class.

Build scripts

Let’s start with the build scripts. We look at both the Windows and UNIX build scripts below.

Windows

For the buildAndRun.ps1 build script, we need to add the usntructured mesh reading class to the compilation and linking stage. In the updated build script below, we have added a compilation statement on line 10 to compile the unstructured mesh reading class, while we have added its created object file on line 14 for the linker.

# 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 .\meshReaderLib\src\readUnstructuredMesh.cpp /Fo".\build\readUnstructuredMesh.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 .\build\readUnstructuredMesh.obj /Fe".\build\cgnsTest.exe" /link /MACHINE:x64 /LIBPATH:"C:\libraries\lib" cgns.lib hdf5.lib msvcrt.lib libcmt.lib
.\build\cgnsTest.exe

UNIX

In a similar manner, we update the UNIX bash script buildAndRun.sh to build our executable. Line 9 was added to compile the unstructured mesh reading class and line 13 was modified to include the compiled object file of the unstructured mesh reading class.

#!/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 ./meshReaderLib/src/readUnstructuredMesh.cpp -o ./build/readUnstructuredMesh.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 ./build/readUnstructuredMesh.o -o ./build/cgnsTest -L~/libs/lib -lcgns

# run executable
./build/cgnsTest

Header files

As we did for the structured grid class, we have to add the unstructured mesh reading class to our header include file, i.e. the meshReaderLib/meshReader.hpp file. This is done on line 3 below.

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

main.cpp for testing

Finally, we create a similar test for reading unstructured grids in the main.cpp file. In the file shown below, we essentially copy and paste the code for reading a structured mesh on lines 9-13 and copy that afterwards. Replacing the class to ReadUnstructuredMesh and reading one of the unstructured CGNS files, we now have a quick test for reading unstructured grids. Any errors during development may then be detected if this call produces an error.

#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;

  auto unstructuredMeshPathNoFamily = std::filesystem::path("mesh/unstructured2DNoFamily.cgns");
  ReadUnstructuredMesh unstructuredMeshNoFamily(unstructuredMeshPathNoFamily);
  unstructuredMeshNoFamily.readMesh();

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

  return 0;
}

What information is stored in an unstructured grid?

Before we jump into the code, I thought it would be good to have a look at the unstructured grid that we are trying to read, so that we have a clear idea of what information an unstructured grid exposes. Below is a schematic drawing of the unstructured grid that we will be reading with our implementation (located in the mesh/ folder) :

There are similarities to the structured grid example that we looked at in our previous article. We have two zones, the bounding box is from 0,0 to 2,1 (i.e. each edge for either of the zones is exactly 1 unit in length, here meters), and we have the same boundary conditions assigned. The interface is at the same location, albeit with fewer vertices that connect it now.

We see a couple of indices in the unstructured grid, the job of the mesh reading library is to read all of that information. The first task is to read all the vertices. Each vertex will have a coordinate for x and y. The vertices shown above show the zone they belong to with a superscript, and the current index within that zone with their subscript (again, indices in CGNS are Fortran-based and start at 1). With the knowledge of the bounding box being 0,0 to 2,1, we can then read the coordinates for the different vertices. For example, we have v_2^1=(0.5,0), v_{11}^1=(0.75,0.25), and v_9^2=(1.5,0.5).

During the coordinate reading, we will read the coordinates of these vertices into 1D arrays. Having the coordinates does not provide any information about the cells they belong to, though, and this information is provided by the element connectivity array (sometimes also called an element/cell lookup table). The element connectivity array will store the indices of the vertices that make up a cell. These cells are shown above using the identifier C_{index}^{zone}. The subscript and superscript are defined the same as for vertices.

An example of element connectivity that makes up a cell is C_8^1=[v_1^1, v_2^1, v_{12}^1] and C_3^2=[v_2^2, v_3^2, v_4^2, v_9^2]. We can see that each cell simply stores the indices of each vertex. Should we now need to construct cell information such as the cell’s centroid, its volume or one of its face area and normal vector, we simply go into the coordinate array and query the coordinates for each of the vertices in the element connectivity array for that specific cell. There are a few more steps here to get the information we need, but that is the nature of an unstructured grid that we have to accept to facilitate easier handling of complex geometries.

Interfaces and boundaries are similar to a structured grid, where we store indices in the interface or boundary, now as a 1D array of indices, rather than 2D arrays of i,j index pairs as we did for structured grids. As long as you feel comfortable with this discussion, the remaining write-up of how to read an unstructured grid should not be too difficult.

Creating a derived class for unstructured grids

With an overview of unstructured grids in general and the updating of existing files out of the way, let’s then have a look at the ReadUnstructuredMesh class which has its interface implemented in the header file at meshReaderLib/include/readUnstructuredMesh.hpp and its source file located at meshReaderLib/src/readUnstructuredMesh.cpp.

Creating the class interface

The interface or class definition is provided below for the ReadUnstructuredMesh class. As we did for the structured grid mesh reading class, we’ll break this into different components and discuss them in turn below.

#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 <string>
#include <cmath>
#include <unordered_map>
#include <algorithm>
#include <stdexcept>
#include <tuple>

#include "cgnslib.h"

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

class MESHREADERLIB_EXPORT ReadUnstructuredMesh : public ReadMeshBase {
public:
  using CoordinateType = typename std::vector<std::vector<std::vector<double>>>;
  using IndexType = std::vector<cgsize_t>;
  using InterfaceConnectivityType = typename std::vector<std::vector<InterfaceConnectivity<IndexType>>>;
  using BoundaryConditionInformationType = typename std::vector<std::vector<BoundaryConditionInformation<IndexType>>>;
  
  using InternalElementConnectivityType = typename std::vector<std::vector<IndexType>>;
  using InterfaceElementConnectivityType = typename std::vector<std::vector<std::vector<IndexType>>>;
  using BoundaryElementConnectivityType = typename std::vector<std::vector<std::vector<IndexType>>>;
  using StartLocationType = typename std::vector<IndexType>;
  using StartLocationForElementConnectivityType = typename std::tuple<std::vector<std::vector<cgsize_t>>,
    std::vector<std::vector<cgsize_t>>, std::vector<std::vector<cgsize_t>>>;

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

  const CoordinateType& getCoordinates() const { return _coordinates; };
  const InternalElementConnectivityType& getInternalCells() const { return _internalElementConnectivity; };
  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:
  void readCellConnectivity();
  StartLocationForElementConnectivityType getStartLocationForElementConnectivity() const;

private:
  CoordinateType _coordinates;
  InterfaceConnectivityType _interfaceConnectivity;
  BoundaryConditionInformationType _boundaryConditions;

  InternalElementConnectivityType _internalElementConnectivity;
  InterfaceElementConnectivityType _interfaceElementConnectivity;
  BoundaryElementConnectivityType _boundaryElementConnectivity;
};

Preprocessor directives

Let’s concentrate on the preprocessor directives first. I’ll also include header include statements here.

#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 <string>
#include <cmath>
#include <unordered_map>
#include <algorithm>
#include <stdexcept>
#include <tuple>

#include "cgnslib.h"

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

We see the same dynamic library import/export statements on lines 3-13 that we have discussed at length previously. However, we see also the include statements on lines 15-28, which are a bit longer than before. We’ll make use of the tuple header here which allows us to return more than one argument from a function. We’ll make use of that during the element connectivity reading.

Specified types

Lines 32-42 define a few types which are listed below.

using CoordinateType = typename std::vector<std::vector<std::vector<double>>>;
using IndexType = std::vector<cgsize_t>;
using InterfaceConnectivityType = typename std::vector<std::vector<InterfaceConnectivity<IndexType>>>;
using BoundaryConditionInformationType = typename std::vector<std::vector<BoundaryConditionInformation<IndexType>>>;

using InternalElementConnectivityType = typename std::vector<std::vector<IndexType>>;
using InterfaceElementConnectivityType = typename std::vector<std::vector<std::vector<IndexType>>>;
using BoundaryElementConnectivityType = typename std::vector<std::vector<std::vector<IndexType>>>;
using StartLocationType = typename std::vector<IndexType>;
using StartLocationForElementConnectivityType = typename std::tuple<std::vector<std::vector<cgsize_t>>,
  std::vector<std::vector<cgsize_t>>, std::vector<std::vector<cgsize_t>>>;

The coordinate type on line 1 is now a 3D array rather than a 4D array as it was for structured grids since we are storing now an unstructured list of coordinates. To index these, we only need to know the zone, vertex, and if it is the x or y coordinate. Hence we only need 3 dimensions. The index type set on line 2 is something we saw for the structured grid as well, where it was set to std::array<unsigned, 2>, as we need an index pair of i,j indices to find a vertex in our grid. For unstructured grids, this is going to be a single index to find a specific vertex, no need for a second index.

The interface connectivity and boundary condition information type specified on lines 3-4 are the same as for the structured grid, but this time we store single vertex indices instead of i,j index pairs in these structs, and this is reflected by the template argument to these two structs, which uses the index type defined on line 2.

Lines 6-8 define element connectivity array types. We discussed above that the element connectivity array will hold information about which vertices make up a cell. We differentiate here between internal, interface, and boundary cells. Internal element connectivity information is stored in a 3D array where the indices index the zone, cell, and then the vertices in that cell. Interface and boundary connectivity data has one more dimension, as we need to be able to index for each zone potentially more than one interface or boundary. Thus, the 4D indices go over zones, interfaces/boundaries, cells, and then vertices within that cell.

Lines 9-10 define additional types that we will need later. The way CGNS stores element connectivity is by placing all information in one large 1D array. We need to figure out where in the 1D array we have to start reading our element connectivity for a specific cell type, for example. The type defined on line 10 will be used to store the start index for the internal, interface, and boundary element connectivity into the global 1D element connectivity array per zone.

Function declarations

The functions defined in the ReadUnstructuredMesh are similar to the ones defined in the structured mesh reading class, with a few additions to read the element connectivity data, as well as to return the element connectivity data. The function declarations are listed below again for convenience:

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

  const CoordinateType& getCoordinates() const { return _coordinates; };
  const InternalElementConnectivityType& getInternalCells() const { return _internalElementConnectivity; };
  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:
  void readCellConnectivity();
  StartLocationForElementConnectivityType getStartLocationForElementConnectivity() const;

We can see on lines 3 and 11-13 that we are implementing here our pure virtual functions defined in the base class. We can also see the different getters on lines 5-8, where we now have one additional getter on line 6 which will return the cell structure of internal cells. We did not have that getter for the structured grid case, as it is not needed (cells can be identified based on i,j index pairs alone). We also see that the type is different and because of these differences, we have not made any getter a pure virtual function in the base class.

There are two additional functions now, which are listed in the private section of the class on lines 16-17, i.e. they are only visible to functions within this class. Both these functions will help us read the entire element connectivity arrays for the internal, interface, and boundary cells. This will be discussed below.

Implementing the interface

Now that we have the interface defined, let’s go through the code and see how we implement the various functions defined in the interface.

Constructor

The constructor is pretty much identical to the structured mesh reading class. We receive the filename and path to the CGNS file, which we then directly pass to the base class constructor. The base class will then store this value and open the CGNS file for us (as well as read the number of bases and assert that we only have a single base).

#include "meshReaderLib/include/readUnstructuredMesh.hpp"

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

Implementing the readMesh() function

The readMesh() function is pretty similar to the structured mesh reading class again, with two small changes; first, we call the readZones() function defined in the base class passing the Unstructured zone type as an argument. Secondly, on line 4, after reading the coordinates, we read the element connectivity as an additional step, which we did not have to do for structured grids.

void ReadUnstructuredMesh::readMesh() {
  readZones(CGNS_ENUMV(Unstructured));
  readCoorinates();
  readCellConnectivity();
  readInterfaceConnectivity();
  readBoundaries();
}

Reading coordinates from a CGNS file

Even reading unstructured grid information is very similar to structured grid information. However, I’d argue that coordinate reading in a CGNS file is always easier for unstructured coordinates, as we always receive our coordinates in a 1D array from the CGNS library, which makes it easier to handle on unstructured grids. In the structured coordinate reading case, we had to then reconstruct a 2D array based on the 1D coordinate array received. We do not need to do this for unstructured grids and can simply copy the coordinates into the global coordinate array. Below is the entire code for reading coordinates and we will discuss it below in greater detail.

void ReadUnstructuredMesh::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 = _zoneSizesMax[zone][0];
    cgsize_t minVertices = _zoneSizesMin[zone][0];

    // based on the max vertices, resize the coordinates array.
    _coordinates[zone].resize(maxVertices);
    for (int i = 0; i < maxVertices; ++i) {
      _coordinates[zone][i].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
      std::vector<double> temp(maxVertices);      
      if (cg_coord_read(_fileIndex, 1, zone + 1, coordname, type, &minVertices, &maxVertices, &temp[0]))
        cg_error_exit();

      // store temporary coordinates in final coordinates array
      for (int vertex = 0; vertex < maxVertices; ++vertex) {
        _coordinates[zone][vertex][dimension] = temp[vertex];
      }
    }
  }
}

Lines 2-6 are identical to the structured mesh reading case, we simply define some variables we need later and make space in our global coordinate array by allocation space for each zone. We then loop over each zone on line 8, where we now simply get the minimum and maximum number of coordinates. For structured grids, we had to get the minimum and maximum for both the x and y direction and thus this was an array with 2 entries.

We resize the global coordinates yet again to make space for all vertices and for both the x and y directions on lines 14-17 and then read the coordinates for both the x and y directions on lines 20-34. Specifically, we read the coordinate information on line 23 (e.g. are coordinates stored as single or double precision, and what is the name of the coordinate array) and then use this information to read the coordinates on line 27 into a temporary array (defined on line 26 for the number of coordinates we are going to read). As alluded to above, we then simply copy the content of the temporary array into the global coordinate array on line 32.

There is not much more to it, we discussed the coordinate reading in more depth for the structured grid class, so again, if this was too quick, I’d suggest having a look at the previous article to make sense of the above code should it not be clear.

Getting start location for elements in connectivity array

We saw in the readMesh() function that after reading the coordinates, we read the element connectivity next. Within that function, we make a call early on to the getStartLocationForElementConnectivity(), in-fact, it is one of the first things we do. Thus, I want to discuss this helper function first and then look at the element connectivity reading itself.

Why we need the start location in the first place

So let’s discuss why we need to find a start location in the element connectivity array first. As with pretty much all other CGNS functions (certainly the ones I have used thus far), whenever we request some form of an array, it is returned to us as a 1D array. If we know the precise structure (e.g. a 2D, 3 by 3 array), then we can also allocate memory for that array, pass it into CGNS, and get it returned with the correct values. We saw that when we were reading zone sizes in the base class.

This works as a 1D, 9-element array will allocate the same chunk of memory compared to a 2D, 3 by 3 array (as long as we use C-style arrays, not C++ containers!). So we just made our lives a bit easier by providing the correct dimensions during the zone reading, but we could have also passed in a 1D, 9-element array. Most of the time, though, we simply pass in a 1D allocated array of the correct size, read the information and then afterwards map it to our C++ data structure or container we want to use to permanently store that information.

With the element connectivity, things get a bit more complicated. When we receive the element connectivity, we need to know 2 things; the number of cells we are reading connectivity information for, and the number of vertices per cell, so that we can correctly size the connectivity array. If we think this through, this means that each cell type has to be stored separately and so when we are reading element connectivity information in the first zone (where we have both triangle and quad elements), we have to read element connectivity twice, one per element type (as they have different number of vertices).

Furthermore, interfaces and boundaries also need element connectivity data, and this information is also stored within each zone. Thus, we need to read several element connectivity arrays per zone. Interfaces and boundaries have helper routines that allow to distinguish interface and boundary connectivity data from other element connectivity arrays, but internal cells don’t. So, we need to come up with a way to be able to distinguish each element connectivity array that we read.

We are not just getting the element connectivity, but a bunch of other information as well when we are reading this data. One of the data we are getting is the start and end index of each cell, as if all of the connectivity information was written into one large connectivity array. In the unstructured grid shown above, we saw that for zone 1, we had 8 triangles and 4 quads. So, the element connectivity may be read for the triangles first, in which case we would have start=1 and end=8. Afterwards, we may read the element connectivity for the quad cells, for which we would have start=9 and end=12. Then we also have interfaces and boundaries, all with their own set of start and end indices.

So you see, if we were to put all element connectivity into one large array per zone, we would need to know the start and end location for each element connectivity data. I follow a slightly different approach here, in that I want to keep all if these connectivity information separately. I want to have on for all internal cells, one for interface cells, and one for boundary cells. I may have more than 1 element type in these but I am keeping internal, interface, and boundary cells separate.

Practically speaking, this means that later when I read the actual connectivity data, I need to know if I have just read internal, interface, or boundary cells. The way for me to distinguish between these is by looking at the start location and then based on the start location infer which type I have just read. Therefore, the helper function we look at below will provide me with that information, i.e. a map between the type of element connectivity and its corresponding start location(s) in the global element connectivity array. So let’s look at the code then.

Implementation of reading the element connectivity start locations

Have a look through the code below and then we’ll discuss it afterwards in greater detail.

std::tuple<std::vector<std::vector<cgsize_t>>, std::vector<std::vector<cgsize_t>>, std::vector<std::vector<cgsize_t>>>
  ReadUnstructuredMesh::getStartLocationForElementConnectivity() const {

  std::vector<std::vector<cgsize_t>> internalElementStart(_numberOfZones);
  std::vector<std::vector<cgsize_t>> interfaceElementStart(_numberOfZones);
  std::vector<std::vector<cgsize_t>> boundaryElementStart(_numberOfZones);

  // for each zone, read the starting location of the element connectivity for internal, interface, and boundary cells
  unsigned totalNumberOfInterfaces = 0;
  for (int zone = 0; zone < _numberOfZones; ++zone) {
    // read interface connectivity 
    int numberOfInterfaces = 0;
    if (cg_nconns(_fileIndex, 1, zone + 1, &numberOfInterfaces)) cg_error_exit();
    for (int interface = 0; interface < numberOfInterfaces; ++interface) {
      char interfaceName[128], donorName[128]; GridLocation_t location; GridConnectivityType_t connectType;
      cgsize_t donorData, numberOfPoints, numberOfDonorCells, start; DataType_t donorDataType; ZoneType_t donorZoneType;
      PointSetType_t pointSet, donorPointSet; 

      if (cg_conn_info(_fileIndex, 1, zone + 1, interface + 1, interfaceName, &location, &connectType, &pointSet,
        &numberOfPoints, donorName, &donorZoneType, &donorPointSet, &donorDataType, &numberOfDonorCells))
        cg_error_exit();

      if (cg_conn_read(_fileIndex, 1, zone + 1, interface + 1, &start, donorDataType, &donorData)) cg_error_exit();
      interfaceElementStart[zone].push_back(start);
      totalNumberOfInterfaces++;
    }

    // read boundary connectivity 
    int numberOfBoundaries = 0;
    if (cg_nbocos(_fileIndex, 1, zone + 1, &numberOfBoundaries)) cg_error_exit();
    for (int boundary = 0; boundary < numberOfBoundaries; ++boundary) {
      int normalList; cgsize_t start;
      if (cg_boco_read(_fileIndex, 1, zone + 1, boundary + 1, &start, &normalList)) cg_error_exit();
      boundaryElementStart[zone].push_back(start);
    }

    // read all connectivity. we can't directly read only connectivity of internal cells but instead have to read all of
    // them, which will include the interface and boundary connectivity, as they are all stored under the same node (we 
    // just happen to have separate functions to read boundary and interface connectivity but not for internal cells).
    // Once we have all connectivity read, we simple remove the interface and boundary connectivity from it and are left
    // with internal connectivity information only.
    int numberOfSections = 0;
    if (cg_nsections(_fileIndex, 1, zone + 1, &numberOfSections)) cg_error_exit();
    for (int section = 0; section < numberOfSections; ++section) {
      char sectionName[128]; ElementType_t elementType; cgsize_t start, end; int nbndry, parentFlag;
      if (cg_section_read(_fileIndex, 1, zone + 1, section + 1, sectionName, &elementType, &start, &end, &nbndry,
        &parentFlag)) cg_error_exit();
      internalElementStart[zone].push_back(start);
    }
    // erase all boundary connectivity start locations from internal elements
    for (const auto &boundary : boundaryElementStart[zone]) {
      internalElementStart[zone].erase(std::remove(internalElementStart[zone].begin(), internalElementStart[zone].end(),
        boundary), internalElementStart[zone].end());
    }
    // erase all interface connectivity start locations from internal elements
    for (const auto &interface : interfaceElementStart[zone]) {
      internalElementStart[zone].erase(std::remove(internalElementStart[zone].begin(), internalElementStart[zone].end(),
        interface), internalElementStart[zone].end());
    }
  }
  assert(totalNumberOfInterfaces % 2 == 0 && "Expected an even number of interface pairs");

  return {internalElementStart, interfaceElementStart, boundaryElementStart};
}

Lines 4-9 are a bit of housekeeping, allocating space for each zone and initialising variables. We define here internalElementStart, interfaceElementStart, and boundaryElementStart, which are 2D std::vectors that store the start location into the element connectivity array for each zone. Between lines 10-60, we loop over all zones and then fill these arrays.

On lines 13-26, we read the start location for interface element connectivity data. Line 13 checks if there is any interface at all in the current zone and if so, we loop over all interfaces in that zone on lines 14-26. Line 19 reads a bunch of interface information that we largely don’t care about (for now), we only care about the donorDataType variable, as we need that on line 23 to read the interface data. We get the start location within the global element connectivity array from that call. We store that start location in the interfaceElementStart array and increase the number of interfaces by 1.

Lines 29-35 are doing the same thing, this time just for boundary connectivity data. On line 33, we read the start location for boundary connectivity data and push that into our boundaryElementStart array.

Similarly, we do the same on lines 42-49 but this time for all sections. Sections will contain not just internal cells, but also interfaces and boundaries. Remember, for the 2D grid we saw above, we had a square domain for each zone, meaning we have 4 sides and thus the number of interfaces and boundaries is going to be 4. With two additional element types on the internal domain (for which we need to store separate element connectivity data), this means we are going to get 6 sections here. We read all the start locations, even if they are for interfaces or boundaries, and store that, for the moment, in the internalElementStart array.

The next thing we want to do is to remove both the interface’s and the boundary’s start locations from the internalElementStart array. This is what we are doing on lines 51-54 and 56-59 for the boundary and interface start locations, respectively. We make use here of the erase/remove idiom and this will require you to have a good grasp of C++’s standard template library (STL), which we discussed previously on this site.

In a nutshell, the STL algorithm std::remove will only identify which items should be removed from the container, but it does not have the authority (or knowledge of how) to do so. Removing elements is the responsibility of the container, but it does not have any specific implementation for that. Remember, the goal of the STL is to write algorithms once and use them for all containers, this means we can’t implement them in container classes themselves. Thus, we first identify the items that need to be removed and then erase them from the container.

Line 61 checks if an even number of interfaces was processed and read. If we have an uneven number, we are missing an interface somewhere. Annoyingly, this behaviour is different from reading structured grids, at least the way that Pointwise is writing CGNS files, where structured grids store global interface information and unstructured grids local interface information, i.e. for each zone. Other mesh generators may differ.

On line 63, we return all start locations using a brace initialised list. In this case, the function returns a std::tuple, whose property it is that it can pack different variables together into a list, essentially. It’s an easy way to return more than a single variable from functions, instead of passing them by reference. In this case, we can create and allocate memory within the function (as we did on lines 4-6) and then simply pass the end product out of the function, no need to allocate memory before entering the function.

Reading cell connectivity

Now that we have the start locations read for internal, interface, and boundary element connectivity data, we are able to go ahead and read the element connectivity data and store that in our class variables. The code for that is given below, we’ll look at it line by line below the code.

void ReadUnstructuredMesh::readCellConnectivity() {
  _internalElementConnectivity.resize(_numberOfZones);
  _interfaceElementConnectivity.resize(_numberOfZones);
  _boundaryElementConnectivity.resize(_numberOfZones);

  // the element connectivity contains the internal, interface, and boundary cells. First, we need to find out what the
  // starting location is for each cell type so that we no which part of the connectivy array to start reading from to
  // get internal, interface, or boudnary cell connectivity.
  auto [internalElementStart, interfaceElementStart, boundaryElementStart] = getStartLocationForElementConnectivity();

  // go through all zones and read the connectivity
  for (int zone = 0; zone < _numberOfZones; ++zone) {
    int numberOfSections;
    if (cg_nsections(_fileIndex, 1, zone + 1, &numberOfSections)) cg_error_exit();

    _interfaceElementConnectivity[zone].resize(interfaceElementStart[zone].size());
    _boundaryElementConnectivity[zone].resize(boundaryElementStart[zone].size());
    
    unsigned internalElementIndex = 0;
    unsigned interfaceElementIndex = 0;
    unsigned boundaryElementIndex = 0;

    // for each zone, go through all sections and read the element connectivity information
    for (int section = 0; section < numberOfSections; ++section) {
      char sectionName[128]; ElementType_t elementType; cgsize_t start, end, parentData, elementDataSize;
      int boundaryFlag, parentFlag, numberOfVerticesPerElement;

      // this routine provides us with the variable 'start', which tells us where to start reading element connectivity
      // from. We use this value later together with the above received start locations for internal, interface, and
      // boundary start locations to figure out what element connectivity we have just read
      if (cg_section_read(_fileIndex, 1, zone + 1, section + 1, sectionName, &elementType, &start, &end, &boundaryFlag,
        &parentFlag)) cg_error_exit();

      // the number of elements in the element connectivity array
      if (cg_ElementDataSize(_fileIndex, 1, zone + 1, section + 1, &elementDataSize)) cg_error_exit();

      // this temporary array will hold all element connectivity information in this section
      std::vector<cgsize_t> elements1D(elementDataSize);
      if (cg_elements_read(_fileIndex, 1, zone + 1, section + 1, &elements1D[0], &parentData)) cg_error_exit();
      if (cg_npe(elementType, &numberOfVerticesPerElement)) cg_error_exit();

      assert(elementDataSize % numberOfVerticesPerElement == 0 &&
        "The number of elements in the section is not a multiple of the number of vertices per element");
      
      unsigned numberOfElements = elementDataSize / numberOfVerticesPerElement;

      // the element connectivity is received from CGNS in a 1D array. It has a total size of
      // size = [numberOfElements * numberOfVerticesPerElement]. What we want, for easier indexing, is a 2D array where
      // each dimension will have a size of [numberOfElements][numberOfVerticesPerElement]. This lambda expression does
      // that for us. We use a lambda here because we want to use this routine 3 time (interal, interface, and boundary)
      // cells. 
      auto map1DElementsto2DVector = [numberOfVerticesPerElement, numberOfElements, &elements1D]
        (std::vector<ReadUnstructuredMesh::IndexType> &target) {
        unsigned counter = 0;
        unsigned elementStartLocation = target.size();
        target.resize(elementStartLocation + numberOfElements);
        for (int element = 0; element < numberOfElements; ++element) {
          unsigned elementIndex = elementStartLocation + element;
          target[elementIndex].resize(numberOfVerticesPerElement);
          for (int vertex = 0; vertex < numberOfVerticesPerElement; ++vertex) {
            target[elementIndex][vertex] = elements1D[counter++] - 1;
          }
        }
      };

      auto &internal = internalElementStart[zone];
      auto &interface = interfaceElementStart[zone];
      auto &boundary = boundaryElementStart[zone];

      // now we make use of the variable 'start' and the starting location for the internal, interface, and boundary
      // cells. We go through all of them and check if they contain the value for 'start'. If so, then we know that
      // we have read the connectivity for the internal, interface, or boundary cells. We then map the 1D connectivity
      // array to the 2D element connectivity array using our above defined lambda expression.
      if (std::find(internal.begin(), internal.end(), start) != internal.end()) {
        map1DElementsto2DVector(_internalElementConnectivity[zone]);
        internalElementIndex++;
      } else if (std::find(interface.begin(), interface.end(), start) != interface.end()) {
        map1DElementsto2DVector(_interfaceElementConnectivity[zone][interfaceElementIndex]);
        interfaceElementIndex++;
      } else if (std::find(boundary.begin(), boundary.end(), start) != boundary.end()) {
        map1DElementsto2DVector(_boundaryElementConnectivity[zone][boundaryElementIndex]);
        boundaryElementIndex++;
      } else {
        throw std::runtime_error("Could not itendify which element connectivity was read");
      }
    }
    int totalSectionsProcessed = internalElementIndex + interfaceElementIndex + boundaryElementIndex;
    assert(numberOfSections == totalSectionsProcessed && "Did not process all sections, some information is lost");
  }
}

Setting up element connectivity reading

Lines 2-4 allocate sufficient memory for our element connectivity arrays and then on line 9 we receive the starting locations for internal, interface, and boundary element connectivity from the getStartLocationForElementConnectivity(), which we discussed just before this section. You can see that we need to use the auto [] syntax here to receive the elements of a std::tuple returned by this function, where we have as many entries within the square brackets as there are variables returned by the tuple. With this information available, we start looping over all zones on line 12.

On line 14, we see that we call cg_nsections(), which we also called in getStartLocationForElementConnectivity(). There will be an overlap of functions we have to call now as trying to get the start location for each element connectivity data will use the same CGNS mid-level library calls.

On lines 16-17, we allocate memory for the interface and boundary element connectivity data, but not for the internal cells. We may have several interfaces and/or boundaries per zone and we want to store the element connectivity data for each interface/boundary separately, hence we need to make space for as many entries as we have in the start location array. Say we have 3 entries in the start location array for the boundary element connectivity data, then we know that this corresponds to 3 different boundary conditions and so we want to create element connectivity data for each separately.

Internal cells only have a single-element connectivity array for each zone, which may contain more than one cell type, so we don’t need to allocate any memory here for the connectivity data. If you look back at how the internal element connectivity data was defined in the class header/interface, you’ll see that it has one dimension less than the interface and boundary element connectivity data.

Reading element connectivity data for internal, interface, and boundary cells

On line 24, we start looping over all sections and remember that each section may store element connectivity data for the internal, interface, or boundary cells. We define some required CGNS variables on lines 25-26 and then start reading the information for the current section on lines 31-32. This call has the start value into the global element connectivity array, so we can compare this value later against the start locations we received on line 9 to see which type of element connectivity we have just read.

Line 35 returns to us the total number of entries in the 1D element connectivity array that we are about to read. We can use this to allocate memory for a temporary array into which we will read the element connectivity array of the current section, which we do on line 38.

On line 39, we then read the element connectivity data into our temporary array and we follow up this call by cg_npe(), which tells us how many vertices there are per cell for the current element. For example, for a triangle, we get 3 and for a quad element, we get 4. If we are on the boundaries or an interface, we get 2 vertices per cell, which is now of type bar (i.e. a line/edge).

Line 42 is a sanity check, we make sure that the size of the 1D element connectivity we read into the temporary array is divisible by the number of vertices per cell. If it is not, something went wrong during mesh generation and export. If this test passes, we can calculate the number of elements using line 45, i.e. take the size of the element connectivity data and divide it by the number of vertices. Since we checked that this is divisible and returns a modulus of 0, we can be assured that this integer division will not spring up any nasty surprises.

A lambda expression to map 1D element connectivity data to 2D

Lines 52-64 require a bit more explanation, so let’s copy the code and put it below so that we can easily refer to it. In the code below, we are defining a lambda expression and this is a good use case for why you want to use them. We looked at lambda expressions before, so you may want to go back to that article first before continuing here if you don’t feel comfortable with them yet.

auto map1DElementsto2DVector = [numberOfVerticesPerElement, numberOfElements, &elements1D]
(std::vector<ReadUnstructuredMesh::IndexType> &target) {
  unsigned counter = 0;
  unsigned elementStartLocation = target.size();
  target.resize(elementStartLocation + numberOfElements);
  for (int element = 0; element < numberOfElements; ++element) {
    unsigned elementIndex = elementStartLocation + element;
    target[elementIndex].resize(numberOfVerticesPerElement);
    for (int vertex = 0; vertex < numberOfVerticesPerElement; ++vertex) {
      target[elementIndex][vertex] = elements1D[counter++] - 1;
    }
  }
};

We create a lambda expression here because we intend to use the same function (lambda expression) for the internal, interface, and boundary cells. The job of the map1DElementsto2DVector[]() lambda expression is to take the 1D element connectivity that we read and then put it into either the internal, interface, or boundary element connectivity array. However, we don’t just want to put the 1D array into our class variables, but rather transform the element connectivity data into a 2D array, so that we can loop over all cells and then query the vertices for that cell easily.

To do so, we first need to capture a few things in the lambda capture block, i.e. between the square brackets []. We need the number of vertices per element, the number of elements, and then the temporary array that holds the 1D connectivity data. Capturing this in the lambda means we don’t have to pass it as an argument. So the only argument we are passing is the internal, interface, or boundary element connectivity array, which we call target here.

On line 4, we define an element start location. This will read the size of either the internal, interface, or boundary element connectivity data we pass into this function. For interfaces and boundaries, this will always be zero, but not necessarily for internal cells. We said that for internal cells we want to store all different types of element types in the same connectivity array.

Thus, say we have already stored element connectivity data for triangles in our internal element connectivity array, by the time we get into this function and want to write quad element connectivity data, we need to know how much data is already in our internal element connectivity data (which is called target here). We then use that size and the number of elements that we calculated before and captured in this lambda to resize our target array (i.e. to allocate additional memory on top of what was already allocated, if anything).

Then, we loop over the number of elements (for which we just made space) on line 6, the following line 7 defines an element index which takes the offset into account that we discussed above (we are skipping over elements that we have already written to the target array), allocate enough memory to hold as many indices as we have vertices for the current element type, and then loop over the number of vertices per element.

Line 10 allows us, then, to store the 1D element connectivity data in the current target as a 2D array. We use the counter variable here defined previously on line 3 to monotonically increase the index we read. We use a little trick here in that we specify the index as counter++. This means that we read the value of counter (which is going to be 0 the first time we read it) and then we increment it by one with the ++ syntax. The next time we read the value for counter, it is going to be 1, after which we increment it again, and the cycle continues.

Since CGNS uses Fortran-based indexing where the first index starts at 1, we have to subtract 1 from all of the element connectivity information we write into our internal, interface, or boundary element connectivity arrays.

Writing element connectivity data into class variables

Returning to the code printed above the lambda expression, we have arrived at lines 66-68, which are just simple helper variables that will help to reduce the code we have to write (i.e. our lines will get shorter using a variable name with fewer characters, that’s all that we are trying to achieve here).

But let’s look at lines 70-85, which I have copied below again for convenience:

// now we make use of the variable 'start' and the starting location for the internal, interface, and boundary
// cells. We go through all of them and check if they contain the value for 'start'. If so, then we know that
// we have read the connectivity for the internal, interface, or boundary cells. We then map the 1D connectivity
// array to the 2D element connectivity array using our above defined lambda expression.
if (std::find(internal.begin(), internal.end(), start) != internal.end()) {
  map1DElementsto2DVector(_internalElementConnectivity[zone]);
  internalElementIndex++;
} else if (std::find(interface.begin(), interface.end(), start) != interface.end()) {
  map1DElementsto2DVector(_interfaceElementConnectivity[zone][interfaceElementIndex]);
  interfaceElementIndex++;
} else if (std::find(boundary.begin(), boundary.end(), start) != boundary.end()) {
  map1DElementsto2DVector(_boundaryElementConnectivity[zone][boundaryElementIndex]);
  boundaryElementIndex++;
} else {
  throw std::runtime_error("Could not itendify which element connectivity was read");
}

We define here an if/else statement and the goal of this if/else statement is to go through all of our element start location arrays (which we call now internal, interface, and boundary on lines 66-68 before) and check if the start location that we read for the current section (some time ago on line 31) does occur in any of the start location arrays we read at the beginning of the function on line 9. If there is a match, we then know which type we have (i.e. element connectivity data for internal, interface, or boundary cells). If there is no match, something went wrong and we terminate with a runtime exception.

The way that we look for the start location we read before is with the std::find algorithm. Since we defined our start location arrays as std::vectors, we can now go through them and find specific entries. The inputs to the std::find algorithm are the start and end iterator of the container, typically obtained with the begin() and end() function if the container has a directional iterator type defined. The third argument is the value we are trying to find, in the case of the start location that we read before.

We look for the start location in all three arrays, i.e. the internal, interface, and boundary start location arrays (lines 5, 8, and 11 above). If we have found the start value in either of them, then the std::find algorithm will return the iterator to the starting location array, which we don’t actually care about. However, if we don’t find the start value in any of these arrays, then the std::find algorithm will return the iterator which is one unit past the last value. This is what the end() function returns as well.

So, all we have to do is use the std::find algorithm and check if the value found is not equal to the last iterator returned by the end() function. If the start value is found in any of the arrays, we return the iterator to that value, which is going to be somewhere within the array and so by definition it is not going to be the same value as the iterator returned by end(). Therefore, if it is not the end() iterator, we know that the start location is in the current array and the if statement will be evaluated as true.

Within either of the if/else statements, we then go ahead and call the lambda expression defined before which will then write the current temporary 1D element connectivity data into the array passed into the lambda expression. You may notice that for internal cells we only pass _internalElementConnectivity[zone] whereas for interfaces and boundaries, we pass _interfaceElementConnectivity[zone][interfaceElementIndex] and _boundaryElementConnectivity[zone][boundaryElementIndex] to the lambda expression.

This has to do, again, with internal cells only storing a single connectivity array per zone, whereas interface and boundary element connectivity data is stored separately for each interface and boundary for each zone. The interfaceElementIndex and boundaryElementIndex were initially set to 0 and we are using a similar trick to what we did with the counter variable in the lambda, where we are just increasing this index every time we visit this if/else statement.

Incidentally, if it is the first time you are hearing about iterators, they allow us to simply loop over data in a data structure, like a std::vector in the above-discussed case. We looked at iterators, containers, and algorithms in my article about the C++ standard template library, or the STL in short, have a look if you need a refresher.

To finish off this section, in the original code posted at the top of this section, on lines 87-88 we check that once we have processed all element connectivity data, we have visited the internal, interface, and boundary element connectivity writing routines as many times as there are sections, or, with other words, that we have not missed reading any data. Should there have been any issues before, we likely would have thrown an exception already in the if/else statement before, but this is one more check we can do to save us from unwanted behaviour.

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

Reading interfaces for unstructured grids is slightly different to reading interfaces on structured grids, or at least if using Pointwise. However, I’d argue that from an educational point of view, this is a good thing. We can see how to read interface information in two different ways and should you use a different mesh generator, you may have either of these forms and should be able to adapt your code. From an end-user point of view, though, it would have been nice to have a consistent interface reading implementation, but there may be good reasons why Pointwise had to implement it this way that I am unaware of.

When we were reading the interface connectivity information on the structured grid, we were reading global interface connectivity, meaning it is stored under the CGNS base node and not for each zone. If you need a refresher, have a look at the structured mesh reading article. For unstructured meshes, Pointwise does store the information under each zone so we can simply loop over all zones and then read the interfaces as we would read boundary conditions. This is shown in the code below, which we will, as usual, discuss afterwards.

void ReadUnstructuredMesh::readInterfaceConnectivity() {
  char interfaceName[128], donorName[128]; GridLocation_t location; GridConnectivityType_t connectType;
  PointSetType_t pointSet, donorPointSet; cgsize_t numberOfPoints, numberOfDonorCells;
  ZoneType_t donorZoneType; DataType_t donorDataType;

  // create map that stores, for each interface (identified by name), the zones that it connects
  std::unordered_map<std::string, std::vector<std::array<int, 2>>> interfaceNameToZone;
  for (int zone = 0; zone < _numberOfZones; ++zone) {
    int numberOfInterfaces;
    if (cg_nconns(_fileIndex, 1, zone + 1, &numberOfInterfaces)) cg_error_exit();
    for (int interface = 0; interface < numberOfInterfaces; ++interface) {
      if (cg_conn_info(_fileIndex, 1, zone + 1, interface + 1, interfaceName, &location, &connectType, &pointSet,
        &numberOfPoints, donorName, &donorZoneType, &donorPointSet, &donorDataType, &numberOfDonorCells))
        cg_error_exit();
      
      if (interfaceNameToZone.find(interfaceName) == interfaceNameToZone.end())
        interfaceNameToZone.emplace(interfaceName, std::vector<std::array<int, 2>>{{zone, interface}});
      else
        interfaceNameToZone[interfaceName].push_back({zone, interface});
    }
  }

  // read interface connectivity for each zone
  _interfaceConnectivity.resize(_numberOfZones);
  for (int zone = 0; zone < _numberOfZones; ++zone) {    
    int numberOfInterfaces;
    if (cg_nconns(_fileIndex, 1, zone + 1, &numberOfInterfaces)) cg_error_exit();
    assert(numberOfInterfaces == _interfaceElementConnectivity[zone].size() && "The number of interfaces should be the"
      " same as the number of interfaces for which elements have been read");

    // for each zone, go through all interfaces stored in this zone and store the connectivity information
    _interfaceConnectivity[zone].resize(numberOfInterfaces);
    for (int interface = 0; interface < numberOfInterfaces; ++interface) {
      if (cg_conn_info(_fileIndex, 1, zone + 1, interface + 1, interfaceName, &location, &connectType, &pointSet,
        &numberOfPoints, donorName, &donorZoneType, &donorPointSet, &donorDataType, &numberOfDonorCells))
        cg_error_exit();
        
      assert(donorZoneType == CGNS_ENUMV(Unstructured) && "Expected unstructured interface");

      auto zones = interfaceNameToZone[interfaceName];
      assert(zones.size() == 2 && "Expected 2 zones per interface");

      const auto &zoneOwner = zones[0][0];
      const auto &interfaceOwner = zones[0][1];
      const auto &zoneNeighbour = zones[1][0];
      const auto &interfaceNeighbour = zones[1][1];

      _interfaceConnectivity[zone][interface].zones = {zoneOwner, zoneNeighbour};
      _interfaceConnectivity[zone][interface].ownerIndex = _interfaceElementConnectivity[zoneOwner][interfaceOwner];
      _interfaceConnectivity[zone][interface].neighbourIndex =
        _interfaceElementConnectivity[zoneNeighbour][interfaceNeighbour];
    }
  }
}

Creating a zone index lookup table

As we have seen before, CGNS is not always providing us with all the information we need. I’m not sure if this is by choice or if the steering committee did overlook some vital information its functions should be returning. It seems odd, given that the CGNS routines return so many different variables, most of which don’t even have data, but are there just in case someone wants to implement and consume it.

When it comes to an interface, the most important part that we want to know is which two zones are connected. Without that information, we don’t have an interface, yet the CGNS standard does not require it. This is why we have to generate yet another look-up table that we will use to store all zones that are connected to the same interface, identified by the interface name. There should only ever be 2, and exactly 2, zones attached to the same interface.

This is what is done between lines 7-21. On line 7, we define a std::unordered_map called interfaceNameToZone, which will be our lookup table, i.e. where we store, for each interface name, all zones that attach to it. We have to store to information, the zone index, and the interface index. We will use these later to uniquely identify interfaces with a zone that belong connect to interfaces on a different zone. We store this informationin a std::vector, which itself stores a std::array with two entries (one for the zone and one for the interface index).

Line 8 loops over all zones, and subsequently we check the total number of interfaces for that zone and loop over all interfaces in that zone on line 11.

On line 12, we call cg_conn_info(), which will provide us with the interface name, i.e. the variable interfaceName. We know the zone that is attached to this interface from our loop variable zone defined in the for loop on line 8.

Line 16 checks if we already have an entry in the lookup table, i.e. the interfaceNameToZone variable. We are using the same std::find() algorithm here discussed above in the element connectivity reading and we are checking again if the algorithm returns the end() iterator, meaning that the interface name is not yet in our lookup table. If that is the case, we create a new entry for it on line 17, where we tie the name of the interface with the zone and interface that we are currently in.

We initialise a std::vector as the second argument and we add one entry using the brace initialiser list. Since the entry into the std::vector is itself a std::array, which also needs to be brace initialised, we have the double brace syntax, i.e. {{zone, interface}}. So the outer braces contain one element for the std::vector, and the inner braces contain two elements for the std::array.

If we have written the interface name to the lookup table and we loop over the remaining zones and interfaces, eventually we will reach the other side of the interface in a different zone. Once that is the case, the interfaceName variable read on line 12 in the cg_conn_info() function is going to be the same as for the interface information we were writing earlier.

When we now get to line 16, we check if the interface name is in our lookup variable, which it is, so the iterator returned will not be equal to the last iterator returned by the end() function. Thus, we jump straight to line 18 and execute the else statement, where we simply add the current zone and interface to the std::vector, which should now hold two std::arrays with the different zone and interface indices.

Reading interface connectivity data

To be fair, quite a lot of the heavy lifting has already been done. We have the lookup table and know which interface name connects to which zones, and we have read the interface element connectivity data, so the remaining steps are rather straightforward with that information at hand. This is done in the code above between lines 24-53. The _interfaceConnectivity variable will contain all of the interface data so we first allocate sufficient memory for all zones, and then find the number of interfaces, allocate memory for each interface per zone, and then loop over all interfaces in the current zone on line 33.

It is perhaps worth pointing out that we do some sanity checks here on line 28, i.e. we are testing the number of element connectivities we have available. Remember, we are storing this information for each interface separately, so if we have 2 interfaces in one zone, we expect 2 separate element connectivity arrays.

Between lines 33-52, we are now reading the interface data and storing the already processed interface data in this part. For simplicity, I have copied the code again below so we don’t have to scroll up and down:

for (int interface = 0; interface < numberOfInterfaces; ++interface) {
  if (cg_conn_info(_fileIndex, 1, zone + 1, interface + 1, interfaceName, &location, &connectType, &pointSet,
    &numberOfPoints, donorName, &donorZoneType, &donorPointSet, &donorDataType, &numberOfDonorCells))
    cg_error_exit();
    
  assert(donorZoneType == CGNS_ENUMV(Unstructured) && "Expected unstructured interface");

  auto zones = interfaceNameToZone[interfaceName];
  assert(zones.size() == 2 && "Expected 2 zones per interface");

  const auto &zoneOwner = zones[0][0];
  const auto &interfaceOwner = zones[0][1];
  const auto &zoneNeighbour = zones[1][0];
  const auto &interfaceNeighbour = zones[1][1];

  _interfaceConnectivity[zone][interface].zones = {zoneOwner, zoneNeighbour};
  _interfaceConnectivity[zone][interface].ownerIndex = _interfaceElementConnectivity[zoneOwner][interfaceOwner];
  _interfaceConnectivity[zone][interface].neighbourIndex =
    _interfaceElementConnectivity[zoneNeighbour][interfaceNeighbour];
}

We loop over all interfaces and then read the information for these interfaces on line 2, as we have done previously when creating the lookup table based on the interface name. Now, we use the donorZoneType, or the neighbouring zone type, to check that the interface data is written as an unstructured grid. It is possible to connect unstructured and structured grids at an interface (and certainly our grid arrangement would allow for that in our small test grid), but this would require a different data structure, essentially a combination of the structured and unstructured mesh reading class. While it is possible, it is uncommon for CFD solvers to do this.

Based on the interface name, we retrieve the owner and neighbour zone and interface indices on line 8 from the lookup table we created earlier. Now we check that we do have indeed only 2 entries for each side of the interface, otherwise it would be undefined. We then go ahead and create some more telling variable names on lines 11-14 and all that’s left is to write all available information into our interface connectivity information array.

Line 16 stores both the owner and neighbouring zone index and for both the owner (line 17) and its neighbour (line 18), we store the element connectivity information. We use the zone and interface index from the lookup table, as well as the element connectivity data that we read in the previous section. With all of that information written to the array, we have finished creating the interface information and can move on to read boundaries, the last remaining step.

Reading boundary conditions from a CGNS file

The boundary condition reading is pretty similar to the structured mesh reading. I’d argue it is actually simpler in this case, as we have already read the boundary element connectivity data. For the structured grid, we received the first and last index as a point range and then had to construct the indices for i and j. The code is given below with explanations to follow.

void ReadUnstructuredMesh::readBoundaries() {
  // go through all zones and check which boundary conditions are stored under each zone
  _boundaryConditions.resize(_numberOfZones);
  for (int zone = 0; zone < _numberOfZones; ++zone) {
    int numberOfBoundaries;
    if (cg_nbocos(_fileIndex, 1, zone + 1, &numberOfBoundaries)) cg_error_exit();
    _boundaryConditions[zone].resize(numberOfBoundaries);

    // for each zone, go through all boundaries stored in this zone and store the boundary information
    for (int boundary = 0; boundary < numberOfBoundaries; ++boundary) {
      char boundaryName[128]; int normalIndex, numberOfDataset; PointSetType_t pointSetType;
      cgsize_t numberOfPointsInBoundary, normalListSize; DataType_t normalDataType; BCType_t boundaryType;

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

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

      _boundaryConditions[zone][boundary].indices = _boundaryElementConnectivity[zone][boundary];

      // Writing the boundary name and type to the boundary structure. We need to check if we need to get this
      // information from a family node or if we can read this directly from the boundary node.
      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);
      }
    }
  }
}

After allocating memory for all zones and then looping over them on lines 3-4, we query the number of boundary conditions on line 6 and then allocate memory for all boundary conditions on the following line. Between lines 10-42, we loop over all boundaries to read their type, name, and write all of this information, together with the element connectivity data read before into our boundary condition array.

On line 14, we receive the boundary name and its type. We check on lines 17-18 that the points received by the CGNS call is of type PointRange, meaning we should only receive the start and end location. This is pretty inconsequential, though, as we are not actually proceeding to read this information. If you encounter issues at this point, you should be able to delete these two lines without consequence. What we do instead is to store the boundary element connectivity data that we already read on line 20. Our boundary condition array now has all the elements available on that boundary.

We check on line 24 whether the boundary type is set to FamilySpecified or not. If it is, then we need to go ahead and read the boundary name and type from the family conditions that we created in the base class. The way to retrieve the information requires pointing the CGNS file to the family node, which we do on line 26 with the cg_goto() function (if we think of nodes in a CGNS file as folders on a hard drive, then cg_goto() is equivalent to changing directories).

Pointing at the right node, we can read the family name on line 32 which we use to look up the information we need from our family lookup table (again, we created that in the base class, see the above-linked article for a description). First, we store the family name, which now is the name of the boundary condition, on line 35 and then proceed to look up the boundary condition stored in the boundary lookup table on line 36 and store its result in our boundary condition array.

If the boundary condition that we read earlier on line 14 is not of type FamilySpecified, then we have read the actual type and do not need to retrieve this information from a family node. In this case, we can directly store the boundary name and its type that we read on line 14 in our boundary condition array, which we do on lines 39-40.

If you compare this implementation with that of the boundary condition reading on the structured mesh, then you will that there are quite a few parallels as discussed above. I went through this description a bit faster as we have already covered most of the functionality in the aforementioned article. If you feel you need a bit more explanation, have a look at that article as well for additional explanations.

Summary

If you made it up to this point, and you read the previous 2 articles, you should have a firm grasp on how to read structured and unstructured grids. I have tried to be as general as possible (i.e. support multi-block/zone mesh reading for both structured and unstructured grids) and chances are you will be able to read pretty much any grid with this implementation.

CGNS is a complex library, allowing you to write the same information in many different ways. Most paths are standardised but some are left to the developers implementing the CGNS routines into their software. I have been using grids generated by Pointwise and thus you should be able to proceed with Pointwise-generated grids.

If you use a different software, there should be enough fail-fast assert statements in place to point you in the right direction. Check which assert is failing and then read through the code to see which modifications you have to do. If there are modifications that I was able to foresee, I mentioned those in the text, so check the writeup here again and see if the solution is not already presented here.

Even though we have limited ourselves to 2D grids only, going to 3D is not that much more complicated. The difference is mainly in the memory allocation and looping over data. If you understand the entire code up until this point, you should be able to extend this to 3D.

This finalised the development of this library. What is left to do for us is to turn this code into an actual library, i.e. compile it as a static or dynamic library, and, we hinted at that in earlier articles as well, create some form of tests that we can run to ensure the mesh reading was successful. This is what we deal with in the next article before closing this series.

---

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.