stupidsimcpp/util/sim/fluidsim.hpp

#include <map>
#include <iostream>
#include <vector>
#include <chrono>
#include <thread>
#include <atomic>
#include <mutex>
#include <cmath>
#include <random>
#include <algorithm>

#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 <GLFW/glfw3.h>
#include "../stb/stb_image.h"

struct fluidParticle {
    Eigen::Matrix<float, 3, 1> velocity;
    Eigen::Matrix<float, 3, 1> acceleration;
    Eigen::Matrix<float, 3, 1> forceAccumulator;
    float density = 0.0f;
    float pressure = 0.0f;
    Eigen::Matrix<float, 3, 1> pressureForce;
    float viscosity = 0.5f;
    Eigen::Matrix<float, 3, 1> viscosityForce;
    float restitution = 5.0f;
    float mass;
};

struct gridConfig {
    float gridSizeCube = 4096;
    const float SMOOTHING_RADIUS = 1024.0f;
    const float REST_DENSITY = 0.00005f;
    const float TIMESTEP = 0.016f;
    const float G_ATTRACTION = 50.0f;
};

Eigen::Matrix<float, 3, 1> posGen() {
    static std::random_device rd;
    static std::mt19937 gen(rd());
    static std::normal_distribution<float> dist(0.0f, 1024.0f);
    
    return Eigen::Matrix<float, 3, 1>(dist(gen),dist(gen),dist(gen));
}

Eigen::Matrix<float, 3, 1> velGen() {
    static std::random_device rd;
    static std::mt19937 gen(rd());
    static std::normal_distribution<float> dist(0.0f, 1.0f);
    
    return Eigen::Matrix<float, 3, 1>(dist(gen),dist(gen),dist(gen));
}

float W_poly6(Eigen::Vector3f rv, float h) {
    float r = rv.squaredNorm();
    
    if (r > h || r < 0) return 0;
    float factor = 315 / (64 * M_PI * pow(h, 9));
    float m = pow(h*h-r, 3);
    return factor * m;
}

Eigen::Vector3f gradW_poly6(Eigen::Vector3f rv, float h) {
    float r = rv.squaredNorm();
    float h2 = h * h;

    if (r > h2 || r < 0) {
        return Eigen::Vector3f::Zero();
    }

    float m = -6 * pow(h*h-r, 2);
    float factor = -945.0f / (32.0f * M_PI * std::pow(h, 9));

    return factor * m * rv;
}

float lapW_poly6(Eigen::Vector3f rv, float h) {
    float r = rv.squaredNorm();
    float h2 = h * h;

    if (r > h || r < 0) return 0;

    float m = h2 - r;
    float term2 = 3.0f * h2 - 7.0f * r;
    float factor = -945.0f / (32.0f * M_PI * std::pow(h, 9));

    return factor * m * term2;
}

float W_spiky(Eigen::Vector3f rv, float h) {
    float r = rv.norm();
    if (r > h || r < 0) return 0;

    float coeff = pow(r-h, 3);
    float factor = 15 / (M_PI * pow(h, 6));

    return factor * coeff;
}

Eigen::Vector3f gradW_spiky(Eigen::Vector3f rv, float h) {
    float r = rv.norm();

    if (r > h || r < 0) {
        return Eigen::Vector3f::Zero();
    }

    float diff = h - r;
    float coeff = -45.0f / (M_PI * std::pow(h, 6));
    
    Eigen::Vector3f direction = rv / r;

    return coeff * std::pow(diff, 2) * direction;
}

float W_visc(Eigen::Vector3f rv, float h) {
    float r = rv.norm();

    if (r > h || r < 0) return 0;

    float r2 = r * r;
    float r3 = r2 * r;
    float h3 = h * h * h;
    
    float coeff = 15.0f / (2.0f * M_PI * h3);
    float term = (-0.5f * r3 / h3) + (r2 / (h * h)) + (h / (2.0f * r)) - 1.0f;

    return coeff * term;
}

float lapW_visc(Eigen::Vector3f rv, float h) {
    float r = rv.norm();

    if (r > h || r < 0) return 0;

    float diff = h - r;
    float coeff = 45.0f / (M_PI * std::pow(h, 6));

    return coeff * diff;
}

