diff --git a/src/serac/numerics/functional/dof_numbering.hpp b/src/serac/numerics/functional/dof_numbering.hpp
index 8dce4e5f1..0c9bec041 100644
--- a/src/serac/numerics/functional/dof_numbering.hpp
+++ b/src/serac/numerics/functional/dof_numbering.hpp
@@ -277,6 +277,9 @@ struct GradientAssemblyLookupTables {
     };
   };
 
+  /// dummy default ctor to enable deferred initialization
+  GradientAssemblyLookupTables() : initialized{false} {};
+
   /**
    * @param block_test_dofs object containing information about dofs for the test space
    * @param block_trial_dofs object containing information about dofs for the trial space
@@ -284,8 +287,8 @@ struct GradientAssemblyLookupTables {
    * @brief create lookup tables describing which degrees of freedom
    * correspond to each domain/boundary element
    */
-  GradientAssemblyLookupTables(const serac::BlockElementRestriction& block_test_dofs,
-                               const serac::BlockElementRestriction& block_trial_dofs)
+  void init(const serac::BlockElementRestriction& block_test_dofs,
+            const serac::BlockElementRestriction& block_trial_dofs)
   {
     // we start by having each element and boundary element emit the (i,j) entry that it
     // touches in the global "stiffness matrix", and also keep track of some metadata about
@@ -342,6 +345,8 @@ struct GradientAssemblyLookupTables {
     }
 
     row_ptr.back() = static_cast<int>(nnz);
+
+    initialized = true;
   }
 
   /**
@@ -368,6 +373,9 @@ struct GradientAssemblyLookupTables {
    * corresponding to the (i,j) entry
    */
   std::unordered_map<Entry, uint32_t, Entry::Hasher> nz_LUT;
+
+  /// @brief specifies if the table has already been initialized or not
+  bool initialized;
 };
 
 }  // namespace serac
diff --git a/src/serac/numerics/functional/functional.hpp b/src/serac/numerics/functional/functional.hpp
index 08826501f..9cd91ed7d 100644
--- a/src/serac/numerics/functional/functional.hpp
+++ b/src/serac/numerics/functional/functional.hpp
@@ -538,7 +538,6 @@ class Functional<test(trials...), exec> {
     Gradient(Functional<test(trials...), exec>& f, uint32_t which = 0)
         : mfem::Operator(f.test_space_->GetTrueVSize(), f.trial_space_[which]->GetTrueVSize()),
           form_(f),
-          lookup_tables(f.G_test_[Domain::Type::Elements], f.G_trial_[Domain::Type::Elements][which]),
           which_argument(which),
           test_space_(f.test_space_),
           trial_space_(f.trial_space_[which]),
@@ -576,6 +575,11 @@ class Functional<test(trials...), exec> {
 
       constexpr bool col_ind_is_sorted = true;
 
+      if (!lookup_tables.initialized) {
+        lookup_tables.init(form_.G_test_[Domain::Type::Elements],
+                           form_.G_trial_[Domain::Type::Elements][which_argument]);
+      }
+
       double* values = new double[lookup_tables.nnz]{};
 
       std::map<mfem::Geometry::Type, ExecArray<double, 3, exec>> element_gradients[Domain::num_types];