FML::TRIANGULATION Namespace Reference

This namespace deals with performin triangulations (Delaunay and Voronoi) More...

Data Structures
class	MPIPeriodicDelaunay
	Create a Periodic Delaunay triangulation with particles spread across tasks. More...
struct	VertexDataWatershed
	The vertex base for the CGAL tesselation needed for the Watershed algorithm. More...
struct	WatershedBasin
	This is a general binning class for doing watershed based on [quantity] (typically volume) When we bin particles to the watershed basins we compute stuff based on what is in this class. More...

Typedefs
typedef CGAL::Exact_predicates_inexact_constructions_kernel	K
typedef CGAL::Periodic_3_Delaunay_triangulation_traits_3< K >	Gt
using	VbDS = CGAL::Periodic_3_triangulation_ds_vertex_base_3<>
using	Vb = CGAL::Triangulation_vertex_base_3< Gt, VbDS >
using	CbDS = CGAL::Periodic_3_triangulation_ds_cell_base_3<>
using	Cb = CGAL::Triangulation_cell_base_3< Gt, CbDS >
template<class T >
using	VbInfo = CGAL::Triangulation_vertex_base_with_info_3< T, Gt, Vb >
template<class T >
using	TDS = CGAL::Triangulation_data_structure_3< VbInfo< T >, Cb >
template<class T >
using	PeriodicDelaunayWithInfo = CGAL::Periodic_3_Delaunay_triangulation_3< Gt, TDS< T > >
using	VD_FIDUCIAL = CGAL::Periodic_3_triangulation_ds_vertex_base_3<>
using	Fb = CGAL::Periodic_2_triangulation_face_base_2< Gt >
using	IDType = int
using	QuantityType = float

Enumerations
enum	PointType { REGULAR_POINT , BOUNDARY_POINT , GUARD_POINT }

Functions
template<class VD , class T >
void	fiducial_assignment_function (VD *v, T *p)
	Fiducial assignment function for when we have info at vertices (the fiducial choice is no info) More...
template<class T , class U >
void	WatershedGeneral (MPIPeriodicDelaunay< T, VertexDataWatershed > &D, T *p, size_t NumPart, std::vector< double > &quantity, std::vector< U > &watershed_groups)
	In this method we perform a tesselation, then we assign [quantity] to vertices and walk the tesselation assigning the particles to Watershed basins Quantity is a vector of size NumPart. More...
template<class T , class U >
void	WatershedDensity (T *p, size_t NumPart, std::vector< U > &watershed_groups, double buffer_fraction=0.30, double random_fraction=0.5, bool do_density_maximum=false)
	In this method we perform a tesselation, compute voronoi volumes and assign density to each cell mass_of_part/volume and then we locate the density minima (or maxima) and walk the tesselation assigning the particles to Watershed basins If density_maximum = true then we assign 1/density to the particles and run the same algorithm so we locate density peaks instead of density minima. More...

Variables
const IDType	NoWatershedID = std::numeric_limits<IDType>::max()
const QuantityType	Infinity = std::numeric_limits<QuantityType>::infinity()

Detailed Description

This namespace deals with performin triangulations (Delaunay and Voronoi)

Typedef Documentation

Cb

using FML::TRIANGULATION::Cb = typedef CGAL::Triangulation_cell_base_3<Gt, CbDS>

Definition at line 56 of file PeriodicDelaunay.h.

CbDS

using FML::TRIANGULATION::CbDS = typedef CGAL::Periodic_3_triangulation_ds_cell_base_3<>

Definition at line 55 of file PeriodicDelaunay.h.

Fb

using FML::TRIANGULATION::Fb = typedef CGAL::Periodic_2_triangulation_face_base_2<Gt>

Definition at line 69 of file PeriodicDelaunay.h.

Gt

typedef CGAL::Periodic_2_Delaunay_triangulation_traits_2< K > FML::TRIANGULATION::Gt

Definition at line 52 of file PeriodicDelaunay.h.

IDType

using FML::TRIANGULATION::IDType = typedef int

Definition at line 695 of file PeriodicDelaunay.h.

K

typedef CGAL::Exact_predicates_inexact_constructions_kernel FML::TRIANGULATION::K

