diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md
index 420af7d17e..b060e18ef1 100644
--- a/RELEASE-NOTES.md
+++ b/RELEASE-NOTES.md
@@ -25,6 +25,7 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/
 - Primal: Adds `NURBSCurve::isLinear()` to check if a curve is (nearly) flat (corresponding to `BezierCurve::isLinear()`)
 - Primal: Adds `NURBSPatch::isTriviallyTrimmed()` to check if the trimming curves for a patch lie on the patch boundaries
 - Quest: Adds support for reading mfem files with variable order NURBS curves (requires mfem>4.9).
+- Quest: Adds OMP support for fast GWN methods for STL/Triangulated STEP input and linearized NURBS Curve input.
 
 ### Removed
 
diff --git a/src/axom/core/numerics/quadrature.cpp b/src/axom/core/numerics/quadrature.cpp
index 7cc0367514..709992ed95 100644
--- a/src/axom/core/numerics/quadrature.cpp
+++ b/src/axom/core/numerics/quadrature.cpp
@@ -15,6 +15,7 @@
 
 #include <cmath>
 #include <cstdint>
+#include <mutex>
 
 namespace axom
 {
@@ -158,7 +159,7 @@ void compute_gauss_legendre_data(int npts,
  * \note If this method has already been called for a given order, it will reuse the same quadrature points
  *  without needing to recompute them
  *
- * \warning The use of a static variable to store cached nodes makes this method not threadsafe.
+ * \note This implementation uses a process-wide cache protected by a mutex for thread safety.
  * 
  * \return The `QuadratureRule` object which contains axom::ArrayView<double>'s of stored nodes and weights
  */
@@ -166,8 +167,10 @@ QuadratureRule get_gauss_legendre(int npts, int allocatorID)
 {
   assert("Quadrature rules must have >= 1 point" && (npts >= 1));
 
-  // Store cached rules keyed by (npts, allocatorID). This cache is not thread-safe.
+  // Store cached rules keyed by (npts, allocatorID).
   static axom::FlatMap<std::uint64_t, GaussLegendreRuleStorage> rule_library(64);
+  static std::mutex rule_library_mutex;
+  const std::lock_guard<std::mutex> lock(rule_library_mutex);
 
   const std::uint64_t key = make_gauss_legendre_key(npts, allocatorID);
 
diff --git a/src/axom/primal/operators/detail/winding_number_2d_memoization.hpp b/src/axom/primal/operators/detail/winding_number_2d_memoization.hpp
index 87084b073a..37ba2257bf 100644
--- a/src/axom/primal/operators/detail/winding_number_2d_memoization.hpp
+++ b/src/axom/primal/operators/detail/winding_number_2d_memoization.hpp
@@ -256,6 +256,117 @@ std::ostream& operator<<(std::ostream& os, const NURBSCurveGWNCache<T>& nCurveCa
 }
 
 }  // namespace detail
+
+/*!
+ * \brief Manage an array of NURBSCurveGWNCache<double>
+ */
+class NURBSCurveCacheManager
+{
+  using NURBSCache = axom::primal::detail::NURBSCurveGWNCache<double>;
+  using NURBSCacheArray = axom::Array<NURBSCache>;
+  using NURBSCacheArrayView = axom::ArrayView<const NURBSCache>;
+
+  using CurveArrayView = axom::ArrayView<const axom::primal::NURBSCurve<double, 2>>;
+
+public:
+  NURBSCurveCacheManager() = default;
+
+  NURBSCurveCacheManager(const CurveArrayView& curves, double bbExpansionAmount = 0.0)
+  {
+    for(const auto& curve : curves)
+    {
+      m_nurbs_caches.push_back(NURBSCache(curve, bbExpansionAmount));
+    }
+  }
+
+  /// A view of the manager object.
+  struct View
+  {
+    NURBSCacheArrayView m_view;
+
+    /// Return the NURBSCacheArrayView.
+    NURBSCacheArrayView caches() const { return m_view; }
+  };
+
+  /// Return a view of this manager to pass into a device function.
+  View view() const { return View {m_nurbs_caches.view()}; }
+
+  /// Return if the underlying array is empty
+  bool empty() const { return m_nurbs_caches.empty(); }
+
+private:
+  NURBSCacheArray m_nurbs_caches;
+};
+
+template <typename ExecSpace>
+struct nurbs_cache_2d_traits
+{
+  using type = NURBSCurveCacheManager;
+};
+
+#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_OPENMP)
+/*!
+ * \brief Manage per-thread arrays of NURBSCurveGWNCache<double>
+ */
+class NURBSCurveCacheManagerOMP
+{
+  using NURBSCache = axom::primal::detail::NURBSCurveGWNCache<double>;
+  using NURBSCachePerThreadArray = axom::Array<axom::Array<NURBSCache>>;
+  using NURBSCachePerThreadArrayView = axom::ArrayView<const axom::Array<NURBSCache>>;
+  using NURBSCacheArrayView = axom::ArrayView<const NURBSCache>;
+
+  using CurveArrayView = axom::ArrayView<const axom::primal::NURBSCurve<double, 2>>;
+
+public:
+  NURBSCurveCacheManagerOMP() = default;
+
+  NURBSCurveCacheManagerOMP(const CurveArrayView& curves, double bbExpansionAmount = 0.0)
+  {
+    const int nt = omp_get_max_threads();
+    m_nurbs_caches.resize(nt);
+    auto nurbs_caches_view = m_nurbs_caches.view();
+
+    // Make the first one
+    nurbs_caches_view[0].resize(curves.size());
+    axom::for_all<axom::OMP_EXEC>(
+      curves.size(),
+      AXOM_LAMBDA(axom::IndexType i) {
+        nurbs_caches_view[0][i] = NURBSCache(curves[i], bbExpansionAmount);
+      });
+
+    // Copy the constructed cache to the other threads' copies (less work than construction)
+    axom::for_all<axom::OMP_EXEC>(
+      1,
+      nt,
+      AXOM_LAMBDA(axom::IndexType t) { nurbs_caches_view[t] = nurbs_caches_view[0]; });
+  }
+
+  /// A view of the manager object.
+  struct View
+  {
+    NURBSCachePerThreadArrayView m_views;
+
+    /// Return the NURBSCacheArrayView for the current OMP thread.
+    NURBSCacheArrayView caches() const { return m_views[omp_get_thread_num()].view(); }
+  };
+
+  /// Return a view of this manager to pass into a device function.
+  View view() const { return View {m_nurbs_caches.view()}; }
+
+  /// Return if the underlying array is empty
+  bool empty() const { return m_nurbs_caches.empty(); }
+
+private:
+  NURBSCachePerThreadArray m_nurbs_caches;
+};
+
+template <>
+struct nurbs_cache_2d_traits<axom::OMP_EXEC>
+{
+  using type = NURBSCurveCacheManagerOMP;
+};
+#endif
+
 }  // namespace primal
 }  // namespace axom
 
