main.cpp

#include "Vector/vector_dist.hpp"
#include "Decomposition/CartDecomposition.hpp"
#include "data_type/aggregate.hpp"

int main(int argc, char* argv[])
{

    openfpm_init(&argc,&argv);
    Vcluster & v_cl = create_vcluster();

    // we will place the particles on a grid like way with 128 particles on each direction
    size_t sz[3] = {128,128,128};

    // The domain
    Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});

    // Boundary conditions
    size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};

    // ghost, big enough to contain the interaction radius
    Ghost<3,float> ghost(1.0/(128-2));



    vector_dist<3,float, aggregate<double,double[3],double[3][3]> > vd(0,box,bc,ghost);



    auto it = vd.getGridIterator(sz);

    // For all the particles
    while (it.isNext())
    {
        vd.add();

        auto key = it.get();

        // The position of the particle is given by the point-index (0,0,0) ; (0,0,1) ... (127,127,127)
        // multiplied by the spacing
        vd.getLastPos()[0] = key.get(0) * it.getSpacing(0);
        vd.getLastPos()[1] = key.get(1) * it.getSpacing(1);
        vd.getLastPos()[2] = key.get(2) * it.getSpacing(2);

        // next point
        ++it;
    }



    vd.ghost_get<0,1,2>();



    float r_cut = 1.0/(128-2);

    // Get cell list
    auto NN = vd.getCellList(r_cut);

    // Get the iterator
    auto it2 = vd.getDomainIterator();

    // For each particle ...
    while (it2.isNext())
    {
        // ... p
        auto p = it2.get();

        // Get the position of the particle p
        Point<3,float> xp = vd.getPos(p);

        // Get an iterator of all the particles neighborhood of p
        auto Np = NN.getNNIterator(NN.getCell(vd.getPos(p)));

        // For each particle near p
        while (Np.isNext())
        {
            // Get the particle q near to p
            auto q = Np.get();

            // Get the position of the particle q
            Point<3,float> xq = vd.getPos(q);

            // Calculate the distance vector between p and q
            Point<3,float> f = (xp - xq);

            // we sum the distance of all the particles
            vd.template getProp<0>(p) += f.norm();;

            // we sum the distance of all the particles
            vd.template getProp<1>(p)[0] += f.get(0);
            vd.template getProp<1>(p)[1] += f.get(1);
            vd.template getProp<1>(p)[2] += f.get(2);

            // Moment of inertia
            vd.template getProp<2>(p)[0][0] += (xp.get(0) - xq.get(0)) * (xp.get(0) - xq.get(0));
            vd.template getProp<2>(p)[0][1] += (xp.get(0) - xq.get(0)) * (xp.get(1) - xq.get(1));
            vd.template getProp<2>(p)[0][2] += (xp.get(0) - xq.get(0)) * (xp.get(2) - xq.get(2));
            vd.template getProp<2>(p)[1][0] += (xp.get(1) - xq.get(1)) * (xp.get(0) - xq.get(0));
            vd.template getProp<2>(p)[1][1] += (xp.get(1) - xq.get(1)) * (xp.get(1) - xq.get(1));
            vd.template getProp<2>(p)[1][2] += (xp.get(1) - xq.get(1)) * (xp.get(2) - xq.get(2));
            vd.template getProp<2>(p)[2][0] += (xp.get(2) - xq.get(2)) * (xp.get(0) - xq.get(0));
            vd.template getProp<2>(p)[2][1] += (xp.get(2) - xq.get(2)) * (xp.get(1) - xq.get(1));
            vd.template getProp<2>(p)[2][2] += (xp.get(2) - xq.get(2)) * (xp.get(2) - xq.get(2));

            ++Np;
        }

        // Next particle p
        ++it2;
    }



    auto verlet = vd.getVerlet(r_cut);

    auto it3 = vd.getDomainIterator();

    // For each particle i verlet.size() == Number of particles
    while (it3.isNext())
    {
        auto p = it3.get();

        // get the position of the particle i
        Point<3,float> xp = vd.getPos(p);

        // get the Neighborhood of p
        auto NNp = verlet.getNNIterator(p.getKey());

        // for each neighborhood j of particle to i
        while (NNp.isNext())
        {
            size_t q = NNp.get();

            // Get the position of the particle neighborhood if i
            Point<3,float> xq = vd.getPos(q);

            // Calculate the distance vector between p and q
            Point<3,float> f = (xp - xq);

            // we sum the distance of all the particles
            vd.template getProp<0>(p) += f.norm();;

            // we sum the distance of all the particles
            vd.template getProp<1>(p)[0] += f.get(0);
            vd.template getProp<1>(p)[1] += f.get(1);
            vd.template getProp<1>(p)[2] += f.get(2);

            // Moment of inertia
            vd.template getProp<2>(p)[0][0] += (xp.get(0) - xq.get(0)) * (xp.get(0) - xq.get(0));
            vd.template getProp<2>(p)[0][1] += (xp.get(0) - xq.get(0)) * (xp.get(1) - xq.get(1));
            vd.template getProp<2>(p)[0][2] += (xp.get(0) - xq.get(0)) * (xp.get(2) - xq.get(2));
            vd.template getProp<2>(p)[1][0] += (xp.get(1) - xq.get(1)) * (xp.get(0) - xq.get(0));
            vd.template getProp<2>(p)[1][1] += (xp.get(1) - xq.get(1)) * (xp.get(1) - xq.get(1));
            vd.template getProp<2>(p)[1][2] += (xp.get(1) - xq.get(1)) * (xp.get(2) - xq.get(2));
            vd.template getProp<2>(p)[2][0] += (xp.get(2) - xq.get(2)) * (xp.get(0) - xq.get(0));
            vd.template getProp<2>(p)[2][1] += (xp.get(2) - xq.get(2)) * (xp.get(1) - xq.get(1));
            vd.template getProp<2>(p)[2][2] += (xp.get(2) - xq.get(2)) * (xp.get(2) - xq.get(2));

            ++NNp;
        }

        ++it3;
    }



    // Get cell list
    auto NN4 = vd.getCellList<CELL_MEMBAL(3,float)>(r_cut);
    auto NN5 = vd.getCellList<CELL_MEMMW(3,float)>(r_cut);



    openfpm_finalize();


}