diff --git a/src/Application.cpp b/src/Application.cpp index fb404b1..a7d4133 100644 --- a/src/Application.cpp +++ b/src/Application.cpp @@ -35,41 +35,16 @@ Application::Application(int width, int height, const char* title) return; } - // Construct dummy velocity field - // TODO: Remove eventually - int fieldSize = 39; - double curl1 = -1.0; - double curl2 = 1.0; - std::vector hori(fieldSize * fieldSize); - std::vector vert(fieldSize * fieldSize); + SDL_SetRenderDrawBlendMode(renderer, SDL_BLENDMODE_BLEND); - for (int y = 0; y < fieldSize; y++) - { - for (int x = 0; x < fieldSize; x++) - { - // map internal coords to "real world" coordinates - double realX = -2.0 + (4.0 / (double)fieldSize) * (double)x; - double realY = -2.0 + (4.0 / (double)fieldSize) * (double)y; - - // hori[y * fieldSize + x] = -(realX + 1.0) / sqrt((realX + 1.0) * (realX + 1.0) + (realY + 1.0) * (realY + 1.0)) + (realX - 1.0) / sqrt((realX - 1.0) * (realX - 1.0) + (realY - 1.0) * (realY - 1.0)); - // vert[y * fieldSize + x] = (realY + 1.0) / sqrt((realX + 1.0) * (realX + 1.0) + (realY + 1.0) * (realY + 1.0)) - (realY - 1.0) / sqrt((realX - 1.0) * (realX - 1.0) + (realY - 1.0) * (realY - 1.0)); - - // hori[y * fieldSize + x] = realY / sqrt(realX * realX + realY * realY); - // vert[y * fieldSize + x] = -realX / sqrt(realX * realX + realY * realY); - - // hori[y * fieldSize + x] = -realY; - // vert[y * fieldSize + x] = realX; - - hori[y * fieldSize + x] = realY; - vert[y * fieldSize + x] = -realX - 1.0*realY; - } - } - - velocity = VectorField(fieldSize, fieldSize, hori, vert); + field = new FluidField(60); + before = std::chrono::steady_clock::now(); } Application::~Application() { + delete field; + SDL_DestroyRenderer(renderer); // Let's just destroy this renderer regardless of other applications in this program :) what could go wrong SDL_DestroyWindow(window); } @@ -109,7 +84,10 @@ void Application::HandleEvents() void Application::Update() { + double frametime = std::chrono::duration_cast>(std::chrono::steady_clock::now() - before).count(); + before = std::chrono::steady_clock::now(); + field->DensityStep(0.0005, frametime); } void Application::Render() @@ -118,7 +96,7 @@ void Application::Render() SDL_RenderClear(renderer); - velocity.Draw(renderer, {10, 10, 790, 790}); + field->Draw(renderer, {10, 10, 990, 990}); SDL_RenderPresent(renderer); } diff --git a/src/Application.hpp b/src/Application.hpp index 4ffbe8b..7e84ab7 100644 --- a/src/Application.hpp +++ b/src/Application.hpp @@ -1,6 +1,7 @@ #pragma once -#include "VectorField.hpp" +#include +#include "FluidField.hpp" // Forward Declarations struct SDL_Window; @@ -24,7 +25,8 @@ private: SDL_Window* window = nullptr; SDL_Renderer* renderer = nullptr; + std::chrono::steady_clock::time_point before; bool shouldClose = false; - VectorField velocity; + FluidField* field; }; \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ac395a3..cc1e8ed 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -4,7 +4,7 @@ cmake_minimum_required (VERSION 3.8) # Add source to this project's executable. -add_executable (EulerFluid "main.cpp" "Application.hpp" "Application.cpp" "VectorField.hpp" "VectorField.cpp") +add_executable (EulerFluid "main.cpp" "Application.hpp" "Application.cpp" "VectorField.hpp" "VectorField.cpp" "FluidField.hpp" "FluidField.cpp") target_include_directories(EulerFluid PUBLIC ${SDL2_INCLUDE_DIRS}) target_link_libraries(EulerFluid PUBLIC ${SDL2_LIBRARIES}) diff --git a/src/FluidField.cpp b/src/FluidField.cpp new file mode 100644 index 0000000..7934ab2 --- /dev/null +++ b/src/FluidField.cpp @@ -0,0 +1,166 @@ +#include "FluidField.hpp" + +#include +#include +#include "VectorField.hpp" + +#define VALUE(arr, x, y) ((arr)[(y) * this->size + (x)]) +#define PVALUE(arr, x, y) ((*(arr))[(y) * this->size + (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); + + std::vector hori(this->size * this->size); + std::vector vert(this->size * this->size); + + for (int y = 1; y < this->size - 1; y++) + { + for (int x = 1; x < this->size - 1; x++) + { + // map internal coords to "real world" coordinates + double realX = -2.0 + (4.0 / (double)size) * (double)x; + double realY = -2.0 + (4.0 / (double)size) * (double)y; + + // hori[y * this->size + x] = 1.0 * realY; + // vert[y * this->size + x] = 1.0 * (- realX - 0.0 * realY); + + hori[y * this->size + x] = 0.2 * realY * realY; + vert[y * this->size + x] = 0.2 * realX * realX; + } + } + + vel = new VectorField(this->size, this->size, hori, vert); + prevVel = new VectorField(this->size, this->size, hori, vert); +} + +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; +} + +void FluidField::ApplyBoundaryConditions(BoundaryCondition condition, std::vector& field) +{ + int N = this->size - 2; + for (int i = 1; i <= N; i++) + { + field[i * (N + 2) + 0] = (condition == BoundaryCondition::InvertHorizontal) ? -field[i * (N + 2) + 1] : field[i * (N + 2) + 1]; + field[i * (N + 2) + N + 1] = (condition == BoundaryCondition::InvertHorizontal) ? -field[i * (N + 2) + N] : field[i * (N + 2) + N]; + field[0 * (N + 2) + i] = (condition == BoundaryCondition::InvertVertical) ? -field[1 * (N + 2) + i] : field[1 * (N + 2) + i]; + field[(N + 1) * (N + 2) + i] = (condition == BoundaryCondition::InvertVertical) ? -field[N * (N + 2) + i] : field[N * (N + 2) + i]; + } + + field[0 * (N + 2) + 0] = 0.5 * (field[0 * (N + 2) + 1] + field[1 * (N + 2) + 0]); + field[(N + 1) * (N + 2) + 0] = 0.5 * (field[(N + 1) * (N + 2) + 1] + field[N * (N + 2) + 0]); + field[0 * (N + 2) + (N + 1)] = 0.5 * (field[0 * (N + 2) + N] + field[1 * (N + 2) + (N + 1)]); + field[(N + 1) * (N + 2) + (N + 1)] = 0.5 * (field[(N + 1) * (N + 2) + N] + field[N * (N + 2) + (N + 1)]); +} + +void FluidField::Diffuse(double diff, double dt) +{ + int N = this->size - 2; + double a = dt * diff * N * N; + + for (int k = 0; k < 20; k++) + { + for (int i = 1; i <= N; i++) + { + 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); + } + } + + ApplyBoundaryConditions(BoundaryCondition::Continuous, *density); + } +} + +void FluidField::Advect(double dt) +{ + int N = this->size - 2; + double dt0 = dt * N; + + for (int i = 1; i <= N; i++) + { + 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); + + if (x < 0.5) x = 0.5; + if (x > N + 0.5) x = N + 0.5; + if (y < 0.5) y = 0.5; + if (y > N + 0.5) y = N + 0.5; + + int i0 = (int)x; + int i1 = i0 + 1; + int j0 = (int)y; + int j1 = j0 + 1; + + double s1 = x - i0; + double s0 = 1 - s1; + 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)); + } + } + + ApplyBoundaryConditions(BoundaryCondition::Continuous, *density); +} + +void FluidField::DensityStep(double diff, double dt) +{ + // AddSource(3, 3, 60.0, dt); + // AddSource(50, 3, 60.0, dt); + // AddSource(3, 50, 60.0, dt); + + int x, y; + SDL_GetMouseState(&x, &y); + int dx = (double)(this->size - 2) / (double)(990 - 10) * (double)(x - 10); + int dy = (double)(this->size - 2) / (double)(990 - 10) * (double)(y - 10); + if(dx > 1 && dx < this->size - 2 && dy > 1 && dy < this->size - 2) + AddSource(dx, dy, 60.0, dt); + + std::swap(prevDensity, density); + Diffuse(diff, dt); + std::swap(prevDensity, density); + Advect(dt); +} + +void FluidField::Draw(SDL_Renderer* renderer, const SDL_Rect& target) +{ + double cellWidth = (double)(target.w - target.x) / (double)this->size; + double cellHeight = (double)(target.h - target.y) / (double)this->size; + + SDL_FRect vectorCenterSquare; + vectorCenterSquare.w = cellWidth; + vectorCenterSquare.h = cellHeight; + + for (int y = 0; y < this->size; y++) + { + for (int x = 0; x < this->size; x++) + { + double densityVal = std::min(PVALUE(density, x, y), 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 + vectorCenterSquare.y = (double)target.y + cellHeight * y; + SDL_RenderFillRectF(renderer, &vectorCenterSquare); + } + } + + vel->Draw(renderer, target); +} diff --git a/src/FluidField.hpp b/src/FluidField.hpp new file mode 100644 index 0000000..c8906cb --- /dev/null +++ b/src/FluidField.hpp @@ -0,0 +1,38 @@ +#pragma once + +#include + +class VectorField; +struct SDL_Renderer; +struct SDL_Rect; + +enum class BoundaryCondition +{ + Continuous, + InvertVertical, + InvertHorizontal +}; + +class FluidField +{ +public: + FluidField(int size); + ~FluidField(); + + void AddSource(int x, int y, double density, double dt); + void ApplyBoundaryConditions(BoundaryCondition condition, std::vector& field); + + void Diffuse(double diff, double dt); + void Advect(double dt); + void DensityStep(double diff, double dt); + + void Draw(SDL_Renderer* renderer, const SDL_Rect& target); + +private: + int size; + + VectorField* vel; + VectorField* prevVel; + std::vector* density; + std::vector* prevDensity; +}; \ No newline at end of file diff --git a/src/VectorField.cpp b/src/VectorField.cpp index 350029e..32693db 100644 --- a/src/VectorField.cpp +++ b/src/VectorField.cpp @@ -45,7 +45,7 @@ void VectorField::Draw(SDL_Renderer* renderer, const SDL_Rect& targetRect) vectorCenterSquare.w = cellWidth / 5.0; vectorCenterSquare.h = cellHeight / 5.0; - SDL_SetRenderDrawColor(renderer, 200, 20, 20, 255); + SDL_SetRenderDrawColor(renderer, 200, 20, 20, 60); for (int y = 0; y < height; y++) { for (int x = 0; x < width; x++) diff --git a/src/VectorField.hpp b/src/VectorField.hpp index e3b1b5d..7f982c3 100644 --- a/src/VectorField.hpp +++ b/src/VectorField.hpp @@ -14,11 +14,12 @@ public: void Draw(SDL_Renderer* renderer, const SDL_Rect& targetRect); -private: - int width, height; - +public: std::vector horizontal; std::vector vertical; +private: + int width, height; + double biggestMagnitude = 0.0; }; \ No newline at end of file diff --git a/src/main.cpp b/src/main.cpp index d533295..08d93ea 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -2,7 +2,7 @@ int main(int argc, char** argv) { - Application* app = new Application(800, 800, "Euler Fluid Simulation"); + Application* app = new Application(1000, 1000, "Euler Fluid Simulation"); app->Launch(); delete app;