diff --git a/mfem b/mfem
index dbbd425a22..6ce18b2005 160000
--- a/mfem
+++ b/mfem
@@ -1 +1 @@
-Subproject commit dbbd425a225b2ad0c3e1dc0e2b3992899e27556a
+Subproject commit 6ce18b200547a606c6121d986c1f59794dac08ec
diff --git a/src/smith/numerics/functional/boundary_integral_kernels.hpp b/src/smith/numerics/functional/boundary_integral_kernels.hpp
index 6ed2c7d362..10663f9098 100644
--- a/src/smith/numerics/functional/boundary_integral_kernels.hpp
+++ b/src/smith/numerics/functional/boundary_integral_kernels.hpp
@@ -174,7 +174,7 @@ void evaluation_kernel_impl(trial_element_type trial_elements, test_element, dou
 
   static constexpr int qpts_per_elem = num_quadrature_points(geom, Q);
 
-  [[maybe_unused]] tuple u = {
+  [[maybe_unused]] tuple u{
       reinterpret_cast<const typename decltype(type<indices>(trial_elements))::dof_type*>(inputs[indices])...};
 
   // for each element in the domain
@@ -184,9 +184,9 @@ void evaluation_kernel_impl(trial_element_type trial_elements, test_element, dou
     auto x_e = x[e];
 
     // batch-calculate values / derivatives of each trial space, at each quadrature point
-    [[maybe_unused]] tuple qf_inputs = {promote_each_to_dual_when < indices ==
-                                        differentiation_index >
-                                            (get<indices>(trial_elements).interpolate(get<indices>(u)[e], rule))...};
+    [[maybe_unused]] tuple qf_inputs{promote_each_to_dual_when < indices ==
+                                     differentiation_index >
+                                         (get<indices>(trial_elements).interpolate(get<indices>(u)[e], rule))...};
 
     // (batch) evalute the q-function at each quadrature point
     auto qf_outputs = batch_apply_qf(qf, t, x_e, J_e, get<indices>(qf_inputs)...);
diff --git a/src/smith/numerics/functional/detail/hexahedron_H1.inl b/src/smith/numerics/functional/detail/hexahedron_H1.inl
index a9122c0366..d540de8d0a 100644
--- a/src/smith/numerics/functional/detail/hexahedron_H1.inl
+++ b/src/smith/numerics/functional/detail/hexahedron_H1.inl
@@ -158,7 +158,8 @@ struct finite_element<mfem::Geometry::CUBE, H1<p, c> > {
           const auto& d10 = get<0>(get<1>(input(Q)));
           const auto& d11 = get<1>(get<1>(input(Q)));
 
-          output[Q] = {d00 * phi_j + dot(d01, dphi_j_dxi), d10 * phi_j + dot(d11, dphi_j_dxi)};
+          output[Q] = std::remove_reference_t<decltype(output[Q])>{d00 * phi_j + dot(d01, dphi_j_dxi),
+                                                                   d10 * phi_j + dot(d11, dphi_j_dxi)};
         }
       }
     }
diff --git a/src/smith/numerics/functional/detail/hexahedron_Hcurl.inl b/src/smith/numerics/functional/detail/hexahedron_Hcurl.inl
index fd32d95a7f..0dafa77601 100644
--- a/src/smith/numerics/functional/detail/hexahedron_Hcurl.inl
+++ b/src/smith/numerics/functional/detail/hexahedron_Hcurl.inl
@@ -324,7 +324,8 @@ struct finite_element<mfem::Geometry::CUBE, Hcurl<p>> {
           const auto& d10 = get<0>(get<1>(input(Q)));
           const auto& d11 = get<1>(get<1>(input(Q)));
 
-          output[Q] = {dot(d00, phi_j) + dot(d01, curl_phi_j), dot(d10, phi_j) + dot(d11, curl_phi_j)};
+          output[Q] = std::remove_reference_t<decltype(output[Q])>{dot(d00, phi_j) + dot(d01, curl_phi_j),
+                                                                   dot(d10, phi_j) + dot(d11, curl_phi_j)};
         }
       }
     }
