VerletList_test.hpp

/*
 * VerletList_test.hpp
 *
 *  Created on: Aug 16, 2016
 *      Author: i-bird
 */

#ifndef OPENFPM_DATA_SRC_NN_VERLETLIST_VERLETLIST_TEST_HPP_
#define OPENFPM_DATA_SRC_NN_VERLETLIST_VERLETLIST_TEST_HPP_

#include "NN/VerletList/VerletList.hpp"
#include "NN/VerletList/VerletListM.hpp"

template<unsigned int dim, typename T> void create_particles_on_grid(grid_sm<dim,void> & g_info, SpaceBox<dim,T> & box, openfpm::vector<Point<dim,T>> & v)
{
    // Create a grid iterator
    grid_key_dx_iterator<dim> g_it(g_info);

    // Usefull definition of points
    Point<dim,T> end = box.getP2() - box.getP1();

    Point<dim,T> middle;
    Point<dim,T> spacing;

    for (size_t i = 0 ; i < dim ; i++)
    {
        middle.get(i) = end.get(i) / g_info.size(i) / 2.0;
        spacing.get(i) = end.get(i) / g_info.size(i);
    }

    while (g_it.isNext())
    {
        // Add 2 particles on each cell

        Point<dim,T> key = Point<dim,T>(g_it.get().toPoint());
        key = pmul(key,spacing) + box.getP1() + middle;

        v.add(key);

        ++g_it;
    }
}

template<unsigned int dim, typename T> void create_particles_on_gridM(grid_sm<dim,void> & g_info,
                                                                      T r_cut,
                                                                      SpaceBox<dim,T> & box,
                                                                      openfpm::vector<pos_v<openfpm::vector<Point<dim,T>>>> & v,
                                                                      CellListM<dim,T,2> & cl)
{
    // Create a grid iterator
    grid_key_dx_iterator<dim> g_it(g_info);

    // Usefull definition of points
    Point<dim,T> end = box.getP2() - box.getP1();

    Point<dim,T> middle;
    Point<dim,T> spacing;

    Point<dim,T> offset[dim];

    size_t div[dim];
    for (size_t i = 0 ; i < dim ; i++)
    {
        offset[i] = middle;
        div[i] = 20 - 1;
    }

    for (size_t i = 0 ; i < dim ; i++)
        offset[i].get(i) += (1.0 / g_info.size(i)) / 32.0;

    for (size_t i = 0 ; i < dim ; i++)
    {
        middle.get(i) = end.get(i) / g_info.size(i) / 2.0;
        spacing.get(i) = end.get(i) / g_info.size(i);
    }

    cl.Initialize(box,div);

    while (g_it.isNext())
    {
        // Add 2 particles on each cell

        Point<dim,T> key = Point<dim,T>(g_it.get().toPoint());
        key = pmul(key,spacing) + box.getP1() + middle;

        v.get(0).pos.add(Point<dim,T>(key + offset[0]));
        v.get(1).pos.add(Point<dim,T>(key + offset[1]));
        v.get(2).pos.add(Point<dim,T>(key + offset[2]));

        cl.add(Point<dim,T>(key + offset[0]),v.get(0).pos.size()-1,0);
        cl.add(Point<dim,T>(key + offset[1]),v.get(1).pos.size()-1,1);
        cl.add(Point<dim,T>(key + offset[2]),v.get(2).pos.size()-1,2);

        ++g_it;
    }
}

