Merge branch 'develop' into feature/han12/cuda_tutorial_image

src/axom/core/Array.hpp

@@ -625,13 +625,34 @@ class Array : public ArrayBase<T, DIM, Array<T, DIM, SPACE>>
    * \param [in] value the value to be added to the back.
    *
    * \note Reallocation is done if the new size will exceed the capacity.
-   * \note If used in a device kernel, the number of push_backs must not exceed
-   *  the capacity, since device-side reallocations aren't supported.
-   * \note Array must be allocated in unified memory if calling on the device.
-   * 
+   *
+   * \pre DIM == 1
+   */
+  void push_back(const T& value);
+
+  /*!
+   * \brief Push a value to the back of the array.
+   *
+   * \param [in] value the value to move to the back.
+   *
+   * \note Reallocation is done if the new size will exceed the capacity.
+   *
+   * \pre DIM == 1
+   */
+  void push_back(T&& value);
+
+  /*!
+   * \brief Inserts new element at the end of the Array.
+   *
+   * \param [in] args the arguments to forward to constructor of the element.
+   *
+   * \note Reallocation is done if the new size will exceed the capacity.
+   * \note The size increases by 1.
+   *
    * \pre DIM == 1
    */
-  AXOM_HOST_DEVICE void push_back(const T& value);
+  template <typename... Args>
+  void emplace_back(Args&&... args);
 
   /*!
    * \brief Push a value to the back of the array.
@@ -642,10 +663,13 @@ class Array : public ArrayBase<T, DIM, Array<T, DIM, SPACE>>
    * \note If used in a device kernel, the number of push_backs must not exceed
    *  the capacity, since device-side reallocations aren't supported.
    * \note Array must be allocated in unified memory if calling on the device.
-   * 
+   *
    * \pre DIM == 1
    */
-  AXOM_HOST_DEVICE void push_back(T&& value);
+  /// @{
+  AXOM_HOST_DEVICE void push_back_device(const T& value);
+  AXOM_HOST_DEVICE void push_back_device(T&& value);
+  /// @}
 
   /*!
    * \brief Inserts new element at the end of the Array.
@@ -657,11 +681,11 @@ class Array : public ArrayBase<T, DIM, Array<T, DIM, SPACE>>
    * \note If used in a device kernel, the number of push_backs must not exceed
    *  the capacity, since device-side reallocations aren't supported.
    * \note Array must be allocated in unified memory if calling on the device.
-   * 
+   *
    * \pre DIM == 1
    */
   template <typename... Args>
-  AXOM_HOST_DEVICE void emplace_back(Args&&... args);
+  AXOM_HOST_DEVICE void emplace_back_device(Args&&... args);
 
   /// @}
 
@@ -1425,15 +1449,15 @@ inline typename Array<T, DIM, SPACE>::ArrayIterator Array<T, DIM, SPACE>::emplac
 
 //------------------------------------------------------------------------------
 template <typename T, int DIM, MemorySpace SPACE>
-AXOM_HOST_DEVICE inline void Array<T, DIM, SPACE>::push_back(const T& value)
+inline void Array<T, DIM, SPACE>::push_back(const T& value)
 {
   static_assert(DIM == 1, "push_back is only supported for 1D arrays");
   emplace_back(value);
 }
 
 //------------------------------------------------------------------------------
 template <typename T, int DIM, MemorySpace SPACE>