diff --git a/src/smith/numerics/functional/detail/hexahedron_L2.inl b/src/smith/numerics/functional/detail/hexahedron_L2.inl
index 028299cd9b..7f880cb05b 100644
--- a/src/smith/numerics/functional/detail/hexahedron_L2.inl
+++ b/src/smith/numerics/functional/detail/hexahedron_L2.inl
@@ -162,7 +162,8 @@ struct finite_element<mfem::Geometry::CUBE, L2<p, c> > {
           const auto& d10 = get<0>(get<1>(input(Q)));
           const auto& d11 = get<1>(get<1>(input(Q)));
 
-          output[Q] = {d00 * phi_j + dot(d01, dphi_j_dxi), d10 * phi_j + dot(d11, dphi_j_dxi)};
+          output[Q] = std::remove_reference_t<decltype(output[Q])>{d00 * phi_j + dot(d01, dphi_j_dxi),
+                                                                   d10 * phi_j + dot(d11, dphi_j_dxi)};
         }
       }
     }
diff --git a/src/smith/numerics/functional/detail/quadrilateral_H1.inl b/src/smith/numerics/functional/detail/quadrilateral_H1.inl
index 467164b7dc..d43da00008 100644
--- a/src/smith/numerics/functional/detail/quadrilateral_H1.inl
+++ b/src/smith/numerics/functional/detail/quadrilateral_H1.inl
@@ -181,7 +181,8 @@ struct finite_element<mfem::Geometry::SQUARE, H1<p, c> > {
         const auto& d10 = get<0>(get<1>(input(Q)));
         const auto& d11 = get<1>(get<1>(input(Q)));
 
-        output[Q] = {d00 * phi_j + dot(d01, dphi_j_dxi), d10 * phi_j + dot(d11, dphi_j_dxi)};
+        output[Q] = std::remove_reference_t<decltype(output[Q])>{d00 * phi_j + dot(d01, dphi_j_dxi),
+                                                                 d10 * phi_j + dot(d11, dphi_j_dxi)};
       }
     }
 
diff --git a/src/smith/numerics/functional/detail/quadrilateral_Hcurl.inl b/src/smith/numerics/functional/detail/quadrilateral_Hcurl.inl
index 61b74c2fdb..1e1e3b20fb 100644
--- a/src/smith/numerics/functional/detail/quadrilateral_Hcurl.inl
+++ b/src/smith/numerics/functional/detail/quadrilateral_Hcurl.inl
@@ -278,7 +278,8 @@ struct finite_element<mfem::Geometry::SQUARE, Hcurl<p> > {
         const auto& d10 = get<0>(get<1>(input(Q)));
         const auto& d11 = get<1>(get<1>(input(Q)));
 
-        output[Q] = {dot(d00, phi_j) + d01 * curl_phi_j, dot(d10, phi_j) + d11 * curl_phi_j};
+        output[Q] = std::remove_reference_t<decltype(output[Q])>{dot(d00, phi_j) + d01 * curl_phi_j,
+                                                                 dot(d10, phi_j) + d11 * curl_phi_j};
       }
     }
 
diff --git a/src/smith/numerics/functional/detail/quadrilateral_L2.inl b/src/smith/numerics/functional/detail/quadrilateral_L2.inl
index c3cb575c26..c7fdd2f0c7 100644
--- a/src/smith/numerics/functional/detail/quadrilateral_L2.inl
+++ b/src/smith/numerics/functional/detail/quadrilateral_L2.inl
@@ -180,7 +180,8 @@ struct finite_element<mfem::Geometry::SQUARE, L2<p, c> > {
         const auto& d10 = get<0>(get<1>(input(Q)));
         const auto& d11 = get<1>(get<1>(input(Q)));
 
-        output[Q] = {d00 * phi_j + dot(d01, dphi_j_dxi), d10 * phi_j + dot(d11, dphi_j_dxi)};
+        output[Q] = std::remove_reference_t<decltype(output[Q])>{d00 * phi_j + dot(d01, dphi_j_dxi),
+                                                                 d10 * phi_j + dot(d11, dphi_j_dxi)};
       }
     }
 
@@ -213,7 +214,8 @@ struct finite_element<mfem::Geometry::SQUARE, L2<p, c> > {
         const auto& d10 = get<0>(get<1>(input(Q)));
         const auto& d11 = get<1>(get<1>(input(Q)));
 
-        output[Q] = {d00 * phi0_j + d01 * phi1_j, d10 * phi0_j + d11 * phi1_j};
+        output[Q] =
+            std::remove_reference_t<decltype(output[Q])>{d00 * phi0_j + d01 * phi1_j, d10 * phi0_j + d11 * phi1_j};
       }
     }
 