template<unsigned int dim, typename T, typename VerS> void Verlet_list_s(SpaceBox<dim,T> & box)
{
    T r_cut = 1.506 * (box.getHigh(0) - box.getLow(0)) / 20;

    for (size_t k = 0 ; k < 5; k++)
    {
        // id Cell list

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

        grid_sm<3,void> ginfo(div);

        openfpm::vector<Point<dim,T>> pos;
        create_particles_on_grid(ginfo,box,pos);

        VerS vl1;

        Box<dim,T> innerBox;
        for (size_t i = 0 ; i < dim ; i++)
        {
            innerBox.setLow(i,r_cut);
            innerBox.setHigh(i,1.0-r_cut);
        }

        vl1.Initialize(box,box,r_cut,pos,pos.size());

        // Check that the verlet is consistent
        for (size_t i = 0 ; i < pos.size() ; i++)
        {
            if (innerBox.isInside(pos.get(i)) == false)
                continue;

            Point<dim,T> p = pos.get(i);

            if (k == 0)
            {BOOST_REQUIRE_EQUAL(vl1.getNNPart(i),19ul);}

            auto NN = vl1.getNNIterator(i);

            bool ret = true;

            while (NN.isNext())
            {
                auto k = NN.get();

                T dist = p.distance(Point<dim,T>(pos.get(k)));

                ret &= (dist < r_cut);

                ++NN;
            }

            BOOST_REQUIRE_EQUAL(ret,true);
        }

        r_cut += 0.075;
    }

    r_cut = 1.506 * (box.getHigh(0) - box.getLow(0)) / 20;

    size_t div[dim];
    for (size_t i = 0 ; i < dim ; i++)
        div[i] = 20-1;

    grid_sm<3,void> ginfo(div);

    openfpm::vector<Point<dim,T>> pos;
    create_particles_on_grid(ginfo,box,pos);


    // Initialize an external cell-list
    CellList<dim,T,Mem_fast<>,shift<3,double>> cli;
    Box<dim,T> bt = box;

    // Calculate the divisions for the Cell-lists
    cl_param_calculate(bt,div,r_cut,Ghost<dim,T>(0.0));

    // Initialize a cell-list
    cli.Initialize(bt,div,5);

    for (size_t i = 0; i < pos.size() ; i++)
        cli.add(pos.get(i), i);


    // Create cell-list

    for (size_t k = 0 ; k < 4; k++)
    {

        VerS vl1;
        vl1.Initialize(cli,r_cut,pos,pos,pos.size());



        VerS vl2;
        vl2.Initialize(box,box,r_cut,pos,pos.size());


        Box<dim,T> innerBox;
        for (size_t i = 0 ; i < dim ; i++)
        {
            innerBox.setLow(i,r_cut);
            innerBox.setHigh(i,1.0-r_cut);
        }

        // Check that the verlet is consistent
        for (size_t i = 0 ; i < pos.size() ; i++)
        {
            if (innerBox.isInside(pos.get(i)) == false)
                continue;

            Point<dim,T> p = pos.get(i);

            if (k == 0)
            {BOOST_REQUIRE_EQUAL(vl1.getNNPart(i),19ul);}
            BOOST_REQUIRE_EQUAL(vl1.getNNPart(i),vl2.getNNPart(i));

            openfpm::vector<size_t> v1;
            openfpm::vector<size_t> v2;


            bool ret = true;
            auto NN = vl1.getNNIterator(i);
            while (NN.isNext())
            {
                auto k = NN.get();

                T dist = p.distance(Point<dim,T>(pos.get(k)));

                ret &= (dist < r_cut);

                v1.add(k);

                ++NN;
            }


            auto NN2 = vl2.getNNIterator(i);
            while (NN2.isNext())
            {
                auto k = NN2.get();

                v2.add(k);

                ++NN2;
            }

            ret = true;
            v1.sort();
            v2.sort();
            for (size_t i = 0 ; i < v1.size(); i++)
            {
                ret &= v1.get(i) == v2.get(i);
            }

            BOOST_REQUIRE_EQUAL(ret,true);
        }

        r_cut += 0.075;
    }
}


