CellList_test.hpp

/*
 * CellList_test.hpp
 *
 *  Created on: Mar 23, 2015
 *      Author: Pietro Incardona
 */

#include "CellList.hpp"
#include "CellListM.hpp"
#include "Grid/grid_sm.hpp"

#ifndef CELLLIST_TEST_HPP_
#define CELLLIST_TEST_HPP_


template<unsigned int dim, typename T, typename CellS> void Test_cell_s(SpaceBox<dim,T> & box)
{
    //Space where is living the Cell list
    //SpaceBox<dim,T> box({0.0f,0.0f,0.0f},{1.0f,1.0f,1.0f});

    // Subdivisions
    size_t div[dim] = {16,16,16};

    // Origin
    Point<dim,T> org({0.0,0.0,0.0});

    // id Cell list
    CellS cl2(box,div);

    // grid info
    grid_sm<dim,void> g_info(div);

    // Test force reallocation in case of Cell list fast
    for (size_t i = 0 ; i < CELL_REALLOC * 3 ; i++)
    {
        cl2.add(org,i);
    }

    // Check the elements
    BOOST_REQUIRE_EQUAL(cl2.getNelements(cl2.getCell(org)),CELL_REALLOC * 3ul);
    for (size_t i = 0 ; i < CELL_REALLOC * 3 ; i++)
    {
        BOOST_REQUIRE_EQUAL(cl2.get(cl2.getCell(org),i),i);
    }


    // id Cell list
    CellS cl1(box,div);

    // Create a grid iterator
    grid_key_dx_iterator<dim> g_it(g_info);

    // Iterate through each element
    // Add 1 element for each cell

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

    Point<dim,T> offset[dim] = {middle,middle,middle};

    // Create offset shift vectors
    for (size_t i = 0 ; i < dim ; i++)
    {
        offset[i].get(i) += (1.0 / div[i]) / 8.0;
    }

    openfpm::vector<Point<dim,T>> pos;
    size_t id = 0;

    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) + offset[0] + box.getP1();
        pos.add(key);

        cl1.add(key,id);
        ++id;

        key = Point<dim,T>(g_it.get().toPoint());
        key = pmul(key,spacing) + offset[1] + box.getP1();
        pos.add(key);

        cl1.add(key,id);
        ++id;

        ++g_it;
    }


    // check the cell are correctly filled

    // reset iterator
    g_it.reset();

    while (g_it.isNext())
    {
        // Check that there are 2 particles on each cell

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

        size_t cell = cl1.getCell(key);
        size_t n_ele = cl1.getNelements(cell);

        BOOST_REQUIRE_EQUAL(n_ele,2ul);
        BOOST_REQUIRE_EQUAL((long int)(cl1.get(cell,1) - cl1.get(cell,0)),1);

        ++g_it;
    }

    // reset itarator
    g_it.reset();


    while (g_it.isNext())
    {
        // remove 1 particle on each cell

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

        auto cell = cl1.getCell(key);

        // Remove the first particle in the cell
        cl1.remove(cell,0);
        ++g_it;
    }


    // Check we have 1 object per cell
    g_it.reset();

    while (g_it.isNext())
    {
        // remove 1 particle on each cell

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

        auto cell = cl1.getCell(key);
        size_t n_ele = cl1.getNelements(cell);

        BOOST_REQUIRE_EQUAL(n_ele,1ul);
        ++g_it;
    }


    // Check we have 1 object per cell

    // Create a grid iterator
    grid_key_dx<dim> p1(1,1,1);
    grid_key_dx<dim> p2(div[0]-2,div[1]-2,div[2]-2);
    grid_key_dx_iterator_sub<dim> g_it_s(g_info,p1,p2);
    id = 0;

    while (g_it_s.isNext())
    {

        Point<dim,T> key = Point<dim,T>(g_it_s.get().toPoint());
        key = pmul(key,spacing) + offset[0] + box.getP1();

        auto NN = cl1.template getNNIterator<NO_CHECK>(cl1.getCell(key));
        size_t total = 0;

        while(NN.isNext())
        {
            // total

            total++;

            ++NN;
        }


        BOOST_REQUIRE_EQUAL(total,(size_t)openfpm::math::pow(3,dim));

        // in SE1_CLASS the cell list consider this construction as an hack
        // disable the test

#ifndef SE_CLASS1

        id = cl1.get(cl1.getCell(key),0);
        auto NNSym = cl1.template getNNIteratorSym<NO_CHECK>(cl1.getCell(key),id,pos);
        total = 0;

        while(NNSym.isNext())
        {
            // total

            total++;

            ++NNSym;
        }

        BOOST_REQUIRE_EQUAL(total,(size_t)openfpm::math::pow(3,dim) / 2 + 1);

#endif

        ++g_it_s;
    }

}


