lauchmelder/NumericalMethods
1
0
Fork 0
Code Issues Pull requests Actions Packages Projects Releases Wiki Activity

added more retentive stuff

Browse source
This commit is contained in:
Lauchmelder 2021-12-16 16:35:37 +01:00
parent 75df9ad96d
commit b947b0b263
8 changed files with 117 additions and 107 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
lib/CMakeLists.txt
View file

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

99
lib/RetentiveEntity.hpp
View file

@ -3,6 +3,76 @@
#include <vector>
#include <array>
#include <functional>
#include <memory>
template<
typename Type,
unsigned int AttentionSpan,
bool IsPointerType>
class _RetentiveEntityBase;
template<
typename Type,
unsigned int AttentionSpan>
class _RetentiveEntityBase<Type, AttentionSpan, false>
{
public:
/**
* @brief Get the entity from `index` generations ago
*
* @param index Amount of generations to go backwards in time
* @return The entity from before `index` generations
*/
Type& operator[](size_t index)
{
return data[index];
}
/**
* @brief Get the most up-to-date entity
*
* @return The entity with generation 0
*/
Type& Current()
{
return data[0];
}
protected:
std::array<Type, AttentionSpan + 1> data;
};
template<
typename Type,
unsigned int AttentionSpan>
class _RetentiveEntityBase<Type, AttentionSpan, true>
{
public:
/**
* @brief Get the entity from `index` generations ago
*
* @param index Amount of generations to go backwards in time
* @return The entity from before `index` generations
*/
Type& operator[](size_t index)
{
return *(data[index]);
}
/**
* @brief Get the most up-to-date entity
*
* @return The entity with generation 0
*/
Type& Current()
{
return *(data[0]);
}
protected:
std::array<std::shared_ptr<Type>, AttentionSpan + 1> data; // Shared for move semantics
};
/**
* @brief Abstract base type for all retentive things
@ -13,8 +83,8 @@
template<
typename Type,
unsigned int AttentionSpan,
typename std::enable_if_t<(AttentionSpan > 0), bool> = true>
class RetentiveEntity
bool IsPointerType = false>
class RetentiveEntity : public _RetentiveEntityBase<Type, AttentionSpan, IsPointerType>
{
public:
RetentiveEntity() {}
@ -23,27 +93,6 @@ public:
RetentiveEntity<Type, AttentionSpan>& operator=(const RetentiveEntity<Type, AttentionSpan>& other) = delete;
RetentiveEntity<Type, AttentionSpan>& operator=(RetentiveEntity<Type, AttentionSpan>&& other) = delete;
/**
* @brief Get the entity from `index` generations ago
*
* @param index Amount of generations to go backwards in time
* @return The entity from before `index` generations
*/
virtual Type& operator[](size_t index)
{
return data[index];
}
/**
* @brief Get the most up-to-date entity
*
* @return The entity with generation 0
*/
virtual Type& Current()
{
return data[0];
}
/**
* @brief Evolve the data in the entity
*
@ -92,8 +141,4 @@ protected:
* @brief Swaps the objects in the array
*/
virtual void CycleGenerations() = 0;
protected:
std::array<Type, AttentionSpan + 1> data;
std::function<void(void)> rule = std::bind([]() {}); // Do nothing by default
};

42
lib/RetentiveObject.hpp
View file

