diff --git a/include/model_variables.h b/include/model_variables.h
index 20d9f7e5a..6f3da0c2e 100644
--- a/include/model_variables.h
+++ b/include/model_variables.h
@@ -34,7 +34,6 @@ class modelResidual
 struct variable_info
 {
   bool                                     is_scalar;
-  unsigned int                             variable_index;
   unsigned int                             global_var_index;
   dealii::EvaluationFlags::EvaluationFlags evaluation_flags;
   dealii::EvaluationFlags::EvaluationFlags residual_flags;
diff --git a/include/variableContainer.h b/include/variableContainer.h
index bfc1518f3..e5a0fb9be 100644
--- a/include/variableContainer.h
+++ b/include/variableContainer.h
@@ -3,11 +3,14 @@
 #ifndef VARIBLECONTAINER_H
 #define VARIBLECONTAINER_H
 
+#include <deal.II/base/exceptions.h>
 #include <deal.II/lac/vector.h>
 #include <deal.II/matrix_free/evaluation_flags.h>
 #include <deal.II/matrix_free/fe_evaluation.h>
 #include <deal.II/matrix_free/matrix_free.h>
 
+#include <boost/unordered_map.hpp>
+
 #include "userInputParameters.h"
 
 template <int dim, int degree, typename T>
@@ -19,15 +22,16 @@ class variableContainer
 
   // Standard contructor, used for most situations
   variableContainer(const dealii::MatrixFree<dim, double> &data,
-                    std::vector<variable_info>             _varInfoList,
-                    std::vector<variable_info>             _varChangeInfoList);
+                    const std::vector<variable_info>      &_varInfoList,
+                    const std::vector<variable_info>      &_varChangeInfoList);
+
   variableContainer(const dealii::MatrixFree<dim, double> &data,
-                    std::vector<variable_info>             _varInfoList);
+                    const std::vector<variable_info>      &_varInfoList);
   // Nonstandard constructor, used when only one index of "data" should be used,
   // use with care!
   variableContainer(const dealii::MatrixFree<dim, double> &data,
-                    std::vector<variable_info>             _varInfoList,
-                    unsigned int                           fixed_index);
+                    const std::vector<variable_info>      &_varInfoList,
+                    const unsigned int                    &fixed_index);
 
   // Methods to get the value/grad/hess in the residual method (this is how the
   // user gets these values in equations.h)
@@ -90,11 +94,6 @@ class variableContainer
   reinit_and_eval_change_in_solution(const vectorType &src,
                                      unsigned int      cell,
                                      unsigned int      var_being_solved);
-  void
-  reinit_and_eval_LHS(const vectorType               &src,
-                      const std::vector<vectorType *> solutionSet,
-                      unsigned int                    cell,
-                      unsigned int                    var_being_solved);
 
   // Only initialize the FEEvaluation object for each variable (used for
   // post-processing)
@@ -111,29 +110,34 @@ class variableContainer
   // The quadrature point index, a method to get the number of quadrature points
   // per cell, and a method to get the xyz coordinates for the quadrature point
   unsigned int q_point;
+
   unsigned int
-  get_num_q_points();
+  get_num_q_points() const;
+
   dealii::Point<dim, T>
-  get_q_point_location();
+  get_q_point_location() const;
 
 private:
-  // The number of variables
-  unsigned int num_var;
-
   // Vectors of the actual FEEvaluation objects for each active variable, split
   // into scalar variables and vector variables for type reasons
-  std::vector<dealii::FEEvaluation<dim, degree, degree + 1, 1, double>>   scalar_vars;
-  std::vector<dealii::FEEvaluation<dim, degree, degree + 1, dim, double>> vector_vars;
+  using scalar_FEEval = dealii::FEEvaluation<dim, degree, degree + 1, 1, double>;
+  using vector_FEEval = dealii::FEEvaluation<dim, degree, degree + 1, dim, double>;
+
+  boost::unordered_map<unsigned int, std::unique_ptr<scalar_FEEval>> scalar_vars_map;
+  boost::unordered_map<unsigned int, std::unique_ptr<vector_FEEval>> vector_vars_map;
 
-  std::vector<dealii::FEEvaluation<dim, degree, degree + 1, 1, double>>
-    scalar_change_in_vars;
-  std::vector<dealii::FEEvaluation<dim, degree, degree + 1, dim, double>>
-    vector_change_in_vars;
+  boost::unordered_map<unsigned int, std::unique_ptr<scalar_FEEval>>
+    scalar_change_in_vars_map;
+  boost::unordered_map<unsigned int, std::unique_ptr<vector_FEEval>>
+    vector_change_in_vars_map;
 
   // Object containing some information about each variable (indices, whether
   // the val/grad/hess is needed, etc)
   std::vector<variable_info> varInfoList;
   std::vector<variable_info> varChangeInfoList;
+
+  // The number of variables
+  unsigned int num_var;
 };
 
 #endif