template<unsigned int dim, typename T, typename CellS> void Test_cell_sM(SpaceBox<dim,T> & box)
{
    //Space where is living the Cell list
    //SpaceBox<dim,T> box({0.0f,0.0f,0.0f},{1.0f,1.0f,1.0f});

    // Subdivisions
    size_t div[dim] = {16,16,16};

    // Origin
    Point<dim,T> org({0.0,0.0,0.0});

    // grid info
    grid_sm<dim,void> g_info(div);


    // CellS = CellListM<dim,T,8>
    CellS cl1(box,div);

    // Create a grid iterator
    grid_key_dx_iterator<dim> g_it(g_info);

    // Iterate through each element
    // Add 1 element for each cell

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

    Point<dim,T> offset[dim] = {middle,middle,middle};

    // Create offset shift vectors
    for (size_t i = 0 ; i < dim ; i++)
    {
        offset[i].get(i) += (1.0 / div[i]) / 8.0;
    }

    openfpm::vector<Point<dim,T>> phase1;
    openfpm::vector<Point<dim,T>> phase2;

    openfpm::vector<pos_v<openfpm::vector<Point<dim,T>>>> phases;
    phases.add(pos_v<openfpm::vector<Point<dim,T>>>(phase1));
    phases.add(pos_v<openfpm::vector<Point<dim,T>>>(phase2));

    size_t id = 0;

    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) + offset[0] + box.getP1();

        phase1.add(key);

        cl1.add(key,id,0);

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

        phase2.add(key);
        cl1.add(key,id,1);
        ++id;

        ++g_it;
    }

    // check the cell are correctly filled

    // reset iterator
    g_it.reset();

    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) + offset[2] + box.getP1();

        size_t cell = cl1.getCell(key);
        size_t n_ele = cl1.getNelements(cell);

        size_t p1 = cl1.getP(cell,1);
        size_t p2 = cl1.getP(cell,0);

        size_t v1 = cl1.getV(cell,1);
        size_t v2 = cl1.getV(cell,0);

        BOOST_REQUIRE_EQUAL(n_ele,2ul);
        BOOST_REQUIRE_EQUAL((long int)(p1 - p2),0);
        BOOST_REQUIRE_EQUAL((long int)(v1 - v2),1);
        ++g_it;
    }
}



template<typename CellList> void Test_CellDecomposer_consistent()
{
    Box<2,float> bx({-1.0/3.0,-1.0/3.0},{1.0/3.0,1.0/3.0});

    size_t div[2] = {36,36};

    CellDecomposer_sm<2,float,shift<2,float>> cd(bx,div,1);

    Box<2,float> bx_sub({-1.0/5.0,-1.0/5.0},{1.0/5.0,1.0/5.0});

    size_t bc[2] = {NON_PERIODIC,NON_PERIODIC};
    CellList cl(cd,bx_sub);
    Box<2,long int> bx_int = cd.convertDomainSpaceIntoGridUnits(bx_sub,bc);

    BOOST_REQUIRE_EQUAL(bx_int.getLow(0),8);
    BOOST_REQUIRE_EQUAL(bx_int.getLow(1),8);

    BOOST_REQUIRE_EQUAL(bx_int.getHigh(0),28);
    BOOST_REQUIRE_EQUAL(bx_int.getHigh(1),28);

    cd.convertCellUnitsIntoDomainSpace(bx_sub);

    BOOST_REQUIRE_EQUAL(bx_sub.getLow(0),-1.0f/5.0f);
    BOOST_REQUIRE_EQUAL(bx_sub.getLow(1),-1.0f/5.0f);

    BOOST_REQUIRE_EQUAL(bx_sub.getHigh(0),1.0f/5.0f);
    BOOST_REQUIRE_EQUAL(bx_sub.getHigh(1),1.0f/5.0f);
}