diff --git a/src/smith/numerics/functional/detail/segment_H1.inl b/src/smith/numerics/functional/detail/segment_H1.inl
index 43721ccece..3f05521a41 100644
--- a/src/smith/numerics/functional/detail/segment_H1.inl
+++ b/src/smith/numerics/functional/detail/segment_H1.inl
@@ -115,7 +115,8 @@ struct finite_element<mfem::Geometry::SEGMENT, H1<p, c> > {
       const auto& d10 = get<0>(get<1>(input(qx)));
       const auto& d11 = get<1>(get<1>(input(qx)));
 
-      output[qx] = {d00 * phi_j + d01 * dphi_j_dxi, d10 * phi_j + d11 * dphi_j_dxi};
+      output[qx] =
+          std::remove_reference_t<decltype(output[qx])>{d00 * phi_j + d01 * dphi_j_dxi, d10 * phi_j + d11 * dphi_j_dxi};
     }
 
     return output;
diff --git a/src/smith/numerics/functional/detail/segment_L2.inl b/src/smith/numerics/functional/detail/segment_L2.inl
index 0437b678ab..524e5b7a84 100644
--- a/src/smith/numerics/functional/detail/segment_L2.inl
+++ b/src/smith/numerics/functional/detail/segment_L2.inl
@@ -115,7 +115,8 @@ struct finite_element<mfem::Geometry::SEGMENT, L2<p, c> > {
       auto& d10 = get<0>(get<1>(input(qx)));
       auto& d11 = get<1>(get<1>(input(qx)));
 
-      output[qx] = {d00 * phi_j + d01 * dphi_j_dxi, d10 * phi_j + d11 * dphi_j_dxi};
+      output[qx] =
+          std::remove_reference_t<decltype(output[qx])>{d00 * phi_j + d01 * dphi_j_dxi, d10 * phi_j + d11 * dphi_j_dxi};
     }
 
     return output;
@@ -144,7 +145,8 @@ struct finite_element<mfem::Geometry::SEGMENT, L2<p, c> > {
       const auto& d10 = get<0>(get<1>(input(qx)));
       const auto& d11 = get<1>(get<1>(input(qx)));
 
-      output[qx] = {d00 * phi0_j + d01 * phi1_j, d10 * phi0_j + d11 * phi1_j};
+      output[qx] =
+          std::remove_reference_t<decltype(output[qx])>{d00 * phi0_j + d01 * phi1_j, d10 * phi0_j + d11 * phi1_j};
     }
 
     return output;
diff --git a/src/smith/numerics/functional/detail/tetrahedron_H1.inl b/src/smith/numerics/functional/detail/tetrahedron_H1.inl
index 05b1b3b508..9480033ec4 100644
--- a/src/smith/numerics/functional/detail/tetrahedron_H1.inl
+++ b/src/smith/numerics/functional/detail/tetrahedron_H1.inl
@@ -354,7 +354,8 @@ struct finite_element<mfem::Geometry::TETRAHEDRON, H1<p, c> > {
       const auto& d10 = get<0>(get<1>(input(i)));
       const auto& d11 = get<1>(get<1>(input(i)));
 
-      output[i] = {d00 * phi_j + dot(d01, dphi_j_dxi), d10 * phi_j + dot(d11, dphi_j_dxi)};
+      output[i] = std::remove_reference_t<decltype(output[i])>{d00 * phi_j + dot(d01, dphi_j_dxi),
+                                                               d10 * phi_j + dot(d11, dphi_j_dxi)};
     }
 
     return output;
diff --git a/src/smith/numerics/functional/detail/tetrahedron_L2.inl b/src/smith/numerics/functional/detail/tetrahedron_L2.inl
index 0f1fd4c1fe..93aa2eb554 100644
--- a/src/smith/numerics/functional/detail/tetrahedron_L2.inl
+++ b/src/smith/numerics/functional/detail/tetrahedron_L2.inl
@@ -359,7 +359,8 @@ struct finite_element<mfem::Geometry::TETRAHEDRON, L2<p, c> > {
       const auto& d10 = get<0>(get<1>(input(i)));
       const auto& d11 = get<1>(get<1>(input(i)));
 
-      output[i] = {d00 * phi_j + dot(d01, dphi_j_dxi), d10 * phi_j + dot(d11, dphi_j_dxi)};
+      output[i] = std::remove_reference_t<decltype(output[i])>{d00 * phi_j + dot(d01, dphi_j_dxi),
+                                                               d10 * phi_j + dot(d11, dphi_j_dxi)};
     }
 
     return output;