diff --git a/src/matrixfree/computeLHS.cc b/src/matrixfree/computeLHS.cc
index 806866ddf..2606b1ade 100644
--- a/src/matrixfree/computeLHS.cc
+++ b/src/matrixfree/computeLHS.cc
@@ -57,7 +57,6 @@ MatrixFreePDE<dim, degree>::getLHS(
   for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
     {
       // Initialize, read DOFs, and set evaulation flags for each variable
-      // variable_list.reinit_and_eval_LHS(src,solutionSet,cell,currentFieldIndex);
       variable_list.reinit_and_eval(solutionSet, cell);
       variable_list.reinit_and_eval_change_in_solution(src, cell, currentFieldIndex);
 
diff --git a/src/userInputParameters/loadVariableAttributes.cc b/src/userInputParameters/loadVariableAttributes.cc
index 4edb12c7e..942936973 100644
--- a/src/userInputParameters/loadVariableAttributes.cc
+++ b/src/userInputParameters/loadVariableAttributes.cc
@@ -50,8 +50,6 @@ userInputParameters<dim>::loadVariableAttributes(
         }
     }
   varInfoListExplicitRHS.reserve(num_var_explicit_RHS);
-  unsigned int scalar_var_index = 0;
-  unsigned int vector_var_index = 0;
   for (unsigned int i = 0; i < number_of_variables; i++)
     {
       variable_info varInfo;
@@ -64,28 +62,9 @@ userInputParameters<dim>::loadVariableAttributes(
 
       varInfo.global_var_index = i;
 
-      !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing)
-        ? varInfo.var_needed = true
-        : varInfo.var_needed = false;
+      varInfo.var_needed = !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing);
 
-      if (var_type[i] == SCALAR)
-        {
-          varInfo.is_scalar = true;
-          if (varInfo.var_needed)
-            {
-              varInfo.variable_index = scalar_var_index;
-              scalar_var_index++;
-            }
-        }
-      else
-        {
-          varInfo.is_scalar = false;
-          if (varInfo.var_needed)
-            {
-              varInfo.variable_index = vector_var_index;
-              vector_var_index++;
-            }
-        }
+      varInfo.is_scalar = var_type[i] == SCALAR;
 
       varInfoListExplicitRHS.push_back(varInfo);
     }
@@ -101,8 +80,6 @@ userInputParameters<dim>::loadVariableAttributes(
         }
     }
   varInfoListNonexplicitRHS.reserve(num_var_nonexplicit_RHS);
-  scalar_var_index = 0;
-  vector_var_index = 0;
   for (unsigned int i = 0; i < number_of_variables; i++)
     {
       variable_info varInfo;
@@ -115,28 +92,9 @@ userInputParameters<dim>::loadVariableAttributes(
 
       varInfo.global_var_index = i;
 
-      !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing)
-        ? varInfo.var_needed = true
-        : varInfo.var_needed = false;
+      varInfo.var_needed = !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing);
 
-      if (var_type[i] == SCALAR)
-        {
-          varInfo.is_scalar = true;
-          if (varInfo.var_needed)
-            {
-              varInfo.variable_index = scalar_var_index;
-              scalar_var_index++;
-            }
-        }
-      else
-        {
-          varInfo.is_scalar = false;
-          if (varInfo.var_needed)
-            {
-              varInfo.variable_index = vector_var_index;
-              vector_var_index++;
-            }
-        }
+      varInfo.is_scalar = var_type[i] == SCALAR;
 
       varInfoListNonexplicitRHS.push_back(varInfo);
     }
@@ -153,8 +111,6 @@ userInputParameters<dim>::loadVariableAttributes(
     }
 
   varInfoListLHS.reserve(num_var_LHS);
-  scalar_var_index = 0;
-  vector_var_index = 0;
   for (unsigned int i = 0; i < number_of_variables; i++)
     {
       variable_info varInfo;
@@ -167,35 +123,14 @@ userInputParameters<dim>::loadVariableAttributes(
 
       varInfo.global_var_index = i;
 
-      !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing)
-        ? varInfo.var_needed = true
-        : varInfo.var_needed = false;
+      varInfo.var_needed = !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing);
 
-      if (var_type[i] == SCALAR)
-        {
-          varInfo.is_scalar = true;
-          if (varInfo.var_needed)
-            {
-              varInfo.variable_index = scalar_var_index;
-              scalar_var_index++;
-            }
-        }
-      else
-        {
-          varInfo.is_scalar = false;
-          if (varInfo.var_needed)
-            {
-              varInfo.variable_index = vector_var_index;
-              vector_var_index++;
-            }
-        }
+      varInfo.is_scalar = var_type[i] == SCALAR;
 
       varInfoListLHS.push_back(varInfo);
     }
 
   varChangeInfoListLHS.reserve(num_var_LHS);
-  scalar_var_index = 0;
-  vector_var_index = 0;
   for (unsigned int i = 0; i < number_of_variables; i++)
     {
       variable_info varInfo;
@@ -209,28 +144,9 @@ userInputParameters<dim>::loadVariableAttributes(
 
       varInfo.global_var_index = i;
 
-      !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing)
-        ? varInfo.var_needed = true
-        : varInfo.var_needed = false;
+      varInfo.var_needed = !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing);
 
