From d6d4a30798cafa2fd303a9abe044846a270a025a Mon Sep 17 00:00:00 2001 From: Yggdrasil75 Date: Fri, 6 Feb 2026 14:31:04 -0500 Subject: [PATCH] better simulation --- tests/fluidsim.cpp | 40 ++++++++++------- util/output/frame.hpp | 102 ++++++++++++++++++++++++++++++++++++++++++ util/sim/fluidsim.hpp | 62 ++++++++++--------------- 3 files changed, 148 insertions(+), 56 deletions(-) diff --git a/tests/fluidsim.cpp b/tests/fluidsim.cpp index c0d7b60..243e6af 100644 --- a/tests/fluidsim.cpp +++ b/tests/fluidsim.cpp @@ -6,8 +6,8 @@ int main() { // Simulation settings const int TOTAL_FRAMES = 10000; - const int WIDTH = 800; - const int HEIGHT = 600; + const int WIDTH = 1024; + const int HEIGHT = 1024; // Setup Camera Camera cam; @@ -26,8 +26,10 @@ int main() { baseParticle.restitution = 200.0f; std::cout << "Starting Fluid Simulation..." << std::endl; + sim.grid.setLODFalloff(0.001); + sim.grid.setLODMinDistance(4096); - for (int frameIdx = 0; frameIdx < TOTAL_FRAMES; ++frameIdx) { + for (int frameIdx = 0; frameIdx < TOTAL_FRAMES; frameIdx++) { if (frameIdx < (TOTAL_FRAMES * 0.1f)) { if (frameIdx % 1 == 0) { @@ -39,41 +41,45 @@ int main() { sim.applyPhysics(); - if (frameIdx % 1 == 0) { - std::cout << "Rendering ultrafast Frame " << frameIdx << " / " << TOTAL_FRAMES << std::endl; + if (frameIdx % 100 == 0) { + std::cout << "Rendering Frame " << frameIdx << " / " << TOTAL_FRAMES << std::endl; - frame renderedFrame = sim.grid.fastRenderFrame(cam, HEIGHT, WIDTH, frame::colormap::RGB); + frame renderedFrame = sim.grid.renderFrame(cam, HEIGHT, WIDTH, frame::colormap::RGB, 4, 5, true, false); std::stringstream ss; - ss << "output/ufframe_" << std::setw(4) << std::setfill('0') << frameIdx << ".bmp"; + ss << "output/frame_" << std::setw(4) << std::setfill('0') << frameIdx << ".bmp"; BMPWriter::saveBMP(ss.str(), renderedFrame); - - } - if (frameIdx % 10 == 0) { + } else if (frameIdx % 25 == 0) { std::cout << "Rendering quick Frame " << frameIdx << " / " << TOTAL_FRAMES << std::endl; frame renderedFrame = sim.grid.renderFrame(cam, HEIGHT, WIDTH, frame::colormap::RGB, 1, 1, true, false); std::stringstream ss; - ss << "output/qframe_" << std::setw(4) << std::setfill('0') << frameIdx << ".bmp"; + ss << "output/frame_" << std::setw(4) << std::setfill('0') << frameIdx << ".bmp"; BMPWriter::saveBMP(ss.str(), renderedFrame); + } else if (frameIdx % 5 == 0) { + std::cout << "Rendering ultrafast Frame " << frameIdx << " / " << TOTAL_FRAMES << std::endl; - } - if (frameIdx % 50 == 0) { - std::cout << "Rendering Frame " << frameIdx << " / " << TOTAL_FRAMES << std::endl; - - frame renderedFrame = sim.grid.renderFrame(cam, HEIGHT, WIDTH, frame::colormap::RGB, 4, 4, true, false); + frame renderedFrame = sim.grid.fastRenderFrame(cam, HEIGHT / 2, WIDTH / 2, frame::colormap::RGB); + renderedFrame.resize(HEIGHT, WIDTH, frame::interpolation::NEAREST); std::stringstream ss; ss << "output/frame_" << std::setw(4) << std::setfill('0') << frameIdx << ".bmp"; BMPWriter::saveBMP(ss.str(), renderedFrame); - } + if (frameIdx % 500 == 0) { sim.grid.printStats(); } } + sim.grid.printStats(); + frame renderedFrame = sim.grid.renderFrame(cam, HEIGHT, WIDTH, frame::colormap::RGB, 5, 7, true, false); + + std::stringstream ss; + ss << "output/frame_" << std::setw(4) << std::setfill('0') << TOTAL_FRAMES << ".bmp"; + BMPWriter::saveBMP(ss.str(), renderedFrame); + std::cout << "Simulation Complete." << std::endl; FunctionTimer::printStats(FunctionTimer::Mode::ENHANCED); return 0; diff --git a/util/output/frame.hpp b/util/output/frame.hpp index 4930dd8..8a61c90 100644 --- a/util/output/frame.hpp +++ b/util/output/frame.hpp @@ -45,6 +45,13 @@ public: HUFFMAN, RAW }; + + enum class interpolation { + NEAREST, + BILINEAR, + AREA, + LANCZOS4 + }; private: size_t getChannels(colormap fmt) const { switch (fmt) { @@ -562,6 +569,35 @@ public: _data.shrink_to_fit(); } + void resize(size_t newWidth, size_t newHeight, interpolation method = interpolation::NEAREST) { + TIME_FUNCTION; + if (cformat != compresstype::RAW) { + throw std::runtime_error("Cannot resize a compressed frame. Decompress first."); + } + if (newWidth == 0 || newHeight == 0) { + throw std::invalid_argument("Target dimensions must be non-zero."); + } + if (width == newWidth && height == newHeight) { + return; + } + + size_t channels = getChannels(colorFormat); + std::vector newData; + newData.resize(newWidth * newHeight * channels); + + if (method == interpolation::NEAREST) { + resizeNearest(newData, newWidth, newHeight, channels); + } else if (method == interpolation::BILINEAR) { + resizeBilinear(newData, newWidth, newHeight, channels); + } + + _data = std::move(newData); + width = newWidth; + height = newHeight; + + resetState(_data.size()); + } + private: std::vector> sortvecs(std::vector> source) { std::sort(source.begin(), source.end(), [](const std::vector & a, const std::vector & b) {return a.size() > b.size();}); @@ -753,6 +789,72 @@ private: throw std::logic_error("Function not yet implemented"); } + void resizeNearest(std::vector& dst, size_t newW, size_t newH, size_t channels) { + const double x_ratio = (double)width / newW; + const double y_ratio = (double)height / newH; + + #pragma omp parallel for + for (size_t y = 0; y < newH; ++y) { + size_t srcY = static_cast(std::floor(y * y_ratio)); + if (srcY >= height) srcY = height - 1; + + size_t destOffsetBase = y * newW * channels; + size_t srcRowOffset = srcY * width * channels; + + for (size_t x = 0; x < newW; ++x) { + size_t srcX = static_cast(std::floor(x * x_ratio)); + if (srcX >= width) srcX = width - 1; + + size_t srcIndex = srcRowOffset + (srcX * channels); + size_t destIndex = destOffsetBase + (x * channels); + + for (size_t c = 0; c < channels; ++c) { + dst[destIndex + c] = _data[srcIndex + c]; + } + } + } + } + + void resizeBilinear(std::vector& dst, size_t newW, size_t newH, size_t channels) { + const float x_ratio = (width > 1) ? static_cast(width - 1) / (newW - 1) : 0; + const float y_ratio = (height > 1) ? static_cast(height - 1) / (newH - 1) : 0; + + for (size_t y = 0; y < newH; ++y) { + float srcY_f = y * y_ratio; + size_t y_l = static_cast(srcY_f); + size_t y_h = (y_l + 1 < height) ? y_l + 1 : y_l; + float y_weight = srcY_f - y_l; + float y_inv = 1.0f - y_weight; + + size_t destOffsetBase = y * newW * channels; + + for (size_t x = 0; x < newW; ++x) { + float srcX_f = x * x_ratio; + size_t x_l = static_cast(srcX_f); + size_t x_h = (x_l + 1 < width) ? x_l + 1 : x_l; + float x_weight = srcX_f - x_l; + float x_inv = 1.0f - x_weight; + + size_t idx_TL = (y_l * width + x_l) * channels; + size_t idx_TR = (y_l * width + x_h) * channels; + size_t idx_BL = (y_h * width + x_l) * channels; + size_t idx_BR = (y_h * width + x_h) * channels; + size_t destIndex = destOffsetBase + (x * channels); + + for (size_t c = 0; c < channels; ++c) { + float val_TL = _data[idx_TL + c]; + float val_TR = _data[idx_TR + c]; + float val_BL = _data[idx_BL + c]; + float val_BR = _data[idx_BR + c]; + + float top = val_TL * x_inv + val_TR * x_weight; + float bottom = val_BL * x_inv + val_BR * x_weight; + float result = top * y_inv + bottom * y_weight; + dst[destIndex + c] = static_cast(result); + } + } + } + } }; inline std::ostream& operator<<(std::ostream& os, frame& f) { diff --git a/util/sim/fluidsim.hpp b/util/sim/fluidsim.hpp index 89d19a6..83b9cad 100644 --- a/util/sim/fluidsim.hpp +++ b/util/sim/fluidsim.hpp @@ -36,7 +36,7 @@ struct fluidParticle { }; struct gridConfig { - float gridSizeCube = 4096; + float gridSizeCube = 8192; const float SMOOTHING_RADIUS = 1024.0f; const float REST_DENSITY = 0.00005f; const float TIMESTEP = 0.016f; @@ -190,13 +190,13 @@ public: void spawnParticles(fluidParticle toSpawn, int count, bool resize = true) { TIME_FUNCTION; toSpawn.mass = newMass; - float size = toSpawn.mass / 10; + //float size = toSpawn.mass / 10; Eigen::Vector3f color = buildGradient(toSpawn.mass / 1000, gradientmap); for (int i = 0; i < count; i++) { Eigen::Matrix pos = posGen(); int id = nextObjectId++; - grid.set(toSpawn, pos, true, color, size, true, id, (toSpawn.mass > 100) ? true : false, 1); + grid.set(toSpawn, pos, true, color, 10, true, id, (toSpawn.mass > 100) ? true : false, 1); idposMap.emplace(id, pos); } if (resize){ @@ -226,27 +226,18 @@ public: } void applyPressure() { - // std::cout << "points to consider: " << idposMap.size(); size_t pointcounter = 0; for (auto& point : idposMap) { - // std::cout << " point: " << pointcounter++; auto node = grid.find(point.second); if (!node) { - // std::cout << " point is missing" << std::endl; continue; } - // std::cout << " nd: " << node->data.density; - // std::cout << " nr: " << node->data.restitution; - // std::cout << " crd: " << config.REST_DENSITY; node->data.pressure = node->data.restitution * (node->data.density - config.REST_DENSITY); - // std::cout << std::endl << std::endl; } for (auto& point : idposMap) { - // std::cout << " point: " << pointcounter++; auto node = grid.find(point.second); if (!node) { - // std::cout << " point is missing" << std::endl; continue; } Eigen::Vector3f pressureForce = Eigen::Vector3f::Zero(); @@ -256,15 +247,11 @@ public: if (node == neighbor) continue; Eigen::Vector3f rv = node->position - neighbor->position; Eigen::Vector3f gradW = gradW_spiky(rv, config.SMOOTHING_RADIUS); - // std::cout << " np: " << node->data.pressure; - // std::cout << " np2: " << neighbor->data.pressure; - // std::cout << " nd: " << neighbor->data.density; float scalarP = (node->data.pressure + neighbor->data.pressure) / (2.0f * neighbor->data.density); pressureForce -= neighbor->data.mass * scalarP * gradW; } node->data.pressureForce = pressureForce; - // std::cout << " pressure: " << pressureForce.norm() << std::endl; } } @@ -280,28 +267,19 @@ public: Eigen::Vector3f rv = node->position - neighbor->position; Eigen::Vector3f velDiff = neighbor->data.velocity - node->data.velocity; float lapW = lapW_visc(rv, config.SMOOTHING_RADIUS); - // std::cout << "nv: " << node->data.viscosity; - // std::cout << "nm: " << neighbor->data.mass; - // std::cout << "veldiff: " << velDiff.norm(); - // std::cout << "nd: " << neighbor->data.density; - // std::cout << "lapl: " << lapW; viscosityForce += node->data.viscosity * neighbor->data.mass * (velDiff / neighbor->data.density) * lapW; } node->data.viscosityForce = viscosityForce; - // std::cout << "viscosity: " << viscosityForce.norm() << std::endl; } } Eigen::Vector3f applyMutualGravity(std::shared_ptr::NodeData>& node) { Eigen::Vector3f gravityForce = Eigen::Vector3f::Zero(); std::vector::NodeData>> neighbors = grid.findInRadius(node->position, config.SMOOTHING_RADIUS); - for (auto& neighbor : neighbors) { if (node == neighbor) continue; - Eigen::Vector3f dir = neighbor->position - node->position; float dist = dir.norm(); - if (dist > EPSILON) { dir.normalize(); float forceMag = (config.G_ATTRACTION * node->data.mass * neighbor->data.mass) / (dist * dist); @@ -310,6 +288,21 @@ public: } return gravityForce; } + + Eigen::Vector3f applyCenterGravity(std::shared_ptr::NodeData>& node) { + Eigen::Vector3f center = Eigen::Vector3f::Zero(); + Eigen::Vector3f direction = center - node->position; + float distSq = direction.squaredNorm(); + float dist = std::sqrt(distSq); + if (dist < 50.0f) { + dist = 50.0f; + distSq = 2500.0f; + } + direction /= dist; + float centerMass = 5000.0f; + float forceMag = (config.G_ATTRACTION * node->data.mass * centerMass) / distSq; + return direction * forceMag; + } void applyForce() { for (auto& point : idposMap) { @@ -318,16 +311,17 @@ public: Eigen::Vector3f internalForces = node->data.pressureForce + node->data.viscosityForce; Eigen::Vector3f acceleration = (internalForces / node->data.mass); Eigen::Vector3f gravity = applyMutualGravity(node); - // std::cout << "gravity: " << gravity.norm(); + gravity += applyCenterGravity(node); node->data.forceAccumulator = (acceleration + gravity) / node->data.mass; - // node->data.forceAccumulator = node->data.pressureForce + node->data.viscosityForce; } } void replaceLost() { std::vector idsToRemove; for (auto& point : idposMap) { - if (!grid.inGrid(point.second)) { + if (config.gridSizeCube / 2 > point.second[0] > config.gridSizeCube / 2 || + config.gridSizeCube / 2 > point.second[1] > config.gridSizeCube / 2 || + config.gridSizeCube / 2 > point.second[2] > config.gridSizeCube / 2) { idsToRemove.push_back(point.first); } } @@ -354,24 +348,14 @@ public: auto node = grid.find(point.second); if (!node) continue; Eigen::Matrix acceleration = node->data.forceAccumulator; - //std::cout << acceleration; node->data.velocity += acceleration * config.TIMESTEP; Eigen::Matrix newPos = point.second + (node->data.velocity); Eigen::Vector3f oldPos = point.second; - // std::cout << "moving " << (newPos - node->position).norm(); - // grid.update(oldPos, newPos, node->data, node->visible, node->color, node->size, - // node->active, node->objectId, node->light, node->emittance, - // node->refraction, node->reflection, node->shape); if (grid.move(oldPos, newPos)) { auto newpoint = grid.find(newPos); - if (!newpoint) std::cout << "failed to move" << std::endl; - // auto oldpoint = grid.find(oldPos); - // if (!oldpoint) std::cout << "still at old position" << std::endl; idposMap[point.first] = newPos; - // std::cout << "moved" << std::endl; - } // else std::cout << "failed to move" << std::endl; + } } - // std::cout << std::endl; replaceLost(); } };