-AXOM_HOST_DEVICE inline void Array<T, DIM, SPACE>::push_back(T&& value)
+inline void Array<T, DIM, SPACE>::push_back(T&& value)
 {
   static_assert(DIM == 1, "push_back is only supported for 1D arrays");
   emplace_back(std::move(value));
@@ -1442,8 +1466,33 @@ AXOM_HOST_DEVICE inline void Array<T, DIM, SPACE>::push_back(T&& value)
 //------------------------------------------------------------------------------
 AXOM_SUPPRESS_HD_WARN
 template <typename T, int DIM, MemorySpace SPACE>
+AXOM_HOST_DEVICE inline void Array<T, DIM, SPACE>::push_back_device(const T& value)
+{
+  static_assert(DIM == 1, "push_back_device is only supported for 1D arrays");
+  emplace_back_device(value);
+}
+
+//------------------------------------------------------------------------------
+template <typename T, int DIM, MemorySpace SPACE>
+AXOM_HOST_DEVICE inline void Array<T, DIM, SPACE>::push_back_device(T&& value)
+{
+  static_assert(DIM == 1, "push_back_device is only supported for 1D arrays");
+  emplace_back_device(std::move(value));
+}
+
+//------------------------------------------------------------------------------
+template <typename T, int DIM, MemorySpace SPACE>
 template <typename... Args>
-AXOM_HOST_DEVICE inline void Array<T, DIM, SPACE>::emplace_back(Args&&... args)
+inline void Array<T, DIM, SPACE>::emplace_back(Args&&... args)
+{
+  static_assert(DIM == 1, "emplace_back is only supported for 1D arrays");
+  emplace(size(), std::forward<Args>(args)...);
+}
+
+//------------------------------------------------------------------------------
+template <typename T, int DIM, MemorySpace SPACE>
+template <typename... Args>
+AXOM_HOST_DEVICE inline void Array<T, DIM, SPACE>::emplace_back_device(Args&&... args)
 {
   static_assert(DIM == 1, "emplace_back is only supported for 1D arrays");
 #ifdef AXOM_DEVICE_CODE
@@ -1685,7 +1734,12 @@ template <typename T, int DIM, MemorySpace SPACE>
 inline void Array<T, DIM, SPACE>::dynamicRealloc(IndexType new_num_elements)
 {
   assert(m_resize_ratio >= 1.0);
-  IndexType new_capacity = new_num_elements * m_resize_ratio + 0.5;
+
+  // Using resize strategy from LLVM libc++ (vector::__recommend()):
+  //   new_capacity = max(capacity() * resize_ratio, new_num_elements)
+  IndexType new_capacity =
+    axom::utilities::max<IndexType>(this->capacity() * m_resize_ratio + 0.5,
+                                    new_num_elements);
   const IndexType block_size = this->blockSize();
   const IndexType remainder = new_capacity % block_size;
   if(remainder != 0)

src/axom/core/MDMapping.hpp

@@ -49,9 +49,9 @@ class MDMapping
     @param [in] fastestStrideLength Stride in the
       fastest-changing direction.
   */
-  MDMapping(const axom::StackArray<T, DIM>& shape,
-            axom::ArrayStrideOrder arrayStrideOrder,
-            int fastestStrideLength = 1)
+  AXOM_HOST_DEVICE MDMapping(const axom::StackArray<T, DIM>& shape,
+                             axom::ArrayStrideOrder arrayStrideOrder,
+                             int fastestStrideLength = 1)
   {
     initializeShape(shape, arrayStrideOrder, fastestStrideLength);
   }
@@ -66,9 +66,9 @@ class MDMapping
       fastest-changing direction.
   */
   template <typename DirType>
-  MDMapping(const axom::StackArray<T, DIM>& shape,
-            const axom::StackArray<DirType, DIM>& slowestDirs,
-            int fastestStrideLength = 1)
+  AXOM_HOST_DEVICE MDMapping(const axom::StackArray<T, DIM>& shape,
+                             const axom::StackArray<DirType, DIM>& slowestDirs,
+                             int fastestStrideLength = 1)
   {
     initializeShape(shape, slowestDirs, fastestStrideLength);
   }
@@ -81,8 +81,8 @@ class MDMapping
     @param [in] orderSource ArrayIndex to copy stride order
       from.
   */
-  MDMapping(const axom::StackArray<T, DIM>& shape,
-            const axom::MDMapping<DIM, T>& orderSource)
+  AXOM_HOST_DEVICE MDMapping(const axom::StackArray<T, DIM>& shape,
+                             const axom::MDMapping<DIM, T>& orderSource)
   {
     initializeShape(shape, orderSource);
   }
@@ -98,7 +98,8 @@ class MDMapping
     clash with the more prevalent usage of constructing from the array's
     shape.
   */
-  MDMapping(const axom::StackArray<T, DIM>& strides) : m_strides(strides)
+  AXOM_HOST_DEVICE MDMapping(const axom::StackArray<T, DIM>& strides)
+    : m_strides(strides)
   {
     initializeStrides(strides);
   }
@@ -117,7 +118,7 @@ class MDMapping
     shape is to be determined.  Initialize the mapping with its own
     slowestDirs() to preserve stride order as its shape changes.
   */
-  MDMapping(ArrayStrideOrder arrayStrideOrder)
+  AXOM_HOST_DEVICE MDMapping(ArrayStrideOrder arrayStrideOrder)
   {
     axom::StackArray<T, DIM> shape;
     for(int d = 0; d < DIM; ++d)
@@ -238,8 +239,8 @@ class MDMapping
                 << "Likely, multi-dim array shape is 1 in some direction.\n"
                 << "Impossible to compute index ordering.\n"
                 << "Please use a different MDMapping initializer.\n";
-#endif
       utilities::processAbort();
+#endif
     }
 
     // 2nd argument doesn't matter because strides are unique.
@@ -275,7 +276,7 @@ class MDMapping
       {
         if(m_strides[m_slowestDirs[s]] < m_strides[m_slowestDirs[d]])
         {
-          std::swap(m_slowestDirs[s], m_slowestDirs[d]);
+          axom::utilities::swap(m_slowestDirs[s], m_slowestDirs[d]);
         }
       }
     }

src/axom/core/tests/core_array.hpp

@@ -25,7 +25,9 @@ axom::IndexType calc_new_capacity(axom::Array<T>& v, axom::IndexType increase)
   axom::IndexType new_num_elements = v.size() + increase;
   if(new_num_elements > v.capacity())
   {
-    return new_num_elements * v.getResizeRatio() + 0.5;
+    axom::IndexType capacity_expanded = v.capacity() * v.getResizeRatio() + 0.5;
+    return axom::utilities::max<axom::IndexType>(capacity_expanded,
+                                                 new_num_elements);
   }
 
   return v.capacity();

src/axom/core/tests/core_array_for_all.hpp

@@ -1180,15 +1180,15 @@ AXOM_TYPED_TEST(core_array_for_all, device_insert)
       {
   #pragma omp critical
         {
-          arr_v[0].emplace_back(3 * idx + 5);
+          arr_v[0].emplace_back_device(3 * idx + 5);
         }
       }
       else
       {
-        arr_v[0].emplace_back(3 * idx + 5);
+        arr_v[0].emplace_back_device(3 * idx + 5);
       }
 #else
-      arr_v[0].emplace_back(3 * idx + 5);
+      arr_v[0].emplace_back_device(3 * idx + 5);
 #endif
     });
 

