From 341066a0cfccea759d6000ef9bf7cb42ddb156ae Mon Sep 17 00:00:00 2001 From: Yggdrasil75 Date: Wed, 4 Feb 2026 15:00:05 -0500 Subject: [PATCH] pushing --- makefile | 2 +- tests/fluidsim.cpp | 75 +++++++++++++++ util/grid/grid3eigen.hpp | 4 + util/sim/fluidsim.cpp | 83 ---------------- util/sim/fluidsim.hpp | 198 +++++++++++++++++++++++++++++++++++++++ 5 files changed, 278 insertions(+), 84 deletions(-) create mode 100644 tests/fluidsim.cpp delete mode 100644 util/sim/fluidsim.cpp create mode 100644 util/sim/fluidsim.hpp diff --git a/makefile b/makefile index d3ef09d..cf3ef30 100644 --- a/makefile +++ b/makefile @@ -34,7 +34,7 @@ endif CXXFLAGS = $(BASE_CXXFLAGS) $(SIMD_CXXFLAGS) # Source files -SRC := $(SRC_DIR)/g3etest.cpp +SRC := $(SRC_DIR)/fluidsim.cpp #SRC := $(SRC_DIR)/g2chromatic2.cpp SRC += $(IMGUI_DIR)/imgui.cpp $(IMGUI_DIR)/imgui_demo.cpp $(IMGUI_DIR)/imgui_draw.cpp $(IMGUI_DIR)/imgui_tables.cpp $(IMGUI_DIR)/imgui_widgets.cpp SRC += $(IMGUI_DIR)/backends/imgui_impl_glfw.cpp $(IMGUI_DIR)/backends/imgui_impl_opengl3.cpp diff --git a/tests/fluidsim.cpp b/tests/fluidsim.cpp new file mode 100644 index 0000000..4fb8d14 --- /dev/null +++ b/tests/fluidsim.cpp @@ -0,0 +1,75 @@ +#include "../util/sim/fluidsim.hpp" +#include "../util/timing_decorator.cpp" + +int main() { + fluidSim sim; + + // Simulation settings + const int TOTAL_FRAMES = 100000; + const int WIDTH = 800; + const int HEIGHT = 600; + + // Setup Camera + Camera cam; + cam.origin = Eigen::Vector3f(0.0f, 400.0f, -1800.0f); // Back and slightly up + cam.direction = (Eigen::Vector3f(0.0f, 0.0f, 0.0f) - cam.origin).normalized(); // Look at center + cam.up = Eigen::Vector3f(0.0f, 1.0f, 0.0f); + cam.fov = 60.0f; + + // Base particle template + fluidParticle baseParticle; + baseParticle.velocity = Eigen::Vector3f::Zero(); + baseParticle.acceleration = Eigen::Vector3f::Zero(); + baseParticle.forceAccumulator = Eigen::Vector3f::Zero(); + baseParticle.pressure = 0.0f; + baseParticle.viscosity = 8.0f; + baseParticle.restitution = 200.0f; + + std::cout << "Starting Fluid Simulation..." << std::endl; + + for (int frameIdx = 0; frameIdx < TOTAL_FRAMES; ++frameIdx) { + + // --- Spawning Logic --- + // Spawn throughout the simulation (e.g., until 90% completion) + if (frameIdx % 10 == 0) { + if (frameIdx < (TOTAL_FRAMES * 0.9f)) { + float t = static_cast(frameIdx) / (TOTAL_FRAMES * 0.9f); + + // Linear interpolation for Mass: Heavy (10.0) -> Light (0.5) + float currentMass = (1.0f - t) * 10.0f + t * 0.5f; + + // Linear interpolation for Size: Large (12.0) -> Small (2.0) + float currentSize = (1.0f - t) * 12.0f + t * 2.0f; + + // Number of particles to spawn this frame (can increase as they get smaller) + int spawnCount = 2 + static_cast(t * 4); + + sim.spawnParticles(baseParticle, spawnCount, currentSize, currentMass); + } + } + + // --- Physics Step --- + sim.applyPhysics(); + + // --- Render & Save --- + if (frameIdx % 50 == 0) { + std::cout << "Rendering Frame " << frameIdx << " / " << TOTAL_FRAMES << std::endl; + + // Generate LODs for faster/better rendering + sim.grid.generateLODs(); + + // Render frame + frame renderedFrame = sim.grid.renderFrame(cam, HEIGHT, WIDTH, frame::colormap::RGB, 4, 4, true, false); + + // Save to BMP + std::stringstream ss; + ss << "output/frame_" << std::setw(4) << std::setfill('0') << frameIdx << ".bmp"; + BMPWriter::saveBMP(ss.str(), renderedFrame); + sim.grid.printStats(); + } + } + + std::cout << "Simulation Complete." << std::endl; + FunctionTimer::printStats(FunctionTimer::Mode::ENHANCED); + return 0; +} \ No newline at end of file diff --git a/util/grid/grid3eigen.hpp b/util/grid/grid3eigen.hpp index 57f29ed..9ceee59 100644 --- a/util/grid/grid3eigen.hpp +++ b/util/grid/grid3eigen.hpp @@ -646,6 +646,10 @@ public: return searchNode(root_.get()); } + bool inGrid(PointType pos) { + return root_.contains(pos); + } + bool remove(const PointType& pos, float tolerance = 0.0001f) { bool found = false; std::function removeNode = [&](OctreeNode* node) -> bool { diff --git a/util/sim/fluidsim.cpp b/util/sim/fluidsim.cpp deleted file mode 100644 index d016746..0000000 --- a/util/sim/fluidsim.cpp +++ /dev/null @@ -1,83 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "../util/grid/grid3eigen.hpp" -#include "../util/output/bmpwriter.hpp" -#include "../util/output/frame.hpp" -#include "../util/timing_decorator.cpp" -#include "../util/noise/pnoise2.hpp" -#include "../util/noise/pnoise.cpp" -#include "../util/output/aviwriter.hpp" - -#include "../imgui/imgui.h" -#include "../imgui/backends/imgui_impl_glfw.h" -#include "../imgui/backends/imgui_impl_opengl3.h" -#include -#include "../stb/stb_image.h" - - -struct fluidParticle { - Eigen::Matrix pos; - Eigen::Matrix velocity; - Eigen::Matrix acceleration; - Eigen::Matrix forceAccumulator; - float pressure; - float viscosity; - float restitution; - float mass; -}; - -struct gridDefaults { - float gridSizeCube = 1024; -}; - -Eigen::Matrix posGen() { - static std::random_device rd; - static std::mt19937 gen(rd()); - static std::normal_distribution dist(0.0f, 1024.0f); - - return Eigen::Matrix( - dist(gen), - dist(gen), - dist(gen) - ); -} -class fluidSim { - std::unordered_map> idposMap; - Octree grid; - - void spawnParticles(fluidParticle toSpawn, int count) { - for (int i = 0; i < count; i++) { - Eigen::Matrix pos = posGen(); - grid.set(toSpawn, pos, true, {toSpawn.mass / 1000, toSpawn.mass / 1000, toSpawn.mass / 1000}, 1, true, 1, false); - - } - } - - void applyPhysics() { - for (auto& meh : idposMap) { - std::shared_ptr::NodeData> mehdata = grid.find(meh.second); - std::vector::NodeData>> neighbors = grid.findInRadius(meh.second, mehdata->data.mass); - ///TODO: - // apply velocity to have particle move - // use mass of neighboring particles to change direction towards them - // ressure pushes particles away if they are too close - // resitution controls how much pressure builds up based on how close particles are - // have to find a way to have a clump use their combined mass instead of individual mass. - //if a particle leaves the grid, destroy it and spawn a new one - - for (auto& neighbor : neighbors) { - - - } - } - } -}; \ No newline at end of file diff --git a/util/sim/fluidsim.hpp b/util/sim/fluidsim.hpp new file mode 100644 index 0000000..504a5d3 --- /dev/null +++ b/util/sim/fluidsim.hpp @@ -0,0 +1,198 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../util/grid/grid3eigen.hpp" +#include "../util/output/bmpwriter.hpp" +#include "../util/output/frame.hpp" +#include "../util/noise/pnoise2.hpp" +#include "../util/noise/pnoise.cpp" +#include "../util/output/aviwriter.hpp" + +#include "../imgui/imgui.h" +#include "../imgui/backends/imgui_impl_glfw.h" +#include "../imgui/backends/imgui_impl_opengl3.h" +#include +#include "../stb/stb_image.h" + +struct fluidParticle { + Eigen::Matrix pos; + Eigen::Matrix velocity; + Eigen::Matrix acceleration; + Eigen::Matrix forceAccumulator; + float density; + float pressure; + float viscosity; + float restitution; + float mass; +}; + +struct gridConfig { + float gridSizeCube = 1024; + const float SMOOTHING_RADIUS = 32.0f; + const float REST_DENSITY = 0.5f; + const float TIMESTEP = 0.016f; + const float G_ATTRACTION = 50.0f; +}; + +Eigen::Matrix posGen() { + static std::random_device rd; + static std::mt19937 gen(rd()); + static std::normal_distribution dist(0.0f, 1024.0f); + + return Eigen::Matrix(dist(gen),dist(gen),dist(gen)); +} + +// Math Helpers for SPH Kernels +float poly6Kernel(float r, float h) { + if (r < 0 || r > h) return 0.0f; + float h2 = h * h; + float r2 = r * r; + float diff = h2 - r2; + float coef = 315.0f / (64.0f * M_PI * pow(h, 9)); + return coef * diff * diff * diff; +} + +Eigen::Vector3f spikyGradientKernel(Eigen::Vector3f r_vec, float r, float h) { + if (r <= 0 || r > h) return Eigen::Vector3f::Zero(); + float diff = h - r; + float coef = -45.0f / (M_PI * pow(h, 6)); + return r_vec.normalized() * coef * diff * diff; +} + +Eigen::Vector3f buildGradient(float value, const std::map& gradientKeys) { + auto exactMatch = gradientKeys.find(value); + if (exactMatch != gradientKeys.end()) { + return exactMatch->second; + } + auto lower = gradientKeys.lower_bound(value); + if (lower == gradientKeys.begin()) { + return gradientKeys.begin()->second; + } + if (lower == gradientKeys.end()) { + return gradientKeys.rbegin()->second; + } + auto upper = lower; + lower = std::prev(lower); + float key1 = lower->first; + float key2 = upper->first; + const Eigen::Vector3f& color1 = lower->second; + const Eigen::Vector3f& color2 = upper->second; + float t = (value - key1) / (key2 - key1); + t = std::clamp(t, 0.0f, 1.0f); + return color1 + t * (color2 - color1); +} + +class fluidSim { +private: + std::unordered_map> idposMap; + float newMass = 1000; + gridConfig config; + int nextObjectId = 0; + std::map gradientmap; +public: + Octree grid; + + fluidSim() : grid({-config.gridSizeCube, -config.gridSizeCube, -config.gridSizeCube}, {config.gridSizeCube, config.gridSizeCube, config.gridSizeCube}) { + grid.setBackgroundColor({0.1f, 0.1f, 0.2f}); + gradientmap.emplace(0.0, Eigen::Vector3f(1, 0, 0)); + gradientmap.emplace(0.5, Eigen::Vector3f(0, 1, 0)); + gradientmap.emplace(1.0, Eigen::Vector3f(0, 0, 1)); + + } + + void spawnParticles(fluidParticle toSpawn, int count) { + TIME_FUNCTION; + for (int i = 0; i < count; i++) { + Eigen::Matrix pos = posGen(); + int id = nextObjectId++; + + Eigen::Vector3f color = buildGradient(toSpawn.mass / 1000, gradientmap); + grid.set(toSpawn, pos, true, color, 1, true, id, false); + idposMap[id] = pos; + } + newMass -= newMass / 10; + } + + void applyPressureDensity() { + #pragma omp parallel for + for (auto& point : idposMap) { + auto node = grid.find(point.second); + if (!node) continue; + std::vector::NodeData>> neighbors = grid.findInRadius(point.second, config.SMOOTHING_RADIUS * node->data.mass); + float density = 0.0f; + for (const auto& neighbor : neighbors) { + if (node == neighbor) continue; + float dist = (point.second-neighbor->position).norm(); + density += neighbor->data.mass * poly6Kernel(dist, config.SMOOTHING_RADIUS * node->data.mass); + } + node->data.density = density; + float pressure = node->data.restitution * (density - config.REST_DENSITY); + node->data.pressure = std::max(pressure, 0.0f); + } + } + + void applyForce() { + #pragma omp parallel for + for (auto& point : idposMap) { + auto node = grid.find(point.second); + if (!node) continue; + Eigen::Vector3f totalForce = -point.second.normalized() * 10.0f * node->data.mass; + std::vector::NodeData>> neighbors = grid.findInRadius(point.second, config.SMOOTHING_RADIUS * node->data.mass); + + for (const auto& neighbor : neighbors) { + if (node == neighbor) continue; + Eigen::Vector3f diff = point.second - neighbor->position; + float dist = diff.norm(); + if (dist < EPSILON) continue; + float densj = neighbor->data.density; + if (densj < EPSILON) densj = EPSILON; + float pressj = neighbor->data.pressure; + + Eigen::Vector3f pressureForceVec = -neighbor->data.mass * ((node->data.pressure + pressj) / (2.0f * densj)) * spikyGradientKernel(diff, dist, config.SMOOTHING_RADIUS); + totalForce += pressureForceVec; + Eigen::Vector3f velDiff = neighbor->data.velocity - node->data.velocity; + float viscositCoef = node->data.viscosity * neighbor->data.mass * (1.0f / densj) * poly6Kernel(dist, config.SMOOTHING_RADIUS); + totalForce += velDiff * viscositCoef; + + float clampDist = std::max(dist, 5.0f); + Eigen::Vector3f dirToNeighbor = (neighbor->position - node->position).normalized(); + float attract = (config.G_ATTRACTION * node->data.mass * neighbor->data.mass) / (clampDist * clampDist); + totalForce += dirToNeighbor * attract; + } + node->data.forceAccumulator = totalForce; + } + } + + void replaceLost() { + int toReplace = 0; + #pragma omp parallel for + for (auto& point : idposMap) { + if (!grid.inGrid(point.second)) { + grid.remove(point.second); + toReplace++; + } + } + fluidParticle newParticles; + newParticles.mass = newMass; + spawnParticles(newParticles, toReplace); + } + + void applyPhysics() { + TIME_FUNCTION; + ///TODO: + // apply velocity to have particle move + // use mass of neighboring particles to change direction towards them + // ressure pushes particles away if they are too close + // resitution controls how much pressure builds up based on how close particles are + // have to find a way to have a clump use their combined mass instead of individual mass. + //if a particle leaves the grid, destroy it and spawn a new one + } +};