diff --git a/src/axom/primal/operators/detail/winding_number_3d_memoization.hpp b/src/axom/primal/operators/detail/winding_number_3d_memoization.hpp
index c4486b81c4..cb568c18ca 100644
--- a/src/axom/primal/operators/detail/winding_number_3d_memoization.hpp
+++ b/src/axom/primal/operators/detail/winding_number_3d_memoization.hpp
@@ -170,6 +170,19 @@ class NURBSPatchGWNCache
       m_averageNormal = m_alteredPatch.calculateTrimmedPatchNormal();
     }
 
+    // Cast direction is set to average normal, unless it is near zero
+    if(m_averageNormal.norm() < 1e-10)
+    {
+      // ...unless the average direction is zero
+      double theta = axom::utilities::random_real(0.0, 2 * M_PI);
+      double u = axom::utilities::random_real(-1.0, 1.0);
+      m_castDirection = Vector<T, 3> {sin(theta) * sqrt(1 - u * u), cos(theta) * sqrt(1 - u * u), u};
+    }
+    else
+    {
+      m_castDirection = m_averageNormal.unitVector();
+    }
+
     m_pboxDiag = m_alteredPatch.getParameterSpaceDiagonal();
 
     // Make a bounding box by doing (trimmed) bezier extraction,
@@ -238,6 +251,7 @@ class NURBSPatchGWNCache
   ///@{
   //! \name Accessors for precomputed data
   const Vector<T, 3>& getAverageNormal() const { return m_averageNormal; }
+  const Vector<T, 3>& getCastDirection() const { return m_castDirection; }
   const BoundingBox<T, 3>& boundingBox() const { return m_bBox; }
   const OrientedBoundingBox<T, 3>& orientedBoundingBox() const { return m_oBox; }
   //@}
@@ -272,7 +286,7 @@ class NURBSPatchGWNCache
   // Per patch data
   BoundingBox<T, 3> m_bBox;
   OrientedBoundingBox<T, 3> m_oBox;
-  Vector<T, 3> m_averageNormal;
+  Vector<T, 3> m_averageNormal, m_castDirection;
   double m_pboxDiag;
 
   // Per trimming curve data, keyed by (whichRefinementLevel, whichRefinementIndex)
@@ -280,6 +294,123 @@ class NURBSPatchGWNCache
 };
 
 }  // namespace detail
+
+/*!
+ * \brief Manage an array of NURBSPatchGWNCache<double>
+ */
+class NURBSPatchCacheManager
+{
+  using NURBSCache = axom::primal::detail::NURBSPatchGWNCache<double>;
+  using NURBSCacheArray = axom::Array<NURBSCache>;
+  using NURBSCacheArrayView = axom::ArrayView<const NURBSCache>;
+
+  using PatchArrayView = axom::ArrayView<const axom::primal::NURBSPatch<double, 3>>;
+
+public:
+  NURBSPatchCacheManager() = default;
+
+  NURBSPatchCacheManager(const PatchArrayView& patchs)
+  {
+    for(const auto& patch : patchs)
+    {
+      m_nurbs_caches.push_back(NURBSCache(patch));
+    }
+  }
+
+  /// A view of the manager object.
+  struct View
+  {
+    NURBSCacheArrayView m_view;
+
+    /// Return the NURBSCacheArrayView.
+    NURBSCacheArrayView caches() const { return m_view; }
+  };
+
+  /// Return a view of this manager to pass into a device function.
+  View view() const { return View {m_nurbs_caches.view()}; }
+
+  /// Return if the underlying array is empty
+  bool empty() const { return m_nurbs_caches.empty(); }
+
+private:
+  NURBSCacheArray m_nurbs_caches;
+};
+
+template <typename ExecSpace>
+struct nurbs_cache_3d_traits
+{
+  using type = NURBSPatchCacheManager;
+};
+
+#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_OPENMP)
+/*!
+ * \brief Manage per-thread arrays of NURBSPatchGWNCache<double>
+ */
+class NURBSPatchCacheManagerOMP
+{
+  using NURBSCache = axom::primal::detail::NURBSPatchGWNCache<double>;
+  using NURBSCachePerThreadArray = axom::Array<axom::Array<NURBSCache>>;
+  using NURBSCachePerThreadArrayView = axom::ArrayView<const axom::Array<NURBSCache>>;
+  using NURBSCacheArrayView = axom::ArrayView<const NURBSCache>;
+
+  using PatchArrayView = axom::ArrayView<const axom::primal::NURBSPatch<double, 3>>;
+
+public:
+  NURBSPatchCacheManagerOMP() = default;
+
+  NURBSPatchCacheManagerOMP(const PatchArrayView& patches)
+  {
+    const int nt = omp_get_max_threads();
+    m_nurbs_caches.resize(nt);
+    auto nurbs_caches_view = m_nurbs_caches.view();
+
+    // Make the first one
+    nurbs_caches_view[0].resize(patches.size());
+    axom::for_all<axom::OMP_EXEC>(
+      patches.size(),
+      AXOM_LAMBDA(axom::IndexType i) { nurbs_caches_view[0][i] = NURBSCache(patches[i]); });
+    SLIC_INFO("Finished the first construction");
+    // Copy the constructed cache to the other threads' copies (less work than construction)
+    axom::for_all<axom::OMP_EXEC>(
+      1,
+      nt,
+      AXOM_LAMBDA(axom::IndexType t) { nurbs_caches_view[t].resize(nurbs_caches_view[0].size()); });
+    axom::for_all<axom::OMP_EXEC>(
+      patches.size(),
+      AXOM_LAMBDA(axom::IndexType i) {
+        for(int t = 0; t < nt; t++)
+        {
+          nurbs_caches_view[t][i] = nurbs_caches_view[0][i];
+        }
+      });
+  }
+
+  /// A view of the manager object.
+  struct View
+  {
+    NURBSCachePerThreadArrayView m_views;
+
+    /// Return the NURBSCacheArrayView for the current OMP thread.
+    NURBSCacheArrayView caches() const { return m_views[omp_get_thread_num()].view(); }
+  };
+
+  /// Return a view of this manager to pass into a device function.
+  View view() const { return View {m_nurbs_caches.view()}; }
+
+  /// Return if the underlying array is empty
+  bool empty() const { return m_nurbs_caches.empty(); }
+
+private:
+  NURBSCachePerThreadArray m_nurbs_caches;
+};
+
+template <>
+struct nurbs_cache_3d_traits<axom::OMP_EXEC>
+{
+  using type = NURBSPatchCacheManagerOMP;
+};
+#endif
+
 }  // namespace primal
 }  // namespace axom
 
