#include <VoronoiMeshSnapshot.hpp>
Public Member Functions | |
VoronoiMeshSnapshot () | |
VoronoiMeshSnapshot (const SimulationItem *item, const Box &extent, const vector< Vec > &sites, bool relax) | |
VoronoiMeshSnapshot (const SimulationItem *item, const Box &extent, SiteListInterface *sli, bool relax) | |
VoronoiMeshSnapshot (const SimulationItem *item, const Box &extent, string filename, bool relax) | |
~VoronoiMeshSnapshot () | |
int | cellIndex (Position bfr) const |
Position | centroidPosition (int m) const |
std::unique_ptr< PathSegmentGenerator > | createPathSegmentGenerator () const |
double | density (int m) const override |
double | density (Position bfr) const override |
Box | extent () const override |
Box | extent (int m) const |
void | foregoVoronoiMesh () |
Position | generatePosition () const override |
Position | generatePosition (int m) const override |
void | getEntities (EntityCollection &entities, Position bfr) const override |
void | getEntities (EntityCollection &entities, Position bfr, Direction bfk) const override |
double | mass () const override |
int | numEntities () const override |
Position | position (int m) const override |
void | readAndClose () override |
void | setExtent (const Box &extent) |
double | volume (int m) const override |
void | writeGridPlotFiles (const SimulationItem *probe) const |
![]() | |
Snapshot () | |
virtual | ~Snapshot () |
double | age (int m) const |
double | bias (int m) const |
void | close () |
double | currentMass (int m) const |
virtual double | density (int m) const =0 |
virtual double | density (Position bfr) const =0 |
virtual Box | extent () const =0 |
virtual Position | generatePosition () const =0 |
virtual Position | generatePosition (int m) const =0 |
virtual void | getEntities (EntityCollection &entities, Position bfr) const =0 |
virtual void | getEntities (EntityCollection &entities, Position bfr, Direction bfk) const =0 |
bool | hasAge () const |
bool | hasBias () const |
bool | hasCurrentMass () const |
bool | hasInitialMass () const |
bool | hasMagneticField () const |
bool | hasMetallicity () const |
bool | hasParameters () const |
bool | hasTemperature () const |
bool | hasVelocity () const |
bool | hasVelocityDispersion () const |
bool | holdsNumber () const |
void | importBias () |
void | importBox () |
void | importCurrentMass () |
void | importMagneticField () |
void | importMass () |
void | importMassDensity () |
void | importMetallicity () |
void | importNumber () |
void | importNumberDensity () |
void | importParameters (const vector< SnapshotParameter > ¶meters) |
void | importPosition () |
void | importSize () |
void | importTemperature () |
void | importVelocity () |
void | importVelocityDispersion () |
double | initialMass (int m) const |
Vec | magneticField (int m) const |
Vec | magneticField (Position bfr) const |
virtual double | mass () const =0 |
double | metallicity (int m) const |
double | metallicity (Position bfr) const |
virtual int | numEntities () const =0 |
void | open (const SimulationItem *item, string filename, string description) |
void | parameters (int m, Array ¶ms) const |
void | parameters (Position bfr, Array ¶ms) const |
virtual Position | position (int m) const =0 |
virtual void | readAndClose () |
void | setCoordinateSystem (CoordinateSystem coordinateSystem) |
void | setMassDensityPolicy (double multiplier, double maxTemperature, bool useMetallicity) |
void | setNeedGetEntities () |
double | SigmaX () const |
double | SigmaY () const |
double | SigmaZ () const |
double | temperature (int m) const |
double | temperature (Position bfr) const |
void | useColumns (string columns) |
Vec | velocity (int m) const |
Vec | velocity (Position bfr) const |
double | velocityDispersion (int m) const |
double | volume () const |
virtual double | volume (int m) const =0 |
Protected Member Functions | |
const Array & | properties (int m) const override |
![]() | |
int | ageIndex () const |
int | biasIndex () const |
int | boxIndex () const |
void | calculateDensityAndMass (Array &rhov, Array &cumrhov, double &mass) |
int | currentMassIndex () const |
int | densityIndex () const |
bool | hasMassDensityPolicy () const |
TextInFile * | infile () |
int | initialMassIndex () const |
Log * | log () const |
void | logMassStatistics (int numIgnored, double totalOriginalMass, double totalMetallicMass, double totalEffectiveMass) |
int | magneticFieldIndex () const |
int | massIndex () const |
double | maxTemperature () const |
int | metallicityIndex () const |
double | multiplier () const |
bool | needGetEntities () const |
int | numParameters () const |
int | parametersIndex () const |
int | positionIndex () const |
virtual const Array & | properties (int m) const =0 |
Random * | random () const |
void | setContext (const SimulationItem *item) |
int | sizeIndex () const |
int | temperatureIndex () const |
Units * | units () const |
bool | useMetallicity () const |
bool | useTemperatureCutoff () const |
int | velocityDispersionIndex () const |
int | velocityIndex () const |
Private Member Functions | |
void | buildMesh (bool relax) |
void | buildSearchPerBlock () |
void | buildSearchSingle () |
Node * | buildTree (vector< int >::iterator first, vector< int >::iterator last, int depth) const |
void | calculateVolume () |
bool | isPointClosestTo (Vec r, int m, const vector< int > &ids) const |
Private Attributes | |
vector< vector< int > > | _blocklists |
vector< Node * > | _blocktrees |
vector< Cell * > | _cells |
Array | _cumrhov |
double | _eps |
Box | _extent |
bool | _foregoVoronoiMesh |
double | _mass |
int | _nb |
int | _nb2 |
int | _nb3 |
Array | _rhov |
Friends | |
class | MySegmentGenerator |
Additional Inherited Members | |
![]() | |
enum class | CoordinateSystem { CARTESIAN , CYLINDRICAL , SPHERICAL } |
A VoronoiMeshSnapshot object represents a Voronoi tessellation or mesh of a cuboidal spatial domain (given a list of generating sites) and offers several facilities related to this mesh. A Voronoi mesh partitions the domain in convex polyhedra. Consider a given set of points in the domain, called sites. For each site there will be a corresponding region consisting of all points in the domain closer to that site than to any other. These regions are called Voronoi cells, and together they form the Voronoi tessellation of the domain.
The VoronoiMeshSnapshot class serves two key purposes: (1) represent snapshot data produced by a hydrodynamical simulation and imported from a column text file, defining a primary source or a transfer medium distribution, and (2) implement an unstructured spatial grid that discretizes the spatial domain as a basis for the radiative transfer simulation itself.
To support the first use case (imported snapshot data), the VoronoiMeshSnapshot class is based on the Snapshot class; it uses the facilities offered there to configure and help read the snapshot data, and it implements all functions in the general Snapshot public interface. In addition it offers functionality that is specific to this snapshot type, such as, for example, the requirement to configure the spatial extent of the domain. In this use case, the client employs the default constructor and then proceeds to configure the snapshot object as described in the Snapshot class header.
To support the second use case (spatial grid for radiative transfer), the class offers specialty constructors that accept a list of the generating sites from various sources, including a programmatically generated list or a user-supplied text column file. Note that this input file includes just the site coordinates, while a snapshot data file would include additional properties for each site. The specialty constructors automatically complete the configuration sequence of the object, so that the getters can be used immediately after construction.
Once an VoronoiMeshSnapshot object has been constructed and fully configured, its data members are no longer modified. Consequently all getters are re-entrant.
To build the Voronoi tessellation, the buildMesh() function in this class uses the code in the voro
subfolder of the SKIRT code hierarchy, which is taken from the Voro++ library written by Chris H. Rycroft (Harvard University/Lawrence Berkeley Laboratory) at https://github.com/chr1shr/voro (git commit 122531f) with minimal changes to avoid compiler warnings.
A distinguishing feature of the Voro++ library is that it carries out cell-based calculations, computing the Voronoi cell for each site individually, rather than computing the Voronoi tessellation as a global network of vertices and edges. It is therefore particularly well-suited for applications that require cell-based properties such as the cell volume, the centroid position or the number of faces or vertices. Equally important in the context of SKIRT, it is easy to distribute the work over parallel execution threads because, after setting up a common search data structure holding all sites, the calculations for the various cells are mutually independent.
The Voro++ approach also has an important drawback. Because cells are handled independently of each other, floating pointing rounding errors sometimes cause inconsistencies where the results for one cell do not properly match the corresponding results for a neighboring cell. This most often happens when the generating sites are very close to each other and/or form certain hard-to-calculate geometries. These problems manifest themselves either as empty cells or as asymmetries in the neighbor lists. To handle these situations, the buildMesh() function removes some sites from the input list as described below.
Firstly, before actually building the Voronoi tessellation, and in addition to removing sites that lie outside of the spatial domain, the function discards sites that lie closer to any previously listed site than
Secondly, after all Voronoi cells have been calculated, the function verifies the cell properties. If any of the cells have a zero volume or if any of the cell neighbors are not mutual, the involved sites are discarded and the Voronoi construction starts anew from scratch with the reduced list of sites. After a maximum of 5 attempts, the function throws a fatal error. In practice, this will hopefully never happen. Discarding incorrectly calculated cells perhaps incurs a slightly higher risk of changing the physcis of the input model. However, in practice it seems that these issues mostly occur in regions of high site density, so that the errors should be fairly limited.
VoronoiMeshSnapshot::VoronoiMeshSnapshot | ( | ) |
The default constructor initializes the snapshot in an invalid state; see the description of the required calling sequence in the Snapshot class header.
VoronoiMeshSnapshot::~VoronoiMeshSnapshot | ( | ) |
The destructor releases any data structures allocated by this class.
VoronoiMeshSnapshot::VoronoiMeshSnapshot | ( | const SimulationItem * | item, |
const Box & | extent, | ||
string | filename, | ||
bool | relax | ||
) |
This constructor reads the site positions from the specified text column file. The input file must contain three columns specifying the x,y,z coordinates. The default unit is parsec, which can be overridden by providing column header info in the file. The constructor completes the configuration for the object (but without importing mass density information or setting a mass density policy) and calls the private buildMesh() and buildSearch() functions to create the relevant data structures.
The item argument specifies a simulation item in the hierarchy of the caller (usually the caller itself) used to retrieve context such as an appropriate logger. The extent argument specifies the extent of the domain as a box lined up with the coordinate axes. Sites located outside of the domain and sites that are too close to another site are discarded. The filename argument specifies the name of the input file, including filename extension but excluding path and simulation prefix. If the relax argument is true, the function performs a single relaxation step on the site positions.
VoronoiMeshSnapshot::VoronoiMeshSnapshot | ( | const SimulationItem * | item, |
const Box & | extent, | ||
SiteListInterface * | sli, | ||
bool | relax | ||
) |
This constructor obtains the site positions from a SiteListInterface instance. The constructor completes the configuration for the object (but without importing mass density information or setting a mass density policy) and calls the private buildMesh() and buildSearch() functions to create the relevant data structures.
The item argument specifies a simulation item in the hierarchy of the caller (usually the caller itself) used to retrieve context such as an appropriate logger. The extent argument specifies the extent of the domain as a box lined up with the coordinate axes. Sites located outside of the domain and sites that are too close to another site are discarded. The sli argument specifies an object that provides the SiteListInterface interface from which to obtain the site positions. If the relax argument is true, the function performs a single relaxation step on the site positions.
VoronoiMeshSnapshot::VoronoiMeshSnapshot | ( | const SimulationItem * | item, |
const Box & | extent, | ||
const vector< Vec > & | sites, | ||
bool | relax | ||
) |
This constructor obtains the site positions from a programmatically prepared list. The constructor completes the configuration for the object (but without importing mass density information or setting a mass density policy) and calls the private buildMesh() and buildSearch() functions to create the relevant data structures.
The item argument specifies a simulation item in the hierarchy of the caller (usually the caller itself) used to retrieve context such as an appropriate logger. The extent argument specifies the extent of the domain as a box lined up with the coordinate axes. Sites located outside of the domain and sites that are too close to another site are discarded. The sites argument specifies the list of site positions. If the relax argument is true, the function performs a single relaxation step on the site positions.
|
private |
Given a list of generating sites (represented as partially initialized Cell objects), this private function builds the Voronoi tessellation and stores the corresponding cell information, including any properties relevant for supporting the interrogation capabilities offered by this class. All other data (such as Voronoi cell vertices, edges and faces) are discarded. In practice, the function adds the sites to a Voro++ container, computes the Voronoi cells one by one, and copies the relevant cell information (such as the list of neighboring cells) from the Voro++ data structures into its own.
If the relax argument is true, the function performs a single relaxation step on the site positions using Lloyd's algorithm (Lloyd 1982; Du, Faber and Gunzburger 1999, SIAM review 41.4, pp 637-676; Dobbels 2017, master thesis). An intermediate Voronoi tessellation is built using the original site positions, and subsequently each site position is replaced by the centroid (mass center) of the corresponding cell. The final tessellation is then constructed with these adjusted site positions, which are distributed more uniformly, thereby avoiding overly elongated cells in the Voronoi tessellation. Relaxation can be quite time-consuming because the Voronoi tessellation must be constructed twice.
|
private |
This private function builds data structures that allow accelerating the operation of the cellIndex() function. It assumes that the Voronoi mesh has already been built.
The domain is partitioned using a linear cubodial grid into cells that are called blocks. For each block, the function builds and stores a list of all Voronoi cells that possibly overlap that block. Retrieving the list of cells possibly overlapping a given point in the domain then comes down to locating the block containing the point (which is trivial since the grid is linear). The current implementation uses a Voronoi cell's enclosing cuboid to test for intersection with a block. Performing a precise intersection test is really slow and the shortened block lists don't substantially accelerate the cellIndex() function.
To further reduce the search time within blocks that overlap with a large number of cells, the function builds a binary search tree on the cell sites for those blocks (see for example en.wikipedia.org/wiki/Kd-tree).
|
private |
This private function builds a data structure that allows accelerating the operation of the cellIndex() function without using the Voronoi mesh. The domain is not partitioned in blocks. The function builds a single binary search tree on all cell sites (see for example en.wikipedia.org/wiki/Kd-tree).
|
private |
Private function to recursively build a binary search tree (see en.wikipedia.org/wiki/Kd-tree)
|
private |
This private function calculates the volumes for all cells without using the Voronoi mesh. It assumes that both mass and mass density columns are being imported.
int VoronoiMeshSnapshot::cellIndex | ( | Position | bfr | ) | const |
This function returns the cell index
The function uses the search data structures created by the private BuildSearch() function to accelerate its operation. It computes the appropriate block index from the coordinates of the specified point, which provides a list of Voronoi cells possibly overlapping the point. If there is a search tree for this block, the function uses it to locate the nearest point in
If the search data structures were not created during construction (which happens when using the default constructor without configuring a mass density policy), invoking the cellIndex() function causes undefined behavior.
Position VoronoiMeshSnapshot::centroidPosition | ( | int | m | ) | const |
This function returns the centroid of the Voronoi cell with index m. If the index is out of range, the behavior is undefined.
std::unique_ptr< PathSegmentGenerator > VoronoiMeshSnapshot::createPathSegmentGenerator | ( | ) | const |
This function creates and hands over ownership of a path segment generator appropriate for the adaptive mesh spatial grid, implemented as a private PathSegmentGenerator subclass. The algorithm used to construct the path is described below.
In the first stage, the function checks whether the start point is inside the domain. If so, the current point is simply initialized to the start point. If not, the function computes the path segment to the first intersection with one of the domain walls and moves the current point inside the domain. Finally the function determines the current cell, i.e. the cell containing the current point.
In the second stage, the function loops over the algorithm that computes the exit point from the current cell, i.e. the intersection of the ray formed by the current point and the path direction with the current cell's boundary. By the nature of Voronoi cells, this algorithm also produces the ID of the neigboring cell without extra cost. If an exit point is found, the loop adds a segment to the output path, updates the current point and the current cell, and continues to the next iteration. If the exit is towards a domain wall, the path is complete and the loop is terminated. If no exit point is found, which shouldn't happen too often, this must be due to computational inaccuracies. In that case, no path segment is added, the current point is advanced by a small amount, and the new current cell is determined by calling the function cellIndex().
The algorithm that computes the exit point has the following input data:
Position of the current point | |
Direction of the ray, normalized | |
ID of the cell containing the current point | |
IDs of the neighboring cells ( | |
Site position for cell | |
Domain boundaries (implicit) |
where the domain wall IDs are defined as follows:
Domain wall ID | Domain wall equation |
-1 | |
-2 | |
-3 | |
-4 | |
-5 | |
-6 |
The line containing the ray can be written as
The algorithm that computes the exit point proceeds as follows:
To calculate
and a point on the plane is given by
The equation of the plane can then be written as
Substituting
If
To calculate
|
overridevirtual |
This function returns the mass density associated with the cell with index m. If no density policy has been set or no mass information is being imported, or if the index is out of range, the behavior is undefined.
Implements Snapshot.
|
overridevirtual |
This function returns the mass density represented by the snapshot at a given point
Implements Snapshot.
|
overridevirtual |
This function returns the extent of the spatial domain as configured through the setExtent() function.
Implements Snapshot.
Box VoronoiMeshSnapshot::extent | ( | int | m | ) | const |
This function returns the bounding box (enclosing cuboid lined up with the coordinate axes) of the Voronoi cell with index m. If the index is out of range, the behavior is undefined.
void VoronoiMeshSnapshot::foregoVoronoiMesh | ( | ) |
This function configures the snapshot to skip construction of the actual Voronoi tessellation and instead use a search tree across all sites. It should be called only if (1) the snapshot has been configured to import both a mass/number density column and a volume-integrated mass/number column, and (2) the snapshot will not be required to generate random positions or trace paths. Violating these conditions will result in undefined behavior.
|
overridevirtual |
This function returns a random position within the spatial domain of the snapshot, drawn from the mass density distribution represented by the snapshot. The function first selects a random cell from the discrete probability distribution formed by the respective cell masses, and then generates a random position uniformly from the volume of that cell. If no density policy has been set or no mass information is being imported, the behavior is undefined.
Implements Snapshot.
|
overridevirtual |
This function returns a random position drawn uniformly from the (polyhedron) volume of the cell with index m. If the index is out of range, the behavior is undefined.
The function generates uniformly distributed random points in the enclosing cuboid until one happens to be inside the cell. The candidate point is inside the cell if it is closer to the cell's site position than to any neighbor cell's site positions.
Implements Snapshot.
|
overridevirtual |
This function sets the specified entity collection to the cell containing the specified point
Implements Snapshot.
|
overridevirtual |
This function replaces the contents of the specified entity collection by the set of cells crossed by the specified path with starting point
Implements Snapshot.
|
private |
This private function returns true if the given point is closer to the site with index m than to the sites with indices ids.
|
overridevirtual |
This function returns the total mass represented by the snapshot, in other words the sum of the masses of all cells. If no density policy has been set or no mass information is being imported, the behavior is undefined.
Implements Snapshot.
|
overridevirtual |
This function returns the number of sites (or, equivalently, cells) in the snapshot.
Implements Snapshot.
|
overridevirtual |
This function returns the position of the site with index m. If the index is out of range, the behavior is undefined.
Implements Snapshot.
|
overrideprotectedvirtual |
This function returns a reference to an array containing the imported properties (in column order) for the cell with index
Implements Snapshot.
|
overridevirtual |
This function reads the snapshot data from the input file, honoring the options set through the configuration functions, stores the data for later use, and closes the file by calling the base class Snapshot::readAndClose() function. Sites located outside of the domain and sites that are too close to another site are discarded. Sites with an associated temperature above the cutoff temperature (if one has been configured) are assigned a density value of zero, so that the corresponding cell has zero mass (regardless of the imported mass/density properties).
The function calls the private buildMesh() function to build the Voronoi mesh based on the imported site positions. If the snapshot configuration requires the ability to determine the density at a given spatial position, the function also calls the private buildSearch() function to create a data structure that accelerates locating the cell containing a given point.
During its operation, the function logs some statistical information about the imported snapshot and the resulting data structures.
Reimplemented from Snapshot.
void VoronoiMeshSnapshot::setExtent | ( | const Box & | extent | ) |
This function sets the extent of the spatial domain for the Voronoi mesh. When using the default constructor, this function must be called during configuration. There is no default; failing to set the extent of the domain results in undefined behavior.
|
overridevirtual |
This function returns the volume of the Voronoi cell with index m. If the index is out of range, the behavior is undefined.
Implements Snapshot.
void VoronoiMeshSnapshot::writeGridPlotFiles | ( | const SimulationItem * | probe | ) | const |
This function outputs grid plot files as described for the SpatialGridPlotProbe. The function reconstructs the Voronoi tesselation in order to produce the coordinates of the Voronoi cell vertices.