-      if (var_type[i] == SCALAR)
-        {
-          varInfo.is_scalar = true;
-          if (varInfo.var_needed)
-            {
-              varInfo.variable_index = scalar_var_index;
-              scalar_var_index++;
-            }
-        }
-      else
-        {
-          varInfo.is_scalar = false;
-          if (varInfo.var_needed)
-            {
-              varInfo.variable_index = vector_var_index;
-              vector_var_index++;
-            }
-        }
+      varInfo.is_scalar = var_type[i] == SCALAR;
 
       varChangeInfoListLHS.push_back(varInfo);
     }
@@ -238,8 +154,6 @@ userInputParameters<dim>::loadVariableAttributes(
   // Load variable information for postprocessing
   // First, the info list for the base field variables
   pp_baseVarInfoList.reserve(number_of_variables);
-  scalar_var_index = 0;
-  vector_var_index = 0;
   for (unsigned int i = 0; i < number_of_variables; i++)
     {
       variable_info varInfo;
@@ -249,28 +163,9 @@ userInputParameters<dim>::loadVariableAttributes(
 
       varInfo.global_var_index = i;
 
-      !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing)
-        ? varInfo.var_needed = true
-        : varInfo.var_needed = false;
+      varInfo.var_needed = !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing);
 
-      if (var_type[i] == SCALAR)
-        {
-          varInfo.is_scalar = true;
-          if (varInfo.var_needed)
-            {
-              varInfo.variable_index = scalar_var_index;
-              scalar_var_index++;
-            }
-        }
-      else
-        {
-          varInfo.is_scalar = false;
-          if (varInfo.var_needed)
-            {
-              varInfo.variable_index = vector_var_index;
-              vector_var_index++;
-            }
-        }
+      varInfo.is_scalar = var_type[i] == SCALAR;
 
       pp_baseVarInfoList.push_back(varInfo);
     }
