颗粒间相互作用设计¶

颗粒之间相互作用搜索¶

颗粒之间相互作用搜索主要由KDTreeFlann类实现，该类调用nanoflann开源代码实现颗粒间相互作用的搜索，关于本开源代码的使用请参考以下链接nanoflann，下面进行介绍：

PointCloud类本类定义在文件PhyStar/Algo/KNN/utils.h中，主要用来构建颗粒云：

数据

struct Point
{
    T x, y, z;
};

std::vector<Point> pts;

方法

返回颗粒云中颗粒的数量：

inline size_t kdtree_get_point_count() const 
{ 
    return pts.size(); 
}

返回颗粒云中第idx个颗粒的第dim个位置分量：

inline T kdtree_get_pt(const size_t idx, const size_t dim) const
    {
        if (dim == 0)
            return pts[idx].x;
        else if (dim == 1)
            return pts[idx].y;
        else
            return pts[idx].z;
    }

暂时理解为占位函数吧：

    template <class BBOX>
    bool kdtree_get_bbox(BBOX & /* bb */) const
    {
        return false;
    }

KDTreeFlann类

数据

PointCloud<Real> cloud;

// construct a kd-tree index:
using my_kd_tree_t = nanoflann::KDTreeSingleIndexAdaptor<
            nanoflann::L2_Simple_Adaptor<Real, PointCloud<Real>>,
            PointCloud<Real>, N_dim>;

my_kd_tree_t *ptr_index = nullptr;
std::vector<nanoflann::ResultItem<uint32_t, Real>> ret_matches;
std::vector<uint32_t> indices;
std::vector<Real> dist;

方法

建立颗粒云：

auto &BuildFromIterator(auto itbegion, auto itend, const size_t maxleaf = 10)
{
    cloud.pts.clear();
    int ip = 0;
    PointCloud<Real>::Point p = {0, 0, 0};
    for (auto it = itbegion; it != itend; it++)
    {
        p.x = (*it).pdata.pos(0);
        if constexpr (N_dim > 1)
        {
            p.y = (*it).pdata.pos(1);
        }
        if constexpr (N_dim > 2)
        {
            p.z = (*it).pdata.pos(2);
        }
        cloud.pts.push_back(p);      // add a particle to the cloud
    }
    if (ptr_index != nullptr)
    {
        delete ptr_index;
        ptr_index = nullptr;
    }
    ptr_index = new my_kd_tree_t(N_dim, cloud, {maxleaf});
    return *this;
};

搜索颗粒附近给定距离的颗粒，返回邻近颗粒的index和对应的距离distance：

std::pair<std::vector<uint32_t>, std::vector<Real>> Search(const Real pos[], const Real radius)
{
    ret_matches.clear();

    auto nMatches = ptr_index->radiusSearch(pos, radius * radius, ret_matches);

    indices.resize(nMatches);
    dist.resize(nMatches);
    for (size_t i = 0; i < nMatches; i++)
    {
        indices[i] = ret_matches[i].first;
        dist[i] = std::sqrt(ret_matches[i].second);
    }
    return {indices, dist};
}

颗粒间的接触设计¶

颗粒间的接触主要由PPContact类进行控制，其记录可能发生碰撞的颗粒和不可能发生碰撞的颗粒，记录接触模型等等，下面进行介绍：

数据

ptr_particle_base存储颗粒的模型列表，在函数Init()进行设置；
contact_comp存储发生接触的颗粒之间的类型，控制碰撞是否只发生在同一类型的颗粒(Single)或不同类型的颗粒(Multi)之间，在函数Init()进行设置；
contact_models存储颗粒的接触模型列表，在函数ContactModelListGenerate()中进行设置；

// storing index for particles which could be contacted
std::vector<int> possible_contact_index;
// storing index for particles which cannot be contacted 
std::vector<int> no_possible_contact_index;

algo::KDTreeFlann<N_dim> tree; // estabilsh  Flann KDtree 

using CM_ptr = std::shared_ptr<ContactModel>;
std::vector<CM_ptr> contact_models;

/**
    * @brief hash of std::pair to establish a hash map
    *
    */