Definition at line 50 of file PeriodicDelaunay.h.

PeriodicDelaunayWithInfo

template<class T >

using FML::TRIANGULATION::PeriodicDelaunayWithInfo = typedef CGAL::Periodic_2_Delaunay_triangulation_2<Gt, TDS<T> >

Definition at line 62 of file PeriodicDelaunay.h.

QuantityType

using FML::TRIANGULATION::QuantityType = typedef float

Definition at line 696 of file PeriodicDelaunay.h.

TDS

template<class T >

using FML::TRIANGULATION::TDS = typedef CGAL::Triangulation_data_structure_2<VbInfo<T>, Fb>

Definition at line 60 of file PeriodicDelaunay.h.

Vb

typedef CGAL::Triangulation_vertex_base_2< Gt > FML::TRIANGULATION::Vb

Definition at line 54 of file PeriodicDelaunay.h.

VbDS

using FML::TRIANGULATION::VbDS = typedef CGAL::Periodic_3_triangulation_ds_vertex_base_3<>

Definition at line 53 of file PeriodicDelaunay.h.

VbInfo

template<class T >

using FML::TRIANGULATION::VbInfo = typedef CGAL::Triangulation_vertex_base_with_info_2<T, Gt, Vb>

Definition at line 58 of file PeriodicDelaunay.h.

VD_FIDUCIAL

typedef CGAL::Periodic_2_triangulation_vertex_base_2< Gt > FML::TRIANGULATION::VD_FIDUCIAL

Definition at line 63 of file PeriodicDelaunay.h.

Enumeration Type Documentation

PointType

enum FML::TRIANGULATION::PointType

Enumerator
REGULAR_POINT
BOUNDARY_POINT
GUARD_POINT

Definition at line 692 of file PeriodicDelaunay.h.

{ REGULAR_POINT, BOUNDARY_POINT, GUARD_POINT };

Function Documentation

fiducial_assignment_function()

template<class VD , class T >

void FML::TRIANGULATION::fiducial_assignment_function	(	VD *	v,
		T *	p
	)

Fiducial assignment function for when we have info at vertices (the fiducial choice is no info)

Template Parameters

VD	The vertex base
The	particle class

Parameters

[in]	v	Pointer to the vertex
[in]	p	Pointer to the particle

Definition at line 88 of file PeriodicDelaunay.h.

{}

WatershedDensity()

template<class T , class U >

void FML::TRIANGULATION::WatershedDensity	(	T *	p,
		size_t	NumPart,
		std::vector< U > &	watershed_groups,
		double	buffer_fraction = 0.30,
		double	random_fraction = 0.5,
		bool	do_density_maximum = false
	)

In this method we perform a tesselation, compute voronoi volumes and assign density to each cell mass_of_part/volume and then we locate the density minima (or maxima) and walk the tesselation assigning the particles to Watershed basins If density_maximum = true then we assign 1/density to the particles and run the same algorithm so we locate density peaks instead of density minima.

The method sets the volume of each particle and the particle must have a set_volume method.

Template Parameters

T	The particle class
U	The watershed binning class

Parameters

[in]	p	Pointer to the particles
[in]	NumPart	Number of local particles.
[out]	watershed_groups	The result of the watershed: a list of watershed groups.
[in]	buffer_fraction	Optional. How big part of the neighbor domain do we include as the buffer.
[in]	random_fraction	Optional. How many (as fraction of the normal particles) random particles do we add (this is to help speed up the tesslation).
[in]	do_density_maximum	Optional. Watershed based on the density (false) or 1/density (true).

