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();


}