In our final article of this series, we will take our project and turn it into a library. This means that we have to create new build scripts that can compile our library code as either a static or a dynamic library. This is what we will focus on in the first half of this article. The second part of the article consists of writing code that tests our implementation. Thus far, we have only ensured that our code does compile and run, but we haven’t done anything in the way of ensuring correctness. This is what we focus on in the second part, where we test the meshes that we read against expected values.
By the end of this article, you will have a complete view of the project structure and all of its code and build scripts, and you will gain some basic understanding of how we can test code systematically. These tests not only ensure that our code works correctly but also showcase how we would use the library code in a real project. From these tests, you will learn and understand how you can use the developed library in your code base.
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-5.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 articles
Where we are and what needs to be done
If you have read all the articles in a row, then you will know that we started developing a base class to read a CGNS-based mesh file, from which we then derived a structured mesh and unstructured mesh reading class. There are so far two limitations that I want to address in this article.
First of all, our code is currently compiled into one executable, but I want to separate the library code from any code that is merely calling the mesh reading classes. For this, similar to what we did with our linear algebra solver library, we are going to write new build scripts for Windows and UNIX that will compile the library code into either a static or dynamic library first. We are then going to use that compiled library to inject it into the main.cpp file, which should be compiled separately from the library code.
Secondly, thus far, within the main.cpp file we are only calling the readMesh() function of either the structured or unstructured grid. This ensures two things; first, it makes sure that we don’t have any compilation and syntax errors in our code, and second, we go through all of our assert statements and we can check if any of them triggers.
So far, the code compiled fine and no asserts evaluated to false, so we can have some confidence that the code is working to a certain degree. However, we haven’t tested the library extensively and we haven’t checked if any of the coordinates, interface, or boundary condition information were read correctly.
Therefore, the second part of this article will look at writing some functions within our main.cpp file that will look at the data we received from the structured and unstructured mesh reading classes and check it for correctness. This requires knowledge of the generated mesh, and this is the reason we used a simplified mesh thus far. It is easier to gather information for a simplified mesh which we can then check against the data received from our mesh reading classes.
With this out of the way, we will now look at the modified build scripts and then look at the test code within the main.cpp file next.
Generating build scripts for static and dynamic libraries
Our build scripts thus far have compiled all library code together with the main.cpp file into a single executable. But as alluded to above, we want to separate the build process so that the library code is compiled into its own output, i.e. a static or dynamic library, which will then be consumed by the main.cpp file. The build scripts are going to create two separate outputs and we can reuse the library code later in other projects.
We typically would also want to include an installation step here, i.e. take the compiled library and all the required header files and copy them into an installation directory, from where you can reuse the code in all other types of projects. We’ll skip this step for now but when we will look at the build process in a separate series we will readdress this issue.
For now, we stick to our manually crafted build scripts. If you are a seasoned programmer, then you will be using build scripts and you may feel that this approach is a bit backward thinking. I would agree that we should not be writing manual build scripts these days anymore, however, I believe there is a lot of value in understanding the build process from start to end, before moving to more complex build scripts.
At some point, we will switch to using build systems, and at that point, all the knowledge and experience we have gathered from these manual build scripts will help us to understand what the build system is doing and we can easily write our own build files or extend existing ones. We will have to debug the output from a build system at which point we need to understand the compilation process. So we stick with the manual approach for now but will go to full build systems soon.
Creating a static library
As discussed in my article on static, dynamic, and header-only libraries, a static library gathers all the compiled object files of the library code and puts that into a convenient archive, which we said can be thought of as a simple zip archive. This static library, or archive, will then be injected into the final executable of any code that wants to use this library. Thus, the executable will have all of the compiled source files available and can resolve function calls to the library.
Windows
Compiling a static library on Windows is rather straight forward, the new build script is given below, but it only changed in two lines compared to our previous build script.
# 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 static library
lib.exe /nologo /OUT:build\meshReader.lib build\readMeshBase.obj build\readStructuredMesh.obj build\readUnstructuredMesh.obj
# link library into main executable
cl.exe /nologo /EHsc /std:c++20 .\build\main.obj /Fe".\build\cgnsTest.exe" /link /MACHINE:x64 /LIBPATH:"C:\libraries\lib" /LIBPATH:.\build meshReader.lib cgns.lib hdf5.lib msvcrt.lib libcmt.lib
# clean up
Remove-Item build\*.obj
# test executable
.\build\cgnsTest.exeWe still start by copying all required dependencies into our build/ directory on lines 2-5. We follow that by compiling all source files (library and non-library code, e.g. the main.cpp file) on lines 8-11 but then go ahead and take only the compiled library code on line 14 and pass that to the linker (lib.exe) on line 14.
The linker works again similarly to a zip archiver, taking the compiled object files and compressing them into what we call here meshReader.lib within the build/ folder, i.e. the argument given after the /OUT: flag.
On line 17, we are now taking the compiled output from the main.cpp file and are creating an executable called cgnsTest.exe in the build/ folder. It has several dependencies, one of which is now our newly compiled library meshReader.lib, along with several other dependencies our meshReader.lib library depends on, i.e. CGNS and its dependencies.
UNIX
The UNIX build process is pretty similar and shown below.
#!/bin/bash
rm -rf build
mkdir -p build
# compile source files into object code
g++ -c -std=c++20 -O3 -I. -I~/libs/include meshReaderLib/src/readMeshBase.cpp -o build/readMeshBase.o
g++ -c -std=c++20 -O3 -I. -I~/libs/include meshReaderLib/src/readStructuredMesh.cpp -o build/readStructuredMesh.o
g++ -c -std=c++20 -O3 -I. -I~/libs/include meshReaderLib/src/readUnstructuredmesh.cpp -o build/readUnstructuredmesh.o
g++ -c -std=c++20 -O3 -I. -I~/libs/include main.cpp -o build/main.o
# link object files into static library
ar rcs build/libMeshReader.a build/readMeshBase.o build/readStructuredMesh.o build/readUnstructuredmesh.o
# compile main function and link library into executable
g++ -O3 -I. build/main.o -o build/cgnsTest -Lbuild -L~/libs/lib -lMeshReader -lcgns
# remove object files
rm build/*.o
# run
./build/cgnsTestWe don’t have to copy the dependencies into our build/ folder, UNIX has more advanced ways to deal with external dependencies, thus, we are only creating a build/ folder for our project on lines 3-4.
We compile all code on lines 7-10 and then pass the compiled library code to the ar executable, which is the archiver on UNIX. On line 16, we take the compiled code (object file) from the main.cpp file and create the cgnsTest executable. We link against the dependencies which here are only our mesh reader library, as well as CGNS.
One thing worth pointing out here, again, is that Windows and UNIX have different naming conventions when it comes to library naming. Windows is pretty relaxed and you can call your libraries whatever you want. On UNIX, though, we use the convention that our library has to start with the prefix lib, followed by the library name. In our case, we have meshReader.lib on Windows but libMeshReader.a on UNIX.
You can see on line 13 that we created our static library as libMeshReader.a but then on line 16, we only specify the library name as meshReader, i.e. the string following the -l flag. UNIX expects the prefix lib but then when it comes to specifying the library name, we should not specify it when we include a library using -l flag, nor should we include the file ending, which again UNIX assumes to be *.a.
For Windows, this is entirely different, where we can name our library anything we want. When it comes to including the library, we have to include the entire library name, including its file ending. The same is true for dynamic libraries discussed below, except that the file ending on Windows now becomes *.dll and *.so on UNIX.
Creating a dynamic library
Returning back to my article on static, dynamic, and header-only libraries, a dynamic library stores compiled code in an archive if you want, similar to how we store compiled code for static libraries. The structure is slightly different, but we don’t have to get held up on semantics here.
This archive is not, however, injected into the executable but instead lives on its own. Whenever we need to use code that is compiled into a dynamic library, we simply look for it at the runtime of our executable and then call functionality within that library. Thus, our executable will always be separated from the library which has its advantages and disadvantages.
The advantage is that we keep our file sizes small. We don’t copy and paste the library into our executable, as we do for static libraries, and so keep file sizes low. The downside is that we have to rely on finding the library at runtime. So even if we compile our executable without errors, which relies on a dynamic library, we may still not be able to execute it if the dynamic libraries it depends on can’t be located.
As an analogy, think of music. There are different ways you can consume it. You can have some MP3 files on your phone. You own your phone and everything on it so whenever you have access to your phone, you can play your mp3 files. This is similar to a static library, i.e. in this scenario, your phone is the executable and the MP3 file is the static library. The executable owns the static library (your phone has access to the MP3) and thus you call the library whenever you need it (you can play the MP3 file whenever you want).
A dynamic library is more like streaming music from the internet. Even if you have access to your phone, you may have run out of data or don’t have WIFI access. In that case, the MP3 file still exists somewhere on the internet but you can’t access it because you can’t obtain/find it. Having a phone doesn’t mean you can stream music. Again, the phone here is our executable and it may compile fine, but it may not be able to find the library at runtime, it needs to know where it is located. If it can’t find it, the executable won’t run.
Windows
The build script for the Windows-based PowerShell build script is given 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 /DCOMPILEDLL .\meshReaderLib\src\readMeshBase.cpp /Fo".\build\readMeshBase.obj"
cl.exe /nologo /EHsc /std:c++20 /I. /I "C:\libraries\include" /c /O2 /DCOMPILEDLL .\meshReaderLib\src\readStructuredMesh.cpp /Fo".\build\readStructuredMesh.obj"
cl.exe /nologo /EHsc /std:c++20 /I. /I "C:\libraries\include" /c /O2 /DCOMPILEDLL .\meshReaderLib\src\readUnstructuredMesh.cpp /Fo".\build\readUnstructuredMesh.obj"
cl.exe /nologo /EHsc /std:c++20 /I. /I "C:\libraries\include" /c /O2 /DCOMPILEDLL .\main.cpp /Fo".\build\main.obj"
# link dynamic library
link.exe /nologo /DLL /OUT:build\meshReader.dll build\readMeshBase.obj build\readStructuredMesh.obj build\readUnstructuredMesh.obj /LIBPATH:"C:\libraries\lib" cgns.lib hdf5.lib msvcrt.lib libcmt.lib
# link library into main executable
cl.exe /nologo /EHsc /std:c++20 .\build\main.obj /Fe".\build\cgnsTest.exe" /link /MACHINE:x64 /LIBPATH:"C:\libraries\lib" /LIBPATH:.\build meshReader.lib cgns.lib hdf5.lib msvcrt.lib libcmt.lib
# clean up
Remove-Item build\*.obj
# test executable
.\build\cgnsTest.exeLines 2-5 are the same as before, i.e. we are simply copying all required dependencies (libraries) into the same folder our executable will live in so it can find all of its dependencies. Lines 8-11 compile the code again, the only difference being that we are passing the pre-processor directive of /DCOMPILEDLL to the source files. This is required in the class header files to mark functions either for import or export. We have talked about this in more length both at the beginning of this series and in our previous series on compiling a linear algebra solver library, which you may want to look through for a refresher.
Next, we pass the compiled library code on line 14 to the linker and specify that we want to create a dynamic library using the /DLL (dynamic linked library) flag. As we discussed in the above-linked previous series as well, even though we create the meshReader.dll dynamic library on line 14, Windows creates both the meshReader.dll and meshReader.lib files, where the meshReader.lib file now acts as a wrapper, delegating any function calls to the meshReader.dll file.
UNIX
UNIX is again slightly different, and the build script using bash is given below.
#!/bin/bash
rm -rf build
mkdir -p build
# compile source files into object code
g++ -c -std=c++20 -fpic -O3 -I. -I~/libs/include meshReaderLib/src/readMeshBase.cpp -o build/readMeshBase.o
g++ -c -std=c++20 -fpic -O3 -I. -I~/libs/include meshReaderLib/src/readStructuredMesh.cpp -o build/readStructuredMesh.o
g++ -c -std=c++20 -fpic -O3 -I. -I~/libs/include meshReaderLib/src/readUnstructuredmesh.cpp -o build/readUnstructuredmesh.o
g++ -c -std=c++20 -O3 -I. -I~/libs/include main.cpp -o build/main.o
# link object files into static library
g++ -shared -O3 -Wall -Wextra -I. -o build/libMeshReader.so build/readMeshBase.o build/readStructuredMesh.o build/readUnstructuredmesh.o
# compile main function and link library into executable
g++ -O3 -Wall -Wextra -I. build/main.o -o build/cgnsTest -Lbuild -L~/libs/lib -Wl,-rpath=build -lMeshReader -lcgns
# remove object files
rm build/*.o
# run
./build/cgnsTestWe do not have to provide the COMPILEDLL flag to the compilation step on lines 7-10, as this is a Microsoft-specific implementation of the C++ standard (making it non-standard?!). We do, though, pass the -fpic compiler flag, which allows us to compile position-independent code. In a nutshell, if we want to load library code at runtime, we face the issue that we need to know, during compilation and linking, where the code is going to be located. With dynamic libraries, we don’t know that and so we use a technique used as position-indepdent code, which overcomes this issue.
On line 13, we call here g++, i.e. the compiler, to create the dynamic library. Internally though, the linker ld is called. This is just a convenient shorthand notation (and by convenient, I really mean convenient, as the linker arguments look really rather cryptic), and by specifying the -shared linker flag, the linker knows that it should create a dynamic library (which it can, since the code was compiled position-indepdent).
On line 16, we compile the main.cpp file, which depdens on the dynamic library we created. Thus, we specify where to find the library (using the -L flag to specify the folder, i.e. build/), as well as the name of the libraries using the -l flag, i.e. meshReader and cgns. The linker will, as explained above, extend both of these names and look for the libmeshReader.so and libcgns.so file instead.
Furthermore, we specify additional linker arguments using the -Wl flag, anything following the comma will be passed to the linker. In this case, we specify the rpath which is a concept only available on UNIX and not Windows. The rpath is the runtime path to check for dynamic libraries, i.e. we can bake that information into the executable. It doesn’t mean that the libraries will actually be in that folder, but it adds this folder to the search path of possible folders in which the dynamic libraries could be located and if it is found in this path, the executable will use the libraries from there.
Creating the test code
What we are seeking to achieve at this point is to modify our main.cpp file to perform some tests on our structured and unstructured mesh reading classes. Since we have implemented two different types of boundary condition reading, i.e. either family or boundary node-based (which we discussed in a previous article in length), we want to test both implementations here. This means we will have two tests for the structured grid, and two tests for the unstructured grid, so four in total.
The test code is given below for the main.cpp file. We discuss a broad overview below the code and then jump into separate sections underneath for testing specific tests.
#include <iostream>
#include <filesystem>
#include <cmath>
#include <string>
#include "meshReaderLib/meshReader.hpp"
void readStructuredMeshCoordinates(ReadStructuredMesh &structuredMesh);
void readStructuredMeshInterfaceConnectivity(ReadStructuredMesh &structuredMesh);
void readStructuredMeshBoundaryConditions(ReadStructuredMesh &structuredMesh);
void readUnstructuredMeshCoordinates(ReadUnstructuredMesh &unstructuredMesh);
void readUnstructuredInternalCells(ReadUnstructuredMesh &unstructuredMesh);
void readUnstructuredMeshInterfaceConnectivity(ReadUnstructuredMesh &unstructuredMesh);
void readUnstructuredMeshBoundaryConditions(ReadUnstructuredMesh &unstructuredMesh);
int main() {
auto structuredMeshPathFamily = std::filesystem::path("mesh/structured2D.cgns");
ReadStructuredMesh structuredMeshFamily(structuredMeshPathFamily);
structuredMeshFamily.readMesh();
readStructuredMeshCoordinates(structuredMeshFamily);
readStructuredMeshInterfaceConnectivity(structuredMeshFamily);
readStructuredMeshBoundaryConditions(structuredMeshFamily);
std::cout << "Structured mesh read successfully from CGNS file (family-based boundary conditions)" << std::endl;
// structured mesh reading test, without family-based boundary condition
auto structuredMeshPathNoFamily = std::filesystem::path("mesh/structured2DNoFamily.cgns");
ReadStructuredMesh structuredMeshNoFamily(structuredMeshPathNoFamily);
structuredMeshNoFamily.readMesh();
readStructuredMeshCoordinates(structuredMeshNoFamily);
readStructuredMeshInterfaceConnectivity(structuredMeshNoFamily);
readStructuredMeshBoundaryConditions(structuredMeshNoFamily);
std::cout << "Structured mesh read successfully from CGNS file (non family-based boundary conditions)" << std::endl;
// unstructured mesh reading test, with family-based boundary condition
auto unstructuredMeshPathFamily = std::filesystem::path("mesh/unstructured2D.cgns");
ReadUnstructuredMesh unstructuredMeshFamily(unstructuredMeshPathFamily);
unstructuredMeshFamily.readMesh();
readUnstructuredMeshCoordinates(unstructuredMeshFamily);
readUnstructuredInternalCells(unstructuredMeshFamily);
readUnstructuredMeshInterfaceConnectivity(unstructuredMeshFamily);
readUnstructuredMeshBoundaryConditions(unstructuredMeshFamily);
std::cout << "Unstructured mesh read successfully from CGNS file (family-based boundary conditions)" << std::endl;
// unstructured mesh reading test, without family-based boundary condition
auto unstructuredMeshPathNoFamily = std::filesystem::path("mesh/unstructured2DNoFamily.cgns");
ReadUnstructuredMesh unstructuredMeshNoFamily(unstructuredMeshPathNoFamily);
unstructuredMeshNoFamily.readMesh();
readUnstructuredMeshCoordinates(unstructuredMeshNoFamily);
readUnstructuredInternalCells(unstructuredMeshNoFamily);
readUnstructuredMeshInterfaceConnectivity(unstructuredMeshNoFamily);
readUnstructuredMeshBoundaryConditions(unstructuredMeshNoFamily);
std::cout << "Unstructured mesh read successfully from CGNS file (non family-based boundary conditions)" << std::endl;
return 0;
}Lines 1-6 import the required header, of which line 6 shows the header we need for our library. This will allow us to make use of the mesh reading library we just developed in the previous articles. Lines 8-10 and 12-15 provide the function definitions for the structured and unstructured grid, respectively, that we will use to test the various parts of the mesh reading process.
Both grids have functions defined to read the coordinates, interfaces, and boundary conditions. We will use these functions to test that the data returned by either class are correct and match the data of the grid. The unstructured grid has one more function, which allows us to check the element connectivity reading.
We define the functions here and then we can call them later below the main() function, which is defined on lines 17-64.
The main() function can be subdivided into 4 categories:
- Lines 18-26: Testing the structured mesh reading, with boundary conditions stored under a family node.
- Lines 28-37: Testing the structured mesh reading, with boundary conditions stored per zone under a boundary condition node.
- Lines 39-49: Testing the unstructured mesh reading, with boundary conditions stored under a family node.
- Lines 51-61: Testing the unstructured mesh reading, with boundary conditions stored per zone under a boundary condition node.
Each of these sections first creates an object of either the ReadStructuredMesh or ReadUnstructuredMesh class, after which each object proceeds by calling the readMesh() function defined within its class. Afterwards, we pass the object to one of the functions defined between lines 8-15, i.e. to test that either the coordinate, interface, boundary condition or element connectivity data was read correctly. The sections below will show the code for each of these test functions and discuss how we ensure that the data is read correctly.
Testing structured mesh reading
Before we jump into the various code sections, let’s remind ourselves what information we are reading. Below is a sketch of the structured grid that we want to read. The index location 1,1 on zone 1 is at x=0 and y=0, and the location 5,5 on zone 2 is at x=2 and y=1. Given that we have an equidistant spacing in the x and y direction, it follows that \Delta x=\Delta y=0.25. Thus, we know the coordinates for our adjacent coordinates, e.g. index 2,1 on zone 1 is at x=0.25 and y=0 and index 3,4 on zone 2 is at x=1.5 and y=0.75.
We can use this knowledge now to read our coordinates and test that we have received the correct coordinates, by comparing the value received from the class with an expected value. Similarly, we know the indices at the interfaces and boundaries, and furthermore we can state that we have one interface and three boundaries per zone. There are two zones in total. All of this information will be checked against what we have read from the structured mesh reading class and compared. Let’s start with the coordinate reading.
Testing coordinate reading
In general, when we test code, we want to test if an expression is true or false. So far, we have used assert statements throughout our code to check if certain preconditions are met before we can continue with the code. Have a look at the following examples:
assert(true); // evaluates to true
assert(false); // evaluates to false
assert(1 == 1); // evaluates to true
assert(1 == 2); // evaluates to false
assert(1 < 2); // evaluates to true
assert(1 > 2); // evaluates to falseAll of these will either return true or false. So when we want to test our code, we read something from the class we are trying to test and then compare that with an expected value (that we can construct, for example, the coordinates based on the structured mesh shown above). If the returned value from the class does meet our expected value, then we move on and read the next value. This is the structure that we will adopt in all subsequent test functions, i.e. read information and compare that against expected values.
Returning to the coordinate reading test, after we have created our mesh reading object and read the grid through the readMesh() function inside the main() function, we pass the mesh object into the readStructuredMeshCoordinates() function, which is shown below. We read the coordinates on line 3, which we receive from the mesh object by calling the getter function getCoordinates() from the class. Afterwards, we test that all data received has been read correctly.
On line 6, we assert that the first index in the coordinate array is of size 2, which is the index looping over all zones. Since we have 2 zones for this structured grid (shown above), we expect the size() to return a value of 2.
We can add some additional documentation using the && "documentation" syntax. The and operator && allows us to test more than one expression. So, for example, we could have assert(1 == 1 && 2 == 2);, which would evaluate to true because both expressions are true. However, assert(1 == 1 && 1 == 2); would evaluate to false, as one of these expressions is false. The and operator && requires both expressions to be true in order for the assert to return true overall.
If you try to cast a string to a boolean value, e.g. bool("documentation"), it will always return true and so we can use this to our advantage to pass some text after the expression to state, in words, what we expect this assert statement to test.
After line 6, we go through each zone and each index pair i,j to test if the read coordinate in x and y is what we expect it to be. Since we are now dealing with floating point numbers, we use another little trick here in that we subtract the expected coordinate value from the actual read coordinate value.
This should provide us with a value of zero, but rounding errors may produce values slightly larger or smaller than zero. Thus, we first take the absolute value using the std::fabs() function (floating point absolute value) to only have positive values. If that value is smaller than a given tolerance (here 0.01), then we can state that we have read to coordinates correctly. There are better ways to compare floating point numbers, and we will cover that in a dedicated series on testing in the future.
We are not going to go through each line, but let’s take an example to see how these lines after line 6 work, and then all subsequent lines should make sense. Since we are zero-based with our indices, we need to add 1 if we want to relate our indices to the schematic above. The coordinate array stores first the zones, then the i indices (x-direction), followed by the j indices (y-direction) and finally the coordinates themselves (either x or y).
Thus, if we look at line 16, for example, we have coordinates[0][3][1][COORDINATE::X]. This means we are reading the x coordinate on zone 1 with index i=4 and j=2 (adding one for each index). Looking at the sketch above, we expect the x value to be 0.75 here and so we subtract the 0.75 value from our coordinate array, which should give us a value of zero (or close to zero accounting for any round-off errors).
The corresponding call for the y coordinate can be found on line 46, and we have an expected value of 0.25 for it. In this way, all code shown on lines 8-67 test each coordinate on zone 1. Subsequently, lines 69-128 test all coordinates on zone 2. The difference here is that all x-coordinates are shifted by a value of 1, since zone 2 starts at x=1. However, the y-coordinates are the same as those on zone 1 as there is no offset.
void readStructuredMeshCoordinates(ReadStructuredMesh &structuredMesh) {
// get the coordinates from the structured mesh
auto coordinates = structuredMesh.getCoordinates();
// test that coordinates were read correctly
assert(coordinates.size() == 2 && "Expected 2 zones");
// test zone 0
assert(std::fabs(coordinates[0][0][0][COORDINATE::X] - 0.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][1][0][COORDINATE::X] - 0.25) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][2][0][COORDINATE::X] - 0.50) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][3][0][COORDINATE::X] - 0.75) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][4][0][COORDINATE::X] - 1.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][0][1][COORDINATE::X] - 0.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][1][1][COORDINATE::X] - 0.25) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][2][1][COORDINATE::X] - 0.50) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][3][1][COORDINATE::X] - 0.75) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][4][1][COORDINATE::X] - 1.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][0][2][COORDINATE::X] - 0.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][1][2][COORDINATE::X] - 0.25) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][2][2][COORDINATE::X] - 0.50) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][3][2][COORDINATE::X] - 0.75) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][4][2][COORDINATE::X] - 1.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][0][3][COORDINATE::X] - 0.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][1][3][COORDINATE::X] - 0.25) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][2][3][COORDINATE::X] - 0.50) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][3][3][COORDINATE::X] - 0.75) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][4][3][COORDINATE::X] - 1.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][0][4][COORDINATE::X] - 0.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][1][4][COORDINATE::X] - 0.25) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][2][4][COORDINATE::X] - 0.50) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][3][4][COORDINATE::X] - 0.75) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][4][4][COORDINATE::X] - 1.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][0][0][COORDINATE::Y] - 0.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][1][0][COORDINATE::Y] - 0.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][2][0][COORDINATE::Y] - 0.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][3][0][COORDINATE::Y] - 0.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][4][0][COORDINATE::Y] - 0.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][0][1][COORDINATE::Y] - 0.25) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][1][1][COORDINATE::Y] - 0.25) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][2][1][COORDINATE::Y] - 0.25) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][3][1][COORDINATE::Y] - 0.25) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][4][1][COORDINATE::Y] - 0.25) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][0][2][COORDINATE::Y] - 0.50) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][1][2][COORDINATE::Y] - 0.50) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][2][2][COORDINATE::Y] - 0.50) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][3][2][COORDINATE::Y] - 0.50) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][4][2][COORDINATE::Y] - 0.50) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][0][3][COORDINATE::Y] - 0.75) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][1][3][COORDINATE::Y] - 0.75) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][2][3][COORDINATE::Y] - 0.75) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][3][3][COORDINATE::Y] - 0.75) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][4][3][COORDINATE::Y] - 0.75) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][0][4][COORDINATE::Y] - 1.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][1][4][COORDINATE::Y] - 1.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][2][4][COORDINATE::Y] - 1.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][3][4][COORDINATE::Y] - 1.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][4][4][COORDINATE::Y] - 1.00) < 0.01 && "Unexpected value for Y-coordinate");
// test zone 1
assert(std::fabs(coordinates[1][0][0][COORDINATE::X] - 1.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][1][0][COORDINATE::X] - 1.25) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][2][0][COORDINATE::X] - 1.50) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][3][0][COORDINATE::X] - 1.75) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][4][0][COORDINATE::X] - 2.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][0][1][COORDINATE::X] - 1.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][1][1][COORDINATE::X] - 1.25) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][2][1][COORDINATE::X] - 1.50) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][3][1][COORDINATE::X] - 1.75) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][4][1][COORDINATE::X] - 2.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][0][2][COORDINATE::X] - 1.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][1][2][COORDINATE::X] - 1.25) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][2][2][COORDINATE::X] - 1.50) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][3][2][COORDINATE::X] - 1.75) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][4][2][COORDINATE::X] - 2.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][0][3][COORDINATE::X] - 1.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][1][3][COORDINATE::X] - 1.25) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][2][3][COORDINATE::X] - 1.50) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][3][3][COORDINATE::X] - 1.75) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][4][3][COORDINATE::X] - 2.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][0][4][COORDINATE::X] - 1.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][1][4][COORDINATE::X] - 1.25) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][2][4][COORDINATE::X] - 1.50) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][3][4][COORDINATE::X] - 1.75) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][4][4][COORDINATE::X] - 2.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][0][0][COORDINATE::Y] - 0.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][1][0][COORDINATE::Y] - 0.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][2][0][COORDINATE::Y] - 0.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][3][0][COORDINATE::Y] - 0.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][4][0][COORDINATE::Y] - 0.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][0][1][COORDINATE::Y] - 0.25) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][1][1][COORDINATE::Y] - 0.25) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][2][1][COORDINATE::Y] - 0.25) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][3][1][COORDINATE::Y] - 0.25) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][4][1][COORDINATE::Y] - 0.25) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][0][2][COORDINATE::Y] - 0.50) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][1][2][COORDINATE::Y] - 0.50) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][2][2][COORDINATE::Y] - 0.50) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][3][2][COORDINATE::Y] - 0.50) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][4][2][COORDINATE::Y] - 0.50) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][0][3][COORDINATE::Y] - 0.75) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][1][3][COORDINATE::Y] - 0.75) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][2][3][COORDINATE::Y] - 0.75) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][3][3][COORDINATE::Y] - 0.75) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][4][3][COORDINATE::Y] - 0.75) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][0][4][COORDINATE::Y] - 1.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][1][4][COORDINATE::Y] - 1.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][2][4][COORDINATE::Y] - 1.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][3][4][COORDINATE::Y] - 1.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][4][4][COORDINATE::Y] - 1.00) < 0.01 && "Unexpected value for Y-coordinate");
}Testing interface reading
When it comes to the interface, we receive the indices of the two zones that share the interface, as well as the vertices for each zone that are located on the interface.
In the code below, we receive the interface information on line 3 and then check that we only have one interface available in total (on line 6). Looking at the mesh sketched above, we see that we only have a single interface (between zones 1 and 2). Furthermore, we expect zone 0 to be the owner and zone 1 to be the neighbour, as shown on lines 7-8. We have to remember again that we use zero-based indexing, hence the offset of one.
Lines 9-10 check that we have exactly 5 vertices in each interface, and we can see that from the above sketch as well. Then, lines 12-22 check that we received the correct indices for the owning zone and likewise, lines 24-34 check the same information for the neighbouring zone.
For example, the owning zone starts at the indices i,j=4,0, which we can see in the sketch above (remember to add 1 for the offset). It goes up to the indices i,j=4,4. We see on lines 12-16 that the i index does not change, but lines 18-22 show that the j index changes, which we are expecting. The same statement can be made about the vertices for the neighbouring zone, except that the i index now remains at 0, as the interface is now located at the left-most index (see the red zone above in the sketch).
void readStructuredMeshInterfaceConnectivity(ReadStructuredMesh &structuredMesh) {
// get the interface connectivity from the structured mesh
auto interfaceConnectivity = structuredMesh.getInterfaceConnectivity();
// test that interface connectivity were read correctly
assert(interfaceConnectivity.size() == 1 && "Expected a single interface");
assert(interfaceConnectivity[0].zones[0] == 0 && "Expecting owning zone to be zone 0");
assert(interfaceConnectivity[0].zones[1] == 1 && "Expecting neighbouring zone to be zone 1");
assert(interfaceConnectivity[0].ownerIndex.size() == 5 && "Expected 5 owner indices");
assert(interfaceConnectivity[0].neighbourIndex.size() == 5 && "Expected 5 neighbour indices");
assert(interfaceConnectivity[0].ownerIndex[0][COORDINATE::X] == 4 && "Unexpected owner index in X");
assert(interfaceConnectivity[0].ownerIndex[1][COORDINATE::X] == 4 && "Unexpected owner index in X");
assert(interfaceConnectivity[0].ownerIndex[2][COORDINATE::X] == 4 && "Unexpected owner index in X");
assert(interfaceConnectivity[0].ownerIndex[3][COORDINATE::X] == 4 && "Unexpected owner index in X");
assert(interfaceConnectivity[0].ownerIndex[4][COORDINATE::X] == 4 && "Unexpected owner index in X");
assert(interfaceConnectivity[0].ownerIndex[0][COORDINATE::Y] == 0 && "Unexpected owner index in Y");
assert(interfaceConnectivity[0].ownerIndex[1][COORDINATE::Y] == 1 && "Unexpected owner index in Y");
assert(interfaceConnectivity[0].ownerIndex[2][COORDINATE::Y] == 2 && "Unexpected owner index in Y");
assert(interfaceConnectivity[0].ownerIndex[3][COORDINATE::Y] == 3 && "Unexpected owner index in Y");
assert(interfaceConnectivity[0].ownerIndex[4][COORDINATE::Y] == 4 && "Unexpected owner index in Y");
assert(interfaceConnectivity[0].neighbourIndex[0][COORDINATE::X] == 0 && "Unexpected owner index in X");
assert(interfaceConnectivity[0].neighbourIndex[1][COORDINATE::X] == 0 && "Unexpected owner index in X");
assert(interfaceConnectivity[0].neighbourIndex[2][COORDINATE::X] == 0 && "Unexpected owner index in X");
assert(interfaceConnectivity[0].neighbourIndex[3][COORDINATE::X] == 0 && "Unexpected owner index in X");
assert(interfaceConnectivity[0].neighbourIndex[4][COORDINATE::X] == 0 && "Unexpected owner index in X");
assert(interfaceConnectivity[0].neighbourIndex[0][COORDINATE::Y] == 0 && "Unexpected owner index in Y");
assert(interfaceConnectivity[0].neighbourIndex[1][COORDINATE::Y] == 1 && "Unexpected owner index in Y");
assert(interfaceConnectivity[0].neighbourIndex[2][COORDINATE::Y] == 2 && "Unexpected owner index in Y");
assert(interfaceConnectivity[0].neighbourIndex[3][COORDINATE::Y] == 3 && "Unexpected owner index in Y");
assert(interfaceConnectivity[0].neighbourIndex[4][COORDINATE::Y] == 4 && "Unexpected owner index in Y");
}Testing boundary condition reading
The boundary conditions are very similar to the interface reading. Instead of receiving vertices for the owner and neighbour, we only receive vertices for the boundary face itself, while we also receive information about the boundary condition type, which we will use to set values on the boundaries.
In the code below, we receive boundary conditions on line 3 and check that we have boundary conditions available for the 2 zones in our mesh on line 5. Then, for each zone, we check that we have exactly 3 boundary conditions on lines 6-7, followed by checking the type of each boundary condition for both zones on lines 9-15.
Similar to the interfaces, we receive the vertices which are located on each boundary. We check this information on lines 17-87 for each boundary. Since we have 3 boundaries per zone (so 6 in total for 2 zones), and per boundary 5 vertices, we have a few more lines to check here but the idea is the same as that for the interface.
For example, the first boundary in zone 0 is of type wall (line 9). We know from the sketch above that the bottom boundary is the wall boundary condition, and there is only a single wall boundary condition in that zone. Thus, we know that the starting point of that boundary is at i,j=0,0 and it will go up to i,j=4,0. This is what we are checking on lines 17-21 for the i index and lines 35-39 for the j index.
The boundaryCondition variable stores as its first index the zone, followed by the index storing the boundary conditions within that zone. Each boundary condition has a property called boundaryType and indices, as we can see from the code below.
void readStructuredMeshBoundaryConditions(ReadStructuredMesh &structuredMesh) {
// get the boundary conditions from the structured mesh
auto boundaryConditions = structuredMesh.getBoundaryConditions();
assert(boundaryConditions.size() == 2 && "Expected 2 zones");
assert(boundaryConditions[0].size() == 3 && "Expected to read 3 boundary conditions for the first zone");
assert(boundaryConditions[1].size() == 3 && "Expected to read 3 boundary conditions for the second zone");
assert(boundaryConditions[0][0].boundaryType == BC::WALL && "Expected wall boundary condition");
assert(boundaryConditions[0][1].boundaryType == BC::SYMMETRY && "Expected symmetry boundary condition");
assert(boundaryConditions[0][2].boundaryType == BC::INLET && "Expected inlet boundary condition");
assert(boundaryConditions[1][0].boundaryType == BC::WALL && "Expected wall boundary condition");
assert(boundaryConditions[1][1].boundaryType == BC::OUTLET && "Expected outlet boundary condition");
assert(boundaryConditions[1][2].boundaryType == BC::SYMMETRY && "Expected symmetry boundary condition");
assert(boundaryConditions[0][0].indices[0][COORDINATE::X] == 0 && "Unexpected index in X");
assert(boundaryConditions[0][0].indices[1][COORDINATE::X] == 1 && "Unexpected index in X");
assert(boundaryConditions[0][0].indices[2][COORDINATE::X] == 2 && "Unexpected index in X");
assert(boundaryConditions[0][0].indices[3][COORDINATE::X] == 3 && "Unexpected index in X");
assert(boundaryConditions[0][0].indices[4][COORDINATE::X] == 4 && "Unexpected index in X");
assert(boundaryConditions[0][1].indices[0][COORDINATE::X] == 0 && "Unexpected index in X");
assert(boundaryConditions[0][1].indices[1][COORDINATE::X] == 1 && "Unexpected index in X");
assert(boundaryConditions[0][1].indices[2][COORDINATE::X] == 2 && "Unexpected index in X");
assert(boundaryConditions[0][1].indices[3][COORDINATE::X] == 3 && "Unexpected index in X");
assert(boundaryConditions[0][1].indices[4][COORDINATE::X] == 4 && "Unexpected index in X");
assert(boundaryConditions[0][2].indices[0][COORDINATE::X] == 0 && "Unexpected index in X");
assert(boundaryConditions[0][2].indices[1][COORDINATE::X] == 0 && "Unexpected index in X");
assert(boundaryConditions[0][2].indices[2][COORDINATE::X] == 0 && "Unexpected index in X");
assert(boundaryConditions[0][2].indices[3][COORDINATE::X] == 0 && "Unexpected index in X");
assert(boundaryConditions[0][2].indices[4][COORDINATE::X] == 0 && "Unexpected index in X");
assert(boundaryConditions[0][0].indices[0][COORDINATE::Y] == 0 && "Unexpected index in Y");
assert(boundaryConditions[0][0].indices[1][COORDINATE::Y] == 0 && "Unexpected index in Y");
assert(boundaryConditions[0][0].indices[2][COORDINATE::Y] == 0 && "Unexpected index in Y");
assert(boundaryConditions[0][0].indices[3][COORDINATE::Y] == 0 && "Unexpected index in Y");
assert(boundaryConditions[0][0].indices[4][COORDINATE::Y] == 0 && "Unexpected index in Y");
assert(boundaryConditions[0][1].indices[0][COORDINATE::Y] == 4 && "Unexpected index in Y");
assert(boundaryConditions[0][1].indices[1][COORDINATE::Y] == 4 && "Unexpected index in Y");
assert(boundaryConditions[0][1].indices[2][COORDINATE::Y] == 4 && "Unexpected index in Y");
assert(boundaryConditions[0][1].indices[3][COORDINATE::Y] == 4 && "Unexpected index in Y");
assert(boundaryConditions[0][1].indices[4][COORDINATE::Y] == 4 && "Unexpected index in Y");
assert(boundaryConditions[0][2].indices[0][COORDINATE::Y] == 0 && "Unexpected index in Y");
assert(boundaryConditions[0][2].indices[1][COORDINATE::Y] == 1 && "Unexpected index in Y");
assert(boundaryConditions[0][2].indices[2][COORDINATE::Y] == 2 && "Unexpected index in Y");
assert(boundaryConditions[0][2].indices[3][COORDINATE::Y] == 3 && "Unexpected index in Y");
assert(boundaryConditions[0][2].indices[4][COORDINATE::Y] == 4 && "Unexpected index in Y");
assert(boundaryConditions[1][0].indices[0][COORDINATE::X] == 0 && "Unexpected index in X");
assert(boundaryConditions[1][0].indices[1][COORDINATE::X] == 1 && "Unexpected index in X");
assert(boundaryConditions[1][0].indices[2][COORDINATE::X] == 2 && "Unexpected index in X");
assert(boundaryConditions[1][0].indices[3][COORDINATE::X] == 3 && "Unexpected index in X");
assert(boundaryConditions[1][0].indices[4][COORDINATE::X] == 4 && "Unexpected index in X");
assert(boundaryConditions[1][1].indices[0][COORDINATE::X] == 4 && "Unexpected index in X");
assert(boundaryConditions[1][1].indices[1][COORDINATE::X] == 4 && "Unexpected index in X");
assert(boundaryConditions[1][1].indices[2][COORDINATE::X] == 4 && "Unexpected index in X");
assert(boundaryConditions[1][1].indices[3][COORDINATE::X] == 4 && "Unexpected index in X");
assert(boundaryConditions[1][1].indices[4][COORDINATE::X] == 4 && "Unexpected index in X");
assert(boundaryConditions[1][2].indices[0][COORDINATE::X] == 0 && "Unexpected index in X");
assert(boundaryConditions[1][2].indices[1][COORDINATE::X] == 1 && "Unexpected index in X");
assert(boundaryConditions[1][2].indices[2][COORDINATE::X] == 2 && "Unexpected index in X");
assert(boundaryConditions[1][2].indices[3][COORDINATE::X] == 3 && "Unexpected index in X");
assert(boundaryConditions[1][2].indices[4][COORDINATE::X] == 4 && "Unexpected index in X");
assert(boundaryConditions[1][0].indices[0][COORDINATE::Y] == 0 && "Unexpected index in Y");
assert(boundaryConditions[1][0].indices[1][COORDINATE::Y] == 0 && "Unexpected index in Y");
assert(boundaryConditions[1][0].indices[2][COORDINATE::Y] == 0 && "Unexpected index in Y");
assert(boundaryConditions[1][0].indices[3][COORDINATE::Y] == 0 && "Unexpected index in Y");
assert(boundaryConditions[1][0].indices[4][COORDINATE::Y] == 0 && "Unexpected index in Y");
assert(boundaryConditions[1][1].indices[0][COORDINATE::Y] == 0 && "Unexpected index in Y");
assert(boundaryConditions[1][1].indices[1][COORDINATE::Y] == 1 && "Unexpected index in Y");
assert(boundaryConditions[1][1].indices[2][COORDINATE::Y] == 2 && "Unexpected index in Y");
assert(boundaryConditions[1][1].indices[3][COORDINATE::Y] == 3 && "Unexpected index in Y");
assert(boundaryConditions[1][1].indices[4][COORDINATE::Y] == 4 && "Unexpected index in Y");
assert(boundaryConditions[1][2].indices[0][COORDINATE::Y] == 4 && "Unexpected index in Y");
assert(boundaryConditions[1][2].indices[1][COORDINATE::Y] == 4 && "Unexpected index in Y");
assert(boundaryConditions[1][2].indices[2][COORDINATE::Y] == 4 && "Unexpected index in Y");
assert(boundaryConditions[1][2].indices[3][COORDINATE::Y] == 4 && "Unexpected index in Y");
assert(boundaryConditions[1][2].indices[4][COORDINATE::Y] == 4 && "Unexpected index in Y");
}Testing unstructured mesh reading
Similar to the tests we have written for the structured grid above, we want to extend these tests to allow for testing of our unstructured grids. The first thing we want to do is to remind ourselves what the unstructured grid looks like, and below is the schematic for it. Similar to the structured grid, it has a bounding box of x=0, y=0 to x=2, y=1, that is, the index location for vertex v_1^1 on zone 1 is at x=0 and y=0, and the index location for vertex v_5^2 on zone 2 is at x=2 and y=1.
We store the coordinates for each vertex in a 1D array and then construct a helper array, the element connectivity array, that states which vertices combine to form a single cell. For example, cell C_5^1 shown above, i.e. the fifth cell in zone 1, contains the vertices v_3^1, v_4^1, and v_11^1. So, instead of just reading the coordinate information, we need to supplement this with the element connectivity information for the internal cells.
We need the information for the interfaces and boundary conditions, and these tests are largely the same as for the structured mesh tests. The difference here is that we are reading 1D element connectivity arrays for interface and boundary faces instead of i,j index pairs pointing us to the location of the vertices. We use the element connectivity of the interfaces and boundary faces to then get the coordinates through the coordinate arrays, since we know the vertex indices stored in the element connectivity. This will become clear once we create these tests. Let’s get started with the coordinate testing first.
Testing coordinate reading
The coordinate testing has many similarities to the structured mesh reading test. First, we receive the coordinates on line 3, and then check that we have read coordinates for two zones on line 6. The first zone should contain 13 vertices and the second zone 9 vertices, as we can see in the above sketch. This is asserted on lines 7-8.
The remaining section follows that of the structured mesh reading test, with the exception that we are not indexing our coordinate array with i,j indices, but instead using a 1D monotonically increasing index for each vertex ID. In other words, the first index of the coordinate array loops over all zones, the second index loops over all vertices, and the third index accesses either the x or y component of the coordinate.
Since we know the bounding box of the domain, we can read of the coordinates from the above sketch. After the mesh was generated in Pointwise, the vertices were moved so that coordinates would have at most 3 significant digits (e.g. 0.25, 0.75, 1.25, etc.).
We make use of the same trick as we did for the structured mesh reading test, where we subtract the expected coordinate value from the value we have read in the coordinate array. This should result in a value of zero if the values are the same, subject to some potential rounding-off errors. We take the absolute value of it and make sure that it is smaller than a specified threshold, e.g. here 0.01.
Let’s look at vertex v_11^1, for example. We expect it to have a value of x=0.75 and y=0.25. Since it has a zone index of 1 and a vertex index of 11, we need to transform that to C++-based zero indexing. Thus, we need to check for the x and y values for the coordinate array entry with zoneID=0 and vertexID=10, i.e. coordinates[0][10][COORDINATE::X] and coordinates[0][10][COORDINATE::Y], which we have on lines 20 and 34, respectively.
Other vertices follow the same pattern and from line 10 onwards, we test both the x and y coordinates in both zones and make sure that the values within the coordinate array match the expected values.
void readUnstructuredMeshCoordinates(ReadUnstructuredMesh &unstructuredMesh) {
// get the coordinates from the unstructured mesh
auto coordinates = unstructuredMesh.getCoordinates();
// test that coordinates were read correctly
assert(coordinates.size() == 2 && "Expected 2 zones");
assert(coordinates[0].size() == 13 && "Expected 13 vertices in zone 0");
assert(coordinates[1].size() == 9 && "Expected 9 vertices in zone 1");
assert(std::fabs(coordinates[0][0][COORDINATE::X] - 0.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][1][COORDINATE::X] - 0.50) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][2][COORDINATE::X] - 1.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][3][COORDINATE::X] - 1.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][4][COORDINATE::X] - 1.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][5][COORDINATE::X] - 0.50) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][6][COORDINATE::X] - 0.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][7][COORDINATE::X] - 0.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][8][COORDINATE::X] - 0.75) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][9][COORDINATE::X] - 0.50) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][10][COORDINATE::X] - 0.75) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][11][COORDINATE::X] - 0.25) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][12][COORDINATE::X] - 0.25) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[0][0][COORDINATE::Y] - 0.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][1][COORDINATE::Y] - 0.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][2][COORDINATE::Y] - 0.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][3][COORDINATE::Y] - 0.50) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][4][COORDINATE::Y] - 1.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][5][COORDINATE::Y] - 1.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][6][COORDINATE::Y] - 1.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][7][COORDINATE::Y] - 0.50) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][8][COORDINATE::Y] - 0.75) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][9][COORDINATE::Y] - 0.50) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][10][COORDINATE::Y] - 0.25) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][11][COORDINATE::Y] - 0.25) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[0][12][COORDINATE::Y] - 0.75) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][0][COORDINATE::X] - 1.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][1][COORDINATE::X] - 1.50) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][2][COORDINATE::X] - 2.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][3][COORDINATE::X] - 2.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][4][COORDINATE::X] - 2.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][5][COORDINATE::X] - 1.50) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][6][COORDINATE::X] - 1.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][7][COORDINATE::X] - 1.00) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][8][COORDINATE::X] - 1.50) < 0.01 && "Unexpected value for X-coordinate");
assert(std::fabs(coordinates[1][0][COORDINATE::Y] - 0.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][1][COORDINATE::Y] - 0.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][2][COORDINATE::Y] - 0.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][3][COORDINATE::Y] - 0.50) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][4][COORDINATE::Y] - 1.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][5][COORDINATE::Y] - 1.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][6][COORDINATE::Y] - 1.00) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][7][COORDINATE::Y] - 0.50) < 0.01 && "Unexpected value for Y-coordinate");
assert(std::fabs(coordinates[1][8][COORDINATE::Y] - 0.50) < 0.01 && "Unexpected value for Y-coordinate");
}Testing element connectivity reading
Once we have read the coordinates, we have to read the element connectivity data next. As alluded to above, the element connectivity data will allow us to map vertex indices to specific cells. The element connectivity test is shown below.
First, we get, as usual, the information from the mesh itself. On line 2, we receive the cells array, which is a 3D array. The first index loops over all zones, the second over all cells within that zone, and the third loops over all vertices within that cell. Thus, we first check that we have read two zones on line 4, and we check that we have 12 cells for the first, as well as 4 cells for the second zone on lines 6-7.
Lines 9-20 check that all cells in the first zone have the current vertex indices, and lines 22-25 do the same check for the second zone. For example, we said that the fifth cell in the first zone, i.e. C_5^1 consists of vertices v_3^1, v_4^1, and v_11^1, see the sketch of the mesh above. Accounting again for zero-based indexing, this means that cell 4 should have 3 vertices with ID 2, 3, and 10. We can see that this is what we are testing on line 13. Note that all vertices are checked on the same line.
We can see that cells 0-7 (or 1-8 using one-based indexing) only contain 3 vertices, i.e. these cells are triangular cells. See lines 9-16 for the first zone. However, there are also 4 quad elements in the first zones, and these we check on lines 17-20. We see that we check here for four vertex IDs. The second zone on lines 22-25 all contain 4 vertices and we know from the mesh shown above that we expect exactly four quad elements here.
void readUnstructuredInternalCells(ReadUnstructuredMesh &unstructuredMesh) {
auto cells = unstructuredMesh.getInternalCells();
assert(cells.size() == 2 && "Expected 2 zones");
assert(cells[0].size() == 12 && "Expected 12 cells in zone 0");
assert(cells[1].size() == 4 && "Expected 4 cells in zone 1");
assert(cells[0][0][0] == 2); assert(cells[0][0][1] == 10); assert(cells[0][0][2] == 1);
assert(cells[0][1][0] == 4); assert(cells[0][1][1] == 8); assert(cells[0][1][2] == 3);
assert(cells[0][2][0] == 12); assert(cells[0][2][1] == 5); assert(cells[0][2][2] == 6);
assert(cells[0][3][0] == 11); assert(cells[0][3][1] == 7); assert(cells[0][3][2] == 0);
assert(cells[0][4][0] == 2); assert(cells[0][4][1] == 3); assert(cells[0][4][2] == 10);
assert(cells[0][5][0] == 12); assert(cells[0][5][1] == 6); assert(cells[0][5][2] == 7);
assert(cells[0][6][0] == 4); assert(cells[0][6][1] == 5); assert(cells[0][6][2] == 8);
assert(cells[0][7][0] == 11); assert(cells[0][7][1] == 0); assert(cells[0][7][2] == 1);
assert(cells[0][8][0] == 8); assert(cells[0][8][1] == 5); assert(cells[0][8][2] == 12); assert(cells[0][8][3] == 9);
assert(cells[0][9][0] == 9); assert(cells[0][9][1] == 12); assert(cells[0][9][2] == 7); assert(cells[0][9][3] == 11);
assert(cells[0][10][0] == 9); assert(cells[0][10][1] == 11); assert(cells[0][10][2] == 1); assert(cells[0][10][3] == 10);
assert(cells[0][11][0] == 8); assert(cells[0][11][1] == 9); assert(cells[0][11][2] == 10); assert(cells[0][11][3] == 3);
assert(cells[1][0][0] == 5); assert(cells[1][0][1] == 6); assert(cells[1][0][2] == 7); assert(cells[1][0][3] == 8);
assert(cells[1][1][0] == 8); assert(cells[1][1][1] == 7); assert(cells[1][1][2] == 0); assert(cells[1][1][3] == 1);
assert(cells[1][2][0] == 8); assert(cells[1][2][1] == 1); assert(cells[1][2][2] == 2); assert(cells[1][2][3] == 3);
assert(cells[1][3][0] == 5); assert(cells[1][3][1] == 8); assert(cells[1][3][2] == 3); assert(cells[1][3][3] == 4);
}Testing interface reading
For the interface, the test is very similar to the structured mesh reading test, albeit with a slightly different structure. Interface information can be either written as a global node in a CGNS file or as a sub-node for each zone. Pointwise has decided to write structured mesh interfaces under a global node and unstructured interfaces under sub-nodes within zones. Thus, we retrieved global interface connectivity data when we tested structured mesh reading but now we get local interface information per zone.
This means we have to slightly change the way we interact with the interface information. We could have also decided to reconstruct the interface information in a consistent way within the class (which would have been probably a better idea, but one I probably overlooked when writing the code). Either way, for the unstructured grid, we now need to deal with interface information per zone.
We receive the interface information on line 3, and since we expect this information to be stored for each zone separately, we check that we have exactly two zones present for this interface information. The structure of the array is such that the first index loops over the zone, and then the second index over all interfaces within that zone. Thus, we are looking at a 2D array.
We expect exactly one interface per zone, which is checked on lines 7-8. We also check that the owning and neighbouring zones have indices 0 and 1 on lines 10-11. There are only two cells in the interface on both sides, which is checked on lines 13-14, and we can verify that with the sketch of the mesh above.
Next, we are checking the vertices on that interface. From the mesh above, we can see that the vertex IDs on the first zone are 3, 4, and 5. Using zero-based indexing, this would translate to 2, 3, and 4. If we wanted to describe the cells on that interface, then we need to prescribe the start and end index for each cell. For the first cell, this would be indices 2 and 3, and, for the second cell, this would be indices 3 and 4. This is what we check on lines 16-19.
The vertices on the second zone, that match the sequence of 2, 3 and 3, 4 are vertices 0, 7 and 7, 6. This is what we check on lines 21-24, which show the element connectivity of the cells on the interface of the neighbour zone. We then repeat the test for the second interface from line 26, which contains the same information. It has to, by definition, as the interfaces connect two zones and so both sides need to have the same information.
We can see that we are storing now the information twice, which may be redundant and wasting storage to some extent. However, the information is locally available within each zone, and so its information is readily available. It is a trade-off between having information where you need it vs. optimising storage, and this is often a call we have to make when writing CFD solvers.
void readUnstructuredMeshInterfaceConnectivity(ReadUnstructuredMesh &unstructuredMesh) {
// get the interface connectivity from the unstructured mesh
auto interfaceConnectivity = unstructuredMesh.getInterfaceConnectivity();
// test that interface connectivity was read correctly
assert(interfaceConnectivity.size() == 2 && "Expected 2 zones");
assert(interfaceConnectivity[0].size() == 1 && "Expected 1 interface in zone 0");
assert(interfaceConnectivity[1].size() == 1 && "Expected 1 interface in zone 1");
assert(interfaceConnectivity[0][0].zones[0] == 0 && "Expecting owning zone to be zone 0");
assert(interfaceConnectivity[0][0].zones[1] == 1 && "Expecting neighbouring zone to be zone 1");
assert(interfaceConnectivity[0][0].ownerIndex.size() == 2 && "Expected 2 owner indices");
assert(interfaceConnectivity[0][0].neighbourIndex.size() == 2 && "Expected 2 neighbour indices");
assert(interfaceConnectivity[0][0].ownerIndex[0][0] == 2 && "Unexpected vertex index");
assert(interfaceConnectivity[0][0].ownerIndex[0][1] == 3 && "Unexpected vertex index");
assert(interfaceConnectivity[0][0].ownerIndex[1][0] == 3 && "Unexpected vertex index");
assert(interfaceConnectivity[0][0].ownerIndex[1][1] == 4 && "Unexpected vertex index");
assert(interfaceConnectivity[0][0].neighbourIndex[0][0] == 0 && "Unexpected vertex index");
assert(interfaceConnectivity[0][0].neighbourIndex[0][1] == 7 && "Unexpected vertex index");
assert(interfaceConnectivity[0][0].neighbourIndex[1][0] == 7 && "Unexpected vertex index");
assert(interfaceConnectivity[0][0].neighbourIndex[1][1] == 6 && "Unexpected vertex index");
assert(interfaceConnectivity[1][0].zones[0] == 0 && "Expecting owning zone to be zone 0");
assert(interfaceConnectivity[1][0].zones[1] == 1 && "Expecting neighbouring zone to be zone 1");
assert(interfaceConnectivity[1][0].ownerIndex.size() == 2 && "Expected 2 owner indices");
assert(interfaceConnectivity[1][0].neighbourIndex.size() == 2 && "Expected 2 neighbour indices");
assert(interfaceConnectivity[1][0].ownerIndex[0][0] == 2 && "Unexpected vertex index");
assert(interfaceConnectivity[1][0].ownerIndex[0][1] == 3 && "Unexpected vertex index");
assert(interfaceConnectivity[1][0].ownerIndex[1][0] == 3 && "Unexpected vertex index");
assert(interfaceConnectivity[1][0].ownerIndex[1][1] == 4 && "Unexpected vertex index");
assert(interfaceConnectivity[1][0].neighbourIndex[0][0] == 0 && "Unexpected vertex index");
assert(interfaceConnectivity[1][0].neighbourIndex[0][1] == 7 && "Unexpected vertex index");
assert(interfaceConnectivity[1][0].neighbourIndex[1][0] == 7 && "Unexpected vertex index");
assert(interfaceConnectivity[1][0].neighbourIndex[1][1] == 6 && "Unexpected vertex index");
}Testing boundary condition reading
The final part is reading the boundary condition, and this is again very similar to the structured boundary condition reading. The only difference, again, is that we do not receive i,j index pairs for coordinates on the boundaries, but rather 1D indices that we can use to index the coordinate array to get the x and y coordinates on the boundaries.
On line 3, we receive the 2D boundary condition array and its composition is such that the first index loops over all zones, and the second over all boundaries within that zone. We check that we have boundary conditions available for both zones on line 6 and that we have 3 boundary conditions per zone on lines 8-9. Each boundary has 2 cells, so we check that each boundary meets this criterion on lines 11-17 and then we check that the correct type was read on lines 19-25, for each boundary condition.
From lines 27 onwards, we check that the correct vertices were read for each boundary condition. For example, we know that the first boundary for the second zone is of type outlet (line 23), so we would expect this boundary to be on the far right of the mesh sketch. We can see that the two cells making up the boundary cells on that boundary have vertices 3, 4 and 4, 5, see the sketch of the mesh above. With zero-based indexing, that becomes 2, 3 and 3, 4. This is what we check on lines 36-37, again with two asserts per line to save some space.
The remaining boundary cells follow a similar pattern and once we have read all boundary information, without any issue, we have finished testing our unstructured mesh reading.
One comment on the unstructured boundary condition reading; even though we can read meshes with boundary conditions stored either under family or boundary nodes within the CGNS file, we would not expect them to produce different results during the testing. The reason is that we want to deal with all of this complexity within the mesh reading class. The end user should not see any of that. The same is true for the structured mesh, where the boundary condition reading also just tests one implementation, and both should work fine without modifications.
void readUnstructuredMeshBoundaryConditions(ReadUnstructuredMesh &unstructuredMesh) {
// get the boundary conditions from the structured mesh
auto boundaryConditions = unstructuredMesh.getBoundaryConditions();
// test that boundary conditions were read correctly
assert(boundaryConditions.size() == 2 && "Expected 2 zones");
assert(boundaryConditions[0].size() == 3 && "Expected 3 boundaries in zone 0");
assert(boundaryConditions[1].size() == 3 && "Expected 3 boundaries in zone 1");
assert(boundaryConditions[0][0].indices.size() == 2 && "Expected 2 boundary indices");
assert(boundaryConditions[0][1].indices.size() == 2 && "Expected 2 boundary indices");
assert(boundaryConditions[0][2].indices.size() == 2 && "Expected 2 boundary indices");
assert(boundaryConditions[1][0].indices.size() == 2 && "Expected 2 boundary indices");
assert(boundaryConditions[1][1].indices.size() == 2 && "Expected 2 boundary indices");
assert(boundaryConditions[1][2].indices.size() == 2 && "Expected 2 boundary indices");
assert(boundaryConditions[0][0].boundaryType == BC::INLET && "Expected inlet boundary condition");
assert(boundaryConditions[0][1].boundaryType == BC::SYMMETRY && "Expected symmetry boundary condition");
assert(boundaryConditions[0][2].boundaryType == BC::WALL && "Expected wall boundary condition");
assert(boundaryConditions[1][0].boundaryType == BC::OUTLET && "Expected outlet boundary condition");
assert(boundaryConditions[1][1].boundaryType == BC::SYMMETRY && "Expected symmetry boundary condition");
assert(boundaryConditions[1][2].boundaryType == BC::WALL && "Expected wall boundary condition");
assert(boundaryConditions[0][0].indices[0][0] == 6); assert(boundaryConditions[0][0].indices[0][1] == 7);
assert(boundaryConditions[0][0].indices[1][0] == 7); assert(boundaryConditions[0][0].indices[1][1] == 0);
assert(boundaryConditions[0][1].indices[0][0] == 4); assert(boundaryConditions[0][1].indices[0][1] == 5);
assert(boundaryConditions[0][1].indices[1][0] == 5); assert(boundaryConditions[0][1].indices[1][1] == 6);
assert(boundaryConditions[0][2].indices[0][0] == 0); assert(boundaryConditions[0][2].indices[0][1] == 1);
assert(boundaryConditions[0][2].indices[1][0] == 1); assert(boundaryConditions[0][2].indices[1][1] == 2);
assert(boundaryConditions[1][0].indices[0][0] == 2); assert(boundaryConditions[1][0].indices[0][1] == 3);
assert(boundaryConditions[1][0].indices[1][0] == 3); assert(boundaryConditions[1][0].indices[1][1] == 4);
assert(boundaryConditions[1][1].indices[0][0] == 4); assert(boundaryConditions[1][1].indices[0][1] == 5);
assert(boundaryConditions[1][1].indices[1][0] == 5); assert(boundaryConditions[1][1].indices[1][1] == 6);
assert(boundaryConditions[1][2].indices[0][0] == 0); assert(boundaryConditions[1][2].indices[0][1] == 1);
assert(boundaryConditions[1][2].indices[1][0] == 1); assert(boundaryConditions[1][2].indices[1][1] == 2);
}Summary
This brings us to the end of not just this article, but also this series. We implemented mesh reading from scratch for both structured and unstructured grids, and we used the CGNS file format to help us achieve that. We turned the code we wrote into a library, which helps us to reuse our code later should we need to use it in another application. We also looked at basic code testing in this article and how we can ensure that the mesh reading is working as expected.
We covered quite some ground in this series but hopefully, you have an idea now how we can achieve mesh reading with the CGNS file format. As I mentioned at the beginning of this series, the CGNS file format comes with its own baggage and disadvantages, but it is the best and complete approach to storing grids and additional CFD data that we may want to exchange between applications. You should aim to understand the basics of CGNS, it will hold you in good stead.
Countless other formats exist that are easier to learn but they are also more restrictive. The CGNS format is the only one that comes to mind that can deal with both structured and unstructured grids. With its latest version of the library, we even have polyhedral elements support, and it does offer support for higher order elements, which is one of the very few formats that has the capabilities to do so.
You see, even if you don’t need these features right now, you may need them in the future, and if you already know how to work with the CGNS format, it’ll become much easier to extend your knowledge to new features than it is to learn CGNS from scratch and trying to get it to work. I hope you have learned something new in this series and are looking forward to the next one.
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.