n_fluids.c

Go to the documentation of this file.

/*
 * Nilorea Library
 * Copyright (C) 2005-2026 Castagnier Mickael
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <https://www.gnu.org/licenses/>.
 */
 
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdbool.h>
#include <strings.h>
#include "nilorea/n_fluids.h"
#include "nilorea/n_common.h"
#include "nilorea/n_thread_pool.h"
 
#define N_FLUID_U_FIELD 0
#define N_FLUID_V_FIELD 1
#define N_FLUID_S_FIELD 2
 
void n_memset(void* dst, const void* val, size_t size, size_t count) {
    char* ptr = (char*)dst;
    while (count-- > 0) {
        memcpy(ptr, val, size);
        ptr += size;
    }
}
 
int destroy_n_fluid(N_FLUID** fluid) {
    __n_assert((*fluid), return FALSE);
 
    list_destroy(&(*fluid)->integrate_chunk_list);
    list_destroy(&(*fluid)->solveIncompressibility_chunk_list);
    list_destroy(&(*fluid)->advectVel_chunk_list);
    list_destroy(&(*fluid)->advectSmoke_chunk_list);
 
    FreeNoLog((*fluid)->u);
    FreeNoLog((*fluid)->newU);
    FreeNoLog((*fluid)->v);
    FreeNoLog((*fluid)->newV);
    FreeNoLog((*fluid)->p);
    FreeNoLog((*fluid)->s);
    FreeNoLog((*fluid)->m);
    FreeNoLog((*fluid)->newM);
    FreeNoLog((*fluid));
 
    return TRUE;
}
 
