Palabos源码：src/particles

h和hh文件
#include "particles/particle3D.hh"
#include "particles/particleIdentifiers3D.hh"
#include "particles/particleField3D.hh"
#include "particles/multiParticleField3D.hh"
#include "particles/particleProcessingFunctional3D.hh"
#include "particles/particleNonLocalTransfer3D.hh"
#include "particles/visualParticle3D.hh"
#include "particles/visualParticleFunctional3D.hh"
#include "particles/visualParticleWrapper3D.hh"
#include "particles/particleVtk3D.hh"

particle3D.h/hh

class Particle3d

Particle3D(plint tag_, Array<T,3> const& position_);

tag并不需要独一无二的存在，仅仅是为了方便使用，可见这里的Particles3D只包含一个tag和一个position。
其定义

template<typename T, template<typename U> class Descriptor>
Particle3D<T,Descriptor>::Particle3D()
    : tag(0),
      position(T(),T(),T())
{ }

template<typename T, template<typename U> class Descriptor>
Particle3D<T,Descriptor>::Particle3D(plint tag_, Array<T,3> const& position_)
    : tag(tag_),
      position(position_)
{ }

virtual void velocityToParticle(TensorField3D<T,3>& velocityField, T scaling=1.) =0;
virtual void velocityToParticle(NTensorField3D<T>& velocityField, T scaling=1.) =0;
virtual void rhoBarJtoParticle(NTensorField3D<T>& rhoBarJfield, bool velIsJ, T scaling=1.) =0;
virtual void fluidToParticle(BlockLattice3D<T,Descriptor>& fluid, T scaling=1.) =0;

耦合particle与fluid，注释里提到流体格子速度需要小于0.25，如果这个条件被违背，那么速度会被削减到0.25 速度的模。

    virtual void advance() =0;
    Array<T,3> const& getPosition() const { return position; }
    Array<T,3>& getPosition() { return position; }
    virtual void reset(Array<T,3> const& position);
    virtual void serialize(HierarchicSerializer& serializer) const;
    virtual void unserialize(HierarchicUnserializer& unserializer);
    virtual int getId() const =0;
    plint getTag() const;
    void setTag(plint tag_);
    virtual Particle3D<T,Descriptor>* clone() const =0;
    virtual bool getScalar(plint whichScalar, T& scalar) const;
    virtual bool setScalar(plint whichScalar, T scalar);
    virtual bool getVector(plint whichVector, Array<T,3>& vector) const;
    virtual bool setVector(plint whichVector, Array<T,3> const& vector);
    virtual bool getTensor(plint whichVector, Array<T,SymmetricTensorImpl<T,3>::n>& tensor) const;
    virtual bool setTensor(plint whichVector, Array<T,SymmetricTensorImpl<T,3>::n> const& tensor);
    virtual bool setScalars(std::vector<T> const& scalars);
    virtual bool setVectors(std::vector<Array<T,3> > const& vectors);
    virtual bool setTensors(std::vector<Array<T,SymmetricTensorImpl<T,3>::n> > const& tensors);
    virtual void rescale(int dxScale, int dtScale);

存在的Particle3d操作，比如getPosition会返回坐标。

PointParticle3D

template<typename T, template<typename U> class Descriptor>
void PointParticle3D<T,Descriptor>::reset(Array<T,3> const& position_)
{
    Particle3D<T,Descriptor>::reset(position_);
    velocity.resetToZero();
}

reset命令重置particle位置，速度清零。

template<typename T, template<typename U> class Descriptor>
void PointParticle3D<T,Descriptor>::velocityToParticle(TensorField3D<T,3>& velocityField, T scaling)
{
    velocity = predictorCorrectorTensorField<T,3>(velocityField, this->getPosition(), scaling);
}

template<typename T, template<typename U> class Descriptor>
void PointParticle3D<T,Descriptor>::velocityToParticle(NTensorField3D<T>& velocityField, T scaling)
{
    velocity = predictorCorrectorNTensorField<T>(velocityField, this->getPosition(), scaling);
}

利用了一个predictor函数计算速度。