@ -21,9 +21,8 @@
*/
template<
typename Type,
unsigned int AttentionSpan,
typename std::enable_if_t<(AttentionSpan > 0), bool> = true>
class RetentiveObject : public RetentiveEntity<Type*, AttentionSpan>
unsigned int AttentionSpan>
class RetentiveObject : public RetentiveEntity<Type, AttentionSpan, true>
{
public:
/**
@ -34,7 +33,7 @@ public:
// Create new vectors for every generation that needs to be remembered
for (int n = 0; n <= AttentionSpan; n++)
{
data[n] = new Type();
data[n] = std::make_shared<Type>();
}
}
@ -45,10 +44,10 @@ public:
*/
RetentiveObject(const Type& initVal)
{
// Create new vectors for every generation that needs to be remembered
// Create new objects for every generation that needs to be remembered
for (int n = 0; n <= AttentionSpan; n++)
{
data[n] = new Type(initVal);
data[n] = std::make_shared<Type>(initVal);
}
}
@ -77,7 +76,6 @@ public:
for (int n = 0; n <= AttentionSpan; n++)
{
data[n] = other.data[n];
other.data[n] = nullptr;
}
return *this;
@ -103,36 +101,6 @@ public:
*this = other;
}
~RetentiveObject()
{
for (int n = 0; n <= AttentionSpan; n++)
{
if (data[n] != nullptr)
delete data[n];
}
}
/**
* @brief Get the object from `index` generations ago
*
* @param index Amount of generations to go backwards in time
* @return The object from before `index` generations
*/
Type& operator[](size_t index) override
{
return *(data[index]);
}
/**
* @brief Get the most up-to-date object
*
* @return The object with generation 0
*/
Type& Current() override
{
return *(data[0]);
}
private:
/**
* @brief Swaps the objects in the array

2
src/EulerFluid/VectorField.cpp → lib/VectorField.cpp
View file

@ -62,7 +62,7 @@ void VectorField::RecalculateMagnitude()
{
double u = horizontal[y * this->width + x];
double v = vertical[y * this->width + x];
double magnitude = u*u + v*v;
double magnitude = u * u + v * v;
biggestMagnitude = std::max(biggestMagnitude, magnitude);
}

0
src/EulerFluid/VectorField.hpp → lib/VectorField.hpp
View file

2
src/EulerFluid/CMakeLists.txt
View file

@ -4,7 +4,7 @@
cmake_minimum_required (VERSION 3.8)
# Add source to this project's executable.
add_executable (EulerFluid "main.cpp" "EulerFluid.hpp" "EulerFluid.cpp" "VectorField.hpp" "VectorField.cpp" "FluidField.hpp" "FluidField.cpp")
add_executable (EulerFluid "main.cpp" "EulerFluid.hpp" "EulerFluid.cpp" "FluidField.hpp" "FluidField.cpp")
target_include_directories(EulerFluid PUBLIC nm_utils)
target_link_libraries(EulerFluid PRIVATE nm_utils)

67
src/EulerFluid/FluidField.cpp
View file

@ -2,6 +2,7 @@
#include <iostream>
#include <SDL.h>
#include "VectorField.hpp"
#define VALUE(arr, x, y) ((arr)[(y) * this->size + (x)])
@ -33,14 +34,12 @@ FluidField::FluidField(int size) :
}
}
vel = new VectorField(this->size, this->size, hori, vert); // gonna do the same inefficient calculations twice, why not
prevVel = new VectorField(this->size, this->size, hori, vert);
velocity = RetentiveObject<VectorField, 1>(VectorField(this->size, this->size, hori, vert));
}
FluidField::~FluidField()
{
delete prevVel;
delete vel;
// Do nothing
}
void FluidField::AddSource(int x, int y, double dens, double dt)
@ -51,8 +50,8 @@ void FluidField::AddSource(int x, int y, double dens, double dt)
void FluidField::AddFlow(int x, int y, double dx, double dy, double dt)
{
VALUE(vel->horizontal, x, y) += dt * dx;
VALUE(vel->vertical, x, y) += dt * dy;
velocity.Current().horizontal[IDX(x, y, size)] += dt * dx;
velocity.Current().vertical[IDX(x, y, size)] += dt * dy;
}
void FluidField::ApplyBoundaryConditions(BoundaryCondition condition, std::vector<double>& field)
@ -101,8 +100,8 @@ void FluidField::Advect(double dt)
{
for (int j = 1; j <= N; j++)
{
double x = i - dt0 * VALUE(vel->horizontal, i, j);
double y = j - dt0 * VALUE(vel->vertical, i, j);
double x = i - dt0 * velocity.Current().horizontal[IDX(i, j, size)];
double y = j - dt0 * velocity.Current().vertical[IDX(i, j, size)];
if (x < 0.5) x = 0.5;
if (x > N + 0.5) x = N + 0.5;
@ -138,13 +137,13 @@ void FluidField::DiffuseVelocity(double visc, double dt)
{
for (int j = 1; j <= N; j++)
{
VALUE(vel->horizontal, i, j) = (VALUE(prevVel->horizontal, i, j) + a * (VALUE(vel->horizontal, i - 1, j) + VALUE(vel->horizontal, i + 1, j) + VALUE(vel->horizontal, i, j - 1) + VALUE(vel->horizontal, i, j + 1))) / (1 + 4 * a);
VALUE(vel->vertical, i, j) = (VALUE(prevVel->vertical, i, j) + a * (VALUE(vel->vertical, i - 1, j) + VALUE(vel->vertical, i + 1, j) + VALUE(vel->vertical, i, j - 1) + VALUE(vel->vertical, i, j + 1))) / (1 + 4 * a);
velocity.Current().horizontal[IDX(i, j, size)] = (velocity[1].horizontal[IDX(i, j, size)] + a * (velocity[1].horizontal[IDX(i - 1, j, size)] + velocity[1].horizontal[IDX(i + 1, j, size)] + velocity[1].horizontal[IDX(i, j - 1, size)] + velocity[1].horizontal[IDX(i, j + 1, size)])) / (1 + 4 * a);
velocity.Current().vertical[IDX(i, j, size)] = (velocity[1].vertical[IDX(i, j, size)] + a * (velocity[1].vertical[IDX(i - 1, j, size)] + velocity[1].vertical[IDX(i + 1, j, size)] + velocity[1].vertical[IDX(i, j - 1, size)] + velocity[1].vertical[IDX(i, j + 1, size)])) / (1 + 4 * a);
}
}
ApplyBoundaryConditions(BoundaryCondition::Continuous, vel->horizontal);
ApplyBoundaryConditions(BoundaryCondition::Continuous, vel->vertical);
ApplyBoundaryConditions(BoundaryCondition::Continuous, velocity.Current().horizontal);
ApplyBoundaryConditions(BoundaryCondition::Continuous, velocity.Current().vertical);
}
}
@ -157,8 +156,8 @@ void FluidField::AdvectVelocity(double dt)
{
for (int j = 1; j <= N; j++)
{
double x = i - dt0 * VALUE(prevVel->horizontal, i, j);
double y = j - dt0 * VALUE(prevVel->vertical, i, j);
double x = i - dt0 * velocity[1].horizontal[IDX(i, j, size)];
double y = j - dt0 * velocity[1].vertical[IDX(i, j, size)];
if (x < 0.5) x = 0.5;
if (x > N + 0.5) x = N + 0.5;
@ -175,16 +174,16 @@ void FluidField::AdvectVelocity(double dt)
double t1 = y - j0;
double t0 = 1 - t1;
VALUE(vel->horizontal, i, j) = s0 * (t0 * VALUE(prevVel->horizontal, i0, j0) + t1 * VALUE(prevVel->horizontal, i0, j1)) +
s1 * (t0 * VALUE(prevVel->horizontal, i1, j0) + t1 * VALUE(prevVel->horizontal, i1, j1));
velocity.Current().horizontal[IDX(i, j, size)] = s0 * (t0 * velocity[1].horizontal[IDX(i0, j0, size)] + t1 * velocity[1].horizontal[IDX(i0, j1, size)]) +
s1 * (t0 * velocity[1].horizontal[IDX(i1, j0, size)] + t1 * velocity[1].horizontal[IDX(i1, j1, size)]);
VALUE(vel->vertical, i, j) = s0 * (t0 * VALUE(prevVel->vertical, i0, j0) + t1 * VALUE(prevVel->vertical, i0, j1)) +
s1 * (t0 * VALUE(prevVel->vertical, i1, j0) + t1 * VALUE(prevVel->vertical, i1, j1));
velocity.Current().vertical[IDX(i, j, size)] = s0 * (t0 * velocity[1].vertical[IDX(i0, j0, size)] + t1 * velocity[1].vertical[IDX(i0, j1, size)]) +
s1 * (t0 * velocity[1].vertical[IDX(i1, j0, size)] + t1 * velocity[1].vertical[IDX(i1, j1, size)]);
}
}
ApplyBoundaryConditions(BoundaryCondition::InvertHorizontal, vel->horizontal);
ApplyBoundaryConditions(BoundaryCondition::InvertVertical, vel->vertical);
ApplyBoundaryConditions(BoundaryCondition::InvertHorizontal, velocity.Current().horizontal);
ApplyBoundaryConditions(BoundaryCondition::InvertVertical, velocity.Current().vertical);
}
void FluidField::VelocityStep(double visc, double dt)
@ -206,11 +205,9 @@ void FluidField::VelocityStep(double visc, double dt)
lastMouseX = dx;
lastMouseY = dy;
std::swap(vel, prevVel);
DiffuseVelocity(visc, dt);
velocity.Evolve(std::bind(&FluidField::DiffuseVelocity, this, visc, dt));
Project();
std::swap(vel, prevVel);
AdvectVelocity(dt);
velocity.Evolve(std::bind(&FluidField::AdvectVelocity, this, dt));
Project();
// vel->RecalculateMagnitude();
@ -225,13 +222,13 @@ void FluidField::Project()
{
for (int j = 1; j <= N; j++)
{
VALUE(prevVel->vertical, i, j) = -0.5 * h * (VALUE(vel->horizontal, i + 1, j) - VALUE(vel->horizontal, i - 1, j) + VALUE(vel->vertical, i, j + 1) - VALUE(vel->vertical, i, j - 1));
VALUE(prevVel->horizontal, i, j) = 0;
velocity[1].vertical[IDX(i, j, size)] = -0.5 * h * (velocity.Current().vertical[IDX(i + 1, j, size)] - velocity.Current().vertical[IDX(i - 1, j, size)] + velocity.Current().vertical[IDX(i, j + 1, size)] - velocity.Current().vertical[IDX(i, j - 1, size)]);
velocity[1].horizontal[IDX(i, j, size)] = 0;
}
}
ApplyBoundaryConditions(BoundaryCondition::Continuous, prevVel->horizontal);
ApplyBoundaryConditions(BoundaryCondition::Continuous, prevVel->vertical);
ApplyBoundaryConditions(BoundaryCondition::Continuous, velocity[1].horizontal);
ApplyBoundaryConditions(BoundaryCondition::Continuous, velocity[1].vertical);
for (int k = 0; k < 20; k++)
{
@ -239,24 +236,24 @@ void FluidField::Project()
{
for (int j = 1; j <= N; j++)
{
VALUE(prevVel->horizontal, i, j) = (VALUE(prevVel->vertical, i, j) + VALUE(prevVel->horizontal, i - 1, j) + VALUE(prevVel->horizontal, i + 1, j) + VALUE(prevVel->vertical, i, j - 1) + VALUE(prevVel->vertical, i, j + 1)) / 4.0;
velocity[1].horizontal[IDX(i, j, size)] = (velocity[1].vertical[IDX(i, j, size)] + velocity[1].vertical[IDX(i - 1, j, size)] + velocity[1].vertical[IDX(i + 1, j, size)] + velocity[1].vertical[IDX(i, j - 1, size)] + velocity[1].vertical[IDX(i, j + 1, size)]) / 4.0;
}
}
ApplyBoundaryConditions(BoundaryCondition::Continuous, prevVel->horizontal);
ApplyBoundaryConditions(BoundaryCondition::Continuous, velocity[1].horizontal);
}
for (int i = 1; i <= N; i++)
{
for (int j = 1; j <= N; j++)
{
VALUE(vel->horizontal, i, j) -= 0.5 * (VALUE(prevVel->horizontal, i + 1, j) - VALUE(prevVel->horizontal, i - 1, j)) / h;
VALUE(vel->vertical, i, j) -= 0.5 * (VALUE(prevVel->horizontal, i, j + 1) - VALUE(prevVel->horizontal, i, j - 1)) / h;
velocity.Current().horizontal[IDX(i, j, size)] -= 0.5 * (velocity[1].horizontal[IDX(i + 1, j, size)] - velocity[1].horizontal[IDX(i - 1, j, size)]) / h;
velocity.Current().vertical[IDX(i, j, size)] -= 0.5 * (velocity[1].horizontal[IDX(i, j + 1, size)] - velocity[1].horizontal[IDX(i, j - 1, size)]) / h;
}
}
ApplyBoundaryConditions(BoundaryCondition::InvertHorizontal, prevVel->horizontal);
ApplyBoundaryConditions(BoundaryCondition::InvertVertical, prevVel->vertical);
ApplyBoundaryConditions(BoundaryCondition::InvertHorizontal, velocity[1].horizontal);
ApplyBoundaryConditions(BoundaryCondition::InvertVertical, velocity[1].vertical);
}
void FluidField::DensityStep(double diff, double dt)
@ -300,5 +297,5 @@ void FluidField::Draw(SDL_Renderer* renderer, const SDL_Rect& target)
}
}
vel->Draw(renderer, target);
velocity.Current().Draw(renderer, target);
}

6
src/EulerFluid/FluidField.hpp
View file

@ -1,9 +1,10 @@
#pragma once
#include <vector>
#include "VectorField.hpp"
#include "RetentiveArray.hpp"
#include "RetentiveObject.hpp"
class VectorField;
struct SDL_Renderer;
struct SDL_Rect;
@ -38,8 +39,7 @@ public:
private:
int size;
VectorField* vel;
VectorField* prevVel;
RetentiveObject<VectorField, 1> velocity;
RetentiveArray<double, 1> density;
int lastMouseX, lastMouseY;