struct Pair_Hash
{
    template <typename T1, typename T2>
    std::size_t operator()(const std::pair<T1, T2>& p) const
    {
        auto h1 = std::hash<T1>{}(p.first);
        auto h2 = std::hash<T2>{}(p.second);
        return h1 ^ h2;
    }
};

std::unordered_map<std::pair<uIntl, uIntl>, ContactStatus, Pair_Hash> contact_pairs;
// the base of pointer list of pointer of PointParticleModel
PointParticleModel** ptr_particle_base = nullptr; 
ContactComponent contact_comp = ContactComponent::Single;
uIntl coll_count = 0;

方法

为每个颗粒搜索邻近的颗粒(不包括虚拟颗粒)，将有可能发生接触作用的颗粒记录在PPContact类的possible_contact_index变量中，将不可能发生接触作用的颗粒记录在no_possible_contact_index变量中，并将颗粒的邻居颗粒记录在ParticleField类中的nbrs_pp_index变量中。

void SearchingParticlePairs(auto& pflist, const Real search_radius)
{
    possible_contact_index.clear();
    no_possible_contact_index.clear();
    tree.BuildFromIterator(pflist.begin(), pflist.vend(), 10);
    for (auto i = 0; i < pflist.MaxLength(); ++i)
    {
        auto [ind, dist] = tree.Search(pflist(i).pdata.pos.data(), search_radius);
        pflist(i).nbrs_pp_index.clear();
        if (ind.size() <= 1)
        {
            no_possible_contact_index.push_back(ind[0]); // the self particle
        }
        else
        {
            possible_contact_index.push_back(ind[0]);
            pflist(i).nbrs_pp_index.resize(ind.size() - 1);
            for (auto ip = 0; ip < ind.size() - 1; ++ip)
            {
                pflist(i).nbrs_pp_index[ip] = ind[ip + 1];
            }
        }
    }
}

清除已经记录的接触的颗粒对：

inline void ClearContactPairs()
{
    contact_pairs.clear();
}

添加已经接触的颗粒对，将其存储在变量contact_pairs中：

inline void AddContactPair(uIntl i, uIntl j)
{
    if (i > j)
    {
        std::swap(i, j);
    }
    std::pair<uIntl, uIntl> ij = { i, j };
    if (contact_pairs.count({ i, j }) < 1)
    {
        contact_pairs.emplace(ij, ContactStatus::Contact);
        // std::cout<<i<<"\t"<<j<<std::endl;
    }
}

记录颗粒碰撞次数：

inline void CollisionCount(uIntl i, uIntl j, const ContactStatus cont_status, const uIntl maxlength)
{
    if (cont_status == ContactStatus::Contact)
    {
        return;
    }
    else
    {
        if (i > j)
        {
            std::swap(i, j);
        }
        uIntl increment = 2;
        if (contact_pairs.count({ i, j }) == 1)
        {
            if (contact_pairs[{i, j}] == ContactStatus::Contact)
            {
                if (j >= maxlength)
                {
                    increment = 1;         // collide with the virtual particle only counts once
                }
                coll_count += increment;   // collide with the particle counts twice
                contact_pairs[{i, j}] = ContactStatus::Do_Not_Contact;
                // std::cout<<i<<"\t"<<j<<std::endl;
            }
        }
    }
}

初始化颗粒碰撞次数，将其置为零：

inline void InitCollisionCount()
{
    coll_count = 0;
}

生成接触模型列表，并存储在contact_models变量中，在颗粒求解器中调用函数CreateParticle()进行设置：

void ContactModelListGenerate(const std::vector<std::string>& strlist) {
    contact_models.clear();
    for (auto& str : strlist)
    {
        if (str == "GhostCollision")
        {
            contact_models.push_back(CM_ptr(new typename PPContactModel<N_dim>::GhostCollision(ptr_particle_base)));
        }
        else if (str == "GlowinskiCollision")
        {
            contact_models.push_back(CM_ptr(new typename PPContactModel<N_dim>::GlowinskiCollision(ptr_particle_base)));
        }
        else
        {
            std::string serr = "[Error At PPContact::ContactModelListGenerate]:" + str + " was not defined\n";
            throw serr;
        }
    }
}

