vector_dist_performance_common.hpp

/*
 * vector_dist_performance_common.hpp
 *
 *  Created on: Dec 25, 2016
 *      Author: i-bird
 */

#ifndef SRC_VECTOR_PERFORMANCE_VECTOR_DIST_PERFORMANCE_COMMON_HPP_
#define SRC_VECTOR_PERFORMANCE_VECTOR_DIST_PERFORMANCE_COMMON_HPP_

#include "Vector/vector_dist.hpp"

template<unsigned int dim, size_t prp = 0, typename T, typename V> void calc_forces(T & NN, V & vd, float r_cut)
{
    auto it_v = vd.getDomainIterator();

    float sum[dim];

    for (size_t i = 0; i < dim; i++)
        sum[i] = 0;

    while (it_v.isNext())
    {
        //key
        vect_dist_key_dx key = it_v.get();

        // Get the position of the particles
        Point<dim,float> p = vd.getPos(key);

        for (size_t i = 0; i < dim; i++)
            sum[i] = 0;

        // Get the neighborhood of the particle
        auto cell_it = NN.template getNNIterator<NO_CHECK>(NN.getCell(p));

        while(cell_it.isNext())
        {
            auto nnp = cell_it.get();

            // p != q
            if (nnp == key.getKey())
            {
                ++cell_it;
                continue;
            }

            Point<dim,float> q = vd.getPos(nnp);

            if (p.distance2(q) < r_cut*r_cut)
            {
                //Calculate the forces
                float num[dim];
                for (size_t i = 0; i < dim; i++)
                    num[i] = vd.getPos(key)[i] - vd.getPos(nnp)[i];

                float denom = 0;
                for (size_t i = 0; i < dim; i++)
                    denom += num[i] * num[i];

                float res[dim];
                for (size_t i = 0; i < dim; i++)
                    res[i] = num[i] / denom;

                for (size_t i = 0; i < dim; i++)
                    sum[i] += res[i];
            }
            //Next particle in a cell
            ++cell_it;
        }

        //Put the forces
        for (size_t i = 0; i < dim; i++)
            vd.template getProp<prp>(key)[i] += sum[i];

        //Next particle in cell list
        ++it_v;
    }
}

template<unsigned int dim, unsigned int prp, typename T, typename V> void cross_calc(T & NN, T & NN2, V & vd, V & vd2)
{
    auto it_v = vd.getDomainIterator();

    while (it_v.isNext())
    {
        //key
        vect_dist_key_dx key = it_v.get();

        // Get the position of the particles
        Point<dim,float> p = vd.getPos(key);

        // Get the neighborhood of the particle
        auto cell_it = NN2.template getNNIterator<NO_CHECK>(NN2.getCell(p));

        double sum = 0.0;

        while(cell_it.isNext())
        {
            auto nnp = cell_it.get();

            Point<dim,float> q = vd2.getPos(nnp);

            sum += norm(p - q);

            //Next particle in a cell
            ++cell_it;
        }

        vd.template getProp<prp>(key) = sum;

        //Next particle in cell list
        ++it_v;
    }
}

template<unsigned int dim, typename v_dist>
void vd_initialize(v_dist & vd, Vcluster<HeapMemory> & v_cl)
{
    // The random generator engine
    std::default_random_engine eg(v_cl.getProcessUnitID()*4313);
    std::uniform_real_distribution<float> ud(0.0f, 1.0f);


    auto it = vd.getIterator();

    while (it.isNext())
    {
        auto key = it.get();

        for (size_t i = 0; i < dim; i++)
        {vd.getPos(key)[i] = ud(eg);}

        ++it;
    }

    vd.map();
}

template<unsigned int dim, typename v_dist>
void vd_initialize_box_nomap(v_dist & vd, const Box<dim,typename v_dist::stype> & box, Vcluster<HeapMemory> & v_cl,int start, int stop)
{
    // The random generator engine
    std::default_random_engine eg(v_cl.getProcessUnitID()*4313);
    std::uniform_real_distribution<float> ud(0.0f, 1.0f);

    auto it = vd.getIterator(start,stop);

    while (it.isNext())
    {
        auto key = it.get();

        for (size_t i = 0; i < dim; i++)
        {vd.getPos(key)[i] = (box.getHigh(i) - box.getLow(i))*ud(eg) + box.getLow(i);}

        ++it;
    }
}


template<unsigned int dim, typename v_dist> void vd_initialize_double(v_dist & vd,v_dist & vd2, Vcluster<> & v_cl, size_t k_int)
{
    // The random generator engine
    std::default_random_engine eg(v_cl.getProcessUnitID()*4313);
    std::uniform_real_distribution<float> ud(0.0f, 1.0f);


    auto it = vd.getIterator();

    while (it.isNext())
    {
        auto key = it.get();

        for (size_t i = 0; i < dim; i++)
        {
            vd.getPos(key)[i] = ud(eg);
            vd2.getPos(key)[i] = vd.getPos(key)[i];
        }

        ++it;
    }

    vd.map();
    vd2.map();
}



template<unsigned int dim, size_t prp = 0, typename T, typename V> void calc_forces_hilb(T & NN, V & vd, float r_cut)
{
    auto it_cl = NN.getIterator();

    float sum[dim];

    for (size_t i = 0; i < dim; i++)
        sum[i] = 0;

    while (it_cl.isNext())
    {
        //key
        auto key = it_cl.get();

        // Get the position of the particles
        Point<dim,float> p = vd.getPos(key);

        for (size_t i = 0; i < dim; i++)
            sum[i] = 0;

        // Get the neighborhood of the particle
        auto cell_it = NN.template getNNIterator<NO_CHECK>(NN.getCell(p));

        while(cell_it.isNext())
        {
            auto nnp = cell_it.get();

            // p != q
            if (nnp == key)
            {
                ++cell_it;
                continue;
            }

            Point<dim,float> q = vd.getPos(nnp);

            if (p.distance2(q) < r_cut*r_cut)
            {
                //Calculate the forces
                float num[dim];
                for (size_t i = 0; i < dim; i++)
                    num[i] = vd.getPos(key)[i] - vd.getPos(nnp)[i];

                float denom = 0;
                for (size_t i = 0; i < dim; i++)
                    denom += num[i] * num[i];

                float res[dim];
                for (size_t i = 0; i < dim; i++)
                    res[i] = num[i] / denom;

                for (size_t i = 0; i < dim; i++)
                    sum[i] += res[i];
            }
            //Next particle in a cell
            ++cell_it;
        }

        //Put the forces
        for (size_t i = 0; i < dim; i++)
            vd.template getProp<prp>(key)[i] += sum[i];

        //Next particle in cell list
        ++it_cl;
    }
}

#endif /* SRC_VECTOR_PERFORMANCE_VECTOR_DIST_PERFORMANCE_COMMON_HPP_ */