template<typename T, template<typename U> class Descriptor>
void PointParticle3D<T,Descriptor>::rhoBarJtoParticle (
        NTensorField3D<T>& rhoBarJfield, bool velIsJ, T scaling )
{
    T rhoBar;
    Array<T,3> j;
    predictorCorrectorRhoBarJ(rhoBarJfield, this->getPosition(), velIsJ, j, rhoBar);
    if (velIsJ) {
        velocity = j*scaling;
    }
    else {
        velocity = j*scaling*Descriptor<T>::invRho(rhoBar);
    }
}

在某些算例里，流场信息是以rhoBarJary形式储存的，则以这种方式计算particle速度。

template<typename T, template<typename U> class Descriptor>
void PointParticle3D<T,Descriptor>::fluidToParticle(BlockLattice3D<T,Descriptor>& fluid, T scaling)
{
    static const T maxVel = 0.25-1.e-6;
    static const T maxVelSqr = maxVel*maxVel;
#ifdef PLB_DEBUG
    Box3D bbox(fluid.getBoundingBox());
#endif
    Dot3D loc(fluid.getLocation());
    Array<T,3> position1(this->getPosition()-Array<T,3>(loc.x,loc.y,loc.z));
    PLB_ASSERT( position1[0] >= bbox.x0+0.5 );
    PLB_ASSERT( position1[0] <= bbox.x1-0.5 );
    PLB_ASSERT( position1[1] >= bbox.y0+0.5 );
    PLB_ASSERT( position1[1] <= bbox.y1-0.5 );
    PLB_ASSERT( position1[2] >= bbox.z0+0.5 );
    PLB_ASSERT( position1[2] <= bbox.z1-0.5 );

    Dot3D intPos( (plint)position1[0], (plint)position1[1], (plint)position1[2] );
    T u = position1[0] - (T)intPos.x;
    T v = position1[1] - (T)intPos.y;
    T w = position1[2] - (T)intPos.z;

    Array<T,3> tmpVel;
    Array<T,3> velocity1; velocity1.resetToZero();
    fluid.get(intPos.x  ,intPos.y  ,intPos.z  ).computeVelocity(tmpVel);
    velocity1 += ((T)1.-u) * ((T)1.-v) * ((T)1.-w) *tmpVel;
    fluid.get(intPos.x  ,intPos.y  ,intPos.z+1).computeVelocity(tmpVel);
    velocity1 += ((T)1.-u) * ((T)1.-v) * (      w) *tmpVel;
    fluid.get(intPos.x  ,intPos.y+1,intPos.z  ).computeVelocity(tmpVel);
    velocity1 += ((T)1.-u) * (      v) * ((T)1.-w) *tmpVel;
    fluid.get(intPos.x  ,intPos.y+1,intPos.z+1).computeVelocity(tmpVel);
    velocity1 += ((T)1.-u) * (      v) * (      w) *tmpVel;
    fluid.get(intPos.x+1,intPos.y  ,intPos.z  ).computeVelocity(tmpVel);
    velocity1 += (      u) * ((T)1.-v) * ((T)1.-w) *tmpVel;
    fluid.get(intPos.x+1,intPos.y  ,intPos.z+1).computeVelocity(tmpVel);
    velocity1 += (      u) * ((T)1.-v) * (      w) *tmpVel;
    fluid.get(intPos.x+1,intPos.y+1,intPos.z  ).computeVelocity(tmpVel);
    velocity1 += (      u) * (      v) * ((T)1.-w) *tmpVel;
    fluid.get(intPos.x+1,intPos.y+1,intPos.z+1).computeVelocity(tmpVel);
    velocity1 += (      u) * (      v) * (      w) *tmpVel;

    velocity1 *= scaling;

    if (normSqr(velocity1)>maxVelSqr) {
        velocity1 /= norm(velocity1);
        velocity1 *= maxVel;
    }

    Array<T,3> position2(position1+velocity1);

    intPos = Dot3D( (plint)position2[0], (plint)position2[1], (plint)position2[2] );
    u = position2[0] - (T)intPos.x;
    v = position2[1] - (T)intPos.y;
    w = position2[2] - (T)intPos.z;

    Array<T,3> velocity2; velocity2.resetToZero();
    fluid.get(intPos.x  ,intPos.y  ,intPos.z  ).computeVelocity(tmpVel);
    velocity2 += ((T)1.-u) * ((T)1.-v) * ((T)1.-w) *tmpVel;
    fluid.get(intPos.x  ,intPos.y  ,intPos.z+1).computeVelocity(tmpVel);
    velocity2 += ((T)1.-u) * ((T)1.-v) * (      w) *tmpVel;
    fluid.get(intPos.x  ,intPos.y+1,intPos.z  ).computeVelocity(tmpVel);
    velocity2 += ((T)1.-u) * (      v) * ((T)1.-w) *tmpVel;
    fluid.get(intPos.x  ,intPos.y+1,intPos.z+1).computeVelocity(tmpVel);
    velocity2 += ((T)1.-u) * (      v) * (      w) *tmpVel;
    fluid.get(intPos.x+1,intPos.y  ,intPos.z  ).computeVelocity(tmpVel);
    velocity2 += (      u) * ((T)1.-v) * ((T)1.-w) *tmpVel;
    fluid.get(intPos.x+1,intPos.y  ,intPos.z+1).computeVelocity(tmpVel);
    velocity2 += (      u) * ((T)1.-v) * (      w) *tmpVel;
    fluid.get(intPos.x+1,intPos.y+1,intPos.z  ).computeVelocity(tmpVel);
    velocity2 += (      u) * (      v) * ((T)1.-w) *tmpVel;
    fluid.get(intPos.x+1,intPos.y+1,intPos.z+1).computeVelocity(tmpVel);
    velocity2 += (      u) * (      v) * (      w) *tmpVel;
    
    velocity2 *= scaling;

    if (normSqr(velocity2)>maxVelSqr) {
        velocity2 /= norm(velocity2);
        velocity2 *= maxVel;
    }

    velocity = (velocity1+velocity2)/(T)2;
}

