From: Wolfgang Bangerth Date: Sat, 12 Dec 2020 17:31:25 +0000 (-0700) Subject: Provide the interface to also store locations of particles in the PropertyPool. X-Git-Tag: v9.3.0-rc1~625^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3e9a7e5e403c66d5dcedad84f4b97998125be444;p=dealii.git Provide the interface to also store locations of particles in the PropertyPool. --- diff --git a/include/deal.II/particles/property_pool.h b/include/deal.II/particles/property_pool.h index 1ce02000a2..b696e206b2 100644 --- a/include/deal.II/particles/property_pool.h +++ b/include/deal.II/particles/property_pool.h @@ -19,6 +19,7 @@ #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -27,7 +28,10 @@ namespace Particles { /** * This class manages a memory space in which all particles associated with - * a ParticleHandler store their properties. The rationale for this class is + * a ParticleHandler store their properties. It also stores the locations + * and reference locations of particles. + * + * The rationale for this class is * that because typically every particle stores the same number of * properties, and because algorithms generally traverse over all particles * doing the same operation on all particles' properties, it is more efficient @@ -96,6 +100,33 @@ namespace Particles void deregister_particle(Handle &handle); + /** + * Return the location of a particle identified by the given `handle`. + */ + const Point & + get_location(const Handle handle) const; + + /** + * Set the location of a particle identified by the given `handle`. + */ + void + set_location(const Handle handle, const Point &new_location); + + /** + * Return the reference_location of a particle identified by the given + * `handle`. + */ + const Point & + get_reference_location(const Handle handle) const; + + /** + * Set the reference location of a particle identified by the given + * `handle`. + */ + void + set_reference_location(const Handle handle, + const Point &new_reference_location); + /** * Return an ArrayView to the properties that correspond to the given * handle @p handle. @@ -122,9 +153,24 @@ namespace Particles */ const unsigned int n_properties; + /** + * A vector that stores the locations of particles. It is indexed in the + * same way as the `reference_locations` and `properties` arrays, i.e., via + * handles. + */ + std::vector> locations; + + /** + * A vector that stores the reference locations of particles. It is indexed + * in the same way as the `locations` and `properties` arrays, i.e., via + * handles. + */ + std::vector> reference_locations; + /** * The currently allocated properties (whether assigned to - * a particle or available for assignment). + * a particle or available for assignment). It is indexed the same way as + * the `locations` and `reference_locations` arrays via handles. */ std::vector properties; @@ -142,6 +188,101 @@ namespace Particles /* ---------------------- inline and template functions ------------------ */ + template + inline const Point & + PropertyPool::get_location(const Handle handle) const + { + const std::vector::size_type data_index = + (handle != invalid_handle) ? handle : 0; + + // Ideally we would need to assert that 'handle' has not been deallocated + // by searching through 'currently_available_handles'. However, this + // is expensive and this function is performance critical, so instead + // just check against the array range, and rely on the fact + // that handles are invalidated when handed over to + // deallocate_properties_array(). + Assert(data_index <= locations.size() - 1, + ExcMessage("Invalid location handle. This can happen if the " + "handle was duplicated and then one copy was deallocated " + "before trying to access the properties.")); + + return locations[data_index]; + } + + + + template + inline void + PropertyPool::set_location(const Handle handle, + const Point &new_location) + { + const std::vector::size_type data_index = + (handle != invalid_handle) ? handle : 0; + + // Ideally we would need to assert that 'handle' has not been deallocated + // by searching through 'currently_available_handles'. However, this + // is expensive and this function is performance critical, so instead + // just check against the array range, and rely on the fact + // that handles are invalidated when handed over to + // deallocate_properties_array(). + Assert(data_index <= locations.size() - 1, + ExcMessage("Invalid location handle. This can happen if the " + "handle was duplicated and then one copy was deallocated " + "before trying to access the properties.")); + + locations[data_index] = new_location; + } + + + + template + inline const Point & + PropertyPool::get_reference_location(const Handle handle) const + { + const std::vector::size_type data_index = + (handle != invalid_handle) ? handle : 0; + + // Ideally we would need to assert that 'handle' has not been deallocated + // by searching through 'currently_available_handles'. However, this + // is expensive and this function is performance critical, so instead + // just check against the array range, and rely on the fact + // that handles are invalidated when handed over to + // deallocate_properties_array(). + Assert(data_index <= reference_locations.size() - 1, + ExcMessage("Invalid location handle. This can happen if the " + "handle was duplicated and then one copy was deallocated " + "before trying to access the properties.")); + + return reference_locations[data_index]; + } + + + + template + inline void + PropertyPool::set_reference_location( + const Handle handle, + const Point &new_reference_location) + { + const std::vector::size_type data_index = + (handle != invalid_handle) ? handle : 0; + + // Ideally we would need to assert that 'handle' has not been deallocated + // by searching through 'currently_available_handles'. However, this + // is expensive and this function is performance critical, so instead + // just check against the array range, and rely on the fact + // that handles are invalidated when handed over to + // deallocate_properties_array(). + Assert(data_index <= locations.size() - 1, + ExcMessage("Invalid location handle. This can happen if the " + "handle was duplicated and then one copy was deallocated " + "before trying to access the properties.")); + + reference_locations[data_index] = new_reference_location; + } + + + template inline ArrayView PropertyPool::get_properties(const Handle handle) @@ -152,7 +293,7 @@ namespace Particles // Ideally we would need to assert that 'handle' has not been deallocated // by searching through 'currently_available_handles'. However, this // is expensive and this function is performance critical, so instead - // just check against the properties range, and rely on the fact + // just check against the array range, and rely on the fact // that handles are invalidated when handed over to // deallocate_properties_array(). Assert(data_index <= properties.size() - n_properties, diff --git a/source/particles/property_pool.cc b/source/particles/property_pool.cc index 1637aafab2..68e2f2f8ff 100644 --- a/source/particles/property_pool.cc +++ b/source/particles/property_pool.cc @@ -14,6 +14,8 @@ // --------------------------------------------------------------------- +#include + #include DEAL_II_NAMESPACE_OPEN @@ -61,6 +63,12 @@ namespace Particles } // Clear vectors and ensure deallocation of memory + locations.clear(); + locations.shrink_to_fit(); + + reference_locations.clear(); + reference_locations.shrink_to_fit(); + properties.clear(); properties.shrink_to_fit(); @@ -84,7 +92,16 @@ namespace Particles } else { - handle = properties.size() / n_properties; + handle = locations.size(); + + // Append new slots to the end of the arrays. Initialize with + // invalid points. + locations.resize(locations.size() + 1, + numbers::signaling_nan>()); + reference_locations.resize(reference_locations.size() + 1, + numbers::signaling_nan>()); + + // In contrast, initialize properties by zero. properties.resize(properties.size() + n_properties, 0.0); } } @@ -122,6 +139,8 @@ namespace Particles void PropertyPool::reserve(const std::size_t size) { + locations.reserve(size); + reference_locations.reserve(size); properties.reserve(size * n_properties); }