main.cpp

Go to the documentation of this file.

/// @file main.cpp
/// @brief 3D SPH fluid simulation  -- hose demo with scene rotation
///
/// Two water jets aimed at the centre of the domain:
///   H   -- toggle hot  hose (top-left corner,  80  degC)
///   C   -- toggle cold hose (top-right corner,  5  degC)
///
/// Arrow keys orbit the camera AND physically tilt the scene so that gravity
/// always points along the camera's "down" direction.  Particles therefore
/// fall toward whichever wall is currently at the bottom of the screen.
///
/// Controls:
///   H / C         -- hot / cold hose on-off
///   Arrow keys    -- rotate scene (gravity tilts with pitch)
///   SPACE         -- pause / resume
///   R             -- reset (clear particles, keep hose state)
///   +  /  -       -- more / fewer substeps per frame
///   ESC           -- quit
 
#include "raylib.h"
#include "raymath.h"
#include "fluid3d.hpp"
#include <cmath>
#include <algorithm>
 
static constexpr float kPI = 3.14159265f;
static constexpr float kG  = 9.81f;
 
// --- colour ------------------------------------------------------------------
 
static Color temp_color(float T, float T_min, float T_max) {
    float t = (T - T_min) / (T_max - T_min + 1e-4f);
    t = (t < 0.0f) ? 0.0f : (t > 1.0f ? 1.0f : t);
    return ColorFromHSV((1.0f - t) * 240.0f, 1.0f, 0.95f);
}
 
// --- solver factory ----------------------------------------------------------
 
static physics::FluidParams3D make_params() {
    physics::FluidParams3D p;
    p.h           = 0.05f;
    p.rho0        = 1000.0f;
    p.gamma       = 7;
    p.c0          = 10.0f;
    p.mu          = 10.0f;
    p.mass        = p.rho0 * (0.8f * p.h) * (0.8f * p.h) * (0.8f * p.h);
    p.gx = 0.0f;  p.gy = -kG;  p.gz = 0.0f;
    p.dt          = 0.001f;
    p.xmin = p.ymin = p.zmin = 0.0f;
    p.xmax = p.ymax = p.zmax = 0.8f;
    p.restitution = 0.01f;
    p.alpha_T     = 0.005f;
    p.h_conv      = 8.0f;
    return p;
}
 
// --- vector helpers (raylib Vector3 doesn't have all we need) ----------------
 
static Vector3 v3_norm(Vector3 v) {
    float len = sqrtf(v.x*v.x + v.y*v.y + v.z*v.z);
    if (len < 1e-8f) return {0.0f, 1.0f, 0.0f};
    return {v.x/len, v.y/len, v.z/len};
}
static Vector3 v3_cross(Vector3 a, Vector3 b) {
    return {a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x};
}
 
// --- Hose --------------------------------------------------------------------
 
struct Hose {
    Vector3 source;
    Vector3 dir;       ///< normalised jet direction
    Vector3 perp_u;    ///< two axes spanning the spray cone
    Vector3 perp_v;
    float   temperature;
    float   speed;
    float   cone_half; ///< half-angle in radians
    bool    active = false;
 
    void init(Vector3 src, Vector3 target,
              float temp, float spd, float half_angle) {
        source      = src;
        temperature = temp;
        speed       = spd;
        cone_half   = half_angle;
 
        dir = v3_norm({target.x - src.x, target.y - src.y, target.z - src.z});
 
        // Build orthonormal basis for the cone
        Vector3 ref = (fabsf(dir.y) < 0.9f) ? Vector3{0,1,0} : Vector3{1,0,0};
        perp_u = v3_norm(v3_cross(dir, ref));
        perp_v = v3_cross(dir, perp_u);   // already unit since dir ⊥ perp_u
    }
 
