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 1ce02000a21ae84d8fd96e0381729a78e88ea4d7..b696e206b22d57e02256cdf43b1738441be78297 100644 (file)
--- a/include/deal.II/particles/property_pool.h
+++ b/include/deal.II/particles/property_pool.h
@@ -19,6 +19,7 @@
 #include <deal.II/base/config.h>
 
 #include <deal.II/base/array_view.h>
+#include <deal.II/base/point.h>
 
 
 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<spacedim> &
+    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<spacedim> &new_location);
+
+    /**
+     * Return the reference_location of a particle identified by the given
+     * `handle`.
+     */
+    const Point<dim> &
+    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<dim> &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<Point<spacedim>> 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<Point<dim>> 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<double> properties;
 
@@ -142,6 +188,101 @@ namespace Particles
 
   /* ---------------------- inline and template functions ------------------ */
 
+  template <int dim, int spacedim>
+  inline const Point<spacedim> &
+  PropertyPool<dim, spacedim>::get_location(const Handle handle) const
+  {
+    const std::vector<double>::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 <int dim, int spacedim>
+  inline void
+  PropertyPool<dim, spacedim>::set_location(const Handle           handle,
+                                            const Point<spacedim> &new_location)
+  {
+    const std::vector<double>::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 <int dim, int spacedim>
+  inline const Point<dim> &
+  PropertyPool<dim, spacedim>::get_reference_location(const Handle handle) const
+  {
+    const std::vector<double>::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 <int dim, int spacedim>
+  inline void
+  PropertyPool<dim, spacedim>::set_reference_location(
+    const Handle      handle,
+    const Point<dim> &new_reference_location)
+  {
+    const std::vector<double>::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 <int dim, int spacedim>
   inline ArrayView<double>
   PropertyPool<dim, spacedim>::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 1637aafab22876856e8fd1d030b2959c85acab36..68e2f2f8ff960eb25681eae3a54a28e9ffef793c 100644 (file)
--- a/source/particles/property_pool.cc
+++ b/source/particles/property_pool.cc
@@ -14,6 +14,8 @@
 // ---------------------------------------------------------------------
 
 
+#include <deal.II/base/signaling_nan.h>
+
 #include <deal.II/particles/property_pool.h>
 
 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<Point<spacedim>>());
+            reference_locations.resize(reference_locations.size() + 1,
+                                       numbers::signaling_nan<Point<dim>>());
+
+            // In contrast, initialize properties by zero.
             properties.resize(properties.size() + n_properties, 0.0);
           }
       }
@@ -122,6 +139,8 @@ namespace Particles
   void
   PropertyPool<dim, spacedim>::reserve(const std::size_t size)
   {
+    locations.reserve(size);
+    reference_locations.reserve(size);
     properties.reserve(size * n_properties);
   }