从最上面开始，fluid.getBoundingBox()得到了整个流体的方形体的坐标，然后后续6行判断大小是拿particle的坐标与方形体对角线两个点的坐标对比，是为了确定particle在流场内部。
其uvw获取也是挺有意思，整型一下坐标，然后再减去它就得到了数值。
接着定义velocity1之后，很快又将其设置为0，我想这一步不可避免，在我之前开发算法的时候，感觉到计算数值不对，似乎就是在声明定义之后没有设置为0，具体原因不得而知，毕竟我没学过c++，也不想学c++。
tmpVel会被赋予整型后坐标那个点的流体速度，然后是8个方向根据uvw，有点求均值的味道，得到了velocity1，scaling应该是与单位转换有关，后续判断速度的模是否大于速度为0.25的模，如果是就强行定义为速度为0.25的模。

后面position2实际上就是预测了以velocity1速度流动一个时间步长的位置，在这里位置也求一个均值速度，再与velocity1平均一下，就得到了particle的速度。

template<typename T, template<typename U> class Descriptor>
void PointParticle3D<T,Descriptor>::advance() {
    PLB_ASSERT( norm(velocity)<1. );
    this->getPosition() += velocity;
}

更新particle的位置。

NormedVelocityParticle3D, RestParticle3D, VerletParticle3D

略

particleIdentifiers3D.h/hh

略

particleField3D.h/hh

class ParticleField3D

ParticleField3D(plint nx, plint ny, plint nz, BlockDataTransfer3D* dataTransfer);

