Aboria

PrevUpHomeNext

Using the spatial data structures

The neighbour-searching functionality within Aboria uses the underlying spatial data structure (i.e. cell list, octtree, kdtree) to perform its task. You can directly access this underlying data structure as well, in order for you to write your own spatial algorithms.

Aboria uses iterators to provide a generic way to interact with the spatial data structure, so that the code that you write is independent of the particlular data structure that you use. An iterator is a generic C++ object that "iterates" through a 1D container. For example, the classic iterator is the loop index i in the code below.

size_t N = 10;
std::vector<double> v(N);
for (size_t i = 0; i < N; ++i) {
  v[i] = i;
}

The STL library abstracts away the particular type of i, and defines a set of iterators for each container. For example, the std::vector class has its own iterator, which you can use as follows.

size_t index = 0;
for (std::vector<double>::iterator i = v.begin(); i != v.end(); ++i) {
  *i = index++;
}

Or in more compact notation

index = 0;
for (auto i = v.begin(); i != v.end(); ++i) {
  *i = index++;
}

The iterators in Aboria are similar to STL iterators in that they can step through a given range of objects using the ++ operator, with some slight differences that we will describe below.

Before we can start using the data structure iterators in Aboria, we need a particle set. Lets create a random set of points in 2D:

N = 100;
typedef Particles<std::tuple<>, 2> Particles_t;
typedef typename Particles_t::position position;
Particles_t particles(N);

std::uniform_real_distribution<double> uniform(0, 1);
for (size_t i = 0; i < N; ++i) {
  auto &gen = get<generator>(particles)[i];
  get<position>(particles)[i] = vdouble2(uniform(gen), uniform(gen));
}

particles.init_neighbour_search(vdouble2(0, 0), vdouble2(1, 1),
                                vdouble2(false, false));

In order to start interacting with the spatial data structures, we need to get its query object from the particle set. A query object is a lightweight object that has all the information neccessary to access either the particle set itself or the underlying spatial data structures. It was designed because the Aboria::Particles container and the neighbour search classes (e.g. Aboria::NeighbourQueryBase) were unsuitable for copying to a gpu in order to perform calculations there, so a simpler class, the query class, was created with the required functionality.

typedef Particles_t::query_type Query_t;
Query_t query = particles.get_query();

The base class for all the query objects is Aboria::NeighbourQueryBase, and all the query classes for the individual data structures are derived from this. Now that we have query, we can create a Aboria::NeighbourQueryBase::child_iterator. This is the lowest level data structure iterator, and allows you to iterate through a set of child nodes attached to a single parent node within a tree structure.

[Note] Note

All the spatial data structures in Aboria are considered trees. For the HyperOctree and kdtree data structures, this description is obvious, but the cell list is also treated as a tree, in this case a tree with one root node having N children, where N is the total number of buckets in the cell list.

[Note] Note

Aboria tends to use the terms nodes and buckets fairly interchangably.

You can create a child_iterator by using the Aboria::NeighbourQueryBase::get_children function. This creates a child_iterator that loops through the children of the root node of the tree. We will use this iterator to loop through all these children and print out the spatial bounds of each node. In order to determine when we have reached the end of the children, we can compare the iterator to false. This pattern is widely used in Aboria, rather than specific end iterators as used in the STL.

for (auto i = query.get_children(); i != false; ++i) {
  std::cout << query.get_bounds(i) << std::endl;
}

Above we use the Aboria::NeighbourQueryBase::get_bounds function to get the bounds of the child iterator. This returns a Aboria::bbox class that contains the minimum and maximum spatial extents of the node pointed to by i.

Note that here we are using the default spatial data structure, a cell list provided by Aboria::CellList, so the "tree" here will only have 2 levels, and the loop above will loop through the second (i.e. non-root) level. We can also create a proper tree structure using the hyper oct-tree data structure given by Aboria::HyperOctree, like so:

typedef Particles<std::tuple<>, 2, std::vector, HyperOctree>
    ParticlesOcttree_t;
ParticlesOcttree_t particles_octtree(N);

for (size_t i = 0; i < N; ++i) {
  auto &gen = get<generator>(particles_octtree)[i];
  get<position>(particles_octtree)[i] =
      vdouble2(uniform(gen), uniform(gen));
}

particles_octtree.init_neighbour_search(vdouble2(0, 0), vdouble2(1, 1),
                                        vdouble2(false, false));

auto query_octtree = particles_octtree.get_query();