初始化，设置颗粒间接触的类型，是否只限定为同种颗粒间的接触还是多种颗粒间的接触，并且传递颗粒模型列表，这个函数在颗粒求解器里面调用函数Prepare()进行设置：

inline void Init(PointParticleModel** ptr_p, const ContactComponent contactcomp = ContactComponent::Single)
{
    ptr_particle_base = ptr_p;
    // cont_model = collmodel;
    contact_comp = contactcomp;
}

调用对应的接触模型进行接触力的计算，并将接触力存储在对应的颗粒的ParticleField类的contact_force变量中，同时计算颗粒间的接触力矩，并将其存储在存储在对应的颗粒的ParticleField类的contact_torque变量中，并返回颗粒之间的接触状态，这个函数在颗粒求解器里面调用函数vSolve()进行设置：

ContactStatus Force(auto& pf_i, auto& pf_j, const Real alpha = 1., const Real beta = 0., const bool whether_only_detect_contact = false)
{
    if ((contact_comp == ContactComponent::Single) & (pf_i.particle_typeid != pf_j.particle_typeid))
    {
        return ContactStatus::Do_Not_Contact;
    }
    if (whether_only_detect_contact) {
        return WhetherContacted(pf_i, pf_j);
    }
    else {
        for (auto& cm : contact_models) {
            cm->Force(pf_i, pf_j, alpha, beta);
        }
        return WhetherContacted(pf_i, pf_j);
    }
}

返回颗粒间的接触状态，是接触还是没有接触：

ContactStatus WhetherContacted(auto& pf_i, auto& pf_j) {
    if (WhichContactShapeType(pf_i, pf_j) == ContactShapeType::SphereToSphere) {
        auto dist = Distance(pf_i, pf_j);
        auto p_part0 = *(ptr_particle_base + pf_i.particle_typeid);
        auto p_part1 = *(ptr_particle_base + pf_j.particle_typeid);
        auto r0 = p_part0->Radius();
        auto r1 = p_part1->Radius();
        auto delta = r0 + r1 - dist;
        if (delta > 0)
        {
            return ContactStatus::Contact;
        }else{
            return ContactStatus::Do_Not_Contact;
        }

    }
    else if (WhichContactShapeType(pf_i, pf_j) == ContactShapeType::EllipsoidToEllipsoid) {
        return ContactStatus::Do_Not_Contact;

    }
    else {
        return ContactStatus::Do_Not_Contact;
    }
}

返回颗粒之间的碰撞类型：

ContactShapeType WhichContactShapeType(auto& pf_i, auto& pf_j)
{
    auto p_part0 = *(ptr_particle_base + pf_i.particle_typeid);
    auto p_part1 = *(ptr_particle_base + pf_j.particle_typeid);
    auto shapeid0 = p_part0->ShapeID();
    auto shapeid1 = p_part1->ShapeID();
    if (shapeid0 == ShapeType::Sphere & shapeid1 == ShapeType::Sphere)
    {
        return ContactShapeType::SphereToSphere;
    }
    else if (shapeid0 == ShapeType::Spheroid | shapeid0 == ShapeType::Ellipsoid |
        shapeid1 == ShapeType::Spheroid | shapeid1 == ShapeType::Ellipsoid)
    {
        return ContactShapeType::EllipsoidToEllipsoid;
    }
    else
    {
        return ContactShapeType::EllipsoidToEllipsoid;
    }
}

颗粒接触模型基类¶

颗粒接触模型¶

颗粒相互作用的一些枚举类型¶

控制颗粒之间的碰撞类型：

enum class ContactShapeType
{
    None = 0,
    SphereToSphere,
    EllipsoidToEllipsoid
};

控制颗粒之间的接触状态：

enum class ContactStatus
{
    Do_Not_Contact = 0, //contact model does not work yet
    Contact             //contact model works
};

控制碰撞是否只发生在同一类型的颗粒(Single)或不同类型的颗粒(Multi)之间
```
enum class ContactComponent
{
    Single = 0,
    Multi,
};
```