There are two cell list data structures within Aboria. Both divide the domain into a regular grid of hypercubes with side length set so that the average number of particles within each box is close to a given value. Each particle in the container is assigned to the cell that contains its position, and neighbourhood queries search within that cell and its neighbours within the given radius.

For example, the following diagram illustrates a cell list data structure in two dimensions, shown as a regular array of grey squares each containing zero or more particles. The user wishes to find all the particles within a given euclidean distance around the red point. To accomplish this query efficiently, Aboria would then search all the red-shaded cells for particles that fall within the red circle.

The first cell list data structure supports serial insertion of particles, and parallel queries. The relevant classes are Aboria::CellList and Aboria::CellListQuery. This data structure can be selected on a per-particle-set basis, by setting the fourth template argument for Aboria::Particles. I.e.

typedef Particles<std::tuple<>, 3, std::vector, CellList>
    particle_bs_serial_t;
particle_bs_serial_t particle_bs_serial;

You will notice that we also need to specify the vector data structure that the particle container uses, which in this case is a std::vector.

The alternative is a cell-list data structure that supports parallel insertion of points, and parallel queries. This constantly re-orders the particles in the particle container so that they are sorted into individual cells, so if particles are changing cells often this can be slower. But theoretically (this hasn't been tested yet) this should speed up neighbourhood search queries as the particles that are local in memory are also local in space. The relevant classes are Aboria::CellListOrdered and Aboria::CellListOrderedQuery, and you can use this data structure like so:

typedef Particles<std::tuple<>, 3, std::vector, CellListOrdered>
    particle_bs_parallel_t;
particle_bs_parallel_t particle_bs_parallel;

The two cell-list datastructures support an alternate neighbour search facility that can be faster than the typical Aboria search iterators described above. The key assumption of this fast search is that the regular cells have a width greater than or equal to two times the search radius, and that the same search (i.e. same radius) is to be performed for every single particle in the set. If we want to use the same search radius radius as before, we can ensure this is true by setting the n_particles_in_leaf argument of the Aboria::Particles::init_neighbour_search function to $n = N\frac{(radius)^D}{V}$, where $N$ is the total number of particles in the set, $V$ is the volume of the domain, and $D$ is the number of spatial dimensions. That is,

const double required_n = N * std::pow(radius, 3) / std::pow(2.0, 3);
particles.init_neighbour_search(min, max, periodic, required_n);

Given this assumption, a fast neighbour search would be to simply look in all the possible pairs of neighbouring cells for possible neighbouring particle pairs. To enable this, Aboria provides the Aboria::get_neighbouring_buckets function, which returns an iterator that steps through all possible pairs of neighbouring buckets. The user can then iterate over each bucket pair, looping through all the particle within in bucket using either the Aboria::CellListQuery::get_bucket_particles or Aboria::CellListOrderedQuery::get_bucket_particles functions. For example, to count up the number of neighbours within a distance of radius, you might write:

for (auto ij = get_neighbouring_buckets(particles.get_query()); ij != false;
     ++ij) {
  auto tpl = *ij;
  auto i = std::get<0>(tpl); // bucket i
  auto j = std::get<1>(tpl); // bucket j
  // position offset to apply to particles in i (for periodic boundaries)
  auto poffset = std::get<2>(tpl);
  for (auto pi = particles.get_query().get_bucket_particles(i); pi != false;
       ++pi) {
    const Vector<double, 3> pi_position = get<position>(*pi) + poffset;
    for (auto pj = particles.get_query().get_bucket_particles(j);
         pj != false; ++pj) {
      if ((pi_position - get<position>(*pj)).squaredNorm() < radius) {
        // each ij bucket pair is counted once, so need to
        // increment neighbour count for pi and pj
        get<neighbours_count>(*pi)++;
        get<neighbours_count>(*pj)++;
      }
    }
  }
}

The above code considers particle pairs within neighbouring buckets, but not those within the same bucket. These pairs can be obtained by simply looping through all the buckets in the cell-list, using the Aboria::CellListQuery::get_subtree or Aboria::CellListOrderedQuery::get_subtree functions.

For example:

for (auto i = particles.get_query().get_subtree(); i != false; ++i) {
  for (auto pi = particles.get_query().get_bucket_particles(*i);
       pi != false; ++pi) {
    get<neighbours_count>(*pi)++; // self is a neighbour
    for (auto pj = pi + 1; pj != false; ++pj) {
      if ((get<position>(*pi) - get<position>(*pj)).squaredNorm() <
          radius) {
        get<neighbours_count>(*pi)++;
        get<neighbours_count>(*pj)++;
      }
    }
  }
}

After the code given above, the variable neighbour_count for each particle will contain the number of neighbouring particles around that particle. Note that this will only be correct if the width of each cell is greater than radius.

A kd-tree builds up a hierarchical tree of cells, with only the leaf cells actually containing particles. It is an efficient data structure to use if you have a high number of dimensions or if your particles are clustered in certain regions of the domain, and so you wish to adapt the size of your cells with the local particle density.

Each level of the tree divides the cells in the parent level along a certain dimension (the dimension is chosen based on the distribution of particles within the cell). Any cells that contain a number of particles that is smaller than a given threshold (set in Aboria::Particles::init_neighbour_search) are marked as leaf cells, and are not divided on subsequent levels.

There are two possible kd-trees in Aboria. The first is a custom implementation, the second wraps the popular NanoFLANN library https://github.com/jlblancoc/nanoflann. However, Aboria's native neighbourhood queries are used instead of those provided with NanoFLANN. The NanoFLANN kd-tree is fastest for running in serial (the construction and updates for this data structure are not done in parallel), where-as the custom kd-tree is best if you are running in parallel, and is the only option if you are using a GPU.

The relevant classes within Aboria are Aboria::Kdtree and Aboria::KdtreeQuery for the custom implementation, and Aboria::KdtreeNanoflann and Aboria::KdtreeNanoflannQuery for the NanoFLANN implementation. You can create a particle set using a kd-tree by setting the Aboria::Particles template arguments accordingly.

typedef Particles<std::tuple<>, 3, std::vector, Kdtree> particle_kdtree_t;
particle_kdtree_t particle_kd_tree;

or

typedef Particles<std::tuple<>, 3, std::vector, KdtreeNanoflann>
    particle_kdtree_nanoflann_t;
particle_kdtree_nanoflann_t particle_kd_tree_nanoflann;

A hyper oct-tree is a generalisation of an oct-tree (in 3 dimensions) to $N$ dimensions. Is also builds up a hierarchical tree of cells, however in this case each level of the tree is split along all dimensions, so that each cell has $2^N$ children. Any cells that contain less that the given number of particles (set in Aboria::Particles::init_neighbour_search) are marked as leaf cells. Empty cells are included in the data structure, but are ignored by any queries.

For example, the diagram below shows the leaf cells of a hyper oct-tree in 2 dimensions (this is the same as a quad-tree). If the user wishes to find all the particles within a given euclidean distance of the red particle, then Aboria will search through all the red-shaded cells for matching particles.

The relevant classes within Aboria are Aboria::octtree and Aboria::HyperOctreeQuery. You can create a particle set using a hyper oct-tree by setting the Aboria::Particles template arguments accordingly.

typedef Particles<std::tuple<>, 3, std::vector, HyperOctree>
    particle_octtree_t;
particle_octtree_t particle_octtree;

Aboria allows the user to search in a transformed coordinate space. This is useful for example in searching in a skew coordinate system, which might be neccessary if you wish to simulate a non-othorhombic periodic lattice system. All of the distance search functions (e.g. Aboria::euclidean_search), take an optional parameter transform, which must be a function object defining two overloaded operator() functions.

For convenience, Aboria provides a ready to use class for linear transformations Aboria::LinearTransform, which is created using Aboria::create_linear_transform. This class only requires the user to define the first operator() function above, and assumes that this point transform is linear with respect to the point position. For example,

struct MyTransform {
  auto operator()(const vdouble3 &v) const {
    return vdouble3(v[0] + 0.3 * v[1], v[1], v[2]);
  }
};

auto skew_x = create_linear_transform<3>(MyTransform());

for (auto i = euclidean_search(particles.get_query(), vdouble3::Constant(0),
                               radius, skew_x);
     i != false; ++i) {
  std::cout << "Found a particle with dx = " << i.dx()
            << " position = " << get<position>(*i)
            << " and id = " << get<id>(*i) << "\n";
}

The call to Aboria::euclidean_search searches for points within a radius of radius from the origin $(0,0,0)$, within new skew coordinate system defined by the transformation given by the matrix $A$

$$ A = \begin{pmatrix} 1 & 0.3 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} $$

It is important to note that the dx vector returned by i.dx() in the code above has been transformed to the new coordinate system, but all the positions stored within the particle set (e.g. get<position>(*i)) will still be in the original coordinate system.