《Fluid Engine Development》 學習筆記2-基礎

斷斷續續花了一個月，終於把這本書的一二兩章啃了下來，理解流體模擬的理論似乎不難，無論是《Fluid Simulation for Computer Graphics》還是《計算流體力學基礎及其應用》都能很好幫助程式設計師去理解這些原理，可在缺乏實踐情況下，這種對原理的理解其實跟死記硬背沒什麼區別。《Fluid Engine Development》提供了一個實現完成的流體模擬引擎以及它的程式設計實現原理，充分幫助程式設計師通過程式設計實現流體動畫引擎，以此完成流體模擬學習的第一步。這不，早在今年一月就嚷嚷研究學習流體模擬卻苦苦掙扎無法入門的我，在抄著作者的程式碼看著作者的書的情況，終於實現了一個流體模擬引擎。我終於可以自信地說自己已經入門了流體模擬o(╥﹏╥)o。

流體模擬

我學習流體模擬的本質原因，是因為迷戀上了複雜的流體動畫，渴望通過程式設計來創造它們。流體模擬主要是指結合流體模擬的物理現象、方程和計算機圖形學的方法來模擬海面、海浪、煙霧等場景。

流體模擬的基本方法可分為三類: 基於紋理變換的流體模擬、基於二維高度場網格的流體模擬以及基於真實物理方程的流體模擬。

基於紋理變換的流體模擬只需要對水面紋理進行法向擾動後、繪製水面的倒影(反射)以及繪製水底的情況(透射)即可繪製出一般的水面效果。 但這種方法由於其根本上沒有水面網格，所以水面起伏的繪製效果不明顯。

基於二維高度場的網格流體模擬方法把水面表示成為一個連續的平面網格，然後生成一系列對應於這張網路的連續的高度紋理-稱為高度圖。接著每個網格頂點對應於一個高度圖的畫素，作為水面高度，從而表示出整個水面。

以上2者方法相對簡單，但真實度不高，本質上只是hack，基於真實物理方程的流體模擬才能製作目前所能呈現在計算機上最真實的流體動畫。儘管這種方法，受限於過高的計算量，一直處於離線渲染領域，很難在實時遊戲中使用，但隨著實時流體模擬技術的進步以及硬體條件的進步，在遊戲中使用這種方法並不遙遠，實際上，這甚至已經成為了現實。看看這個視訊，實時流體模擬已經能做得很好了。

流體模擬的基本方法主要有三種，一種使用粒子（拉格朗日方法），一種使用網格（尤拉方法），第三種則是2者的混合。本文主要講基於粒子的流體模擬。

基於物理的動畫

流體動畫首先是動畫，動畫的本質是在給定時間序列播放一系列影象，由此我們可以編寫幀的結構體和動畫的抽象類。我們建立新的動畫時，只要編寫類繼承Animation，再重寫onUpdate函式，在指定幀更新影象即可。

struct Frame final {
    public:
        int index = 0;
        double timeIntervalInSeconds = 1.0 / 60.0;
        double GetTimeInSeconds() const;
        void Advance();
        void Advance(int delta);
};
    
class Animation{
public:
    void Update(const Frame& frame);
protected:

    //**********************************************
    //the function is called from Animation:Update();
    //this function and implement its logic for updating the animation state.
    //**********************************************
    virtual void onUpdate(const Frame& frame) = 0;
};

void Animation::Update(const Frame & frame){
    onUpdate(frame);
}

double Frame::GetTimeInSeconds() const{
    return index * timeIntervalInSeconds;
}

void Frame::Advance(){
    ++index;
}

void Frame::Advance(int delta){
    index += delta;
}

流體動畫屬於基於物理的動畫，它不同於正常的動畫，正常的動畫並不依賴外界的輸入以及其他幀的資料和狀態，它們大多是根據時間狀態的改變來播放不同的影象，或是通過以時間為輸入的函式（例如sin(t)影象）來切換狀態。基於物理的動畫正好相反，當前幀的資料和狀態完全是由上一幀決定的。由此編寫了PhysicsAnimation類

class PhysicsAnimation : public Animation{
    public:
        double GetCurrentTime() const { return _currentTime; }
    private:
        void initialize();
        void timeStep(double timeIntervalInSeconds);