diff --git a/src/axom/primal/operators/winding_number.hpp b/src/axom/primal/operators/winding_number.hpp
index 3d0d97b6a0..69db83afb9 100644
--- a/src/axom/primal/operators/winding_number.hpp
+++ b/src/axom/primal/operators/winding_number.hpp
@@ -347,6 +347,23 @@ double winding_number(const Point<T, 2>& query,
   return gwn;
 }
 
+/// \brief Overload for views
+template <typename T>
+double winding_number(const Point<T, 2>& query,
+                      const axom::ArrayView<const detail::NURBSCurveGWNCache<T>>& nurbs_curve_arr,
+                      double edge_tol = 1e-8,
+                      double EPS = 1e-8)
+{
+  double gwn = 0;
+  bool dummy_isOnCurve = false;
+  for(int i = 0; i < nurbs_curve_arr.size(); ++i)
+  {
+    gwn += winding_number(query, nurbs_curve_arr[i], dummy_isOnCurve, edge_tol, EPS);
+  }
+
+  return gwn;
+}
+
 //! \brief Overload without optional return parameter
 template <typename T>
 double winding_number(const Point<T, 2>& query,
@@ -742,23 +759,9 @@ double winding_number(const Point<T, 3>& query,
                       const double disk_size = 0.01,
                       const double EPS = 1e-8)
 {
-  // Select the cast direction as an average normal of the untrimmed surface
-  auto cast_direction = nurbs.getAverageNormal();
-  if(cast_direction.norm() < 1e-10)
-  {
-    // ...unless the average direction is zero
-    double theta = axom::utilities::random_real(0.0, 2 * M_PI);
-    double u = axom::utilities::random_real(-1.0, 1.0);
-    cast_direction = Vector<T, 3> {sin(theta) * sqrt(1 - u * u), cos(theta) * sqrt(1 - u * u), u};
-  }
-  else
-  {
-    cast_direction = cast_direction.unitVector();
-  }
-
   return detail::nurbs_winding_number(query,
                                       nurbs,
-                                      cast_direction,
+                                      nurbs.getCastDirection(),
                                       edge_tol,
                                       ls_tol,
                                       quad_tol,
@@ -766,6 +769,32 @@ double winding_number(const Point<T, 3>& query,
                                       EPS);
 }
 
+/// \brief Overload for a single query and an ArrayView
+template <typename T>
+double winding_number(const Point<T, 3>& query,
+                      const axom::ArrayView<const detail::NURBSPatchGWNCache<T>>& nurbs_arr,
+                      const double edge_tol = 1e-8,
+                      const double ls_tol = 1e-8,
+                      const double quad_tol = 1e-8,
+                      const double disk_size = 0.01,
+                      const double EPS = 1e-8)
+{
+  double ret_val = 0.0;
+  for(int i = 0; i < nurbs_arr.size(); ++i)
+  {
+    ret_val += detail::nurbs_winding_number(query,
+                                            nurbs_arr[i],
+                                            nurbs_arr[i].getCastDirection(),
+                                            edge_tol,
+                                            ls_tol,
+                                            quad_tol,
+                                            disk_size,
+                                            EPS);
+  }
+
+  return ret_val;
+}
+
 /*!
  * \brief Computes the GWN for a 3D point wrt a generic 3D surface object
  *
@@ -830,26 +859,6 @@ axom::Array<double> winding_number(const axom::Array<Point<T, 3>>& query_arr,
                                    const double disk_size = 0.01,
                                    const double EPS = 1e-8)
 {
-  // Pull precomputed cast directions for each patch
-  axom::Array<Vector<T, 3>> cast_direction_arr(0, nurbs_arr.size());
-  for(int i = 0; i < nurbs_arr.size(); ++i)
-  {
-    // Select the cast direction as an average normal of the untrimmed surface
-    cast_direction_arr.emplace_back(nurbs_arr[i].getAverageNormal());
-    if(cast_direction_arr[i].norm() < 1e-10)
-    {
-      // ...unless the average direction is zero
-      double theta = axom::utilities::random_real(0.0, 2 * M_PI);
-      double u = axom::utilities::random_real(-1.0, 1.0);
-      cast_direction_arr[i] =
-        Vector<T, 3> {sin(theta) * sqrt(1 - u * u), cos(theta) * sqrt(1 - u * u), u};
-    }
-    else
-    {
-      cast_direction_arr[i] = cast_direction_arr[i].unitVector();
-    }
-  }
-
   axom::Array<double> ret_val(query_arr.size());
   for(int n = 0; n < query_arr.size(); ++n)
   {
@@ -859,7 +868,7 @@ axom::Array<double> winding_number(const axom::Array<Point<T, 3>>& query_arr,
     {
       ret_val[n] += detail::nurbs_winding_number(query_arr[n],
                                                  nurbs_arr[i],
-                                                 cast_direction_arr[i],
+                                                 nurbs_arr[i].getCastDirection(),
                                                  edge_tol,
                                                  ls_tol,
                                                  quad_tol,
diff --git a/src/axom/quest/GWNMethods.hpp b/src/axom/quest/GWNMethods.hpp
index 0f6a7b052b..0c16c0c5ec 100644
--- a/src/axom/quest/GWNMethods.hpp
+++ b/src/axom/quest/GWNMethods.hpp
@@ -128,11 +128,12 @@ void generate_gwn_query_mesh(mfem::DataCollection& dc,
 ///@{
 /// \name Query methods for 2D GWN applications
 
+template <typename ExecSpace>
 class DirectGWN2D
 {
 public:
   using CurveArrayType = axom::Array<axom::primal::NURBSCurve<double, 2>>;
-  using NURBSCacheArray = axom::Array<axom::primal::detail::NURBSCurveGWNCache<double>>;
+  using NURBSCacheManager = typename axom::primal::nurbs_cache_2d_traits<ExecSpace>::type;
 
   DirectGWN2D() = default;
 
@@ -152,11 +153,7 @@ class DirectGWN2D
       AXOM_ANNOTATE_SCOPE("preprocessing");
       if(use_memoization)
       {
-        m_nurbs_caches.reserve(input_curves.size());
-        for(const auto& curv : input_curves)
-        {
-          m_nurbs_caches.emplace_back(curv);
-        }
+        m_nurbs_cache_mgr = NURBSCacheManager(input_curves);
       }
     }
     timer.stop();
@@ -203,17 +200,16 @@ class DirectGWN2D
     {
       AXOM_ANNOTATE_SCOPE("query");
       const primal::WindingTolerances tol_copy = tol;
-      const auto input_curves_view = m_input_curves_view;
 
       // Use non-memoized form
-      if(m_nurbs_caches.empty())
+      if(m_nurbs_cache_mgr.empty())
       {
-        axom::for_all<axom::SEQ_EXEC>(num_query_points, [=, &winding, &inout](axom::IndexType nidx) {
+        axom::for_all<ExecSpace>(num_query_points, [=, &winding, &inout](axom::IndexType nidx) {
           const auto q = query_point(static_cast<int>(nidx));
           double wn {};
-          for(const auto& cache : input_curves_view)
+          for(const auto& curve : m_input_curves_view)
           {
-            wn += axom::primal::winding_number(q, cache, tol_copy.edge_tol, tol_copy.EPS);
+            wn += axom::primal::winding_number(q, curve, tol_copy.edge_tol, tol_copy.EPS);
           }
           winding[static_cast<int>(nidx)] = wn;
           inout[static_cast<int>(nidx)] = std::lround(wn);
@@ -221,14 +217,14 @@ class DirectGWN2D
       }
       else  // Use memoized form
       {
-        const auto nurbs_caches_view = m_nurbs_caches.view();
-        axom::for_all<axom::SEQ_EXEC>(num_query_points, [=, &winding, &inout](axom::IndexType nidx) {
+        const auto cache_mgr_view = m_nurbs_cache_mgr.view();
+        axom::for_all<ExecSpace>(num_query_points, [=, &winding, &inout](axom::IndexType nidx) {
           const auto q = query_point(static_cast<int>(nidx));
-          double wn {};
-          for(const auto& cache : nurbs_caches_view)
-          {
-            wn += axom::primal::winding_number(q, cache, tol_copy.edge_tol, tol_copy.EPS);
-          }
+          const auto caches_view = cache_mgr_view.caches();
+
+          const double wn =
+            axom::primal::winding_number(q, caches_view, tol_copy.edge_tol, tol_copy.EPS);
+
           winding[static_cast<int>(nidx)] = wn;
           inout[static_cast<int>(nidx)] = std::lround(wn);
         });
@@ -243,7 +239,7 @@ class DirectGWN2D
       "Querying {:L} samples in winding number field with{} memoization took {:.3Lf} seconds"
       " (@ {:.0Lf} queries per second; {:.6Lf} ms per query)",
       num_query_points,
-      m_nurbs_caches.empty() ? "out" : "",
+      m_nurbs_cache_mgr.empty() ? "out" : "",
       query_time_s,
       num_query_points / query_time_s,
       ms_per_query));
@@ -253,10 +249,10 @@ class DirectGWN2D
 
 private:
   axom::ArrayView<const axom::primal::NURBSCurve<double, 2>> m_input_curves_view;
-  NURBSCacheArray m_nurbs_caches;
+  NURBSCacheManager m_nurbs_cache_mgr;
 };
 
-template <int ORDER>
+template <typename ExecSpace, int ORDER>
 class PolylineGWN2D
 {
 public:
@@ -291,7 +287,7 @@ class PolylineGWN2D
       m_segments.resize(poly_mesh->getNumberOfCells());
       auto segments_view = m_segments.view();
 
-      axom::mint::for_all_cells<axom::SEQ_EXEC, axom::mint::xargs::coords>(
+      axom::mint::for_all_cells<ExecSpace, axom::mint::xargs::coords>(
         poly_mesh,
         AXOM_LAMBDA(axom::IndexType cellIdx,
                     const axom::numerics::Matrix<double>& coords,
@@ -316,7 +312,7 @@ class PolylineGWN2D
         auto aabbs_view = aabbs.view();
         const auto segments_view = m_segments.view();
 
-        axom::for_all<axom::SEQ_EXEC>(
+        axom::for_all<ExecSpace>(
           nlines,
           AXOM_LAMBDA(axom::IndexType i) {
             aabbs_view[i] = BoxType {segments_view[i].source(), segments_view[i].target()};
@@ -340,7 +336,7 @@ class PolylineGWN2D
         };
 
         const auto traverser = m_bvh.getTraverser();
-        m_internal_moments = traverser.reduce_tree<axom::SEQ_EXEC, GWNMoments>(compute_moments);
+        m_internal_moments = traverser.template reduce_tree<ExecSpace, GWNMoments>(compute_moments);
       }
       stage_timer.stop();
       SLIC_INFO(
@@ -396,7 +392,7 @@ class PolylineGWN2D
         const auto traverser = m_bvh.getTraverser();
         const auto internal_moments_view = m_internal_moments.view();
 
-        axom::for_all<axom::SEQ_EXEC>(num_query_points, [=, &winding, &inout](axom::IndexType index) {
+        axom::for_all<ExecSpace>(num_query_points, [=, &winding, &inout](axom::IndexType index) {
           const double wn = axom::quest::fast_approximate_winding_number(query_point(index),
                                                                          traverser,
                                                                          segments_view,
@@ -410,7 +406,7 @@ class PolylineGWN2D
       // Use direct formula
       else
       {
-        axom::for_all<axom::SEQ_EXEC>(num_query_points, [=, &winding, &inout](axom::IndexType index) {
+        axom::for_all<ExecSpace>(num_query_points, [=, &winding, &inout](axom::IndexType index) {
           const axom::primal::Point<double, 2> q = query_point(static_cast<int>(index));
           double wn {};
           for(const auto& seg : segments_view)
@@ -444,17 +440,20 @@ class PolylineGWN2D
 
   // Only needed for fast approximation method
   axom::Array<GWNMoments> m_internal_moments;
-  axom::spin::BVH<2, axom::SEQ_EXEC> m_bvh;
+  axom::spin::BVH<2, ExecSpace> m_bvh;
 };
 ///@}
 
 ///@{
 /// \name Query methods for 3D GWN applications
+
+template <typename ExecSpace>
 class DirectGWN3D
 {
 public:
   using PatchArrayType = axom::Array<axom::primal::NURBSPatch<double, 3>>;
   using NURBSCacheArray = axom::Array<axom::primal::detail::NURBSPatchGWNCache<double>>;
+  using NURBSCacheManager = typename axom::primal::nurbs_cache_3d_traits<ExecSpace>::type;
 
   DirectGWN3D() = default;
 
@@ -474,11 +473,7 @@ class DirectGWN3D
       AXOM_ANNOTATE_SCOPE("preprocessing");
       if(use_memoization)
       {
-        m_nurbs_caches.reserve(input_patches.size());
-        for(const auto& patch : input_patches)
-        {
-          m_nurbs_caches.emplace_back(patch);
-        }
+        m_nurbs_cache_mgr = NURBSCacheManager(input_patches);
       }
     }
     timer.stop();
@@ -532,9 +527,9 @@ class DirectGWN3D
       const auto input_patches_view = m_input_patches_view;
 
       // Use non-memoized form
-      if(m_nurbs_caches.empty())
+      if(m_nurbs_cache_mgr.empty())
       {
-        axom::for_all<axom::SEQ_EXEC>(num_query_points, [=, &winding, &inout](axom::IndexType nidx) {
+        axom::for_all<ExecSpace>(num_query_points, [=, &winding, &inout](axom::IndexType nidx) {
           const auto q = query_point(static_cast<int>(nidx));
           double wn {};
           for(const auto& patch : input_patches_view)
@@ -553,20 +548,19 @@ class DirectGWN3D
       }
       else  // Use memoized form
       {
-        const auto nurbs_patches_view = m_nurbs_caches.view();
-        axom::for_all<axom::SEQ_EXEC>(num_query_points, [=, &winding, &inout](axom::IndexType nidx) {
+        const auto cache_mgr_view = m_nurbs_cache_mgr.view();
+        axom::for_all<ExecSpace>(num_query_points, [=, &winding, &inout](axom::IndexType nidx) {
           const auto q = query_point(static_cast<int>(nidx));
-          double wn {};
-          for(const auto& cache : nurbs_patches_view)
-          {
-            wn += axom::primal::winding_number(q,
-                                               cache,
-                                               tol_copy.edge_tol,
-                                               tol_copy.ls_tol,
-                                               tol_copy.quad_tol,
-                                               tol_copy.disk_size,
-                                               tol_copy.EPS);
-          }
+          const auto caches_view = cache_mgr_view.caches();
+
+          const double wn = axom::primal::winding_number(q,
+                                                         caches_view,
+                                                         tol_copy.edge_tol,
+                                                         tol_copy.ls_tol,
+                                                         tol_copy.quad_tol,
+                                                         tol_copy.disk_size,
+                                                         tol_copy.EPS);
+
           winding[static_cast<int>(nidx)] = wn;
           inout[static_cast<int>(nidx)] = std::lround(wn);
         });
@@ -581,7 +575,7 @@ class DirectGWN3D
       "Querying {:L} samples in winding number field with{} memoization took {:.3Lf} seconds"
       " (@ {:.0Lf} queries per second; {:.6Lf} ms per query)",
       num_query_points,
-      m_nurbs_caches.empty() ? "out" : "",
+      m_nurbs_cache_mgr.empty() ? "out" : "",
       query_time_s,
       num_query_points / query_time_s,
       ms_per_query));
@@ -591,10 +585,10 @@ class DirectGWN3D
 
 private:
   axom::ArrayView<const axom::primal::NURBSPatch<double, 3>> m_input_patches_view;
-  NURBSCacheArray m_nurbs_caches;
+  NURBSCacheManager m_nurbs_cache_mgr;
 };
 
-template <int ORDER>
+template <typename ExecSpace, int ORDER>
 class TriangleGWN3D
 {
 public:
@@ -628,7 +622,7 @@ class TriangleGWN3D
       // Iterate over mesh nodes and get a bounding box for the shape
       BoxType shape_bbox;
       BoxType* shape_bbox_ptr = &shape_bbox;
-      axom::mint::for_all_nodes<axom::SEQ_EXEC, axom::mint::xargs::xyz>(
+      axom::mint::for_all_nodes<ExecSpace, axom::mint::xargs::xyz>(
         tri_mesh,
         AXOM_LAMBDA(axom::IndexType, double x, double y, double z) {
           shape_bbox_ptr->addPoint(Point3D {x, y, z});
@@ -645,7 +639,7 @@ class TriangleGWN3D
 
       m_triangles.resize(ntris);
       auto triangles_view = m_triangles.view();
-      axom::mint::for_all_cells<axom::SEQ_EXEC, axom::mint::xargs::coords>(
+      axom::mint::for_all_cells<ExecSpace, axom::mint::xargs::coords>(
         tri_mesh,
         AXOM_LAMBDA(axom::IndexType cellIdx,
                     const axom::numerics::Matrix<double>& coords,
@@ -677,7 +671,7 @@ class TriangleGWN3D
         auto aabbs_view = aabbs.view();
         const auto triangles_view = m_triangles.view();
 
-        axom::for_all<axom::SEQ_EXEC>(
+        axom::for_all<ExecSpace>(
           ntris,
           AXOM_LAMBDA(axom::IndexType i) {
             aabbs_view[i] =
@@ -702,7 +696,7 @@ class TriangleGWN3D
         };
 
         const auto traverser = m_bvh.getTraverser();
-        m_internal_moments = traverser.reduce_tree<axom::SEQ_EXEC, GWNMoments>(compute_moments);
+        m_internal_moments = traverser.template reduce_tree<ExecSpace, GWNMoments>(compute_moments);
       }
       stage_timer.stop();
       SLIC_INFO(
@@ -766,7 +760,7 @@ class TriangleGWN3D
         const auto traverser = m_bvh.getTraverser();
         const auto internal_moments_view = m_internal_moments.view();
 
-        axom::for_all<axom::SEQ_EXEC>(num_query_points, [=, &winding, &inout](axom::IndexType index) {
+        axom::for_all<ExecSpace>(num_query_points, [=, &winding, &inout](axom::IndexType index) {
           const double wn = axom::quest::fast_approximate_winding_number(scaled_query_point(index),
                                                                          traverser,
                                                                          triangles_view,
@@ -780,7 +774,7 @@ class TriangleGWN3D
       // Use direct formula
       else
       {
-        axom::for_all<axom::SEQ_EXEC>(num_query_points, [=, &winding, &inout](axom::IndexType index) {
+        axom::for_all<ExecSpace>(num_query_points, [=, &winding, &inout](axom::IndexType index) {
           const auto q = scaled_query_point(static_cast<int>(index));
           double wn {};
           for(const auto& tri : triangles_view)
@@ -814,7 +808,7 @@ class TriangleGWN3D
 
   // Only needed for fast approximation method
   axom::Array<GWNMoments> m_internal_moments;
-  axom::spin::BVH<3, axom::SEQ_EXEC> m_bvh;
+  axom::spin::BVH<3, ExecSpace> m_bvh;
 
   // Parameters for normalization
   axom::primal::Point<double, 3> m_shape_center;
@@ -837,23 +831,23 @@ enum class GWNInputType
 template <typename GWNQueryType>
 struct gwn_input_traits;
 
-template <int ORDER>
-struct gwn_input_traits<axom::quest::PolylineGWN2D<ORDER>>
+template <typename ExecSpace, int ORDER>
+struct gwn_input_traits<axom::quest::PolylineGWN2D<ExecSpace, ORDER>>
   : std::integral_constant<GWNInputType, GWNInputType::Polyline>
 { };
 
-template <>
-struct gwn_input_traits<axom::quest::DirectGWN2D>
+template <typename ExecSpace>
+struct gwn_input_traits<axom::quest::DirectGWN2D<ExecSpace>>
   : std::integral_constant<GWNInputType, GWNInputType::Curve>
 { };
 
-template <int ORDER>
-struct gwn_input_traits<axom::quest::TriangleGWN3D<ORDER>>
+template <typename ExecSpace, int ORDER>
+struct gwn_input_traits<axom::quest::TriangleGWN3D<ExecSpace, ORDER>>
   : std::integral_constant<GWNInputType, GWNInputType::Triangulation>
 { };
 
-template <>
-struct gwn_input_traits<axom::quest::DirectGWN3D>
+template <typename ExecSpace>
+struct gwn_input_traits<axom::quest::DirectGWN3D<ExecSpace>>
   : std::integral_constant<GWNInputType, GWNInputType::Surface>
 { };
 
diff --git a/src/axom/quest/examples/CMakeLists.txt b/src/axom/quest/examples/CMakeLists.txt
index b9574c2cee..b6431e547e 100644
--- a/src/axom/quest/examples/CMakeLists.txt
+++ b/src/axom/quest/examples/CMakeLists.txt
@@ -719,6 +719,7 @@ if(MFEM_FOUND AND OPENCASCADE_FOUND)
                         --output-prefix wn3d_${_model}_res${_res_i}
                         --no-vis
                         --stats
+                        --policy seq
                         query_mesh --res ${_res_i} ${_res_j} ${_res_k} --order 1)
             set_tests_properties(${_testname} PROPERTIES
                 PASS_REGULAR_EXPRESSION "${_pass_regex}")
diff --git a/src/axom/quest/examples/quest_winding_number_2d.cpp b/src/axom/quest/examples/quest_winding_number_2d.cpp
index 89d227f328..5e29371646 100644
--- a/src/axom/quest/examples/quest_winding_number_2d.cpp
+++ b/src/axom/quest/examples/quest_winding_number_2d.cpp
@@ -5,7 +5,7 @@
 // SPDX-License-Identifier: (BSD-3-Clause)
 
 /*!
- * \file quest_winding_number2d.cpp
+ * \file quest_winding_number_2d.cpp
  * \brief Example that computes the winding number of a grid of points
  * against a collection of 2D parametric rational curves.
  * Supports MFEM meshes in the cubic positive Bernstein basis or the (rational)
@@ -42,16 +42,18 @@ using BoundingBox2D = primal::BoundingBox<double, 2>;
 
 using NURBSCurve2D = primal::NURBSCurve<double, 2>;
 
+using RuntimePolicy = axom::runtime_policy::Policy;
+
 namespace
 {
 /**
- * This helper class takes in an mfem mesh (potentially with variable order curves from mfem>4.9)
- * and writes out a version that is compatible with mfem@4.9
- * In particular, this allows us to visualize it with current versions of VisIt which do not yet support this feature.
- *
- * We can remove this class once downstream applications (such as VisIt) are updated to a version of mfem
- * that supports the NURBS patches format.
- */
+   * This helper class takes in an mfem mesh (potentially with variable order curves from mfem>4.9)
+   * and writes out a version that is compatible with mfem@4.9
+   * In particular, this allows us to visualize it with current versions of VisIt which do not yet support this feature.
+   *
+   * We can remove this class once downstream applications (such as VisIt) are updated to a version of mfem
+   * that supports the NURBS patches format.
+   */
 class MFEM49ElevatedNURBSMeshWriter
 {
 public:
@@ -270,6 +272,8 @@ class Input
   bool stats {false};
   std::string elevatedMeshFile;
 
+  axom::runtime_policy::Policy policy = RuntimePolicy::seq;
+
   const std::array<std::string, 2> valid_algorithms {"direct", "fast-approximation"};
   std::string algorithm {valid_algorithms[1]};  // fast-approximation
 
@@ -333,6 +337,16 @@ class Input
       ->capture_default_str()
       ->check(axom::utilities::ValidCaliperMode);
 #endif
+    std::stringstream pol_sstr;
+    pol_sstr << "Set runtime policy method.";
+    pol_sstr << "\nSet to 'seq' or 0 to use the RAJA sequential policy.";
+#ifdef AXOM_RUNTIME_POLICY_USE_OPENMP
+    pol_sstr << "\nSet to 'omp' or 1 to use the RAJA OpenMP policy.";
+#endif
+
+    app.add_option("-p, --policy", policy, pol_sstr.str())
+      ->capture_default_str()
+      ->transform(axom::CLI::CheckedTransformer(axom::runtime_policy::s_nameToPolicy));
 
     // Options for triangulation of the input STEP file
     auto* linearize_curves_subcommand =
@@ -396,30 +410,55 @@ class Input
   }
 };
 
-using GWNQueryType = std::variant<axom::quest::DirectGWN2D,
-                                  axom::quest::PolylineGWN2D<0>,
-                                  axom::quest::PolylineGWN2D<1>,
-                                  axom::quest::PolylineGWN2D<2>>;
+using GWNQueryType = std::variant<axom::quest::DirectGWN2D<axom::SEQ_EXEC>,
+                                  axom::quest::PolylineGWN2D<axom::SEQ_EXEC, 0>,
+                                  axom::quest::PolylineGWN2D<axom::SEQ_EXEC, 1>,
+                                  axom::quest::PolylineGWN2D<axom::SEQ_EXEC, 2>
+#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_OPENMP)
+                                  ,
+                                  axom::quest::DirectGWN2D<axom::OMP_EXEC>,
+                                  axom::quest::PolylineGWN2D<axom::OMP_EXEC, 0>,
+                                  axom::quest::PolylineGWN2D<axom::OMP_EXEC, 1>,
+                                  axom::quest::PolylineGWN2D<axom::OMP_EXEC, 2>
+#endif
+                                  >;
 
-GWNQueryType make_gwn_query(bool linearize_curves, int approximation_order)
+template <typename ExecSpace>
+GWNQueryType pick_gwn_method(bool linearize_curves, int approximation_order)
 {
   if(linearize_curves)
   {
     if(approximation_order == 0)
     {
-      return axom::quest::PolylineGWN2D<0> {};
+      return axom::quest::PolylineGWN2D<ExecSpace, 0> {};
     }
     else if(approximation_order == 1)
     {
-      return axom::quest::PolylineGWN2D<1> {};
+      return axom::quest::PolylineGWN2D<ExecSpace, 1> {};
     }
-    else
+    else  // approximation_order == 2
     {
-      return axom::quest::PolylineGWN2D<2> {};
+      return axom::quest::PolylineGWN2D<ExecSpace, 2> {};
     }
   }
 
-  return axom::quest::DirectGWN2D {};
+  return axom::quest::DirectGWN2D<ExecSpace> {};
+}
+
+GWNQueryType make_gwn_query(axom::runtime_policy::Policy policy,
+                            bool linearize_curves,
+                            int approximation_order)
+{
+#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_OPENMP)
+  if(policy == RuntimePolicy::omp)
+  {
+    SLIC_INFO(axom::fmt::format("Using policy omp with {} threads", omp_get_max_threads()));
+    return pick_gwn_method<axom::OMP_EXEC>(linearize_curves, approximation_order);
+  }
+#endif
+
+  SLIC_INFO("Using policy seq");
+  return pick_gwn_method<axom::SEQ_EXEC>(linearize_curves, approximation_order);
 }
 
 int main(int argc, char** argv)
@@ -442,7 +481,7 @@ int main(int argc, char** argv)
   }
 
   axom::utilities::raii::AnnotationsWrapper annotation_raii_wrapper(input.annotationMode);