N_FLUID* new_n_fluid(double density, double gravity, size_t numIters, double dt, double overRelaxation, size_t sx, size_t sy) {
    N_FLUID* fluid = NULL;
 
    Malloc(fluid, N_FLUID, 1);
    __n_assert(fluid, return NULL);
 
    fluid->density = density;
    fluid->gravity = gravity;
    fluid->numIters = numIters;
    fluid->dt = dt;
    fluid->h = 1.0 / 100.0;
    fluid->overRelaxation = overRelaxation;
    fluid->numX = sx + 2;
    fluid->numY = sy + 2;
    fluid->numZ = 1; /* Z expansion not yet implemented */
    fluid->numCells = fluid->numX * fluid->numY * fluid->numZ;
 
    Malloc(fluid->u, double, fluid->numCells);
    if (!fluid->u) goto cleanup_fluid;
    Malloc(fluid->newU, double, fluid->numCells);
    if (!fluid->newU) goto cleanup_fluid;
    Malloc(fluid->v, double, fluid->numCells);
    if (!fluid->v) goto cleanup_fluid;
    Malloc(fluid->newV, double, fluid->numCells);
    if (!fluid->newV) goto cleanup_fluid;
    Malloc(fluid->p, double, fluid->numCells);
    if (!fluid->p) goto cleanup_fluid;
    Malloc(fluid->s, double, fluid->numCells);
    if (!fluid->s) goto cleanup_fluid;
    Malloc(fluid->m, double, fluid->numCells);
    if (!fluid->m) goto cleanup_fluid;
    Malloc(fluid->newM, double, fluid->numCells);
    if (!fluid->newM) goto cleanup_fluid;
 
    fluid->showSmoke = 1;
    fluid->showPressure = 0;
    fluid->showPaint = 0;
    fluid->fluid_production_percentage = 0.1;
    fluid->cScale = 16.0;
    // double precision. Taking 'value', if( fabs( value ) < float_tolerance )  value is considered as zero
    fluid->negative_float_tolerance = -0.00001;
    fluid->positive_float_tolerance = 0.00001;
 
    double d_val = 1.0;
    n_memset(fluid->m, &d_val, sizeof(d_val), fluid->numCells);
 
    long int nb_cores_l = get_nb_cpu_cores();
    if (nb_cores_l <= 0) {
        nb_cores_l = 1;
    }
    size_t nb_cores = (size_t)nb_cores_l;
 
    // precalculate and allocated N_PROC_PARAMS lists for threaded computing
    fluid->integrate_chunk_list = new_generic_list(nb_cores);
    fluid->solveIncompressibility_chunk_list = new_generic_list(nb_cores);
    fluid->advectVel_chunk_list = new_generic_list(nb_cores);
    fluid->advectSmoke_chunk_list = new_generic_list(nb_cores);
 
    size_t steps = (fluid->numX / nb_cores > 0) ? fluid->numX / nb_cores : 1;
 
    N_FLUID_THREAD_PARAMS* params = NULL;
    /* integrate */
    for (size_t i = 1; i < fluid->numX - 1; i += steps) {
        Malloc(params, N_FLUID_THREAD_PARAMS, 1);
        params->ptr = fluid;
        params->x_start = i;
        params->x_end = i + steps;
        params->y_start = 1;
        params->y_end = fluid->numY - 1;
        list_push(fluid->integrate_chunk_list, params, &free);
    }
    /* set the last batch to cover the full x range */
    if (fluid->integrate_chunk_list->end) {
        params = (N_FLUID_THREAD_PARAMS*)fluid->integrate_chunk_list->end->ptr;
        params->x_end = fluid->numX;
    }
 
    /* solveIncompressibility */
    for (size_t i = 1; i < fluid->numX - 1; i += steps) {
        Malloc(params, N_FLUID_THREAD_PARAMS, 1);
        params->ptr = fluid;
        params->x_start = i;
        params->x_end = i + steps;
        params->y_start = 1;
        params->y_end = fluid->numY - 1;
        list_push(fluid->solveIncompressibility_chunk_list, params, &free);
    }
    /* set the last batch to cover the full x range */
    if (fluid->solveIncompressibility_chunk_list->end) {
        params = (N_FLUID_THREAD_PARAMS*)fluid->solveIncompressibility_chunk_list->end->ptr;
        params->x_end = fluid->numX - 1;
    }
 
    /* advectVel */
    for (size_t i = 1; i < fluid->numX - 1; i += steps) {
        Malloc(params, N_FLUID_THREAD_PARAMS, 1);
        params->ptr = fluid;
        params->x_start = i;
        params->x_end = i + steps;
        params->y_start = 1;
        params->y_end = fluid->numY;
        list_push(fluid->advectVel_chunk_list, params, &free);
    }
    /* set the last batch to cover the full x range */
    if (fluid->advectVel_chunk_list->end) {
        params = (N_FLUID_THREAD_PARAMS*)fluid->advectVel_chunk_list->end->ptr;
        params->x_end = fluid->numX;
    }
 
    /* advectSmoke */
    for (size_t i = 1; i < fluid->numX - 1; i += steps) {
        Malloc(params, N_FLUID_THREAD_PARAMS, 1);
        params->ptr = fluid;
        params->x_start = i;
        params->x_end = i + steps;
        params->y_start = 1;
        params->y_end = fluid->numY - 1;
        list_push(fluid->advectSmoke_chunk_list, params, &free);
    }
    /* set the last batch to cover the full x range */
    if (fluid->advectSmoke_chunk_list->end) {
        params = (N_FLUID_THREAD_PARAMS*)fluid->advectSmoke_chunk_list->end->ptr;
        params->x_end = fluid->numX - 1;
    }
 
    return fluid;
 
cleanup_fluid:
    FreeNoLog(fluid->u);
    FreeNoLog(fluid->newU);
    FreeNoLog(fluid->v);
    FreeNoLog(fluid->newV);
    FreeNoLog(fluid->p);
    FreeNoLog(fluid->s);
    FreeNoLog(fluid->m);
    FreeNoLog(fluid->newM);
    FreeNoLog(fluid);
    return NULL;
} /* new_n_fluid */
 
void* n_fluid_integrate_proc(void* ptr) {
    N_FLUID_THREAD_PARAMS* params = (N_FLUID_THREAD_PARAMS*)ptr;
    N_FLUID* fluid = (N_FLUID*)params->ptr;
 
    size_t n = fluid->numY;
    for (size_t i = params->x_start; i < params->x_end; i++) {
        for (size_t j = params->y_start; j < params->y_end; j++) {
            if (!_z(fluid, s[i * n + j]) && !_z(fluid, s[i * n + j - 1]))
                fluid->v[i * n + j] += fluid->gravity * fluid->dt;
        }
    }
    return NULL;
}
 