        Frame _currentFrame;
        double _currentTime = 0.0;
protected:
    virtual void onUpdate(const Frame& frame) override final;

    //**********************************************
    //the function is called from Animation:timeStep(double);
    //This function is called for each time-step;
    //**********************************************
    virtual void onTimeStep(double timeIntervalInSeconds) = 0;

    //**********************************************
    //the function is called from Animation:initialize();
    //Called at frame 0 to initialize the physics state.;
    //**********************************************
    virtual void onInitialize() = 0;
};

void PhysicsAnimation::initialize(){
    onInitialize();
}


void PhysicsAnimation::onUpdate(const Frame & frame){
    if (frame.index > _currentFrame.index) {
        if (_currentFrame.index < 0) {
            initialize();
        }

        int numberOfFrames = frame.index - _currentFrame.index;

        for (int i = 0; i < numberOfFrames; ++i) {
            timeStep(frame.timeIntervalInSeconds);
        }

        _currentFrame = frame;
    }
}

我們可以看到PhysicsAnimation重寫了onUpdate函式，採用漸近更新每一幀

粒子動畫系統

我們使用粒子模擬流體，那就不可避免的要使用粒子系統

以下是我個人流體引擎粒子系統動畫系統的部分標頭檔案

class ParticleSystemData3
{
    public:
    typedef std::vector<double> ScalarData;
    typedef std::vector<Vector3D> VectorData;

    ParticleSystemData3();
    virtual ~ParticleSystemData3();
    explicit ParticleSystemData3(size_t numberOfParticles);
    void Resize(size_t newNumberOfParticles);
    size_t GetNumberOfParticles() const;
    size_t AddVectorData(const Vector3D& initialVal = Vector3D::zero);
    size_t AddScalarData(double initialVal = 0.0);
    const std::vector<Vector3D>& GetPositions() const;
    std::vector<Vector3D>& GetPositions();
    const std::vector<Vector3D>& GetVelocities() const;
    std::vector<Vector3D>& GetVelocities();
    const std::vector<Vector3D>& GetForces() const;
    std::vector<Vector3D>& GetForces();
    void AddParticle(
        const Vector3D& newPosition,
        const Vector3D& newVelocity,
        const Vector3D& newForce = Vector3D::zero);
    void AddParticles(
        const std::vector<Vector3D>& newPositions,
        const std::vector<Vector3D>& newVelocities,
        const std::vector<Vector3D>& newForces);
    const std::vector<double>& ScalarDataAt(size_t idx) const;
    std::vector<double>& ScalarDataAt(size_t idx);
    const std::vector<Vector3D>& VectorDataAt(size_t idx) const;
    std::vector<Vector3D>& VectorDataAt(size_t idx);

    double GetParticleRadius() const;
    void SetParticleRadius(double newRadius);
    double GetParticleMass() const;
    virtual void SetParticleMass(double newMass);

    void BuildNeighborSearcher(double maxSearchRadius);
    void BuildNeighborLists(double maxSearchRadius);

    const std::shared_ptr<PointNeighborSearcher3> & GetNeighborSearcher() const;
    const std::vector<std::vector<size_t>>& GetNeighborLists() const;
    protected:
    size_t _positionIdx;
    size_t _velocityIdx;
    size_t _forceIdx;
    size_t _numberOfParticles = 0;
    double _radius = 1e-3;
    double _mass = 1e-3;

    std::vector<ScalarData> _scalarDataList;
    std::vector<VectorData> _vectorDataList;
    std::shared_ptr<PointNeighborSearcher3> _neighborSearcher;
    std::vector<std::vector<size_t>>_neighborLists;
};

class ParticleSystemSolver3 : public PhysicsAnimation{
private:
    void timeStepStart(double timeStepInSeconds);
    void timeStepEnd(double timeStepInSeconds);
    void timeIntegration(double timeIntervalInSeconds);
    void resolveCollision();
    void updateCollider(double timeStepInSeconds);
    void updateEmitter(double timeStepInSeconds);
    ParticleSystemData3::VectorData _newPositions;
    ParticleSystemData3::VectorData _newVelocities;
protected:
    void onTimeStep(double timeIntervalInSeconds) override;
    virtual void onInitialize() override;
    //**********************************************
    //the function is called in ParticleSystemSolver3:timeStepStart(double);
    // Called when a time-step is about to begin;
    //**********************************************
    virtual void onTimeStepStart(double timeStepInSeconds);

