main.cpp

Go to the documentation of this file.

/// @file apps/em_demo/main.cpp
///
/// DC current through an aluminum rod + movable permanent magnet.
///
/// Physics:
///   Rod B field  -- precomputed once via div(sigmagradphi)=0 -> J -> laplacianA=-mu_0J -> B=curlA
///   Magnet       -- magnetic dipole field, evaluated analytically every frame:
///                 B_dip = mu_0/4pi * [3(m*r_hat)r_hat - m] / r^3   (m = M0*y_hat)
///   Total B      -- B_rod (sampled from grid) + B_dip (analytical)
///
/// Controls:
///   LMB drag           -- orbit camera
///   RMB drag           -- pan
///   Scroll             -- zoom
///   Shift + LMB drag   -- move magnet (horizontal -> XZ, vertical -> Y)
#include "field.hpp"
 
#include <raylib.h>
#include <raymath.h>
 
#include <chrono>
#include <cmath>
#include <cstdio>
#include <vector>
 
// Scene parameters
static constexpr int   NX = 32, NY = 32, NZ = 32;
static constexpr float DX = 0.05f;
static constexpr float SIM_DOMAIN = NX * DX;   // 1.6 m
 
static constexpr int   ROD_CX  = NX / 2;
static constexpr int   ROD_CZ  = NZ / 2;
static constexpr int   ROD_R   = 4;        // cells (0.2 m)
static constexpr float ROD_SIG = 1.0f;
static constexpr float BG_SIG  = 1e-6f;
static constexpr float V_TOP   = 1.0f;
static constexpr float V_BOT   = 0.0f;
 
// Magnet dipole moment magnitude [A*m^2]  -- tuned so dipole is
// visible relative to the rod field at ~0.3 m distance
static constexpr float MAG_M0  = 0.01f;
static constexpr float MAG_VIS_R = 0.07f;  // drawn sphere radius
 
static constexpr int   STRIDE  = 4;          // sample every N cells in x,y,z
static constexpr float MAX_LEN = DX * 2.2f;  // longest arrow in world units
 
// Helpers
static bool inside_rod(int i, int k) {
    int di = i - ROD_CX, dk = k - ROD_CZ;
    return di*di + dk*dk <= ROD_R*ROD_R;
}
static int flat(int i, int j, int k) { return k*NY*NX + j*NX + i; }
 
static Color field_color(float t) {
    t = t < 0.0f ? 0.0f : (t > 1.0f ? 1.0f : t);
    return ColorFromHSV((1.0f - t) * 240.0f, 1.0f, 0.95f);
}
 
static void draw_arrow(Vector3 origin, Vector3 dir, float len, Color col) {
    Vector3 tip = { origin.x + dir.x*len,
                    origin.y + dir.y*len,
                    origin.z + dir.z*len };
    DrawLine3D(origin, tip, col);
    DrawSphere(tip, len * 0.12f, col);
}
 
// Dipole field:  B_dip = mu_0/4pi * [3(m*r)r - m*r^2] / r^5
// with m = M0*y_hat  ->  m*r = M0*ry
static void dipole_B(float wx, float wy, float wz,
                     float mx, float my, float mz,
                     float& bx, float& by, float& bz) {
    float rx = wx - mx, ry = wy - my, rz = wz - mz;
    float r2 = rx*rx + ry*ry + rz*rz;
    if (r2 < 1e-5f) { bx = by = bz = 0.0f; return; }
    float r5   = r2 * r2 * sqrtf(r2);
    float mdr  = MAG_M0 * ry;          // m*r  (m = M0*y_hat)
    constexpr float C = 1e-7f;         // mu_0/4pi
    bx = C * 3.0f * mdr * rx / r5;
    by = C * MAG_M0 * (3.0f*ry*ry - r2) / r5;
    bz = C * 3.0f * mdr * rz / r5;
}
 