int n_fluid_integrate(N_FLUID* fluid) {
    __n_assert(fluid, return FALSE);
 
    size_t n = fluid->numY;
    for (size_t i = 1; i < fluid->numX; i++) {
        for (size_t j = 1; j < fluid->numY - 1; j++) {
            if (!_z(fluid, s[i * n + j]) && !_z(fluid, s[i * n + j - 1]))
                fluid->v[i * n + j] += fluid->gravity * fluid->dt;
        }
    }
    return TRUE;
}
 
void* n_fluid_solveIncompressibility_proc(void* ptr) {
    N_FLUID_THREAD_PARAMS* params = (N_FLUID_THREAD_PARAMS*)ptr;
    N_FLUID* fluid = (N_FLUID*)params->ptr;
 
    double cp = (fluid->density * fluid->h) / fluid->dt;
 
    size_t n = fluid->numY;
    for (size_t i = params->x_start; i < params->x_end; i++) {
        for (size_t j = params->y_start; j < params->y_end; j++) {
            if (_z(fluid, s[i * n + j]))
                continue;
 
            double sx0 = fluid->s[(i - 1) * n + j];
            double sx1 = fluid->s[(i + 1) * n + j];
            double sy0 = fluid->s[i * n + j - 1];
            double sy1 = fluid->s[i * n + j + 1];
            double s = sx0 + sx1 + sy0 + sy1;
            if (_zd(fluid, s))
                continue;
 
            double div = fluid->u[(i + 1) * n + j] - fluid->u[i * n + j] + fluid->v[i * n + j + 1] - fluid->v[i * n + j];
            double p = (-div * fluid->overRelaxation) / s;
            fluid->p[i * n + j] += cp * p;
            fluid->u[i * n + j] -= sx0 * p;
            fluid->u[(i + 1) * n + j] += sx1 * p;
            fluid->v[i * n + j] -= sy0 * p;
            fluid->v[i * n + j + 1] += sy1 * p;
        }
    }
    return NULL;
}
 
int n_fluid_solveIncompressibility(N_FLUID* fluid) {
    __n_assert(fluid, return FALSE);
    size_t n = fluid->numY;
 
    double cp = (fluid->density * fluid->h) / fluid->dt;
 
    for (size_t iter = 0; iter < fluid->numIters; iter++) {
        for (size_t i = 1; i < fluid->numX - 1; i++) {
            for (size_t j = 1; j < fluid->numY - 1; j++) {
                if (_z(fluid, s[i * n + j]))
                    continue;
 
                double sx0 = fluid->s[(i - 1) * n + j];
                double sx1 = fluid->s[(i + 1) * n + j];
                double sy0 = fluid->s[i * n + j - 1];
                double sy1 = fluid->s[i * n + j + 1];
                double s = sx0 + sx1 + sy0 + sy1;
                if (_zd(fluid, s))
                    continue;
 
                double div = fluid->u[(i + 1) * n + j] - fluid->u[i * n + j] + fluid->v[i * n + j + 1] - fluid->v[i * n + j];
                double p = (-div * fluid->overRelaxation) / s;
                fluid->p[i * n + j] += cp * p;
                fluid->u[i * n + j] -= sx0 * p;
                fluid->u[(i + 1) * n + j] += sx1 * p;
                fluid->v[i * n + j] -= sy0 * p;
                fluid->v[i * n + j + 1] += sy1 * p;
            }
        }
    }
    return TRUE;
}
 
int n_fluid_extrapolate(N_FLUID* fluid) {
    __n_assert(fluid, return FALSE);
    size_t n = fluid->numY;
    for (size_t i = 0; i < fluid->numX; i++) {
        fluid->u[i * n + 0] = fluid->u[i * n + 1];
        fluid->u[i * n + fluid->numY - 1] = fluid->u[i * n + fluid->numY - 2];
    }
    for (size_t j = 0; j < fluid->numY; j++) {
        fluid->v[0 * n + j] = fluid->v[1 * n + j];
        fluid->v[(fluid->numX - 1) * n + j] = fluid->v[(fluid->numX - 2) * n + j];
    }
    return TRUE;
}
 