template<unsigned int dim, typename T, typename CellS> void Test_NN_iterator_radius(SpaceBox<dim,T> & box)
{
    //Space where is living the Cell list
    //SpaceBox<dim,T> box({0.0f,0.0f,0.0f},{1.0f,1.0f,1.0f});

    // Subdivisions
    size_t div[dim];
    size_t div2[dim];

    for (size_t i = 0 ; i < dim ; i++)
    {
        div[i] = 17;
        div2[i] = 34;
    }

    // grid info
    grid_sm<dim,void> g_info(div);
    grid_sm<dim,void> g_info2(div2);


    // CellS = CellListM<dim,T,8>
    CellS cl1(box,div,2);
    CellS cl2(box,div2,3);

    T radius = (box.getHigh(0) - box.getLow(0))/div[0];

    cl2.setRadius( radius );

    // create a vector of random point

    openfpm::vector<Point<dim,float>> vrp;

    for (size_t j = 0 ; j < 10000 ; j++)
    {
        vrp.add();

        for (size_t i = 0 ; i < dim ; i++)
        {
            vrp.template get<0>(j)[i] = ((float)rand() / (float)RAND_MAX)*(box.getHigh(i) - box.getLow(i)) + box.getLow(i);
        }
    }

    auto g_it = vrp.getIterator();

    while (g_it.isNext())
    {
        Point<dim,T> xp = vrp.get(g_it.get());

        size_t debug = cl1.getCell(xp);

        cl1.add(xp,g_it.get());
        cl2.add(xp,g_it.get());

        ++g_it;
    }

    // Get the neighborhood of each particle and compare the numbers of cell

    bool match = true;
    auto g_it2 = vrp.getIterator();

    size_t number_of_nn = 0;
    size_t number_of_nn2 = 0;

    while (g_it2.isNext())
    {
        Point<dim,T> xp = vrp.get(g_it2.get());

        openfpm::vector<size_t> ids1;
        openfpm::vector<size_t> ids2;

        size_t local1 = 0;
        size_t local2 = 0;

        auto NNit = cl1.getNNIterator(cl1.getCell(xp));

        while (NNit.isNext())
        {
            auto q = NNit.get();

            // calculate distance

            Point<dim,T> xq = vrp.get(q);
            Point<dim,T> r = xq - xp;

            if (r.norm() <= radius)
            {ids1.add(q);}

            number_of_nn++;
            local1++;

            ++NNit;
        }

        auto NN2it = cl2.getNNIteratorRadius(cl2.getCell(xp));

        while (NN2it.isNext())
        {
            auto q = NN2it.get();

            Point<dim,T> xq = vrp.get(q);
            Point<dim,T> r = xq - xp;

            if (r.norm() <= radius)
            {ids2.add(q);}

            number_of_nn2++;
            local2++;

            ++NN2it;
        }

        // Sort ids1
        ids1.sort();
        ids2.sort();

        match &= ids1.size() == ids2.size();

        for (size_t i = 0 ; i < ids1.size() ; i++)
        {match &= ids1.get(i) == ids2.get(i);}

        if (match == false)
        {break;}

        ++g_it2;
    }

    BOOST_REQUIRE_EQUAL(match,true);
    BOOST_REQUIRE(number_of_nn2 < number_of_nn);
}

BOOST_AUTO_TEST_SUITE( CellList_test )

BOOST_AUTO_TEST_CASE ( NN_radius_check )
{
    SpaceBox<2,float> box1({0.1,0.1},{0.3,0.5});
    SpaceBox<3,float> box2({0.0,0.1,0.2},{0.3,0.7,0.5});
    SpaceBox<3,float> box3({0.0,0.0,0.0},{1.0,1.0,1.0});

    std::cout << "Test cell list radius" << "\n";

    Test_NN_iterator_radius<2,float,CellList<2,float,Mem_fast<>,shift<2,float>>>(box1);
    Test_NN_iterator_radius<3,float,CellList<3,float,Mem_fast<>,shift<3,float>>>(box2);
    Test_NN_iterator_radius<3,float,CellList<3,float,Mem_fast<>,shift<3,float>>>(box3);

    std::cout << "End cell list" << "\n";
}