    //**********************************************
    //the function is called in ParticleSystemSolver3:timeStepEnd(double);
    // Called when a time-step is about to end;
    //**********************************************
    virtual void onTimeStepEnd(double timeStepInSeconds);

    //**********************************************
    //the function is called in ParticleSystemSolver3:onTimeStep(double);
    //accumulate forces
    //**********************************************
    virtual void accumulateForces(double timeIntervalInSeconds);

    std::shared_ptr<ParticleSystemData3> _particleSystemData;
    std::shared_ptr<VectorField3> _wind;
    Vector3D _gravity = Vector3D(0.0, -9.8, 0.0);
    double _dragCoefficient = 1e-4;
    std::shared_ptr<Collider3> _collider;
    double _restitutionCoefficient = 0.0;
    std::shared_ptr<ParticleEmitter3> _emitter;
};

主要看粒子系統重寫的onTimeStep函式，它用以更新粒子系統，它依次呼叫了5個函式。

accumulateForces用以累加作用於粒子上的力，timeIntegration用以更新粒子速度和位置.resolveCollision處理碰撞，timeStepStart，timeStepEnd分別用作每幀邏輯更新前後的事件呼叫以及內部資料更新

void ParticleSystemSolver3::onTimeStep(double timeIntervalInSeconds)
{
    timeStepStart(timeIntervalInSeconds);

    accumulateForces(timeIntervalInSeconds);
    timeIntegration(timeIntervalInSeconds);
    resolveCollision();

    timeStepEnd(timeIntervalInSeconds);
}

具體實現如下

void CalfFluidEngine::ParticleSystemSolver3::timeStepStart(double timeStepInSeconds)
{
    auto& forces = _particleSystemData->GetForces();
    tbb::parallel_for(
        tbb::blocked_range<size_t>(0, forces.size()),
        [&](const tbb::blocked_range<size_t> & b) {
        for (size_t i = b.begin(); i != b.end(); ++i)
            forces[i] = Vector3D::zero;
    });

    updateCollider(timeStepInSeconds);
    updateEmitter(timeStepInSeconds);

    size_t n = _particleSystemData->GetNumberOfParticles();
    _newPositions.resize(n);
    _newVelocities.resize(n);

    onTimeStepStart(timeStepInSeconds);
}

void CalfFluidEngine::ParticleSystemSolver3::accumulateForces(double timeIntervalInSeconds)
{
    size_t n = _particleSystemData->GetNumberOfParticles();
    auto& forces = _particleSystemData->GetForces();
    auto& velocities = _particleSystemData->GetVelocities();
    auto& positions = _particleSystemData->GetPositions();
    const double mass = _particleSystemData->GetParticleMass();

    tbb::parallel_for(
        tbb::blocked_range<size_t>(0, n),
        [&](const tbb::blocked_range<size_t> & b) {
            for (size_t i = b.begin(); i != b.end(); ++i){
                // Gravity
                Vector3D force = mass * _gravity;

                // Wind forces
                Vector3D relativeVel = velocities[i] - _wind->Sample(positions[i]);
                force += -_dragCoefficient * relativeVel;

                forces[i] += force;
        }
        
    });
}

void CalfFluidEngine::ParticleSystemSolver3::timeIntegration(double timeIntervalInSeconds)
{
    size_t n = _particleSystemData->GetNumberOfParticles();
    auto& forces = _particleSystemData->GetForces();
    auto& velocities = _particleSystemData->GetVelocities();
    auto& positions = _particleSystemData->GetPositions();
    const double mass = _particleSystemData->GetParticleMass();

    tbb::parallel_for(
        tbb::blocked_range<size_t>(0, n),
        [&](const tbb::blocked_range<size_t> & b) {
        for (size_t i = b.begin(); i != b.end(); ++i) {
            Vector3D& newVelocity = _newVelocities[i];
            newVelocity = velocities[i];
            newVelocity.AddScaledVector(forces[i] / mass, timeIntervalInSeconds);

            Vector3D& newPosition = _newPositions[i];
            newPosition = positions[i];
            newPosition.AddScaledVector(newVelocity, timeIntervalInSeconds);
        }

    });
}

