diff --git a/src/Application.cpp b/src/Application.cpp index 253e3cc..3b6e516 100644 --- a/src/Application.cpp +++ b/src/Application.cpp @@ -37,7 +37,7 @@ Application::Application(int width, int height, const char* title) SDL_SetRenderDrawBlendMode(renderer, SDL_BLENDMODE_BLEND); - field = new FluidField(100); + field = new FluidField(60); before = std::chrono::steady_clock::now(); } @@ -87,6 +87,7 @@ void Application::Update() double frametime = std::chrono::duration_cast>(std::chrono::steady_clock::now() - before).count(); before = std::chrono::steady_clock::now(); + field->VelocityStep(0.002, frametime); field->DensityStep(0.001, frametime); } diff --git a/src/FluidField.cpp b/src/FluidField.cpp index a2f5ef0..2f49300 100644 --- a/src/FluidField.cpp +++ b/src/FluidField.cpp @@ -24,17 +24,14 @@ FluidField::FluidField(int size) : 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.25 * realY); + // hori[y * this->size + x] = 1.0 * realY; + // vert[y * this->size + x] = 1.0 * (- realX - 0.25 * realY); // hori[y * this->size + x] = 0.2 * realY * realY; // vert[y * this->size + x] = 0.2 * realX * realX; } } - ApplyBoundaryConditions(BoundaryCondition::InvertHorizontal, hori); - ApplyBoundaryConditions(BoundaryCondition::InvertVertical, vert); - 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); } @@ -54,6 +51,12 @@ void FluidField::AddSource(int x, int y, double dens, double dt) PVALUE(density, x, y) = std::max(PVALUE(density, x, y), 0.0); } +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; +} + void FluidField::ApplyBoundaryConditions(BoundaryCondition condition, std::vector& field) { int N = this->size - 2; @@ -126,9 +129,140 @@ void FluidField::Advect(double dt) ApplyBoundaryConditions(BoundaryCondition::Continuous, *density); } +void FluidField::DiffuseVelocity(double visc, double dt) +{ + int N = this->size - 2; + double a = dt * visc * N * N; + + for (int k = 0; k < 20; k++) + { + for (int i = 1; i <= N; i++) + { + 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); + } + } + + ApplyBoundaryConditions(BoundaryCondition::Continuous, vel->horizontal); + ApplyBoundaryConditions(BoundaryCondition::Continuous, vel->vertical); + } +} + +void FluidField::AdvectVelocity(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(prevVel->horizontal, i, j); + double y = j - dt0 * VALUE(prevVel->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; + + 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)); + + 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)); + } + } + + ApplyBoundaryConditions(BoundaryCondition::InvertHorizontal, vel->horizontal); + ApplyBoundaryConditions(BoundaryCondition::InvertVertical, vel->vertical); +} + +void FluidField::VelocityStep(double visc, double dt) +{ + // AddFlow(15, 30, 3000.0, 0.0, dt); + // AddFlow(45, 30, -3000.0, 0.0, dt); + // AddFlow(30, 15, 0.0, 3000.0, dt); + + int x, y; + Uint32 buttons = 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 (buttons & SDL_BUTTON_RMASK) + { + AddFlow(lastMouseX, lastMouseY, (dx - lastMouseX) * 500.0, (dy - lastMouseY) * 500.0, dt); + } + + lastMouseX = dx; + lastMouseY = dy; + + std::swap(vel, prevVel); + DiffuseVelocity(visc, dt); + Project(); + std::swap(vel, prevVel); + AdvectVelocity(dt); + Project(); + + vel->RecalculateMagnitude(); +} + +void FluidField::Project() +{ + int N = this->size - 2; + double h = 1.0 / (double)N; + + for (int i = 1; i <= N; i++) + { + 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; + } + } + + ApplyBoundaryConditions(BoundaryCondition::Continuous, prevVel->horizontal); + ApplyBoundaryConditions(BoundaryCondition::Continuous, prevVel->vertical); + + for (int k = 0; k < 20; k++) + { + for (int i = 1; i <= N; i++) + { + 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; + } + } + + ApplyBoundaryConditions(BoundaryCondition::Continuous, prevVel->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; + } + } + + ApplyBoundaryConditions(BoundaryCondition::InvertHorizontal, prevVel->horizontal); + ApplyBoundaryConditions(BoundaryCondition::InvertVertical, prevVel->vertical); +} + void FluidField::DensityStep(double diff, double dt) { - AddSource(this->size / 2, this->size / 2, -10000000.0, dt); // AddSource(50, 3, 60.0, dt); // AddSource(3, 50, 60.0, dt); @@ -137,15 +271,10 @@ void FluidField::DensityStep(double diff, double dt) int dx = (double)(this->size - 2) / (double)(990 - 10) * (double)(x - 10); int dy = (double)(this->size - 2) / (double)(990 - 10) * (double)(y - 10); - int factor = 0; if (buttons & SDL_BUTTON_LMASK) - factor = 1; - else if (buttons & SDL_BUTTON_RMASK) - factor = -1; - - if(dx > 0 && dx < this->size - 1 && dy > 0 && dy < this->size - 1) - AddSource(dx, dy, 300.0 * factor, 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); diff --git a/src/FluidField.hpp b/src/FluidField.hpp index c8906cb..a33ac70 100644 --- a/src/FluidField.hpp +++ b/src/FluidField.hpp @@ -20,12 +20,18 @@ public: ~FluidField(); void AddSource(int x, int y, double density, double dt); + void AddFlow(int x, int y, double dx, double dy, 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 DiffuseVelocity(double visc, double dt); + void AdvectVelocity(double dt); + void VelocityStep(double visc, double dt); + void Project(); + void Draw(SDL_Renderer* renderer, const SDL_Rect& target); private: @@ -35,4 +41,6 @@ private: VectorField* prevVel; std::vector* density; std::vector* prevDensity; + + int lastMouseX, lastMouseY; }; \ No newline at end of file diff --git a/src/VectorField.cpp b/src/VectorField.cpp index ef702f2..7cadd41 100644 --- a/src/VectorField.cpp +++ b/src/VectorField.cpp @@ -24,22 +24,7 @@ VectorField::VectorField(int width, int height, const std::vector& hori, horizontal = hori; vertical = vert; - for (int y = 0; y < this->height; y++) - { - for (int x = 0; x < this->width; x++) - { - double u = horizontal[y * this->width + x]; - double v = vertical[y * this->width + x]; - double magnitude = u + v; - - biggestMagnitude += (biggestMagnitude < magnitude) * (magnitude - biggestMagnitude); - } - } - - if (biggestMagnitude == 0.0) // should use an epsilon probably - biggestMagnitude = 1.0; - - biggestMagnitude = sqrt(biggestMagnitude * 5.0); + RecalculateMagnitude(); } void VectorField::Draw(SDL_Renderer* renderer, const SDL_Rect& targetRect) @@ -68,3 +53,23 @@ void VectorField::Draw(SDL_Renderer* renderer, const SDL_Rect& targetRect) } } } + +void VectorField::RecalculateMagnitude() +{ + for (int y = 0; y < this->height; y++) + { + for (int x = 0; x < this->width; x++) + { + double u = horizontal[y * this->width + x]; + double v = vertical[y * this->width + x]; + double magnitude = u + v; + + biggestMagnitude += (biggestMagnitude < magnitude) * (magnitude - biggestMagnitude); + } + } + + if (biggestMagnitude == 0.0) // should use an epsilon probably + biggestMagnitude = 1.0; + + biggestMagnitude = sqrt(biggestMagnitude * 5.0); +} diff --git a/src/VectorField.hpp b/src/VectorField.hpp index 7f982c3..cbf16fe 100644 --- a/src/VectorField.hpp +++ b/src/VectorField.hpp @@ -13,6 +13,7 @@ public: VectorField(int width, int height, const std::vector& hori, const std::vector& vert); void Draw(SDL_Renderer* renderer, const SDL_Rect& targetRect); + void RecalculateMagnitude(); public: std::vector horizontal;