diff --git a/include/numerics/wrapped_functor.h b/include/numerics/wrapped_functor.h
index dba9b90285..f2a6490a11 100644
--- a/include/numerics/wrapped_functor.h
+++ b/include/numerics/wrapped_functor.h
@@ -123,14 +123,28 @@ class WrappedFunctor : public FEMFunctionBase<Output>
 template <typename Output>
 void WrappedFunctor<Output>::init_context (const FEMContext & c)
 {
-  for (auto dim : c.elem_dimensions())
+  // If the context was constructed with an active var subset (e.g. projector
+  // projecting only selected variables), then only those FE objects are
+  // guaranteed to exist.
+  const std::vector<unsigned int> * active = c.active_vars();
+
+  auto check_and_init = [&c](unsigned int v, unsigned short dim)
+    {
+      FEAbstract * fe = nullptr;
+      c.get_element_fe(v, fe, dim);
+      libmesh_assert_msg(fe, "FEAbstract pointer is null for variable "
+                             << v << " in dimension " << dim);
+      fe->get_nothing();
+    };
+
+  for (const auto dim : c.elem_dimensions())
     {
-      for (auto v : make_range(c.n_vars()))
-        {
-          FEAbstract * fe;
-          c.get_element_fe(v, fe, dim);
-          fe->get_nothing();
-        }
+      if (active)
+        for (const auto v : *active)
+          check_and_init(v, dim);
+      else
+        for (const auto v : make_range(c.n_vars()))
+          check_and_init(v, dim);
     }
 }
 
diff --git a/include/systems/system.h b/include/systems/system.h
index f2208f1269..d4b9c50aaa 100644
--- a/include/systems/system.h
+++ b/include/systems/system.h
@@ -40,6 +40,7 @@
 #include <string_view>
 #include <vector>
 #include <memory>
+#include <optional>
 
 // This define may be useful in --disable-optional builds when it is
 // possible that libmesh will not have any solvers available.
@@ -514,9 +515,15 @@ class System : public ReferenceCountedObject<System>,
    * user-provided cloneable functors.
    * A gradient \p g is only required/used for projecting onto finite
    * element spaces with continuous derivatives.