diff --git a/src/smith/numerics/functional/detail/triangle_H1.inl b/src/smith/numerics/functional/detail/triangle_H1.inl
index 26143df89c..600fa82922 100644
--- a/src/smith/numerics/functional/detail/triangle_H1.inl
+++ b/src/smith/numerics/functional/detail/triangle_H1.inl
@@ -260,7 +260,8 @@ struct finite_element<mfem::Geometry::TRIANGLE, H1<p, c> > {
       const auto& d10 = get<0>(get<1>(input(i)));
       const auto& d11 = get<1>(get<1>(input(i)));
 
-      output[i] = {d00 * phi_j + dot(d01, dphi_j_dxi), d10 * phi_j + dot(d11, dphi_j_dxi)};
+      output[i] = std::remove_reference_t<decltype(output[i])>{d00 * phi_j + dot(d01, dphi_j_dxi),
+                                                               d10 * phi_j + dot(d11, dphi_j_dxi)};
     }
 
     return output;
diff --git a/src/smith/numerics/functional/detail/triangle_L2.inl b/src/smith/numerics/functional/detail/triangle_L2.inl
index ba52595913..bad9895df7 100644
--- a/src/smith/numerics/functional/detail/triangle_L2.inl
+++ b/src/smith/numerics/functional/detail/triangle_L2.inl
@@ -270,7 +270,8 @@ struct finite_element<mfem::Geometry::TRIANGLE, L2<p, c> > {
       const auto& d10 = get<0>(get<1>(input(i)));
       const auto& d11 = get<1>(get<1>(input(i)));
 
-      output[i] = {d00 * phi_j + dot(d01, dphi_j_dxi), d10 * phi_j + dot(d11, dphi_j_dxi)};
+      output[i] = std::remove_reference_t<decltype(output[i])>{d00 * phi_j + dot(d01, dphi_j_dxi),
+                                                               d10 * phi_j + dot(d11, dphi_j_dxi)};
     }
 
     return output;
@@ -300,7 +301,8 @@ struct finite_element<mfem::Geometry::TRIANGLE, L2<p, c> > {
       const auto& d10 = get<0>(get<1>(input(i)));
       const auto& d11 = get<1>(get<1>(input(i)));
 
-      output[i] = {d00 * phi0_j + d01 * phi1_j, d10 * phi0_j + d11 * phi1_j};
+      output[i] =
+          std::remove_reference_t<decltype(output[i])>{d00 * phi0_j + d01 * phi1_j, d10 * phi0_j + d11 * phi1_j};
     }
 
     return output;