Definition at line 1304 of file PeriodicDelaunay.h.

                                                               {
 
            static_assert(FML::PARTICLE::has_set_volume<T>(),
                          "[WatershedDensity] We require the particle to have a set_volume / get_volume method");
 
            if (FML::ThisTask == 0) {
                if (do_density_maximum)
                    std::cout << "[WatershedDensity] We will locate minima of 1/(voronoid density) of particles\n";
                else
                    std::cout << "[WatershedDensity] We will locate minima of the voronoi density of particles\n";
                std::cout << "[WatershedDensity] Buffer fraction " << buffer_fraction << " of neighbor domains\n";
            }
 
            // Vertex assignement function
            std::function<void(VertexDataWatershed *, T *)> vertex_assignment_function =
                [](VertexDataWatershed * v, T * p) { v->part_ptr = p ? (void *)p : nullptr; };
 
            // Create tesselation
            if (FML::ThisTask == 0)
                std::cout << "[WatershedDensity] Computing tesselation\n";
            MPIPeriodicDelaunay<T, VertexDataWatershed> D;
            D.create(p, NumPart, buffer_fraction, random_fraction, vertex_assignment_function);
            // double dx_buffer = D.get_dx_buffer();
 
            // Compute Voronoi volumes
            if (FML::ThisTask == 0)
                std::cout << "[WatershedDensity] Computing voronoi volumes from tesselation\n";
            std::vector<double> volumes;
            D.VoronoiVolume(volumes);
 
            // Set the volume to the particles
            if constexpr (FML::PARTICLE::has_set_volume<T>()) {
                for (size_t i = 0; i < NumPart; i++) {
                    FML::PARTICLE::SetVolume(p[i], volumes[i]);
                }
            }
 
            // Convert volumes to mass density
            // If you want to tesselated based on an other property then this is all we
            // need to change
            std::vector<double> & density = volumes;
            double mass = 1.0;
            for (size_t i = 0; i < NumPart; i++) {
                if constexpr (FML::PARTICLE::has_get_mass<T>()) {
                    mass = FML::PARTICLE::GetMass(p[i]);
                }
                if (do_density_maximum)
                    // To find clusters use 1/density as the quantity to watershed
                    density[i] = volumes[i] / mass;
                else
                    // To find voids use density as watershed quantity
                    density[i] = mass / volumes[i];
            }
 
            // Run general algorithm
            WatershedGeneral<T, U>(D, p, NumPart, density, watershed_groups);
        }

WatershedGeneral()

template<class T , class U >

void FML::TRIANGULATION::WatershedGeneral	(	MPIPeriodicDelaunay< T, VertexDataWatershed > &	D,
		T *	p,
		size_t	NumPart,
		std::vector< double > &	quantity,
		std::vector< U > &	watershed_groups
	)

In this method we perform a tesselation, then we assign [quantity] to vertices and walk the tesselation assigning the particles to Watershed basins Quantity is a vector of size NumPart.

When quantity is the density or inverse density of the particle then we get a Void or Cluster finder

If the buffer is too small the results will not be perfect, we give warnings when this is the case (instead of just throwing) as some times this is fine. To be sure of the result its a good idea to try with a smaller number of CPUs or large buffer and check that the results match

Template Parameters

T	The particle class
U	The watershed binning class

Parameters

[in]	D	MPIPeriodicDelaunay tesselation (already created).
[in]	p	Pointer to the particles.
[in]	NumPart	Number of local particles.
[in]	quantity	Vector with the quantity to watershed on (e.g. the density of the particles). quantity[i] corresponds to the quantity for particle i.
[out]	watershed_groups	The result of the watershed: a list of watershed groups.