-  AXOM_ANNOTATE_SCOPE("winding number example");
+  AXOM_ANNOTATE_SCOPE("2D winding number example");
 
   if(!input.elevatedMeshFile.empty())
   {
@@ -515,7 +554,8 @@ int main(int argc, char** argv)
   mfem::DataCollection dc("winding_query");
   {
     // Create the desired winding number query instance
-    auto wn_query = make_gwn_query(app.got_subcommand("linearize_curves"), input.approximation_order);
+    auto wn_query =
+      make_gwn_query(input.policy, app.got_subcommand("linearize_curves"), input.approximation_order);
 
     // Generate the query grid and fields
     quest::generate_gwn_query_mesh(dc,
diff --git a/src/axom/quest/examples/quest_winding_number_3d.cpp b/src/axom/quest/examples/quest_winding_number_3d.cpp
index 5dda1c53eb..e043e0c2a5 100644
--- a/src/axom/quest/examples/quest_winding_number_3d.cpp
+++ b/src/axom/quest/examples/quest_winding_number_3d.cpp
@@ -40,6 +40,8 @@ using BoundingBox3D = primal::BoundingBox<double, 3>;
 using NURBSPatch3D = quest::STEPReader::NURBSPatch;
 using Triangle3D = primal::Triangle<double, 3>;
 
+using RuntimePolicy = axom::runtime_policy::Policy;
+
 //------------------------------------------------------------------------------
 // CLI input
 //------------------------------------------------------------------------------
@@ -57,6 +59,8 @@ class Input
   bool validate {false};
   bool stats {false};
 
+  axom::runtime_policy::Policy policy = RuntimePolicy::seq;
+
   const std::array<std::string, 2> valid_algorithms {"direct", "fast-approximation"};
   std::string algorithm {valid_algorithms[1]};  // fast-approximation
 
@@ -160,6 +164,16 @@ class Input
       ->capture_default_str()
       ->check(axom::utilities::ValidCaliperMode);
 #endif
+    std::stringstream pol_sstr;
+    pol_sstr << "Set runtime policy method.";
+    pol_sstr << "\nSet to 'seq' or 0 to use the RAJA sequential policy.";
+#ifdef AXOM_RUNTIME_POLICY_USE_OPENMP
+    pol_sstr << "\nSet to 'omp' or 1 to use the RAJA OpenMP policy.";
+#endif
+
+    app.add_option("-p, --policy", policy, pol_sstr.str())
+      ->capture_default_str()
+      ->transform(axom::CLI::CheckedTransformer(axom::runtime_policy::s_nameToPolicy));
 
     auto* query_mesh_subcommand =
       app.add_subcommand("query_mesh")->description("Options for setting up a query mesh")->fallthrough();
@@ -231,30 +245,55 @@ class Input
   }
 };
 
-using GWNQueryType = std::variant<axom::quest::DirectGWN3D,
-                                  axom::quest::TriangleGWN3D<0>,
-                                  axom::quest::TriangleGWN3D<1>,
-                                  axom::quest::TriangleGWN3D<2>>;
+using GWNQueryType = std::variant<axom::quest::DirectGWN3D<axom::SEQ_EXEC>,
+                                  axom::quest::TriangleGWN3D<axom::SEQ_EXEC, 0>,
+                                  axom::quest::TriangleGWN3D<axom::SEQ_EXEC, 1>,
+                                  axom::quest::TriangleGWN3D<axom::SEQ_EXEC, 2>
+#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_OPENMP)
+                                  ,
+                                  axom::quest::DirectGWN3D<axom::OMP_EXEC>,
+                                  axom::quest::TriangleGWN3D<axom::OMP_EXEC, 0>,
+                                  axom::quest::TriangleGWN3D<axom::OMP_EXEC, 1>,
+                                  axom::quest::TriangleGWN3D<axom::OMP_EXEC, 2>
+#endif
+                                  >;
 
-GWNQueryType make_gwn_query(Input input)
+template <typename ExecSpace>
+GWNQueryType pick_gwn_method(bool triangulate, int approximation_order)
 {
-  if(input.triangulate)
+  if(triangulate)
   {
-    if(input.approximation_order == 0)
+    if(approximation_order == 0)
     {
-      return axom::quest::TriangleGWN3D<0> {};
+      return axom::quest::TriangleGWN3D<ExecSpace, 0> {};
     }
-    else if(input.approximation_order == 1)
+    else if(approximation_order == 1)
     {
-      return axom::quest::TriangleGWN3D<1> {};
+      return axom::quest::TriangleGWN3D<ExecSpace, 1> {};
     }
-    else
+    else  // approximation_order == 2
     {
-      return axom::quest::TriangleGWN3D<2> {};
+      return axom::quest::TriangleGWN3D<ExecSpace, 2> {};
     }
   }
 
-  return axom::quest::DirectGWN3D {};
+  return axom::quest::DirectGWN3D<ExecSpace> {};
+}
+
+GWNQueryType make_gwn_query(axom::runtime_policy::Policy policy,
+                            bool triangulate,
+                            int approximation_order)
+{
+#if defined(AXOM_USE_RAJA) && defined(AXOM_USE_OPENMP)
+  if(policy == RuntimePolicy::omp)
+  {
+    SLIC_INFO(axom::fmt::format("Using policy omp with {} threads", omp_get_max_threads()));
+    return pick_gwn_method<axom::OMP_EXEC>(triangulate, approximation_order);
+  }
+#endif
+
+  SLIC_INFO("Using policy seq");
+  return pick_gwn_method<axom::SEQ_EXEC>(triangulate, approximation_order);
 }
 
 int main(int argc, char** argv)
@@ -347,7 +386,7 @@ int main(int argc, char** argv)
     SLIC_INFO(axom::fmt::format(
       axom::utilities::locale(),
       "Loaded {} trimmed NURBS patches (with {} trimming curves) in {:.3Lf} seconds",
-      patches.size(),
+      step_reader.numPatches(),
       num_trimming_curves,
       read_timer.elapsed()));
 
@@ -398,7 +437,7 @@ int main(int argc, char** argv)
   mfem::DataCollection dc("winding_query");
   {
     // Create the desired winding number query instance
-    auto wn_query = make_gwn_query(input);
+    auto wn_query = make_gwn_query(input.policy, input.triangulate, input.approximation_order);
 
     // Generate the query grid and fields
     quest::generate_gwn_query_mesh(dc,
diff --git a/src/axom/quest/tests/quest_gwn_methods.cpp b/src/axom/quest/tests/quest_gwn_methods.cpp
index 886594caaf..3238eb825c 100644
--- a/src/axom/quest/tests/quest_gwn_methods.cpp
+++ b/src/axom/quest/tests/quest_gwn_methods.cpp
@@ -154,7 +154,8 @@ TEST(quest_gwn_methods, gwn_moment_data_triangle)
 }
 
 //------------------------------------------------------------------------------
-TEST(quest_gwn_methods, mfem_mesh_linearization)
+template <typename ExecSpace>
+void check_mfem_mesh_linearization()
 {
   using NURBSCurve2D = axom::primal::NURBSCurve<double, 2>;
   const std::string fileName = pjoin(AXOM_DATA_DIR, "contours", "svg", "mfem_logo_simp.mesh");
@@ -212,19 +213,19 @@ TEST(quest_gwn_methods, mfem_mesh_linearization)
 
   // Direct
   SLIC_INFO("Testing Direct Evaluation");
-  axom::quest::DirectGWN2D gwn_direct {};
+  axom::quest::DirectGWN2D<ExecSpace> gwn_direct {};
   gwn_direct.preprocess(curves);
   gwn_direct.query(dc[0], tol);
 
   // Linearized
   SLIC_INFO("Testing Direct Evaluation of Triangulation");
-  axom::quest::PolylineGWN2D<0> gwn_polyline {};
+  axom::quest::PolylineGWN2D<ExecSpace, 0> gwn_polyline {};
   gwn_polyline.preprocess(&poly_mesh, useDirectPolyline);
   gwn_polyline.query(dc[1], tol);
 
   // Linearized, fast approximation
   SLIC_INFO("Testing Fast-Approximate Evaluation of Triangulation");
-  axom::quest::PolylineGWN2D<0> gwn_polyline_fast {};
+  axom::quest::PolylineGWN2D<ExecSpace, 0> gwn_polyline_fast {};
   gwn_polyline_fast.preprocess(&poly_mesh, !useDirectPolyline);
   gwn_polyline_fast.query(dc[2], tol);
 
@@ -245,7 +246,8 @@ TEST(quest_gwn_methods, mfem_mesh_linearization)
 
 #ifdef AXOM_USE_OPENCASCADE
 //------------------------------------------------------------------------------
-TEST(quest_gwn_methods, step_file_triangulation)
+template <typename ExecSpace>
+void check_step_file_triangulation()
 {
   const std::string fileName = pjoin(AXOM_DATA_DIR, "quest", "step", "nut.step");
 
@@ -298,19 +300,19 @@ TEST(quest_gwn_methods, step_file_triangulation)
 
   // Direct
   SLIC_INFO("Testing Direct Evaluation");
-  axom::quest::DirectGWN3D gwn_direct {};
+  axom::quest::DirectGWN3D<ExecSpace> gwn_direct {};
   gwn_direct.preprocess(patches);
   gwn_direct.query(dc[0], tol);
 
   // Triangulated
   SLIC_INFO("Testing Direct Evaluation of Polyline");
-  axom::quest::TriangleGWN3D<0> gwn_tri {};
+  axom::quest::TriangleGWN3D<ExecSpace, 0> gwn_tri {};
   gwn_tri.preprocess(&tri_mesh, useDirectTriangle);
   gwn_tri.query(dc[1], tol);
 
   // Triangulated, fast approximation
   SLIC_INFO("Testing Fast-Approximate Evaluation of Polyline");
-  axom::quest::TriangleGWN3D<0> gwn_tri_fast {};
+  axom::quest::TriangleGWN3D<ExecSpace, 0> gwn_tri_fast {};
   gwn_tri_fast.preprocess(&tri_mesh, !useDirectTriangle);
   gwn_tri_fast.query(dc[2], tol);
 
@@ -330,6 +332,33 @@ TEST(quest_gwn_methods, step_file_triangulation)
 }
 #endif
 
+//------------------------------------------------------------------------------
+TEST(quest_gwn_methods, mfem_mesh_linearization)
+{
+  check_mfem_mesh_linearization<axom::SEQ_EXEC>();
+}
+
+#if defined AXOM_USE_OPENMP && defined(AXOM_USE_RAJA)
+TEST(quest_gwn_methods, mfem_mesh_linearization_omp)
+{
+  check_mfem_mesh_linearization<axom::OMP_EXEC>();
+}
+#endif
+
+#ifdef AXOM_USE_OPENCASCADE
+TEST(quest_gwn_methods, step_file_triangulation)
+{
+  check_step_file_triangulation<axom::SEQ_EXEC>();
+}
+#endif
+
+#if defined(AXOM_USE_OPENCASCADE) && defined(AXOM_USE_OPENMP) && defined(AXOM_USE_RAJA)
+TEST(quest_gwn_methods, step_file_triangulation_omp)
+{
+  check_step_file_triangulation<axom::OMP_EXEC>();
+}
+#endif
+
 //------------------------------------------------------------------------------
 int main(int argc, char *argv[])
 {