main.cpp

#include "Vector/vector_dist.hpp"
#include "data_type/aggregate.hpp"
#include "Plot/GoogleChart.hpp"
#include "Plot/util.hpp"
#include "timer.hpp"
 
int main(int argc, char* argv[])
{
 
    openfpm_init(&argc,&argv);
 
    constexpr int gid = 0;
    constexpr int nn_norm = 1;
    constexpr int nn_sym = 2;
 
    // used to define the domain
    float L = 1000.0;
 
    // Domain
    Box<3,float> box({-L,-L,-L},{L,L,L});
 
    // Boundary conditions
    size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
 
    // cut-off radius
    float r_cut = 100.0;
 
    // ghost
    Ghost<3,float> ghost(r_cut);
 
    // Point and global id
    struct point_and_gid
    {
        // global id of the particle
        size_t id;
 
        // Position of the neighborhood particle
        Point<3,float> xq;
 
        // Used to reorder the neighborhood particles by id
        bool operator<(const struct point_and_gid & pag) const
        {
            return (id < pag.id);
        }
    };
 
 
 
    // Particle properties list
    typedef  aggregate<size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;
 
 
 
    // Distributed vector
    vector_dist<3,float, part_prop > vd(4096,box,bc,ghost);
 
    // used to calculate a global index for each particle
    size_t start = vd.accum();
 
    auto it = vd.getDomainIterator();
 
    while (it.isNext())
    {
        auto key = it.get();
 
        vd.getPos(key)[0] = 2.0*L*((float)rand()/RAND_MAX) - L;
        vd.getPos(key)[1] = 2.0*L*((float)rand()/RAND_MAX) - L;
        vd.getPos(key)[2] = 2.0*L*((float)rand()/RAND_MAX) - L;
 
        vd.getProp<gid>(key) = key.getKey() + start;
 
        ++it;
    }
 
 
 
    vd.map();
 
    // sync the ghost
    vd.ghost_get<gid>();
 
 
 
    auto NN = vd.getCellList(r_cut);
    auto p_it = vd.getDomainIterator();
 
    while (p_it.isNext())
    {
        auto p = p_it.get();
 
        Point<3,float> xp = vd.getPos(p);
 
        auto Np = NN.getNNIterator(NN.getCell(vd.getPos(p)));
 
        while (Np.isNext())
        {
            auto q = Np.get();
 
            // skip self interaction
            if (p.getKey() == q)
            {
                ++Np;
                continue;
            }
 
            Point<3,float> xq = vd.getPos(q);
            Point<3,float> f = (xp - xq);
 
            float distance = f.norm();
 
            // if the distance smalle than the cut-off radius add it to the neighborhood list
            if (distance < r_cut )
            {
                vd.getProp<nn_norm>(p).add();
                vd.getProp<nn_norm>(p).last().xq = xq;
                vd.getProp<nn_norm>(p).last().id = vd.getProp<0>(q);
            }
 
            // Next neighborhood
            ++Np;
        }
 
        // Next particle
        ++p_it;
    }
 
 
 
    auto NN2 = vd.getCellListSym(r_cut);
 
    auto p_it2 = vd.getDomainIterator();
 
    while (p_it2.isNext())
    {
        auto p = p_it2.get();
 
        Point<3,float> xp = vd.getPos(p);
 
        auto Np = NN2.template getNNIteratorSym<NO_CHECK>(NN2.getCell(vd.getPos(p)),p.getKey(),vd.getPosVector());
 
        while (Np.isNext())
        {
            auto q = Np.get();
 
            if (p.getKey() == q)
            {
                ++Np;
                continue;
            }
 
            // repulsive
 
            Point<3,float> xq = vd.getPos(q);
            Point<3,float> f = (xp - xq);
 
            float distance = f.norm();
 
            // Particle should be inside r_cut range
 
            if (distance < r_cut )
            {
                vd.getProp<nn_sym>(p).add();
                vd.getProp<nn_sym>(q).add();
 
                vd.getProp<nn_sym>(p).last().xq = xq;
                vd.getProp<nn_sym>(q).last().xq = xp;
                vd.getProp<nn_sym>(p).last().id = vd.getProp<0>(q);
                vd.getProp<nn_sym>(q).last().id = vd.getProp<0>(p);
            }
 
            ++Np;
        }
 
        ++p_it2;
    }
 
    vd.ghost_put<merge_,2>();
 
 
 
    auto p_it3 = vd.getDomainIterator();
 
    bool ret = true;
    while (p_it3.isNext())
    {
        auto p = p_it3.get();
 
        vd.getProp<nn_norm>(p).sort();
        vd.getProp<nn_sym>(p).sort();
 
        ret &= vd.getProp<nn_norm>(p).size() == vd.getProp<nn_sym>(p).size();
 
        for (size_t i = 0 ; i < vd.getProp<1>(p).size() ; i++)
        {
            ret &= vd.getProp<nn_norm>(p).get(i).id == vd.getProp<nn_sym>(p).get(i).id;
 
            if (box.isInside(vd.getProp<nn_norm>(p).get(i).xq) == true)
            {
                ret &= vd.getProp<nn_norm>(p).get(i).xq == vd.getProp<nn_sym>(p).get(i).xq;
            }
        }
 
        ++p_it3;
    }
 
    if (ret != true)
    {
        std::cout << "ERROR" << std::endl;
        exit(1);
    }
 
 
 
    openfpm_finalize();
 
 
}