Commit 8deb38db authored by incardon's avatar incardon

Fixing multiphase Cell-list and verlet

parent cdd40c06
......@@ -3,7 +3,7 @@
AC_PREREQ(2.59)
AC_INIT(FULL-PACKAGE-NAME, VERSION, BUG-REPORT-ADDRESS)
AC_INIT(OpenFPM_data, 0.8.0, BUG-REPORT-ADDRESS)
AC_CANONICAL_SYSTEM
AC_CONFIG_SRCDIR([src/main.cpp])
......
......@@ -47,6 +47,41 @@ class grid_key_dx_iterator_sub_bc : public grid_key_dx_iterator_sub<dim,warn>
public:
grid_key_dx_iterator_sub_bc(const grid_key_dx_iterator_sub_bc & tmp)
{
this->operator=(tmp);
}
grid_key_dx_iterator_sub_bc(grid_key_dx_iterator_sub_bc && tmp)
{
this->operator=(tmp);
}
grid_key_dx_iterator_sub_bc & operator=(const grid_key_dx_iterator_sub_bc & tmp)
{
for (size_t i = 0 ; i < dim ; i++)
bc[i] = tmp.bc[i];
act = tmp.act;
boxes = tmp.boxes;
return *this;
}
grid_key_dx_iterator_sub_bc & operator=(grid_key_dx_iterator_sub_bc && tmp)
{
for (size_t i = 0 ; i < dim ; i++)
bc[i] = tmp.bc[i];
act = tmp.act;
boxes.swap(tmp.boxes);
return *this;
}
/*! \brief Constructor
*
*
......
......@@ -38,6 +38,8 @@
\endverbatim
*
* \warning A MUST BE fully contained into B
*
* The skin is the part in between the two boxes A and B
*
......@@ -51,9 +53,12 @@
template<unsigned int dim>
class grid_skin_iterator_bc
{
protected:
//! Internal iterator for each faces
grid_key_dx_iterator_sub_bc<dim> sub_it[2*dim];
private:
//! Actual iterator
size_t act;
......@@ -67,7 +72,7 @@ public:
* \param bc boundary conditions
*
*/
template <typename T> grid_skin_iterator_bc(grid_sm<dim,T> & g_sm, const Box<dim,size_t> & A, const Box<dim,size_t> & B, const size_t (& bc)[dim])
template <typename T> grid_skin_iterator_bc(const grid_sm<dim,T> & g_sm, const Box<dim,size_t> & A, const Box<dim,size_t> & B, const size_t (& bc)[dim])
:act(0)
{
for (size_t i = 0 ; i < dim ; i++)
......
......@@ -18,7 +18,7 @@ struct is_grid: std::false_type {};
*
* ### Example
*
* \snippet util.hpp Check is_grid
* \snippet util_test.hpp Check is_grid
*
* return true if T is a grid
*
......
......@@ -677,6 +677,7 @@ public:
ACTION_ON_ERROR(CELL_DECOMPOSER);
}
#endif
/* coverity[dead_error_line] */
key.set_d(s,ConvertToID(pos,s));
}
......@@ -1298,6 +1299,7 @@ Box "b" <-----------------+ | | | | | | Grid (7, 6)
/////////// Low boundary
// we are at the positive border (We are assuming that there are not rounding error in the decomposition)
/* coverity[copy_paste_error] */
if (b_d.getHigh(i) == box.getLow(i))
g_box.setHigh(i,off[i]);
......@@ -1418,6 +1420,7 @@ Box "b" <-----------------+ | | | | | | Grid (7, 6)
/////////// Low boundary
// we are at the positive border (We are assuming that there are not rounding error in the decomposition)
/* coverity[copy_paste_error] */
if (b_d.getHigh(i) == box.getLow(i))
g_box.setHigh(i,off[i]);
......
......@@ -339,6 +339,8 @@ private:
public:
typedef CellNNIteratorSym<dim,CellList<dim,T,FAST,transform,base>,RUNTIME,NO_CHECK> SymNNIterator;
//! Object type that the structure store
typedef typename base::value_type value_type;
......@@ -434,6 +436,7 @@ public:
// create the array that store the number of particle on each cell and se it to 0
InitializeStructures(this->gr_cell.getSize(),this->gr_cell.size());
from_cd = false;
}
//! Default Constructor
......
......@@ -154,8 +154,6 @@ BOOST_AUTO_TEST_CASE( ParticleIt_Cells_iterator )
while (it_cl.isNext())
{
auto p_key = it_cl.get();
count++;
++it_cl;
}
......@@ -211,8 +209,6 @@ BOOST_AUTO_TEST_CASE( ParticleIt_Cells_iterator )
while (it_cl2.isNext())
{
auto p_key = it_cl2.get();
count++;
++it_cl2;
}
......@@ -247,8 +243,6 @@ BOOST_AUTO_TEST_CASE( ParticleIt_Cells_NN_iterator )
NN.Initialize(box,div,1);
float pos[dim];
grid_key_dx_iterator<3> it(gs2);
float spacing[dim];
......@@ -321,7 +315,6 @@ BOOST_AUTO_TEST_CASE( ParticleIt_Cells_NN_iterator )
while (it_cl.isNext())
{
auto p_key = it_cl.get();
auto NN_it = it_cl.getNNIteratorCSR(vp);
size_t size_NN = 0;
......
......@@ -10,6 +10,7 @@
#include "CellNNIteratorM.hpp"
struct PV_cl
{
//! particle id
......@@ -59,6 +60,8 @@ class CellListM : public CellBase
public:
typedef CellNNIteratorSymM<dim,CellListM<dim,T,sh_byte,CellBase>,sh_byte,RUNTIME,NO_CHECK> SymNNIterator;
//! Default Constructor
CellListM()
{};
......@@ -260,36 +263,13 @@ public:
* \return Cell-list structure
*
*/
template<unsigned int impl> inline CellNNIteratorM<dim,CellListM<dim,T,sh_byte,CellBase>,sh_byte,SYM,impl> getNNIteratorSym(size_t cell)
{
CellNNIteratorM<dim,CellListM<dim,T,sh_byte,CellBase>,sh_byte,SYM,impl> cln(cell,CellListM<dim,T,sh_byte,CellBase>::NNc_sym,*this);
return cln;
}
/*! \brief Get the Neighborhood iterator
*
* It iterate across all the element of the selected cell and the near cells
*
* \verbatim
* *
x *
\endverbatim
*
* * x is the selected cell
* * * are the near cell
*
* \param cell cell id
*
* \return Cell-list structure
*
*/
template<unsigned int impl> inline CellNNIteratorM<dim,CellListM<dim,T,sh_byte,CellBase>,sh_byte,CRS,impl> getNNIteratorCross(size_t cell)
template<unsigned int impl> inline CellNNIteratorSymM<dim,CellListM<dim,T,sh_byte,CellBase>,sh_byte,SYM,impl>
getNNIteratorSym(size_t cell,
size_t p,
const openfpm::vector<Point<dim,typename CellBase::stype>> & pos,
const openfpm::vector<pos_v<dim,typename CellBase::stype>> & v)
{
CellNNIteratorM<dim,CellListM<dim,T,sh_byte,CellBase>,sh_byte,CRS,impl> cln(cell,CellListM<dim,T,sh_byte,CellBase>::NNc_cr,*this);
CellNNIteratorSymM<dim,CellListM<dim,T,sh_byte,CellBase>,sh_byte,SYM,impl> cln(cell,p,CellListM<dim,T,sh_byte,CellBase>::NNc_sym,*this,pos,v);
return cln;
}
......
......@@ -257,6 +257,13 @@ template<unsigned int dim, typename T, typename CellS> void Test_cell_sM(SpaceBo
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<dim,T>> phases;
phases.add(pos_v<dim,T>(phase1));
phases.add(pos_v<dim,T>(phase2));
size_t id = 0;
while (g_it.isNext())
......@@ -266,13 +273,15 @@ template<unsigned int dim, typename T, typename CellS> void Test_cell_sM(SpaceBo
Point<dim,T> key = Point<dim,T>(g_it.get().toPoint());
key = pmul(key,spacing) + offset[0] + box.getP1();
cl1.add(key,id,1);
++id;
phase1.add(key);
cl1.add(key,id,0);
key = Point<dim,T>(g_it.get().toPoint());
key = pmul(key,spacing) + offset[1] + box.getP1();
cl1.add(key,id,2);
phase2.add(key);
cl1.add(key,id,1);
++id;
++g_it;
......@@ -300,7 +309,7 @@ template<unsigned int dim, typename T, typename CellS> void Test_cell_sM(SpaceBo
size_t v2 = cl1.getV(cell,0);
BOOST_REQUIRE_EQUAL(n_ele,2ul);
BOOST_REQUIRE_EQUAL((long int)(p1 - p2),1);
BOOST_REQUIRE_EQUAL((long int)(p1 - p2),0);
BOOST_REQUIRE_EQUAL((long int)(v1 - v2),1);
++g_it;
}
......@@ -316,6 +325,8 @@ template<unsigned int dim, typename T, typename CellS> void Test_cell_sM(SpaceBo
Point<dim,T> key = Point<dim,T>(g_it_s.get().toPoint());
key = pmul(key,spacing) + offset[0] + box.getP1();
size_t i = g_info.LinId(g_it_s.get());
auto NN = cl1.template getNNIterator<NO_CHECK>(cl1.getCell(key));
size_t total1 = 0;
size_t total2 = 0;
......@@ -336,7 +347,7 @@ template<unsigned int dim, typename T, typename CellS> void Test_cell_sM(SpaceBo
BOOST_REQUIRE_EQUAL(total2,(size_t)openfpm::math::pow(3,dim));
auto NNSym = cl1.template getNNIteratorSym<NO_CHECK>(cl1.getCell(key));
auto NNSym = cl1.template getNNIteratorSym<NO_CHECK>(cl1.getCell(key),i,phase1,phases);
total1 = 0;
total2 = 0;
......@@ -344,7 +355,13 @@ template<unsigned int dim, typename T, typename CellS> void Test_cell_sM(SpaceBo
{
// total
if (NNSym.getV() == 1)
if (NNSym.getV() == 0 && NNSym.getP() == i)
{
++NNSym;
continue;
}
if (NNSym.getV() == 0)
total1++;
else
total2++;
......@@ -352,9 +369,39 @@ template<unsigned int dim, typename T, typename CellS> void Test_cell_sM(SpaceBo
++NNSym;
}
BOOST_REQUIRE_EQUAL(total1,(size_t)openfpm::math::pow(3,dim) / 2 + 1);
BOOST_REQUIRE_EQUAL(total1,(size_t)openfpm::math::pow(3,dim) / 2);
BOOST_REQUIRE_EQUAL(total2,(size_t)openfpm::math::pow(3,dim) / 2 + 1);
//////////////////////////////////////////////////////////////////////////
auto NNSym2 = cl1.template getNNIteratorSym<NO_CHECK>(cl1.getCell(key),i,phase2,phases);
total1 = 0;
total2 = 0;
while(NNSym2.isNext())
{
// total
if (NNSym2.getV() == 1 && NNSym2.getP() == i)
{
++NNSym2;
continue;
}
if (NNSym2.getV() == 0)
total1++;
else
total2++;
++NNSym2;
}
BOOST_REQUIRE_EQUAL(total1,(size_t)openfpm::math::pow(3,dim) / 2);
BOOST_REQUIRE_EQUAL(total2,(size_t)openfpm::math::pow(3,dim) / 2);
++i;
++g_it_s;
}
}
......
......@@ -78,4 +78,18 @@ template<unsigned int dim, typename T, typename CellList> void populate_cell_lis
populate_cell_list_sym(pos,cli,g_m);
}
/*! \brief Structure that contain a reference to a vector of particles
*
*
*/
template<unsigned int dim, typename T>
struct pos_v
{
openfpm::vector<Point<dim,T>> & pos;
pos_v(openfpm::vector<Point<dim,T>> & pos)
:pos(pos)
{}
};
#endif /* OPENFPM_DATA_SRC_NN_CELLLIST_CELLLIST_UTIL_HPP_ */
......@@ -309,4 +309,6 @@ public:
}
};
#include "CellNNIteratorRuntime.hpp"
#endif /* CELLNNITERATOR_FULL_HPP_ */
......@@ -13,6 +13,116 @@
#include <boost/integer/integer_mask.hpp>
/*! \brief Iterator for the neighborhood of the cell structures
*
* In general you never create it directly but you get it from the CellList structures
*
* It iterate across all the element of the selected cell and the near cells accordingly yo the
* symmetric scheme
*
* \tparam dim dimensionality of the space where the cell live
* \tparam Cell cell type on which the iterator is working
* \tparam NNc_size neighborhood size
* \tparam impl implementation specific options NO_CHECK do not do check on access, SAFE do check on access
*
*/
template<unsigned int dim, typename Cell, unsigned int sh_byte,unsigned int NNc_size, unsigned int impl>
class CellNNIteratorSymM : public CellNNIterator<dim,Cell,NNc_size,impl>
{
typedef boost::low_bits_mask_t<sizeof(size_t)*8-sh_byte> mask_low;
//! index of the particle p
size_t p;
//! Position of the particles p
const openfpm::vector<Point<dim,typename Cell::stype>> & pos;
//! Position of the particle p
const openfpm::vector<pos_v<dim,typename Cell::stype>> & ps;
/*! Select the next valid element
*
*/
inline void selectValid()
{
if (this->NNc[this->NNc_id] == 0)
{
while (this->start_id < this->stop_id)
{
size_t q = this->cl.get_lin(this->start_id);
for (long int i = dim-1 ; i >= 0 ; i--)
{
if (pos.template get<0>(p)[i] < ps.get(q >> (sizeof(size_t)*8-sh_byte)).pos.template get<0>(q & mask_low::sig_bits_fast)[i])
return;
else if (pos.template get<0>(p)[i] > ps.get(q >> (sizeof(size_t)*8-sh_byte)).pos.template get<0>(q & mask_low::sig_bits_fast)[i])
goto next;
}
if (q >= p) return;
next:
this->start_id++;
}
CellNNIterator<dim,Cell,NNc_size,impl>::selectValid();
}
else
{
CellNNIterator<dim,Cell,NNc_size,impl>::selectValid();
}
}
public:
/*! \brief
*
* Cell NN iterator
*
* \param cell Cell id
* \param NNc Cell neighborhood indexes (relative)
* \param cl Cell structure
*
*/
CellNNIteratorSymM(size_t cell, size_t p, const long int (&NNc)[NNc_size], Cell & cl, const openfpm::vector<Point<dim,typename Cell::stype>> & pos, const openfpm::vector<pos_v<dim,typename Cell::stype>> & ps)
:CellNNIterator<dim,Cell,NNc_size,impl>(cell,NNc,cl),p(p),pos(pos),ps(ps)
{
selectValid();
}
/*! \brief Get the value of the cell
*
* \return the next element object
*
*/
inline size_t getP()
{
return CellNNIterator<dim,Cell,NNc_size,impl>::get() & mask_low::sig_bits_fast;
}
/*! \brief Get the value of the cell
*
* \return the next element object
*
*/
inline size_t getV()
{
return (CellNNIterator<dim,Cell,NNc_size,impl>::get()) >> (sizeof(size_t)*8-sh_byte);
}
/*! \brief take the next element
*
* \return itself
*
*/
inline CellNNIteratorSymM<dim,Cell,sh_byte,NNc_size,impl> & operator++()
{
this->start_id++;
selectValid();
return *this;
}
};
/*! \brief Iterator for the neighborhood of the cell structures
*
* In general you never create it directly but you get it from the CellList structures
......@@ -111,5 +221,6 @@ public:
}
};
#include "CellNNIteratorRuntimeM.hpp"
#endif /* OPENFPM_DATA_SRC_NN_CELLLIST_CELLNNITERATORM_HPP_ */
/*
* CellNNIteratorRuntimeM.hpp
*
* Created on: Jan 26, 2017
* Author: i-bird
*/
#ifndef OPENFPM_DATA_SRC_NN_CELLLIST_CELLNNITERATORRUNTIMEM_HPP_
#define OPENFPM_DATA_SRC_NN_CELLLIST_CELLNNITERATORRUNTIMEM_HPP_
/*! \brief Iterator for the neighborhood of the cell structures
*
* In general you never create it directly but you get it from the CellList structures
*
* It iterate across all the element of the selected cell and the near cells
*
* \note to calculate quantities that involve a total reduction (like energies) use the CellIteratorSymRed
*
* \tparam dim dimensionality of the space where the cell live
* \tparam Cell cell type on which the iterator is working
* \tparam impl implementation specific options NO_CHECK do not do check on access, SAFE do check on access
*
*/
template<unsigned int dim, typename Cell, unsigned int sh_byte, unsigned int impl>
class CellNNIteratorM<dim,Cell,sh_byte,RUNTIME,impl> : public CellNNIterator<dim,Cell,RUNTIME,impl>
{
typedef boost::low_bits_mask_t<sizeof(size_t)*8-sh_byte> mask_low;
public:
/*! \brief
*
* Cell NN iterator
*
* \param cell Cell id
* \param NNc Cell neighborhood indexes (relative)
* \param cl Cell structure
*
*/
CellNNIteratorM(size_t cell, const long int * NNc, size_t NNc_size, Cell & cl)
:CellNNIterator<dim,Cell,RUNTIME,impl>(cell,NNc,NNc_size,cl)
{}
/*! \brief Get the value of the cell
*
* \return the next element object
*
*/
inline size_t getP()
{
return CellNNIterator<dim,Cell,RUNTIME,impl>::get() & mask_low::sig_bits_fast;
}
/*! \brief Get the value of the cell
*
* \return the next element object
*
*/
inline size_t getV()
{
return (CellNNIterator<dim,Cell,RUNTIME,impl>::get()) >> (sizeof(size_t)*8-sh_byte);
}
};
/*! \brief Symmetric iterator for the neighborhood of the cell structures
*
* In general you never create it directly but you get it from the CellList structures
*
* It iterate across all the element of the selected cell and the near cells.
*
* \note if we query the neighborhood of p and q is the neighborhood of p
* when we will query the neighborhood of q p is not present. This is
* useful to implement formula like \f$ \sum_{q = neighborhood(p) and p <= q} \f$
*
* \tparam dim dimensionality of the space where the cell live
* \tparam Cell cell type on which the iterator is working
* \tparam NNc_size neighborhood size
* \tparam impl implementation specific options NO_CHECK do not do check on access, SAFE do check on access
*
*/
template<unsigned int dim, typename Cell, unsigned int sh_byte, unsigned int impl>
class CellNNIteratorSymM<dim,Cell,sh_byte,RUNTIME,impl> : public CellNNIterator<dim,Cell,RUNTIME,impl>
{
typedef boost::low_bits_mask_t<sizeof(size_t)*8-sh_byte> mask_low;
//! index of the particle p
size_t p;
const openfpm::vector<Point<dim,typename Cell::stype>> & pos;
//! Position of the particles in the phases
const openfpm::vector<pos_v<dim,typename Cell::stype>> & ps;
/*! Select the next valid element
*
*/
inline void selectValid()
{
if (this->NNc[this->NNc_id] == 0)
{
while (this->start_id < this->stop_id)
{
size_t q = this->cl.get_lin(this->start_id);
for (long int i = dim-1 ; i >= 0 ; i--)
{
if (pos.template get<0>(p)[i] < ps.get(q >> (sizeof(size_t)*8-sh_byte)).pos.template get<0>(q & mask_low::sig_bits_fast)[i])
return;
else if (pos.template get<0>(p)[i] > ps.get(q >> (sizeof(size_t)*8-sh_byte)).pos.template get<0>(q & mask_low::sig_bits_fast)[i])
goto next;
}
if (q >= p) return;
next:
this->start_id++;
}
CellNNIterator<dim,Cell,RUNTIME,impl>::selectValid();
}
else
{
CellNNIterator<dim,Cell,RUNTIME,impl>::selectValid();
}
}
public:
/*! \brief
*
* Cell NN iterator
*
* \param cell Cell id
* \param NNc Cell neighborhood indexes (relative)
* \param cl Cell structure
*
*/
CellNNIteratorSymM(size_t cell,