void CalfFluidEngine::ParticleSystemSolver3::resolveCollision()
{
    if (_collider != nullptr) {
        size_t numberOfParticles = _particleSystemData->GetNumberOfParticles();
        const double radius = _particleSystemData->GetParticleRadius();

        tbb::parallel_for(
            tbb::blocked_range<size_t>(size_t(0), numberOfParticles),
            [&](const tbb::blocked_range<size_t> & b) {
            for (size_t i = b.begin(); i != b.end(); ++i)
            {
                _collider->ResolveCollision(
                    radius,
                    _restitutionCoefficient,
                    &_newPositions[i],
                    &_newVelocities[i]);
            }
                
        });
    }
}

void CalfFluidEngine::ParticleSystemSolver3::timeStepEnd(double timeStepInSeconds)
{
    size_t n = _particleSystemData->GetNumberOfParticles();
    auto& positions = _particleSystemData->GetPositions();
    auto& velocities = _particleSystemData->GetVelocities();

    tbb::parallel_for(
        tbb::blocked_range<size_t>(0, n),
        [&](const tbb::blocked_range<size_t> & b) {
        for (size_t i = b.begin(); i != b.end(); ++i)
        {
            positions[i] = _newPositions[i];
            velocities[i] = _newVelocities[i];
        }
            
    });

    onTimeStepEnd(timeStepInSeconds);
}

碰撞檢測和處理

上一節我們提到了碰撞，這個需要特別列出一節來講。

碰撞在物理模擬乃至物理引擎領域都是一個非常複雜的問題，《Fluid Engine Development》對於這一塊並沒有進行深入，碰撞檢測只是簡單檢測碰撞點是否在幾何表面的內部，以及碰撞點，粒子的中心間的距離是否小於粒子的半徑。

至於碰撞的處理，主要是改變粒子的速度，和把粒子修正到正確的位置。

粒子的正確位置為碰撞點+碰撞點法線 * 粒子半徑,也就是緊貼碰撞面表面的位置

至於粒子的速度，這裡把粒子相對碰撞面的速度拆解切向速度和法向速度

經過碰撞，法向速度變更為 \(- R V_{n}\) （R是恢復係數，Vn是法向速度）

切向速度變更為 \(max\left( 1 - u \frac{\left|\Delta V_{n}\right|}{V_{t}} \right) V_{t}\)

切向速度的變更大家或許會有疑惑，這裡列出推導過程

假定\(\Delta V_{n} = V_{n}^{new} - V_{n} = (- R - 1)V_{n}\)

已知\(\Delta V_{t} = a_{t}\Delta t = F_{f}/m \Delta t = u F_{n}/m \Delta t\)

又上類推得\(\Delta V_{n} = a_{n}\Delta t = F_{n}/m \Delta t\)

所以\(\Delta V_{t} = u \Delta V_{n}\)

因為新切向速度必然小於原速度，所以\(\Delta V_{t} = min(u\Delta V_{n},V_{t})\)

可得 \(V_{t}^{new} =max\left( 1 - u \frac{\left|\Delta V_{n}\right|}{V_{t}} \right) V_{t}\)

具體程式碼實現如下

struct ColliderQueryResult final {
            double distance;
            Vector3D point;
            Vector3D normal;
            Vector3D velocity;
};