    void emit(physics::FluidSolver3D& solver, int count,
              int max_particles = 1800) const {
        if (!active) return;
        if ((int)solver.particles().size() >= max_particles) return;
 
        const float max_r = tanf(cone_half);
        for (int i = 0; i < count; ++i) {
            float angle = (float)GetRandomValue(0, 10000) / 10000.0f * 2.0f * kPI;
            float r     = (float)GetRandomValue(0, 10000) / 10000.0f * max_r;
            float su = cosf(angle) * r;
            float sv = sinf(angle) * r;
 
            // Spray direction with random cone spread
            Vector3 d = v3_norm({
                dir.x + su * perp_u.x + sv * perp_v.x,
                dir.y + su * perp_u.y + sv * perp_v.y,
                dir.z + su * perp_u.z + sv * perp_v.z
            });
 
            // Small position jitter so particles don't all stack
            float jx = (float)(GetRandomValue(-100, 100)) * 0.00008f;
            float jy = (float)(GetRandomValue(-100, 100)) * 0.00008f;
            float jz = (float)(GetRandomValue(-100, 100)) * 0.00008f;
 
            solver.add_particle(source.x + jx, source.y + jy, source.z + jz,
                                d.x * speed, d.y * speed, d.z * speed,
                                temperature);
        }
    }
 
    // Visual nozzle position
    Vector3 pos() const { return source; }
};
 
// --- main --------------------------------------------------------------------
 