Now particles_octtree contains a full oct-tree, dividing the spatial domain into a hierarchical set of boxes that make up our tree data structure. The simplest iteration we might want to do on the tree is a depth-first iteration, which is easiest achieved by recursion. The Aboria::NeighbourQueryBase::get_children function can be used to get the children of a Aboria::NeighbourQueryBase::child_iterator, and using a C++ lambda function to provide the recursion we can implement a depth-first iteration like so

std::cout << "recursive depth-first" << std::endl;
for (auto i = query_octtree.get_children(); i != false; ++i) {
  std::function<void(decltype(i) &)> depth_first;
  depth_first = [&](const auto &parent) {
    std::cout << query_octtree.get_bounds(parent) << std::endl;
    for (auto i = query_octtree.get_children(parent); i != false; ++i) {
      depth_first(i);
    }
  };
  depth_first(i);
}

This construction might be a bit clumsy to use in practice however, so Aboria provides a special depth-first iterator Aboria::NeighbourQueryBase::all_iterator to allow you to write a loop equivalent to the recursive depth-first code given above.

The Aboria::NeighbourQueryBase::get_subtree function returns a Aboria::NeighbourQueryBase::all_iterator that performs a depth-first iteration over the tree. Note that you can also pass in a child_iterator to Aboria::NeighbourQueryBase::get_subtree to iterate over the sub-tree below a particular node of the tree.

std::cout << "subtree depth-first" << std::endl;
for (auto i = query_octtree.get_subtree(); i != false; ++i) {
  std::cout << query_octtree.get_bounds(i.get_child_iterator())
            << std::endl;
}

You might also want to distinguish between leaf nodes (nodes with no children) and non-leaf nodes. You can do this with the Aboria::NeighbourQueryBase::is_leaf_node function, which takes a reference to a node (rather than an iterator), and can be used like so

std::cout << "subtree depth-first showing leaf nodes" << std::endl;
for (auto i = query_octtree.get_subtree(); i != false; ++i) {
  auto ci = i.get_child_iterator();
  if (query_octtree.is_leaf_node(*ci)) {
    std::cout << "leaf node with bounds = " << query_octtree.get_bounds(ci)
              << std::endl;
  } else {
    std::cout << "non-leaf node with bounds = "
              << query_octtree.get_bounds(ci) << std::endl;
  }
}

Leaf nodes in the tree are the only nodes that contain particles. You can loop through all the particles in a given leaf node using the Aboria::NeighbourQueryBase::get_bucket_particles function, which returns an iterator. Note for non-leaf nodes, the Aboria::NeighbourQueryBase::get_bucket_particles will return an iterator that is immediatelly false, so this loop is safe even for non-leaf nodes.

std::cout << "subtree depth-first showing leaf nodes and particles"
          << std::endl;
for (auto i = query_octtree.get_subtree(); i != false; ++i) {
  auto ci = i.get_child_iterator();
  if (query_octtree.is_leaf_node(*ci)) {
    std::cout << "leaf node with bounds = " << query_octtree.get_bounds(ci)
              << std::endl;
    for (auto j = query_octtree.get_bucket_particles(*ci); j != false;
         ++j) {
      std::cout << "\t has particle with position" << get<position>(*j)
                << std::endl;
    }
  } else {
    std::cout << "non-leaf node with bounds = "
              << query_octtree.get_bounds(ci) << std::endl;
  }
}

Aboria also provides functions to query leaf nodes, or buckets, within a certain distance of a point, and these are used internally for the neighbour search functionality discussed in earlier sections. You can use the Aboria::NeighbourQueryBase::get_buckets_near_point function, which returns a Aboria::NeighbourQueryBase::query_iterator of all the buckets with a given distance of a point. This function also takes a template argument P, which refers to the p-norm distance that it uses (i.e. P=2 is the standard euclidean distance).

[Caution] Caution

The distance search provided by Aboria::NeighbourQueryBase::get_buckets_near_point does not respect the periodicity of the domain, so if you did a search near a lhs edge of a periodic domain, it would not pick up buckets on the neighbouring periodic rhs edge.

const int P = 2;
const vdouble2 search_point = vdouble2(0.5, 0.5);
const double search_radius = 0.1;
std::cout << "searching within " << search_point << " of point "
          << search_point << std::endl;
for (auto i = query_octtree.get_buckets_near_point<P>(search_point,
                                                      search_radius);
     i != false; ++i) {
  auto ci = i.get_child_iterator();
  std::cout << "\t found bucket at " << query_octtree.get_bounds(ci)
            << std::endl;
  for (auto j = query_octtree.get_bucket_particles(*ci); j != false; ++j) {
    std::cout << "\t\t found particle at " << get<position>(*j)
              << std::endl;
  }
}

PrevUpHomeNext