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

added velocity evolution

Browse source
This commit is contained in:
Lauchmelder 2021-12-10 22:18:20 +01:00
parent aa8393ee4c
commit f8e8cc64b4
5 changed files with 174 additions and 30 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
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(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::duration<double>>(std::chrono::steady_clock::now() - before).count();
before = std::chrono::steady_clock::now();
field->VelocityStep(0.002, frametime);
field->DensityStep(0.001, frametime);
}

155
src/FluidField.cpp
View file

@ -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<double>& 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);

8
src/FluidField.hpp
View file

@ -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<double>& 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<double>* density;
std::vector<double>* prevDensity;
int lastMouseX, lastMouseY;
};

37
src/VectorField.cpp
View file

@ -24,22 +24,7 @@ VectorField::VectorField(int width, int height, const std::vector<double>& 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);
}

1
src/VectorField.hpp
View file

@ -13,6 +13,7 @@ public:
VectorField(int width, int height, const std::vector<double>& hori, const std::vector<double>& vert);
void Draw(SDL_Renderer* renderer, const SDL_Rect& targetRect);
void RecalculateMagnitude();
public:
std::vector<double> horizontal;