可见颗粒区域由一个3d区域构成，

 virtual void addParticle(Box3D domain, Particle3D<T,Descriptor>* particle) =0;
    /// Remove all particles found in the indicated domain.
    virtual void removeParticles(Box3D domain) =0;
    /// Remove all particles of a certain tag found in the indicated domain.
    virtual void removeParticles(Box3D domain, plint tag) =0;
    /// Return all particles found in the indicated domain.
    virtual void findParticles(Box3D domain,
                               std::vector<Particle3D<T,Descriptor>*>& found) =0;
    /// Return all particles found in the indicated domain; const version
    virtual void findParticles(Box3D domain,
                               std::vector<Particle3D<T,Descriptor> const*>& found) const =0;
    /// Execute velocity-particle interaction for all particles contained in the domain.
    virtual void velocityToParticleCoupling(Box3D domain, TensorField3D<T,3>& velocity, T scaling=0.) =0;
    /// Execute velocity-particle interaction for all particles contained in the domain.
    virtual void velocityToParticleCoupling(Box3D domain, NTensorField3D<T>& velocity, T scaling=0.) =0;
    /// Execute velocity-particle interaction for all particles contained in the domain.
    virtual void rhoBarJtoParticleCoupling(Box3D domain, NTensorField3D<T>& rhoBarJ, bool velIsJ, T scaling=0.) =0;
    /// Execute fluid-particle interaction for all particles contained in the domain.
    virtual void fluidToParticleCoupling(Box3D domain, BlockLattice3D<T,Descriptor>& lattice, T scaling=0.) =0;
    /// Advance all particles contained in the domain. When the speed of a particle drops
    ///   below sqrt(cutOffSpeedSqr), the particle is eliminated. Negative cutOffSpeedSqr means
    ///   no cutoff.
    virtual void advanceParticles(Box3D domain, T cutOffSpeedSqr=-1.) =0;
    virtual identifiers::BlockId getBlockId() const { return identifiers::ParticleId; }

与上面的只针对particle不同，这里是针对一个particle field。

template<typename T, template<typename U> class Descriptor>
void DenseParticleField3D<T,Descriptor>::velocityToParticleCoupling (
        Box3D domain, TensorField3D<T,3>& velocityField, T scaling )
{
    Box3D finalDomain;
    if( intersect(domain, particleGrid.getBoundingBox(), finalDomain) )
    {
        for (plint iX=finalDomain.x0; iX<=finalDomain.x1; ++iX) {
            for (plint iY=finalDomain.y0; iY<=finalDomain.y1; ++iY) {
                for (plint iZ=finalDomain.z0; iZ<=finalDomain.z1; ++iZ) {
                    for (pluint iParticle=0; iParticle<particleGrid.get(iX,iY,iZ).size(); ++iParticle) {
                        particleGrid.get(iX,iY,iZ)[iParticle]->velocityToParticle(velocityField, scaling);
                    }
                }
            }
        }
    }
}

template<typename T, template<typename U> class Descriptor>
void DenseParticleField3D<T,Descriptor>::fluidToParticleCoupling (
        Box3D domain, BlockLattice3D<T,Descriptor>& lattice, T scaling )
{
    Box3D finalDomain;
    if( intersect(domain, particleGrid.getBoundingBox(), finalDomain) )
    {
        for (plint iX=finalDomain.x0; iX<=finalDomain.x1; ++iX) {
            for (plint iY=finalDomain.y0; iY<=finalDomain.y1; ++iY) {
                for (plint iZ=finalDomain.z0; iZ<=finalDomain.z1; ++iZ) {
                    for (pluint iParticle=0; iParticle<particleGrid.get(iX,iY,iZ).size(); ++iParticle) {
                        particleGrid.get(iX,iY,iZ)[iParticle]->fluidToParticle(lattice, scaling);
                    }
                }
            }
        }
    }
}

而其hh文件中的代码表明是对particlefield上的每一个particle进行fluidToParticle，velocityToParticle等操作。

DenseParticleField3D，LightParticleDataTransfer3D，等

略

multiParticleField3D

节选如下：
MultiParticleField3D (
            MultiBlockManagement3D const& multiBlockManagement_,
            CombinedStatistics* combinedStatistics_ );
    MultiParticleField3D(plint nx_, plint ny_, plint nz_);

在算例里，实际应用是这样的：