double n_fluid_sampleField(N_FLUID* fluid, double x, double y, uint32_t field) {
    __n_assert(fluid, return FALSE);
    size_t n = fluid->numY;
    double h1 = 1.0 / fluid->h;
    double h2 = 0.5 * fluid->h;
 
    x = MAX(MIN(x, (double)fluid->numX * fluid->h), fluid->h);
    y = MAX(MIN(y, (double)fluid->numY * fluid->h), fluid->h);
 
    double dx = 0.0;
    double dy = 0.0;
 
    const double* f = NULL;
    switch (field) {
        case N_FLUID_U_FIELD:
            f = fluid->u;
            dy = h2;
            break;
        case N_FLUID_V_FIELD:
            f = fluid->v;
            dx = h2;
            break;
        case N_FLUID_S_FIELD:
            f = fluid->m;
            dx = h2;
            dy = h2;
            break;
        default:
            f = NULL;
            return 0.0;
    }
 
    if (!f) return 0.0;
 
    double x0 = MIN(floor((x - dx) * h1), (double)(fluid->numX - 1));
    double tx = ((x - dx) - x0 * fluid->h) * h1;
    double x1 = MIN(x0 + 1, (double)(fluid->numX - 1));
 
    double y0 = MIN(floor((y - dy) * h1), (double)(fluid->numY - 1));
    double ty = ((y - dy) - y0 * fluid->h) * h1;
    double y1 = MIN(y0 + 1, (double)(fluid->numY - 1));
 
    double sx = 1.0 - tx;
    double sy = 1.0 - ty;
 
    double val = sx * sy * f[(size_t)(x0 * (double)n + y0)] +
                 tx * sy * f[(size_t)(x1 * (double)n + y0)] +
                 tx * ty * f[(size_t)(x1 * (double)n + y1)] +
                 sx * ty * f[(size_t)(x0 * (double)n + y1)];
 
    return val;
}
 
double n_fluid_avgU(N_FLUID* fluid, size_t i, size_t j) {
    __n_assert(fluid, return FALSE);
    if (j == 0) return 0.0;
    size_t n = fluid->numY;
    double u = (fluid->u[i * n + j - 1] + fluid->u[i * n + j] +
                fluid->u[(i + 1) * n + j - 1] + fluid->u[(i + 1) * n + j]) *
               0.25;
 
    return u;
}
 
double n_fluid_avgV(N_FLUID* fluid, size_t i, size_t j) {
    __n_assert(fluid, return FALSE);
    if (i == 0) return 0.0;
    size_t n = fluid->numY;
    double v = (fluid->v[(i - 1) * n + j] + fluid->v[i * n + j] +
                fluid->v[(i - 1) * n + j + 1] + fluid->v[i * n + j + 1]) *
               0.25;
    return v;
}
 
void* n_fluid_advectVel_proc(void* ptr) {
    N_FLUID_THREAD_PARAMS* params = (N_FLUID_THREAD_PARAMS*)ptr;
    N_FLUID* fluid = (N_FLUID*)params->ptr;
 
    size_t n = fluid->numY;
    double h2 = 0.5 * fluid->h;
    for (size_t i = params->x_start; i < params->x_end; i++) {
        for (size_t j = params->y_start; j < params->y_end; j++) {
            size_t index = i * n + j;
            // u component
            if (!_z(fluid, s[index]) && !_z(fluid, s[(i - 1) * n + j]) && j < fluid->numY - 1) {
                double x = (double)i * fluid->h;
                double y = (double)j * fluid->h + h2;
                double u = fluid->u[index];
                double v = n_fluid_avgV(fluid, i, j);
                // double v = n_fluid_sampleField( fluid , x , y , N_FLUID_V_FIELD );
                x = x - fluid->dt * u;
                y = y - fluid->dt * v;
                u = n_fluid_sampleField(fluid, x, y, N_FLUID_U_FIELD);
                fluid->newU[index] = u;
            }
            // v component
            if (!_z(fluid, s[index]) && !_z(fluid, s[index - 1]) && i < fluid->numX - 1) {
                double x = (double)i * fluid->h + h2;
                double y = (double)j * fluid->h;
                double u = n_fluid_avgU(fluid, i, j);
                // double u = n_fluid_sampleField( fluid , x , y , N_FLUID_U_FIELD );
                double v = fluid->v[index];
                x = x - fluid->dt * u;
                y = y - fluid->dt * v;
                v = n_fluid_sampleField(fluid, x, y, N_FLUID_V_FIELD);
                fluid->newV[index] = v;
            }
        }
    }
    return NULL;
}
 