@@ -292,8 +187,6 @@ userInputParameters<dim>::loadVariableAttributes(
 
   // The info list for the postprocessing field variables
   pp_varInfoList.reserve(pp_number_of_variables);
-  scalar_var_index = 0;
-  vector_var_index = 0;
   for (unsigned int i = 0; i < pp_number_of_variables; i++)
     {
       variable_info varInfo;
@@ -303,18 +196,9 @@ userInputParameters<dim>::loadVariableAttributes(
         variable_attributes.equation_dependency_parser.eval_flags_residual_postprocess[i];
 
       varInfo.global_var_index = i;
-      if (pp_var_type[i] == SCALAR)
-        {
-          varInfo.is_scalar      = true;
-          varInfo.variable_index = scalar_var_index;
-          scalar_var_index++;
-        }
-      else
-        {
-          varInfo.is_scalar      = false;
-          varInfo.variable_index = vector_var_index;
-          vector_var_index++;
-        }
+
+      varInfo.is_scalar = pp_var_type[i] == SCALAR;
+
       pp_varInfoList.push_back(varInfo);
     }
 }
diff --git a/src/variableContainer/variableContainer.cc b/src/variableContainer/variableContainer.cc
index 491a449d8..14ac824e8 100644
--- a/src/variableContainer/variableContainer.cc
+++ b/src/variableContainer/variableContainer.cc
@@ -1,44 +1,48 @@
-// All of the methods for the 'variableContainer' class
 #include "../../include/variableContainer.h"
 
 template <int dim, int degree, typename T>
 variableContainer<dim, degree, T>::variableContainer(
   const dealii::MatrixFree<dim, double> &data,
-  std::vector<variable_info>             _varInfoList,
-  std::vector<variable_info>             _varChangeInfoList)
+  const std::vector<variable_info>      &_varInfoList,
+  const std::vector<variable_info>      &_varChangeInfoList)
+  : varInfoList(std::move(_varInfoList))
+  , varChangeInfoList(std::move(_varChangeInfoList))
+  , num_var(varInfoList.size())
 {
-  varInfoList       = _varInfoList;
-  varChangeInfoList = _varChangeInfoList;
-
-  num_var = varInfoList.size();
-
   for (unsigned int i = 0; i < num_var; i++)
     {
-      if (varInfoList[i].var_needed)
+      const auto &var_info        = varInfoList[i];
+      const auto &var_change_info = varChangeInfoList[i];
+
+      if (var_info.var_needed)
         {
-          if (varInfoList[i].is_scalar)
+          const unsigned int var_index = var_info.global_var_index;
+
+          if (var_info.is_scalar)
             {
-              dealii::FEEvaluation<dim, degree, degree + 1, 1, double> var(data, i);
-              scalar_vars.push_back(var);
+              scalar_vars_map.emplace(var_index,
+                                      std::make_unique<scalar_FEEval>(data, i));
             }
           else
             {
-              dealii::FEEvaluation<dim, degree, degree + 1, dim, double> var(data, i);
-              vector_vars.push_back(var);
+              vector_vars_map.emplace(var_index,
+                                      std::make_unique<vector_FEEval>(data, i));
             }
         }
 
-      if (varChangeInfoList[i].var_needed)
+      if (var_change_info.var_needed)
         {
-          if (varChangeInfoList[i].is_scalar)
+          const unsigned int var_index = var_change_info.global_var_index;
+
+          if (var_change_info.is_scalar)
             {
-              dealii::FEEvaluation<dim, degree, degree + 1, 1, double> var(data, i);
-              scalar_change_in_vars.push_back(var);
+              scalar_change_in_vars_map.emplace(var_index,
+                                                std::make_unique<scalar_FEEval>(data, i));
             }
           else
             {
-              dealii::FEEvaluation<dim, degree, degree + 1, dim, double> var(data, i);
-              vector_change_in_vars.push_back(var);
+              vector_change_in_vars_map.emplace(var_index,
+                                                std::make_unique<vector_FEEval>(data, i));
             }
         }
     }
@@ -47,26 +51,28 @@ variableContainer<dim, degree, T>::variableContainer(
 template <int dim, int degree, typename T>
 variableContainer<dim, degree, T>::variableContainer(
   const dealii::MatrixFree<dim, double> &data,
-  std::vector<variable_info>             _varInfoList)
+  const std::vector<variable_info>      &_varInfoList)
+  : varInfoList(std::move(_varInfoList))
+  , num_var(varInfoList.size())
 {
-  varInfoList = _varInfoList;
-
-  num_var = varInfoList.size();
-
   for (unsigned int i = 0; i < num_var; i++)
     {
-      if (varInfoList[i].var_needed)
+      const auto &var_info = varInfoList[i];
+
+      if (!var_info.var_needed)
         {
-          if (varInfoList[i].is_scalar)
-            {
-              dealii::FEEvaluation<dim, degree, degree + 1, 1, double> var(data, i);
-              scalar_vars.push_back(var);
-            }
-          else
-            {
-              dealii::FEEvaluation<dim, degree, degree + 1, dim, double> var(data, i);
-              vector_vars.push_back(var);
-            }
+          continue;
+        }
+
+      const unsigned int var_index = var_info.global_var_index;
+
+      if (var_info.is_scalar)
+        {
+          scalar_vars_map.emplace(var_index, std::make_unique<scalar_FEEval>(data, i));
+        }
+      else
+        {
+          vector_vars_map.emplace(var_index, std::make_unique<vector_FEEval>(data, i));
         }
     }
 }
@@ -76,74 +82,90 @@ variableContainer<dim, degree, T>::variableContainer(
 template <int dim, int degree, typename T>
 variableContainer<dim, degree, T>::variableContainer(
   const dealii::MatrixFree<dim, double> &data,
-  std::vector<variable_info>             _varInfoList,
-  unsigned int                           fixed_index)
+  const std::vector<variable_info>      &_varInfoList,
+  const unsigned int                    &fixed_index)
+  : varInfoList(std::move(_varInfoList))
+  , num_var(varInfoList.size())
 {
-  varInfoList = _varInfoList;
-
-  num_var = varInfoList.size();
-
   for (unsigned int i = 0; i < num_var; i++)
     {
-      if (varInfoList[i].var_needed)
+      const auto &var_info = varInfoList[i];
+
+      if (!var_info.var_needed)
         {
-          if (varInfoList[i].is_scalar)
-            {
-              dealii::FEEvaluation<dim, degree, degree + 1, 1, double> var(data,
-                                                                           fixed_index);
-              scalar_vars.push_back(var);
-            }
-          else
-            {
-              dealii::FEEvaluation<dim, degree, degree + 1, dim, double> var(data,
-                                                                             fixed_index);
-              vector_vars.push_back(var);
-            }
+          continue;
+        }
+
+      const unsigned int var_index = var_info.global_var_index;
+
+      if (var_info.is_scalar)
+        {
+          scalar_vars_map.emplace(var_index,
+                                  std::make_unique<scalar_FEEval>(data, fixed_index));
+        }
+      else
+        {
+          vector_vars_map.emplace(var_index,
+                                  std::make_unique<vector_FEEval>(data, fixed_index));
         }
     }
 }
 
 template <int dim, int degree, typename T>
 unsigned int
-variableContainer<dim, degree, T>::get_num_q_points()
+variableContainer<dim, degree, T>::get_num_q_points() const
 {
-  if (scalar_vars.size() > 0)
+  if (!scalar_vars_map.empty())
+    {
+      return scalar_vars_map.begin()->second->n_q_points;
+    }
+  else if (!vector_vars_map.empty())
     {
-      return scalar_vars[0].n_q_points;
+      return vector_vars_map.begin()->second->n_q_points;
     }
-  else if (vector_vars.size() > 0)
+  else if (!scalar_change_in_vars_map.empty())
     {
-      return vector_vars[0].n_q_points;
+      return scalar_change_in_vars_map.begin()->second->n_q_points;
     }
-  else if (scalar_change_in_vars.size() > 0)
+  else if (!vector_change_in_vars_map.empty())
     {
-      return scalar_change_in_vars[0].n_q_points;
+      return vector_change_in_vars_map.begin()->second->n_q_points;
     }
   else
     {
-      return vector_change_in_vars[0].n_q_points;
+      AssertThrow(false,
+                  dealii::ExcMessage(
+                    "PRISMS-PF Error: When trying to access the number of quadrature "
+                    "points, all FEEvaluation object containers were empty."));
     }
 }
 
 template <int dim, int degree, typename T>
 dealii::Point<dim, T>
-variableContainer<dim, degree, T>::get_q_point_location()
+variableContainer<dim, degree, T>::get_q_point_location() const
 {
-  if (scalar_vars.size() > 0)
+  if (!scalar_vars_map.empty())
+    {
+      return scalar_vars_map.begin()->second->quadrature_point(q_point);
+    }
+  else if (!vector_vars_map.empty())
     {
-      return scalar_vars[0].quadrature_point(q_point);
+      return vector_vars_map.begin()->second->quadrature_point(q_point);
     }
-  else if (vector_vars.size() > 0)
+  else if (!scalar_change_in_vars_map.empty())
     {
-      return vector_vars[0].quadrature_point(q_point);
+      return scalar_change_in_vars_map.begin()->second->quadrature_point(q_point);
     }
-  else if (scalar_change_in_vars.size() > 0)
+  else if (!vector_change_in_vars_map.empty())
     {
-      return scalar_change_in_vars[0].quadrature_point(q_point);
+      return vector_change_in_vars_map.begin()->second->quadrature_point(q_point);
     }
   else
     {
-      return vector_change_in_vars[0].quadrature_point(q_point);
+      AssertThrow(false,
+                  dealii::ExcMessage(
+                    "PRISMS-PF Error: When trying to access the quadrature point "
+                    "location, all FEEvaluation object containers were empty."));
     }
 }
 
@@ -154,29 +176,35 @@ variableContainer<dim, degree, T>::reinit_and_eval(const std::vector<vectorType
 {
   for (unsigned int i = 0; i < num_var; i++)
     {
-      if (varInfoList[i].var_needed)
+      const auto &var_info = varInfoList[i];
+
+      if (!var_info.var_needed)
         {
-          if (varInfoList[i].is_scalar)
-            {
-              scalar_vars[varInfoList[i].variable_index].reinit(cell);
-              scalar_vars[varInfoList[i].variable_index].read_dof_values(*src[i]);
-              scalar_vars[varInfoList[i].variable_index].evaluate(
-                varInfoList[i].evaluation_flags);
-            }
-          else
-            {
-              vector_vars[varInfoList[i].variable_index].reinit(cell);
-              vector_vars[varInfoList[i].variable_index].read_dof_values(*src[i]);
-              vector_vars[varInfoList[i].variable_index].evaluate(
-                varInfoList[i].evaluation_flags);
-            }
+          continue;
+        }
+
+      const unsigned int var_index = var_info.global_var_index;
+
+      if (var_info.is_scalar)
+        {
+          auto *scalar_FEEval_ptr = scalar_vars_map[var_index].get();
+          scalar_FEEval_ptr->reinit(cell);
+          scalar_FEEval_ptr->read_dof_values(*src[i]);
+          scalar_FEEval_ptr->evaluate(var_info.evaluation_flags);
+        }
+      else
+        {
+          auto *vector_FEEval_ptr = vector_vars_map[var_index].get();
+          vector_FEEval_ptr->reinit(cell);
+          vector_FEEval_ptr->read_dof_values(*src[i]);
+          vector_FEEval_ptr->evaluate(var_info.evaluation_flags);
         }
     }
 }
 
 /**
- * This is specialized for the LHS where there will be only one change in the
- * solution needed. The RHS method takes the src as a vector of vectorTypes.
+ * This is specialized for the LHS where a change in solution is needed. The RHS method
+ * takes the src as a vector of vectorTypes.
  */
 template <int dim, int degree, typename T>
 void
@@ -187,82 +215,42 @@ variableContainer<dim, degree, T>::reinit_and_eval_change_in_solution(
 {
   if (varChangeInfoList[var_being_solved].is_scalar)
     {
-      scalar_change_in_vars[0].reinit(cell);
-      scalar_change_in_vars[0].read_dof_values(src);
-      scalar_change_in_vars[0].evaluate(
-        varChangeInfoList[var_being_solved].evaluation_flags);
+      auto *scalar_FEEval_ptr = scalar_change_in_vars_map[var_being_solved].get();
+      scalar_FEEval_ptr->reinit(cell);
+      scalar_FEEval_ptr->read_dof_values(src);
+      scalar_FEEval_ptr->evaluate(varChangeInfoList[var_being_solved].evaluation_flags);
     }
   else
     {
-      vector_change_in_vars[0].reinit(cell);
-      vector_change_in_vars[0].read_dof_values(src);
-      vector_change_in_vars[0].evaluate(
-        varChangeInfoList[var_being_solved].evaluation_flags);
+      auto *vector_FEEval_ptr = vector_change_in_vars_map[var_being_solved].get();
+      vector_FEEval_ptr->reinit(cell);
+      vector_FEEval_ptr->read_dof_values(src);
+      vector_FEEval_ptr->evaluate(varChangeInfoList[var_being_solved].evaluation_flags);
     }
 }
 
 template <int dim, int degree, typename T>
 void
-variableContainer<dim, degree, T>::reinit_and_eval_LHS(
-  const vectorType               &src,
-  const std::vector<vectorType *> solutionSet,
-  unsigned int                    cell,
-  unsigned int                    var_being_solved)
+variableContainer<dim, degree, T>::reinit(unsigned int cell)
 {
   for (unsigned int i = 0; i < num_var; i++)
     {
-      if (varInfoList[i].var_needed)
+      const auto &var_info = varInfoList[i];
+
+      if (!var_info.var_needed)
         {
-          if (varInfoList[i].is_scalar)
-            {
-              scalar_vars[varInfoList[i].variable_index].reinit(cell);
-              if (i == var_being_solved)
-                {
-                  scalar_vars[varInfoList[i].variable_index].read_dof_values(src);
-                }
-              else
-                {
-                  scalar_vars[varInfoList[i].variable_index].read_dof_values(
-                    *solutionSet[i]);
-                }
-              scalar_vars[varInfoList[i].variable_index].evaluate(
-                varInfoList[i].evaluation_flags);
-            }
-          else
-            {
-              vector_vars[varInfoList[i].variable_index].reinit(cell);
-              if (i == var_being_solved)
-                {
-                  vector_vars[varInfoList[i].variable_index].read_dof_values(src);
-                }
-              else
-                {
-                  vector_vars[varInfoList[i].variable_index].read_dof_values(
-                    *solutionSet[i]);
-                }
-              vector_vars[varInfoList[i].variable_index].evaluate(
-                varInfoList[i].evaluation_flags);
-            }
+          continue;
         }
-    }
-}
 
-template <int dim, int degree, typename T>
-void
-variableContainer<dim, degree, T>::reinit(unsigned int cell)
-{
-  for (unsigned int i = 0; i < num_var; i++)
-    {
-      if (varInfoList[i].var_needed)
+      const unsigned int var_index = var_info.global_var_index;
+
+      if (var_info.is_scalar)
         {
-          if (varInfoList[i].is_scalar)
-            {
-              scalar_vars[varInfoList[i].variable_index].reinit(cell);
-            }
-          else
-            {
-              vector_vars[varInfoList[i].variable_index].reinit(cell);
-            }
+          scalar_vars_map[var_index]->reinit(cell);
+        }
+      else
+        {
+          vector_vars_map[var_index]->reinit(cell);
         }
     }
 }
@@ -274,22 +262,26 @@ variableContainer<dim, degree, T>::integrate_and_distribute(
 {
   for (unsigned int i = 0; i < num_var; i++)
     {
-      if (!(varInfoList[i].residual_flags & dealii::EvaluationFlags::nothing))
+      const auto &var_info = varInfoList[i];
+
+      if (var_info.residual_flags & dealii::EvaluationFlags::nothing)
         {
-          if (varInfoList[i].is_scalar)
-            {
-              scalar_vars[varInfoList[i].variable_index].integrate(
-                varInfoList[i].residual_flags);
-              scalar_vars[varInfoList[i].variable_index].distribute_local_to_global(
-                *dst[i]);
-            }
-          else
-            {
-              vector_vars[varInfoList[i].variable_index].integrate(
-                varInfoList[i].residual_flags);
-              vector_vars[varInfoList[i].variable_index].distribute_local_to_global(
-                *dst[i]);
-            }
+          continue;
+        }
+
+      const unsigned int var_index = var_info.global_var_index;
+
+      if (var_info.is_scalar)
+        {
+          auto *scalar_FEEval_ptr = scalar_vars_map[var_index].get();
+          scalar_FEEval_ptr->integrate(var_info.residual_flags);
+          scalar_FEEval_ptr->distribute_local_to_global(*dst[i]);
+        }
+      else
+        {
+          auto *vector_FEEval_ptr = vector_vars_map[var_index].get();
+          vector_FEEval_ptr->integrate(var_info.residual_flags);
+          vector_FEEval_ptr->distribute_local_to_global(*dst[i]);
         }
     }
 }
@@ -303,15 +295,15 @@ variableContainer<dim, degree, T>::integrate_and_distribute_change_in_solution_L
   // integrate
   if (varChangeInfoList[var_being_solved].is_scalar)
     {
-      scalar_change_in_vars[0].integrate(
-        varChangeInfoList[var_being_solved].residual_flags);
-      scalar_change_in_vars[0].distribute_local_to_global(dst);
+      auto *scalar_FEEval_ptr = scalar_change_in_vars_map[var_being_solved].get();
+      scalar_FEEval_ptr->integrate(varChangeInfoList[var_being_solved].residual_flags);
+      scalar_FEEval_ptr->distribute_local_to_global(dst);
     }
   else
     {
-      vector_change_in_vars[0].integrate(
-        varChangeInfoList[var_being_solved].residual_flags);
-      vector_change_in_vars[0].distribute_local_to_global(dst);
+      auto *vector_FEEval_ptr = vector_change_in_vars_map[var_being_solved].get();
+      vector_FEEval_ptr->integrate(varChangeInfoList[var_being_solved].residual_flags);
+      vector_FEEval_ptr->distribute_local_to_global(dst);
     }
 }
 
@@ -325,8 +317,7 @@ variableContainer<dim, degree, T>::get_scalar_value(
   if (varInfoList[global_variable_index].evaluation_flags &
       dealii::EvaluationFlags::values)
     {
-      return scalar_vars[varInfoList[global_variable_index].variable_index].get_value(
-        q_point);
+      return scalar_vars_map.at(global_variable_index)->get_value(q_point);
     }
   else
     {
@@ -346,8 +337,7 @@ variableContainer<dim, degree, T>::get_scalar_gradient(
   if (varInfoList[global_variable_index].evaluation_flags &
       dealii::EvaluationFlags::gradients)
     {
-      return scalar_vars[varInfoList[global_variable_index].variable_index].get_gradient(
-        q_point);
+      return scalar_vars_map.at(global_variable_index)->get_gradient(q_point);
     }
   else
     {
@@ -367,8 +357,7 @@ variableContainer<dim, degree, T>::get_scalar_hessian(
   if (varInfoList[global_variable_index].evaluation_flags &
       dealii::EvaluationFlags::hessians)
     {
-      return scalar_vars[varInfoList[global_variable_index].variable_index].get_hessian(
-        q_point);
+      return scalar_vars_map.at(global_variable_index)->get_hessian(q_point);
     }
   else
     {
@@ -388,8 +377,7 @@ variableContainer<dim, degree, T>::get_vector_value(
   if (varInfoList[global_variable_index].evaluation_flags &
       dealii::EvaluationFlags::values)
     {
-      return vector_vars[varInfoList[global_variable_index].variable_index].get_value(
-        q_point);
+      return vector_vars_map.at(global_variable_index)->get_value(q_point);
     }
   else
     {
@@ -409,8 +397,7 @@ variableContainer<dim, degree, T>::get_vector_gradient(
   if (varInfoList[global_variable_index].evaluation_flags &
       dealii::EvaluationFlags::gradients)
     {
-      return vector_vars[varInfoList[global_variable_index].variable_index].get_gradient(
-        q_point);
+      return vector_vars_map.at(global_variable_index)->get_gradient(q_point);
     }
   else
     {
@@ -430,8 +417,7 @@ variableContainer<dim, degree, T>::get_vector_hessian(
   if (varInfoList[global_variable_index].evaluation_flags &
       dealii::EvaluationFlags::hessians)
     {
-      return vector_vars[varInfoList[global_variable_index].variable_index].get_hessian(
-        q_point);
+      return vector_vars_map.at(global_variable_index)->get_hessian(q_point);
     }
   else
     {
@@ -453,9 +439,7 @@ variableContainer<dim, degree, T>::get_change_in_scalar_value(
   if (varChangeInfoList[global_variable_index].evaluation_flags &
       dealii::EvaluationFlags::values)
     {
-      return scalar_change_in_vars[varChangeInfoList[global_variable_index]
-                                     .variable_index]
-        .get_value(q_point);
+      return scalar_change_in_vars_map.at(global_variable_index)->get_value(q_point);
     }
   else
     {
@@ -475,9 +459,7 @@ variableContainer<dim, degree, T>::get_change_in_scalar_gradient(
   if (varChangeInfoList[global_variable_index].evaluation_flags &
       dealii::EvaluationFlags::gradients)
     {
-      return scalar_change_in_vars[varChangeInfoList[global_variable_index]
-                                     .variable_index]
-        .get_gradient(q_point);
+      return scalar_change_in_vars_map.at(global_variable_index)->get_gradient(q_point);
     }
   else
     {
@@ -497,9 +479,7 @@ variableContainer<dim, degree, T>::get_change_in_scalar_hessian(
   if (varChangeInfoList[global_variable_index].evaluation_flags &
       dealii::EvaluationFlags::hessians)
     {
-      return scalar_change_in_vars[varChangeInfoList[global_variable_index]
-                                     .variable_index]
-        .get_hessian(q_point);
+      return scalar_change_in_vars_map.at(global_variable_index)->get_hessian(q_point);
     }
   else
     {
@@ -519,9 +499,7 @@ variableContainer<dim, degree, T>::get_change_in_vector_value(
   if (varChangeInfoList[global_variable_index].evaluation_flags &
       dealii::EvaluationFlags::values)
     {
-      return vector_change_in_vars[varChangeInfoList[global_variable_index]
-                                     .variable_index]
-        .get_value(q_point);
+      return vector_change_in_vars_map.at(global_variable_index)->get_value(q_point);
     }
   else
     {
@@ -541,9 +519,7 @@ variableContainer<dim, degree, T>::get_change_in_vector_gradient(
   if (varChangeInfoList[global_variable_index].evaluation_flags &
       dealii::EvaluationFlags::gradients)
     {
-      return vector_change_in_vars[varChangeInfoList[global_variable_index]
-                                     .variable_index]
-        .get_gradient(q_point);
+      return vector_change_in_vars_map.at(global_variable_index)->get_gradient(q_point);
     }
   else
     {
@@ -563,9 +539,7 @@ variableContainer<dim, degree, T>::get_change_in_vector_hessian(
   if (varChangeInfoList[global_variable_index].evaluation_flags &
       dealii::EvaluationFlags::hessians)
     {
-      return vector_change_in_vars[varChangeInfoList[global_variable_index]
-                                     .variable_index]
-        .get_hessian(q_point);
+      return vector_change_in_vars_map.at(global_variable_index)->get_hessian(q_point);
     }
   else
     {
@@ -584,8 +558,7 @@ variableContainer<dim, degree, T>::set_scalar_value_term_RHS(
   unsigned int global_variable_index,
   T            val)
 {
-  scalar_vars[varInfoList[global_variable_index].variable_index].submit_value(val,
-                                                                              q_point);
+  scalar_vars_map[global_variable_index]->submit_value(val, q_point);
 }
 
 template <int dim, int degree, typename T>
@@ -594,8 +567,7 @@ variableContainer<dim, degree, T>::set_scalar_gradient_term_RHS(
   unsigned int              global_variable_index,
   dealii::Tensor<1, dim, T> grad)
 {
-  scalar_vars[varInfoList[global_variable_index].variable_index].submit_gradient(grad,
-                                                                                 q_point);
+  scalar_vars_map[global_variable_index]->submit_gradient(grad, q_point);
 }
 
 template <int dim, int degree, typename T>
@@ -604,8 +576,7 @@ variableContainer<dim, degree, T>::set_vector_value_term_RHS(
   unsigned int              global_variable_index,
   dealii::Tensor<1, dim, T> val)
 {
-  vector_vars[varInfoList[global_variable_index].variable_index].submit_value(val,
-                                                                              q_point);
+  vector_vars_map[global_variable_index]->submit_value(val, q_point);
 }
 
 template <int dim, int degree, typename T>
@@ -614,8 +585,7 @@ variableContainer<dim, degree, T>::set_vector_gradient_term_RHS(
   unsigned int              global_variable_index,
   dealii::Tensor<2, dim, T> grad)
 {
-  vector_vars[varInfoList[global_variable_index].variable_index].submit_gradient(grad,
-                                                                                 q_point);
+  vector_vars_map[global_variable_index]->submit_gradient(grad, q_point);
 }
 
 template <int dim, int degree, typename T>
@@ -624,8 +594,7 @@ variableContainer<dim, degree, T>::set_scalar_value_term_LHS(
   unsigned int global_variable_index,
   T            val)
 {
-  scalar_change_in_vars[varChangeInfoList[global_variable_index].variable_index]
-    .submit_value(val, q_point);
+  scalar_change_in_vars_map[global_variable_index]->submit_value(val, q_point);
 }
 
 template <int dim, int degree, typename T>
@@ -634,8 +603,7 @@ variableContainer<dim, degree, T>::set_scalar_gradient_term_LHS(
   unsigned int              global_variable_index,
   dealii::Tensor<1, dim, T> grad)
 {
-  scalar_change_in_vars[varChangeInfoList[global_variable_index].variable_index]
-    .submit_gradient(grad, q_point);
+  scalar_change_in_vars_map[global_variable_index]->submit_gradient(grad, q_point);
 }
 
 template <int dim, int degree, typename T>
@@ -644,8 +612,7 @@ variableContainer<dim, degree, T>::set_vector_value_term_LHS(
   unsigned int              global_variable_index,
   dealii::Tensor<1, dim, T> val)
 {
-  vector_change_in_vars[varChangeInfoList[global_variable_index].variable_index]
-    .submit_value(val, q_point);
+  vector_change_in_vars_map[global_variable_index]->submit_value(val, q_point);
 }
 
 template <int dim, int degree, typename T>
@@ -654,8 +621,7 @@ variableContainer<dim, degree, T>::set_vector_gradient_term_LHS(
   unsigned int              global_variable_index,
   dealii::Tensor<2, dim, T> grad)
 {
-  vector_change_in_vars[varChangeInfoList[global_variable_index].variable_index]
-    .submit_gradient(grad, q_point);
+  vector_change_in_vars_map[global_variable_index]->submit_gradient(grad, q_point);
 }
 
 template class variableContainer<2, 1, dealii::VectorizedArray<double>>;