Put on the kettle and get a cup of tea (or two), this is going to be a bit of a longer read. We will look at the standard template library, or STL, and we will look at a few example how we could use it in our every day CFD solver development process. The STL provides us with common data structures and algorithms that operate on the data structures, which can take a lot of boilerplate code out of your codebase. You never have to write a sorting algorithm again, and the good news is, that the algorithms and data structures that are provided, are implemented by C++ veterans that know how to write efficient and performant code. Chances are, that if you use the STL, you will notice a performance increase, especially if you used a lot of your own implementations that are already covered by the STL. In this article, then, I want to look at the most important data structures (or containers, as they are known in the STL) and algorithms, and then finally explain briefly how we glue both of these together with so-called iterators. The STL is yet another important tool we need to know, so if you are serious about C++ development and writing efficient CFD solver and tools, you need to know about the STL (and then, ideally use it).
In this series
- Part 1: Choosing the right programming language for CFD development
- Part 2: Why you should use C++ for CFD development
- Part 3: The complete guide to memory management in C++ for CFD
- Part 4: Object-orientated programming in CFD
- Part 5: How to handle inheritance and class hierarchies in C++
- Part 6: Templates in C++: Boost your CFD solver performance
- Part 7: Enhance readability with operator overloading in C++
- Part 8: The power of the standard template library (STL) in C++
- Part 9: Understanding Lambda expressions and how to use them in C++
- Part 10: Reduce memory bugs with smart pointers in C++
In this article
Introduction
Now then, with all of the knowledge that we have built up until now, you should be able to understand the standard template library, or STL for short. The STL is a collection of the most common data structures, algorithms, and iterators. We will look into each of them but suffice it to say at this stage, this is a powerful set of tools you should be using in your code, everywhere. Most people shy away from learning the STL, as some concepts in C++ are already difficult enough to master, but the STL helps us to write less code and make the code that we write more effective and optimised for speed. The STL is written by C++ veterans who understand memory management and who have written the code in a way that extracts the most performance for us. We would be fools not to make use of that. The only thing we need to get right is to choose the right data structure and algorithms for the task, but thankfully it is not as difficult as it may sound.
STL ingredient 1/3: Containers (data structures)
Data structures, or containers called in the STL, are the fundamental building blocks of any code. They hold data in a way that can be easily processed for a specific task while optimising memory layout and access. There are a handful of containers available in C++ and thankfully, we only really need to understand a few of them. I would argue that if you can master only 3 containers, you will be able to write efficient CFD solvers. In the following section, I want to shed some light on the most important containers that we should add to our programming arsenal and show how to use them in the context of real CFD solver development tasks.
Your new best friends: std::vector
We have already looked into std::vector and std::array throughout this series, but I wanted to spend more time on them and introduce some more features. We said that std::vector allocates its memory on the heap, while the std::array allocates memory on the stack. Thus, we use std::vectors for large arrays, such as coordinates x, y, and z, solution vectors for velocity, pressure, temperature, etc., while we use std::array for smaller arrays such as a normal vector, a transformation matrix, etc. std::array needs to know at compile time the size of the array, so we cannot use it in situations where memory needs to be dynamically allocated at runtime. std::vectors, on the other hand, can be resized to any length dynamically and so can be modified as we go along.
We could now look at the documentation of the std::vector class, but instead, let’s jump straight into code and see some of the most useful and used functions that the std::vector class exposes to us:
#include <iostream>
#include <vector>
int main()
{
// declaration
std::vector<double> vector1;
std::vector<double> vector2(5);
// managing memory
vector1.reserve(1000);
std::cout << vector1.size() << ", " << vector2.size() << std::endl;
std::cout << &vector1[0] << ", " << &vector2[0] << std::endl;
vector1.resize(50);
std::cout << vector1.size() << ", " << vector2.size() << std::endl;
std::cout << &vector1[0] << ", " << &vector2[0] << std::endl;
vector2.resize(65);
std::cout << vector1.size() << ", " << vector2.size() << std::endl;
std::cout << &vector1[0] << ", " << &vector2[0] << std::endl;
// adding data
vector1[0] = 3.1415;
vector1[1] = 1.618;
vector2.resize(0);
vector2.push_back(3.1415);
vector2.push_back(1.618);
vector2.pop_back();
vector2.shrink_to_fit();
return 0;
}The code above will produce the following output:
0, 5
0xa9b2e0, 0xa9b2b0
50, 5
0xa9b2e0, 0xa9b2b0
50, 65
0xa9b2e0, 0xa9e240Let’s look at the code and see what we are doing here. On lines 7-8, we declare two vectors, the first without any specific size, the second with space for 5 doubles. Then, on line 11, we are reserving space for 1000 doubles for vector1, but when we print the size of each vector on the next line, vector1 still has a size of 0 (see the first line of the output). reserve only reserves space, but does not actually allocate it. This is useful if we know that our vector will change in size frequently.
Say you are writing a CFD solver and you want to support adaptive mesh refinement. You know that your mesh will change over time, and so will your solution vector. So, you need to have some capacity to let your std::vector grow. Hence, we reserve some memory in advance. Only once we call resize, as we do on line 15, does the size of vector1 change. So you may ask yourself then, if I know that my size may change over time, why don’t I just allocate that much memory and just forget about it?
Well, memory allocation is costly, and guess what, at the beginning of this section I mentioned that the STL is written by people who do understand memory management and who can optimise our code to minimise performance overhead due to bad memory management. Well, that’s exactly what the std::vector is doing for us. When we call resize, the std::vector will actually allocate more memory than we need, but just a small amount. So next time, when you want to add a few elements, you don’t need to allocate more memory. There is a trade-off between how often you want to allocate memory and the size of your std::vector, but the way it is implemented is quite clever and better than what you and I probably can come up with ourselves.
OK, so we do understand the resize function now, but why bother with reserve? Well, if you know your std::vector will not change over time (e.g. no adaptive mesh refinement), then you can just make a single call to resize and be done with it, no need to first call reserve. This function simply blocks off memory for us, which we can use later for memory allocation. The reason that this is important is that a std::vector lays out its elements contiguously in memory. So, if we want to grow a vector, but there already is some data in the memory address just next to our last element in the vector, then we have to copy the whole vector into a different place in memory, and that is costly. We want to avoid that at all costs and so by reserving space, we ensure that neither our application nor the operating system or any other application can take that space away from us. We haven’t allocated yet, as this is time-consuming and wasteful of resources, but we can grow into it if we need to.
On line 19, we see that we can resize a vector that previously had already allocated some memory. But what we discussed above now may become a problem. vector2 was only initially allocated five doubles, now we request 65 doubles. If the next 480 bytes of memory are free (60 additional doubles, each of size 8 bytes), then we simply allocate memory and extend our vector, if it isn’t then it is first copied elsewhere and then more memory is allocated. And if we inspect the output of the program, we can see that this is exactly what is happening. By printing the memory address of the first entry in the vector using either &vector1[0] or &vector2[0], we see that after resizing vector2 on line 19, we see that the address of the first entry has changed, indicating that the entire vector was now copied to a different memory location. This is wasteful memory copying and so if we know that our vector will increase slowly, we should make sure that we call reserve first.
Accessing data is done via the usual square brackets, as seen on lines 24-25. We can use this to both access and retrieve data. This type of adding data is only possible, though, if we have called resize before, or, we have specified explicitly during the object creation how many elements we want to allocate memory for. If we do not have any size allocated, as we indicate by calling resize(0) on line 27, then we can’t add elements at specified locations and, instead, need to use the push_back() function, which will allocate memory for at least one element and add the value provided in parentheses as the value stored in the vector (again, the STL implementation is clever and adds a bit more than we need to make subsequent calls more efficient). Similarly, we can remove the last element using the pop_back() function, and, we can also call the shrink_to_fit() function to minimise memory usage for our vector, if we have frequently called resize or push_back and pop_back. If you find yourself using push_back() in your code and know the size of your vector, use reserve() first to minimise your memory allocation and copying costs.
At this point, we can’t avoid talking about cppreference. It is one of those websites I have a love-hate relationship with. I love it because it is up-to-date with all of the latest implementations and proposals by the C++ steering committee, which releases a new standard every 3 years that C++ compilers have to implement. It is the most complete and comprehensive resource available when it comes to C++ development. I hate it because half of the time I go through it I don’t understand what I am reading. I am proud to to say that I do understand at least half of what I am reading, which used to be a lot less and thus I did not frequent the website too often in the past. But as proficient C++ developers, we need to be comfortable enough to explore the website. It has a really good summary for us on all the functions of the different containers, and the first site I want you to take in is that of the std::vector, where you will find more information about the different functions it exposes to us. Most of the functions come with example code to show how to use them, but some examples are not always the clearest in my opinion.
std::vector’s little sister: The std::array
The std::array is very similar to the std::vector, with the already discussed main difference that the size of the array is declared at compile time. So there are no functions we can call to dynamically resize our array. To specify the size, you need to provide a second template argument, which is similar to how we would define an array in C:
#include <array>
int main() {
// array definition in C
double cArray[10];
// array definition in C++
std::array<double, 10> cppArray;
cArray[0] = 3.141592;
cppArray[0] = 3.141592;
std::cout << cArray[0] << ", " << cppArray[0] << std::endl;
return 0;
}Both cArray and cppArray can store the same amount of data and they can be accessed in the same way using the square brackets, as shown on lines 10 and 11. Remember that you can always replace a call to cppArray[0] = 3.141592; mentally with cppArray.operator[](0) = 3.141592;, at least in C++. We discussed this in our previous article on operator overloading, give it a read if you need a refresher.
In C++, you should always use std::array instead of the older C syntax to declare an array. If you do, you will be able to use all of the algorithms the STL provides us which we will look at in a second. The older C syntax of arrays requires you to implement any fundamental algorithm yourself from scratch, like sorting your data, for which the STL provides highly optimised codes already. Remember, we want to be DRY (Don’t repeat yourself) programmers, reusing existing code (even if it is not written by us) will help us to do just that.
Accessing data and assigning values is done in the same way as with std::vector, i.e. with the square brackets operator, and I would say that probably 90% of all the data structures you will be using are std::vectors and std::arrays. So if you understand all of what was discussed thus far in this article, you are pretty much halfway through understanding the STL and only need to master algorithms and iterators now, though iterators just need to be understood conceptually; you never really have to deal with them as they are hidden from you, but we are getting ahead of ourselves. Check out the std::array reference on cppreference to see what functions you can call. You will realise that it feels very similar to a std::vector, minus the memory management part, and so we won’t repeat the example code from above.
Collecting named parameters at their value: std::map and associates
Before continuing, it is worth pointing out that you can find an overview of all possible data structures (containers) on cppreference. Of all of these containers, the only other data structure you really ought to know are the std::map and std::unordered_map. The idea behind both of them is that you store a collection of so-called key-value pairs. Other languages such as Python may call it a dictionary. Essentially, for each key, which can be of any data type, you store a corresponding value. If you want to get mathematical, this is called a bijective mapping. When is this useful? Typically in situations where we want to store a named parameter (i.e. using a std::string or char to identify it by name) and you want to store any associated data with that parameter. If you know that you have more than one parameter that follows this scheme, then you can put all of them into a std::map, which allows you to grow this list of named parameters as your solver is executed. A good example is that of reading boundary conditions, i.e. you want to store the name of the boundary patch itself, and then the associated boundary condition, e.g. is it an inlet, an outlet, or a wall? This is shown in the following example:
#include <iostream>
#include <map>
#include <utility>
#include <string>
int main()
{
enum BC_TYPE {WALL = 0, INLET, OUTLET};
std::map<std::string, int> boundaryConditions;
boundaryConditions.emplace(std::make_pair(std::string("inlet-1"), BC_TYPE::INLET));
boundaryConditions.emplace(std::make_pair(std::string("inlet-2"), BC_TYPE::INLET));
boundaryConditions.emplace(std::make_pair(std::string("outlet"), BC_TYPE::OUTLET));
boundaryConditions.emplace(std::make_pair(std::string("aerofoil"), BC_TYPE::WALL));
boundaryConditions.emplace(std::make_pair(std::string("flap"), BC_TYPE::WALL));
return 0;
}In this example, we are storing the name of the boundary patch as the key, here as a std::string, and then we store the associated boundary condition as an int, which we specify as an enum type (which by default is of type int). Later, when we are applying our boundary conditions, and we want to know what type our “aerofoil” boundary patch has, then we can query the map, similar to a std::vector and std::array with the square bracket operator, i.e.:
std::cout << boundaryConditions["aerfoil"] << std::endl; // prints 0, i.e. first entry in enum BY_TYPE
std::cout << boundaryConditions["inlet-1"] << std::endl; // prints 1, i.e. second entry in enum BY_TYPE
std::cout << boundaryConditions["outlet"] << std::endl; // prints 2, i.e. third entry in enum BY_TYPEThis way we can store as many key-value pairs as we want. The only difference between a std::map and std::unordered_map is the way elements are stored, and thus, this impacts performance, i.e., the time it takes to insert, delete, and access data. A std::map is stored in an orderly fashion based on its keys, and this is important if you want to traverse the std::map ordered by its keys. Typically, we do not have this requirement if we are only using the std::map to look up parameters, such as the boundary condition type above. For this, the std::unordered_map is better suited, which sacrifices the ordering of the keys (which is now random, or, unordered) for speed in terms of inserting, deleting, and accessing data. The good news is, that you simply have to change the data type from std::map to std::unordered_map, and the rest of the code remains the same. Try this with the example above (and remember to include the right header, for std::unordered_map it is #include <unordered_map>).
Another good example is storing of simulation parameters, for example, the number of inner and outer iterations, which turbulence model to use, what schemes to use to discretize your equations in time and space, etc. We may then define a few different std::maps to hold ints, floats, and other data types, based on what we want to store. Again, an std::unordered_map would be most performant, as we would only ever access specific information here, rather than insisting on walking through the map ordered by its keys. It has to be said, though, that for smaller objects (i.e. those which are only storing boundary conditions and parameters), it is unlikely that the performance difference between std::map and std::unordered_map is going to be noticeable.
Using std::map for cleaner code design
At this point, let’s try to inject some good code design while talking about parameter storage. Why is it a good decision to store parameters in a std::map (or std::unordered_map), why don’t we just create a struct and declare variables in there that we want to set to a certain type? OK, let’s create that struct:
#include <string>
struct parameters {
int numberOfIterations = 1000;
float CFL = 0.7;
std::string meshFile = "mesh/airfoil.cgns";
};We can add more parameters as needed, but for this example, let’s stick to these three parameters, which have different types, i.e. int, float, and std::string. If we wanted to store parameters using a std::map, we need a bit more boilerplate code, as we need to support adding
#include <iostream>
#include <type_traits>
#include <string>
#include <map>
class parameters {
public:
template<typename Type>
void add(std::string name, Type data) {
if constexpr (std::is_same<Type, int>::value)
_intParameter.emplace(std::make_pair(name, data));
else if constexpr (std::is_same<Type, float>::value)
_floatParameter.emplace(std::make_pair(name, data));
else if constexpr (std::is_same<Type, std::string>::value)
_stringParameter.emplace(std::make_pair(name, data));
}
template<typename Type>
Type get(std::string name) {
if constexpr (std::is_same<Type, int>::value)
return _intParameter[name];
else if constexpr (std::is_same<Type, float>::value)
return _floatParameter[name];
else if constexpr (std::is_same<Type, std::string>::value)
return _stringParameter[name];
}
private:
std::map<std::string, int> _intParameter;
std::map<std::string, float> _floatParameter;
std::map<std::string, std::string> _stringParameter;
};
int main() {
parameters p;
p.add<int>("numberOfIterations", 1000);
p.add<float>("CFL", 0.7);
p.add<std::string>("meshFile", "mesh/airfoil.cgns");
std::cout << p.get<int>("numberOfIterations") << std::endl;
std::cout << p.get<float>("CFL") << std::endl;
std::cout << p.get<std::string>("meshFile") << std::endl;
return 0;
}We define the parameter class between lines 6 and 32, which contains two functions; add() and get(). We use these to add either an int, float, or std::string parameter, and we achieve this through the use of templates. We could have defined several variations of these functions, i.e. addInt(), addFloat(), addString(), etc., but we want to make use of our newfound knowledge on templates that we gained in our article on templates. Using the keyword constexpr, we can evaluate the if/else statement at compile time, which we need to do when dealing with templates (everything gets evaluated at compile time with templates!). The type_traits header helps us to compare types, rather than the data these types hold. We create three separate std::maps on line 29-31, each of which either stores an int, float, or std::string as the parameter type. The add() and get() function then checks which of these maps to query based on the type provided to the function. Finally, on lines 36-38 we see how we can now add data to this class, and on line 40-42 how to retrieve data.
So why is this example better than the simple struct we provided in the example before? After all, we created more code and had to tease our brain with templates, why spend this extra effort? Because the version with std::maps follows the open-closed principle, and the struct version doesn’t. In a nutshell, the open-closed principle advocates class designs that need to be written once and then can be extended without changing the code (the class is open for extension by additional data types, but closed to modifications, i.e. the need to change the underlying code). In our example, if we wanted to add a fourth parameter turbulenceModel with a value of "kEpsilon", then we can simply call our add() function in the std::map example, and we would add this parameter (extension). If we used our simple struct, however, we would need to change the class, and thus it wouldn’t be closed (modification).
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
STL ingredient 2/3: Algorithms
Right, so at this point you should feel at least acquainted with the most useful containers. Again, there are a few more, but the ones we looked at are the most commonly used containers for CFD development (or, in general, for situations where we deal with large amounts of data and need to operate on them frequently, as we do typically for any physics-based simulation). While having specific data structures in C++ is a nice addition compared to C, the standard template library does not stop there. One of the most powerful aspects of it (and C++ as a whole) are the available algorithms, and there are so many that it is difficult to keep track of all of them. For the most part, these are defined in the <algorithm> header file, and this is the one we need to include if we want to make use of them. I would suggest at this point to have a quick look at the available algorithms on cppreference to get an idea of what is already implemented and available for you to use. Chances are that you can use parts of these algorithms to create your own, more complex algorithms that are based on these already defined ones.
We spent some time familiarising ourselves with containers in the STL, and now we are ready to reap the awards for it. If we use a container, we all of a sudden have access to a lot of different algorithms that we can apply to it. This is why we don’t want to use C-style arrays as we saw in a previous example. The STL provides us with a lot of boilerplate code that is optimised, highly efficient, and there really isn’t a good excuse for us not to make use of it. So in the following, I want to look through some of the more common examples that you may find.
A gentle introduction: std::fill
The first example I want to look at is the std::fill algorithm. As the name (and documentation) suggests, it is used to fill containers with content. Let’s look at an example how to use it with a std::vector:
#include <iostream>
#include <vector>
#include <algorithm>
void printVector(const std::vector<double> &vector) {
for (const auto &entry : vector)
std::cout << entry << " ";
std::cout << std::endl;
}
int main() {
std::vector<double> vec(5);
std::fill(vec.begin(), vec.end(), 3.141592);
printVector(vec);
return 0;
}Notice how we now include the <algorithm> header on line 3, to make use of the STL algorithms. Lines 5-9 define a convenience function which we just use to print out the content of our vector. Within the main function on lines 11-16, we define first a std::vector with a size of 5 doubles, and then we use std::fill to set each value within that std::vector equal to 3.141592. If we print the content of the vector, as we do on line 14, we get:
3.14159 3.14159 3.14159 3.14159 3.14159 So what are the first two arguments to the std::fill function? Essentially we need to specify the range over which we want to fill our vector. In our case, we want to fill the entire vector, and each suitable container has a begin() and end() function implemented to indicate what is the first and last element in that data structure. The point is that the algorithms implemented in the STL are generic and want to cater for a broad range of data structures. You may want to write your own data structure (container) because you decide that the available ones do not cover your needs exactly. If you provide the two functions begin() and end(), then you’ll be able to use the std::fill algorithm on your data structure as well! (you also need to implement a forward iterator, but we haven’t looked at iterators yet, so let’s not get ahead of ourselves).
So I said the STL algorithms are generic and want to cover a large range of applications. Say, for example, you do not want to change the value of the first and last element in the std::vector, in this case, you could simply replace the call to std::fill in the following way:
std::fill(vec.begin() + 1, vec.end() - 1, 3.141592);For completeness, the code now prints:
0 3.14159 3.14159 3.14159 0 You see the first and the last elements are now zero because we did not fill these elements. They could have any value, really, as we did not explicitly initialise these values beforehand, and you should always make sure that all your values are initialised to a value you want, don’t rely on the compiler to fill an empty vector with zeros if you do not specify their value explicitly. If you don’t, you may get a non-zero number which can lead to code that is difficult to debug.
So, to recap, an algorithm in the STL typically requires you to provide a start and end point between which you want to apply the algorithm to the data of the underlying container. Some algorithms, like the std::fill algorithm, also require you to provide additional inputs, all of which you can find by looking up the definition for the algorithm on cppreference.
Get your memory in order: std::sort
Sorting should be self-explanatory, if we have a vector with different entries, sometimes we need to sort the vector in ascending or descending order. Take a simple example, we have a 1D mesh (i.e. just a line) and we apply uniform spacing to it:
#include <iostream>
#include <vector>
#include <algorithm>
void printVector(const std::vector<double> &vector) {
for (const auto &entry : vector)
std::cout << entry << " ";
std::cout << std::endl;
}
int main() {
int numberOfVertices = 11;
double domainStart = 0.0;
double domainEnd = 1.0;
double dx = (domainEnd - domainStart) / (numberOfVertices - 1);
std::vector<double> x(numberOfVertices);
for (int i = 0; i < numberOfVertices; ++i)
x[i] = i * dx;
printVector(x);
return 0;
}This is the world’s simplest mesh generator (probably)! We can ignore the definition again at the top on lines 5-9, which we simply use to inspect what values our vector holds. We define a few parameters about the mesh, i.e. the number of vertices we want to use, and the start and end of the domain. We use that to calculate the mesh spacing dx, after which we define the coordinate array x which will print its content on line 23 as follows:
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Now let’s say we want to implement adaptive mesh refinement, and we may want to add a few points around x=0.5. If we wanted to add the points 0.45, 0.475, 0.525 and 0.55, we would achieve this in C++, using a std::vector, in the following way:
x.push_back(0.475);
x.push_back(0.55);
x.push_back(0.45);
x.push_back(0.525);I have not used any particular order here, just to show that it doesn’t matter. If you include this code now on line 22 in the example code above, then you will now print the following content:
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.475 0.55 0.45 0.525All of the newly added entries are appended to the back of the vector. Surely, we would like to have some sorted range here and this is exactly where the std::sort algorithm comes in. It only requires you to specify the range in which you want to sort the vector, i.e. using the begin() and end() function again. There are no further arguments for the std::sort algorithm. If we add this call to our code after we have added our new entries, we end up with the following code (reproduced in full for completeness):
#include <iostream>
#include <vector>
#include <algorithm>
void printVector(const std::vector<double> &vector) {
for (const auto &entry : vector)
std::cout << entry << " ";
std::cout << std::endl;
}
int main() {
int numberOfVertices = 11;
double domainStart = 0.0;
double domainEnd = 1.0;
double dx = (domainEnd - domainStart) / (numberOfVertices - 1);
std::vector<double> x(numberOfVertices);
for (int i = 0; i < numberOfVertices; ++i)
x[i] = i * dx;
x.push_back(0.475);
x.push_back(0.55);
x.push_back(0.45);
x.push_back(0.525);
std::sort(x.begin(), x.end());
printVector(x);
return 0;
}We can see here the call to std::sort on line 28, which by now should look familiar and predictable. For completeness, line 30 will now produce the following output:
0 0.1 0.2 0.3 0.4 0.45 0.475 0.5 0.525 0.55 0.6 0.7 0.8 0.9 1There is not much more to it than that. Whenever you find yourself in the need to sort your std::vector, don’t ask ChatGPT for un-optimised code, use the STL std::fill algorithm and enjoy highly optimised code for free.
Turning everything on its head: std::reverse
Sometimes you find yourself in the need to reverse your data. This is particularly useful if you try to apply a function which can only work from the beginning of your data structure (container), but it can’t be used from the back of your container. The std::string type is one such example. It allows you to find the first character which is not of a specific value. When would this be useful? When you are reading a string from a text file and you want to get rid of any leading and trailing whitespace or tabs. Say, for example, you have the following code:
#include <iostream>
#include <string>
int main() {
std::string text("\t a string with leading and trailing whitespaces and tabs \t ");
std::cout << "\"" << text << "\"" << std::endl;
return 0;
}I am using here \” to escape the quotation mark so that we can see where the string starts and ends when we print it on line 6. When we do, we get the following output:
" a string with leading and trailing whitespaces and tabs "If we want to remove the whitespace, then we can use the inbuild function find_first_not_of() of the std::string class. This function will look for a set of characters specified as function arguments, and if the current character in the string is not equal to any of the provided arguments, then it has found the first entry not of that type and it will return this position in the string. The full example is shown in the following:
#include <iostream>
#include <string>
#include <algorithm>
void removeWhitespace(std::string &text) {
int start = text.find_first_not_of("\t ");
text = text.substr(start, text.size() - start);
std::cout << "\"" << text << "\"" << std::endl;
}
int main() {
std::string text("\t a string with leading and trailing whitespaces and tabs \t ");
std::cout << "\"" << text << "\"" << std::endl;
removeWhitespace(text);
std::reverse(text.begin(), text.end());
removeWhitespace(text);
std::reverse(text.begin(), text.end());
std::cout << "\"" << text << "\"" << std::endl;
return 0;
}This code will print the following output:
" a string with leading and trailing whitespaces and tabs "
"a string with leading and trailing whitespaces and tabs "
"sbat dna secapsetihw gniliart dna gnidael htiw gnirts a"
"a string with leading and trailing whitespaces and tabs"In the removeWhitespace() function, we find the first character which is not a space a tab (\t), which will return the position in the string for which character that is the case here. Then, we create a substring of the original one, where we start the string at the position which does not contain any whitespace anymore. The length of the substring, i.e. the second parameter, is the length of the original string minus the length at the front that we are chopping off.
On line 15, we remove the trailing whitespace and at the end of the function on line 8 we print the new line, which produces the second print statement in the output above. The leading whitespace is indeed gone. Then, we reverse the string on line 16 using std::reverse, so that the trailing whitespace becomes a leading whitespace, we can then remove it with the same function, and finally reverse the string once more to have it in the correct order. We can see that the second call to removeWhitespace() does indeed remove the trailing whitespace, but the string is still in reverse order, as seen on the third output line. The second call call to std::reverse on line 19 bring the string back in its correct form, and the final print statement on line 21 confirms this, i.e. the fourth and final line of the output shows the string in the correct order without any leading or trailing whitespace.
When is this useful? Consider the following text file:
1 4 1 2 3 5 5 5 5 1 1
2 4 5 3 3 6 6 6 6 1 1
3 9 8 7 7 10 11 12 12 1 1
4 4 1 9 7 3 2 10 12 1 1There are leading and trailing whitespace in this text file. This represents part of an unstructured mesh file, using the StarCD format. If we wanted to read this file, we would first get rid of the whitespace, then split the string into individual components (which may then have either leading or trailing whitespace), and then remove the whitespace from them again. So this is useful whenever we are dealing with ASCII-based text files, and we want to remove spaces between entries.
Clone yourself: std::copy
std::copy is probably self-explanatory, but I just wanted to show here how to use it, especially for partial copies. So let’s look at an example straight away:
#include <iostream>
#include <vector>
#include <algorithm>
void printVector(const std::vector<int> &vector) {
for (const auto &entry : vector)
std::cout << entry << " ";
std::cout << std::endl;
}
int main() {
std::vector<int> vec = {0, 1, 2, 3, 4, 5};
std::vector<int> vec_copy(vec.size());
printVector(vec);
printVector(vec_copy);
std::cout << std::endl;
// copy entire vector
std::copy(vec.begin(), vec.end(), vec_copy.begin());
printVector(vec);
printVector(vec_copy);
std::cout << std::endl;
// reset vec_copy to zero everywhere
std::fill(vec_copy.begin(), vec_copy.end(), 0);
// copy vector partially, ignore first and last element
std::copy(vec.begin() + 1, vec.end() - 1, vec_copy.begin() + 1);
printVector(vec);
printVector(vec_copy);
std::cout << std::endl;
return 0;
}The above code will produce the following output:
0 1 2 3 4 5
0 0 0 0 0 0
0 1 2 3 4 5
0 1 2 3 4 5
0 1 2 3 4 5
0 1 2 3 4 0We have made use here of the printVector() function on lines 5-9 that we saw before, except now we are dealing with std::vectors of type int. We define two vectors on lines 12 and 13, where we use the aggregate initialisation to specify the entries in the first vector (i.e. we set the right-hand side equal to a list going from 0 to 5 and this will become the content of the vector. Initially, the vec_copy does not contain any values, but when we use std::copy on line 20 and then print the content again on lines 22 and 23, we see that the copied vector now contains the same values as the original vector. At this point I should point out that if this is what you want to achieve, then you could have also simply written vec_copy = vec;. Why? Because the std::vector has a function implemented that looks something like the following (the exact implementation will depend on the compiler):
template<typename Type>
std::vector<Type> operator=(const std::vector<Type> &rhs) {
std::vector<Type> copiedVector;
std::copy(rhs.begin(), rhs.end(), copiedVector.begin());
return copiedVector;
}I.e. the std::vector will have an overloaded operator= which works for any type the vector is defined for as we use templates here. If this looks unfamiliar, then you may want to have a look at my previous articles on operator overloading and templates to familiarise yourself with the concepts. If it does look familiar, then you will have spotted the std::copy buried on line 4, so this would achieve exactly the same thing.
Back to our previous example, we can reset the vec_copy back to having zero elements everywhere, as we do on line 27 using std::fill what we already saw. We are also able to partially copy a vector, as we do on line 30, where we ignore the first and last elements. It uses a similar syntax to what we used when we discussed std::fill, so hopefully nothing surprising here for you.
We make copies all the time, so there is not a particular CFD example that would serve as a good case study. But more complex algorithms may be composed of several calls to the STL, of which std::copy may be one of them.
Shed some light on your data: std::find
So this is a pretty common use case. If you have a std::vector filled with entries, you may want to figure out at what position a certain entry is located so that you can index the std::vector later at this index. Say, for example, you have a std::vector, that stores ints, which represent the face ID of boundary patches. You want to find a specific face for which you know the ID, but you don’t know where in the vector it is, exactly (for example, if we use an unstructured grid). In that case, we can use std::find to figure out where this face is located within the face ID vector. We can then later use that information, to access another vector wich stores the actual information about the face, i.e. the normal vector, integration point, surface area, etc., which will be located at the same index that we just found in the face ID vector.
To understand the usage, consider the following simplified example:
#include <iostream>
#include <vector>
#include <algorithm>
int main() {
std::vector<int> vec = {0, 7, 3, 5, 1, 11};
int value = 7;
auto iterator = std::find(vec.begin(), vec.end(), value);
auto index = std::distance(vec.begin(), iterator);
std::cout << value << " is located at index " << index << std::endl;
value = 12; // does not exist!
iterator = std::find(vec.begin(), vec.end(), value);
index = std::distance(vec.begin(), iterator);
std::cout << value << " is located at index " << index << std::endl;
// check if found value is actually within the vector
if (iterator != vec.end())
std::cout << value << " is located at index " << index << std::endl;
else
std::cout << value << " is not located in vecor" << std::endl;
return 0;
}This will produce the following output:
7 is located at index 1
12 is located at index 6
12 is not located in vecorThis demonstrates pretty much all there is to std::find. We define a vector on line 6 which does not have its entries ordered. Then we specify a value we want to find in the vector, here 7. On line 9, we use std::find to get an iterator to where the element is stored. We have not looked into iterators (they will be the final ingredient in our STL discussion and we will get to them shortly), but suffice it to say they provide us with some positional information as to where we are in the vector. On line 10 we use the std::distance algorithm which calculates the distance between the iterator (i.e. the position in the vector) and the beginning of the vector, which will return the index for us. We print that on line 11 and is shown as the first line of output above. Remember that C++ is zero-based with its indices, so if we say that the value 7 is at index 1, it is the second entry in the std::vector, the first entry would be at index 0.
If we repeat this step for an element that does not exist in the vector, as we do on line 13 where we now look for an element of value 12 (which is not in the vector), repeat the steps above and then print the result (lines 15-17), we still get something printed. However, if you look carefully, index 6 is not part of the vector, i.e. the last index in the vector is 5. C++ says that if an element is not within the vector, the index should be one index past the end of the vector. This is not just the case for std::find, but any other algorithm where you are looking for data. To overcome this issue, we need to check if the iterator is not equal to the end() of our vector, as we do on line 20. This line may seem strange, but the end() of a container is defined as the element which is one index past the last index in the container. So here, the last index is 5 in our example, so end() would point to index 6. So, as long as our iterator is not equal to the end() of the container, we have found an element, otherwise, we haven’t, and so the if/else construct of lines 20 to 23 shows how to check for this in code. The corresponding final print statement in the output shown above correctly states that the value 12 is not contained in the std::vector.
At this point, you have probably gotten an idea of how to use the STL algorithms, they typically follow the routine in this section, but they are actually much more powerful than what we saw this far. But for that, we need to first understand lambda functions, which we will cover later (and then return to the STL algorithms). The final ingredient is iterators which we have already seen a few times, see let’s look into them next.
STL Ingredient 3/3: Iterators
We have come a long way, we now know about the different data structures (containers) that are available in C++, and we have also looked at how to apply certain algorithms to them. But the reality is that not all data structures can be used with all algorithms. So how does C++ know which can be used and which can’t? In short: Iterators.
An iterator is a special type we have already come across a few times in the example above because it was unavoidable. Iterators provide some mechanism to step through our data structure. In the case of the std::vector or std::array, this is simple enough. We saw that by using an algorithm such as std::copy or std::fill, we specify a start and end location, and then C++ magically know how to iterate over all elements and fill or copy them. If we wanted to access the next element in a std::vector or std::array, we know that memory is laid out contiguously, so simply going to the next memory address adjacent to the current one allows us to step through (or iterate over) our containers. For other containers, such as std::unordered_map, this is not necessarily possible, as it does not have any specific order, and so we can’t say what element is next.
The idea is the following then: Define a couple of iterators that allow to access the data in a specific order, for example, going through a data structure (container) from start to end, or in reverse order from end to start. Then, for each data structure, provide mechanisms to step through your containers with all possible iterators. If a certain container can’t possibly implement a certain iterator, then all algorithms that rely on this particular iterator to step over data can’t be used on this container. So which iterators do we have then in C++? There are 5, and they are given below:
- Input Iterator
- Output Iterator
- Forward Iterator
- Bidirectional Iterator
- Random-access iterator
Let’s look briefly at them. The input and output iterator are pretty much self-explanatory; they allow us to read and write data at a specific location in the data structure. An iterator is nothing else than a pointer with special functionalities, so the input and output iterator point at a specific address in memory (which is associated with the data structure), and then it has the functionality to either read or write data. The forward iterator allows us to go through our data structure from start to end (and we have seen this plenty of times in action now, e.g. using std::fill or std::copy, where we started at begin() and iterated until we reached end() of a container), while the bidirectional iterator also allows us to go from the back to the beginning of the data structure as well as forward. Finally, the random-access iterator allows us to arbitrarily access values at any point directly. This is required for data structures where we can’t iterate over the data with a forward or bidirectional iterator, e.g. std::unordered_map.
The std::vector has implementations for all of these 5 iterators, and thus it will work with any algorithm provided by the STL. But let’s introduce another data structure, the linked list. It is not very important in CFD applications, so we don’t tend to use it that much, also because it is very similar to a std::vector but less performant. A linked list can store any number of values and like a std::vector grow in size. The difference is in its memory layout; while a std::vector requires contiguous allocation of memory addresses, a linked list can allocate memory anywhere in memory and thus never has to be copied, even when resizing the list. This is advantageous in applications where we know our data structure is always growing in size (outside of CFD, think of a comment function on social media, the comments are never decreasing but always increasing (ignoring deletion here), so you don’t want to be copying your data all the time when you are resizing your list (adding comments)). The way a list then is able to access its data is to store the address of the first entry and then for each entry, store the value and the pointer to the next element. This means we can use a forward iterator, but not the random access iterator, because we always need to go through all entries in order to find the next entry. If we also store a pointer back to the last element, then we can also use a bidirectional iterator.
Let’s look at some iterators in code to see how they would look like. The good news is, we are going to look at them here once, and then we will forget about them because the way that they are implemented in C++ means that we never explicitly have to use them, unless we are writing our own data structure class, but even then, we can reuse existing ones from other data structures.
#include <iostream>
#include <vector>
int main() {
std::vector<int> vec = {1, 2, 3, 4, 5};
// forward iterator
for (std::vector<int>::iterator it = vec.begin(); it != vec.end(); ++it) {
std::cout << *it << " ";
}
std::cout << std::endl;
// backward iterator
for (std::vector<int>::iterator it = vec.end() - 1; it != vec.begin() - 1; --it) {
std::cout << *it << " ";
}
std::cout << std::endl;
// forward iterator, normal use case
for (const auto &element : vec)
std::cout << element << " ";
std::cout << std::endl;
return 0;
}We define a vector here again on line 5 and then loop over it on line 7 using a forward iterator. Why is it forward? You can see that the first statement is a call to the iterator, which is initialized to vec.begin(). We say we want to iterate over this vector until we reach the end of the vector, as indicated by the vec.end() statement. Finally, we increment the iterator, using the ++it statement. Since we increment the iterator, it is a forward iterator. Remember, an iterator is just a pointer to the next value, so incrementing it will point to the next value in memory. Also, on line 8, we are printing the iterator. Since it is a pointer, we have to dereference it with the leading asterisk, i.e. *.
On line 13, we are using a backwards (bidirectional) iterator to loop from the back of the vector to the beginning of it. Notice how the final argument is --it, indicating a decrement in the iterator. But I also mentioned that we basically can ignore iterators altogether. The range-based for loop is one such example, and this is shown on lines 19-20. Here, we are looping over all the elements that do exist within vec. Notice how we do not need to dereference the element on line 20; in this case, we already get the dereferenced value back, which is more convenient and usually what we want. Range-based for loops have only two arguments; the first is the dereferenced value, and if we don’t care about its type, we can simply use auto here to automatically detect the type (the compiler will do this for us). In this example, we are getting a const reference, i.e. we are not copying the value into element but rather element is pointing to the value stored within the vector and we can use it directly. The second argument is the container over which we want to iterate. Simple enough, if you don’t care about the index, ranged-based for loops are the way to go in C++, all thanks to iterators.
Finally, one last word on iterators. Whenever you see calls to begin() and end(), an iterator is used in the background to loop over your data structure. If you are using the STL (and hopefully by now you are convinced that you need to), you will see a lot of calls to these two functions; whenever you do, now you know that in the background one of the 5 iterators is being used to walk over your data in an abstract manner (i.e. a manner which can be done for any data structure that implements these iterators).
Summary
You made it, this time we looked at the STL and we did so in much more detail. But we needed to look at three different concepts in order to understand the STL, and there is no point in cutting this article into three pieces. Let’s review the three main building blocks in the STL that help us to write efficient, highly-optimised code:
Containers: The are the underlying data structures that are defined in C++ and ready for us to use. You will find most classical data structures available in the STL so looking through the available data structures and reading up on when to use which is a good idea.
Algorithms: Just storing data is important, but we are not gaining anything if we do not operate on the data. Algorithms provide us with a way to perform common tasks, irrespective of the underlying data structure (container).
Iterators: These act as the glue between containers and algorithms, and provide a mechanism for algorithms to loop over containers. Different algorithms require different types of access to data, so long as the relevant iterator is implemented within a container class, we are able to use the algorithm, regardless of the underlying container.
And, the STL does not stop there. You can extend the STL by providing your very own container and then make use of all of the algorithms implemented in the STL (assuming you provide implementations for all necessary iterators). For CFD applications we probably want to define our own scalarField or vectorField class, where it may make sense to provide a definition for the operator+, operator-, and operator* to be able to add, subtract, and multiply scalar or vector fields together. If we wanted to parallelise our code, we could hide the parallel implementation within these overloaded operators and still have the same simple syntax of simply adding our data structures together. And we could make this much more complex (e.g. providing a class template which provides the parallel implementation, for example, if we wanted to provide parallelisation both on the CPU and GPU), but the syntax would always remain the same. We could still make use of the STL algorithm and enjoy their advantages, while doing all our case specific implementations in our custom container.
Be brave, use the STL. I waited far too long to embrace it, these days I can’t program without it anymore. It is useful, it is highly optimised, and neither you or me could write data structures or algorithms that are as performant as they are in the STL!
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.