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

improved magnitude finding algorithm

Browse source
This commit is contained in:
Lauchmelder 2021-12-10 18:55:04 +01:00
parent bd061968dc
commit aa8393ee4c
3 changed files with 33 additions and 23 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

6
src/Application.cpp
View file

@ -37,7 +37,7 @@ Application::Application(int width, int height, const char* title)
SDL_SetRenderDrawBlendMode(renderer, SDL_BLENDMODE_BLEND);
field = new FluidField(60);
field = new FluidField(100);
before = std::chrono::steady_clock::now();
}
@ -87,7 +87,7 @@ void Application::Update()
double frametime = std::chrono::duration_cast<std::chrono::duration<double>>(std::chrono::steady_clock::now() - before).count();
before = std::chrono::steady_clock::now();
field->DensityStep(0.0005, frametime);
field->DensityStep(0.001, frametime);
}
void Application::Render()
@ -96,7 +96,7 @@ void Application::Render()
SDL_RenderClear(renderer);
field->Draw(renderer, {10, 10, 990, 990});
field->Draw(renderer, {0, 0, 1000, 1000});
SDL_RenderPresent(renderer);
}

34
src/FluidField.cpp
View file

@ -24,15 +24,18 @@ 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.0 * 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;
// 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);
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);
}
@ -56,16 +59,17 @@ void FluidField::ApplyBoundaryConditions(BoundaryCondition condition, std::vecto
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];
VALUE(field, 0 , i) = (condition == BoundaryCondition::InvertHorizontal) ? -VALUE(field, 1, i) : VALUE(field, 1, i); // VALUE(field, N, i);
VALUE(field, N + 1 , i) = (condition == BoundaryCondition::InvertHorizontal) ? -VALUE(field, N, i) : VALUE(field, N, i); // VALUE(field, 1, i);
VALUE(field, i , 0) = (condition == BoundaryCondition::InvertVertical) ? -VALUE(field, i, 1) : VALUE(field, i, 1); // VALUE(field, i, N);
VALUE(field, i , N + 1) = (condition == BoundaryCondition::InvertVertical) ? -VALUE(field, i, N) : VALUE(field, i, N); // VALUE(field, i, 1);
}
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)]);
VALUE(field, 0 , 0 ) = 0.5 * (VALUE(field, 1, 0 ) + VALUE(field, 0, 1 ));
VALUE(field, 0 , N + 1 ) = 0.5 * (VALUE(field, 1, N + 1) + VALUE(field, 0, N ));
VALUE(field, N + 1 , 0 ) = 0.5 * (VALUE(field, N, 0 ) + VALUE(field, N + 1, 1));
VALUE(field, N + 1 , N + 1 ) = 0.5 * (VALUE(field, N, N + 1) + VALUE(field, N + 1, N));
}
void FluidField::Diffuse(double diff, double dt)
@ -124,7 +128,7 @@ void FluidField::Advect(double dt)
void FluidField::DensityStep(double diff, double dt)
{
// AddSource(3, 3, 60.0, dt);
AddSource(this->size / 2, this->size / 2, -10000000.0, dt);
// AddSource(50, 3, 60.0, dt);
// AddSource(3, 50, 60.0, dt);
@ -141,7 +145,7 @@ void FluidField::DensityStep(double diff, double dt)
factor = -1;
if(dx > 0 && dx < this->size - 1 && dy > 0 && dy < this->size - 1)
AddSource(dx, dy, 60.0 * factor, dt);
AddSource(dx, dy, 300.0 * factor, dt);
std::swap(prevDensity, density);
Diffuse(diff, dt);

16
src/VectorField.cpp
View file

@ -1,5 +1,7 @@
#include "VectorField.hpp"
#include <chrono>
#include <iostream>
#include <SDL.h>
VectorField::VectorField() :
@ -22,18 +24,22 @@ VectorField::VectorField(int width, int height, const std::vector<double>& hori,
horizontal = hori;
vertical = vert;
for (double u : horizontal)
for (int y = 0; y < this->height; y++)
{
for (double v : vertical)
for (int x = 0; x < this->width; x++)
{
biggestMagnitude = std::max(biggestMagnitude, u * u + v * v);
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);
biggestMagnitude = sqrt(biggestMagnitude * 5.0);
}
void VectorField::Draw(SDL_Renderer* renderer, const SDL_Rect& targetRect)
@ -45,7 +51,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, 60);
SDL_SetRenderDrawColor(renderer, 200, 20, 20, 100);
for (int y = 0; y < height; y++)
{
for (int x = 0; x < width; x++)