int main() {
    InitWindow(1200, 800, "SPH 3D  -- Hose Demo");
    SetTargetFPS(60);
 
    physics::FluidParams3D params = make_params();
    physics::FluidSolver3D solver(params);
    // Start empty  -- hoses fill the scene
 
    // Scene/camera orbit angles
    //   yaw:   horizontal rotation around world-Y (left/right arrows)
    //   pitch: elevation angle (up/down arrows)
    //
    // Gravity follows camera "down": g = G * (sinY*sinP, -cosP, cosY*sinP)
    // so pitching physically tilts the container  -- pure yaw (vertical rotation)
    // leaves gravity unchanged, which is correct.
    float cam_yaw   =  kPI * 0.75f;  // 135 deg  -- starts from +X+Z corner (matches old view)
    float cam_pitch =  0.45f;         // ~26 deg elevation
    const float CAM_DIST  = 2.1f;
    const float ROT_SPEED = 1.6f;    // rad/s
    const Vector3 center  = {0.4f, 0.4f, 0.4f};
 
    // Cold hose: top-right-front corner -> centre
    Hose cold_hose, hot_hose;
    cold_hose.init({0.77f, 0.76f, 0.02f}, center,  5.0f, 1.8f, kPI / 36.0f); // 5 deg
    // Hot hose:  top-left-front corner -> centre
    hot_hose .init({0.03f, 0.76f, 0.02f}, center, 80.0f, 1.8f, kPI / 36.0f);
 
    Camera3D camera{};
    camera.fovy       = 50.0f;
    camera.projection = CAMERA_PERSPECTIVE;
    camera.target     = center;
 
    auto apply_camera = [&]() {
        float cp = cosf(cam_pitch), sp = sinf(cam_pitch);
        float cy = cosf(cam_yaw),   sy = sinf(cam_yaw);
        camera.position = {
            center.x + CAM_DIST * cp * sy,
            center.y + CAM_DIST * sp,
            center.z + CAM_DIST * cp * cy
        };
        // Camera up = scene's local Y axis (derived from orbit angles)
        camera.up = { -sy * sp, cp, -cy * sp };
    };
 
    auto apply_gravity = [&]() {
        float cp = cosf(cam_pitch), sp = sinf(cam_pitch);
        float cy = cosf(cam_yaw),   sy = sinf(cam_yaw);
        auto& pm = solver.params_mut();
        pm.gx =  kG * sy * sp;
        pm.gy = -kG * cp;
        pm.gz =  kG * cy * sp;
    };
 
    apply_camera();
    apply_gravity();
 
    bool paused   = false;
    int  substeps = 4;
 
    while (!WindowShouldClose()) {
        const float dt_frame = GetFrameTime();
 
        // -- input ---------------------------------------------------------
        if (IsKeyPressed(KEY_SPACE))  paused = !paused;
        if (IsKeyPressed(KEY_R))      solver.clear();
        if (IsKeyPressed(KEY_H))      hot_hose.active  = !hot_hose.active;
        if (IsKeyPressed(KEY_C))      cold_hose.active = !cold_hose.active;
 
        if (IsKeyPressed(KEY_EQUAL)  || IsKeyPressed(KEY_KP_ADD))
            substeps = std::min(substeps + 1, 20);
        if (IsKeyPressed(KEY_MINUS)  || IsKeyPressed(KEY_KP_SUBTRACT))
            substeps = std::max(substeps - 1, 1);
 
        // Arrow keys: rotate scene
        bool rotated = false;
        if (IsKeyDown(KEY_LEFT))  { cam_yaw   -= ROT_SPEED * dt_frame; rotated = true; }
        if (IsKeyDown(KEY_RIGHT)) { cam_yaw   += ROT_SPEED * dt_frame; rotated = true; }
        if (IsKeyDown(KEY_UP))    { cam_pitch += ROT_SPEED * dt_frame; rotated = true; }
        if (IsKeyDown(KEY_DOWN))  { cam_pitch -= ROT_SPEED * dt_frame; rotated = true; }
        cam_pitch = (cam_pitch < -kPI * 0.44f) ? -kPI * 0.44f
                  : (cam_pitch >  kPI * 0.44f) ?  kPI * 0.44f : cam_pitch;
 
        if (rotated) { apply_camera(); apply_gravity(); }
 
        // -- physics -------------------------------------------------------
        if (!paused) {
            cold_hose.emit(solver, 2);
            hot_hose .emit(solver, 2);
            for (int s = 0; s < substeps; ++s)
                solver.step();
        }
 
        // -- render --------------------------------------------------------
        const float T_min = solver.min_temp();
        const float T_max = solver.max_temp();
 
        BeginDrawing();
        ClearBackground(Color{18, 18, 28, 255});
        BeginMode3D(camera);
 
        // Domain box
        DrawCubeWires({0.4f, 0.4f, 0.4f}, 0.8f, 0.8f, 0.8f, Color{60, 60, 80, 200});
 
        // Nozzles (always drawn so user can see where they are)
        DrawSphere(cold_hose.pos(), 0.022f, cold_hose.active ? SKYBLUE : Color{40,40,100,180});
        DrawSphere(hot_hose .pos(), 0.022f, hot_hose.active  ? ORANGE  : Color{100,40,20,180});
 
        // Particles
        for (const physics::Particle3D& p : solver.particles()) {
            Color c = temp_color(p.temperature, T_min, T_max);
            DrawSphere({p.x, p.y, p.z}, 0.016f, c);
        }
 
        EndMode3D();
 
        // -- HUD -----------------------------------------------------------
        DrawFPS(10, 10);
        DrawText(TextFormat("Particles: %d / 1800", (int)solver.particles().size()),
                 10, 36, 20, WHITE);
        DrawText(TextFormat("Substeps/frame: %d  (%.2fx realtime)",
                            substeps, substeps * solver.params().dt * 60.0f),
                 10, 60, 18, LIGHTGRAY);
        DrawText(TextFormat("T range: %.1f - %.1f  degC", T_min, T_max),
                 10, 84, 18, LIGHTGRAY);
 
        // Gravity direction indicator
        const auto& pm = solver.params();
        DrawText(TextFormat("g = (%.2f, %.2f, %.2f) m/s^2", pm.gx, pm.gy, pm.gz),
                 10, 108, 16, Color{160, 160, 160, 255});
 
        // Hose status
        DrawText("H = HOT  hose:",  10, 136, 18, LIGHTGRAY);
        DrawText(hot_hose.active  ? "ON"  : "off",
                 170, 136, 18, hot_hose.active  ? ORANGE  : GRAY);
        DrawText("C = COLD hose:", 10, 158, 18, LIGHTGRAY);
        DrawText(cold_hose.active ? "ON"  : "off",
                 170, 158, 18, cold_hose.active ? SKYBLUE : GRAY);
 
        if (paused) {
            DrawRectangle(0, 0, 1200, 800, Color{0, 0, 0, 100});
            DrawText("PAUSED", 1200/2 - 60, 800/2 - 20, 40, YELLOW);
        }
 
        // Controls bar
        DrawRectangle(0, 770, 1200, 30, Color{0, 0, 0, 180});
        DrawText("H=hot  C=cold  Arrows=rotate scene  SPACE=pause  R=reset  +/-=speed  ESC=quit",
                 8, 778, 14, Color{180, 180, 180, 255});
 
        EndDrawing();
    }
 
    CloseWindow();
    return 0;
}