Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Properly sample all properties from imported particle media #228

Merged
merged 2 commits into from
Dec 10, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 0 additions & 9 deletions SKIRT/core/AdaptiveMeshSnapshot.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -11,9 +11,7 @@
#include "PathSegmentGenerator.hpp"
#include "Random.hpp"
#include "SpatialGridPath.hpp"
#include "StringUtils.hpp"
#include "TextInFile.hpp"
#include "Units.hpp"

////////////////////////////////////////////////////////////////////

Expand Down Expand Up @@ -373,13 +371,6 @@ const Array& AdaptiveMeshSnapshot::properties(int m) const

////////////////////////////////////////////////////////////////////

int AdaptiveMeshSnapshot::nearestEntity(Position bfr) const
{
return cellIndex(bfr);
}

////////////////////////////////////////////////////////////////////

void AdaptiveMeshSnapshot::getEntities(EntityCollection& entities, Position bfr) const
{
entities.addSingle(cellIndex(bfr));
Expand Down
4 changes: 0 additions & 4 deletions SKIRT/core/AdaptiveMeshSnapshot.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -228,10 +228,6 @@ class AdaptiveMeshSnapshot : public Snapshot
range, the behavior is undefined. */
const Array& properties(int m) const override;

/** This function returns the index \f$0\le m \le N_\mathrm{ent}-1\f$ of the cell containing
the specified point \f${\bf{r}}\f$, or -1 if the point is outside the domain. */
int nearestEntity(Position bfr) const override;

public:
/** This function sets the specified entity collection to the cell containing the specified
point \f${\bf{r}}\f$, or to the empty collection if the point is outside the domain. */
Expand Down
8 changes: 0 additions & 8 deletions SKIRT/core/CellSnapshot.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,6 @@
#include "Random.hpp"
#include "StringUtils.hpp"
#include "TextInFile.hpp"
#include "Units.hpp"

////////////////////////////////////////////////////////////////////

Expand Down Expand Up @@ -404,13 +403,6 @@ const Array& CellSnapshot::properties(int m) const

////////////////////////////////////////////////////////////////////

int CellSnapshot::nearestEntity(Position bfr) const
{
return _grid ? _grid->cellIndexFor(bfr) : -1;
}

////////////////////////////////////////////////////////////////////

void CellSnapshot::getEntities(EntityCollection& entities, Position bfr) const
{
entities.addSingle(_grid->cellIndexFor(bfr));
Expand Down
5 changes: 0 additions & 5 deletions SKIRT/core/CellSnapshot.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -109,11 +109,6 @@ class CellSnapshot : public Snapshot
range, the behavior is undefined. */
const Array& properties(int m) const override;

/** This function returns the index \f$0\le m \le N_\mathrm{ent}-1\f$ of the cell containing
the specified point \f${\bf{r}}\f$, or -1 if the point is outside the domain, if there are
no cells in the snapshot, or if the search data structures were not created. */
int nearestEntity(Position bfr) const override;

public:
/** This function sets the specified entity collection to the cell containing the specified
point \f${\bf{r}}\f$, or to the empty collection if the point is outside the domain or if
Expand Down
27 changes: 0 additions & 27 deletions SKIRT/core/ParticleSnapshot.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,6 @@
#include "SmoothingKernel.hpp"
#include "StringUtils.hpp"
#include "TextInFile.hpp"
#include "Units.hpp"

////////////////////////////////////////////////////////////////////

Expand Down Expand Up @@ -212,24 +211,6 @@ class ParticleSnapshot::ParticleGrid : public Box
return _listv[index(_m, i, j, k)];
}

/** This function returns a pointer to the particle centered nearest to the specified position,
or the null pointer if the specified position is outside of the grid. */
const Particle* nearestParticle(Vec r) const
{
const Particle* nearestParticle = nullptr;
double nearestSquaredDistance = std::numeric_limits<double>::infinity();
for (const Particle* particle : particlesFor(r))
{
double d2 = (r - particle->center()).norm2();
if (d2 < nearestSquaredDistance)
{
nearestParticle = particle;
nearestSquaredDistance = d2;
}
}
return nearestParticle;
}

/** This function replaces the contents of the specified entity collection by the set of
particles with a smoothing kernel that overlaps the specified point \f${\bf{r}}\f$. The
weight corresponding to each particle is set to the particle's smoothing kernel value at
Expand Down Expand Up @@ -555,14 +536,6 @@ const Array& ParticleSnapshot::properties(int m) const

////////////////////////////////////////////////////////////////////

int ParticleSnapshot::nearestEntity(Position bfr) const
{
const Particle* nearestParticle = _grid ? _grid->nearestParticle(bfr) : nullptr;
return nearestParticle ? nearestParticle->index() : -1;
}

////////////////////////////////////////////////////////////////////

void ParticleSnapshot::getEntities(EntityCollection& entities, Position bfr) const
{
_grid->getEntities(entities, bfr, _kernel);
Expand Down
6 changes: 0 additions & 6 deletions SKIRT/core/ParticleSnapshot.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -108,12 +108,6 @@ class ParticleSnapshot : public Snapshot
of range, the behavior is undefined. */
const Array& properties(int m) const override;

/** This function returns the index \f$0\le m \le N_\mathrm{ent}-1\f$ of the particle whose
center is nearest to the specified point \f${\bf{r}}\f$, or -1 if the point is outside the
domain, if there are no particles in the snapshot, or if the search data structures were
not created. */
int nearestEntity(Position bfr) const override;

public:
/** This function replaces the contents of the specified entity collection by the set of
particles with a smoothing kernel that overlaps the specified point \f${\bf{r}}\f$. The
Expand Down
70 changes: 29 additions & 41 deletions SKIRT/core/Snapshot.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
///////////////////////////////////////////////////////////////// */

#include "Snapshot.hpp"
#include "EntityCollection.hpp"
#include "Log.hpp"
#include "Random.hpp"
#include "StringUtils.hpp"
Expand Down Expand Up @@ -242,29 +243,13 @@ double Snapshot::initialMass(int m) const

////////////////////////////////////////////////////////////////////

double Snapshot::initialMass(Position bfr) const
{
int m = nearestEntity(bfr);
return m >= 0 ? initialMass(m) : 0.;
}

////////////////////////////////////////////////////////////////////

double Snapshot::currentMass(int m) const
{
return properties(m)[currentMassIndex()];
}

////////////////////////////////////////////////////////////////////

double Snapshot::currentMass(Position bfr) const
{
int m = nearestEntity(bfr);
return m >= 0 ? currentMass(m) : 0.;
}

////////////////////////////////////////////////////////////////////

double Snapshot::metallicity(int m) const
{
return properties(m)[metallicityIndex()];
Expand All @@ -274,8 +259,9 @@ double Snapshot::metallicity(int m) const

double Snapshot::metallicity(Position bfr) const
{
int m = nearestEntity(bfr);
return m >= 0 ? metallicity(m) : 0.;
thread_local EntityCollection entities; // can be reused for all queries in a given execution thread
getEntities(entities, bfr);
return entities.averageValue([this](int m) { return metallicity(m); }, [this](int m) { return currentMass(m); });
}

////////////////////////////////////////////////////////////////////
Expand All @@ -287,14 +273,6 @@ double Snapshot::age(int m) const

////////////////////////////////////////////////////////////////////

double Snapshot::age(Position bfr) const
{
int m = nearestEntity(bfr);
return m >= 0 ? age(m) : 0.;
}

////////////////////////////////////////////////////////////////////

double Snapshot::temperature(int m) const
{
return properties(m)[temperatureIndex()];
Expand All @@ -304,8 +282,9 @@ double Snapshot::temperature(int m) const

double Snapshot::temperature(Position bfr) const
{
int m = nearestEntity(bfr);
return m >= 0 ? temperature(m) : 0.;
thread_local EntityCollection entities; // can be reused for all queries in a given execution thread
getEntities(entities, bfr);
return entities.averageValue([this](int m) { return temperature(m); }, [this](int m) { return currentMass(m); });
}

////////////////////////////////////////////////////////////////////
Expand All @@ -320,8 +299,9 @@ Vec Snapshot::velocity(int m) const

Vec Snapshot::velocity(Position bfr) const
{
int m = nearestEntity(bfr);
return m >= 0 ? velocity(m) : Vec();
thread_local EntityCollection entities; // can be reused for all queries in a given execution thread
getEntities(entities, bfr);
return entities.averageValue([this](int m) { return velocity(m); }, [this](int m) { return currentMass(m); });
}

////////////////////////////////////////////////////////////////////
Expand All @@ -333,14 +313,6 @@ double Snapshot::velocityDispersion(int m) const

////////////////////////////////////////////////////////////////////

double Snapshot::velocityDispersion(Position bfr) const
{
int m = nearestEntity(bfr);
return m >= 0 ? velocityDispersion(m) : 0.;
}

////////////////////////////////////////////////////////////////////

Vec Snapshot::magneticField(int m) const
{
const auto& propv = properties(m);
Expand All @@ -351,8 +323,9 @@ Vec Snapshot::magneticField(int m) const

Vec Snapshot::magneticField(Position bfr) const
{
int m = nearestEntity(bfr);
return m >= 0 ? magneticField(m) : Vec();
thread_local EntityCollection entities; // can be reused for all queries in a given execution thread
getEntities(entities, bfr);
return entities.averageValue([this](int m) { return magneticField(m); }, [this](int m) { return currentMass(m); });
}

////////////////////////////////////////////////////////////////////
Expand All @@ -369,7 +342,22 @@ void Snapshot::parameters(int m, Array& params) const

void Snapshot::parameters(Position bfr, Array& params) const
{
int m = nearestEntity(bfr);
thread_local EntityCollection entities; // can be reused for all queries in a given execution thread
getEntities(entities, bfr);

// look for the entity with the highest weight
double wmax = 0.;
int m = -1;
for (const auto& entity : entities)
{
if (entity.second > wmax)
{
wmax = entity.second;
m = entity.first;
}
}

// if found, use it
if (m >= 0)
parameters(m, params);
else
Expand Down
38 changes: 5 additions & 33 deletions SKIRT/core/Snapshot.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -359,13 +359,6 @@ class Snapshot
range, the behavior is undefined. */
virtual const Array& properties(int m) const = 0;

/** This function returns the index \f$0\le m \le N_\mathrm{ent}-1\f$ for the entity at or
nearest to the specified point \f${\bf{r}}\f$, or -1 if the point is outside the domain or
if there are no entities in the snapshot. For a cell-based snapshot, the function returns
the index of the cell containing the given point. For a particle-based snapshot, the
function returns the index of the particle whose center is nearest to the given point. */
virtual int nearestEntity(Position bfr) const = 0;

public:
/** This function replaces the contents of the specified entity collection by the set of
entities that overlap the specified point \f${\bf{r}}\f$, with their corresponding weights.
Expand Down Expand Up @@ -407,11 +400,6 @@ class Snapshot
range, the behavior is undefined. */
double initialMass(int m) const;

/** This function returns the initial mass of the entity nearest to (or at) the specified point
\f${\bf{r}}\f$. If the point is outside the domain, the function returns zero. If the
initial mass is not being imported, the behavior is undefined. */
double initialMass(Position bfr) const;

/** This function returns true if the current mass is being imported, and false otherwise. */
bool hasCurrentMass() const { return _currentMassIndex >= 0; }

Expand All @@ -420,11 +408,6 @@ class Snapshot
range, the behavior is undefined. */
double currentMass(int m) const;

/** This function returns the current mass of the entity nearest to (or at) the specified point
\f${\bf{r}}\f$. If the point is outside the domain, the function returns zero. If the
current mass is not being imported, the behavior is undefined. */
double currentMass(Position bfr) const;

/** This function returns true if the metallicity is being imported, and false otherwise. */
bool hasMetallicity() const { return _metallicityIndex >= 0; }

Expand All @@ -433,7 +416,7 @@ class Snapshot
range, the behavior is undefined. */
double metallicity(int m) const;

/** This function returns the metallicity of the entity nearest to (or at) the specified point
/** This function returns the metallicity at the specified point
\f${\bf{r}}\f$. If the point is outside the domain, the function returns zero. If the
metallicity is not being imported, the behavior is undefined. */
double metallicity(Position bfr) const;
Expand All @@ -446,11 +429,6 @@ class Snapshot
*/
double age(int m) const;

/** This function returns the age of the entity nearest to (or at) the specified point
\f${\bf{r}}\f$. If the point is outside the domain, the function returns zero. If the age
is not being imported, the behavior is undefined. */
double age(Position bfr) const;

/** This function returns true if the temperature is being imported, and false otherwise. */
bool hasTemperature() const { return _temperatureIndex >= 0; }

Expand All @@ -459,7 +437,7 @@ class Snapshot
range, the behavior is undefined. */
double temperature(int m) const;

/** This function returns the temperature of the entity nearest to (or at) the specified point
/** This function returns the temperature at the specified point
\f${\bf{r}}\f$. If the point is outside the domain, the function returns zero. If the
temperature is not being imported, the behavior is undefined. */
double temperature(Position bfr) const;
Expand All @@ -472,7 +450,7 @@ class Snapshot
the behavior is undefined. */
Vec velocity(int m) const;

/** This function returns the velocity of the entity nearest to (or at) the specified point
/** This function returns the velocity at the specified point
\f${\bf{r}}\f$. If the point is outside the domain, the function returns zero velocity. If
the velocity is not being imported, the behavior is undefined. */
Vec velocity(Position bfr) const;
Expand All @@ -486,12 +464,6 @@ class Snapshot
of range, the behavior is undefined. */
double velocityDispersion(int m) const;

/** This function returns the velocity dispersion of the entity nearest to (or at) the
specified point \f${\bf{r}}\f$. If the point is outside the domain, the function returns
zero dispersion. If the velocity dispersion is not being imported, the behavior is
undefined. */
double velocityDispersion(Position bfr) const;

/** This function returns true if the magnetic field is being imported, and false otherwise. */
bool hasMagneticField() const { return _magneticFieldIndex >= 0; }

Expand All @@ -500,7 +472,7 @@ class Snapshot
range, the behavior is undefined. */
Vec magneticField(int m) const;

/** This function returns the magnetic field vector of the entity nearest to (or at) the
/** This function returns the magnetic field vector at the
specified point \f${\bf{r}}\f$. If the point is outside the domain, the function returns a
zero magnetic field. If the magnetic field is not being imported, the behavior is
undefined. */
Expand All @@ -515,7 +487,7 @@ class Snapshot
index is out of range, the behavior is undefined. */
void parameters(int m, Array& params) const;

/** This function stores the parameters of the entity nearest to (or at) the specified point
/** This function stores the parameters at the specified point
\f${\bf{r}}\f$ into the given array. If the point is outside the domain, the function
returns the appropriate number of zero parameter values. If parameters are not being
imported, the behavior is undefined. */
Expand Down
8 changes: 0 additions & 8 deletions SKIRT/core/VoronoiMeshSnapshot.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,6 @@
#include "StringUtils.hpp"
#include "Table.hpp"
#include "TextInFile.hpp"
#include "Units.hpp"
#include <set>
#include "container.hh"

Expand Down Expand Up @@ -1108,13 +1107,6 @@ const Array& VoronoiMeshSnapshot::properties(int m) const

////////////////////////////////////////////////////////////////////

int VoronoiMeshSnapshot::nearestEntity(Position bfr) const
{
return _blocktrees.size() ? cellIndex(bfr) : -1;
}

////////////////////////////////////////////////////////////////////

void VoronoiMeshSnapshot::getEntities(EntityCollection& entities, Position bfr) const
{
entities.addSingle(cellIndex(bfr));
Expand Down
Loading