void CalfFluidEngine::Collider3::ResolveCollision(
    double radius, 
    double restitutionCoefficient, 
    Vector3D * position, 
    Vector3D * velocity)
{
    ColliderQueryResult res;
    getClosestPoint(_surface, *position, &res);

    if (isPenetrating(res, *position, radius))
    {
        Vector3D targetNormal = res.normal;
        Vector3D targetPoint = res.point + radius * targetNormal;
        Vector3D colliderVelAtTargetPoint = res.velocity;

        Vector3D relativeVel = *velocity - colliderVelAtTargetPoint;
        double normalDotRelativeVel = Vector3D::Dot(targetNormal, relativeVel);
        Vector3D relativeVelN = normalDotRelativeVel * targetNormal;
        Vector3D relativeVelT = relativeVel - relativeVelN;

        if (normalDotRelativeVel < 0.0) {
            Vector3D deltaRelativeVelN =
                (-restitutionCoefficient - 1.0) * relativeVelN;
            relativeVelN *= -restitutionCoefficient;

            // Apply friction to the tangential component of the velocity
            // From Bridson et al., Robust Treatment of Collisions, Contact and
            // Friction for Cloth Animation, 2002
            // http://graphics.stanford.edu/papers/cloth-sig02/cloth.pdf
            if (relativeVelT.SquareMagnitude() > 0.0) {
                double frictionScale = std::max(
                    1.0 - _frictionCoeffient * deltaRelativeVelN.Magnitude() /
                    relativeVelT.Magnitude(),
                    0.0);
                relativeVelT *= frictionScale;
            }

            *velocity =
                relativeVelN + relativeVelT + colliderVelAtTargetPoint;
        }

        *position = targetPoint;
    }
}

bool CalfFluidEngine::Collider3::isPenetrating(
    const ColliderQueryResult & colliderPoint, 
    const Vector3D & position, 
    double radius)
{
    return _surface->IsInside(position) ||
        colliderPoint.distance < radius;
}

粒子發射器

粒子系統有一個名為updateEmitter的函式需要注意一下，它在timeStepStart函式和onInitialize()函式中被呼叫，用於隨機生成粒子系統的粒子資料。不過在例項demo中，粒子只生成了一次，畢竟每幀重新生成粒子是非常耗時的。粒子發射器的核心程式碼如下。

void CalfFluidEngine::VolumeParticleEmitter3::onUpdate(double currentTimeInSeconds, double timeIntervalInSeconds)
{
    auto particles = GetTarget();

    if (particles == nullptr) {
        return;
    }

    if (_numberOfEmittedParticles > 0 && _isOneShot) {
        return;
    }

    std::vector<Vector3D> newPositions;
    std::vector<Vector3D> newVelocities;
    std::vector<Vector3D> Forces;
    emit(particles, &newPositions, &newVelocities);

    particles->AddParticles(newPositions, newVelocities, Forces);
}

void CalfFluidEngine::VolumeParticleEmitter3::emit(const std::shared_ptr<ParticleSystemData3>& particles, std::vector<Vector3D>* newPositions, std::vector<Vector3D>* newVelocities)
{
    if (_surface == nullptr) return;
    _surface->Update();

    BoundingBox3D region = _bounds;
    if (_surface->IsBounded()) {
        BoundingBox3D surfaceBox = _surface->GetBoundingBox();
        region.lowerCorner = max(region.lowerCorner, surfaceBox.lowerCorner);
        region.upperCorner = min(region.upperCorner, surfaceBox.upperCorner);
    }

    const double j = GetJitter();
    const double maxJitterDist = 0.5 * j * _spacing;

    if (_allowOverlapping || _isOneShot) 
    {
        _pointGenerator->ForEachPoint(region, _spacing, [&](const Vector3D& point) {
            Vector3D randomDir = uniformSampleSphere(random(_rng), random(_rng));
            Vector3D offset = maxJitterDist * randomDir;
            Vector3D candidate = point + offset;
            if (_surface->SignedDistance(candidate) <= 0.0) {
                if (_numberOfEmittedParticles < _maxNumberOfParticles) {
                    newPositions->push_back(candidate);
                    ++_numberOfEmittedParticles;
                }
                else {
                    return false;
                }
            }

            return true;
        });
    }
    else 
    {
        PointHashGridSearcher3 neighborSearcher(
            Vector3<size_t>(
                kDefaultHashGridResolution, 
                kDefaultHashGridResolution,
                kDefaultHashGridResolution),
            2.0 * _spacing);
        if (!_allowOverlapping) {
            neighborSearcher.Build(particles->GetPositions());
        }

        _pointGenerator->ForEachPoint(region, _spacing, [&](const Vector3D& point) {
            Vector3D randomDir = uniformSampleSphere(random(_rng), random(_rng));
            Vector3D offset = maxJitterDist * randomDir;
            Vector3D candidate = point + offset;
            if (_surface->SignedDistance(candidate) <= 0.0 &&
                (!_allowOverlapping &&!neighborSearcher.HasNearbyPoint(candidate, _spacing))) {
                if (_numberOfEmittedParticles < _maxNumberOfParticles) {
                    newPositions->push_back(candidate);
                    neighborSearcher.Add(candidate);
                    ++_numberOfEmittedParticles;
                }
                else {
                    return false;
                }
            }

            return true;
        });
    }

    newVelocities->resize(newPositions->size());

    tbb::parallel_for(
        tbb::blocked_range<size_t>(0, newVelocities->size()),
        [&](const tbb::blocked_range<size_t> & b) {
        for (size_t i = b.begin(); i != b.end(); ++i)
            (*newVelocities)[i] = velocityAt((*newPositions)[i]);
    });
}