BOOST_AUTO_TEST_CASE( CellList_use)
{
    std::cout << "Test cell 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});
    Test_cell_s<3,double,CellList<3,double,Mem_fast<>>>(box);
    Test_cell_s<3,double,CellList<3,double,Mem_fast<>,shift<3,double>> >(box2);
    Test_cell_sM<3,double,CellListM<3,double,8>>(box);
    Test_cell_sM<3,double,CellListM<3,double,8>>(box2);


    Test_cell_s<3,double,CellList<3,double,Mem_bal<>>>(box);
    Test_cell_s<3,double,CellList<3,double,Mem_mw<>>>(box);

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

    // Test the cell list
}

BOOST_AUTO_TEST_CASE( CellList_consistent )
{
    Test_CellDecomposer_consistent<CellList<2,float,Mem_fast<>,shift<2,float>>>();
}

BOOST_AUTO_TEST_CASE( CellList_NNc_csr_calc )
{
    openfpm::vector<std::pair<grid_key_dx<3>,grid_key_dx<3>>> cNN;

    NNcalc_csr(cNN);

    BOOST_REQUIRE_EQUAL(cNN.size(),14ul);

    BOOST_REQUIRE(cNN.get(0).first == grid_key_dx<3>(0,0,0));
    BOOST_REQUIRE(cNN.get(0).second == grid_key_dx<3>(0,0,0));

    BOOST_REQUIRE(cNN.get(1).first == grid_key_dx<3>(0,0,0));
    BOOST_REQUIRE(cNN.get(1).second == grid_key_dx<3>(0,0,1));

    BOOST_REQUIRE(cNN.get(2).first == grid_key_dx<3>(0,0,1));
    BOOST_REQUIRE(cNN.get(2).second == grid_key_dx<3>(0,1,0));

    BOOST_REQUIRE(cNN.get(3).first == grid_key_dx<3>(0,0,0));
    BOOST_REQUIRE(cNN.get(3).second == grid_key_dx<3>(0,1,0));

    BOOST_REQUIRE(cNN.get(4).first == grid_key_dx<3>(0,0,0));
    BOOST_REQUIRE(cNN.get(4).second == grid_key_dx<3>(0,1,1));

    BOOST_REQUIRE(cNN.get(5).first == grid_key_dx<3>(0,1,1));
    BOOST_REQUIRE(cNN.get(5).second == grid_key_dx<3>(1,0,0));

    BOOST_REQUIRE(cNN.get(6).first == grid_key_dx<3>(0,1,0));
    BOOST_REQUIRE(cNN.get(6).second == grid_key_dx<3>(1,0,0));

    BOOST_REQUIRE(cNN.get(7).first == grid_key_dx<3>(0,1,0));
    BOOST_REQUIRE(cNN.get(7).second == grid_key_dx<3>(1,0,1));

    BOOST_REQUIRE(cNN.get(8).first == grid_key_dx<3>(0,0,1));
    BOOST_REQUIRE(cNN.get(8).second == grid_key_dx<3>(1,0,0));

    BOOST_REQUIRE(cNN.get(9).first == grid_key_dx<3>(0,0,0));
    BOOST_REQUIRE(cNN.get(9).second == grid_key_dx<3>(1,0,0));

    BOOST_REQUIRE(cNN.get(10).first == grid_key_dx<3>(0,0,0));
    BOOST_REQUIRE(cNN.get(10).second == grid_key_dx<3>(1,0,1));

    BOOST_REQUIRE(cNN.get(11).first == grid_key_dx<3>(0,0,1));
    BOOST_REQUIRE(cNN.get(11).second == grid_key_dx<3>(1,1,0));

    BOOST_REQUIRE(cNN.get(12).first == grid_key_dx<3>(0,0,0));
    BOOST_REQUIRE(cNN.get(12).second == grid_key_dx<3>(1,1,0));

    BOOST_REQUIRE(cNN.get(13).first == grid_key_dx<3>(0,0,0));
    BOOST_REQUIRE(cNN.get(13).second == grid_key_dx<3>(1,1,1));

}



BOOST_AUTO_TEST_SUITE_END()

#endif /* CELLLIST_TEST_HPP_ */