int n_fluid_advectVel(N_FLUID* fluid) {
    __n_assert(fluid, return FALSE);
 
    memcpy(fluid->newU, fluid->u, fluid->numCells * sizeof(double));
    memcpy(fluid->newV, fluid->v, fluid->numCells * sizeof(double));
 
    size_t n = fluid->numY;
    double h2 = 0.5 * fluid->h;
    for (size_t i = 1; i < fluid->numX; i++) {
        for (size_t j = 1; j < fluid->numY; j++) {
            size_t index = i * n + j;
            // u component
            if (!_z(fluid, s[index]) && !_z(fluid, s[(i - 1) * n + j]) && j < fluid->numY - 1) {
                double x = (double)i * fluid->h;
                double y = (double)j * fluid->h + h2;
                double u = fluid->u[index];
                double v = n_fluid_avgV(fluid, i, j);
                // double v = n_fluid_sampleField( fluid , x , y , N_FLUID_V_FIELD );
                x = x - fluid->dt * u;
                y = y - fluid->dt * v;
                u = n_fluid_sampleField(fluid, x, y, N_FLUID_U_FIELD);
                fluid->newU[index] = u;
            }
            // v component
            if (!_z(fluid, s[index]) && !_z(fluid, s[index - 1]) && i < fluid->numX - 1) {
                double x = (double)i * fluid->h + h2;
                double y = (double)j * fluid->h;
                double u = n_fluid_avgU(fluid, i, j);
                // double u = n_fluid_sampleField( fluid , x , y , N_FLUID_U_FIELD );
                double v = fluid->v[index];
                x = x - fluid->dt * u;
                y = y - fluid->dt * v;
                v = n_fluid_sampleField(fluid, x, y, N_FLUID_V_FIELD);
                fluid->newV[index] = v;
            }
        }
    }
    double* ptr = fluid->u;
    fluid->u = fluid->newU;
    fluid->newU = ptr;
 
    ptr = fluid->v;
    fluid->v = fluid->newV;
    fluid->newV = ptr;
 
    return TRUE;
}
 
void* n_fluid_advectSmoke_proc(void* ptr) {
    N_FLUID_THREAD_PARAMS* params = (N_FLUID_THREAD_PARAMS*)ptr;
    N_FLUID* fluid = (N_FLUID*)params->ptr;
 
    size_t n = fluid->numY;
    double h2 = 0.5 * fluid->h;
    for (size_t i = params->x_start; i < params->x_end; i++) {
        for (size_t j = params->y_start; j < params->y_end; j++) {
            size_t index = i * n + j;
            if (!_z(fluid, s[index])) {
                double u = (fluid->u[index] + fluid->u[(i + 1) * n + j]) * 0.5;
                double v = (fluid->v[index] + fluid->v[index + 1]) * 0.5;
                double x = (double)i * fluid->h + h2 - fluid->dt * u;
                double y = (double)j * fluid->h + h2 - fluid->dt * v;
 
                fluid->newM[index] = n_fluid_sampleField(fluid, x, y, N_FLUID_S_FIELD);
            }
        }
    }
    return NULL;
}
 
int n_fluid_advectSmoke(N_FLUID* fluid) {
    __n_assert(fluid, return FALSE);
 
    size_t n = fluid->numY;
    double h2 = 0.5 * fluid->h;
 
    memcpy(fluid->newM, fluid->m, fluid->numCells * sizeof(double));
 
    for (size_t i = 1; i < fluid->numX - 1; i++) {
        for (size_t j = 1; j < fluid->numY - 1; j++) {
            size_t index = i * n + j;
            if (!_z(fluid, s[index])) {
                double u = (fluid->u[index] + fluid->u[(i + 1) * n + j]) * 0.5;
                double v = (fluid->v[index] + fluid->v[index + 1]) * 0.5;
                double x = (double)i * fluid->h + h2 - fluid->dt * u;
                double y = (double)j * fluid->h + h2 - fluid->dt * v;
 
                fluid->newM[index] = n_fluid_sampleField(fluid, x, y, N_FLUID_S_FIELD);
            }
        }
    }
    double* ptr = fluid->m;
    fluid->m = fluid->newM;
    fluid->newM = ptr;
 
    return TRUE;
}
 
