main.cpp

#include "Vector/vector_dist.hpp"
#include "Decomposition/CartDecomposition.hpp"
#include "data_type/aggregate.hpp"
#include "NN/CellList/CellListM.hpp"
#include "Vector/vector_dist_multiphase_functions.hpp"

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

    openfpm_init(&argc,&argv);

    // Vcluster for general usage
    Vcluster & v_cl = create_vcluster();

    // 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};

    // cut-off radius
    float r_cut = 0.05;

    // ghost, big enough to contain the interaction radius
    Ghost<3,float> ghost(r_cut);

    // The set of phases
    openfpm::vector< vector_dist<3,float, aggregate<double,double>> > phases;

    // first phase
    phases.add( vector_dist<3,float, aggregate<double,double>>(4096,box,bc,ghost) );

    // The other 3 phases
    phases.add( vector_dist<3,float, aggregate<double,double>>(phases.get(0).getDecomposition(),4096) );
    phases.add( vector_dist<3,float, aggregate<double,double>>(phases.get(0).getDecomposition(),4096) );
    phases.add( vector_dist<3,float, aggregate<double,double>>(phases.get(0).getDecomposition(),4096) );





    // An iterator over the particles of the phase0
    auto it = phases.get(0).getDomainIterator();

    // For all the particles of phase0
    while (it.isNext())
    {
        // particle p
        auto p = it.get();

        // Assign the position of the particles to each phase
        for (size_t i = 0 ; i < phases.size(); i++)
        {
            // Assign random position
            phases.get(i).getPos(p)[0] = (float)rand() / RAND_MAX;
            phases.get(i).getPos(p)[1] = (float)rand() / RAND_MAX;
            phases.get(i).getPos(p)[2] = (float)rand() / RAND_MAX;
        }

        // Next particle
        ++it;
    }


    // Redistribute and sync all the phases
    for (size_t i = 0 ; i < 4 ; i++)
    {
        phases.get(i).map();
        phases.get(i).ghost_get<>();
    }



    {

    // Get the cell list of the phase1
    auto CL_phase1 = phases.get(1).getCellList(r_cut);

    // This function create a Verlet-list between phases 0 and 1
    auto NN_ver01 = createVerlet(phases.get(0),phases.get(1),CL_phase1,r_cut);



    // Get an iterator of the particles of the phase0
    it = phases.get(0).getDomainIterator();

    // For each particle of the phase0
    while (it.isNext())
    {
        // Get the particle p
        auto p = it.get();

        // reset the counter
        phases.get(0).getProp<0>(p) = 0;

        // Get an iterator over the neighborhood particles for the particle p
        auto Np = NN_ver01.template getNNIterator<NO_CHECK>(p.getKey());

        // For each neighborhood of the particle p
        while (Np.isNext())
        {
            // Neighborhood particle q
            auto q = Np.get();

            // Count the number of particles
            // Here in general we can do our computation
            phases.get(0).getProp<0>(p)++;

            // Next particle
            ++Np;
        }

        // Next particle
        ++it;
    }

    }



    // This function create an "Empty" Multiphase Cell List from all the phases

    // We just create a cell list with slightly bigget r_cut to avoid error in Verlet creation
    // because of round off errors
    float r_cut2 = r_cut*1.00001;
    auto CL_all = createCellListM<2>(phases,r_cut2);

    // This create a Verlet-list between phase0 and all the other phases
    auto NNver0_all = createVerletM<2>(0,phases.get(0),phases,CL_all,r_cut);



    // Get an iterator for the phase0 particles
    it = phases.get(0).getDomainIterator();

    while (it.isNext())
    {
        // Get the particle p
        auto p = it.get();

        // Get an interator over the neighborhood of the particle p
        auto Np = NNver0_all.template getNNIterator<NO_CHECK>(p.getKey());

        // reset the counter
        phases.get(0).getProp<0>(p) = 0;

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

            // Get from which phase it come from
            auto ph_q = Np.getV();

            // count
            phases.get(0).getProp<0>(p)++;

            // Next particle
            ++Np;
        }

        // Next particle
        ++it;
    }


    {

    // Get the cell list of the phase1
    auto CL_phase1 = phases.get(1).getCellListSym(r_cut);

    // This function create a Verlet-list between phases 0 and 1
    auto NN_ver01 = createVerletSym(phases.get(0),phases.get(1),CL_phase1,r_cut);

    // Get an iterator over the real and ghost particles
    it = phases.get(0).getDomainAndGhostIterator();

    // For each particles
    while (it.isNext())
    {
        // Get the particle p
        auto p = it.get();

        // Reset the counter
        phases.get(0).getProp<0>(p) = 0;

        // Next particle
        ++it;
    }

    // Compute interaction from phase0 to phase1

    // Get an iterator over the real particles of phase0
    it = phases.get(0).getDomainIterator();

    // For each particle of the phase0
    while (it.isNext())
    {
        // get particle p
        auto p = it.get();

        // Get the neighborhood of particle p
        auto Np = NN_ver01.template getNNIterator<NO_CHECK>(p.getKey());

        // For each neighborhood of the particle p
        while (Np.isNext())
        {
            // Neighborhood particle q of p
            auto q = Np.get();

            // increment the counter for the phase 0 and 1
            phases.get(0).getProp<0>(p)++;
            phases.get(1).getProp<0>(q)++;

            // Next particle
            ++Np;
        }

        // Next particle
        ++it;
    }

    // Merge the information of the ghost with the real particles
    phases.get(0).ghost_put<add_,0>();
    phases.get(1).ghost_put<add_,0>();

    }



    // Get an iterator over the phase0
    it = phases.get(0).getDomainAndGhostIterator();

    // For each particle of the phase 0
    while (it.isNext())
    {
        // get the particle p
        auto p = it.get();

        // Reset the counter for the domain and ghost particles
        phases.get(0).getProp<0>(p) = 0;
        phases.get(1).getProp<0>(p) = 0;
        phases.get(2).getProp<0>(p) = 0;
        phases.get(3).getProp<0>(p) = 0;

        // next particle
        ++it;
    }

    // This function create an "Empty" Multiphase Cell List
    CL_all = createCellListSymM<2>(phases,r_cut);

    // Type of the multiphase Verlet-list
    typedef decltype(createVerletSymM<2>(0,phases.get(0),phases,CL_all,r_cut)) verlet_type;

    // for each phase we create one Verlet-list that contain the neighborhood
    // from all the phases
    verlet_type NNver_all[4];

    // Here we create a Verlet-list between each phases

    // 0 to all
    NNver_all[0] = createVerletSymM<2>(0,phases.get(0),phases,CL_all,r_cut);

    // 1 to all
    NNver_all[1] = createVerletSymM<2>(1,phases.get(1),phases,CL_all,r_cut);

    // 2 to all
    NNver_all[2] = createVerletSymM<2>(2,phases.get(2),phases,CL_all,r_cut);

    // 3 to all
    NNver_all[3] = createVerletSymM<2>(3,phases.get(3),phases,CL_all,r_cut);

    // For each phase
    for (size_t i = 0 ; i < phases.size() ; i++)
    {
        // Get an iterator over the particles of that phase
        it = phases.get(i).getDomainIterator();

        // for each particle of the phase
        while (it.isNext())
        {
            // Get the particle p
            auto p = it.get();

            // Get an iterator for neighborhood of the particle p
            auto Np = NNver_all[i].template getNNIterator<NO_CHECK>(p.getKey());

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

                // Get from which phase it come from
                auto ph_q = Np.getV();

                // increment the counter on both p and q
                phases.get(i).getProp<0>(p)++;
                phases.get(ph_q).getProp<0>(q)++;

                // Next particle
                ++Np;
            }

            // Next particle
            ++it;
        }
    }

    // Merge the ghost part with the real particles
    phases.get(0).ghost_put<add_,0>();
    phases.get(1).ghost_put<add_,0>();
    phases.get(2).ghost_put<add_,0>();
    phases.get(3).ghost_put<add_,0>();



    // we create a reference to phase0 particle for convenience
    vector_dist<3,float, aggregate<double,double> > & current_phase = phases.get(0);

    // Get the iterator of the particles of phase 0
    auto it2 = current_phase.getIterator();

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

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

        // Reset to zero the propety 0
        current_phase.getProp<0>(p) = 0.0;

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

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

            // Get from which phase it come from
            auto ph_q = Np.getV();

            Point<3,float> xq = phases.get(ph_q).getPos(q);

            // we accumulate all the distances
            current_phase.getProp<0>(p) = norm(xp - xq);

            ++Np;
        }

        // Next particle p
        ++it2;
    }



    openfpm_finalize();


}