Vector3D CalfFluidEngine::VolumeParticleEmitter3::velocityAt(const Vector3D & point) const
{
    Vector3D r = point - _surface->transform.GetTranslation();
    return _linearVel + Vector3D::Cross(_angularVel,r) + _initialVel;
}

void CalfFluidEngine::BccLatticePointGenerator::ForEachPoint(const BoundingBox3D & boundingBox, double spacing, const std::function<bool(const Vector3D&)>& callback) const
{
    double halfSpacing = spacing / 2.0;
    double boxWidth = boundingBox.GetWidth();
    double boxHeight = boundingBox.GetHeight();
    double boxDepth = boundingBox.GetDepth();

    Vector3D position;
    bool hasOffset = false;
    bool shouldQuit = false;
    for (int k = 0; k * halfSpacing <= boxDepth && !shouldQuit; ++k) {
        position.z = k * halfSpacing + boundingBox.lowerCorner.z;

        double offset = (hasOffset) ? halfSpacing : 0.0;

        for (int j = 0; j * spacing + offset <= boxHeight && !shouldQuit; ++j) {
            position.y = j * spacing + offset + boundingBox.lowerCorner.y;

            for (int i = 0; i * spacing + offset <= boxWidth; ++i) {
                position.x = i * spacing + offset + boundingBox.lowerCorner.x;
                if (!callback(position)) {
                    shouldQuit = true;
                    break;
                }
            }
        }

        hasOffset = !hasOffset;
    }
}

鄰域粒子搜尋

我們看到上節中粒子發射器的程式碼實現裡使用了PointHashGridSearcher3這個類，這是我為什麼說粒子發生器需要注意的原因，這是個非常精妙的演算法，它將3D點對映到整數，它被這樣設計的意義是為了加速搜尋。

它將長方體的包圍盒區域劃分成多個長方形塊，將3D點的座標除以長方形塊的邊長，將它轉換成三維整形數（Vector3），再將三維整數轉換成普通整形數作為3D座標的雜湊值。

程式碼如下

size_t CalfFluidEngine::PointHashGridSearcher3::GetHashKeyFromBucketIndex(const Vector3<size_t>& bucketIndex) const {
    Vector3<size_t> wrappedIndex = bucketIndex;
    wrappedIndex.x = bucketIndex.x % _resolution.x;
    wrappedIndex.y = bucketIndex.y % _resolution.y;
    wrappedIndex.z = bucketIndex.z % _resolution.z;
    return (wrappedIndex.z 
        * _resolution.y + wrappedIndex.y) 
        * _resolution.x + wrappedIndex.x;
}

Vector3<size_t> CalfFluidEngine::PointHashGridSearcher3::GetBucketIndex(const Vector3D & position) const
{
    Vector3<size_t> bucketIndex;
    bucketIndex.x = static_cast<size_t>(
        std::floor(position.x / _gridSpacing));
    bucketIndex.y = static_cast<size_t>(
        std::floor(position.y / _gridSpacing));
    bucketIndex.z = static_cast<size_t>(
        std::floor(position.z / _gridSpacing));
    return bucketIndex;
}

