在一个solver或者util运行时，会自动加载链接器链接到它身上的库，但想在运行时动态地通过字典文件灵活加载动态库该怎么办呢？可以通过这样的操作：

const_cast<Time&>(io.time()).libs().open
(
    dict,
    "dynamicFvMeshLibs",
    IOobjectConstructorTablePtr_
);

操作会加载dict中dynamicFvMeshLibs声明的库，并通过与IOobjectConstructorTablePtr_ 这个列表对比是否为已经加载过的类。跳转到libs()函数返回的dlLibraryTable对象的open函数

template<class TablePtr>
 bool Foam::dlLibraryTable::open
 (
     const dictionary& dict,
     const word& libsEntry,
     const TablePtr& tablePtr
 )
 {
     fileNameList libNames;
     dict.readIfPresent(libsEntry, libNames);
  
     label nOpen = 0;
  
     for (const fileName& libName : libNames)
     {
         const label nEntries = (tablePtr ? tablePtr->size() : 0);
  
         if (dlLibraryTable::open(libName))
         {
             ++nOpen;
  
             if (debug && (!tablePtr || tablePtr->size() <= nEntries))
             {
                 WarningInFunction
                     << "library " << libName
                     << " did not introduce any new entries"
                     << nl << endl;
             }
         }
         ...

注意传入的tablePtr是const& ，意味着在函数体内不能修改。
dict.readIfPresent(libsEntry, libNames);通过libsEntry读取一系列动态库名，并在for循环中遍历，并通过open函数真正加载动态库。需要注意的是nEntries在open函数前等于tablePtr->size()，而在之后又与tablePtr->size() <= nEntries进行对比，只能说明tablePtr在open中被修改了，但tablePtr的是const的呀！很好理解，传入的IOobjectConstructorTablePtr_是全局静态变量，在函数体外修改是没得关系的。具体就是在open动态库的时候addToRunTimeSelectionTable(dynamicFvMesh, dynamicRefineFvMesh, IOobject);这一句被调用，它会把自己的typeName插入到IOobjectConstructorTablePtr_，从而更新此列表。

你可能要自己写个multiRegionReactingFoam这样的求解器来读取多个fvMesh

你要精确预测轨迹建议还是用DEM吧。

你要的MPPIC大概都在这里了 MPPIC

当然是固相应力呀！Harris & Crighton模型

const labelList face = mesh.faces()[index];
forAll(face, pointi)
{
    auto point = mesh.points()[face[pointi]];
}

修正下

value += this->operator[](i).calcCoupled(p, td, dt, mass, Re, muc);

cloud.drag()[this->cell()] += np0*drag;

cloudVolDrag.correctBoundaryConditions();

这个题我会
第一步，在ParticleForce类中添加虚函数

virtual forceSuSp calcDragForce
(
      const typename CloudType::parcelType& p,
      const typename CloudType::parcelType::trackingData& td,
      const scalar dt,
      const scalar mass,
      const scalar Re,
      const scalar muc
 ) const
{
     return forceSuSp(Zero);
}

第二步，在ErgunWenYuDragForce类添加此函数，并返回

...
{
    return this->calcCoupled
   (
       p,
       td,
       dt,
       mass,
       Re,
       muc
   );

第三步，在ParticleForceList类中添加函数

forceSuSp calcDragForce
(
      const typename CloudType::parcelType& p,
      const typename CloudType::parcelType::trackingData& td,
      const scalar dt,
      const scalar mass,
      const scalar Re,
      const scalar muc
 ) const
{
    forceSuSp value(Zero);

    if (calcCoupled_)
    {
        forAll(*this, i)
        {
            value += this->operator[](i).calcDragForce(p, td, dt, mass, Re, muc);
        }
    }

    return value;
}

当然也可以这么写

...
{
    forceSuSp value(Zero);

    if (calcCoupled_)
    {
        forAll(*this, i)
        {
            if(isA<ErgunWenYuDragForce>(this->operator[](i)))
                 value = this->operator[](i).calcDragForce(p, td, dt, mass, Re, muc);
        }
    }

    return value;
}

这么写的话，前面两步就可以不用做了
第四步，在KinematicParcel类的calcVelocity函数中添加

template<class ParcelType>
template<class TrackCloudType>
const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
(
    TrackCloudType& cloud,
    trackingData& td,
    const scalar dt,
    const scalar Re,
    const scalar mu,
    const scalar mass,
    const vector& Su,
    vector& dUTrans,
    vector& drag,    //添加返回参数
    scalar& Spu
) const
{
    ...
    const forceSuSp FDrag = forces.calcDragForce(p, ttd, dt, mass, Re, mu);
    ...
    drag = -(FDrag.Sp()*td.Uc() + FDrag.Su());
    ...
}

在calc函数中添加

template<class ParcelType>
template<class TrackCloudType>
void Foam::KinematicParcel<ParcelType>::calc
(
    TrackCloudType& cloud,
    trackingData& td,
    const scalar dt
)
{
    ...
    vector drag = Zero;

    // Calculate new particle velocity
    this->U_ =
        calcVelocity(cloud, td, dt, Re, td.muc(), mass0, Su, dUTrans, darg, Spu);

    if (cloud.solution().coupled())
    {
        // Update momentum transfer
        cloud.UTrans()[this->cell()] += np0*dUTrans;

        cloud.dragForce()[this->cell()] += np0*drag;

        // Update momentum transfer coefficient
        cloud.UCoeff()[this->cell()] += np0*Spu;
    }

第五步，在KinematicCloud类中添加变量

autoPtr<volVectorField::Internal> drag_;

构造中添加

    drag_
    (
        new volVectorField::Internal
        (
            IOobject
            (
                this->name() + ":drag",
                this->db().time().timeName(),
                this->db(),
                IOobject::READ_IF_PRESENT,
                IOobject::AUTO_WRITE
            ),
            mesh_,
            dimensionedVector(dimForce, Zero)
        )
    ),

以及函数

template<class CloudType>
inline Foam::DimensionedField<Foam::vector, Foam::volMesh>&
Foam::KinematicCloud<CloudType>::drag()
{
    return *drag_;
}


template<class CloudType>
inline const Foam::DimensionedField<Foam::vector, Foam::volMesh>&
Foam::KinematicCloud<CloudType>::drag() const
{
    return *drag_;
}

第六步，主程序添加

        auto cloudDrag(kinematicCloud.drag());
        volVectorField cloudVolDrag
        (
            IOobject
            (
                "cloudVolDrag",
                runTime.timeName(),
                mesh
            ),
            mesh,
            dimensionedVector(cloudDrag.dimensions()/dimVolume, Zero),
            zeroGradientFvPatchVectorField::typeName
        );

        cloudVolDrag.primitiveFieldRef() = -cloudDrag/mesh.V();
        cloudVolSUSu.correctBoundaryConditions();

大概是这样，没编译过，新鲜的。里面的量纲和曳力应该还需要调整。

...

其实大部分情况下，一般都只会用到曳力，所以可以直接用kinematicCloud.SU(Uc)就行了，那这样的话，前面所有的都不用做了

class base
{
private:
    bool ownedByRegistry_;
public:
    template<class Type>
    static void store(Type* p)
    {
	//p->base::ownedByRegistry_ =true;
	p->ownedByRegistry_ = true;
	p->printInfo();
    }
};

class derived
:public base
{
public:
    derived() {}
    void printInfo()
    {
	Info << "derived class" << endl;
    }
};

int main(int argc, char *argv[])
{
...
    derived d;
    base::store(&d);
...
}

wordHashSet unMatchedkeys(dict.toc());
forAll(names, nameI)
{
    const word& name = names[nameI];
    const entry* ePtr = dict.findEntry(name, keyType::REGEX);
    if(ePtr)
       unMatchedKeys.erase(ePtr->keyWord());
}

valueInternalCoeffs, valueBoundaryCoeffs, gradientInternalCoeffs, gradientBoundaryCoeffs are all key coeffs for boundary correct.

In calculatedFvPatchField, those are all:

tmp<Field<Type> > calculatedFvPatchField<Type>::valueInternalCoeffs
(
    const tmp<scalarField>&
) const
{
    FatalErrorIn
    (
        "calculatedFvPatchField<Type>::"
        "valueInternalCoeffs(const tmp<scalarField>&) const"
    )   << "\n    "
...

when construct convection tern

fvm.internalCoeffs()[patchI] = patchFlux*psf.valueInternalCoeffs(pw);
fvm.boundaryCoeffs()[patchI] = -patchFlux*psf.valueBoundaryCoeffs(pw);

and laplacian term

fvm.internalCoeffs()[patchI] = patchGamma*psf.gradientInternalCoeffs();
fvm.boundaryCoeffs()[patchI] = -patchGamma*psf.gradientBoundaryCoeffs();

so...