int n_fluid_simulate(N_FLUID* fluid) {
    __n_assert(fluid, return FALSE);
    n_fluid_integrate(fluid);
    memset(fluid->p, 0, fluid->numCells * sizeof(double));
    n_fluid_solveIncompressibility(fluid);
    n_fluid_extrapolate(fluid);
    n_fluid_advectVel(fluid);
    n_fluid_advectSmoke(fluid);
    return TRUE;
}
 
int n_fluid_simulate_threaded(N_FLUID* fluid, THREAD_POOL* thread_pool) {
    __n_assert(fluid, return FALSE);
 
    // n_fluid_integrate( fluid );
    list_foreach(node, fluid->integrate_chunk_list) {
        add_threaded_process(thread_pool, &n_fluid_integrate_proc, (void*)node->ptr, SYNCED_PROC);
    }
    start_threaded_pool(thread_pool);
    wait_for_synced_threaded_pool(thread_pool);
    refresh_thread_pool(thread_pool);
 
    // set pressure to 0
    memset(fluid->p, 0, fluid->numCells * sizeof(double));
 
    // n_fluid_solveIncompressibility( fluid );
    for (size_t iter = 0; iter < fluid->numIters; iter++) {
        list_foreach(node, fluid->solveIncompressibility_chunk_list) {
            add_threaded_process(thread_pool, &n_fluid_solveIncompressibility_proc, (void*)node->ptr, SYNCED_PROC);
        }
        start_threaded_pool(thread_pool);
        wait_for_synced_threaded_pool(thread_pool);
        refresh_thread_pool(thread_pool);
    }
 
    // extrapolate
    n_fluid_extrapolate(fluid);
 
    // n_fluid_advectVel( fluid );
    memcpy(fluid->newU, fluid->u, fluid->numCells * sizeof(double));
    memcpy(fluid->newV, fluid->v, fluid->numCells * sizeof(double));
 
    list_foreach(node, fluid->advectVel_chunk_list) {
        add_threaded_process(thread_pool, &n_fluid_advectVel_proc, (void*)node->ptr, SYNCED_PROC);
    }
    start_threaded_pool(thread_pool);
    wait_for_synced_threaded_pool(thread_pool);
    refresh_thread_pool(thread_pool);
 
    double* ptr = fluid->u;
    fluid->u = fluid->newU;
    fluid->newU = ptr;
 
    ptr = fluid->v;
    fluid->v = fluid->newV;
    fluid->newV = ptr;
 
    // n_fluid_advectSmoke( fluid );
    memcpy(fluid->newM, fluid->m, fluid->numCells * sizeof(double));
    list_foreach(node, fluid->advectSmoke_chunk_list) {
        add_threaded_process(thread_pool, &n_fluid_advectSmoke_proc, (void*)node->ptr, SYNCED_PROC);
    }
    start_threaded_pool(thread_pool);
    wait_for_synced_threaded_pool(thread_pool);
    refresh_thread_pool(thread_pool);
 
    ptr = fluid->m;
    fluid->m = fluid->newM;
    fluid->newM = ptr;
 
    return TRUE;
}
 
int n_fluid_setObstacle(N_FLUID* fluid, double x, double y, double vx, double vy, double r) {
    __n_assert(fluid, return FALSE);
 
    size_t n = fluid->numY;
    for (size_t i = 1; i < fluid->numX - 2; i++) {
        for (size_t j = 1; j < fluid->numY - 2; j++) {
            double dx = ((double)i + 0.5) - x;
            double dy = ((double)j + 0.5) - y;
 
            if (i > 7 && (dx * dx + dy * dy < r * r)) {
                fluid->s[i * n + j] = 0.0;
                fluid->m[i * n + j] = 1.0;
                fluid->u[i * n + j] = vx;
                fluid->u[(i + 1) * n + j] = vx;
                fluid->v[i * n + j] = vy;
                fluid->v[i * n + j + 1] = vy;
            }
        }
    }
    return TRUE;
}
 
int n_fluid_resetObstacles(N_FLUID* fluid) {
    __n_assert(fluid, return FALSE);
 
    size_t n = fluid->numY;
    for (size_t i = 1; i < fluid->numX - 2; i++) {
        for (size_t j = 1; j < fluid->numY - 2; j++) {
            fluid->s[i * n + j] = 1.0;
        }
    }
 
    return TRUE;
}
 