PointHashGridSearcher3內部有一個二維陣列，或者可以稱它為桶陣列，桶陣列的鍵值就是由3D座標轉換而來的key，代表一個包圍盒區域內的長方形塊。桶陣列儲存所對應長方塊內的所有粒子在粒子陣列中的索引。

當需要查詢鄰域粒子時，就可以通過粒子座標所對應的key值查詢獲取粒子所在長方塊附近其他長方形塊的key，從而獲取粒子所在長方形塊和粒子附近長方形塊內其他粒子。

獲取附近key值程式碼實現如下

void CalfFluidEngine::PointHashGridSearcher3::getNearbyKeys(const Vector3D & position, size_t * nearbyKeys) const
{
    Vector3<size_t> originIndex = GetBucketIndex(position), nearbyBucketIndices[8];

    for (int i = 0; i < 8; i++) {
        nearbyBucketIndices[i] = originIndex;
    }

    if ((originIndex.x + 0.5f) * _gridSpacing <= position.x) {
        nearbyBucketIndices[4].x += 1;
        nearbyBucketIndices[5].x += 1;
        nearbyBucketIndices[6].x += 1;
        nearbyBucketIndices[7].x += 1;
    }
    else {
        nearbyBucketIndices[4].x -= 1;
        nearbyBucketIndices[5].x -= 1;
        nearbyBucketIndices[6].x -= 1;
        nearbyBucketIndices[7].x -= 1;
    }

    if ((originIndex.y + 0.5f) * _gridSpacing <= position.y) {
        nearbyBucketIndices[2].y += 1;
        nearbyBucketIndices[3].y += 1;
        nearbyBucketIndices[6].y += 1;
        nearbyBucketIndices[7].y += 1;
    }
    else {
        nearbyBucketIndices[2].y -= 1;
        nearbyBucketIndices[3].y -= 1;
        nearbyBucketIndices[6].y -= 1;
        nearbyBucketIndices[7].y -= 1;
    }

    if ((originIndex.z + 0.5f) * _gridSpacing <= position.z) {
        nearbyBucketIndices[1].z += 1;
        nearbyBucketIndices[3].z += 1;
        nearbyBucketIndices[5].z += 1;
        nearbyBucketIndices[7].z += 1;
    }
    else {
        nearbyBucketIndices[1].z -= 1;
        nearbyBucketIndices[3].z -= 1;
        nearbyBucketIndices[5].z -= 1;
        nearbyBucketIndices[7].z -= 1;
    }

    for (int i = 0; i < 8; i++) {
        nearbyKeys[i] = GetHashKeyFromBucketIndex(nearbyBucketIndices[i]);
    }

}

獲取粒子附近其他粒子程式碼如下

size_t nearbyKeys[8];
getNearbyKeys(origin, nearbyKeys);

const double queryRadiusSquared = radius * radius;

for (int i = 0; i < 8; i++) {
    const auto& bucket = _buckets[nearbyKeys[i]];
    size_t numberOfPointsInBucket = bucket.size();

    for (size_t j = 0; j < numberOfPointsInBucket; ++j) {
        size_t pointIndex = bucket[j];
        double rSquared = (_points[pointIndex] - origin).SquareMagnitude();
        if (rSquared <= queryRadiusSquared) {
            ......
        }
    }
}

向量運算元

下一節將進入流體力學的部分，在那之前先重溫或學習一下關於矢算場的數學知識，如果不瞭解的話，會對Navier-Stocks方程的一些數學符號感到疑惑。

《Fluid Engine Development》有介紹這部分的知識，更詳細的大家可以閱讀另一位博主推薦的《失算場論札記》,不過這本書我淘寶也買不到，在學校圖書館借了看，後來花9塊錢買了向量分析與場論-工程數學-第四版

偏導數

\(\frac{\partial f}{\Delta } = \frac{f(x + \Delta,y,z) - f(x,y,z)}{ \Delta}\)

梯度

\(\nabla f(x,y,z) = (\frac{\partial f}{\partial x} , \frac{\partial f}{\partial y} , \frac{\partial f}{\partial z} )\)

\(\nabla\) 是梯度運算元，用以測量標量的變化率和變換方向，

梯度是標量場增長最快的方向，梯度的長度是這個最大的變化率

