created retentive array datastructure

lib/CMakeLists.txt

@@ -1,7 +1,7 @@
 
 add_library(nm_utils STATIC
 	"Window.cpp"
-)
+ "RetentiveArray.hpp")
 
 target_include_directories(nm_utils PUBLIC ${SDL2_INCLUDE_DIRS} ${CMAKE_CURRENT_SOURCE_DIR})
 target_link_libraries(nm_utils PUBLIC ${SDL2_LIBRARIES})

lib/RetentiveArray.hpp

@@ -0,0 +1,180 @@
+#pragma once
+
+#include <vector>
+#include <array>
+#include <functional>
+
+/**
+* @brief An array that remembers its past.
+* 
+ * A retentive array is a kind of buffered array that "remembers"
+ * the state of the array from the last n generations.
+ * This is extremely useful for numerically solving PDEs, as the values
+ * of the current array elements usually depend on the previous few states
+ * of the simulation.
+ * 
+ * This structure handles the evolution of these kinds of arrays and takes
+ * care of the memory allocation/freeing as well as properly swapping the 
+ * arrays.
+ * 
+ * @tparam Type The type of the elements in the array
+ * @tparam AttentionSpan How many generations of the data will be remembered by the array
+ */
+template<
+	typename Type, 
+	unsigned int AttentionSpan, 
+	typename std::enable_if_t<(AttentionSpan > 0), bool> = true>
+class RetentiveArray
+{
+public:
+	/**
+	 * @brief Constructs a new, empty retentive array 
+	 */
+	RetentiveArray()
+	{
+		// Do nothing
+	}
+
+	/**
+	 * @brief Constructs a new, empty retentive array
+	 * 
+	 * @param size The size of the data contained in the array
+	 */
+	RetentiveArray(size_t size)
+	{
+		// Create new vectors for every generation that needs to be remembered
+		for (int n = 0; n <= AttentionSpan; n++)
+		{
+			data[n] = std::vector<Type>(size, Type());
+		}
+	}
+
+	/**
+	 * @brief Constructs a new retentive array with the data of another array
+	 *
+	 * @param other The retentive array to copy from
+	 */
+	RetentiveArray<Type, AttentionSpan>& operator=(const RetentiveArray<Type, AttentionSpan>& other)
+	{
+		for (int n = 0; n <= AttentionSpan; n++)
+		{
+			data[n] = std::vector<Type>(other.data[n]);
+		}
+
+		return *this;
+	}
+
+	/**
+	 * @brief Constructs a new retentive array pointing to the data of an rvalue
+	 *
+	 * @param other The retentive array to set its data to
+	 */
+	RetentiveArray<Type, AttentionSpan>& operator=(RetentiveArray<Type, AttentionSpan>&& other)
+	{
+		for (int n = 0; n <= AttentionSpan; n++)
+		{
+			data[n] = std::vector<Type>(std::move(other.data[n]));
+		}
+
+		return *this;
+	}
+
+	/**
+	 * @brief Constructs a new retentive array with the data of another array
+	 *
+	 * @param other The retentive array to copy from
+	 */
+	RetentiveArray(const RetentiveArray<Type, AttentionSpan>& other)
+	{
+		*this = other;
+	}
+
+	/**
+	 * @brief Constructs a new retentive array pointing to the data of an rvalue
+	 *
+	 * @param other The retentive array to set its data to
+	 */
+	RetentiveArray(RetentiveArray<Type, AttentionSpan>&& other)
+	{
+		*this = other;
+	}
+
+	~RetentiveArray()
+	{
+		// Do nothing
+	}
+
+	/**
+	 * @brief Get the array from `index` generations ago
+	 * 
+	 * @param index Amount of generations to go backwards in time
+	 * @return The array from before `index` generations
+	 */
+	std::vector<Type>& operator[](size_t index)
+	{
+		return data[index];
+	}
+
+	/**
+	 * @brief Get the most up-to-date array
+	 * 
+	 * @return The array with generation 0
+	 */
+	std::vector<Type>& Current()
+	{
+		return data[0];
+	}
+
+	/**
+	 * @brief Evolve the data in the array
+	 * 
+	 * Simply calls the member function pointer. The called function is then
+	 * supposed to perform whatever is needed to compute the new array contents.
+	 * 
+	 * The data is swapped BEFORE calling this function, so the evolution function should
+	 * modify the data in the Current() vector (index 0)
+	 */
+	void Evolve()
+	{
+		// Evolve the array
+		Evolve(rule);
+	}
+
+	/**
+	 * @brief Evolve the data in the array
+	 *
+	 * Simply calls the provided function pointer. The called function is then
+	 * supposed to perform whatever is needed to compute the new array contents.
+	 *
+	 * The data is swapped BEFORE calling this function, so the evolution function should
+	 * modify the data in the Current() vector (index 0)
+	 * 
+	 * @param rule A function that evolves the data in the array
+	 */
+	void Evolve(std::function<void(void)> rule)
+	{
+		// Cycle the array contents, so that the previous "current" array is now the 2nd. 
+		// The 2nd becomes the 3rd etc, and the former last array becomes the new current
+		for (int n = 1; n <= AttentionSpan; n++)
+		{
+			data[0].swap(data[1]);
+		}
+
+		rule();
+	}
+
+	/**
+	 * @brief Set a general rule for evolving the data
+	 * 
+	 * @param rule The new evolution rule
+	 */
+	void SetEvolutionRule(std::function<void(void)> rule)
+	{
+		this->rule = rule;
+	}
+
+
+private:
+	std::array<std::vector<Type>, AttentionSpan + 1> data;
+	std::function<void(void)> rule = std::bind([]() {});	// Do nothing by default
+};