diff --git a/src/smith/numerics/functional/domain_integral_kernels.hpp b/src/smith/numerics/functional/domain_integral_kernels.hpp
index 3b2121dfaa..9b0f14853f 100644
--- a/src/smith/numerics/functional/domain_integral_kernels.hpp
+++ b/src/smith/numerics/functional/domain_integral_kernels.hpp
@@ -170,7 +170,7 @@ void evaluation_kernel_impl(trial_element_tuple trial_elements, test_element, do
 
   [[maybe_unused]] auto qpts_per_elem = num_quadrature_points(geom, Q);
 
-  [[maybe_unused]] tuple u = {
+  [[maybe_unused]] tuple u{
       reinterpret_cast<const typename decltype(type<indices>(trial_elements))::dof_type*>(inputs[indices])...};
 
   // for each element in the domain
@@ -181,9 +181,9 @@ void evaluation_kernel_impl(trial_element_tuple trial_elements, test_element, do
 
     //[[maybe_unused]] static constexpr trial_element_tuple trial_element_tuple{};
     // batch-calculate values / derivatives of each trial space, at each quadrature point
-    [[maybe_unused]] tuple qf_inputs = {promote_each_to_dual_when < indices ==
-                                        differentiation_index >
-                                            (get<indices>(trial_elements).interpolate(get<indices>(u)[e], rule))...};
+    [[maybe_unused]] tuple qf_inputs{promote_each_to_dual_when < indices ==
+                                     differentiation_index >
+                                         (get<indices>(trial_elements).interpolate(get<indices>(u)[e], rule))...};
 
     // use J_e to transform values / derivatives on the parent element
     // to the to the corresponding values / derivatives on the physical element
diff --git a/src/smith/numerics/functional/functional.hpp b/src/smith/numerics/functional/functional.hpp
index f5a40d8259..ef56f9646a 100644
--- a/src/smith/numerics/functional/functional.hpp
+++ b/src/smith/numerics/functional/functional.hpp
@@ -681,7 +681,7 @@ class Functional<test(trials...), exec> {
       // mfem::Vector arg0 = ...;
       // mfem::Vector arg1 = ...;
       // e.g. auto [value, gradient_wrt_arg1] = my_functional(arg0, differentiate_wrt(arg1));
-      return {output_T_, grad_[wrt]};
+      return smith::tuple<mfem::Vector&, Gradient&>{output_T_, grad_[wrt]};
     }
 
     if constexpr (wrt == NO_DIFFERENTIATION) {
diff --git a/src/smith/numerics/functional/functional_qoi.inl b/src/smith/numerics/functional/functional_qoi.inl
index 8946924693..4c03c992c8 100644
--- a/src/smith/numerics/functional/functional_qoi.inl
+++ b/src/smith/numerics/functional/functional_qoi.inl
@@ -405,7 +405,7 @@ class Functional<double(trials...), exec> {
       // mfem::Vector arg0 = ...;
       // mfem::Vector arg1 = ...;
       // e.g. auto [value, gradient_wrt_arg1] = my_functional(arg0, differentiate_wrt(arg1));
-      return {output_T_[0], grad_[wrt]};
+      return smith::tuple<double&, Gradient&>{output_T_[0], grad_[wrt]};
     }
 
     if constexpr (wrt == NO_DIFFERENTIATION) {
diff --git a/src/smith/numerics/functional/interior_face_integral_kernels.hpp b/src/smith/numerics/functional/interior_face_integral_kernels.hpp
index 5744793dd3..5142060379 100644
--- a/src/smith/numerics/functional/interior_face_integral_kernels.hpp
+++ b/src/smith/numerics/functional/interior_face_integral_kernels.hpp
@@ -162,7 +162,7 @@ void evaluation_kernel_impl(trial_element_type trial_elements, test_element, dou
 
   static constexpr int qpts_per_elem = num_quadrature_points(geom, Q);
 
-  [[maybe_unused]] tuple u = {
+  [[maybe_unused]] tuple u{
       reinterpret_cast<const typename decltype(type<indices>(trial_elements))::dof_type_if*>(inputs[indices])...};
 
   // for each element in the domain
@@ -172,9 +172,9 @@ void evaluation_kernel_impl(trial_element_type trial_elements, test_element, dou
     auto x_e = x[e];
 
     // batch-calculate values / derivatives of each trial space, at each quadrature point
-    [[maybe_unused]] tuple qf_inputs = {promote_each_to_dual_when < indices ==
-                                        differentiation_index >
-                                            (get<indices>(trial_elements).interpolate(get<indices>(u)[e], rule))...};
+    [[maybe_unused]] tuple qf_inputs{promote_each_to_dual_when < indices ==
+                                     differentiation_index >
+                                         (get<indices>(trial_elements).interpolate(get<indices>(u)[e], rule))...};
 
     // (batch) evalute the q-function at each quadrature point
     auto qf_outputs = batch_apply_qf(qf, t, x_e, J_e, get<indices>(qf_inputs)...);
diff --git a/src/smith/numerics/functional/tuple.hpp b/src/smith/numerics/functional/tuple.hpp
index 9331b8ef07..f6fd605826 100644
--- a/src/smith/numerics/functional/tuple.hpp
+++ b/src/smith/numerics/functional/tuple.hpp
@@ -13,227 +13,6 @@
 
 #include "mfem.hpp"
 
-#include "smith/infrastructure/accelerator.hpp"
-
-namespace mfem::future {
-
-/**
- * @brief Smith compatibility extension for the MFEM tuple implementation.
- *
- * MFEM's tuple copy currently stores up to nine elements, while Smith's
- * historical tuple API supports ten and eleven elements.
- */
-template <typename T0, typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7,
-          typename T8, typename T9>
-struct tuple<T0, T1, T2, T3, T4, T5, T6, T7, T8, T9> {
-  T0 v0;  ///< The first member of the tuple
-  T1 v1;  ///< The second member of the tuple
-  T2 v2;  ///< The third member of the tuple
-  T3 v3;  ///< The fourth member of the tuple
-  T4 v4;  ///< The fifth member of the tuple
-  T5 v5;  ///< The sixth member of the tuple
-  T6 v6;  ///< The seventh member of the tuple
-  T7 v7;  ///< The eighth member of the tuple
-  T8 v8;  ///< The ninth member of the tuple
-  T9 v9;  ///< The tenth member of the tuple
-};
-
-/**
- * @brief Smith compatibility extension for the MFEM tuple implementation.
- */
-template <typename T0, typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7,
-          typename T8, typename T9, typename T10>
-struct tuple<T0, T1, T2, T3, T4, T5, T6, T7, T8, T9, T10> {
-  T0 v0;    ///< The first member of the tuple
-  T1 v1;    ///< The second member of the tuple
-  T2 v2;    ///< The third member of the tuple
-  T3 v3;    ///< The fourth member of the tuple
-  T4 v4;    ///< The fifth member of the tuple
-  T5 v5;    ///< The sixth member of the tuple
-  T6 v6;    ///< The seventh member of the tuple
-  T7 v7;    ///< The eighth member of the tuple
-  T8 v8;    ///< The ninth member of the tuple
-  T9 v9;    ///< The tenth member of the tuple
-  T10 v10;  ///< The eleventh member of the tuple
-};
-
-namespace detail {
-
-/// @brief Return element @p i from a 10- or 11-element tuple.
-template <int i, typename Tuple>
-MFEM_HOST_DEVICE constexpr auto& tuple_get_extended(Tuple& values)
-{
-  if constexpr (i == 0) {
-    return values.v0;
-  }
-  if constexpr (i == 1) {
-    return values.v1;
-  }
-  if constexpr (i == 2) {
-    return values.v2;
-  }
-  if constexpr (i == 3) {
-    return values.v3;
-  }
-  if constexpr (i == 4) {
-    return values.v4;
-  }
-  if constexpr (i == 5) {
-    return values.v5;
-  }
-  if constexpr (i == 6) {
-    return values.v6;
-  }
-  if constexpr (i == 7) {
-    return values.v7;
-  }
-  if constexpr (i == 8) {
-    return values.v8;
-  }
-  if constexpr (i == 9) {
-    return values.v9;
-  }
-  if constexpr (i == 10) {
-    return values.v10;
-  }
-  MFEM_UNREACHABLE();
-}
-
-/// @brief Return const element @p i from a 10- or 11-element tuple.
-template <int i, typename Tuple>
-MFEM_HOST_DEVICE constexpr const auto& tuple_get_extended(const Tuple& values)
-{
-  if constexpr (i == 0) {
-    return values.v0;
-  }
-  if constexpr (i == 1) {
-    return values.v1;
-  }
-  if constexpr (i == 2) {
-    return values.v2;
-  }
-  if constexpr (i == 3) {
-    return values.v3;
-  }
-  if constexpr (i == 4) {
-    return values.v4;
-  }
-  if constexpr (i == 5) {
-    return values.v5;
-  }
-  if constexpr (i == 6) {
-    return values.v6;
-  }
-  if constexpr (i == 7) {
-    return values.v7;
-  }
-  if constexpr (i == 8) {
-    return values.v8;
-  }
-  if constexpr (i == 9) {
-    return values.v9;
-  }
-  if constexpr (i == 10) {
-    return values.v10;
-  }
-  MFEM_UNREACHABLE();
-}
-
-}  // namespace detail
-
-/// @brief Return mutable element @p i from a 10-element tuple.
-template <int i, typename T0, typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7,
-          typename T8, typename T9>
-MFEM_HOST_DEVICE constexpr auto& get(tuple<T0, T1, T2, T3, T4, T5, T6, T7, T8, T9>& values)
-{
-  static_assert(i < 10);
-  return detail::tuple_get_extended<i>(values);
-}
-
-/// @brief Return mutable element @p i from an 11-element tuple.
-template <int i, typename T0, typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7,
-          typename T8, typename T9, typename T10>
-MFEM_HOST_DEVICE constexpr auto& get(tuple<T0, T1, T2, T3, T4, T5, T6, T7, T8, T9, T10>& values)
-{
-  static_assert(i < 11);
-  return detail::tuple_get_extended<i>(values);
-}
-
-/// @brief Return const element @p i from a 10-element tuple.
-template <int i, typename T0, typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7,
-          typename T8, typename T9>
-MFEM_HOST_DEVICE constexpr const auto& get(const tuple<T0, T1, T2, T3, T4, T5, T6, T7, T8, T9>& values)
-{
-  static_assert(i < 10);
-  return detail::tuple_get_extended<i>(values);
-}
-
-/// @brief Return const element @p i from an 11-element tuple.
-template <int i, typename T0, typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7,
-          typename T8, typename T9, typename T10>
-MFEM_HOST_DEVICE constexpr const auto& get(const tuple<T0, T1, T2, T3, T4, T5, T6, T7, T8, T9, T10>& values)
-{
-  static_assert(i < 11);
-  return detail::tuple_get_extended<i>(values);
-}
-
-/**
- * @brief a function intended to be used for extracting the ith type from a tuple.
- *
- * @note type<i>(my_tuple) returns a value, whereas get<i>(my_tuple) returns a reference
- *
- * @tparam i the index of the tuple to query
- * @tparam T0 The first type stored in the tuple
- * @tparam T1 The second type stored in the tuple
- * @tparam T2 The third type stored in the tuple
- * @tparam T3 The fourth type stored in the tuple
- * @tparam T4 The fifth type stored in the tuple
- * @tparam T5 The sixth type stored in the tuple
- * @tparam T6 The seventh type stored in the tuple
- * @tparam T7 The eighth type stored in the tuple
- * @tparam T8 The ninth type stored in the tuple
- * @tparam T9 The tenth type stored in the tuple
- * @param values the tuple of values
- * @return a copy of the ith entry of the input
- */
-template <int i, typename T0, typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7,
-          typename T8, typename T9>
-MFEM_HOST_DEVICE constexpr auto type(const tuple<T0, T1, T2, T3, T4, T5, T6, T7, T8, T9>& values)
-{
-  static_assert(i < 10);
-  return detail::tuple_get_extended<i>(values);
-}
-
-/**
- * @brief a function intended to be used for extracting the ith type from a tuple.
- *
- * @note type<i>(my_tuple) returns a value, whereas get<i>(my_tuple) returns a reference
- *
- * @tparam i the index of the tuple to query
- * @tparam T0 The first type stored in the tuple
- * @tparam T1 The second type stored in the tuple
- * @tparam T2 The third type stored in the tuple
- * @tparam T3 The fourth type stored in the tuple
- * @tparam T4 The fifth type stored in the tuple
- * @tparam T5 The sixth type stored in the tuple
- * @tparam T6 The seventh type stored in the tuple
- * @tparam T7 The eighth type stored in the tuple
- * @tparam T8 The ninth type stored in the tuple
- * @tparam T9 The tenth type stored in the tuple
- * @tparam T10 The eleventh type stored in the tuple
- * @param values the tuple of values
- * @return a copy of the ith entry of the input
- */
-template <int i, typename T0, typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7,
-          typename T8, typename T9, typename T10>
-MFEM_HOST_DEVICE constexpr auto type(const tuple<T0, T1, T2, T3, T4, T5, T6, T7, T8, T9, T10>& values)
-{
-  static_assert(i < 11);
-  return detail::tuple_get_extended<i>(values);
-}
-
-}  // namespace mfem::future
-
 namespace smith {
 
 /// @brief Expose MFEM tuple in the Smith namespace.
diff --git a/src/smith/numerics/functional/tuple_tensor_dual_functions.hpp b/src/smith/numerics/functional/tuple_tensor_dual_functions.hpp
index 0b80bf3977..869e313b67 100644
--- a/src/smith/numerics/functional/tuple_tensor_dual_functions.hpp
+++ b/src/smith/numerics/functional/tuple_tensor_dual_functions.hpp
@@ -913,7 +913,7 @@ inline SMITH_HOST_DEVICE tuple<vec3, mat3> eig_symm(const mat3& A)
                 {Q[2][order[0]], Q[2][order[1]], Q[2][order[2]]}}};
   // clang-format on
 
-  return {eigvals, eigvecs};
+  return tuple<vec3, mat3>{eigvals, eigvecs};
 }
 
 /*