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