In this article, we will look at implementing a base class from which we will derive all future mesh reading classes (for structured and unstructured grids). This implementation will already contain some functionality to interact with the CGNS file and read some basic information, so we will get a feel for how to work with CGNS files and what functions to call to read the required information. Along with the base class, we also set up the project structure, including the build scripts for Windows and UNIX, as well as a first implementation of the main.cpp file, which will serve as a test file to read our different types of grids.
By the end of this article, you should have a good grasp of the class structure required to read either a structured or unstructured mesh file and through the exposure of the CGNS mid-level library calls that we will carry out, you should also feel comfortable reading basic information from a CGNS file such as reading the number of bases, zones, families, as well as reading boundary conditions directly from he family nodes.
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-2.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
Introduction
In the previous two articles, I introduced you to the general structure of a CGNS file, what it can (and currently still) can’t do, and I also laid out my motivation for why we ought to learn this format if we are serious about CFD development. We used the utility cgnslist, which is compiled as part of the CGNS library and sits among various other CGNS tools that we haven’t looked at for now. With cgnslist, we were able to inspect 4 different structured and unstructured grids that we want to read with our mesh reading library that we will create as part of this series.
In today’s article, we will look at the basic library structure, the build scripts for Windows and UNIX (Linux, macOS), and then develop a base class from which we will derive a structured and unstructured mesh reading version. Both these classes will be developed in the next two articles. We are also going to start to write some code to test our library, which we will put, for now, in the main.cpp file. There is a more formal process for testing software and this will be introduced in a different series. For now, we just use this file to test and make sure that we are happy that these files have been read correctly.
Project structure
Let’s start with a high-level overview of our project/library structure. We will develop an initial structure in this article, and then build upon this in the next articles that follow. The project structure that we will be working with looks like the one shown below (and you can download this initial project using the zip file linked at the beginning of this article):
meshReadingLibRoot
├── mesh
│ ├── structured2D.cgns
│ ├── structured2DNoFamily.cgns
│ ├── unstructured2D.cgns
│ └── unstructured2DNoFamily.cgns
├── meshReaderLib
│ ├── include
│ │ ├── readMeshBase.hpp
│ │ └── types.hpp
│ ├── src
│ │ └── readMeshBase.cpp
│ └── meshReader.hpp
├── buildAndRun.ps1
├── buildAndRun.sh
└── main.cppWe have already looked in detail at the mesh/ folder and its content in the previous article. The meshReaderLib/ contains the actual library that we want to develop. As usual, we split the interface/API (application programming interface) and the implementation into two separate folders. Within the meshReaderLib/include/ folder, we collect all header files that tell us what functions need to be implemented and also single to a user of this library which functions to call. We then implement these functions in files stored within the meshReaderLib/src/ folder.
For the moment, there are only the readMeshBase.hpp (header) and readMeshBase.cpp (source) files, but as stated above, we will fill these folders with structured and unstructured mesh reading classes as we go along in this series. You’ll also notice a file called types.hpp, which is located in the include/ folder. This is a personal preference; I like to have a dedicated place in my libraries where I can define custom types and, in particular, enums. We will start populating it today and then use it throughout this series.
When we want to use our library, we need to first compile it as a static/dynamic library (we will develop the build scripts for that towards the end of this series). If you are unsure what the difference is between static/dynamic libraries, we looked at that in much greater depth before in my article on static, dynamic, and header-only libraries. The next thing we need to do is to provide an interface to the library. This is what the meshReaderLib/meshReader.hpp file is for; we include this file in the code where we want to use the mesh reading library, and it will provide all the information to the compiler to use this library. The linker will then get the remaining information from the compiled libraries.
The two build scripts that we will use for now are called buildAndRun.ps1 (Windows) and buildAndRun.sh (UNIX). If you are an astute observer, you’ll notice that we used batch files on Windows before (i.e. buildAndRun.bat), see my previous series on how to build a linear algebra library – basic project structure. Batch files are pretty much outdated by now so we will be using the more up-to-date PowerShell scripting language, and those files end in *.ps1. The build script we developed before should still run as a PowerShell script, it is just the file extension that has changed.
As alluded to above, the main.cpp file, for now, will just be a dumping ground to test that the library reading was successful and that no error messages were encountered during the mesh reading. We will extend this function a bit more in the final part of this series, where we will do some more serious testing to ensure our implementation is actually doing what it is supposed to.
This website exists to create a community of like-minded CFD enthusiasts and I’d love to start a discussion with you. If you would like to be part of it, sign up using the link below and you will receive my OpenFOAM quick reference guide, as well as my guide on Tools every CFD developer needs for free.
Join now
Cross-platform Build scripts
Let’s have a look then at the build scripts we will be using to compile this library on Windows or UNIX (Linux, macOS). We will create an initial structure for these build scripts and then update them as we add more files to our library.
Compiling on Windows
The build script for Windows (buildAndRun.ps1) is shown below:
# 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 .\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 /Fe".\build\cgnsTest.exe" /link /MACHINE:x64 /LIBPATH:"C:\libraries\lib" cgns.lib hdf5.lib msvcrt.lib libcmt.lib
.\build\cgnsTest.exeThe first 5 lines prepare the project. First, we remove the build/ folder and then afterwards re-create it (we just want to get rid of any previously created object files which we dump inside the build/ folder). We then copy the hdf5.dll and zlib.dll library inside the build/ folder. As we discussed in our overview on static, dynamic, and header-only libraries (if you haven’t checked it out, it would be good to go through this first), DLL stands for dynamic-linked libraries and is Windows terminology for dynamic libraries. Why do we need to do this?
The way dynamic libraries are located when an executable is launched is platform-specific and, well, Windows doesn’t have a sensible way of dealing with them. On Windows, you need to dump all of your dynamic libraries in the same folder where you are planning to keep your executable. If you are on Windows right now, go to C:\Program Files\ and pick any folder (ideally, a somewhat larger application). Then, find the executable of that application and look at files ending with *.dll. Just as an example, looking at Blender on my PC, this is what I can see:
I count 21 different *.dll files, and the executable is blender.exe (file extension not shown). For the executable to be able to find all of these dependencies (libraries), they need to be next to the executable. This gives us a horrible folder structure, but, well, this is the price to pay to develop on or for Windows. (Yes, I know we can use environmental variables to specify search paths, but that approach is rarely used. We don’t have an easy way to set the library search path as easily and cleanly as on UNIX, so dll-dumping is a common exercise on Windows).
Thus, we, too have to copy our dependencies into the folder where we keep our executable, which is the build/ folder. I am assuming here that libraries are installed in C:\libraries, so that we can find C:\libraries\bin\hdf5.dll and C:\libraries\bin\zlib.dll. If you have them installed at a different location, you will need to change the copy location. We don’t have to copy the CGNS library because we will use the static version of the library (for no particular reason, but I thought it would be good to show how to deal with both static and dynamic libraries). Static libraries can be located anywhere on the disk.
With all the required dependencies out of the way, we can go ahead and compile the readMeshBase.cpp and main.cpp files on lines 8-9. We are using the C++ 2020 standard, as we will make use of the filesystem header, which was only introduced in 2020. Any compiler modern compiler should offer you C++ 2020 support. We include the current (.) directory, as well as C:\libraries\include, i.e. where we expect header files that define the interface for the various libraries we use (primarily CGNS, but it, in turn, requires knowledge of HDF5 as well).
We define a preprocessor directive called COMPILELIB, this is again a requirement on Windows (we can call it whatever we want) if we are planning to compile our library into a dynamic library. I discussed this at great length in my article on How to write a CFD library: Compiling and testing which you may want to check out. It describes the purpose of this pre-processor directive and how this is used in the code later. We will come across it later when we look at the base class implementation as well.
The rest is pretty much straightforward, we use /O2 as a compiler optimisation, but we could have used /O0 here with the /Zi flag for debugging purposes as well (and during development, this is advisable). In this case, we simply compile the library with optimisation already turned on.
On line 12, we link both object files (meshReaderBase.obj and main.obj) into the executable cgnsTest.exe, which is located inside the build/ folder. We specify that we want to compile a 64-bit version (because, apparently, Microsoft thinks that this is necessary) and tell the linker where to find additional (static) libraries. In this case, in C:\libraries\lib. In this folder, the linker will find the cgns.lib file, which it will then copy into our executable. Confusingly on Windows, we also have to link against hdf5.lib, despite already linking against the *.dll file.
Here, the *.lib file will simply provide some additional information but no actual implementation, that is provided by the *.dll file. Both hdf5.lib and hdf5.dll are needed, but not zlib.lib, oh no, we can’t have consistency on Windows! We also throw in some Windows specific libraries (msvcrt.lib and libcmt.lib, god knows why) and this should produce then an executable we can use. If you wonder how to figure out which libraries you have to link against, it is a process of trial and error. You do the best you can to second guess what Windows need and then you have to plow through the compiler/linker error messages and see what else you need.
If you are getting error messages of the kind of undefined reference, that means that you are correctly importing header files, but the dynamic library can’t be found. You then need to figure out where the undefined reference is coming from (e.g. which library) and then you have to located that library and copy it into your build/ folder, or wherever you keep your executable.
The joy of developing on Windows …
Compiling on UNIX
Developing on UNIX is much better for our sanity and once you understand the basics, it is really rather straightforward and also somewhat fool-proof. Starting out with code development, especially linking against external libraries, is always easier on UNIX and thus probably the best starting point. If you want to follow along on a UNIX distribution or macOS, then you need the below-shown buildAndRun.sh script:
#!/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 ./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 -o ./build/cgnsTest -L~/libs/lib -lcgns
# run executable
./build/cgnsTestLines 3-4 simply remove and recreate the build/ folder, but there is no need to copy and dynamic libraries around, which on UNIX have a file ending of *.so (Linux) or *.dylib (macOS). The reason is simple, UNIX provides default locations where libraries and header-include files are stored (either within the /usr/ or /usr/local/ folder). So, if I simply installed CGNS, HDF5 and zlib without specifying an install directory, they would end up in /usr/include/ and /usr/lib/ and compilers are instructed to always look into these folders when looking for libraries that do not live next to the executable (we can still use this Windows hack if we want to).
In our case, we install the libraries in a non-default location, i.e. ~/libs, thus the libraries are located in ~/libs/lib/ and the header-include files are in ~/libs/include/. We can see that we use the include path on lines 7-8 when we compile the meshReaderBase.cpp and main.cpp source files, and then we link against the CGNS library on line 11 and specify the location of the library with the -L flag.
If you looked inside the ~/libs/lib folder, you would find a file called libcgns.a, where a stands for archive and is Linux terminology for a static library. The UNIX convention is that in this case, the prefix lib and the file ending .a is removed from the library name, so we only link against cgns with the -l flag (-lcgns) instead of writing -llibcgns.a on line 11.
Notice the absence of any UNIX-specific libraries or even the need to link against HDF5. If we had to link against any dynamic library that wasn’t in the default /usr/lib/ folder, we could provide the linker with the runtime path to check for any dynamic libraries. We would achieve that with the -Wl,-rpath ~/libs/lib flag, where -Wl is a flag saying that anything that follows after the comma should be given to the linker. The linker then receives the path variable (runtime path) and we can specify where the dynamic libraries are located. Windows does not provide this feature, which is annoying.
One thing to note here: if we manually compile and link our code, make sure that during linking, you provide object files first (meshReaderBase.o, main.o) and afterwards, you provide library information with the -L and -l flags. The order is important. Otherwise, the linker may not find the required libraries by the time it needs to know about them.
The last line (line 14) simply executes our executable. Now, this is much cleaner than the Windows build script, but it is good to know the differences between operating systems. I’ll leave you to decide which one you prefer for development.
The mesh reading base class
So then, with all of the prerequisites out of the way, let us dive into the first class, the mesh reading base class. We’ll look at some helper files first and then look at the base class in more detail.
Custom types
The first thing I want to look at is the meshReaderLib/include/types.hpp file. As mentioned previously, I like to have a dedicated file in which I can define custom types and enums. For this library, we will only create a few enums but at other times it may be necessary to define some custom types. I provided some examples in my article on C++ templates, which you may want to check out.
The content of the meshReaderLib/include/types.hpp file is given below:
#pragma once
enum COORDINATE {X = 0, Y};
enum INTERFACE {OWNER = 0, NEIGHBOUR};
enum BC {WALL = 0, INLET, OUTLET, SYMMETRY};We want to ensure we only include the file once; hence we use the #pragma once statement at the top (and we’ll do that for all remaining header files). The reason is to avoid cyclic dependencies which, granted, will not happen here even if we remove the #pragma once statement, but I am just consistently applying it to my library code. There is a detailed description of header include guards over at learncpp.com, which covers why we need #pragma once statements.
The enums that are defined will help us to index arrays in a much more readable way. They are just syntactic sugar, and allow us to transform a less readable loops like
for (int i = 0; i < 1000; ++i)
std::cout << coordinates[i][0] << ", " << coordinates[i][1] << std::endl;into something more meaningful and easier to understand like
for (int vertex = 0; vertex < numberOfVertices; ++vertex)
std::cout << coordinates[vertex][COORDINATE::X] << ", " << coordinates[vertex][COORDINATE::Y] << std::endl;Comparing the two loops, I think we can all come to the conclusion that the second variant is easier to understand. If you need some additional help with enums, learncpp got you covered.
The INTERFACE enum will help us to identify zones on the left/right of the interface, and the BC enum is there to allow us easy assignment of boundary conditions that we want to support.
The library interface
The next file we have a brief look at is the meshReaderLib/meshReader.hpp file. This file will contain links to all other header files so that when we include this file somewhere in a different code which wants to consume our mesh reading library, then it will have access to all function signatures and the compiler can go ahead and compile the code. This is the only header file where we don’t need the header include guards (#pragma once, and yes I know, I said that we are going to apply them across all header files, but since this header file is not implementing anything but rather just collecting header files, it is exempt).
#include "meshReaderLib/include/readMeshBase.hpp"
#include "meshReaderLib/include/types.hpp"We will add to this file as we add more header files to our library, namely the structured and unstructured mesh reader classes.
The base class interface
So then, this is the first bit of the mesh reading library. The base class actually already implements quite a bit of logic which is shared across both structured and unstructured mesh reading. Thus, we can collect this in the base class. The file meshReaderLib/include/readMeshBase.hpp is shown 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 <iostream>
#include <filesystem>
#include <cassert>
#include <vector>
#include <array>
#include <unordered_map>
#include "cgnslib.h"
#include "meshReaderLib/include/types.hpp"
template<typename IndexType>
struct MESHREADERLIB_EXPORT InterfaceConnectivity {
std::array<int, 2> zones;
std::vector<IndexType> ownerIndex;
std::vector<IndexType> neighbourIndex;
};
template<typename IndexType>
struct MESHREADERLIB_EXPORT BoundaryConditionInformation {
std::string boundaryName;
unsigned boundaryType;
std::vector<IndexType> indices;
};
class MESHREADERLIB_EXPORT ReadMeshBase {
public:
using ZoneType = typename std::vector<std::vector<cgsize_t>>;
public:
ReadMeshBase(std::filesystem::path cgnsFilePath);
~ReadMeshBase();
// interface / API
public:
virtual void readMesh() = 0;
// internal implementation hidden from the user
protected:
void readBases();
void readZones(ZoneType_t zoneType);
void createBoundaryNameLookupTable();
// internal implementations that must be implemented by derived classes
protected:
virtual void readCoorinates() = 0;
virtual void readInterfaceConnectivity() = 0;
virtual void readBoundaries() = 0;
unsigned getBoundaryConditionID(BCType_t bcType);
protected:
std::filesystem::path _cgnsFilePath;
int _numberOfZones = 0;
int _numberOfBoundaries = 0;
int _fileIndex = 0;
ZoneType _zoneSizesMax, _zoneSizesMin;
std::unordered_map<std::string, cgsize_t> _zoneNamesMap;
std::unordered_map<std::string, BCType_t> _boundaryConditionMap;
};There is quite a bit to look through, so let’s split this file up into its separate parts and look at that in more detail.
Required pre-processor directives for dynamic library linkage
At the top, we define some preprocessor directives (lines 3-13). We saw in our previous series on writing a linear algebra library, that we need to define these directives to allow for dynamic libraries to be compiled on Windows. That’s all, if you are only targeting UNIX, then you can get rid of them. I won’t repeat the reasoning or a more detailed description here as well have already discussed this in much greater depth in the above-linked article.
Custom types for interface and boundary condition information
After defining some header includes that we’ll need to implement the library, there are two additional structures we define, repeated below for convenience.
template<typename IndexType>
struct MESHREADERLIB_EXPORT InterfaceConnectivity {
std::array<int, 2> zones;
std::vector<IndexType> ownerIndex;
std::vector<IndexType> neighbourIndex;
};The InterfaceConnectivity type (defined on lines 26-31 in the header file above) defines interface connectivity information. When we store an interface, we want to know the two zones that are connected by the interface (which are stored on the std::array<int, 2> zones array). There will be two zones, one on each side of the interface, and thus we need to define the array of type int with 2 entries.
Next, we store the ownerIndex and neighbourIndex arrays within a std::vector. For each vertex in the interface, we want to know the index in the owning and the neighbouring zone. Since we don’t know how many vertices there will be per zone, we have to declare them as std::vector and then at runtime determine how many indices we want to store on the owning and neighbouring zones. Take again how structured and unstructured grid as an example:
We said that there is a single interface between both black and red zones. If we concentrate on the structured grid for the moment, and say that the origin is at the bottom left corner, then we can make the following statements:
The size of the ownerIndex and neighbourIndex std::vector will be in both cases 5 (5 vertices across the interface). If we start at 0, 0 at the bottom left, and the x-axis is going to the right, the y-axis to the top, then the ownerIndex std::vector will look like the following:
std::vector<IndexType> ownerIndex{{4, 0}, {4, 1}, {4, 2}, {4, 3}, {4, 4}};I am assuming that the interface starts at the bottom and goes to the top, but it may be stored the other way around. The x index is always 4, as we are at the right-most part of the zone and only the index in y is changing. The point is, we know all vertices through the ownerIndex. Similarly, the neighbourIndex would be stored in the following way:
std::vector<IndexType> neighbourIndex{{0, 0}, {0, 1}, {0, 2}, {0, 3}, {0, 4}};In this case, we are at the left-most part, x is again constant and 0, while the index only changes in the y-direction.
For both indices, we define the type of the std::vector through a template. We do that because the index type will change from structured to unstructured. For structured meshes, as we saw above, we will have the same number of indices as dimensions, i.e. in our 2D grid we have 2 indices for the x and y directions. For unstructured grids, though, we will always store information as 1D arrays. We store a 1D array of indices, through which we can find out what the coordinates are, for example, by indexing either the x or y coordinate array with the indices stored on the interface. Thus, using a template here allows us to use the struct for both structured and unstructured grids.
There is a second struct for boundary condition information, repeated below for convenience and given on lines 33-38 in the header file given above:
template<typename IndexType>
struct MESHREADERLIB_EXPORT BoundaryConditionInformation {
std::string boundaryName;
unsigned boundaryType;
std::vector<IndexType> indices;
};This struct is similar to the interface connectivity struct in that it stores information for each boundary condition such as its name, its type and again the indices that allow us to find the location of the boundary. For the structured mesh shown above, we may have
std::vector<IndexType> indices{{0, 0}, {1, 0}, {2, 0}, {3, 0}, {4, 0}};For the boundary at the bottom (south) face, which has a type of wall and a name of wall. We store this information for all boundaries and thus can then later determine whether to use a Dirichlet or Neumann condition for the flow quantities solved for at that boundary.
Base class for the mesh reading library
The base class definition is given below again for convenience:
class MESHREADERLIB_EXPORT ReadMeshBase {
public:
using ZoneType = typename std::vector<std::vector<cgsize_t>>;
public:
ReadMeshBase(std::filesystem::path cgnsFilePath);
~ReadMeshBase();
// interface / API
public:
virtual void readMesh() = 0;
// internal implementation hidden from the user
protected:
void readBases();
void readZones(ZoneType_t zoneType);
void createBoundaryNameLookupTable();
// internal implementations that must be implemented by derived classes
protected:
virtual void readCoorinates() = 0;
virtual void readInterfaceConnectivity() = 0;
virtual void readBoundaries() = 0;
unsigned getBoundaryConditionID(BCType_t bcType);
protected:
std::filesystem::path _cgnsFilePath;
int _numberOfZones = 0;
int _numberOfBoundaries = 0;
int _fileIndex = 0;
ZoneType _zoneSizesMax, _zoneSizesMin;
std::unordered_map<std::string, cgsize_t> _zoneNamesMap;
std::unordered_map<std::string, BCType_t> _boundaryConditionMap;
};We define a ZoneType on line 3, which will store the sizes of each zone (i.e. number of vertices per zone). We use this type on line 32. Lines 6-7 define the constructor and destructor, nothing too interesting happening here, we just take the path to a CGNS file as a parameter to the constructor.
Then, there is only a single pure virtual function for the public interface, defined on line 11. This is the only function we can call from outside the base class and will go through the motion of reading the entire mesh from the CGNS file. We looked at what virtual functions are and how virtual and pure virtual functions differ. If you need a refresher, then have a look at my article on how to handle inheritance in C++.
On lines 21-23, there are a few more pure virtual functions that need to be implemented by both the structured and unstructured mesh reading class (that we will implement in the next articles). These will read the coordinates, interfaces and boundary conditions. Since these give us slightly different information (and, in the case of the interface, they may even be stored under different CGNS nodes and thus require completely different functions to read the information required), we delegate the responsibility to implement these functions to the child classes (i.e. the structured and unstructured mesh reading classes).
Line 24 provides a function definition to map CGNS boundary conditions IDs to the boundary IDs we defined above in the meshReaderLib/include/types.hpp file, specifically, the BC enum. We redefine the IDs here, as we do not want to store all the boundary IDs and only use a subset of all available boundary conditions. You’ll also notice that this function, along with the ones defined on lines 21-23, are all marked as protected (line 20). This is to limit their scope to child classes which will derive from this base class.
We looked at access modifiers (public, protected, private) in my article on object-orientated programming. They are commonly used together with class inheritance so if you do not feel comfortable yet with either of these concepts, the two articles linked above from my series on what every CFD developer needs to know about C++ will get you started.
Lines 27-30 store a few self-explanatory variables, except for the _fileIndex variable, perhaps. This is just a file ID used to identify a CGNS file, in case we have more than one open. In this case, we have a mechanism to distinguish between several CGNS files that may be open. Lines 33-34 define two helper variables that we will need to identify the zone and boundary condition we are currently on. We will populate them in the source file (below) and then make use of them later in the child classes.
I really should have provided a custom type for the variables on lines 33-34, this must have slipped through my review before starting to write up this series. It’s annoying me endlessly but it’ll do its job.
Before we close this section and move on to the actual implementation, it may be worthwhile to mention a few of the CGNS-specific types. You’ll notice types such as ZoneType_t (line 16), BCType_t (line 24) and cgsize_t (line 33). These are typically some form of int, unsigned, long, short or combinations thereof, see for example a collection of possible C data types here: C Data Types. To figure out which type needs to be used, the CGNS mid-level library provides all types required to store the right information in the correct type, and this has influenced which variable types I have selected here.
The implementation (source file)
Below is the entire content of the meshReaderLib/src/readMeshBase.cpp file, i.e. the implementation of some functions from the base class that can already be implemented. We’ll chop them off into smaller chunks afterwards again to make them easier to discuss in their own subsections.
#include "meshReaderLib/include/readMeshBase.hpp"
ReadMeshBase::ReadMeshBase(std::filesystem::path cgnsFilePath) : _cgnsFilePath(cgnsFilePath) {
if (cg_open(_cgnsFilePath.string().c_str(), CG_MODE_READ, &_fileIndex)) cg_error_exit();
readBases();
createBoundaryNameLookupTable();
}
ReadMeshBase::~ReadMeshBase() {
cg_close(_fileIndex);
}
void ReadMeshBase::readBases() {
// check how many bases there are. Only a single base per mesh is supported, which should always be true.
int numberOfBases;
if (cg_nbases(_fileIndex, &numberOfBases)) cg_error_exit();
assert(numberOfBases == 1 && "Mesh reading only supported for a mesh with a single base");
}
void ReadMeshBase::readZones(ZoneType_t zoneTypeRequired) {
char zonename[64];
cgsize_t sizes[3][3];
// check how many zones there are in the CGNS file
if (cg_nzones(_fileIndex, 1, &_numberOfZones)) cg_error_exit();
_zoneSizesMax.resize(_numberOfZones);
_zoneSizesMin.resize(_numberOfZones);
for (int zone = 0; zone < _numberOfZones; ++zone) {
// check that the zone within the mesh file contains an structured/unstructured mesh, depending on function argument
ZoneType_t zoneType;
if (cg_zone_type(_fileIndex, 1, zone + 1, &zoneType)) cg_error_exit();
assert(zoneType == zoneTypeRequired && "Wrong mesh type read, expected either structured or unstrucutred");
// ensure that we only have 1 grid per zone. This should always be the case but we check to avoid nasty surprises
int numberOfGrids;
if (cg_ngrids(_fileIndex, 1, zone + 1, &numberOfGrids)) cg_error_exit();
assert(numberOfGrids == 1 && "Mesh reading only supported for a mesh with a single grid per zone");
// get the name of the zone and store it in a temporary map. We'll need that for reading boundary conditions later
if (cg_zone_read(_fileIndex, 1, zone + 1, zonename, sizes[0])) cg_error_exit();
_zoneNamesMap.emplace(std::string(zonename), zone);
// we only support 2D meshes for now, make sure this is the case
int numberOfCoordinatesToRead;
if (cg_ncoords(_fileIndex, 1, zone + 1, &numberOfCoordinatesToRead)) cg_error_exit();
assert(numberOfCoordinatesToRead == 2 && "Mesh reading only supported for 2D meshes");
// store the number of vertices (start and end) for each coordiante (e.g. x and y). Start is always 1
_zoneSizesMin[zone].resize(numberOfCoordinatesToRead);
_zoneSizesMax[zone].resize(numberOfCoordinatesToRead);
for (int coordinate = 0; coordinate < numberOfCoordinatesToRead; ++coordinate) {
_zoneSizesMin[zone][coordinate] = 1;
_zoneSizesMax[zone][coordinate] = sizes[0][coordinate];
}
}
}
void ReadMeshBase::createBoundaryNameLookupTable() {
// In CGNS, we can read boundary conditions in two ways; either by placing them directly onto the grid, or by
// inserting an additional layer (called a family) which connects the grid with a boundary condition. The
// advantage of using a family node is that the boundary condition remains intact, even if the mesh is changed.
// The disadvantage is that we have to do a bit more work to read the boundary condition information. First, we
// need to read the family information and then get the boundary condition name and type from the family. Later,
// we need to to map that information to the correct boundary condition. In the following, we check if the
// boundary conditions are stored under a family (i.e. the bc type is FamilySpecified) and then we read the family
// node with all of its boundary information, which we'll store in the boundaryConditionMap for later lookup.
int numberOfFamilies = 0;
if (cg_nfamilies(_fileIndex, 1, &numberOfFamilies)) cg_error_exit();
for (int family = 0; family < numberOfFamilies; ++family) {
char familyName[128];
int isBoundaryCondition, nGeo;
if (cg_family_read(_fileIndex, 1, family + 1, familyName, &isBoundaryCondition, &nGeo)) cg_error_exit();
char familyBoundaryName[128];
BCType_t familyBoundaryType;
// only if isBoundaryCondition is equal to 1 are boundary conditions stored. For a value other than 1, other
// family information may be stored e.g. different volume conditions (useful for separating zones/volumes to
// apply MRF or sliding grids, for example). We are only interested in reading boundary conditions, though.
if (isBoundaryCondition == 1) {
if (cg_fambc_read(_fileIndex, 1, family + 1, 1, familyBoundaryName, &familyBoundaryType)) cg_error_exit();
_boundaryConditionMap.emplace(familyName, familyBoundaryType);
}
}
}
unsigned ReadMeshBase::getBoundaryConditionID(BCType_t bcType) {
if (bcType == 20 || bcType == 21 || bcType == 22 || bcType == 23 || bcType == 24)
return BC::WALL;
else if (bcType == 9 || bcType == 10 || bcType == 11 || bcType == 18)
return BC::INLET;
else if (bcType == 13 || bcType == 14 || bcType == 15 || bcType == 19)
return BC::OUTLET;
else if (bcType == 16 || bcType == 17)
return BC::SYMMETRY;
else {
throw std::runtime_error("Unsupported boundary condition type specified");
exit(-1);
}
}The constructor and destructor
The constructor and destructor of the ReadMeshBase class are given below.
ReadMeshBase::ReadMeshBase(std::filesystem::path cgnsFilePath) : _cgnsFilePath(cgnsFilePath) {
if (cg_open(_cgnsFilePath.string().c_str(), CG_MODE_READ, &_fileIndex)) cg_error_exit();
readBases();
createBoundaryNameLookupTable();
}
ReadMeshBase::~ReadMeshBase() {
cg_close(_fileIndex);
}Let’s spend some time on the constructor, lines 1-5, as it we will see this function calling pattern repeating. Each CGNS function call returns an error code. We could simply store this error code in a variable and then check if it is any other value than zero. However, there is an accepted shorthand notation, and that is the one you see on line 2.
We place each CGNS function call inside an if statement. Should the if statement evaluate to true, then we call cg_error_exit(), which you can see at the end of line 2. We make use of a little trick here; both true and false evaluate to numerical values, where we have true = 1 and false = 0. So, if the CGNS function call returns an error code of 0, meaning the function call was successful, then the if statement evaluates to false and the cg_error_exit() function is not called. On the other hand, if there is an error, the return code will be anything except 0 and thus evaluate to true, so we do call cg_error_exit().
The cg_error_exit() function remembers the return code of the last called CGNS function, and there are a few error codes defined in the library. So, if there was an error, the function will give you a brief description of what went wrong, which should help in debugging your code.
Let’s then look at the CGNS function itself. On line 2, we make a call to cg_open() to open a CGNS file and, on line 8, inside the destructor, we close the CGNS file again with a call to cg_close(). The cg_open() function expects a path and filename, which we provide as the first argument.
We use here the C++ filesystem header, which makes it easier to store path and filename information across different operating systems. But the CGNS library is written in plain C and so expects a char for its filename. A file path from the filesystem library can be converted to a C++ std::string using the string() method, and, a C++ std::string can be converted to a C-style string using the c_str() method. This is why you see this chain of functions coupled together, we are downcasting from a filesystem path variable to a C-style char.
The second argument defines whether to open the file in read-only or edit mode, essentially. We want to read from the file only, without making any modifications, so we use there the CGNS-defined variable CG_MODE_READ. Finally, we pass the variable _fileIndex as the last argument. The importance here is that we pass the variable by reference (using the & symbol in front of the variable) and so the function will modify this value for us. If you are not comfortable with references, I took a deep dive into them in my article on memory management.
This is a common and repeating pattern in the CGNS library. We pass all inputs first and follow them by outputs. In this example, _cgnsFilePath and CG_MODE_READ (i.e. the first two arguments) are inputs, and the last variable (_fileIndex) is an output of this function. By output I mean that the value of the variable is changed. If you look up the function definitions for cg_open(), for example, in the CGNS mid-level library, you’ll find the following representation:
You can see that the first two arguments are shown in blue, and this means they are inputs, whereas the last argument is given in red(-ish) colour. Whenever you look up the function names, this is how you differentiate between inputs and outputs.
There are two more functions the constructor calls; the readBases() and createBoundaryNameLookupTable(). Both of these functions are implemented below and we make sure that both functions are read before we proceed with calling any other function within the class.
Reading the number of bases
The code for the readBases() function is given below. But before we look at the function, remember that a base node sits just underneath the root node of a CGNS file. We looked at base nodes and the general structure in my first article of this series. But just as a quick recap, we need several bases to store either a deforming mesh, or, keeping the mesh the same, different flow solutions (if we decide to store the flow solution inside a CGNS file, in which case, unsteady data could be written as a series of different bases). For mesh reading, we don’t have different flow solutions to worry about and general of static meshes, so we would expect only a single bases.
void ReadMeshBase::readBases() {
// check how many bases there are. Only a single base per mesh is supported, which should always be true.
int numberOfBases;
if (cg_nbases(_fileIndex, &numberOfBases)) cg_error_exit();
assert(numberOfBases == 1 && "Mesh reading only supported for a mesh with a single base");
}First, we define a local variable numberOfBases here, which we pass on line 4 to the cg_nbases() function. This function will check, for a given file index (remember, the file index is used to identify a CGNS file that was opened before and we received the _fileIndex back as an output from our call to cg_open()), how many bases there are. Our file CGNS reading is based on the assumption that there is only 1 base, which is not a very strong assumption to make, so this should be true in general. We make sure this is true with the assert statement on line 5.
The reason we use assert statements here is the same as with our linear algebra solver library we looked at in the previous series. We should use as many assert statements as we can in our code as they come for free. What do I mean by that? During development, we want to make sure that the code is behaving as expected, and assert statements help us to verify this automatically at various checkpoints. Here, for example, we expect only to get a single base back, otherwise, the rest of our code will behave unexpectedly and we are unlikely to succeed in reading the mesh successfully.
We can always turn assert statements off (using the /DNDEBUG (Windows) or -DNDEBUG (UNIX) compiler arguments, NDEBUG = no debug), and during the pre-processor state, all assert statements will be removed from our file before they are compiled. In this way, once we are happy that the code is working, we remove all debugging statements and don’t lose performance (hence assert statements come for free, in terms of computational overhead). This philosophy is known as the fail-fast method.
Reading all zones and zone information
Next, we want to read, for each base (of which we only have one), the number of zones and then store some content for each zone. The full code for reading zones is given below:
void ReadMeshBase::readZones(ZoneType_t zoneTypeRequired) {
char zonename[64];
cgsize_t sizes[3][3];
// check how many zones there are in the CGNS file
if (cg_nzones(_fileIndex, 1, &_numberOfZones)) cg_error_exit();
_zoneSizesMax.resize(_numberOfZones);
_zoneSizesMin.resize(_numberOfZones);
for (int zone = 0; zone < _numberOfZones; ++zone) {
// check that the zone within the mesh file contains an structured/unstructured mesh, depending on function argument
ZoneType_t zoneType;
if (cg_zone_type(_fileIndex, 1, zone + 1, &zoneType)) cg_error_exit();
assert(zoneType == zoneTypeRequired && "Wrong mesh type read, expected either structured or unstrucutred");
// ensure that we only have 1 grid per zone. This should always be the case but we check to avoid nasty surprises
int numberOfGrids;
if (cg_ngrids(_fileIndex, 1, zone + 1, &numberOfGrids)) cg_error_exit();
assert(numberOfGrids == 1 && "Mesh reading only supported for a mesh with a single grid per zone");
// get the name of the zone and store it in a temporary map. We'll need that for reading boundary conditions later
if (cg_zone_read(_fileIndex, 1, zone + 1, zonename, sizes[0])) cg_error_exit();
_zoneNamesMap.emplace(std::string(zonename), zone);
// we only support 2D meshes for now, make sure this is the case
int numberOfCoordinatesToRead;
if (cg_ncoords(_fileIndex, 1, zone + 1, &numberOfCoordinatesToRead)) cg_error_exit();
assert(numberOfCoordinatesToRead == 2 && "Mesh reading only supported for 2D meshes");
// store the number of vertices (start and end) for each coordiante (e.g. x and y). Start is always 1
_zoneSizesMin[zone].resize(numberOfCoordinatesToRead);
_zoneSizesMax[zone].resize(numberOfCoordinatesToRead);
for (int coordinate = 0; coordinate < numberOfCoordinatesToRead; ++coordinate) {
_zoneSizesMin[zone][coordinate] = 1;
_zoneSizesMax[zone][coordinate] = sizes[0][coordinate];
}
}
}There are two things to notice here; first, while we do call the function readBases() in the constructor, we do not call the function readZones() (above). Secondly, the function expects one function argument (zoneTypeRequired) which is either going to be Structured or Unstructured. These two values are defined in the CGNS header file within an enum. Thus, later when we derive the structured and unstructured mesh reading class, we have to make a call to this function to get all zone information.
We are primarily interested in reading the sizes of the zones, i.e. how many vertices there are for each zone. This will allow us to allocate sufficient memory before we read coordinates. First, we read the number of zones on line 6, which follows a similar pattern to the reading of the number of bases. Notice, though, that we had cg_nbases(_fileIndex, &numberOfBases) before and now we have cg_nzones(_fileIndex, 1, &_numberOfZones). When we read the number of zones, we also have to specify for which base (index) we want to read zone information.
If we had several bases, then we would have to loop over all bases and then read the number of zones for each base, e.g.:
int numberOfBases;
if (cg_nbases(_fileIndex, &numberOfBases)) cg_error_exit();
for (int base = 0; base < numberOfBases; ++base) {
int numberOfZones;
if (cg_nzones(_fileIndex, base + 1, &_numberOfZones)) cg_error_exit();
// read rest of the zone information ...
}CGNS indices are Fortran-based, meaning that the lowest index is 1, not 0, as in C/C++, so we need to specify an offset by using base + 1. We will see this pattern repeating for other CGNS functions as well, as we go deeper into the zones.
Back to the code above, after we have read the number of zones, we loop over all of them on line 10. On line 13, we check that we are reading the expected zone type (i.e. either a structured or unstructured grid) which we assert on line 14. Line 18 checks that there is only a single grid for this zone, just to make sure that the mesh was written correctly by the mesh generator (we expect to have only a single collection of grid information per zone).
Line 22, then, reads the actual information that we are interested in, specifically the sizes array, which we specified on line 3. We pass in the first argument, and, since C-arrays are contiguous in memory, even if they are 2-dimensional (as in our case, i.e. we have sizes[3][3]), the function can write information into the other allocated memory slots as well. Also notice that since we loop over all zones now, we have to specify for each function call which zone to read information for, e.g. cg_function(_fileIndex, 1, zone + 1, ...). We hard code the index for the bases to be 1 (as we only have one) but we need to read all zone information separately as we loop over all zones.
We’ll later need a lookup table that tells us the zone index for a given zone name. We store that on line 23 and will see that in a later article in action. Line 27 enforces the fact that we are currently only able to read 2D grids. When I originally planned this series, I wanted to support both 2D and 3D but that got quickly out of hand in terms of the code required, it can be done, but look at how long this article is already, and we haven’t even gotten into any of the mesh reading, just some helper routines. And, there would have been a lot more templates, which may have been difficult to follow then.
Lines 31-36 then basically store the size of each zone. We specify the minimum index, which is always 1 (remember, again, CGNS uses Fortran-based indexing, 1 is the smallest, not 0). Then, we store the maximum index for each zone and coordinate (i.e. the x or y direction). This allows us, again, to correctly size our arrays before we pass them back to CGNS to read our coordinates. The reason we have to specify the smallest index as well is generality. We could decide to only read a portion of a grid if that is what we wanted to do. There may be some case-specific reasons why you would want to do that, but during mesh reading, we typically want to read everything in one go.
Creating a boundary condition lookup table
Next, we will create a boundary condition lookup table. We need this in case the boundary conditions are written under a family node as discussed above. If boundary conditions are stored in the zone, then there will be no family nodes and so the loop on line 13 will never execute. There is a longer comment at the top which summarises this function as well.
void ReadMeshBase::createBoundaryNameLookupTable() {
// In CGNS, we can read boundary conditions in two ways; either by placing them directly onto the grid, or by
// inserting an additional layer (called a family) which connects the grid with a boundary condition. The
// advantage of using a family node is that the boundary condition remains intact, even if the mesh is changed.
// The disadvantage is that we have to do a bit more work to read the boundary condition information. First, we
// need to read the family information and then get the boundary condition name and type from the family. Later,
// we need to to map that information to the correct boundary condition. In the following, we check if the
// boundary conditions are stored under a family (i.e. the bc type is FamilySpecified) and then we read the family
// node with all of its boundary information, which we'll store in the boundaryConditionMap for later lookup.
int numberOfFamilies = 0;
if (cg_nfamilies(_fileIndex, 1, &numberOfFamilies)) cg_error_exit();
for (int family = 0; family < numberOfFamilies; ++family) {
char familyName[128];
int isBoundaryCondition, nGeo;
if (cg_family_read(_fileIndex, 1, family + 1, familyName, &isBoundaryCondition, &nGeo)) cg_error_exit();
char familyBoundaryName[128];
BCType_t familyBoundaryType;
// only if isBoundaryCondition is equal to 1 are boundary conditions stored. For a value other than 1, other
// family information may be stored e.g. different volume conditions (useful for separating zones/volumes to
// apply MRF or sliding grids, for example). We are only interested in reading boundary conditions, though.
if (isBoundaryCondition == 1) {
if (cg_fambc_read(_fileIndex, 1, family + 1, 1, familyBoundaryName, &familyBoundaryType)) cg_error_exit();
_boundaryConditionMap.emplace(familyName, familyBoundaryType);
}
}
}First, we check the number of families (lines 10-11), which, like the zones, also sit under the base node. We loop over all families, if there are any, on line 13, and then read the family information on line 16. Of particular interest is the flag isBoundaryCondition. If this is set to true (even though the CGNS mid-level library requires you to specify these values as an int, but we also saw above that true = 1 and false = 0, i.e. they are just integers, essentially), then we can expect boundary conditions to be present.
In my previous article we saw that we may have volume conditions as well stored under the family nodes so we don’t want to read them as part of the boundary condition reading task. So, only if isBoundaryCondition == true (or equal to 1, as we check on line 24), we want to read boundary information, which we do on line 25. In particular, we are interested in the boundary name and the type. The name is what we gave during mesh generation, and the type is one of the 21 different CGNS types we can choose from.
Mapping CGNS boundary conditions ID
Remember a bit further up we defined within the meshReaderLib/include/types.hpp file a few enums? One of them was given as
enum BC {WALL = 0, INLET, OUTLET, SYMMETRY};These are the boundary conditions we want to support, for now (and who knows, maybe we get periodic boundary conditions one day in a sensible way, and we may add one more to the list). I come very much from an incompressible background, with these 4 boundary conditions you get quite far, if you are dealing with compressible flows, you may need more than that (or maybe you just use farfield condition, in which case you actually need less). The point is, I want to specify the boundary conditions that my mesh should have, and then I need to map all available CGNS boundary conditions to my types. This is done in the code below.
unsigned ReadMeshBase::getBoundaryConditionID(BCType_t bcType) {
if (bcType == 20 || bcType == 21 || bcType == 22 || bcType == 23 || bcType == 24)
return BC::WALL;
else if (bcType == 9 || bcType == 10 || bcType == 11 || bcType == 18)
return BC::INLET;
else if (bcType == 13 || bcType == 14 || bcType == 15 || bcType == 19)
return BC::OUTLET;
else if (bcType == 16 || bcType == 17)
return BC::SYMMETRY;
else {
throw std::runtime_error("Unsupported boundary condition type specified");
exit(-1);
}
}We get one function argument, bcType, which is a number pointing to one of the 21 boundary conditions (you notice that we check up to a bcType value of 24, there are a few boundary conditions that I discard, such as general and user-defined, hence, the actual number of boundary conditions stored in the CGNS header file are more than 21, but there are 21 usable boundary conditions.
We then simply check what value bcType has and then map it to the best approximation. I am using here the numbers of the boundary conditions as this saves some space (but for educational purposes I should have used the types instead, my bad). If you want to figure out which boundary condition corresponds to which index, you will need to inspect the CGNS header file, i.e. cgnslib.h. Below is a screenshot of the relevant section:
With that information, let’s rewrite the first line where we check for wall boundary conditions. So, instead of writing
if (bcType == 20 || bcType == 21 || bcType == 22 || bcType == 23 || bcType == 24)
return BC::WALL;using the information above, we could have also written this line as
if (bcType == BCWall || bcType == BCWallInviscid || bcType == BCWallViscous || bcType == BCWallViscousHeatFlux || bcType == BCWallViscousIsothermal)
return BC::WALL;A bit longer, but definitely more readable. I must of on a mental vacation while writing this function, but well, you get the idea, hopefully. Using these enum types is probably a better way to map the boundary conditions, especially should additional boundary conditions be added in the future, in which case the numbers may not be the same. So, it is not just more readable but also less error prone.
Initial library testing
Well, we start by populating some code into our main.cpp file, which is where we want to start testing our library. For the moment, we can’t actually do that much. We can’t create an object from our base class, because it contains pure virtual functions, and a class containing pure virtual functions can not be instantiated. So we need to wait until we derive from hat class, implement all pure virtual functions and then test those derived classes instead. For now, we simply specify the path to one of the mesh files and print it.
#include <iostream>
#include <filesystem>
#include <cmath>
#include <string>
#include "meshReaderLib/meshReader.hpp"
int main() {
auto structuredMeshPathFamily = std::filesystem::path("mesh/structured2D.cgns");
std::cout << "Mesh path: " << structuredMeshPathFamily << std::endl;
return 0;
}Summary
Right, so this was quite a long write-up, and we haven’t even started writing much in the way of structured od unstructured mesh reading code. This will happen in the next two articles. But for now, we looked at how to set up the basic project structure, which closely follows that of the previous library we wrote. It makes sense to keep these developments in a library form, as this allows us later to rescue the code we are writing now.
We looked at the basic build scripts and the same caveats as for the last series apply; we will eventually have to replace them with proper build scripts and we will do so in the future using CMake. For now, we looked at the raw compiler commands, which should help us understand how to call our compiler with which flags to get the final executable we want.
In this article, then, we looked at the first class, which serves as our base class for mesh reading. It defines a few functions as pure virtual functions which will need to be implemented by the structured and unstructured mesh reading class (e.g. coordinate, interface, and boundary condition reading). The function calls will be very similar, but the information that is stored is different, so there will be some differences in code and hence it makes sense to split this code into their own child classes.
We looked at some common functionality as well, specifically how to read bases and zones, as well as how to extract boundary conditions from family nodes, should there be any. So we already saw how to interact with a CGNS file through the mid-level library calls and we will continue to use this API in upcoming articles.
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.