template<unsigned int dim, typename T, typename VerS> void Verlet_list_sM(SpaceBox<dim,T> & box)
{
    T r_cut = (box.getHigh(0) - box.getLow(0)) / 20;

    for (size_t k = 0 ; k <= 0 ; k++)
    {
        // id Cell list

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

        grid_sm<3,void> ginfo(div);

        openfpm::vector<Point<dim,T>> ps1;
        openfpm::vector<Point<dim,T>> ps2;
        openfpm::vector<Point<dim,T>> ps3;

        openfpm::vector<pos_v<openfpm::vector<Point<dim,T>>>> pos2;
        pos2.add(pos_v<openfpm::vector<Point<dim,T>>>(ps1));
        pos2.add(pos_v<openfpm::vector<Point<dim,T>>>(ps2));
        pos2.add(pos_v<openfpm::vector<Point<dim,T>>>(ps3));

        CellListM<dim,T,2> cl;
        create_particles_on_gridM(ginfo,r_cut,box,pos2,cl);

        VerS vl1;

        Box<dim,T> innerBox;
        for (size_t i = 0 ; i < dim ; i++)
        {
            innerBox.setLow(i,r_cut);
            innerBox.setHigh(i,1.0-r_cut);
        }

        vl1.Initialize(cl,0,r_cut,pos2.get(0).pos,pos2,pos2.get(0).pos.size());

        // Check that the verlet is consistent
        for (size_t i = 0 ; i < pos2.get(0).pos.size() ; i++)
        {
            if (innerBox.isInside(pos2.get(0).pos.get(i)) == false)
                continue;

            Point<dim,T> p = pos2.get(0).pos.get(i);

            if (k == 0)
            {BOOST_REQUIRE_EQUAL(vl1.getNNPart(i),57ul);}

            auto NN = vl1.getNNIterator(i);

            bool ret = true;

            while (NN.isNext())
            {
                auto k = NN.getP();
                auto v = NN.getV();

                T dist = p.distance(Point<dim,T>(pos2.get(v).pos.get(k)));

                ret &= (dist < r_cut);

                ++NN;
            }

            BOOST_REQUIRE_EQUAL(ret,true);
        }

        r_cut += 0.075;
    }

    r_cut = (box.getHigh(0) - box.getLow(0)) / 20;

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

    grid_sm<3,void> ginfo(div);

    // Initialize an external cell-list
    CellListM<dim,T,2> cli;


    openfpm::vector<Point<dim,T>> ps1;
    openfpm::vector<Point<dim,T>> ps2;
    openfpm::vector<Point<dim,T>> ps3;

    openfpm::vector<pos_v<openfpm::vector<Point<dim,T>>>> pos2;
    pos2.add(pos_v<openfpm::vector<Point<dim,T>>>(ps1));
    pos2.add(pos_v<openfpm::vector<Point<dim,T>>>(ps2));
    pos2.add(pos_v<openfpm::vector<Point<dim,T>>>(ps3));

    create_particles_on_gridM(ginfo,r_cut,box,pos2,cli);



    // Create cell-list

    for (size_t k = 0 ; k <= 0 ; k++)
    {

        VerS vl1;
        vl1.Initialize(cli,0,r_cut,pos2.get(0).pos,pos2,pos2.get(0).pos.size());



        VerS vl2;
        vl2.Initialize(cli,0,r_cut,pos2.get(0).pos,pos2,pos2.get(0).pos.size());


        Box<dim,T> innerBox;
        for (size_t i = 0 ; i < dim ; i++)
        {
            innerBox.setLow(i,r_cut);
            innerBox.setHigh(i,1.0-r_cut);
        }

        // Check that the verlet is consistent
        for (size_t i = 0 ; i < pos2.get(0).pos.size() ; i++)
        {
            if (innerBox.isInside(pos2.get(0).pos.get(i)) == false)
                continue;

            Point<dim,T> p = pos2.get(0).pos.get(i);

            if (k == 0)
            {BOOST_REQUIRE_EQUAL(vl1.getNNPart(i),57ul);}
            BOOST_REQUIRE_EQUAL(vl1.getNNPart(i),vl2.getNNPart(i));

            openfpm::vector<size_t> v1;
            openfpm::vector<size_t> v2;


            bool ret = true;
            auto NN = vl1.getNNIterator(i);
            while (NN.isNext())
            {
                auto k = NN.getP();
                auto v = NN.getV();

                T dist = p.distance(Point<dim,T>(pos2.get(v).pos.get(k)));

                ret &= (dist < r_cut);

                v1.add(NN.get());

                ++NN;
            }


            auto NN2 = vl2.getNNIterator(i);
            while (NN2.isNext())
            {
                auto k = NN2.get();

                v2.add(k);

                ++NN2;
            }

            ret = true;
            v1.sort();
            v2.sort();
            for (size_t i = 0 ; i < v1.size(); i++)
            {
                ret &= v1.get(i) == v2.get(i);
            }

            BOOST_REQUIRE_EQUAL(ret,true);
        }

        r_cut += 0.075;
    }
}


BOOST_AUTO_TEST_SUITE( VerletList_test )

BOOST_AUTO_TEST_CASE( VerletList_use)
{
    std::cout << "Test verlet list" << "\n";

    SpaceBox<3,double> box({0.0f,0.0f,0.0f},{1.0f,1.0f,1.0f});
    SpaceBox<3,double> box2({-1.0f,-1.0f,-1.0f},{1.0f,1.0f,1.0f});
//  Verlet_list_s<3,double,VerletList<3,double,FAST,shift<3,double>>>(box);
//  Verlet_list_s<3,double,VerletList<3,double,FAST,shift<3,double>> >(box2);
    Verlet_list_sM<3,double,VerletListM<3,double,2>>(box);
//  Verlet_list_sM<3,double,CellListM<3,double,8>>(box2);

    std::cout << "End verlet list" << "\n";

    // Test the cell list
}

BOOST_AUTO_TEST_SUITE_END()



#endif /* OPENFPM_DATA_SRC_NN_VERLETLIST_VERLETLIST_TEST_HPP_ */