梯度本意是一個向量（向量），當某一函式在某點處沿著該方向的方向導數取得該點處的最大值，即函式在該點處沿方向變化最快，變化率最大（為該梯度的模）

對程式開發者而言，它是由一階偏導組成的向量

當我們需要知道函式在特定方向上的變化率，通過將梯度與該方向的向量點乘，我們就可以計算出函式在這一點的方向導數 $\frac{\partial f}{\partial n} =\nabla f \cdot n $

散度

\(\nabla \cdot f(x) = (\frac{\partial }{\partial x} , \frac{\partial }{\partial y} , \frac{\partial }{\partial z} ) \cdot f(x)\)

$\nabla \cdot $ 是散度運算元，散度就是通量密度，通量即通過一個面的某物理量

在流體力學中，速度場的散度是體積膨脹率，散度(divergence)運算元用於描述向量場中在某點的聚集或發散程度

對程式開發者，這就是梯度運算元與向量場的點乘

流體模擬中，流體往往有著難以壓縮的性質，通常認為速度場的散度為零

旋度

\(\nabla\times f(x) = (\frac{\partial }{\partial x} , \frac{\partial }{\partial y} , \frac{\partial }{\partial z} ) \times f(x)\)

\(\nabla\times\) 是旋度運算元，旋度就是環量密度，以水旋渦為例，沿著圓周關於速度的線積分就是環量，環量除以圓面積，然後讓面積無限小，比值就是旋度，旋度用以描述某個點的旋轉程度

對程式開發者，這就是梯度運算元與向量場的叉積

拉普拉斯運算元

$\nabla^{2} f(x) = \nabla \cdot \nabla f(x) =\frac{\partial^{2} f}{\partial x^{2} } +\frac{\partial^{2} f}{\partial y^{2} } +\frac{\partial^{2} f}{\partial z^{2} } $

$\nabla^{2} $是拉普拉斯運算元，又稱拉氏運算元，它的本質是梯度的散度

拉普拉斯運算元通常用於描述一個函式中某點的值與它鄰域的平均值之差，或者說用於描述某點 在它的鄰域中是“多麼與眾不同”

拉普拉斯運算元可以檢測影象的邊緣和峰值，也可以用以模糊向量場的尖銳部分（拉普拉斯運算元 加上原始向量場）

流體運動的基本原理

現在開始進入正題，前面講的都是程式設計相關的細節，現在開始流體力學的部分。

流體運動的本質是源於多種力的組合。

在保證密度不變的情況下，流體受到重力，壓力，粘度力三種力的作用。

重力不需要我太多介紹，主要是壓力和粘度力

壓力和壓力梯度力

風是高壓區吹到低壓區的，同樣的規則適用於流體。

壓力其實是壓力梯度力，當你潛水時，隨著深度的增加，收到的壓力也隨之增加。這種壓力差會沿著深度在與重力相反的方向產生向上的壓力。正因為這種力，水可以保持其狀態而不是收縮

假設一個立方體大小的流體，單位大小是L，密度為 \(\rho\)， 假設壓力沿x軸方向，左右之間的壓力差是\(\Delta p\)，壓力大小就等於\(\Delta p L^{2}\)，所以 $F_{p} = -\Delta p L^{2}= ma_{p} = L^{3}\rho a_{p} $

可推導得 \(a_{p} = \frac{-\Delta p}{\rho L}\)，假定立方體非常小，可得 $a_{p} = -\frac{\partial p}{\partial x} \frac{1}{\rho} $

接著我們將它推導到三維上 ，可得 $a_{p} = -\frac{\nabla p}{\rho} $

所以$ a = g -\frac{\nabla p}{\rho} + \cdots$

粘度

流體模擬中粘度類似於阻尼，它試圖模糊2點間的速度差異，它導致液體有了一定的粘稠感。正好，拉普拉斯運算元能做到模糊向量場。

\(V_{new} = V + \Delta t u \nabla^{2}V\) u 是粘度係數

由上式可以得到粘度影響的加速度$ a_{V} =u \nabla^{2}V$

所以$ a = g -\frac{\nabla p}{\rho} + u \nabla^{2}V$ 這就是 著名的 納維斯托克斯(Navier--Stokes, NS)方程