Eigen::Vector3f buildGradient(float value, const std::map<float, Eigen::Vector3f>& 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<size_t, Eigen::Matrix<float, 3, 1>> idposMap;
    float newMass = 1000;
    gridConfig config;
    int nextObjectId = 0;
    std::map<float, Eigen::Vector3f> gradientmap;
    float closeThresh = 0.01f * config.SMOOTHING_RADIUS;
public:
    Octree<fluidParticle> 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, bool resize = true) {
        TIME_FUNCTION;
        toSpawn.mass = newMass;
        float size = toSpawn.mass / 10;
        Eigen::Vector3f color = buildGradient(toSpawn.mass / 1000, gradientmap);
        for (int i = 0; i < count; i++) {
            Eigen::Matrix<float, 3, 1> pos = posGen();
            int id = nextObjectId++;
            
            grid.set(toSpawn, pos, true, color, size, true, id, (toSpawn.mass > 100) ? true : false, 1);
            idposMap.emplace(id, pos);
        }
        if (resize){
            newMass *= 0.999f;
        }
    }

    float sphKernel(Eigen::Vector3f rv) {
        float r = rv.norm();
        if (r < closeThresh) return W_spiky(rv, config.G_ATTRACTION);
        else return W_poly6(rv,config.G_ATTRACTION);
    }

    void computeDensities() {
        for (auto& point : idposMap) {
            float densSum = 0;
            auto node = grid.find(point.second);
            if (!node) continue;
            std::vector<std::shared_ptr<Octree<fluidParticle>::NodeData>> neighbors = grid.findInRadius(point.second, config.SMOOTHING_RADIUS);
            for (auto& neighbor : neighbors) {
                Eigen::Vector3f rv = node->position - neighbor->position;
                float w = sphKernel(rv);
                densSum += neighbor->data.mass * w;
            }
            node->data.density = densSum;
        }
    }

    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();

            std::vector<std::shared_ptr<Octree<fluidParticle>::NodeData>> neighbors = grid.findInRadius(point.second, config.SMOOTHING_RADIUS);
            for (auto& neighbor : neighbors) {
                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;
        }
    }

    void applyViscosity() {
        for (auto& point : idposMap) {
            auto node = grid.find(point.second);
            if (!node) continue;
            Eigen::Vector3f viscosityForce = Eigen::Vector3f::Zero();
            if (node->data.velocity == Eigen::Vector3f::Zero()) node->data.velocity = velGen();
            
            std::vector<std::shared_ptr<Octree<fluidParticle>::NodeData>> neighbors = grid.findInRadius(point.second, config.SMOOTHING_RADIUS);
            for (auto& neighbor : neighbors) {
                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<Octree<fluidParticle>::NodeData>& node) {
        Eigen::Vector3f gravityForce = Eigen::Vector3f::Zero();
        std::vector<std::shared_ptr<Octree<fluidParticle>::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);
                gravityForce += dir * forceMag;
            }
        }
        return gravityForce;
    }

    void applyForce() {
        for (auto& point : idposMap) {
            auto node = grid.find(point.second);
            if (!node) continue;
            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();
            node->data.forceAccumulator = (acceleration + gravity) / node->data.mass;
            // node->data.forceAccumulator = node->data.pressureForce + node->data.viscosityForce;
        }
    }

    void replaceLost() {
        std::vector<size_t> idsToRemove;
        for (auto& point : idposMap) { 
            if (!grid.inGrid(point.second)) {
                idsToRemove.push_back(point.first);
            }
        }

        for (size_t id : idsToRemove) {
            grid.remove(idposMap[id]);
            idposMap.erase(id);
        }

        if (!idsToRemove.empty()) {
            fluidParticle newParticles;
            spawnParticles(newParticles, idsToRemove.size(), false);
        }
    }

    void applyPhysics() {
        TIME_FUNCTION;
        computeDensities();
        applyPressure();
        applyViscosity();
        applyForce();

        for (auto& point : idposMap) {
            auto node = grid.find(point.second);
            if (!node) continue;
            Eigen::Matrix<float, 3, 1> acceleration = node->data.forceAccumulator;
            //std::cout << acceleration;
            node->data.velocity += acceleration * config.TIMESTEP;
            Eigen::Matrix<float, 3, 1> 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();
    }
};