diff --git a/src/axom/primal/geometry/Point.hpp b/src/axom/primal/geometry/Point.hpp
index e48096dd62..c8033877ce 100644
--- a/src/axom/primal/geometry/Point.hpp
+++ b/src/axom/primal/geometry/Point.hpp
@@ -37,7 +37,8 @@ AXOM_HOST_DEVICE bool operator==(const Point<T, NDIMS>& lhs,
  * \brief Inequality comparison operator for points
  */
 template <typename T, int NDIMS>
-bool operator!=(const Point<T, NDIMS>& lhs, const Point<T, NDIMS>& rhs);
+AXOM_HOST_DEVICE bool operator!=(const Point<T, NDIMS>& lhs,
+                                 const Point<T, NDIMS>& rhs);
 
 /*!
  * \brief Overloaded output operator for points
@@ -113,7 +114,8 @@ class Point
    * \return d the dimension of the point.
    * \post d >= 1.
    */
-  static int dimension() { return NDIMS; };
+  AXOM_HOST_DEVICE
+  static constexpr int dimension() { return NDIMS; };
 
   /// \name Overloaded [] operator methods
   ///@{
diff --git a/src/axom/primal/geometry/Segment.hpp b/src/axom/primal/geometry/Segment.hpp
index 9c787a2889..f363b311a5 100644
--- a/src/axom/primal/geometry/Segment.hpp
+++ b/src/axom/primal/geometry/Segment.hpp
@@ -28,13 +28,15 @@ class Segment;
  * \brief Equality comparison operator for Segment
  */
 template <typename T, int NDIMS>
-bool operator==(const Segment<T, NDIMS>& lhs, const Segment<T, NDIMS>& rhs);
+AXOM_HOST_DEVICE bool operator==(const Segment<T, NDIMS>& lhs,
+                                 const Segment<T, NDIMS>& rhs);
 
 /*!
  * \brief Inequality comparison operator for Segment
  */
 template <typename T, int NDIMS>
-bool operator!=(const Segment<T, NDIMS>& lhs, const Segment<T, NDIMS>& rhs);
+AXOM_HOST_DEVICE bool operator!=(const Segment<T, NDIMS>& lhs,
+                                 const Segment<T, NDIMS>& rhs);
 
 /*!
  * \brief Overloaded output operator for Segment
@@ -130,7 +132,10 @@ class Segment
   /*!
    * \brief Returns the length of the segment
    */
-  double length() const { return VectorType(m_source, m_target).norm(); }
+  AXOM_HOST_DEVICE double length() const
+  {
+    return VectorType(m_source, m_target).norm();
+  }
 
   /*!
    * \brief Returns a vector normal to the segment
@@ -138,7 +143,7 @@ class Segment
    * \note Only available in 2D
    */
   template <int TDIM = NDIMS>
-  typename std::enable_if<TDIM == 2, VectorType>::type normal() const
+  AXOM_HOST_DEVICE typename std::enable_if<TDIM == 2, VectorType>::type normal() const
   {
     return VectorType {m_target[1] - m_source[1], m_source[0] - m_target[0]};
   }
@@ -146,6 +151,7 @@ class Segment
   /*!
    * \brief Equality comparison operator for segments
    */
+  AXOM_HOST_DEVICE
   friend inline bool operator==(const Segment& lhs, const Segment& rhs)
   {
     return lhs.m_source == rhs.m_source && lhs.m_target == rhs.m_target;
@@ -154,6 +160,7 @@ class Segment
   /*!
    * \brief Inequality operator for segments
    */
+  AXOM_HOST_DEVICE
   friend inline bool operator!=(const Segment& lhs, const Segment& rhs)
   {
     return !(lhs == rhs);
diff --git a/src/axom/primal/geometry/Vector.hpp b/src/axom/primal/geometry/Vector.hpp
index 9c00380b7f..baea921f8e 100644
--- a/src/axom/primal/geometry/Vector.hpp
+++ b/src/axom/primal/geometry/Vector.hpp
@@ -101,7 +101,7 @@ AXOM_HOST_DEVICE Vector<T, NDIMS> operator-(const Point<T, NDIMS>& h,
  * \return C resulting vector from unary negation.
  */
 template <typename T, int NDIMS>
-Vector<T, NDIMS> operator-(const Vector<T, NDIMS>& vec1);
+AXOM_HOST_DEVICE Vector<T, NDIMS> operator-(const Vector<T, NDIMS>& vec1);
 
 /*!
  * \brief Scalar multiplication of vector; Scalar on rhs.
@@ -131,7 +131,8 @@ AXOM_HOST_DEVICE Vector<T, NDIMS> operator*(const T scalar,
  * \pre scalar != 0.0
  */
 template <typename T, int NDIMS>
-Vector<T, NDIMS> operator/(const Vector<T, NDIMS>& vec, const T scalar);
+AXOM_HOST_DEVICE Vector<T, NDIMS> operator/(const Vector<T, NDIMS>& vec,
+                                            const T scalar);
 
 /*!
  * \brief Overloaded output operator for vectors
@@ -229,7 +230,7 @@ class Vector
    * \post d >= 1.
    */
   AXOM_HOST_DEVICE
-  int dimension() const { return NDIMS; };
+  static constexpr int dimension() { return NDIMS; };
 
   /*!
    * \brief Access operator for individual components.
@@ -273,6 +274,7 @@ class Vector
   /*!
    * \brief Inequality operator for points
    */
+  AXOM_HOST_DEVICE
   friend bool operator!=(const Vector& lhs, const Vector& rhs)
   {
     return !(lhs == rhs);
@@ -413,6 +415,7 @@ class Vector
    * \param [in] z the z--coordinate of the vector. Default is 0.0.
    * \return v a Vector instance with the given coordinates.
    */
+  AXOM_HOST_DEVICE
   static Vector make_vector(const T& x, const T& y, const T& z = 0.0);
 
 private:
@@ -656,7 +659,8 @@ AXOM_HOST_DEVICE Point<T, NDIMS> operator-(const Point<T, NDIMS>& P,
 
 //------------------------------------------------------------------------------
 template <typename T, int NDIMS>
-inline Vector<T, NDIMS> operator/(const Vector<T, NDIMS>& vec, const T scalar)
+AXOM_HOST_DEVICE inline Vector<T, NDIMS> operator/(const Vector<T, NDIMS>& vec,
+                                                   const T scalar)
 {
   Vector<T, NDIMS> result(vec);
   result /= scalar;
@@ -683,7 +687,7 @@ AXOM_HOST_DEVICE Vector<T, NDIMS> operator-(const Point<T, NDIMS>& h,
 
 //------------------------------------------------------------------------------
 template <typename T, int NDIMS>
-inline Vector<T, NDIMS> operator-(const Vector<T, NDIMS>& vec1)
+AXOM_HOST_DEVICE inline Vector<T, NDIMS> operator-(const Vector<T, NDIMS>& vec1)
 {
   Vector<T, NDIMS> result(vec1);
   result.negate();
@@ -700,9 +704,9 @@ std::ostream& operator<<(std::ostream& os, const Vector<T, NDIMS>& vec)
 
 //------------------------------------------------------------------------------
 template <typename T, int NDIMS>
-inline Vector<T, NDIMS> Vector<T, NDIMS>::make_vector(const T& x,
-                                                      const T& y,
-                                                      const T& z)
+AXOM_HOST_DEVICE inline Vector<T, NDIMS> Vector<T, NDIMS>::make_vector(const T& x,
+                                                                       const T& y,
+                                                                       const T& z)
 {
   T tmp_array[3] = {x, y, z};
   return Vector(tmp_array, NDIMS);