Definition at line 733 of file PeriodicDelaunay.h.

                                                               {
 
            assert(quantity.size() == NumPart);
            [[maybe_unused]] const double dx_buffer = D.get_dx_buffer();
 
            // Fetch tesselation
            auto dt = D.get_delaunay_triangulation();
            using Vertex_handle = typename decltype(dt)::Vertex_handle;
 
            // For boundaries and communication
            [[maybe_unused]] size_t count_left = 0;
            [[maybe_unused]] size_t count_right = 0;
            [[maybe_unused]] size_t recv_left = 0;
            [[maybe_unused]] size_t recv_right = 0;
            [[maybe_unused]] size_t nboundary = 0;
 
            // Assign quantity to particles
            // Here we fetch the vertex handles from the tesselation
            auto vs = D.get_vertex_handles_regular();
            assert(vs.size() == NumPart);
            for (size_t i = 0; i < NumPart; i++) {
                vs[i]->info().quantity = quantity[i];
                vs[i]->info().point_type = REGULAR_POINT;
            }
 
            // Deal with the boundary
#ifdef USE_MPI
            int LeftTask = (FML::ThisTask - 1 + FML::NTasks) % FML::NTasks;
            int RightTask = (FML::ThisTask + 1) % FML::NTasks;
            if (FML::NTasks > 1) {
 
                // Count how many particles we have on the left and right
                // NB: this relies on we do exactly the same thing in MPIPeriodicDelaunay
                // A more roboust way is to communicate index we want, gather, send etc.
                count_left = 0;
                count_right = 0;
                for (size_t i = 0; i < NumPart; i++) {
                    auto * pos = FML::PARTICLE::GetPos(p[i]);
                    if (pos[0] < FML::xmin_domain + dx_buffer)
                        count_left++;
                    if (pos[0] > FML::xmax_domain - dx_buffer)
                        count_right++;
                }
 
                // Comunicate how many to send
                MPI_Status status;
                MPI_Sendrecv(&count_left,
                             sizeof(count_left),
                             MPI_CHAR,
                             LeftTask,
                             0,
                             &recv_right,
                             sizeof(recv_right),
                             MPI_CHAR,
                             RightTask,
                             0,
                             MPI_COMM_WORLD,
                             &status);
                MPI_Sendrecv(&count_right,
                             sizeof(count_right),
                             MPI_CHAR,
                             RightTask,
                             0,
                             &recv_left,
                             sizeof(recv_left),
                             MPI_CHAR,
                             LeftTask,
                             0,
                             MPI_COMM_WORLD,
                             &status);
 
                // Gather quantity
                std::vector<QuantityType> quantity_to_send_left(count_left);
                std::vector<QuantityType> quantity_to_send_right(count_right);
                std::vector<QuantityType> quantity_to_recv_left(recv_left);
                std::vector<QuantityType> quantity_to_recv_right(recv_right);
                count_left = 0;
                count_right = 0;
                for (size_t i = 0; i < NumPart; i++) {
                    auto * pos = FML::PARTICLE::GetPos(p[i]);
                    if (pos[0] < FML::xmin_domain + dx_buffer) {
                        quantity_to_send_left[count_left] = quantity[i];
                        count_left++;
                    }
                    if (pos[0] > FML::xmax_domain - dx_buffer) {
                        quantity_to_send_right[count_right] = quantity[i];
                        count_right++;
                    }
                }
 
                // Send quantity
                MPI_Sendrecv(quantity_to_send_left.data(),
                             sizeof(QuantityType) * count_left,
                             MPI_CHAR,
                             LeftTask,
                             0,
                             quantity_to_recv_right.data(),
                             sizeof(QuantityType) * recv_right,
                             MPI_CHAR,
                             RightTask,
                             0,
                             MPI_COMM_WORLD,
                             &status);
                MPI_Sendrecv(quantity_to_send_right.data(),
                             sizeof(QuantityType) * count_right,
                             MPI_CHAR,
                             RightTask,
                             0,
                             quantity_to_recv_left.data(),
                             sizeof(QuantityType) * recv_left,
                             MPI_CHAR,
                             LeftTask,
                             0,
                             MPI_COMM_WORLD,
                             &status);
 
                // Assign quantity to particles
                // Here we fetch the vertex handles from the tesselation
                auto vs_boundary = D.get_vertex_handles_boundary();
                nboundary = recv_left + recv_right;
                assert(vs_boundary.size() == nboundary);
                for (size_t i = 0; i < recv_left; i++) {
                    vs_boundary[i]->info().quantity = quantity_to_recv_left[i];
                    vs_boundary[i]->info().point_type = BOUNDARY_POINT;
                }
                for (size_t i = 0; i < recv_right; i++) {
                    vs_boundary[i + recv_left]->info().quantity = quantity_to_recv_right[i];
                    vs_boundary[i + recv_left]->info().point_type = BOUNDARY_POINT;
                }
 
                // Assign minimum neighbor
                for (size_t i = 0; i < nboundary; i++) {
                    auto v = vs_boundary[i];
                    std::vector<Vertex_handle> vertices;
                    dt.adjacent_vertices(v, std::back_inserter(vertices));
                    auto min_quantity = Infinity;
                    for (auto & vnew : vertices) {
                        if (vnew->info().point_type == GUARD_POINT)
                            continue;
                        if (vnew->info().quantity < min_quantity)
                            min_quantity = vnew->info().quantity;
                    }
                    v->info().min_quantity_nbor = min_quantity;
                    if (v->info().min_quantity_nbor == Infinity) {
                        std::cout << "[WatershedGeneral] Warning boundary particle has no normal nbors " << v->point()
                                  << " " << FML::ThisTask << "\n";
                    }
                }
            }
#endif
 
            // Find quantity minima
            std::vector<char> is_minimum(NumPart, 0);
            for (size_t i = 0; i < NumPart; i++) {
                std::vector<Vertex_handle> vertices;
                dt.adjacent_vertices(vs[i], std::back_inserter(vertices));
                auto quantity = vs[i]->info().quantity;
                is_minimum[i] = 1;
                for (auto & v : vertices) {
                    if (v == vs[i])
                        continue;
                    if (v->info().quantity < quantity) {
                        is_minimum[i] = 0;
                        break;
                    }
                }
            }
 
            // Count number of minima
            size_t nminima = 0;
            for (size_t i = 0; i < NumPart; i++) {
                if (is_minimum[i] == 1) {
                    nminima++;
                }
            }
 
            // Ensure unique task id
            std::vector<int> nminima_per_task(FML::NTasks, 0);
            nminima_per_task[FML::ThisTask] = nminima;
#ifdef USE_MPI
            MPI_Allreduce(MPI_IN_PLACE, nminima_per_task.data(), FML::NTasks, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
#endif
            int id_start = 0;
            for (int i = 0; i < FML::ThisTask; i++)
                id_start += nminima_per_task[i];
            size_t ntotal_minima = 0;
            for (int i = 0; i < FML::NTasks; i++)
                ntotal_minima += size_t(nminima_per_task[i]);
 
            if (FML::ThisTask == 0)
                std::cout << "[WatershedGeneral] We found " << ntotal_minima << " minimum points\n";
 
            // Assign unique wathershed ID - this ID is unique across tasks
            // Gather vertex handles to all the minimas
            decltype(vs) vminima(nminima);
            nminima = 0;
            for (size_t i = 0; i < NumPart; i++) {
                if (is_minimum[i]) {
                    auto v = vs[i];
                    vminima[nminima] = v;
                    v->info().WatershedID = id_start + nminima;
                    nminima++;
                }
            }
            is_minimum.clear();
            is_minimum.shrink_to_fit();
 
            // Find the smallest quantity in the nbors
            for (size_t i = 0; i < NumPart; i++) {
                auto v = vs[i];
                std::vector<Vertex_handle> vertices;
                dt.adjacent_vertices(v, std::back_inserter(vertices));
                // auto quantity = v->info().quantity;
                auto min_quantity = Infinity;
                for (auto & vnew : vertices) {
                    if (vnew->info().point_type == GUARD_POINT)
                        continue;
                    if (vnew->info().quantity < min_quantity)
                        min_quantity = vnew->info().quantity;
                }
                assert(min_quantity < Infinity);
                v->info().min_quantity_nbor = min_quantity;
            }
 
            // Recursive call to find the watershed basins we set all neighbors with higher values of quantity as
            // part of the basin and then do a recursive call with each of these neighbors as the starting point
            std::function<int(Vertex_handle, IDType)> recursive_through_vertices = [&](Vertex_handle v, IDType id) {
                if (v->info().point_type == GUARD_POINT)
                    return 0;
 
                // Assign current vertex
                int count = 0;
                v->info().WatershedID = id;
                if (v->info().point_type == REGULAR_POINT)
                    count++;
 
                // Loop over neighbors
                std::vector<Vertex_handle> vertices;
                dt.adjacent_vertices(v, std::back_inserter(vertices));
                auto quantity = v->info().quantity;
                for (auto & vnew : vertices) {
                    if (vnew == v)
                        continue;
                    if (vnew->info().WatershedID == id)
                        continue;
                    if (vnew->info().point_type == GUARD_POINT)
                        continue;
                    if (quantity > vnew->info().min_quantity_nbor)
                        continue;
                    if (vnew->info().quantity > quantity) {
                        if (vnew->info().point_type == REGULAR_POINT) {
                            count++;
                        }
                        vnew->info().WatershedID = id;
                        count += recursive_through_vertices(vnew, id);
                    }
                }
                return count;
            };
 
            // We can now do this step with OpenMP if we want as no two basins should be able to overlap
            size_t ntot = 0;
            for (size_t i = 0; i < nminima; i++) {
                auto v = vminima[i];
                auto id = v->info().WatershedID;
                ntot += recursive_through_vertices(v, id);
            }
#ifdef DEBUG_TESSELATION
            std::cout << "We have " << nminima << " minima with " << ntot << " of " << NumPart
                      << " particles tied up on task " << FML::ThisTask << "\n";
#endif
 
            // Now we need to merge across tasks
            // We communicate boundary particles and their WatershedID to the left and right and do the merging
            // If some of the particles have id > 0 then they will be merged
            // We should really do this thing one by one (i,i+1) to ensure structures that span many tasks gets
            // assigned correctly. This is rarely an issue unless we have extremely few particles and many tasks
#ifdef USE_MPI
            if (FML::NTasks > 1) {
 
                std::vector<IDType> watershed_id_to_send_left(count_left);
                std::vector<IDType> watershed_id_to_send_right(count_right);
                std::vector<IDType> watershed_id_to_recv_left(recv_left);
                std::vector<IDType> watershed_id_to_recv_right(recv_right);
 
                // We merge back and forth until we are done
                // or maximum 3 times
                for (int s = 0; s < 3; s++) {
                    auto vs = D.get_vertex_handles_regular();
                    assert(vs.size() == NumPart);
                    count_left = 0;
                    count_right = 0;
                    for (size_t i = 0; i < NumPart; i++) {
                        auto id = vs[i]->info().WatershedID;
                        auto * pos = FML::PARTICLE::GetPos(p[i]);
                        if (pos[0] < FML::xmin_domain + dx_buffer) {
                            watershed_id_to_send_left[count_left] = id;
                            count_left++;
                        }
                        if (pos[0] > FML::xmax_domain - dx_buffer) {
                            watershed_id_to_send_right[count_right] = id;
                            count_right++;
                        }
                    }
                    assert(watershed_id_to_send_left.size() == count_left);
                    assert(watershed_id_to_send_right.size() == count_right);
 
                    MPI_Status status;
                    MPI_Sendrecv(watershed_id_to_send_left.data(),
                                 sizeof(IDType) * count_left,
                                 MPI_CHAR,
                                 LeftTask,
                                 0,
                                 watershed_id_to_recv_right.data(),
                                 sizeof(IDType) * recv_right,
                                 MPI_CHAR,
                                 RightTask,
                                 0,
                                 MPI_COMM_WORLD,
                                 &status);
                    MPI_Sendrecv(watershed_id_to_send_right.data(),
                                 sizeof(IDType) * count_right,
                                 MPI_CHAR,
                                 RightTask,
                                 0,
                                 watershed_id_to_recv_left.data(),
                                 sizeof(IDType) * recv_left,
                                 MPI_CHAR,
                                 LeftTask,
                                 0,
                                 MPI_COMM_WORLD,
                                 &status);
 
                    // Update Watershed IDs
                    auto vs_boundary = D.get_vertex_handles_boundary();
                    for (size_t i = 0; i < recv_left; i++) {
                        auto v = vs_boundary[i];
                        auto id = watershed_id_to_recv_left[i];
                        v->info().WatershedID = id;
                    }
                    for (size_t i = 0; i < recv_right; i++) {
                        auto v = vs_boundary[i + recv_left];
                        auto id = watershed_id_to_recv_right[i];
                        v->info().WatershedID = id;
                    }
 
                    // Assign IDs to vertices recursively
                    int ntot = 0;
                    for (size_t i = 0; i < nboundary; i++) {
                        auto v = vs_boundary[i];
                        auto id = v->info().WatershedID;
                        if (id == NoWatershedID)
                            continue;
                        ntot += recursive_through_vertices(v, id);
                    }
#ifdef DEBUG_TESSELATION
                    std::cout << "In merging on " << FML::ThisTask << " we assigned " << ntot << "\n";
#endif
                    // If no more particles gets assigned we stop
                    FML::SumOverTasks(&ntot);
                    if (ntot == 0)
                        break;
                    if (s == 2 and FML::ThisTask == 0)
                        std::cout << "[WatershedGeneral] Warning we merged 3 times and still not done, buffer should "
                                     "be larger, cannot guarantee the results are perfect\n";
                }
            }
#endif
 
            // Count how many regular particles we have processed
            long long int assigned = 0;
            for (size_t i = 0; i < NumPart; i++) {
                if (vs[i]->info().WatershedID != NoWatershedID) {
                    assigned++;
                } else {
                    std::cout << "[WatershedGeneral] Warning: on task " << FML::ThisTask << " particle  " << i
                              << " pos: " << vs[i]->point() << " has not been assigned!\n";
                }
            }
#ifdef DEBUG_TESSELATION
            std::cout << "Processed on task " << FML::ThisTask << " is " << assigned << " / " << NumPart << "\n";
#endif
 
            long long int npartotal = NumPart;
            FML::SumOverTasks(&assigned);
            FML::SumOverTasks(&npartotal);
 
            if (FML::ThisTask == 0) {
                std::cout << "[WatershedGeneral] Total particles assigned " << assigned << " / " << npartotal << "\n";
                if (assigned < npartotal) {
                    std::cout << "[WatershedGeneral] Warning we have not mange to assign all particles to their "
                                 "respective minima, cannot guarantee the results are perfect\n";
                }
            }
 
            // Compute how many particles there are per group and communicate this
            std::vector<int> n_per_group(ntotal_minima, 0);
            for (size_t i = 0; i < NumPart; i++) {
                if (vs[i]->info().WatershedID != NoWatershedID) {
                    n_per_group[vs[i]->info().WatershedID]++;
                }
            }
#ifdef USE_MPI
            MPI_Allreduce(MPI_IN_PLACE, n_per_group.data(), ntotal_minima, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
#endif
 
            // Compute the minimum position and make sure all tasks have this
            std::vector<double> x_minima(ntotal_minima, 0.0);
            std::vector<double> y_minima(ntotal_minima, 0.0);
            std::vector<double> z_minima(ntotal_minima, 0.0);
            for (size_t i = 0; i < nminima; i++) {
                auto v = vminima[i];
 
                // Get points in dual polyhedra, i.e. voronoi cell
                // std::vector< decltype(p) > pts;
                // auto tetra = dt.dual(v, std::back_inserter(pts));
                // std::cout << "Dual polyhedra: " << pts.size() << "\n";
 
                // XXX Add compute circumcenter of core particle and 3 closest low density region
 
                // Minimum density particle
                auto pos = v->point();
                x_minima[v->info().WatershedID] = pos[0];
                y_minima[v->info().WatershedID] = pos[1];
#if CGAL_NDIM == 3
                z_minima[v->info().WatershedID] = pos[2];
#endif
            }
#ifdef USE_MPI
            MPI_Allreduce(MPI_IN_PLACE, x_minima.data(), ntotal_minima, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
            MPI_Allreduce(MPI_IN_PLACE, y_minima.data(), ntotal_minima, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
            MPI_Allreduce(MPI_IN_PLACE, z_minima.data(), ntotal_minima, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
#endif
 
            // Time to compile up the results. Make the result array
            watershed_groups.resize(ntotal_minima);
 
            // Initialize groups on all tasks
            for (size_t i = 0; i < ntotal_minima; i++) {
#if CGAL_NDIM == 2
                double pos[] = {x_minima[i], y_minima[i]};
#elif CGAL_NDIM == 3
                double pos[] = {x_minima[i], y_minima[i], z_minima[i]};
#endif
                watershed_groups[i].init(pos);
            }
 
            // Loop through all particles and assign data to groups
            for (size_t i = 0; i < NumPart; i++) {
                auto v = vs[i];
                auto id = v->info().WatershedID;
                if (id != NoWatershedID) {
                    watershed_groups[id].add_particle((T *)v->info().part_ptr, v->info().quantity);
                }
            }
 
            // Communicate stuff so that task 0 has all the data
            // NB: assumes watershed_groups is a simple type so that sizeof works, i.e. no dynamic allocated
            // objects in the class
#ifdef USE_MPI
            for (int i = 1; i < FML::NTasks; i++) {
                size_t bytes = sizeof(U) * ntotal_minima;
                if (FML::ThisTask == 0) {
                    std::vector<U> watershed_groups_from_other_task(ntotal_minima);
                    MPI_Status status;
                    MPI_Recv(watershed_groups_from_other_task.data(), bytes, MPI_BYTE, i, 0, MPI_COMM_WORLD, &status);
 
                    // Merge in the groups
                    for (size_t j = 0; j < ntotal_minima; j++) {
                        watershed_groups[j].merge(watershed_groups_from_other_task[j]);
                    }
                } else if (FML::ThisTask == i) {
                    MPI_Send(watershed_groups.data(), bytes, MPI_BYTE, 0, 0, MPI_COMM_WORLD);
                }
                MPI_Barrier(MPI_COMM_WORLD);
            }
#endif
 
            // Finalize the binning
            if (FML::ThisTask == 0) {
                for (size_t i = 0; i < ntotal_minima; i++) {
                    watershed_groups[i].finalize();
                }
            } else {
                // We only keep the data on task 0 so clear the rest to avoid any confusion
                watershed_groups.clear();
                watershed_groups.shrink_to_fit();
            }
 
            // XXX The stuff below should be in WatershedDensity or we should store quantity_min above!
            // If we sort we need to keep the id
            // Sort them by minimum density from large to small
            // std::sort(watershed_groups.begin(), watershed_groups.end(), [](const auto i, const auto j){ return
            // i.density_min > j.density_min; } );
 
            /*
            // NB: only task 0 has groups
            // Second watershed merging. Compile up list of links between groups and send it all to task 0
            std::vector< std::set<IDType> > set_of_group_links(ntotal_minima);
            for(size_t i = 0; i < NumPart; i++){
            auto v = vs[i];
            auto id = v->info().WatershedID;
            auto & curset = set_of_group_links[id];
            std::vector<Vertex_handle> vertices;
            dt.adjacent_vertices(v, std::back_inserter(vertices));
            std::vector<IDType> links;
            for(auto &vnew : vertices){
            auto nborid = vnew->info().WatershedID;
            if(nborid != id and nborid != NoWatershedID){
            curset.insert(nborid); // Add to list
            }
            }
            }
 
            // We now have a list of sets containing all the links between groups. Compile this up
            std::vector<int> num_links(ntotal_minima,0);
            std::vector<std::vector<IDType>> links(ntotal_minima);
            for(int i = 0; i < ntotal_minima; i++){
            num_links[i] = set_of_group_links[i].size();
            links[i].reserve(num_links[i]);
            for(auto val : set_of_group_links[i]){
            links[i].push_back(val);
            }
            std::cout << FML::ThisTask << " group " << i << " has " << num_links[i] << " neighbors\n";
            set_of_group_links.clear();
            set_of_group_links.shrink_to_fit();
            }
 
            // We now have all the links. We should really do a communication here
            // to ensure we have all the same links. We typically have 10-50 links.
            // Now we sort ntotal_minima group ids by density and merge them. We only
            // keep the groups that are higher than the minimum density
 
            for(int i = 0; i < ntotal_minima; i++){
            for(auto nborid : links){
            // Check if density is higher if so then merge
            }
            }
 
             */
 
            // Compile up struct with {id, quantity_min} for each group. Then sort by quantity_min
            // Now loop over this struct. index tells us the location in watershed_groups
        }

Variable Documentation

Infinity

const QuantityType FML::TRIANGULATION::Infinity = std::numeric_limits<QuantityType>::infinity()

Definition at line 698 of file PeriodicDelaunay.h.

NoWatershedID

const IDType FML::TRIANGULATION::NoWatershedID = std::numeric_limits<IDType>::max()

Definition at line 697 of file PeriodicDelaunay.h.