src/EulerFluid/FluidField.cpp

@@ -7,11 +7,12 @@
 #define VALUE(arr, x, y) ((arr)[(y) * this->size + (x)])
 #define PVALUE(arr, x, y) ((*(arr))[(y) * this->size + (x)])
 
+#define IDX(x, y, w) ((y) * (w) + (x))
+
 FluidField::FluidField(int size) :
 	size(size + 2)
 {
-	density		= new std::vector<double>(this->size * this->size, 0.0);
-	prevDensity = new std::vector<double>(this->size * this->size, 0.0);
+	density = RetentiveArray<double, 1>(this->size * this->size);
 
 	std::vector<double> hori(this->size * this->size);
 	std::vector<double> vert(this->size * this->size);
@@ -38,17 +39,14 @@ FluidField::FluidField(int size) :
 
 FluidField::~FluidField()
 {
-	delete prevDensity;
-	delete density;
-
 	delete prevVel;
 	delete vel;
 }
 
 void FluidField::AddSource(int x, int y, double dens, double dt)
 {
-	PVALUE(density, x, y) += dt * dens;
-	PVALUE(density, x, y) = std::max(PVALUE(density, x, y), 0.0);
+	density.Current()[IDX(x, y, size)] = dt * dens;
+	density.Current()[IDX(x, y, size)] = std::max(density[0][IDX(x, y, size)], 0.0);
 }
 
 void FluidField::AddFlow(int x, int y, double dx, double dy, double dt)
@@ -86,11 +84,11 @@ void FluidField::Diffuse(double diff, double dt)
 		{
 			for (int j = 1; j <= N; j++)
 			{
-				PVALUE(density, i, j) = (PVALUE(prevDensity, i, j) + a * (PVALUE(density, i - 1, j) + PVALUE(density, i + 1, j) + PVALUE(density, i, j - 1) + PVALUE(density, i, j + 1))) / (1 + 4 * a);
+				density[0][IDX(i, j, size)] = (density[1][IDX(i, j, size)] + a * (density[0][IDX(i - 1, j, size)] + density[0][IDX(i + 1, j, size)] + density[0][IDX(i, j - 1, size)] + density[0][IDX(i, j + 1, size)])) / (1 + 4 * a);
 			}
 		}
 
-		ApplyBoundaryConditions(BoundaryCondition::Continuous, *density);
+		ApplyBoundaryConditions(BoundaryCondition::Continuous, density[0]);
 	}
 }
 
@@ -121,12 +119,12 @@ void FluidField::Advect(double dt)
 			double t1 = y - j0;
 			double t0 = 1 - t1;
 
-			PVALUE(density, i, j) = s0 * (t0 * PVALUE(prevDensity, i0, j0) + t1 * PVALUE(prevDensity, i0, j1)) +
-									s1 * (t0 * PVALUE(prevDensity, i1, j0) + t1 * PVALUE(prevDensity, i1, j1));
+			density.Current()[IDX(i, j, size)] = s0 * (t0 * density[1][IDX(i0, j0, size)] + t1 * density[1][IDX(i0, j1, size)]) +
+									s1 * (t0 * density[1][IDX(i1, j0, size)] + t1 * density[1][IDX(i1, j1, size)]);
 		}
 	}
 
-	ApplyBoundaryConditions(BoundaryCondition::Continuous, *density);
+	ApplyBoundaryConditions(BoundaryCondition::Continuous, density[0]);
 }
 
 void FluidField::DiffuseVelocity(double visc, double dt)
@@ -276,10 +274,8 @@ void FluidField::DensityStep(double diff, double dt)
 		if(dx > 0 && dx < this->size - 1 && dy > 0 && dy < this->size - 1)
 			AddSource(dx, dy, 100.0, dt);
 
-	std::swap(prevDensity, density);	
-	Diffuse(diff, dt);
-	std::swap(prevDensity, density);
-	Advect(dt);
+	density.Evolve(std::bind(&FluidField::Diffuse, this, diff, dt));
+	density.Evolve(std::bind(&FluidField::Advect, this, dt));
 }
 
 void FluidField::Draw(SDL_Renderer* renderer, const SDL_Rect& target)
@@ -295,7 +291,7 @@ void FluidField::Draw(SDL_Renderer* renderer, const SDL_Rect& target)
 	{
 		for (int x = 0; x < this->size; x++)
 		{
-			double densityVal = std::min(PVALUE(density, x, y), 1.0);
+			double densityVal = std::min(density.Current()[IDX(x, y, size)], 1.0);
 			SDL_SetRenderDrawColor(renderer, densityVal * 255, densityVal * 255, densityVal * 255, 255);
 
 			vectorCenterSquare.x = (double)target.x + cellWidth * x;	// cellWidth * x + cellWidth / 2 - cellWidth / 10

src/EulerFluid/FluidField.hpp

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <vector>
+#include "RetentiveArray.hpp"
 
 class VectorField;
 struct SDL_Renderer;
@@ -39,8 +40,7 @@ private:
 
 	VectorField* vel;
 	VectorField* prevVel;
-	std::vector<double>* density;
-	std::vector<double>* prevDensity;
+	RetentiveArray<double, 1> density;
 
 	int lastMouseX, lastMouseY;
 };

src/EulerFluid/main.cpp

@@ -1,4 +1,5 @@
 #include "EulerFluid.hpp"
+#include "RetentiveArray.hpp"
 
 #include <thread>