// Main
int main() {
    physics::ScalarField3D sigma(NX, NY, NZ, DX);
    for (int k = 0; k < NZ; ++k)
    for (int j = 0; j < NY; ++j)
    for (int i = 0; i < NX; ++i)
        sigma.grid().set(i, j, k, inside_rod(i,k) ? ROD_SIG : BG_SIG);
 
    std::vector<physics::ElectrodeBC> bcs;
    for (int k = 0; k < NZ; ++k)
    for (int i = 0; i < NX; ++i) {
        if (!inside_rod(i, k)) continue;
        bcs.push_back({ flat(i, 0,    k), V_BOT });
        bcs.push_back({ flat(i, NY-1, k), V_TOP });
    }
 
    printf("Solving... ");
    fflush(stdout);
    auto t0 = std::chrono::steady_clock::now();
 
    physics::ScalarField3D phi(NX, NY, NZ, DX);
    auto r_phi = physics::ElectricSolver::solve_potential(phi, sigma, bcs);
 
    physics::VectorField3D J = physics::MagneticSolver::current_density(sigma, phi);
    physics::VectorField3D B = physics::MagneticSolver::solve_magnetic_field(J);
 
    double solve_ms = std::chrono::duration<double,std::milli>(
                          std::chrono::steady_clock::now() - t0).count();
    printf("done in %.0f ms\n", solve_ms);
 
    const int W = 1280, H = 800;
    InitWindow(W, H, "EM Demo  -- rod + magnet");
    SetTargetFPS(60);
 
    Camera3D camera  = {};
    camera.up        = { 0.0f, 1.0f, 0.0f };
    camera.fovy      = 45.0f;
    camera.projection = CAMERA_PERSPECTIVE;
 
    Vector3 orbit_target = { SIM_DOMAIN*0.5f, SIM_DOMAIN*0.5f, SIM_DOMAIN*0.5f };
    float   azimuth   =  0.785f;
    float   elevation =  0.524f;
    float   distance  = SIM_DOMAIN * 2.6f;
 
    auto update_cam = [&]() {
        camera.position = {
            orbit_target.x + distance * cosf(elevation) * sinf(azimuth),
            orbit_target.y + distance * sinf(elevation),
            orbit_target.z + distance * cosf(elevation) * cosf(azimuth)
        };
        camera.target = orbit_target;
    };
 
    float mag_x = ROD_CX * DX + 0.45f;
    float mag_y = SIM_DOMAIN * 0.5f;
    float mag_z = ROD_CZ * DX;
 
    // Render loop
    while (!WindowShouldClose()) {
        bool shift = IsKeyDown(KEY_LEFT_SHIFT) || IsKeyDown(KEY_RIGHT_SHIFT);
 
        // --- Mouse drag ---
        if (IsMouseButtonDown(MOUSE_LEFT_BUTTON)) {
            Vector2 d = GetMouseDelta();
            if (shift) {
                // Move magnet: horizontal -> along camera's right (XZ), vertical -> Y
                Vector3 view  = Vector3Normalize(
                    Vector3Subtract(camera.target, camera.position));
                Vector3 right = Vector3Normalize(
                    Vector3CrossProduct(view, camera.up));
                // Project right into XZ so magnet stays on a horizontal plane
                // unless the user drags vertically
                float spd = distance * 0.0012f;
                mag_x += right.x * d.x * spd;
                mag_z += right.z * d.x * spd;
                mag_y -= d.y * spd;
            } else {
                azimuth   -= d.x * 0.005f;
                elevation += d.y * 0.005f;
                if (elevation >  1.55f) elevation =  1.55f;
                if (elevation < -1.55f) elevation = -1.55f;
            }
        }
 
        // --- Scroll zoom ---
        float scroll = GetMouseWheelMove();
        if (scroll != 0.0f) {
            distance -= scroll * distance * 0.08f;
            if (distance < 0.05f) distance = 0.05f;
        }
 
        // --- Right drag pan ---
        if (IsMouseButtonDown(MOUSE_RIGHT_BUTTON)) {
            Vector2 d = GetMouseDelta();
            float spd = distance * 0.001f;
            float raz = azimuth + 1.5708f;
            orbit_target.x += (cosf(raz) * d.x
                - sinf(elevation) * sinf(azimuth) * d.y) * spd;
            orbit_target.y -= cosf(elevation) * d.y * spd;
            orbit_target.z += (sinf(raz) * d.x
                - sinf(elevation) * cosf(azimuth) * d.y) * spd;
        }
 
        update_cam();
 
        // Compute B_max over 3-D sample grid (rod + dipole combined)
        float B_max = 1e-30f;
        for (int j = 0; j < NY; j += STRIDE)
        for (int k = 0; k < NZ; k += STRIDE)
        for (int i = 0; i < NX; i += STRIDE) {
            float wx = i*DX, wy = j*DX, wz = k*DX;
            float bx = (float)B.x.grid()(i, j, k);
            float by = (float)B.y.grid()(i, j, k);
            float bz = (float)B.z.grid()(i, j, k);
            float dbx, dby, dbz;
            dipole_B(wx, wy, wz, mag_x, mag_y, mag_z, dbx, dby, dbz);
            bx += dbx; by += dby; bz += dbz;
            float mag = sqrtf(bx*bx + by*by + bz*bz);
            if (mag > B_max) B_max = mag;
        }
 
        // Draw
        BeginDrawing();
        ClearBackground({ 15, 15, 20, 255 });
 
        BeginMode3D(camera);
 
        // Domain wireframe
        DrawCubeWires({ SIM_DOMAIN*0.5f, SIM_DOMAIN*0.5f, SIM_DOMAIN*0.5f },
                      SIM_DOMAIN, SIM_DOMAIN, SIM_DOMAIN, { 60, 60, 80, 255 });
 
        // Aluminum rod
        float rod_wr = ROD_R * DX;
        DrawCylinder({ ROD_CX*DX, 0.0f, ROD_CZ*DX },
                     rod_wr, rod_wr, SIM_DOMAIN, 24, { 180, 185, 195, 220 });
        DrawCylinderWires({ ROD_CX*DX, 0.0f, ROD_CZ*DX },
                          rod_wr, rod_wr, SIM_DOMAIN, 24, { 100, 110, 130, 255 });
        // Electrode caps
        DrawCylinder({ ROD_CX*DX, SIM_DOMAIN - DX*0.5f, ROD_CZ*DX },
                     rod_wr, rod_wr, DX*0.5f, 24, { 220, 60, 60, 255 });
        DrawCylinder({ ROD_CX*DX, 0.0f, ROD_CZ*DX },
                     rod_wr, rod_wr, DX*0.5f, 24, { 60, 80, 220, 255 });
 
        // Magnet sphere  -- N pole (top, red) and S pole (bottom, blue)
        DrawSphere({ mag_x, mag_y + MAG_VIS_R*0.5f, mag_z }, MAG_VIS_R, { 210, 60, 60, 255 });
        DrawSphere({ mag_x, mag_y - MAG_VIS_R*0.5f, mag_z }, MAG_VIS_R, { 60, 80, 210, 255 });
        // Dipole moment axis (white line)
        DrawLine3D({ mag_x, mag_y - MAG_VIS_R*1.6f, mag_z },
                   { mag_x, mag_y + MAG_VIS_R*1.6f, mag_z }, WHITE);
 
        // B field arrows  -- full 3D volume sample grid
        // Length uses log scale so weak far-field arrows are still visible
        // but strong near-field arrows don't dominate. Skip near-zero.
        static constexpr float B_THRESH = 0.002f;  // fraction of B_max to skip
        for (int j = 0; j < NY; j += STRIDE)
        for (int k = 0; k < NZ; k += STRIDE)
        for (int i = 0; i < NX; i += STRIDE) {
            float wx = i*DX, wy = j*DX, wz = k*DX;
            float bx = (float)B.x.grid()(i, j, k);
            float by = (float)B.y.grid()(i, j, k);
            float bz = (float)B.z.grid()(i, j, k);
            float dbx, dby, dbz;
            dipole_B(wx, wy, wz, mag_x, mag_y, mag_z, dbx, dby, dbz);
            bx += dbx; by += dby; bz += dbz;
 
            float mag = sqrtf(bx*bx + by*by + bz*bz);
            float t   = mag / B_max;
            if (t < B_THRESH) continue;
 
            // log scale: strong fields -> long arrows, weak -> short but visible
            float t_log = logf(1.0f + (t / B_THRESH - 1.0f)) /
                          logf(1.0f + (1.0f / B_THRESH - 1.0f));
            float len = MAX_LEN * (0.15f + 0.85f * t_log);
            draw_arrow({ wx, wy, wz }, { bx/mag, by/mag, bz/mag }, len,
                       field_color(t));
        }
 
        EndMode3D();
 
        // HUD
        DrawRectangle(8, 8, 440, 148, { 0, 0, 0, 160 });
        DrawText("Aluminum Rod + Permanent Magnet", 16, 14, 18, WHITE);
        char buf[128];
        snprintf(buf, sizeof(buf), "Solve: %.0f ms | phi iters=%llu res=%.1e",
                 solve_ms, (unsigned long long)r_phi.iterations, r_phi.residual);
        DrawText(buf, 16, 38, 14, LIGHTGRAY);
        snprintf(buf, sizeof(buf), "Magnet (%.2f, %.2f, %.2f)  M0=%.3f A*m^2",
                 (double)mag_x, (double)mag_y, (double)mag_z, (double)MAG_M0);
        DrawText(buf, 16, 58, 14, LIGHTGRAY);
        snprintf(buf, sizeof(buf), "B_max = %.3e T (rod+dipole)", (double)B_max);
        DrawText(buf, 16, 78, 14, LIGHTGRAY);
        DrawText("LMB: orbit  RMB: pan  Scroll: zoom", 16, 100, 14, GRAY);
        DrawText("Shift+LMB: move magnet (XZ + Y)", 16, 118, 14,
                 shift ? YELLOW : GRAY);
 
        // Colour legend
        for (int px = 0; px < 200; ++px) {
            Color c = field_color(px / 199.0f);
            DrawRectangle(W - 220 + px, H - 36, 1, 18, c);
        }
        DrawText("0", W - 224, H - 36, 14, LIGHTGRAY);
        DrawText("B_max", W - 16, H - 36, 14, LIGHTGRAY);
 
        EndDrawing();
    }
 
    CloseWindow();
    return 0;
}