Back to the main page Back to the documentation 1. Overview: Tesselation 2. Creating a triangulation 3. Voronoi volume of the particles 4. Adding a vertex base

Delaunay triangulation is loosely speaking a way of linking together particles to form triangles (tetrahedra) that are 'as equilateral as possible'. The dual graph of a Delaunay triangulation is called the Voronoi diagram and provides a great way of assigning volumes to particles: a Voronoi cell contains the region in space that is closer to its center point than any other point. This has many applications like for example computing volume (instead of mass) averaged fields (DTFE / VTFE), locating halos and voids (ZOBOV) and many others. We use the CGAL library to produce a periodic tringulation in parallel. CGAL (for good reasons) don't have a MPI implementation so what we do is to produce this triangulation on each task with buffer particles and tie it together. We also perform several checks to ensure that the resulting tesselation is correct. A good advice on using this is to use it on as few tasks as you possible can.

Figure 1: A Delaunay triangulation on a set of points (black dots) and the corresponding Voronoi diagram (red cells).

The class is templated on the particle class and on the vertex base (what data do you want to store at the vertices (particles positions) ) and the create method takes in a function that assigns this data to the vertices:

template<class T, class VD = VD_FIDUCIAL>
  MPIPeriodicDelaunay<T, VD> PdT;

// This class has the create function:
void create{
  T *part, 
  size_t NumPart, 
  double buffer_fraction = 0.25, 
  double random_fraction = 0.3, 
  std::function<void(VD *, T *p)> assignment_function = fiducial_assignment_function<VD, T>);

// Create a simple triangulation with the fiducial settings
MPIPeriodicDelaunay<Particle> PdT;
PdT.create(part, NumPart);

Above buffer_fraction is the fraction of the neighboring domain that we communicate. Creating a triangulation is fastest if the points are randomly distributed in the box so shuffle the particles and also include some random particles to speed it up (yes, doing more work is actually faster!). The random_fraction is how many randoms we use compared to the real particles. Something like $10-30\%$ seems to be optimal in the tests I have done. If you don't want any data assigned to the vertices then don't provide the vertex base or the assignment function and we will use the fiducial one.

When creating the triangulation we store a handle to the vertices which you can get by calling the class methods:

std::vector<Vertex_handle> & get_vertex_handles_regular();
std::vector<Vertex_handle> & get_vertex_handles_boundary();

// Example
auto vs = PdT.get_vertex_handles_regular();

// Get the vertex of particle i
int i = 0;
auto v = vs[i];

// Get the point of particle i
auto point = v->point();

// If we created it with a vertex base, say struct VertexBase{ double quantity; } 
// fetch this for particle i
auto quantity = v->info().quantity;

Once we have created the triangulation from a set of particles we can extract the volume as follows:

std::vector<double> volumes;
PdT.VoronoiVolume(volumes);

This method fills the vector volumes with the volumes of the particles used to create the triangulation (and stored in the same order as the particles are in).

The CGAL triangulation is just a collection of points and links between them and has no info about the particles that they came from. We can deal with this by providing a vertex base and a vertex base assignement function to add data to the vertices. This can be anything you want and is really useful for computing stuff on the triangulation (though this does cost a lot of memory so be aware of this!). Here is an example of adding a ptr to the particle at the vertices:

// The vertex base (we store a pointer to the particle at the vertices)
struct VertexData{  void *part_ptr; };

// Set up the class with this vertex base
FML::TRIANGULATION::MPIPeriodicDelaunay<struct VertexData> PdT;

// Set up the vertex assignment function (NB: some vertices are randoms so they don't correspond to a real particle)
std::function<void(struct VertexData *, Particle *)> vertex_assignment_function = [](VertexData *v, Particle *p){
   if(p) v->part_ptr = p;
   else v->part_ptr = nullptr;
};

// Create the triangulation
PdT.create(part, NumPart, buffer_fraction, random_fraction, vertex_assignment_function);