From 4d59b8a53d1dcc6338370b50ab6db26ea3604812 Mon Sep 17 00:00:00 2001 From: Lauchmelder Date: Thu, 16 Dec 2021 15:07:11 +0100 Subject: [PATCH] created retentive array datastructure --- lib/CMakeLists.txt | 2 +- lib/RetentiveArray.hpp | 180 ++++++++++++++++++++++++++++++++++ src/EulerFluid/FluidField.cpp | 30 +++--- src/EulerFluid/FluidField.hpp | 4 +- src/EulerFluid/main.cpp | 1 + 5 files changed, 197 insertions(+), 20 deletions(-) create mode 100644 lib/RetentiveArray.hpp diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index 6f1ed48..d4c3b4a 100644 --- a/lib/CMakeLists.txt +++ b/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}) diff --git a/lib/RetentiveArray.hpp b/lib/RetentiveArray.hpp new file mode 100644 index 0000000..bba74dd --- /dev/null +++ b/lib/RetentiveArray.hpp @@ -0,0 +1,180 @@ +#pragma once + +#include +#include +#include + +/** +* @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(size, Type()); + } + } + + /** + * @brief Constructs a new retentive array with the data of another array + * + * @param other The retentive array to copy from + */ + RetentiveArray& operator=(const RetentiveArray& other) + { + for (int n = 0; n <= AttentionSpan; n++) + { + data[n] = std::vector(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& operator=(RetentiveArray&& other) + { + for (int n = 0; n <= AttentionSpan; n++) + { + data[n] = std::vector(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& 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&& 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& operator[](size_t index) + { + return data[index]; + } + + /** + * @brief Get the most up-to-date array + * + * @return The array with generation 0 + */ + std::vector& 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 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 rule) + { + this->rule = rule; + } + + +private: + std::array, AttentionSpan + 1> data; + std::function rule = std::bind([]() {}); // Do nothing by default +}; diff --git a/src/EulerFluid/FluidField.cpp b/src/EulerFluid/FluidField.cpp index d9bd10b..c24e368 100644 --- a/src/EulerFluid/FluidField.cpp +++ b/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(this->size * this->size, 0.0); - prevDensity = new std::vector(this->size * this->size, 0.0); + density = RetentiveArray(this->size * this->size); std::vector hori(this->size * this->size); std::vector 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 diff --git a/src/EulerFluid/FluidField.hpp b/src/EulerFluid/FluidField.hpp index a33ac70..0e40600 100644 --- a/src/EulerFluid/FluidField.hpp +++ b/src/EulerFluid/FluidField.hpp @@ -1,6 +1,7 @@ #pragma once #include +#include "RetentiveArray.hpp" class VectorField; struct SDL_Renderer; @@ -39,8 +40,7 @@ private: VectorField* vel; VectorField* prevVel; - std::vector* density; - std::vector* prevDensity; + RetentiveArray density; int lastMouseX, lastMouseY; }; \ No newline at end of file diff --git a/src/EulerFluid/main.cpp b/src/EulerFluid/main.cpp index c905b3f..6325798 100644 --- a/src/EulerFluid/main.cpp +++ b/src/EulerFluid/main.cpp @@ -1,4 +1,5 @@ #include "EulerFluid.hpp" +#include "RetentiveArray.hpp" #include