+   * elem_range \p active_local_range, if provided, indicates the range of elements
+   * over which to perform the projection.
+   * variable_numbers \p variable_numbers, if provided, indicates the variable numbers
+   * onto which to project.
    */
   void project_solution (FunctionBase<Number> * f,
-                         FunctionBase<Gradient> * g = nullptr) const;
+                         FunctionBase<Gradient> * g = nullptr,
+                         std::optional<ConstElemRange> active_local_range = std::nullopt,
+                         std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;
 
   /**
    * Projects arbitrary functions onto the current solution.
@@ -524,9 +531,15 @@ class System : public ReferenceCountedObject<System>,
    * user-provided cloneable functors.
    * A gradient \p g is only required/used for projecting onto finite
    * element spaces with continuous derivatives.
+   * elem_range \p active_local_range, if provided, indicates the range of elements
+   * over which to perform the projection.
+   * variable_numbers \p variable_numbers, if provided, indicates the variable numbers
+   * onto which to project.
    */
   void project_solution (FEMFunctionBase<Number> * f,
-                         FEMFunctionBase<Gradient> * g = nullptr) const;
+                         FEMFunctionBase<Gradient> * g = nullptr,
+                         std::optional<ConstElemRange> active_local_range = std::nullopt,
+                         std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;
 
   /**
    * Projects arbitrary functions onto the current solution.
@@ -545,7 +558,9 @@ class System : public ReferenceCountedObject<System>,
                                               const std::string & unknown_name);
   void project_solution (ValueFunctionPointer fptr,
                          GradientFunctionPointer gptr,
-                         const Parameters & parameters) const;
+                         const Parameters & parameters,
+                         std::optional<ConstElemRange> active_local_range = std::nullopt,
+                         std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;
 
   /**
    * Projects arbitrary functions onto a vector of degree of freedom
@@ -554,6 +569,10 @@ class System : public ReferenceCountedObject<System>,
    * user-provided cloneable functors.
    * A gradient \p g is only required/used for projecting onto finite
    * element spaces with continuous derivatives.
+   * elem_range \p active_local_range, if provided, indicates the range of elements
+   * over which to perform the projection.
+   * variable_numbers \p variable_numbers, if provided, indicates the variable numbers
+   * onto which to project.
    *
    * Constrain the new vector using the requested adjoint rather than
    * primal constraints if is_adjoint is non-negative.
@@ -561,7 +580,9 @@ class System : public ReferenceCountedObject<System>,
   void project_vector (NumericVector<Number> & new_vector,
                        FunctionBase<Number> * f,
                        FunctionBase<Gradient> * g = nullptr,
-                       int is_adjoint = -1) const;
+                       int is_adjoint = -1,
+                       std::optional<ConstElemRange> active_local_range = std::nullopt,
+                       std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;
 
   /**
    * Projects arbitrary functions onto a vector of degree of freedom
@@ -570,6 +591,10 @@ class System : public ReferenceCountedObject<System>,
    * user-provided cloneable functors.
    * A gradient \p g is only required/used for projecting onto finite
    * element spaces with continuous derivatives.
+   * elem_range \p active_local_range, if provided, indicates the range of elements
+   * over which to perform the projection.
+   * variable_numbers \p variable_numbers, if provided, indicates the variable numbers
+   * onto which to project.
    *
    * Constrain the new vector using the requested adjoint rather than
    * primal constraints if is_adjoint is non-negative.
@@ -577,7 +602,9 @@ class System : public ReferenceCountedObject<System>,
   void project_vector (NumericVector<Number> & new_vector,
                        FEMFunctionBase<Number> * f,
                        FEMFunctionBase<Gradient> * g = nullptr,
-                       int is_adjoint = -1) const;
+                       int is_adjoint = -1,
+                       std::optional<ConstElemRange> active_local_range = std::nullopt,
+                       std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;
 
   /**
    * Projects arbitrary functions onto a vector of degree of freedom
@@ -586,6 +613,10 @@ class System : public ReferenceCountedObject<System>,
    * represented by function pointers.
    * A gradient \p gptr is only required/used for projecting onto
    * finite element spaces with continuous derivatives.
+   * elem_range \p active_local_range, if provided, indicates the range of elements
+   * over which to perform the projection.
+   * variable_numbers \p variable_numbers, if provided, indicates the variable numbers
+   * onto which to project.
    *
    * Constrain the new vector using the requested adjoint rather than
    * primal constraints if is_adjoint is non-negative.
@@ -594,7 +625,9 @@ class System : public ReferenceCountedObject<System>,
                        GradientFunctionPointer gptr,
                        const Parameters & parameters,
                        NumericVector<Number> & new_vector,
-                       int is_adjoint = -1) const;
+                       int is_adjoint = -1,
+                       std::optional<ConstElemRange> active_local_range = std::nullopt,
+                       std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;
 
   /**
    * Projects arbitrary boundary functions onto a vector of degree of
@@ -607,11 +640,14 @@ class System : public ReferenceCountedObject<System>,
    * user-provided cloneable functors.
    * A gradient \p g is only required/used for projecting onto finite
    * element spaces with continuous derivatives.
+   * elem_range \p active_local_range, if provided, indicates the range of elements
+   * over which to perform the projection.
    */
   void boundary_project_solution (const std::set<boundary_id_type> & b,
                                   const std::vector<unsigned int> & variables,
                                   FunctionBase<Number> * f,
-                                  FunctionBase<Gradient> * g = nullptr);
+                                  FunctionBase<Gradient> * g = nullptr,
+                                  std::optional<ConstElemRange> active_local_range = std::nullopt);
 
   /**
    * Projects arbitrary boundary functions onto a vector of degree of
@@ -624,12 +660,15 @@ class System : public ReferenceCountedObject<System>,
    * represented by function pointers.
    * A gradient \p gptr is only required/used for projecting onto
    * finite element spaces with continuous derivatives.
+   * elem_range \p active_local_range, if provided, indicates the range of elements
+   * over which to perform the projection.
    */
   void boundary_project_solution (const std::set<boundary_id_type> & b,
                                   const std::vector<unsigned int> & variables,
                                   ValueFunctionPointer fptr,
                                   GradientFunctionPointer gptr,
-                                  const Parameters & parameters);
+                                  const Parameters & parameters,
+                                  std::optional<ConstElemRange> active_local_range = std::nullopt);
 
   /**
    * Projects arbitrary boundary functions onto a vector of degree of
@@ -642,6 +681,8 @@ class System : public ReferenceCountedObject<System>,
    * user-provided cloneable functors.
    * A gradient \p g is only required/used for projecting onto finite
    * element spaces with continuous derivatives.
+   * elem_range \p active_local_range, if provided, indicates the range of elements
+   * over which to perform the projection.
    *
    * Constrain the new vector using the requested adjoint rather than
    * primal constraints if is_adjoint is non-negative.
@@ -651,7 +692,8 @@ class System : public ReferenceCountedObject<System>,
                                 NumericVector<Number> & new_vector,
                                 FunctionBase<Number> * f,
                                 FunctionBase<Gradient> * g = nullptr,
-                                int is_adjoint = -1) const;
+                                int is_adjoint = -1,
+                                std::optional<ConstElemRange> active_local_range = std::nullopt) const;
 
   /**
    * Projects arbitrary boundary functions onto a vector of degree of
@@ -664,6 +706,8 @@ class System : public ReferenceCountedObject<System>,
    * represented by function pointers.
    * A gradient \p gptr is only required/used for projecting onto
    * finite element spaces with continuous derivatives.
+   * elem_range \p active_local_range, if provided, indicates the range of elements
+   * over which to perform the projection.
    *
    * Constrain the new vector using the requested adjoint rather than
    * primal constraints if is_adjoint is non-negative.
@@ -674,7 +718,8 @@ class System : public ReferenceCountedObject<System>,
                                 GradientFunctionPointer gptr,
                                 const Parameters & parameters,
                                 NumericVector<Number> & new_vector,
-                                int is_adjoint = -1) const;
+                                int is_adjoint = -1,
+                                std::optional<ConstElemRange> active_local_range = std::nullopt) const;
 
   /**
    * \returns The system number.
@@ -1984,7 +2029,9 @@ class System : public ReferenceCountedObject<System>,
    * primal constraints if is_adjoint is non-negative.
    */
   void project_vector (NumericVector<Number> &,
-                       int is_adjoint = -1) const;
+                       int is_adjoint = -1,
+                       std::optional<ConstElemRange> active_local_range = std::nullopt,
+                       std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;
 
   /**
    * Projects the vector defined on the old mesh onto the
@@ -1996,7 +2043,9 @@ class System : public ReferenceCountedObject<System>,
    */
   void project_vector (const NumericVector<Number> &,
                        NumericVector<Number> &,
-                       int is_adjoint = -1) const;
+                       int is_adjoint = -1,
+                       std::optional<ConstElemRange> active_local_range = std::nullopt,
+                       std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;
 
   /**
    * Whether this object should condense out constrained degrees of freedom
diff --git a/src/systems/system_projection.C b/src/systems/system_projection.C
index 3970fd4a44..06290fe0cf 100644
--- a/src/systems/system_projection.C
+++ b/src/systems/system_projection.C
@@ -245,7 +245,9 @@ public:
 // ------------------------------------------------------------
 // System implementation
 void System::project_vector (NumericVector<Number> & vector,
-                             int is_adjoint) const
+                             int is_adjoint,
+                             std::optional<ConstElemRange> active_local_range,
+                             std::optional<std::vector<unsigned int>> variable_numbers) const
 {
   // Create a copy of the vector, which currently
   // contains the old data.
@@ -253,7 +255,7 @@ void System::project_vector (NumericVector<Number> & vector,
     old_vector (vector.clone());
 
   // Project the old vector to the new vector
-  this->project_vector (*old_vector, vector, is_adjoint);
+  this->project_vector (*old_vector, vector, is_adjoint, active_local_range, variable_numbers);
 }
 
 
@@ -264,7 +266,9 @@ void System::project_vector (NumericVector<Number> & vector,
  */
 void System::project_vector (const NumericVector<Number> & old_v,
                              NumericVector<Number> & new_v,
-                             int is_adjoint) const
+                             int is_adjoint,
+                             std::optional<ConstElemRange> active_local_range,
+                             std::optional<std::vector<unsigned int>> variable_numbers) const
 {
   LOG_SCOPE ("project_vector(old,new)", "System");
 
@@ -285,9 +289,12 @@ void System::project_vector (const NumericVector<Number> & old_v,
   std::unique_ptr<NumericVector<Number>> local_old_vector_built;
   const NumericVector<Number> * old_vector_ptr = nullptr;
 
-  ConstElemRange active_local_elem_range
-    (this->get_mesh().active_local_elements_begin(),
-     this->get_mesh().active_local_elements_end());
+  if (!active_local_range)
+    {
+      active_local_range.emplace
+        (this->get_mesh().active_local_elements_begin(),
+         this->get_mesh().active_local_elements_end());
+    }
 
   // If the old vector was uniprocessor, make the new
   // vector uniprocessor
@@ -304,7 +311,7 @@ void System::project_vector (const NumericVector<Number> & old_v,
     {
       // Build a send list for efficient localization
       BuildProjectionList projection_list(*this);
-      Threads::parallel_reduce (active_local_elem_range,
+      Threads::parallel_reduce (active_local_range.value(),
                                 projection_list);
 
       // Create a sorted, unique send_list
@@ -328,7 +335,7 @@ void System::project_vector (const NumericVector<Number> & old_v,
     {
       // Build a send list for efficient localization
       BuildProjectionList projection_list(*this);
-      Threads::parallel_reduce (active_local_elem_range,
+      Threads::parallel_reduce (active_local_range.value(),
                                 projection_list);
 
       // Create a sorted, unique send_list
@@ -362,8 +369,22 @@ void System::project_vector (const NumericVector<Number> & old_v,
 
   if (n_variables)
     {
-      std::vector<unsigned int> vars(n_variables);
-      std::iota(vars.begin(), vars.end(), 0);
+      std::vector<unsigned int> vars;
+      if (variable_numbers)
+        {
+          vars = *variable_numbers;
+          for (auto v : vars)
+            if (v >= n_variables)
+              libmesh_error_msg("ERROR: variable number " << v <<
+                               " out of range for system with " <<
+                               n_variables << " variables.");
+        }
+      else
+        {
+          vars.resize(n_variables);
+          std::iota(vars.begin(), vars.end(), 0);
+        }
+
       std::vector<unsigned int> regular_vars, vector_vars;
       for (auto var : vars)
       {
@@ -389,7 +410,7 @@ void System::project_vector (const NumericVector<Number> & old_v,
             g(*this, old_vector, &regular_vars);
 
           FEMProjector projector(*this, f, &g, setter, regular_vars);
-          projector.project(active_local_elem_range);
+          projector.project(active_local_range.value());
         }
 
       if (!vector_vars.empty())
@@ -403,7 +424,7 @@ void System::project_vector (const NumericVector<Number> & old_v,
           OldSolutionValue<Tensor, &FEMContext::point_gradient> g_vector(*this, old_vector, &vector_vars);
 
           FEMVectorProjector vector_projector(*this, f_vector, &g_vector, setter, vector_vars);
-          vector_projector.project(active_local_elem_range);
+          vector_projector.project(active_local_range.value());
         }
 
       // Copy the SCALAR dofs from old_vector to new_vector
@@ -412,7 +433,7 @@ void System::project_vector (const NumericVector<Number> & old_v,
       if (this->processor_id() == (this->n_processors()-1))
         {
           const DofMap & dof_map = this->get_dof_map();
-          for (auto var : make_range(this->n_vars()))
+          for (auto var : vars)
             if (this->variable(var).type().family == SCALAR)
               {
                 // We can just map SCALAR dofs directly across
@@ -1026,11 +1047,13 @@ void System::projection_matrix (SparseMatrix<Number> & proj_mat) const
  */
 void System::project_solution (ValueFunctionPointer fptr,
                                GradientFunctionPointer gptr,
-                               const Parameters & function_parameters) const
+                               const Parameters & function_parameters,
+                               std::optional<ConstElemRange> active_local_range,
+                               std::optional<std::vector<unsigned int>> variable_numbers) const
 {
   WrappedFunction<Number> f(*this, fptr, &function_parameters);
   WrappedFunction<Gradient> g(*this, gptr, &function_parameters);
-  this->project_solution(&f, &g);
+  this->project_solution(&f, &g, active_local_range, variable_numbers);
 }
 
 
@@ -1039,9 +1062,11 @@ void System::project_solution (ValueFunctionPointer fptr,
  * projections and nodal interpolations on each element.
  */
 void System::project_solution (FunctionBase<Number> * f,
-                               FunctionBase<Gradient> * g) const
+                               FunctionBase<Gradient> * g,
+                               std::optional<ConstElemRange> active_local_range,
+                               std::optional<std::vector<unsigned int>> variable_numbers) const
 {
-  this->project_vector(*solution, f, g);
+  this->project_vector(*solution, f, g, /*is_adjoint=*/-1, active_local_range, variable_numbers);
 
   solution->localize(*current_local_solution, _dof_map->get_send_list());
 }
@@ -1052,9 +1077,11 @@ void System::project_solution (FunctionBase<Number> * f,
  * projections and nodal interpolations on each element.
  */
 void System::project_solution (FEMFunctionBase<Number> * f,
-                               FEMFunctionBase<Gradient> * g) const
+                               FEMFunctionBase<Gradient> * g,
+                               std::optional<ConstElemRange> active_local_range,
+                               std::optional<std::vector<unsigned int>> variable_numbers) const
 {
-  this->project_vector(*solution, f, g);
+  this->project_vector(*solution, f, g, /*is_adjoint=*/-1, active_local_range, variable_numbers);
 
   solution->localize(*current_local_solution, _dof_map->get_send_list());
 }
@@ -1068,11 +1095,13 @@ void System::project_vector (ValueFunctionPointer fptr,
                              GradientFunctionPointer gptr,
                              const Parameters & function_parameters,
                              NumericVector<Number> & new_vector,
-                             int is_adjoint) const
+                             int is_adjoint,
+                             std::optional<ConstElemRange> active_local_range,
+                             std::optional<std::vector<unsigned int>> variable_numbers) const
 {
   WrappedFunction<Number> f(*this, fptr, &function_parameters);
   WrappedFunction<Gradient> g(*this, gptr, &function_parameters);
-  this->project_vector(new_vector, &f, &g, is_adjoint);
+  this->project_vector(new_vector, &f, &g, is_adjoint, active_local_range, variable_numbers);
 }
 
 /**
@@ -1082,7 +1111,9 @@ void System::project_vector (ValueFunctionPointer fptr,
 void System::project_vector (NumericVector<Number> & new_vector,
                              FunctionBase<Number> * f,
                              FunctionBase<Gradient> * g,
-                             int is_adjoint) const
+                             int is_adjoint,
+                             std::optional<ConstElemRange> active_local_range,
+                             std::optional<std::vector<unsigned int>> variable_numbers) const
 {
   LOG_SCOPE ("project_vector(FunctionBase)", "System");
 
@@ -1094,10 +1125,10 @@ void System::project_vector (NumericVector<Number> & new_vector,
     {
       WrappedFunctor<Gradient> g_fem(*g);
 
-      this->project_vector(new_vector, &f_fem, &g_fem, is_adjoint);
+      this->project_vector(new_vector, &f_fem, &g_fem, is_adjoint, active_local_range, variable_numbers);
     }
   else
-    this->project_vector(new_vector, &f_fem, nullptr, is_adjoint);
+    this->project_vector(new_vector, &f_fem, nullptr, is_adjoint, active_local_range, variable_numbers);
 }
 
 
@@ -1108,22 +1139,41 @@ void System::project_vector (NumericVector<Number> & new_vector,
 void System::project_vector (NumericVector<Number> & new_vector,
                              FEMFunctionBase<Number> * f,
                              FEMFunctionBase<Gradient> * g,
-                             int is_adjoint) const
+                             int is_adjoint,
+                             std::optional<ConstElemRange> active_local_range,
+                             std::optional<std::vector<unsigned int>> variable_numbers) const
 {
   LOG_SCOPE ("project_fem_vector()", "System");
 
   libmesh_assert (f);
 
-  ConstElemRange active_local_range
-    (this->get_mesh().active_local_elements_begin(),
-     this->get_mesh().active_local_elements_end() );
+  if (!active_local_range)
+    {
+      active_local_range.emplace
+        (this->get_mesh().active_local_elements_begin(),
+         this->get_mesh().active_local_elements_end());
+    }
 
   VectorSetAction<Number> setter(new_vector);
 
   const unsigned int n_variables = this->n_vars();
 
-  std::vector<unsigned int> vars(n_variables);
-  std::iota(vars.begin(), vars.end(), 0);
+  std::vector<unsigned int> vars;
+  if (variable_numbers)
+    {
+      vars = *variable_numbers;
+      for (auto v : vars)
+        if (v >= n_variables)
+          libmesh_error_msg("ERROR: variable number " << v <<
+                           " out of range for system with " <<
+                           n_variables << " variables.");
+    }
+  else
+    {
+      vars.resize(n_variables);
+      std::iota(vars.begin(), vars.end(), 0);
+    }
+
 
   // Use a typedef to make the calling sequence for parallel_for() a bit more readable
   typedef
@@ -1137,12 +1187,12 @@ void System::project_vector (NumericVector<Number> & new_vector,
       FEMFunctionWrapper<Gradient> gw(*g);
 
       FEMProjector projector(*this, fw, &gw, setter, vars);
-      projector.project(active_local_range);
+      projector.project(active_local_range.value());
     }
   else
     {
       FEMProjector projector(*this, fw, nullptr, setter, vars);
-      projector.project(active_local_range);
+      projector.project(active_local_range.value());
     }
 
   // Also, load values into the SCALAR dofs
@@ -1154,7 +1204,7 @@ void System::project_vector (NumericVector<Number> & new_vector,
       FEMContext context( *this );
 
       const DofMap & dof_map = this->get_dof_map();
-      for (auto var : make_range(this->n_vars()))
+      for (auto var : vars)
         if (this->variable(var).type().family == SCALAR)
           {
             // FIXME: We reinit with an arbitrary element in case the user
@@ -1197,7 +1247,7 @@ void System::project_vector (NumericVector<Number> & new_vector,
     {
       // Look for a spline node: a NodeElem with a rational variable
       // on it.
-      for (auto & elem : active_local_range)
+      for (auto & elem : active_local_range.value())
         if (elem->type() == NODEELEM)
           for (auto rational_var : rational_vars)
             if (rational_var->active_on_subdomain(elem->subdomain_id()))
@@ -1237,11 +1287,13 @@ void System::boundary_project_solution (const std::set<boundary_id_type> & b,
                                         const std::vector<unsigned int> & variables,
                                         ValueFunctionPointer fptr,
                                         GradientFunctionPointer gptr,
-                                        const Parameters & function_parameters)
+                                        const Parameters & function_parameters,
+                                        std::optional<ConstElemRange> active_local_range)
+
 {
   WrappedFunction<Number> f(*this, fptr, &function_parameters);
   WrappedFunction<Gradient> g(*this, gptr, &function_parameters);
-  this->boundary_project_solution(b, variables, &f, &g);
+  this->boundary_project_solution(b, variables, &f, &g, active_local_range);
 }
 
 
@@ -1253,9 +1305,10 @@ void System::boundary_project_solution (const std::set<boundary_id_type> & b,
 void System::boundary_project_solution (const std::set<boundary_id_type> & b,
                                         const std::vector<unsigned int> & variables,
                                         FunctionBase<Number> * f,
-                                        FunctionBase<Gradient> * g)
+                                        FunctionBase<Gradient> * g,
+                                        std::optional<ConstElemRange> active_local_range)
 {
-  this->boundary_project_vector(b, variables, *solution, f, g);
+  this->boundary_project_vector(b, variables, *solution, f, g, -1 /*is_adjoint*/, active_local_range);
 
   solution->localize(*current_local_solution);
 }
@@ -1274,12 +1327,13 @@ void System::boundary_project_vector (const std::set<boundary_id_type> & b,
                                       GradientFunctionPointer gptr,
                                       const Parameters & function_parameters,
                                       NumericVector<Number> & new_vector,
-                                      int is_adjoint) const
+                                      int is_adjoint,
+                                      std::optional<ConstElemRange> active_local_range) const
 {
   WrappedFunction<Number> f(*this, fptr, &function_parameters);
   WrappedFunction<Gradient> g(*this, gptr, &function_parameters);
   this->boundary_project_vector(b, variables, new_vector, &f, &g,
-                                is_adjoint);
+                                is_adjoint, active_local_range);
 }
 
 /**
@@ -1291,13 +1345,20 @@ void System::boundary_project_vector (const std::set<boundary_id_type> & b,
                                       NumericVector<Number> & new_vector,
                                       FunctionBase<Number> * f,
                                       FunctionBase<Gradient> * g,
-                                      int is_adjoint) const
+                                      int is_adjoint,
+                                      std::optional<ConstElemRange> active_local_range) const
 {
   LOG_SCOPE ("boundary_project_vector()", "System");
 
+  if (!active_local_range)
+  {
+    active_local_range.emplace
+      (this->get_mesh().active_local_elements_begin(),
+       this->get_mesh().active_local_elements_end());
+  }
+
   Threads::parallel_for
-    (ConstElemRange (this->get_mesh().active_local_elements_begin(),
-                     this->get_mesh().active_local_elements_end() ),
+    (active_local_range.value(),
      BoundaryProjectSolution(b, variables, *this, f, g,
                              this->get_equation_systems().parameters,
                              new_vector)
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 6f27e83ae1..0c4a85fea3 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -94,6 +94,7 @@ unit_tests_sources = \
   mesh/write_nodeset_data.C \
   mesh/write_edgeset_data.C \
   mesh/write_vec_and_scalar.C \
+  mesh/project_solution_test.C \
   numerics/composite_function_test.C \
   numerics/coupling_matrix_test.C \
   numerics/distributed_vector_test.C \
diff --git a/tests/Makefile.in b/tests/Makefile.in
index 741fc6a7fc..72024a3bdc 100644
--- a/tests/Makefile.in
+++ b/tests/Makefile.in
@@ -224,7 +224,8 @@ am__unit_tests_dbg_SOURCES_DIST = driver.C libmesh_cppunit.h \
 	mesh/mapped_subdomain_partitioner_test.C \
 	mesh/write_elemset_data.C mesh/write_sideset_data.C \
 	mesh/write_nodeset_data.C mesh/write_edgeset_data.C \
-	mesh/write_vec_and_scalar.C numerics/composite_function_test.C \
+	mesh/write_vec_and_scalar.C mesh/project_solution_test.C \
+	numerics/composite_function_test.C \
 	numerics/coupling_matrix_test.C \
 	numerics/distributed_vector_test.C \
 	numerics/eigen_sparse_vector_test.C \
@@ -333,6 +334,7 @@ am__objects_2 = unit_tests_dbg-driver.$(OBJEXT) \
 	mesh/unit_tests_dbg-write_nodeset_data.$(OBJEXT) \
 	mesh/unit_tests_dbg-write_edgeset_data.$(OBJEXT) \
 	mesh/unit_tests_dbg-write_vec_and_scalar.$(OBJEXT) \
+	mesh/unit_tests_dbg-project_solution_test.$(OBJEXT) \
 	numerics/unit_tests_dbg-composite_function_test.$(OBJEXT) \
 	numerics/unit_tests_dbg-coupling_matrix_test.$(OBJEXT) \
 	numerics/unit_tests_dbg-distributed_vector_test.$(OBJEXT) \
@@ -427,7 +429,8 @@ am__unit_tests_devel_SOURCES_DIST = driver.C libmesh_cppunit.h \
 	mesh/mapped_subdomain_partitioner_test.C \
 	mesh/write_elemset_data.C mesh/write_sideset_data.C \
 	mesh/write_nodeset_data.C mesh/write_edgeset_data.C \
-	mesh/write_vec_and_scalar.C numerics/composite_function_test.C \
+	mesh/write_vec_and_scalar.C mesh/project_solution_test.C \
+	numerics/composite_function_test.C \
 	numerics/coupling_matrix_test.C \
 	numerics/distributed_vector_test.C \
 	numerics/eigen_sparse_vector_test.C \
@@ -535,6 +538,7 @@ am__objects_4 = unit_tests_devel-driver.$(OBJEXT) \
 	mesh/unit_tests_devel-write_nodeset_data.$(OBJEXT) \
 	mesh/unit_tests_devel-write_edgeset_data.$(OBJEXT) \
 	mesh/unit_tests_devel-write_vec_and_scalar.$(OBJEXT) \
+	mesh/unit_tests_devel-project_solution_test.$(OBJEXT) \
 	numerics/unit_tests_devel-composite_function_test.$(OBJEXT) \
 	numerics/unit_tests_devel-coupling_matrix_test.$(OBJEXT) \
 	numerics/unit_tests_devel-distributed_vector_test.$(OBJEXT) \
@@ -625,7 +629,8 @@ am__unit_tests_oprof_SOURCES_DIST = driver.C libmesh_cppunit.h \
 	mesh/mapped_subdomain_partitioner_test.C \
 	mesh/write_elemset_data.C mesh/write_sideset_data.C \
 	mesh/write_nodeset_data.C mesh/write_edgeset_data.C \
-	mesh/write_vec_and_scalar.C numerics/composite_function_test.C \
+	mesh/write_vec_and_scalar.C mesh/project_solution_test.C \
+	numerics/composite_function_test.C \
 	numerics/coupling_matrix_test.C \
 	numerics/distributed_vector_test.C \
 	numerics/eigen_sparse_vector_test.C \
@@ -733,6 +738,7 @@ am__objects_6 = unit_tests_oprof-driver.$(OBJEXT) \
 	mesh/unit_tests_oprof-write_nodeset_data.$(OBJEXT) \
 	mesh/unit_tests_oprof-write_edgeset_data.$(OBJEXT) \
 	mesh/unit_tests_oprof-write_vec_and_scalar.$(OBJEXT) \
+	mesh/unit_tests_oprof-project_solution_test.$(OBJEXT) \
 	numerics/unit_tests_oprof-composite_function_test.$(OBJEXT) \
 	numerics/unit_tests_oprof-coupling_matrix_test.$(OBJEXT) \
 	numerics/unit_tests_oprof-distributed_vector_test.$(OBJEXT) \
@@ -823,7 +829,8 @@ am__unit_tests_opt_SOURCES_DIST = driver.C libmesh_cppunit.h \
 	mesh/mapped_subdomain_partitioner_test.C \
 	mesh/write_elemset_data.C mesh/write_sideset_data.C \
 	mesh/write_nodeset_data.C mesh/write_edgeset_data.C \
-	mesh/write_vec_and_scalar.C numerics/composite_function_test.C \
+	mesh/write_vec_and_scalar.C mesh/project_solution_test.C \
+	numerics/composite_function_test.C \
 	numerics/coupling_matrix_test.C \
 	numerics/distributed_vector_test.C \
 	numerics/eigen_sparse_vector_test.C \
@@ -931,6 +938,7 @@ am__objects_8 = unit_tests_opt-driver.$(OBJEXT) \
 	mesh/unit_tests_opt-write_nodeset_data.$(OBJEXT) \
 	mesh/unit_tests_opt-write_edgeset_data.$(OBJEXT) \
 	mesh/unit_tests_opt-write_vec_and_scalar.$(OBJEXT) \
+	mesh/unit_tests_opt-project_solution_test.$(OBJEXT) \
 	numerics/unit_tests_opt-composite_function_test.$(OBJEXT) \
 	numerics/unit_tests_opt-coupling_matrix_test.$(OBJEXT) \
 	numerics/unit_tests_opt-distributed_vector_test.$(OBJEXT) \
@@ -1021,7 +1029,8 @@ am__unit_tests_prof_SOURCES_DIST = driver.C libmesh_cppunit.h \
 	mesh/mapped_subdomain_partitioner_test.C \
 	mesh/write_elemset_data.C mesh/write_sideset_data.C \
 	mesh/write_nodeset_data.C mesh/write_edgeset_data.C \
-	mesh/write_vec_and_scalar.C numerics/composite_function_test.C \
+	mesh/write_vec_and_scalar.C mesh/project_solution_test.C \
+	numerics/composite_function_test.C \
 	numerics/coupling_matrix_test.C \
 	numerics/distributed_vector_test.C \
 	numerics/eigen_sparse_vector_test.C \
@@ -1129,6 +1138,7 @@ am__objects_10 = unit_tests_prof-driver.$(OBJEXT) \
 	mesh/unit_tests_prof-write_nodeset_data.$(OBJEXT) \
 	mesh/unit_tests_prof-write_edgeset_data.$(OBJEXT) \
 	mesh/unit_tests_prof-write_vec_and_scalar.$(OBJEXT) \
+	mesh/unit_tests_prof-project_solution_test.$(OBJEXT) \
 	numerics/unit_tests_prof-composite_function_test.$(OBJEXT) \
 	numerics/unit_tests_prof-coupling_matrix_test.$(OBJEXT) \
 	numerics/unit_tests_prof-distributed_vector_test.$(OBJEXT) \
@@ -1394,6 +1404,7 @@ am__depfiles_remade = ./$(DEPDIR)/unit_tests_dbg-driver.Po \
 	mesh/$(DEPDIR)/unit_tests_dbg-mixed_dim_mesh_test.Po \
 	mesh/$(DEPDIR)/unit_tests_dbg-mixed_order_test.Po \
 	mesh/$(DEPDIR)/unit_tests_dbg-nodal_neighbors.Po \
+	mesh/$(DEPDIR)/unit_tests_dbg-project_solution_test.Po \
 	mesh/$(DEPDIR)/unit_tests_dbg-simplex_refinement_test.Po \
 	mesh/$(DEPDIR)/unit_tests_dbg-slit_mesh_test.Po \
 	mesh/$(DEPDIR)/unit_tests_dbg-spatial_dimension_test.Po \
@@ -1433,6 +1444,7 @@ am__depfiles_remade = ./$(DEPDIR)/unit_tests_dbg-driver.Po \
 	mesh/$(DEPDIR)/unit_tests_devel-mixed_dim_mesh_test.Po \
 	mesh/$(DEPDIR)/unit_tests_devel-mixed_order_test.Po \
 	mesh/$(DEPDIR)/unit_tests_devel-nodal_neighbors.Po \
+	mesh/$(DEPDIR)/unit_tests_devel-project_solution_test.Po \
 	mesh/$(DEPDIR)/unit_tests_devel-simplex_refinement_test.Po \
 	mesh/$(DEPDIR)/unit_tests_devel-slit_mesh_test.Po \
 	mesh/$(DEPDIR)/unit_tests_devel-spatial_dimension_test.Po \
@@ -1472,6 +1484,7 @@ am__depfiles_remade = ./$(DEPDIR)/unit_tests_dbg-driver.Po \
 	mesh/$(DEPDIR)/unit_tests_oprof-mixed_dim_mesh_test.Po \
 	mesh/$(DEPDIR)/unit_tests_oprof-mixed_order_test.Po \
 	mesh/$(DEPDIR)/unit_tests_oprof-nodal_neighbors.Po \
+	mesh/$(DEPDIR)/unit_tests_oprof-project_solution_test.Po \
 	mesh/$(DEPDIR)/unit_tests_oprof-simplex_refinement_test.Po \
 	mesh/$(DEPDIR)/unit_tests_oprof-slit_mesh_test.Po \
 	mesh/$(DEPDIR)/unit_tests_oprof-spatial_dimension_test.Po \
@@ -1511,6 +1524,7 @@ am__depfiles_remade = ./$(DEPDIR)/unit_tests_dbg-driver.Po \
 	mesh/$(DEPDIR)/unit_tests_opt-mixed_dim_mesh_test.Po \
 	mesh/$(DEPDIR)/unit_tests_opt-mixed_order_test.Po \
 	mesh/$(DEPDIR)/unit_tests_opt-nodal_neighbors.Po \
+	mesh/$(DEPDIR)/unit_tests_opt-project_solution_test.Po \
 	mesh/$(DEPDIR)/unit_tests_opt-simplex_refinement_test.Po \
 	mesh/$(DEPDIR)/unit_tests_opt-slit_mesh_test.Po \
 	mesh/$(DEPDIR)/unit_tests_opt-spatial_dimension_test.Po \
@@ -1550,6 +1564,7 @@ am__depfiles_remade = ./$(DEPDIR)/unit_tests_dbg-driver.Po \
 	mesh/$(DEPDIR)/unit_tests_prof-mixed_dim_mesh_test.Po \
 	mesh/$(DEPDIR)/unit_tests_prof-mixed_order_test.Po \
 	mesh/$(DEPDIR)/unit_tests_prof-nodal_neighbors.Po \
+	mesh/$(DEPDIR)/unit_tests_prof-project_solution_test.Po \
 	mesh/$(DEPDIR)/unit_tests_prof-simplex_refinement_test.Po \
 	mesh/$(DEPDIR)/unit_tests_prof-slit_mesh_test.Po \
 	mesh/$(DEPDIR)/unit_tests_prof-spatial_dimension_test.Po \
@@ -2293,7 +2308,8 @@ unit_tests_sources = driver.C libmesh_cppunit.h stream_redirector.h \
 	mesh/mapped_subdomain_partitioner_test.C \
 	mesh/write_elemset_data.C mesh/write_sideset_data.C \
 	mesh/write_nodeset_data.C mesh/write_edgeset_data.C \
-	mesh/write_vec_and_scalar.C numerics/composite_function_test.C \
+	mesh/write_vec_and_scalar.C mesh/project_solution_test.C \
+	numerics/composite_function_test.C \
 	numerics/coupling_matrix_test.C \
 	numerics/distributed_vector_test.C \
 	numerics/eigen_sparse_vector_test.C \
@@ -2690,6 +2706,8 @@ mesh/unit_tests_dbg-write_edgeset_data.$(OBJEXT):  \
 	mesh/$(am__dirstamp) mesh/$(DEPDIR)/$(am__dirstamp)
 mesh/unit_tests_dbg-write_vec_and_scalar.$(OBJEXT):  \
 	mesh/$(am__dirstamp) mesh/$(DEPDIR)/$(am__dirstamp)
+mesh/unit_tests_dbg-project_solution_test.$(OBJEXT):  \
+	mesh/$(am__dirstamp) mesh/$(DEPDIR)/$(am__dirstamp)
 numerics/$(am__dirstamp):
 	@$(MKDIR_P) numerics
 	@: >>numerics/$(am__dirstamp)
@@ -2986,6 +3004,8 @@ mesh/unit_tests_devel-write_edgeset_data.$(OBJEXT):  \
 	mesh/$(am__dirstamp) mesh/$(DEPDIR)/$(am__dirstamp)
 mesh/unit_tests_devel-write_vec_and_scalar.$(OBJEXT):  \
 	mesh/$(am__dirstamp) mesh/$(DEPDIR)/$(am__dirstamp)
+mesh/unit_tests_devel-project_solution_test.$(OBJEXT):  \
+	mesh/$(am__dirstamp) mesh/$(DEPDIR)/$(am__dirstamp)
 numerics/unit_tests_devel-composite_function_test.$(OBJEXT):  \
 	numerics/$(am__dirstamp) numerics/$(DEPDIR)/$(am__dirstamp)
 numerics/unit_tests_devel-coupling_matrix_test.$(OBJEXT):  \
@@ -3234,6 +3254,8 @@ mesh/unit_tests_oprof-write_edgeset_data.$(OBJEXT):  \
 	mesh/$(am__dirstamp) mesh/$(DEPDIR)/$(am__dirstamp)
 mesh/unit_tests_oprof-write_vec_and_scalar.$(OBJEXT):  \
 	mesh/$(am__dirstamp) mesh/$(DEPDIR)/$(am__dirstamp)
+mesh/unit_tests_oprof-project_solution_test.$(OBJEXT):  \
+	mesh/$(am__dirstamp) mesh/$(DEPDIR)/$(am__dirstamp)
 numerics/unit_tests_oprof-composite_function_test.$(OBJEXT):  \
 	numerics/$(am__dirstamp) numerics/$(DEPDIR)/$(am__dirstamp)
 numerics/unit_tests_oprof-coupling_matrix_test.$(OBJEXT):  \
@@ -3482,6 +3504,8 @@ mesh/unit_tests_opt-write_edgeset_data.$(OBJEXT):  \
 	mesh/$(am__dirstamp) mesh/$(DEPDIR)/$(am__dirstamp)
 mesh/unit_tests_opt-write_vec_and_scalar.$(OBJEXT):  \
 	mesh/$(am__dirstamp) mesh/$(DEPDIR)/$(am__dirstamp)
+mesh/unit_tests_opt-project_solution_test.$(OBJEXT):  \
+	mesh/$(am__dirstamp) mesh/$(DEPDIR)/$(am__dirstamp)
 numerics/unit_tests_opt-composite_function_test.$(OBJEXT):  \
 	numerics/$(am__dirstamp) numerics/$(DEPDIR)/$(am__dirstamp)
 numerics/unit_tests_opt-coupling_matrix_test.$(OBJEXT):  \
@@ -3730,6 +3754,8 @@ mesh/unit_tests_prof-write_edgeset_data.$(OBJEXT):  \
 	mesh/$(am__dirstamp) mesh/$(DEPDIR)/$(am__dirstamp)
 mesh/unit_tests_prof-write_vec_and_scalar.$(OBJEXT):  \
 	mesh/$(am__dirstamp) mesh/$(DEPDIR)/$(am__dirstamp)
+mesh/unit_tests_prof-project_solution_test.$(OBJEXT):  \
+	mesh/$(am__dirstamp) mesh/$(DEPDIR)/$(am__dirstamp)
 numerics/unit_tests_prof-composite_function_test.$(OBJEXT):  \
 	numerics/$(am__dirstamp) numerics/$(DEPDIR)/$(am__dirstamp)
 numerics/unit_tests_prof-coupling_matrix_test.$(OBJEXT):  \
@@ -4053,6 +4079,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_dbg-mixed_dim_mesh_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_dbg-mixed_order_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_dbg-nodal_neighbors.Po@am__quote@ # am--include-marker
+@AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_dbg-project_solution_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_dbg-simplex_refinement_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_dbg-slit_mesh_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_dbg-spatial_dimension_test.Po@am__quote@ # am--include-marker
@@ -4092,6 +4119,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_devel-mixed_dim_mesh_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_devel-mixed_order_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_devel-nodal_neighbors.Po@am__quote@ # am--include-marker
+@AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_devel-project_solution_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_devel-simplex_refinement_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_devel-slit_mesh_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_devel-spatial_dimension_test.Po@am__quote@ # am--include-marker
@@ -4131,6 +4159,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_oprof-mixed_dim_mesh_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_oprof-mixed_order_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_oprof-nodal_neighbors.Po@am__quote@ # am--include-marker
+@AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_oprof-project_solution_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_oprof-simplex_refinement_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_oprof-slit_mesh_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_oprof-spatial_dimension_test.Po@am__quote@ # am--include-marker
@@ -4170,6 +4199,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_opt-mixed_dim_mesh_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_opt-mixed_order_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_opt-nodal_neighbors.Po@am__quote@ # am--include-marker
+@AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_opt-project_solution_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_opt-simplex_refinement_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_opt-slit_mesh_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_opt-spatial_dimension_test.Po@am__quote@ # am--include-marker
@@ -4209,6 +4239,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_prof-mixed_dim_mesh_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_prof-mixed_order_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_prof-nodal_neighbors.Po@am__quote@ # am--include-marker
+@AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_prof-project_solution_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_prof-simplex_refinement_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_prof-slit_mesh_test.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@mesh/$(DEPDIR)/unit_tests_prof-spatial_dimension_test.Po@am__quote@ # am--include-marker
@@ -5477,6 +5508,20 @@ mesh/unit_tests_dbg-write_vec_and_scalar.obj: mesh/write_vec_and_scalar.C
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_dbg_CPPFLAGS) $(CPPFLAGS) $(unit_tests_dbg_CXXFLAGS) $(CXXFLAGS) -c -o mesh/unit_tests_dbg-write_vec_and_scalar.obj `if test -f 'mesh/write_vec_and_scalar.C'; then $(CYGPATH_W) 'mesh/write_vec_and_scalar.C'; else $(CYGPATH_W) '$(srcdir)/mesh/write_vec_and_scalar.C'; fi`
 
+mesh/unit_tests_dbg-project_solution_test.o: mesh/project_solution_test.C
+@am__fastdepCXX_TRUE@	$(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_dbg_CPPFLAGS) $(CPPFLAGS) $(unit_tests_dbg_CXXFLAGS) $(CXXFLAGS) -MT mesh/unit_tests_dbg-project_solution_test.o -MD -MP -MF mesh/$(DEPDIR)/unit_tests_dbg-project_solution_test.Tpo -c -o mesh/unit_tests_dbg-project_solution_test.o `test -f 'mesh/project_solution_test.C' || echo '$(srcdir)/'`mesh/project_solution_test.C
+@am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) mesh/$(DEPDIR)/unit_tests_dbg-project_solution_test.Tpo mesh/$(DEPDIR)/unit_tests_dbg-project_solution_test.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	$(AM_V_CXX)source='mesh/project_solution_test.C' object='mesh/unit_tests_dbg-project_solution_test.o' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_dbg_CPPFLAGS) $(CPPFLAGS) $(unit_tests_dbg_CXXFLAGS) $(CXXFLAGS) -c -o mesh/unit_tests_dbg-project_solution_test.o `test -f 'mesh/project_solution_test.C' || echo '$(srcdir)/'`mesh/project_solution_test.C
+
+mesh/unit_tests_dbg-project_solution_test.obj: mesh/project_solution_test.C
+@am__fastdepCXX_TRUE@	$(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_dbg_CPPFLAGS) $(CPPFLAGS) $(unit_tests_dbg_CXXFLAGS) $(CXXFLAGS) -MT mesh/unit_tests_dbg-project_solution_test.obj -MD -MP -MF mesh/$(DEPDIR)/unit_tests_dbg-project_solution_test.Tpo -c -o mesh/unit_tests_dbg-project_solution_test.obj `if test -f 'mesh/project_solution_test.C'; then $(CYGPATH_W) 'mesh/project_solution_test.C'; else $(CYGPATH_W) '$(srcdir)/mesh/project_solution_test.C'; fi`
+@am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) mesh/$(DEPDIR)/unit_tests_dbg-project_solution_test.Tpo mesh/$(DEPDIR)/unit_tests_dbg-project_solution_test.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	$(AM_V_CXX)source='mesh/project_solution_test.C' object='mesh/unit_tests_dbg-project_solution_test.obj' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_dbg_CPPFLAGS) $(CPPFLAGS) $(unit_tests_dbg_CXXFLAGS) $(CXXFLAGS) -c -o mesh/unit_tests_dbg-project_solution_test.obj `if test -f 'mesh/project_solution_test.C'; then $(CYGPATH_W) 'mesh/project_solution_test.C'; else $(CYGPATH_W) '$(srcdir)/mesh/project_solution_test.C'; fi`
+
 numerics/unit_tests_dbg-composite_function_test.o: numerics/composite_function_test.C
 @am__fastdepCXX_TRUE@	$(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_dbg_CPPFLAGS) $(CPPFLAGS) $(unit_tests_dbg_CXXFLAGS) $(CXXFLAGS) -MT numerics/unit_tests_dbg-composite_function_test.o -MD -MP -MF numerics/$(DEPDIR)/unit_tests_dbg-composite_function_test.Tpo -c -o numerics/unit_tests_dbg-composite_function_test.o `test -f 'numerics/composite_function_test.C' || echo '$(srcdir)/'`numerics/composite_function_test.C
 @am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) numerics/$(DEPDIR)/unit_tests_dbg-composite_function_test.Tpo numerics/$(DEPDIR)/unit_tests_dbg-composite_function_test.Po
@@ -7143,6 +7188,20 @@ mesh/unit_tests_devel-write_vec_and_scalar.obj: mesh/write_vec_and_scalar.C
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_devel_CPPFLAGS) $(CPPFLAGS) $(unit_tests_devel_CXXFLAGS) $(CXXFLAGS) -c -o mesh/unit_tests_devel-write_vec_and_scalar.obj `if test -f 'mesh/write_vec_and_scalar.C'; then $(CYGPATH_W) 'mesh/write_vec_and_scalar.C'; else $(CYGPATH_W) '$(srcdir)/mesh/write_vec_and_scalar.C'; fi`
 
+mesh/unit_tests_devel-project_solution_test.o: mesh/project_solution_test.C
+@am__fastdepCXX_TRUE@	$(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_devel_CPPFLAGS) $(CPPFLAGS) $(unit_tests_devel_CXXFLAGS) $(CXXFLAGS) -MT mesh/unit_tests_devel-project_solution_test.o -MD -MP -MF mesh/$(DEPDIR)/unit_tests_devel-project_solution_test.Tpo -c -o mesh/unit_tests_devel-project_solution_test.o `test -f 'mesh/project_solution_test.C' || echo '$(srcdir)/'`mesh/project_solution_test.C
+@am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) mesh/$(DEPDIR)/unit_tests_devel-project_solution_test.Tpo mesh/$(DEPDIR)/unit_tests_devel-project_solution_test.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	$(AM_V_CXX)source='mesh/project_solution_test.C' object='mesh/unit_tests_devel-project_solution_test.o' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_devel_CPPFLAGS) $(CPPFLAGS) $(unit_tests_devel_CXXFLAGS) $(CXXFLAGS) -c -o mesh/unit_tests_devel-project_solution_test.o `test -f 'mesh/project_solution_test.C' || echo '$(srcdir)/'`mesh/project_solution_test.C
+
+mesh/unit_tests_devel-project_solution_test.obj: mesh/project_solution_test.C
+@am__fastdepCXX_TRUE@	$(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_devel_CPPFLAGS) $(CPPFLAGS) $(unit_tests_devel_CXXFLAGS) $(CXXFLAGS) -MT mesh/unit_tests_devel-project_solution_test.obj -MD -MP -MF mesh/$(DEPDIR)/unit_tests_devel-project_solution_test.Tpo -c -o mesh/unit_tests_devel-project_solution_test.obj `if test -f 'mesh/project_solution_test.C'; then $(CYGPATH_W) 'mesh/project_solution_test.C'; else $(CYGPATH_W) '$(srcdir)/mesh/project_solution_test.C'; fi`
+@am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) mesh/$(DEPDIR)/unit_tests_devel-project_solution_test.Tpo mesh/$(DEPDIR)/unit_tests_devel-project_solution_test.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	$(AM_V_CXX)source='mesh/project_solution_test.C' object='mesh/unit_tests_devel-project_solution_test.obj' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_devel_CPPFLAGS) $(CPPFLAGS) $(unit_tests_devel_CXXFLAGS) $(CXXFLAGS) -c -o mesh/unit_tests_devel-project_solution_test.obj `if test -f 'mesh/project_solution_test.C'; then $(CYGPATH_W) 'mesh/project_solution_test.C'; else $(CYGPATH_W) '$(srcdir)/mesh/project_solution_test.C'; fi`
+
 numerics/unit_tests_devel-composite_function_test.o: numerics/composite_function_test.C
 @am__fastdepCXX_TRUE@	$(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_devel_CPPFLAGS) $(CPPFLAGS) $(unit_tests_devel_CXXFLAGS) $(CXXFLAGS) -MT numerics/unit_tests_devel-composite_function_test.o -MD -MP -MF numerics/$(DEPDIR)/unit_tests_devel-composite_function_test.Tpo -c -o numerics/unit_tests_devel-composite_function_test.o `test -f 'numerics/composite_function_test.C' || echo '$(srcdir)/'`numerics/composite_function_test.C
 @am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) numerics/$(DEPDIR)/unit_tests_devel-composite_function_test.Tpo numerics/$(DEPDIR)/unit_tests_devel-composite_function_test.Po
@@ -8809,6 +8868,20 @@ mesh/unit_tests_oprof-write_vec_and_scalar.obj: mesh/write_vec_and_scalar.C
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_oprof_CPPFLAGS) $(CPPFLAGS) $(unit_tests_oprof_CXXFLAGS) $(CXXFLAGS) -c -o mesh/unit_tests_oprof-write_vec_and_scalar.obj `if test -f 'mesh/write_vec_and_scalar.C'; then $(CYGPATH_W) 'mesh/write_vec_and_scalar.C'; else $(CYGPATH_W) '$(srcdir)/mesh/write_vec_and_scalar.C'; fi`
 
+mesh/unit_tests_oprof-project_solution_test.o: mesh/project_solution_test.C
+@am__fastdepCXX_TRUE@	$(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_oprof_CPPFLAGS) $(CPPFLAGS) $(unit_tests_oprof_CXXFLAGS) $(CXXFLAGS) -MT mesh/unit_tests_oprof-project_solution_test.o -MD -MP -MF mesh/$(DEPDIR)/unit_tests_oprof-project_solution_test.Tpo -c -o mesh/unit_tests_oprof-project_solution_test.o `test -f 'mesh/project_solution_test.C' || echo '$(srcdir)/'`mesh/project_solution_test.C
+@am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) mesh/$(DEPDIR)/unit_tests_oprof-project_solution_test.Tpo mesh/$(DEPDIR)/unit_tests_oprof-project_solution_test.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	$(AM_V_CXX)source='mesh/project_solution_test.C' object='mesh/unit_tests_oprof-project_solution_test.o' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_oprof_CPPFLAGS) $(CPPFLAGS) $(unit_tests_oprof_CXXFLAGS) $(CXXFLAGS) -c -o mesh/unit_tests_oprof-project_solution_test.o `test -f 'mesh/project_solution_test.C' || echo '$(srcdir)/'`mesh/project_solution_test.C
+
+mesh/unit_tests_oprof-project_solution_test.obj: mesh/project_solution_test.C
+@am__fastdepCXX_TRUE@	$(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_oprof_CPPFLAGS) $(CPPFLAGS) $(unit_tests_oprof_CXXFLAGS) $(CXXFLAGS) -MT mesh/unit_tests_oprof-project_solution_test.obj -MD -MP -MF mesh/$(DEPDIR)/unit_tests_oprof-project_solution_test.Tpo -c -o mesh/unit_tests_oprof-project_solution_test.obj `if test -f 'mesh/project_solution_test.C'; then $(CYGPATH_W) 'mesh/project_solution_test.C'; else $(CYGPATH_W) '$(srcdir)/mesh/project_solution_test.C'; fi`
+@am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) mesh/$(DEPDIR)/unit_tests_oprof-project_solution_test.Tpo mesh/$(DEPDIR)/unit_tests_oprof-project_solution_test.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	$(AM_V_CXX)source='mesh/project_solution_test.C' object='mesh/unit_tests_oprof-project_solution_test.obj' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_oprof_CPPFLAGS) $(CPPFLAGS) $(unit_tests_oprof_CXXFLAGS) $(CXXFLAGS) -c -o mesh/unit_tests_oprof-project_solution_test.obj `if test -f 'mesh/project_solution_test.C'; then $(CYGPATH_W) 'mesh/project_solution_test.C'; else $(CYGPATH_W) '$(srcdir)/mesh/project_solution_test.C'; fi`
+
 numerics/unit_tests_oprof-composite_function_test.o: numerics/composite_function_test.C
 @am__fastdepCXX_TRUE@	$(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_oprof_CPPFLAGS) $(CPPFLAGS) $(unit_tests_oprof_CXXFLAGS) $(CXXFLAGS) -MT numerics/unit_tests_oprof-composite_function_test.o -MD -MP -MF numerics/$(DEPDIR)/unit_tests_oprof-composite_function_test.Tpo -c -o numerics/unit_tests_oprof-composite_function_test.o `test -f 'numerics/composite_function_test.C' || echo '$(srcdir)/'`numerics/composite_function_test.C
 @am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) numerics/$(DEPDIR)/unit_tests_oprof-composite_function_test.Tpo numerics/$(DEPDIR)/unit_tests_oprof-composite_function_test.Po
@@ -10475,6 +10548,20 @@ mesh/unit_tests_opt-write_vec_and_scalar.obj: mesh/write_vec_and_scalar.C
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_opt_CPPFLAGS) $(CPPFLAGS) $(unit_tests_opt_CXXFLAGS) $(CXXFLAGS) -c -o mesh/unit_tests_opt-write_vec_and_scalar.obj `if test -f 'mesh/write_vec_and_scalar.C'; then $(CYGPATH_W) 'mesh/write_vec_and_scalar.C'; else $(CYGPATH_W) '$(srcdir)/mesh/write_vec_and_scalar.C'; fi`
 
+mesh/unit_tests_opt-project_solution_test.o: mesh/project_solution_test.C
+@am__fastdepCXX_TRUE@	$(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_opt_CPPFLAGS) $(CPPFLAGS) $(unit_tests_opt_CXXFLAGS) $(CXXFLAGS) -MT mesh/unit_tests_opt-project_solution_test.o -MD -MP -MF mesh/$(DEPDIR)/unit_tests_opt-project_solution_test.Tpo -c -o mesh/unit_tests_opt-project_solution_test.o `test -f 'mesh/project_solution_test.C' || echo '$(srcdir)/'`mesh/project_solution_test.C
+@am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) mesh/$(DEPDIR)/unit_tests_opt-project_solution_test.Tpo mesh/$(DEPDIR)/unit_tests_opt-project_solution_test.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	$(AM_V_CXX)source='mesh/project_solution_test.C' object='mesh/unit_tests_opt-project_solution_test.o' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_opt_CPPFLAGS) $(CPPFLAGS) $(unit_tests_opt_CXXFLAGS) $(CXXFLAGS) -c -o mesh/unit_tests_opt-project_solution_test.o `test -f 'mesh/project_solution_test.C' || echo '$(srcdir)/'`mesh/project_solution_test.C
+
+mesh/unit_tests_opt-project_solution_test.obj: mesh/project_solution_test.C
+@am__fastdepCXX_TRUE@	$(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_opt_CPPFLAGS) $(CPPFLAGS) $(unit_tests_opt_CXXFLAGS) $(CXXFLAGS) -MT mesh/unit_tests_opt-project_solution_test.obj -MD -MP -MF mesh/$(DEPDIR)/unit_tests_opt-project_solution_test.Tpo -c -o mesh/unit_tests_opt-project_solution_test.obj `if test -f 'mesh/project_solution_test.C'; then $(CYGPATH_W) 'mesh/project_solution_test.C'; else $(CYGPATH_W) '$(srcdir)/mesh/project_solution_test.C'; fi`
+@am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) mesh/$(DEPDIR)/unit_tests_opt-project_solution_test.Tpo mesh/$(DEPDIR)/unit_tests_opt-project_solution_test.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	$(AM_V_CXX)source='mesh/project_solution_test.C' object='mesh/unit_tests_opt-project_solution_test.obj' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_opt_CPPFLAGS) $(CPPFLAGS) $(unit_tests_opt_CXXFLAGS) $(CXXFLAGS) -c -o mesh/unit_tests_opt-project_solution_test.obj `if test -f 'mesh/project_solution_test.C'; then $(CYGPATH_W) 'mesh/project_solution_test.C'; else $(CYGPATH_W) '$(srcdir)/mesh/project_solution_test.C'; fi`
+
 numerics/unit_tests_opt-composite_function_test.o: numerics/composite_function_test.C
 @am__fastdepCXX_TRUE@	$(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_opt_CPPFLAGS) $(CPPFLAGS) $(unit_tests_opt_CXXFLAGS) $(CXXFLAGS) -MT numerics/unit_tests_opt-composite_function_test.o -MD -MP -MF numerics/$(DEPDIR)/unit_tests_opt-composite_function_test.Tpo -c -o numerics/unit_tests_opt-composite_function_test.o `test -f 'numerics/composite_function_test.C' || echo '$(srcdir)/'`numerics/composite_function_test.C
 @am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) numerics/$(DEPDIR)/unit_tests_opt-composite_function_test.Tpo numerics/$(DEPDIR)/unit_tests_opt-composite_function_test.Po
@@ -12141,6 +12228,20 @@ mesh/unit_tests_prof-write_vec_and_scalar.obj: mesh/write_vec_and_scalar.C
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_prof_CPPFLAGS) $(CPPFLAGS) $(unit_tests_prof_CXXFLAGS) $(CXXFLAGS) -c -o mesh/unit_tests_prof-write_vec_and_scalar.obj `if test -f 'mesh/write_vec_and_scalar.C'; then $(CYGPATH_W) 'mesh/write_vec_and_scalar.C'; else $(CYGPATH_W) '$(srcdir)/mesh/write_vec_and_scalar.C'; fi`
 
+mesh/unit_tests_prof-project_solution_test.o: mesh/project_solution_test.C
+@am__fastdepCXX_TRUE@	$(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_prof_CPPFLAGS) $(CPPFLAGS) $(unit_tests_prof_CXXFLAGS) $(CXXFLAGS) -MT mesh/unit_tests_prof-project_solution_test.o -MD -MP -MF mesh/$(DEPDIR)/unit_tests_prof-project_solution_test.Tpo -c -o mesh/unit_tests_prof-project_solution_test.o `test -f 'mesh/project_solution_test.C' || echo '$(srcdir)/'`mesh/project_solution_test.C
+@am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) mesh/$(DEPDIR)/unit_tests_prof-project_solution_test.Tpo mesh/$(DEPDIR)/unit_tests_prof-project_solution_test.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	$(AM_V_CXX)source='mesh/project_solution_test.C' object='mesh/unit_tests_prof-project_solution_test.o' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_prof_CPPFLAGS) $(CPPFLAGS) $(unit_tests_prof_CXXFLAGS) $(CXXFLAGS) -c -o mesh/unit_tests_prof-project_solution_test.o `test -f 'mesh/project_solution_test.C' || echo '$(srcdir)/'`mesh/project_solution_test.C
+
+mesh/unit_tests_prof-project_solution_test.obj: mesh/project_solution_test.C
+@am__fastdepCXX_TRUE@	$(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_prof_CPPFLAGS) $(CPPFLAGS) $(unit_tests_prof_CXXFLAGS) $(CXXFLAGS) -MT mesh/unit_tests_prof-project_solution_test.obj -MD -MP -MF mesh/$(DEPDIR)/unit_tests_prof-project_solution_test.Tpo -c -o mesh/unit_tests_prof-project_solution_test.obj `if test -f 'mesh/project_solution_test.C'; then $(CYGPATH_W) 'mesh/project_solution_test.C'; else $(CYGPATH_W) '$(srcdir)/mesh/project_solution_test.C'; fi`
+@am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) mesh/$(DEPDIR)/unit_tests_prof-project_solution_test.Tpo mesh/$(DEPDIR)/unit_tests_prof-project_solution_test.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	$(AM_V_CXX)source='mesh/project_solution_test.C' object='mesh/unit_tests_prof-project_solution_test.obj' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_prof_CPPFLAGS) $(CPPFLAGS) $(unit_tests_prof_CXXFLAGS) $(CXXFLAGS) -c -o mesh/unit_tests_prof-project_solution_test.obj `if test -f 'mesh/project_solution_test.C'; then $(CYGPATH_W) 'mesh/project_solution_test.C'; else $(CYGPATH_W) '$(srcdir)/mesh/project_solution_test.C'; fi`
+
 numerics/unit_tests_prof-composite_function_test.o: numerics/composite_function_test.C
 @am__fastdepCXX_TRUE@	$(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(unit_tests_prof_CPPFLAGS) $(CPPFLAGS) $(unit_tests_prof_CXXFLAGS) $(CXXFLAGS) -MT numerics/unit_tests_prof-composite_function_test.o -MD -MP -MF numerics/$(DEPDIR)/unit_tests_prof-composite_function_test.Tpo -c -o numerics/unit_tests_prof-composite_function_test.o `test -f 'numerics/composite_function_test.C' || echo '$(srcdir)/'`numerics/composite_function_test.C
 @am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) numerics/$(DEPDIR)/unit_tests_prof-composite_function_test.Tpo numerics/$(DEPDIR)/unit_tests_prof-composite_function_test.Po
@@ -13372,6 +13473,7 @@ distclean: distclean-am
 	-rm -f mesh/$(DEPDIR)/unit_tests_dbg-mixed_dim_mesh_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_dbg-mixed_order_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_dbg-nodal_neighbors.Po
+	-rm -f mesh/$(DEPDIR)/unit_tests_dbg-project_solution_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_dbg-simplex_refinement_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_dbg-slit_mesh_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_dbg-spatial_dimension_test.Po
@@ -13411,6 +13513,7 @@ distclean: distclean-am
 	-rm -f mesh/$(DEPDIR)/unit_tests_devel-mixed_dim_mesh_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_devel-mixed_order_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_devel-nodal_neighbors.Po
+	-rm -f mesh/$(DEPDIR)/unit_tests_devel-project_solution_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_devel-simplex_refinement_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_devel-slit_mesh_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_devel-spatial_dimension_test.Po
@@ -13450,6 +13553,7 @@ distclean: distclean-am
 	-rm -f mesh/$(DEPDIR)/unit_tests_oprof-mixed_dim_mesh_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_oprof-mixed_order_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_oprof-nodal_neighbors.Po
+	-rm -f mesh/$(DEPDIR)/unit_tests_oprof-project_solution_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_oprof-simplex_refinement_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_oprof-slit_mesh_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_oprof-spatial_dimension_test.Po
@@ -13489,6 +13593,7 @@ distclean: distclean-am
 	-rm -f mesh/$(DEPDIR)/unit_tests_opt-mixed_dim_mesh_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_opt-mixed_order_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_opt-nodal_neighbors.Po
+	-rm -f mesh/$(DEPDIR)/unit_tests_opt-project_solution_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_opt-simplex_refinement_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_opt-slit_mesh_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_opt-spatial_dimension_test.Po
@@ -13528,6 +13633,7 @@ distclean: distclean-am
 	-rm -f mesh/$(DEPDIR)/unit_tests_prof-mixed_dim_mesh_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_prof-mixed_order_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_prof-nodal_neighbors.Po
+	-rm -f mesh/$(DEPDIR)/unit_tests_prof-project_solution_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_prof-simplex_refinement_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_prof-slit_mesh_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_prof-spatial_dimension_test.Po
@@ -14014,6 +14120,7 @@ maintainer-clean: maintainer-clean-am
 	-rm -f mesh/$(DEPDIR)/unit_tests_dbg-mixed_dim_mesh_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_dbg-mixed_order_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_dbg-nodal_neighbors.Po
+	-rm -f mesh/$(DEPDIR)/unit_tests_dbg-project_solution_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_dbg-simplex_refinement_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_dbg-slit_mesh_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_dbg-spatial_dimension_test.Po
@@ -14053,6 +14160,7 @@ maintainer-clean: maintainer-clean-am
 	-rm -f mesh/$(DEPDIR)/unit_tests_devel-mixed_dim_mesh_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_devel-mixed_order_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_devel-nodal_neighbors.Po
+	-rm -f mesh/$(DEPDIR)/unit_tests_devel-project_solution_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_devel-simplex_refinement_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_devel-slit_mesh_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_devel-spatial_dimension_test.Po
@@ -14092,6 +14200,7 @@ maintainer-clean: maintainer-clean-am
 	-rm -f mesh/$(DEPDIR)/unit_tests_oprof-mixed_dim_mesh_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_oprof-mixed_order_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_oprof-nodal_neighbors.Po
+	-rm -f mesh/$(DEPDIR)/unit_tests_oprof-project_solution_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_oprof-simplex_refinement_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_oprof-slit_mesh_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_oprof-spatial_dimension_test.Po
@@ -14131,6 +14240,7 @@ maintainer-clean: maintainer-clean-am
 	-rm -f mesh/$(DEPDIR)/unit_tests_opt-mixed_dim_mesh_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_opt-mixed_order_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_opt-nodal_neighbors.Po
+	-rm -f mesh/$(DEPDIR)/unit_tests_opt-project_solution_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_opt-simplex_refinement_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_opt-slit_mesh_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_opt-spatial_dimension_test.Po
@@ -14170,6 +14280,7 @@ maintainer-clean: maintainer-clean-am
 	-rm -f mesh/$(DEPDIR)/unit_tests_prof-mixed_dim_mesh_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_prof-mixed_order_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_prof-nodal_neighbors.Po
+	-rm -f mesh/$(DEPDIR)/unit_tests_prof-project_solution_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_prof-simplex_refinement_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_prof-slit_mesh_test.Po
 	-rm -f mesh/$(DEPDIR)/unit_tests_prof-spatial_dimension_test.Po
diff --git a/tests/mesh/project_solution_test.C b/tests/mesh/project_solution_test.C
new file mode 100644
index 0000000000..b549302e80
--- /dev/null
+++ b/tests/mesh/project_solution_test.C
@@ -0,0 +1,132 @@
+#include <libmesh/equation_systems.h>
+#include <libmesh/replicated_mesh.h>
+#include <libmesh/mesh_generation.h>
+#include <libmesh/dof_map.h>
+#include <libmesh/system.h>
+#include <libmesh/elem.h>
+#include <libmesh/numeric_vector.h>
+
+#include "test_comm.h"
+#include "libmesh_cppunit.h"
+
+#include <cmath>
+#include <vector>
+
+#include <libmesh/exodusII_io.h>
+
+using namespace libMesh;
+
+namespace
+{
+
+  Number baseline_linear_function(const Point &p,
+                                  const Parameters &,
+                                  const std::string &,
+                                  const std::string &)
+    {
+      return 0.5 * p(0) + 0.25 * p(1);
+    }
+
+  Number circle_projection_function(const Point &p,
+                                    const Parameters &,
+                                    const std::string &,
+                                    const std::string &)
+    {
+      return std::cos(0.5 * libMesh::pi * p(0)) *
+            std::sin(0.5 * libMesh::pi * p(1)) *
+            std::cos(0.5 * libMesh::pi * p(2));
+    }
+
+}
+
+class ProjectSolutionTest : public CppUnit::TestCase
+{
+public:
+  LIBMESH_CPPUNIT_TEST_SUITE(ProjectSolutionTest);
+  CPPUNIT_TEST(test_partial_project_solution);
+  CPPUNIT_TEST_SUITE_END();
+
+private:
+  void test_partial_project_solution()
+  {
+    LOG_UNIT_TEST;
+
+#if LIBMESH_DIM < 2
+    return;
+#else
+    ReplicatedMesh mesh(*TestCommWorld, /*dim=*/2);
+
+    MeshTools::Generation::build_square(mesh,
+                                        /*nx=*/5, /*ny=*/5,
+                                        /*xmin=*/-1., /*xmax=*/1.,
+                                        /*ymin=*/-1., /*ymax=*/1.,
+                                        QUAD4);
+
+    EquationSystems es(mesh);
+    System &sys = es.add_system<System>("ProjSys");
+
+    // one dof per element
+    const unsigned int u_var = sys.add_variable("u", CONSTANT, MONOMIAL);
+    const unsigned int v_var = sys.add_variable("v", CONSTANT, MONOMIAL);
+
+    es.init();
+
+    // baseline
+    sys.project_solution(baseline_linear_function, /*gptr=*/nullptr, es.parameters);
+    // ExodusII_IO(mesh).write_equation_systems("before_project.e", es);
+
+    const DofMap &dof_map = sys.get_dof_map();
+
+    // collect elements whose vertex_average lies inside the circle
+    const Real r2 = 0.5 * 0.5;
+
+    std::vector<const Elem *> selected_elems;
+    for (const auto &e : mesh.active_local_element_ptr_range())
+      {
+        const Point c = e->vertex_average();
+        if (c(0) * c(0) + c(1) * c(1) <= r2)
+          selected_elems.push_back(e);
+      }
+
+    // ConstElemRange expects mesh element iterators OR a vec_type*
+    // Here we use the vector-backed constructor.
+    ConstElemRange circle_range(&selected_elems);
+
+    // project only on selected elements and only on u variable
+    std::vector<unsigned int> vars_to_project = {u_var};
+    sys.project_solution(circle_projection_function, /*gptr=*/nullptr, es.parameters, circle_range, vars_to_project);
+    // ExodusII_IO(mesh).write_equation_systems("after_project.e", es);
+
+    // Check that restricted projection only overwrites elements inside the circle.
+    // Outside elements must keep the baseline projection.
+    const Real tol = 1e-3;
+
+    for (const auto &e : mesh.active_local_element_ptr_range())
+      {
+        const Point c = e->vertex_average();
+        const bool inside = (c(0) * c(0) + c(1) * c(1) <= r2);
+
+        std::vector<dof_id_type> u_indices, v_indices;
+        dof_map.dof_indices(e, u_indices, u_var);
+        dof_map.dof_indices(e, v_indices, v_var);
+
+        const Real u_val = (*sys.solution)(u_indices[0]);
+        const Real v_val = (*sys.solution)(v_indices[0]);
+
+        const Real val_baseline = baseline_linear_function(c, es.parameters, "", "");
+        const Real val_projected = circle_projection_function(c, es.parameters, "", "");
+        if (inside)
+          // Inside the circle, values should be equal to the projected value.
+          LIBMESH_ASSERT_FP_EQUAL(u_val, val_projected, tol);
+        else
+          // Outside the circle, values must remain unchanged.
+          LIBMESH_ASSERT_FP_EQUAL(u_val, val_baseline, tol);
+
+        // v variable was never projected, should always equal baseline
+        LIBMESH_ASSERT_FP_EQUAL(v_val, val_baseline, tol);
+      }
+#endif
+  }
+};
+
+CPPUNIT_TEST_SUITE_REGISTRATION(ProjectSolutionTest);