// Definition of a particle field.
    MultiParticleField3D<DenseParticleField3D<T,DESCRIPTOR> >* particles=0;
    particles = new MultiParticleField3D<DenseParticleField3D<T,DESCRIPTOR> > (
        lattice->getMultiBlockManagement(),
        defaultMultiBlockPolicy3D().getCombinedStatistics() );

particleProcessingFunctional3D.h/hh

class CountParticlesFunctional3D

计算所有种类particle数量。
这里用到了一个数据处理器，我们来看看数据处理器是怎么样的。

    virtual void processGenericBlocks(Box3D domain, std::vector<AtomicBlock3D*> fields);
    virtual CountParticlesFunctional3D<T,Descriptor>* clone() const;
    plint getNumParticles() const;

void CountParticlesFunctional3D<T,Descriptor>::processGenericBlocks (
        Box3D domain, std::vector<AtomicBlock3D*> blocks )
{
    PLB_PRECONDITION( blocks.size()==1 );
    ParticleField3D<T,Descriptor>& particleField
        = *dynamic_cast<ParticleField3D<T,Descriptor>*>(blocks[0]);
    std::vector<Particle3D<T,Descriptor>*> particles;
    particleField.findParticles(domain, particles);
    this->getStatistics().gatherIntSum(numParticlesId, (plint)particles.size());
}

可见，主要的代码是particleField.findParticles(domain, particles);。

在算例中，主要的应用则是依靠如下代码使用：

 countParticles(*particles, particles->getBoundingBox())

class AverageParticleVelocityFunctional3D

返回particles平均速度。

CopySelectParticles3D

略

InjectParticlesFunctional3D

略

InjectRandomParticlesFunctional3D

基本同上，但是以概率形式生成particles。

MaskedInjectRandomParticlesFunctional3D

基本同上，但生成particles的坐标是以布尔模板定义的。

N_MaskedInjectRandomParticlesFunctional3D， AnalyticalInjectRandomParticlesFunctional3D， MaskedAnalyticalInjectRandomParticlesFunctional3D

略

InjectEquallySpacedParticlesFunctional3D

每个cell生成x方向生成nx个颗粒，y方向生成ny个颗粒，z方向生成nz的课题，即每个cell生成nxxyxz个颗粒。

MaskedInjectEquallySpacedParticlesFunctional3D

略

略

AbsorbParticlesFunctional3D

在这个区域，particles会被移除。

略

RemoveParticlesFromWall3D

将避面内部的particles移除

PushParticlesAwayFromWall3D

将壁面内部的particles推回流场里面。

FluidToParticleCoupling3D

执行流体颗粒耦合，这个过程颗粒不会移动，流场不变。

    FluidToParticleCoupling3D(T scaling_);

particleField.fluidToParticleCoupling(domain, fluid, scaling);

即执行particleField的流体和颗粒耦合，在之前的particleField里，实际上最后是对所有的particle逐个进行计算的。

VelocityToParticleCoupling3D

    VelocityToParticleCoupling3D(T scaling_);

particleField.velocityToParticleCoupling(domain, velocity, scaling);

即执行particleField的速度和颗粒耦合，在之前的particleField里，实际上最后是对所有的particle逐个进行计算的。

N_VelocityToParticleCoupling3D， RhoBarJtoParticleCoupling3D

略

AdvanceParticlesFunctional3D

AdvanceParticlesFunctional3D(T cutOffValue_ = -1.);

主要代码如上，如果particle速度低于cutOffValue_的平方根就会被移除，但如果是负数，则意味着不会移除。

AdvanceParticlesEveryWhereFunctional3D， VerletUpdateVelocity3D， VerletUpdateVelocitySelective3D， VerletParticleInteractionForce3D, CountAndAccumulateTaggedParticlesRefined3D, CountTaggedParticles3D, ComputeConcentrationOfTaggedParticles3D, MaskedComputeConcentrationOfTaggedParticles3D, N_MaskedComputeConcentrationOfTaggedParticles3D

略

particleNonLocalTransfer3D.hh, particles/visualParticle3D.hh, particles/visualParticleFunctional3D.hh, particles/visualParticleWrapper3D.hh, particles/particleVtk3D.hh

略