int n_fluid_setObstacleFromBitmap(N_FLUID* fluid, ALLEGRO_BITMAP* bitmap, double x, double y, double vx, double vy, double r) {
    __n_assert(fluid, return FALSE);
 
    al_lock_bitmap(bitmap, al_get_bitmap_format(bitmap), ALLEGRO_LOCK_READONLY);
 
    size_t n = fluid->numY;
    for (size_t i = 1; i < fluid->numX - 2; i++) {
        for (size_t j = 1; j < fluid->numY - 2; j++) {
            double dx = ((double)i + 0.5) - x;
            double dy = ((double)j + 0.5) - y;
 
            if (i > 7 && (dx * dx + dy * dy < r * r)) {
                fluid->s[i * n + j] = 0.0;
                fluid->m[i * n + j] = 1.0;
                fluid->u[i * n + j] = vx;
                fluid->v[i * n + j] = vy;
 
                fluid->u[(i + 1) * n + j] = vx;
                fluid->v[i * n + j + 1] = vy;
            }
        }
    }
 
    al_unlock_bitmap(bitmap);
 
    return TRUE;
}
 
ALLEGRO_COLOR n_fluid_getSciColor(const N_FLUID* fluid, double val, double minVal, double maxVal) {
    val = MIN(MAX(val, minVal), maxVal - 0.0001);
    double d = maxVal - minVal;
    if (_zd(fluid, d)) {
        val = 0.5;
    } else {
        val = (val - minVal) / d;
    }
    double m = 0.25;
    size_t num = (size_t)floor(val / m);
    double s = (val - (double)num * m) / m;
    double r = 0.0, g = 0.0, b = 0.0;
    switch (num) {
        case 0:
            r = 0.0;
            g = s;
            b = 1.0;
            break;
        case 1:
            r = 0.0;
            g = 1.0;
            b = 1.0 - s;
            break;
        case 2:
            r = s;
            g = 1.0;
            b = 0.0;
            break;
        case 3:
            r = 1.0;
            g = 1.0 - s;
            b = 0.0;
            break;
    }
    // return[255*r,255*g,255*b, 255]
    return al_map_rgb_f((float)r, (float)g, (float)b);
}
 
int n_fluid_draw(N_FLUID* fluid) {
    __n_assert(fluid, return FALSE);
 
    size_t n = fluid->numY;
 
    double minP = fluid->p[0];
    double maxP = fluid->p[0];
 
    if (fluid->showPressure) {
        for (size_t i = 0; i < fluid->numCells; i++) {
            minP = MIN(minP, fluid->p[i]);
            maxP = MAX(maxP, fluid->p[i]);
        }
    }
 
    ALLEGRO_COLOR color;
    double cScale = fluid->cScale;
 
    for (size_t i = 0; i < fluid->numX; i++) {
        for (size_t j = 0; j < fluid->numY; j++) {
            double xd = (double)i * cScale;
            double yd = (double)j * cScale;
            double cxd = xd + cScale;
            double cyd = yd + cScale;
 
            double s = fluid->m[i * n + j];
 
            if (fluid->showPaint) {
                color = n_fluid_getSciColor(fluid, s, 0.0, 1.0);
            } else if (fluid->showPressure) {
                double p = fluid->p[i * n + j];
                color = n_fluid_getSciColor(fluid, p, minP, maxP);
                if (fluid->showSmoke) {
                    float color_vec_f[3] = {0.0, 0.0, 0.0};
                    al_unmap_rgb_f(color, color_vec_f, color_vec_f + 1, color_vec_f + 2);
                    color = al_map_rgb_f((float)MAX(0.0, color_vec_f[0] - s), (float)MAX(0.0, color_vec_f[1] - s), (float)MAX(0.0, color_vec_f[2] - s));
                }
            } else if (fluid->showSmoke) {
                color = al_map_rgb_f((float)(1.0 - s), 0.0f, 0.0f);
            } else {
                color = al_map_rgb_f((float)s, (float)s, (float)s);
            }
            al_draw_filled_rectangle((float)xd, (float)yd, (float)cxd, (float)cyd, color);
        }
    }
    return TRUE;
}