src/axom/mint/execution/internal/for_all_faces.hpp

@@ -993,7 +993,7 @@ inline void for_all_faces_impl(xargs::coords,
           coords[3 * i + 2] = z[nodeID];
         }
 
-        numerics::Matrix<double> coordsMatrix(dimension, 4, coords, NO_COPY);
+        numerics::Matrix<double> coordsMatrix(dimension, numNodes, coords, NO_COPY);
         kernel(faceID, coordsMatrix, nodeIDs);
       });
   }

src/axom/mint/tests/mint_mesh_connectivity_array.cpp

@@ -41,8 +41,9 @@ IndexType calc_ID_capacity(const ConnectivityArray<NO_INDIRECTION>& connec,
   if(new_n_IDs > connec.getIDCapacity())
   {
     IndexType stride = connec.getNumberOfValuesForID();
-    IndexType newCapacity =
-      static_cast<IndexType>(new_n_IDs * stride * connec.getResizeRatio() + 0.5);
+    IndexType newCapacity = axom::utilities::max<IndexType>(
+      connec.getValueCapacity() * connec.getResizeRatio() + 0.5,
+      new_n_IDs * stride);
     IndexType remainder = newCapacity % stride;
     if(remainder != 0)
     {
@@ -67,7 +68,9 @@ IndexType calc_ID_capacity(const ConnectivityArray<TYPED_INDIRECTION>& connec,
   IndexType new_n_IDs = connec.getNumberOfIDs() + increase;
   if(new_n_IDs > connec.getIDCapacity())
   {
-    return static_cast<IndexType>(new_n_IDs * connec.getResizeRatio() + 0.5);
+    return axom::utilities::max<IndexType>(
+      connec.getIDCapacity() * connec.getResizeRatio() + 0.5,
+      new_n_IDs);
   }
 
   return connec.getIDCapacity();
@@ -101,7 +104,9 @@ IndexType calc_value_capacity(const ConnectivityArray<TYPE>& connec,
   IndexType new_n_values = connec.getNumberOfValues() + increase;
   if(new_n_values > connec.getValueCapacity())
   {
-    return static_cast<IndexType>(new_n_values * connec.getResizeRatio() + 0.5);
+    return axom::utilities::max<IndexType>(
+      connec.getValueCapacity() * connec.getResizeRatio() + 0.5,
+      new_n_values);
   }
 
   return connec.getValueCapacity();

src/axom/mint/tests/mint_mesh_unstructured_mesh.cpp

@@ -2200,7 +2200,8 @@ void resize_cells(UnstructuredMesh<TOPO>* mesh)
   /* Append one more, should trigger a resize. */
   append_cell_single(mesh, 1);
   n_cells++;
-  cell_capacity = static_cast<IndexType>(n_cells * resize_ratio + 0.5);
+  cell_capacity =
+    axom::utilities::max<IndexType>(cell_capacity * resize_ratio + 0.5, n_cells);
   connec_capacity = mesh->getCellNodesCapacity();
   ASSERT_EQ(n_cells, mesh->getNumberOfCells());
   ASSERT_EQ(cell_capacity, mesh->getCellCapacity());
@@ -2217,7 +2218,8 @@ void resize_cells(UnstructuredMesh<TOPO>* mesh)
   mesh->setCellResizeRatio(resize_ratio);
   append_cell_multiple(mesh, 100);
   n_cells += 100;
-  cell_capacity = static_cast<IndexType>(resize_ratio * n_cells + 0.5);
+  cell_capacity =
+    axom::utilities::max<IndexType>(cell_capacity * resize_ratio + 0.5, n_cells);
   ASSERT_EQ(n_cells, mesh->getNumberOfCells());
   ASSERT_EQ(cell_capacity, mesh->getCellCapacity());
 

src/axom/primal/geometry/Polygon.hpp

@@ -193,8 +193,25 @@ class Polygon
    * \sa axom::StaticArray::push_back() for behavior when array type is static
    *     and the list of vertices is full.
    */
-  AXOM_HOST_DEVICE
-  void addVertex(const PointType& pt) { m_vertices.push_back(pt); }
+  /// @{
+  template <PolygonArray P_ARRAY_TYPE = ARRAY_TYPE,
+            std::enable_if_t<P_ARRAY_TYPE == PolygonArray::Static, int> = 0>
+  AXOM_HOST_DEVICE void addVertex(const PointType& pt)
+  {
+    m_vertices.push_back(pt);
+  }
+
+  template <PolygonArray P_ARRAY_TYPE = ARRAY_TYPE,
+            std::enable_if_t<P_ARRAY_TYPE == PolygonArray::Dynamic, int> = 0>
+  AXOM_HOST_DEVICE void addVertex(const PointType& pt)
+  {
+#ifdef AXOM_DEVICE_CODE
+    m_vertices.push_back_device(pt);
+#else
+    m_vertices.push_back(pt);
+#endif
+  }
+  /// @}
 
   /// Clears the list of vertices (dynamic array specialization).
   /// Specializations are necessary to remove __host__ __device__ warning for

src/axom/sidre/core/Array.hpp

@@ -571,7 +571,9 @@ inline void Array<T, DIM>::dynamicRealloc(axom::IndexType new_num_elements)
                 "Resize ratio of " << this->m_resize_ratio
                                    << " doesn't support dynamic resizing");
 
-  IndexType new_capacity = new_num_elements * this->m_resize_ratio + 0.5;
+  IndexType new_capacity = axom::utilities::max<IndexType>(
+    this->capacity() * this->getResizeRatio() + 0.5,
+    new_num_elements);
   const IndexType block_size = this->blockSize();
   const IndexType remainder = new_capacity % block_size;
   if(remainder != 0)

src/axom/sidre/tests/sidre_mcarray.cpp

@@ -34,7 +34,14 @@ axom::IndexType calc_new_capacity(MCArray<T>& v, axom::IndexType increase)
   axom::IndexType new_num_tuples = (v.size() / num_components) + increase;
   if((new_num_tuples * num_components) > v.capacity())
   {
-    return new_num_tuples * v.getResizeRatio() + 0.5;
+    axom::IndexType new_capacity = v.capacity() * v.getResizeRatio() + 0.5;
+    axom::IndexType remainder = new_capacity % num_components;
+    if(remainder > 0)
+    {
+      new_capacity += num_components - remainder;
+    }
+    return axom::utilities::max<axom::IndexType>(new_capacity / num_components,
+                                                 new_num_tuples);
   }
 
   return v.capacity() / num_components;