diff --git a/examples/vector_fe/vector_fe_ex10/grad_div_exact_solution.h b/examples/vector_fe/vector_fe_ex10/grad_div_exact_solution.h
index 2858a1ba6f..a572e51143 100644
--- a/examples/vector_fe/vector_fe_ex10/grad_div_exact_solution.h
+++ b/examples/vector_fe/vector_fe_ex10/grad_div_exact_solution.h
@@ -29,34 +29,42 @@ class GradDivExactSolution
   GradDivExactSolution() = default;
   ~GradDivExactSolution() = default;
 
-  RealGradient operator() (Real x, Real y, Real z)
+  RealGradient operator() (Point p)
   {
-    libmesh_ignore(z);
+    Point pp = R.transpose()*p;
+    Real x = pp(0), y = pp(1);
 
     const Real ux = cos(k*x)*sin(k*y);
     const Real uy = sin(k*x)*cos(k*y);
 
-    return RealGradient(ux, uy);
+    return R*RealGradient(ux, uy)/R.det();
   }
 
-  RealTensor grad(Real x, Real y, Real z)
+  RealTensor grad(Point p)
   {
-    libmesh_ignore(z);
+    Point pp = R.transpose()*p;
+    Real x = pp(0), y = pp(1);
 
     const Real dux_dx = -k*sin(k*x)*sin(k*y);
     const Real dux_dy =  k*cos(k*x)*cos(k*y);
     const Real duy_dx = dux_dy;
     const Real duy_dy = dux_dx;
 
-    return RealTensor(dux_dx, dux_dy, Real(0), duy_dx, duy_dy);
+    return R*RealTensor(dux_dx, dux_dy, Real(0), duy_dx, duy_dy)*R.transpose()/R.det();
   }
 
-  RealGradient forcing(Real x, Real y, Real z)
+  RealGradient forcing(Point p)
   {
-    return (2*k*k + 1)*operator()(x, y, z);
+    return (2*k*k + 1)*operator()(p);
+  }
+
+  static void RM(RealTensor T)
+  {
+    R = T;
   }
 
 private:
+  static RealTensor R;
   const Real k = pi;
 };
 
diff --git a/examples/vector_fe/vector_fe_ex10/solution_function.h b/examples/vector_fe/vector_fe_ex10/solution_function.h
index d24eddd0ad..d2365c6eb4 100644
--- a/examples/vector_fe/vector_fe_ex10/solution_function.h
+++ b/examples/vector_fe/vector_fe_ex10/solution_function.h
@@ -45,10 +45,9 @@ class SolutionFunction : public FunctionBase<Number>
                            DenseVector<Number> & output)
   {
     output.zero();
-    const Real x=p(0), y=p(1), z=p(2);
-    output(0) = soln(x, y, z)(0);
-    output(1) = soln(x, y, z)(1);
-    output(2) = soln(x, y, z)(2);
+    output(0) = soln(p)(0);
+    output(1) = soln(p)(1);
+    output(2) = soln(p)(2);
   }
 
   virtual Number component(unsigned int component_in,
@@ -83,10 +82,9 @@ class SolutionGradient : public FunctionBase<Gradient>
                            DenseVector<Gradient> & output)
   {
     output.zero();
-    const Real x=p(0), y=p(1), z=p(2);
-    output(0) = soln.grad(x, y, z).row(0);
-    output(1) = soln.grad(x, y, z).row(1);
-    output(2) = soln.grad(x, y, z).row(2);
+    output(0) = soln.grad(p).row(0);
+    output(1) = soln.grad(p).row(1);
+    output(2) = soln.grad(p).row(2);
   }
 
   virtual Gradient component(unsigned int component_in,
diff --git a/examples/vector_fe/vector_fe_ex10/vector_fe_ex10.C b/examples/vector_fe/vector_fe_ex10/vector_fe_ex10.C
index 72edb80b30..304e7ae009 100644
--- a/examples/vector_fe/vector_fe_ex10/vector_fe_ex10.C
+++ b/examples/vector_fe/vector_fe_ex10/vector_fe_ex10.C
@@ -65,10 +65,12 @@
 #include "libmesh/getpot.h"
 #include "libmesh/exodusII_io.h"
 
-
 // Bring in everything from the libMesh namespace.
 using namespace libMesh;
 
+// Define static data member holding (optional) rotation matrix
+RealTensor GradDivExactSolution::R;
+
 // Function prototype. This is the function that will assemble
 // the linear system for our grad-div problem. Note that the
 // function will take the EquationSystems object and the
@@ -135,6 +137,12 @@ int main (int argc, char ** argv)
   // Make sure the code is robust against nodal reorderings.
   MeshTools::Modification::permute_elements(mesh);
 
+  // Make sure the code is robust against solves on 2d meshes rotated out of
+  // the xy plane. By default, all Euler angles are zero, the rotation matrix
+  // is the identity, and the mesh stays in place.
+  const Real phi = infile("phi", 0.), theta = infile("theta", 0.), psi = infile("psi", 0.);
+  GradDivExactSolution::RM(MeshTools::Modification::rotate(mesh, phi, theta, psi));
+
   // Print information about the mesh to the screen.
   mesh.print_info();
 
@@ -376,14 +384,9 @@ void assemble_graddiv(EquationSystems & es,
           // This involves a single loop in which we integrate the "forcing
           // function" in the PDE against the vector test functions (k).
           {
-            // The location of the current quadrature point.
-            const Real x = q_point[qp](0);
-            const Real y = q_point[qp](1);
-            const Real z = q_point[qp](2);
-
             // "f" is the forcing function, given by the well-known
             // "method of manufactured solutions".
-            RealGradient f = GradDivExactSolution().forcing(x, y, z);
+            RealGradient f = GradDivExactSolution().forcing(q_point[qp]);
 
             // Loop to integrate the vector test functions (k) against the
             // forcing function.
@@ -428,14 +431,8 @@ void assemble_graddiv(EquationSystems & es,
               // Loop over the face quadrature points for integration.
               for (unsigned int qp=0; qp<qface.n_points(); qp++)
                 {
-                  // The location on the boundary of the current
-                  // face quadrature point.
-                  const Real xf = qface_point[qp](0);
-                  const Real yf = qface_point[qp](1);
-                  const Real zf = qface_point[qp](2);
-
                   // The boundary value for the vector variable.
-                  RealGradient vector_value = GradDivExactSolution()(xf, yf, zf);
+                  RealGradient vector_value = GradDivExactSolution()(qface_point[qp]);
 
                   // We use the penalty method to set the flux of the vector
                   // variable at the boundary, i.e. the RT vector boundary dof.
diff --git a/examples/vector_fe/vector_fe_ex3/curl_curl_exact_solution.h b/examples/vector_fe/vector_fe_ex3/curl_curl_exact_solution.h
index 7890266474..52be21227f 100644
--- a/examples/vector_fe/vector_fe_ex3/curl_curl_exact_solution.h
+++ b/examples/vector_fe/vector_fe_ex3/curl_curl_exact_solution.h
@@ -29,39 +29,43 @@ class CurlCurlExactSolution
   CurlCurlExactSolution() = default;
   ~CurlCurlExactSolution() = default;
 
-  RealGradient operator() (Real x, Real y)
+  RealGradient operator() (Point p)
   {
-    const Real ux =  cos(pi*x)*sin(pi*y);
-    const Real uy = -sin(pi*x)*cos(pi*y);
+    Point pp = R.transpose()*p;
+    Real x = pp(0), y = pp(1);
 
-    return RealGradient(ux, uy);
+    const Real ux =  cos(k*x)*sin(k*y);
+    const Real uy = -sin(k*x)*cos(k*y);
+
+    return R*RealGradient(ux, uy);
   }
 
-  RealTensor grad(Real x, Real y)
+  RealTensor grad(Point p)
   {
-    const Real dux_dx = -pi*sin(pi*x)*sin(pi*y);
-    const Real dux_dy = pi*cos(pi*x)*cos(pi*y);
-    const Real duy_dx = -pi*cos(pi*x)*cos(pi*y);
-    const Real duy_dy = pi*sin(pi*x)*sin(pi*y);
+    Point pp = R.transpose()*p;
+    Real x = pp(0), y = pp(1);
+
+    const Real dux_dx = -k*sin(k*x)*sin(k*y);
+    const Real dux_dy =  k*cos(k*x)*cos(k*y);
+    const Real duy_dx = -k*cos(k*x)*cos(k*y);
+    const Real duy_dy =  k*sin(k*x)*sin(k*y);
 
-    return RealTensor(dux_dx, dux_dy, Real(0), duy_dx, duy_dy);
+    return R*RealTensor(dux_dx, dux_dy, Real(0), duy_dx, duy_dy)*R.transpose();
   }
 
-  RealGradient curl(Real x, Real y)
+  RealGradient forcing(Point p)
   {
-    const Real dux_dy =  pi*cos(pi*x)*cos(pi*y);
-    const Real duy_dx = -pi*cos(pi*x)*cos(pi*y);
-
-    return RealGradient(Real(0), Real(0), duy_dx - dux_dy);
+    return (2*k*k + 1)*operator()(p);
   }
 
-  RealGradient forcing(Real x, Real y)
+  static void RM(RealTensor T)
   {
-    const Real fx =  (2*pi*pi + 1)*cos(pi*x)*sin(pi*y);
-    const Real fy = -(2*pi*pi + 1)*sin(pi*x)*cos(pi*y);
-
-    return RealGradient(fx, fy);
+    R = T;
   }
+
+private:
+  static RealTensor R;
+  const Real k = pi;
 };
 
 #endif // CURL_CURL_EXACT_SOLUTION_H
diff --git a/examples/vector_fe/vector_fe_ex3/curl_curl_system.C b/examples/vector_fe/vector_fe_ex3/curl_curl_system.C
index 1ccda00377..b42fb68aca 100644
--- a/examples/vector_fe/vector_fe_ex3/curl_curl_system.C
+++ b/examples/vector_fe/vector_fe_ex3/curl_curl_system.C
@@ -151,7 +151,7 @@ bool CurlCurlSystem::element_time_derivative (bool request_jacobian,
       c.interior_curl(u_var, qp, curl_u);
 
       // Value of the forcing function at this quadrature point
-      RealGradient f = this->forcing(qpoint[qp]);
+      RealGradient f = soln.forcing(qpoint[qp]);
 
       // First, an i-loop over the degrees of freedom.
       for (unsigned int i=0; i != n_u_dofs; i++)
@@ -209,9 +209,9 @@ bool CurlCurlSystem::side_time_derivative (bool request_jacobian,
       Gradient u;
       c.side_value(u_var, qp, u);
 
-      RealGradient N(normals[qp](0), normals[qp](1));
+      RealGradient N(normals[qp](0), normals[qp](1), normals[qp](2));
 
-      Gradient u_exact = this->exact_solution(qpoint[qp]);
+      Gradient u_exact = soln(qpoint[qp]);
 
       Gradient Ncu = (u - u_exact).cross(N);
 
@@ -227,17 +227,3 @@ bool CurlCurlSystem::side_time_derivative (bool request_jacobian,
 
   return request_jacobian;
 }
-
-RealGradient CurlCurlSystem::forcing(const Point & p)
-{
-  Real x = p(0); Real y = p(1);
-
-  return soln.forcing(x, y);
-}
-
-RealGradient CurlCurlSystem::exact_solution(const Point & p)
-{
-  Real x = p(0); Real y = p(1);
-
-  return soln(x, y);
-}
diff --git a/examples/vector_fe/vector_fe_ex3/solution_function.h b/examples/vector_fe/vector_fe_ex3/solution_function.h
index 9fc8efe689..8f62e09994 100644
--- a/examples/vector_fe/vector_fe_ex3/solution_function.h
+++ b/examples/vector_fe/vector_fe_ex3/solution_function.h
@@ -47,19 +47,18 @@ class SolutionFunction : public FunctionBase<Number>
                            DenseVector<Number> & output)
   {
     output.zero();
-    const Real x=p(0), y=p(1);
     // libMesh assumes each component of the vector-valued variable is stored
     // contiguously.
-    output(_u_var)   = soln(x, y)(0);
-    output(_u_var+1) = soln(x, y)(1);
+    output(_u_var)   = soln(p)(0);
+    output(_u_var+1) = soln(p)(1);
+    output(_u_var+2) = soln(p)(2);
   }
 
   virtual Number component(unsigned int component_in,
                            const Point & p,
                            const Real)
   {
-    const Real x=p(0), y=p(1);
-    return soln(x, y)(component_in);
+    return soln(p)(component_in);
   }
 
   virtual std::unique_ptr<FunctionBase<Number>> clone() const
@@ -88,16 +87,15 @@ class SolutionGradient : public FunctionBase<Gradient>
                            DenseVector<Gradient> & output)
   {
     output.zero();
-    const Real x=p(0), y=p(1);
-    output(_u_var)   = soln.grad(x, y).row(0);
-    output(_u_var+1) = soln.grad(x, y).row(1);
+    output(_u_var)   = soln.grad(p).row(0);
+    output(_u_var+1) = soln.grad(p).row(1);
+    output(_u_var+2) = soln.grad(p).row(2);
   }
 
   virtual Gradient component(unsigned int component_in, const Point & p,
                              const Real)
   {
-    const Real x=p(0), y=p(1);
-    return soln.grad(x, y).row(component_in);
+    return soln.grad(p).row(component_in);
   }
 
   virtual std::unique_ptr<FunctionBase<Gradient>> clone() const
diff --git a/examples/vector_fe/vector_fe_ex3/vector_fe_ex3.C b/examples/vector_fe/vector_fe_ex3/vector_fe_ex3.C
index 3b257beba3..fe2cb4d9e2 100644
--- a/examples/vector_fe/vector_fe_ex3/vector_fe_ex3.C
+++ b/examples/vector_fe/vector_fe_ex3/vector_fe_ex3.C
@@ -42,10 +42,12 @@
 #include "libmesh/steady_solver.h"
 #include "solution_function.h"
 
-
 // Bring in everything from the libMesh namespace
 using namespace libMesh;
 
+// Define static data member holding (optional) rotation matrix
+RealTensor CurlCurlExactSolution::R;
+
 // The main program.
 int main (int argc, char ** argv)
 {
@@ -91,6 +93,12 @@ int main (int argc, char ** argv)
   // Make sure the code is robust against nodal reorderings.
   MeshTools::Modification::permute_elements(mesh);
 
+  // Make sure the code is robust against solves on 2d meshes rotated out of
+  // the xy plane. By default, all Euler angles are zero, the rotation matrix
+  // is the identity, and the mesh stays in place.
+  const Real phi = infile("phi", 0.), theta = infile("theta", 0.), psi = infile("psi", 0.);
+  CurlCurlExactSolution::RM(MeshTools::Modification::rotate(mesh, phi, theta, psi));
+
   // Print information about the mesh to the screen.
   mesh.print_info();
 
@@ -132,6 +140,7 @@ int main (int argc, char ** argv)
   // Print information about the system to the screen.
   equation_systems.print_info();
 
+  // Solve the system.
   system.solve();
 
   ExactSolution exact_sol(equation_systems);
diff --git a/include/fe/fe.h b/include/fe/fe.h
index 8922fd4701..ab88722c1a 100644
--- a/include/fe/fe.h
+++ b/include/fe/fe.h
@@ -389,7 +389,8 @@ class FE : public FEGenericBase<typename FEOutputType<T>::type>
   static void nodal_soln(const Elem * elem, const Order o,
                          const std::vector<Number> & elem_soln,
                          std::vector<Number> & nodal_soln,
-                         bool add_p_level = true);
+                         bool add_p_level = true,
+                         const unsigned vdim = 1);
 
   /**
    * Build the nodal soln on one side from the (full) element soln.
@@ -401,7 +402,8 @@ class FE : public FEGenericBase<typename FEOutputType<T>::type>
                               const unsigned int side,
                               const std::vector<Number> & elem_soln,
                               std::vector<Number> & nodal_soln_on_side,
-                              bool add_p_level = true);
+                              bool add_p_level = true,
+                              const unsigned vdim = 1);
 
   /**
    * \returns The number of shape functions associated with
@@ -791,7 +793,8 @@ class FE : public FEGenericBase<typename FEOutputType<T>::type>
                                       const unsigned int side,
                                       const std::vector<Number> & elem_soln,
                                       std::vector<Number> & nodal_soln_on_side,
-                                      bool add_p_level = true);
+                                      bool add_p_level = true,
+                                      const unsigned vdim = 1);
 
   /**
    * An array of the node locations on the last
diff --git a/include/fe/fe_interface.h b/include/fe/fe_interface.h
index 622f2d7f5e..b88497bd01 100644
--- a/include/fe/fe_interface.h
+++ b/include/fe/fe_interface.h
@@ -289,7 +289,8 @@ class FEInterface
                          const Elem * elem,
                          const std::vector<Number> & elem_soln,
                          std::vector<Number> & nodal_soln,
-                         const bool add_p_level = true);
+                         const bool add_p_level = true,
+                         const unsigned int vdim = 1);
 
 
   /**
@@ -305,7 +306,8 @@ class FEInterface
                               const unsigned int side,
                               const std::vector<Number> & elem_soln,
                               std::vector<Number> & nodal_soln,
-                              const bool add_p_level = true);
+                              const bool add_p_level = true,
+                              const unsigned int vdim = 1);
 
 #ifdef LIBMESH_ENABLE_DEPRECATED
   /**
diff --git a/include/fe/fe_macro.h b/include/fe/fe_macro.h
index 64806dbcd1..5f2beb29d4 100644
--- a/include/fe/fe_macro.h
+++ b/include/fe/fe_macro.h
@@ -51,7 +51,7 @@
   template LIBMESH_EXPORT void         FE<2,SUBDIVISION>::init_shape_functions(const std::vector<Point> &, const Elem *); \
   template LIBMESH_EXPORT void         FE<2,SUBDIVISION>::init_dual_shape_functions(unsigned int, unsigned int); \
   template LIBMESH_EXPORT void         FE<2,SUBDIVISION>::default_all_shape_derivs (const Elem * elem, const Order o, const std::vector<Point> & p, std::vector<std::vector<Real>> * comps[3], const bool add_p_level); \
-  template LIBMESH_EXPORT void         FE<2,SUBDIVISION>::default_side_nodal_soln(const Elem * elem, const Order o, const unsigned int side, const std::vector<Number> & elem_soln, std::vector<Number> & nodal_soln_on_side, bool add_p_level)
+  template LIBMESH_EXPORT void         FE<2,SUBDIVISION>::default_side_nodal_soln(const Elem * elem, const Order o, const unsigned int side, const std::vector<Number> & elem_soln, std::vector<Number> & nodal_soln_on_side, bool add_p_level, const unsigned)
 
 #else // LIBMESH_ENABLE_INFINITE_ELEMENTS
 
@@ -64,7 +64,7 @@
   template LIBMESH_EXPORT void         FE<2,SUBDIVISION>::init_shape_functions(const std::vector<Point> &, const Elem *); \
   template LIBMESH_EXPORT void         FE<2,SUBDIVISION>::init_dual_shape_functions(unsigned int, unsigned int); \
   template LIBMESH_EXPORT void         FE<2,SUBDIVISION>::default_all_shape_derivs (const Elem * elem, const Order o, const std::vector<Point> & p, std::vector<std::vector<Real>> * comps[3], const bool add_p_level); \
-  template LIBMESH_EXPORT void         FE<2,SUBDIVISION>::default_side_nodal_soln(const Elem * elem, const Order o, const unsigned int side, const std::vector<Number> & elem_soln, std::vector<Number> & nodal_soln_on_side, bool add_p_level)
+  template LIBMESH_EXPORT void         FE<2,SUBDIVISION>::default_side_nodal_soln(const Elem * elem, const Order o, const unsigned int side, const std::vector<Number> & elem_soln, std::vector<Number> & nodal_soln_on_side, bool add_p_level, const unsigned)
 
 #endif // LIBMESH_ENABLE_INFINITE_ELEMENTS
 
@@ -126,7 +126,8 @@ void FE<_dim,_fetype>::nodal_soln(const Elem * elem,                    \
                                   const Order order,                    \
                                   const std::vector<Number> & elem_soln,\
                                   std::vector<Number> & nodal_soln,     \
-                                  const bool add_p_level)               \
+                                  const bool add_p_level,               \
+                                  const unsigned)                       \
 { UNPACK _funcname(elem, order, elem_soln, nodal_soln, add_p_level); }
 
 #define LIBMESH_FE_NODAL_SOLN(fetype, _funcname)                        \
@@ -136,15 +137,16 @@ LIBMESH_FE_NODAL_SOLN_DIM(fetype, (_funcname), 2)                       \
 LIBMESH_FE_NODAL_SOLN_DIM(fetype, (_funcname), 3)
 
 
-#define LIBMESH_FE_SIDE_NODAL_SOLN_DIM(_fetype, _dim)                   \
-template <>                                                             \
-void FE<_dim,_fetype>::side_nodal_soln(const Elem * elem,               \
-                                       const Order order,               \
-                                       const unsigned int side,         \
-                                       const std::vector<Number> & elem_soln,\
-                                       std::vector<Number> & nodal_soln, \
-                                       const bool add_p_level)          \
-{ default_side_nodal_soln(elem, order, side, elem_soln, nodal_soln, add_p_level); }
+#define LIBMESH_FE_SIDE_NODAL_SOLN_DIM(_fetype, _dim)                         \
+template <>                                                                   \
+void FE<_dim,_fetype>::side_nodal_soln(const Elem * elem,                     \
+                                       const Order order,                     \
+                                       const unsigned int side,               \
+                                       const std::vector<Number> & elem_soln, \
+                                       std::vector<Number> & nodal_soln,      \
+                                       const bool add_p_level,                \
+                                       const unsigned vdim)                   \
+{ default_side_nodal_soln(elem, order, side, elem_soln, nodal_soln, add_p_level, vdim); }
 
 #define LIBMESH_FE_SIDE_NODAL_SOLN(fetype)                              \
 LIBMESH_FE_SIDE_NODAL_SOLN_DIM(fetype, 0)                               \
diff --git a/src/error_estimation/exact_solution.C b/src/error_estimation/exact_solution.C
index 426fa36caa..c36bfdcb85 100644
--- a/src/error_estimation/exact_solution.C
+++ b/src/error_estimation/exact_solution.C
@@ -805,7 +805,7 @@ void ExactSolution::_compute_error(std::string_view sys_name,
 
           // Compute the value of the error at this quadrature point
           typename FEGenericBase<OutputShape>::OutputNumber exact_val(0);
-          RawAccessor<typename FEGenericBase<OutputShape>::OutputNumber> exact_val_accessor( exact_val, dim );
+          RawAccessor<typename FEGenericBase<OutputShape>::OutputNumber> exact_val_accessor( exact_val, n_vec_dim );
           if (_exact_values.size() > sys_num && _exact_values[sys_num])
             {
               for (unsigned int c = 0; c < n_vec_dim; c++)
diff --git a/src/fe/fe.C b/src/fe/fe.C
index ef48b4eca2..788a4e6aa9 100644
--- a/src/fe/fe.C
+++ b/src/fe/fe.C
@@ -710,10 +710,11 @@ FE<Dim,T>::default_side_nodal_soln(const Elem * elem, const Order o,
                                    const unsigned int side,
                                    const std::vector<Number> & elem_soln,
                                    std::vector<Number> & nodal_soln_on_side,
-                                   const bool add_p_level)
+                                   const bool add_p_level,
+                                   const unsigned vdim)
 {
   std::vector<Number> full_nodal_soln;
-  nodal_soln(elem, o, elem_soln, full_nodal_soln, add_p_level);
+  nodal_soln(elem, o, elem_soln, full_nodal_soln, add_p_level, vdim);
   const std::vector<unsigned int> side_nodes =
     elem->nodes_on_side(side);
 
diff --git a/src/fe/fe_hierarchic_vec.C b/src/fe/fe_hierarchic_vec.C
index b1511a976b..2fb2899a95 100644
--- a/src/fe/fe_hierarchic_vec.C
+++ b/src/fe/fe_hierarchic_vec.C
@@ -101,7 +101,8 @@ void FE<2,HIERARCHIC_VEC>::nodal_soln(const Elem * elem,
                                       const Order order,
                                       const std::vector<Number> & elem_soln,
                                       std::vector<Number> & nodal_soln,
-                                      const bool add_p_level)
+                                      const bool add_p_level,
+                                      const unsigned)
 { hierarchic_vec_nodal_soln(elem, order, elem_soln, 2 /*dim*/, nodal_soln, add_p_level); }
 
 template <>
@@ -109,7 +110,8 @@ void FE<3,HIERARCHIC_VEC>::nodal_soln(const Elem * elem,
                                       const Order order,
                                       const std::vector<Number> & elem_soln,
                                       std::vector<Number> & nodal_soln,
-                                      const bool add_p_level)
+                                      const bool add_p_level,
+                                      const unsigned)
 { hierarchic_vec_nodal_soln(elem, order, elem_soln, 3 /*dim*/, nodal_soln, add_p_level); }
 
 LIBMESH_FE_SIDE_NODAL_SOLN(HIERARCHIC_VEC)
diff --git a/src/fe/fe_interface.C b/src/fe/fe_interface.C
index 4e5caa9b98..0b797b4cf0 100644
--- a/src/fe/fe_interface.C
+++ b/src/fe/fe_interface.C
@@ -622,13 +622,13 @@ void FEInterface::dofs_on_edge(const Elem * const elem,
 
 
 
-
 void FEInterface::nodal_soln(const unsigned int dim,
                              const FEType & fe_t,
                              const Elem * elem,
                              const std::vector<Number> & elem_soln,
                              std::vector<Number> &       nodal_soln,
-                             const bool add_p_level)
+                             const bool add_p_level,
+                             const unsigned int vdim)
 {
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
@@ -642,7 +642,7 @@ void FEInterface::nodal_soln(const unsigned int dim,
 
   const Order order = fe_t.order;
 
-  void_fe_with_vec_switch(nodal_soln(elem, order, elem_soln, nodal_soln, add_p_level));
+  void_fe_with_vec_switch(nodal_soln(elem, order, elem_soln, nodal_soln, add_p_level, vdim));
 }
 
 
@@ -652,7 +652,8 @@ void FEInterface::side_nodal_soln(const FEType & fe_t,
                                   const unsigned int side,
                                   const std::vector<Number> & elem_soln,
                                   std::vector<Number> &       nodal_soln,
-                                  const bool add_p_level)
+                                  const bool add_p_level,
+                                  const unsigned int vdim)
 {
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
@@ -667,7 +668,7 @@ void FEInterface::side_nodal_soln(const FEType & fe_t,
   const Order order = fe_t.order;
   const unsigned int dim = elem->dim();
 
-  void_fe_with_vec_switch(side_nodal_soln(elem, order, side, elem_soln, nodal_soln, add_p_level));
+  void_fe_with_vec_switch(side_nodal_soln(elem, order, side, elem_soln, nodal_soln, add_p_level, vdim));
 }
 
 
@@ -2638,9 +2639,8 @@ FEFieldType FEInterface::field_type (const FEFamily & fe_family)
 unsigned int FEInterface::n_vec_dim (const MeshBase & mesh,
                                      const FEType & fe_type)
 {
-  //FIXME: We currently assume that the number of vector components is tied
-  //       to the mesh dimension. This will break for mixed-dimension meshes.
-  return field_type(fe_type.family) == TYPE_VECTOR ? mesh.mesh_dimension() : 1;
+  // We assume the number of vector components is the mesh spatial dimension.
+  return field_type(fe_type.family) == TYPE_VECTOR ? mesh.spatial_dimension() : 1;
 }
 
 
diff --git a/src/fe/fe_lagrange_vec.C b/src/fe/fe_lagrange_vec.C
index b10d53ee19..6bc57eb62a 100644
--- a/src/fe/fe_lagrange_vec.C
+++ b/src/fe/fe_lagrange_vec.C
@@ -641,7 +641,8 @@ void FE<2,LAGRANGE_VEC>::nodal_soln(const Elem * elem,
                                     const Order order,
                                     const std::vector<Number> & elem_soln,
                                     std::vector<Number> & nodal_soln,
-                                    const bool add_p_level)
+                                    const bool add_p_level,
+                                    const unsigned)
 { lagrange_vec_nodal_soln(elem, order, elem_soln, 2 /*dim*/, nodal_soln, add_p_level); }
 
 template <>
@@ -649,7 +650,8 @@ void FE<3,LAGRANGE_VEC>::nodal_soln(const Elem * elem,
                                     const Order order,
                                     const std::vector<Number> & elem_soln,
                                     std::vector<Number> & nodal_soln,
-                                    const bool add_p_level)
+                                    const bool add_p_level,
+                                    const unsigned)
 { lagrange_vec_nodal_soln(elem, order, elem_soln, 3 /*dim*/, nodal_soln, add_p_level); }
 
 LIBMESH_FE_SIDE_NODAL_SOLN(LAGRANGE_VEC)
diff --git a/src/fe/fe_monomial_vec.C b/src/fe/fe_monomial_vec.C
index 3a5cd27602..2647a6792f 100644
--- a/src/fe/fe_monomial_vec.C
+++ b/src/fe/fe_monomial_vec.C
@@ -119,7 +119,8 @@ FE<2, MONOMIAL_VEC>::nodal_soln(const Elem * elem,
                                 const Order order,
                                 const std::vector<Number> & elem_soln,
                                 std::vector<Number> & nodal_soln,
-                                const bool add_p_level)
+                                const bool add_p_level,
+                                const unsigned)
 {
   monomial_vec_nodal_soln(elem, order, elem_soln, 2 /*dim*/, nodal_soln, add_p_level);
 }
@@ -130,7 +131,8 @@ FE<3, MONOMIAL_VEC>::nodal_soln(const Elem * elem,
                                 const Order order,
                                 const std::vector<Number> & elem_soln,
                                 std::vector<Number> & nodal_soln,
-                                const bool add_p_level)
+                                const bool add_p_level,
+                                const unsigned)
 {
   monomial_vec_nodal_soln(elem, order, elem_soln, 3 /*dim*/, nodal_soln, add_p_level);
 }
diff --git a/src/fe/fe_nedelec_one.C b/src/fe/fe_nedelec_one.C
index e04fa146f8..3a3429bb05 100644
--- a/src/fe/fe_nedelec_one.C
+++ b/src/fe/fe_nedelec_one.C
@@ -43,6 +43,7 @@ void nedelec_one_nodal_soln(const Elem * elem,
                             const Order order,
                             const std::vector<Number> & elem_soln,
                             const int dim,
+                            const int vdim,
                             std::vector<Number> & nodal_soln,
                             const bool add_p_level)
 {
@@ -51,7 +52,7 @@ void nedelec_one_nodal_soln(const Elem * elem,
 
   const Order totalorder = order + add_p_level*elem->p_level();
 
-  nodal_soln.resize(n_nodes*dim);
+  nodal_soln.resize(n_nodes*vdim);
 
   FEType p_refined_fe_type(totalorder, NEDELEC_ONE);
 
@@ -82,8 +83,8 @@ void nedelec_one_nodal_soln(const Elem * elem,
   for (unsigned int n = 0; n < n_nodes; n++)
     // u = Sum (u_i phi_i)
     for (unsigned int i=0; i<n_sf; i++)
-      for (int d = 0; d < dim; d++)
-        nodal_soln[dim*n+d] += elem_soln[i]*(vis_phi[i][n](d));
+      for (int d = 0; d < vdim; d++)
+        nodal_soln[vdim*n+d] += elem_soln[i]*(vis_phi[i][n](d));
 
   return;
 } // nedelec_one_nodal_soln
@@ -337,7 +338,8 @@ void FE<0,NEDELEC_ONE>::nodal_soln(const Elem *,
                                    const Order,
                                    const std::vector<Number> &,
                                    std::vector<Number> &,
-                                   bool)
+                                   bool,
+                                   const unsigned)
 { NEDELEC_LOW_D_ERROR_MESSAGE }
 
 template <>
@@ -345,7 +347,8 @@ void FE<1,NEDELEC_ONE>::nodal_soln(const Elem *,
                                    const Order,
                                    const std::vector<Number> &,
                                    std::vector<Number> &,
-                                   bool)
+                                   bool,
+                                   const unsigned)
 { NEDELEC_LOW_D_ERROR_MESSAGE }
 
 template <>
@@ -353,16 +356,18 @@ void FE<2,NEDELEC_ONE>::nodal_soln(const Elem * elem,
                                    const Order order,
                                    const std::vector<Number> & elem_soln,
                                    std::vector<Number> & nodal_soln,
-                                   const bool add_p_level)
-{ nedelec_one_nodal_soln(elem, order, elem_soln, 2 /*dim*/, nodal_soln, add_p_level); }
+                                   const bool add_p_level,
+                                   const unsigned vdim)
+{ nedelec_one_nodal_soln(elem, order, elem_soln, 2 /*dim*/, vdim, nodal_soln, add_p_level); }
 
 template <>
 void FE<3,NEDELEC_ONE>::nodal_soln(const Elem * elem,
                                    const Order order,
                                    const std::vector<Number> & elem_soln,
                                    std::vector<Number> & nodal_soln,
-                                   const bool add_p_level)
-{ nedelec_one_nodal_soln(elem, order, elem_soln, 3 /*dim*/, nodal_soln, add_p_level); }
+                                   const bool add_p_level,
+                                   const unsigned)
+{ nedelec_one_nodal_soln(elem, order, elem_soln, 3 /*dim*/, 3 /*vdim*/, nodal_soln, add_p_level); }
 
 LIBMESH_FE_SIDE_NODAL_SOLN(NEDELEC_ONE)
 
diff --git a/src/fe/fe_raviart.C b/src/fe/fe_raviart.C
index ce2b775606..e05e9ac301 100644
--- a/src/fe/fe_raviart.C
+++ b/src/fe/fe_raviart.C
@@ -47,6 +47,7 @@ void raviart_thomas_nodal_soln(const Elem * elem,
                                const Order order,
                                const std::vector<Number> & elem_soln,
                                const int dim,
+                               const int vdim,
                                std::vector<Number> & nodal_soln,
                                const bool add_p_level)
 {
@@ -55,7 +56,7 @@ void raviart_thomas_nodal_soln(const Elem * elem,
 
   const Order totalorder = order + add_p_level*elem->p_level();
 
-  nodal_soln.resize(n_nodes*dim);
+  nodal_soln.resize(n_nodes*vdim);
 
   FEType p_refined_fe_type(totalorder, RAVIART_THOMAS);
 
@@ -85,8 +86,8 @@ void raviart_thomas_nodal_soln(const Elem * elem,
   for (unsigned int n = 0; n < n_nodes; n++)
     // u = Sum (u_i phi_i)
     for (unsigned int i=0; i<n_sf; i++)
-      for (int d = 0; d < dim; d++)
-        nodal_soln[dim*n+d] += elem_soln[i]*(vis_phi[i][n](d));
+      for (int d = 0; d < vdim; d++)
+        nodal_soln[vdim*n+d] += elem_soln[i]*(vis_phi[i][n](d));
 
   return;
 } // raviart_thomas_nodal_soln
@@ -316,7 +317,8 @@ void FE<0,RAVIART_THOMAS>::nodal_soln(const Elem *,
                                       const Order,
                                       const std::vector<Number> &,
                                       std::vector<Number> &,
-                                      bool)
+                                      bool,
+                                      const unsigned)
 { RAVIART_LOW_D_ERROR_MESSAGE }
 
 template <>
@@ -324,7 +326,8 @@ void FE<1,RAVIART_THOMAS>::nodal_soln(const Elem *,
                                       const Order,
                                       const std::vector<Number> &,
                                       std::vector<Number> &,
-                                      bool)
+                                      bool,
+                                      const unsigned)
 { RAVIART_LOW_D_ERROR_MESSAGE }
 
 template <>
@@ -332,23 +335,26 @@ void FE<2,RAVIART_THOMAS>::nodal_soln(const Elem * elem,
                                       const Order order,
                                       const std::vector<Number> & elem_soln,
                                       std::vector<Number> & nodal_soln,
-                                      const bool add_p_level)
-{ raviart_thomas_nodal_soln(elem, order, elem_soln, 2 /*dim*/, nodal_soln, add_p_level); }
+                                      const bool add_p_level,
+                                      const unsigned vdim)
+{ raviart_thomas_nodal_soln(elem, order, elem_soln, 2 /*dim*/, vdim, nodal_soln, add_p_level); }
 
 template <>
 void FE<3,RAVIART_THOMAS>::nodal_soln(const Elem * elem,
                                       const Order order,
                                       const std::vector<Number> & elem_soln,
                                       std::vector<Number> & nodal_soln,
-                                      const bool add_p_level)
-{ raviart_thomas_nodal_soln(elem, order, elem_soln, 3 /*dim*/, nodal_soln, add_p_level); }
+                                      const bool add_p_level,
+                                      const unsigned)
+{ raviart_thomas_nodal_soln(elem, order, elem_soln, 3 /*dim*/, 3 /*vdim*/, nodal_soln, add_p_level); }
 
 template <>
 void FE<0,L2_RAVIART_THOMAS>::nodal_soln(const Elem *,
                                          const Order,
                                          const std::vector<Number> &,
                                          std::vector<Number> &,
-                                         bool)
+                                         bool,
+                                         const unsigned)
 { RAVIART_LOW_D_ERROR_MESSAGE }
 
 template <>
@@ -356,7 +362,8 @@ void FE<1,L2_RAVIART_THOMAS>::nodal_soln(const Elem *,
                                          const Order,
                                          const std::vector<Number> &,
                                          std::vector<Number> &,
-                                         bool)
+                                         bool,
+                                         const unsigned)
 { RAVIART_LOW_D_ERROR_MESSAGE }
 
 template <>
@@ -364,16 +371,18 @@ void FE<2,L2_RAVIART_THOMAS>::nodal_soln(const Elem * elem,
                                          const Order order,
                                          const std::vector<Number> & elem_soln,
                                          std::vector<Number> & nodal_soln,
-                                         const bool add_p_level)
-{ raviart_thomas_nodal_soln(elem, order, elem_soln, 2 /*dim*/, nodal_soln, add_p_level); }
+                                         const bool add_p_level,
+                                         const unsigned vdim)
+{ raviart_thomas_nodal_soln(elem, order, elem_soln, 2 /*dim*/, vdim, nodal_soln, add_p_level); }
 
 template <>
 void FE<3,L2_RAVIART_THOMAS>::nodal_soln(const Elem * elem,
                                          const Order order,
                                          const std::vector<Number> & elem_soln,
                                          std::vector<Number> & nodal_soln,
-                                         const bool add_p_level)
-{ raviart_thomas_nodal_soln(elem, order, elem_soln, 3 /*dim*/, nodal_soln, add_p_level); }
+                                         const bool add_p_level,
+                                         const unsigned)
+{ raviart_thomas_nodal_soln(elem, order, elem_soln, 3 /*dim*/, 3 /*vdim*/, nodal_soln, add_p_level); }
 
 LIBMESH_FE_SIDE_NODAL_SOLN(RAVIART_THOMAS)
 LIBMESH_FE_SIDE_NODAL_SOLN(L2_RAVIART_THOMAS)
diff --git a/src/fe/fe_side_hierarchic.C b/src/fe/fe_side_hierarchic.C
index 87c4a7f5a6..6ed616b1cc 100644
--- a/src/fe/fe_side_hierarchic.C
+++ b/src/fe/fe_side_hierarchic.C
@@ -208,7 +208,8 @@ void FE<0, SIDE_HIERARCHIC>::side_nodal_soln
    const unsigned int,
    const std::vector<Number> &,
    std::vector<Number> &,
-   bool)
+   bool,
+   const unsigned)
 {
   libmesh_error_msg("No side variables in 0D!");
 }
@@ -219,7 +220,8 @@ void FE<1, SIDE_HIERARCHIC>::side_nodal_soln
    const unsigned int side,
    const std::vector<Number> & elem_soln,
    std::vector<Number> & nodal_soln_on_side,
-   const bool /*add_p_level*/)
+   const bool /*add_p_level*/,
+   const unsigned)
 {
   libmesh_assert_less(side, 2);
   nodal_soln_on_side.resize(1);
@@ -233,7 +235,8 @@ void FE<2, SIDE_HIERARCHIC>::side_nodal_soln
    const unsigned int side,
    const std::vector<Number> & elem_soln,
    std::vector<Number> & nodal_soln_on_side,
-   const bool add_p_level)
+   const bool add_p_level,
+   const unsigned)
 {
   libmesh_assert_equal_to(elem->dim(), 2);
   side_hierarchic_side_nodal_soln(elem, o, side, elem_soln,
@@ -248,7 +251,8 @@ void FE<3, SIDE_HIERARCHIC>::side_nodal_soln
    const unsigned int side,
    const std::vector<Number> & elem_soln,
    std::vector<Number> & nodal_soln_on_side,
-   const bool add_p_level)
+   const bool add_p_level,
+   const unsigned)
 {
   libmesh_assert_equal_to(elem->dim(), 3);
   side_hierarchic_side_nodal_soln(elem, o, side, elem_soln,
diff --git a/src/fe/fe_subdivision_2D.C b/src/fe/fe_subdivision_2D.C
index 56ec42cf5d..4f8d939456 100644
--- a/src/fe/fe_subdivision_2D.C
+++ b/src/fe/fe_subdivision_2D.C
@@ -877,7 +877,8 @@ void FE<2,SUBDIVISION>::nodal_soln(const Elem * elem,
                                    const Order,
                                    const std::vector<Number> & elem_soln,
                                    std::vector<Number> & nodal_soln,
-                                   const bool /*add_p_level*/)
+                                   const bool /*add_p_level*/,
+                                   const unsigned)
 {
   libmesh_assert(elem);
   libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION);
diff --git a/src/mesh/mesh_modification.C b/src/mesh/mesh_modification.C
index 3eaa297420..aa1bbb397b 100644
--- a/src/mesh/mesh_modification.C
+++ b/src/mesh/mesh_modification.C
@@ -349,6 +349,9 @@ MeshTools::Modification::rotate (MeshBase & mesh,
 #if LIBMESH_DIM == 3
   const auto R = RealTensorValue::intrinsic_rotation_matrix(phi, theta, psi);
 
+  if (theta)
+    mesh.set_spatial_dimension(3);
+
   for (auto & node : mesh.node_ptr_range())
     {
       Point & pt = *node;
diff --git a/src/mesh/nemesis_io_helper.C b/src/mesh/nemesis_io_helper.C
index c00e4c1a3e..94b79fd376 100644
--- a/src/mesh/nemesis_io_helper.C
+++ b/src/mesh/nemesis_io_helper.C
@@ -2644,11 +2644,11 @@ Nemesis_IO_Helper::write_element_values(const MeshBase & mesh,
       const System & system = es.get_system(sys_num);
 
       // We need to check if the constant monomial is a scalar or a vector and set the number of
-      // components as the mesh dimension for the latter case as per es.find_variable_numbers().
+      // components as the mesh spatial dimension for the latter as per es.find_variable_numbers().
       // Even for the case where a variable is not active on any subdomain belonging to the
       // processor, we still need to know this number to update 'var_ctr'.
       const unsigned int n_comps =
-        (system.variable_type(var) == FEType(CONSTANT, MONOMIAL_VEC)) ? mesh.mesh_dimension() : 1;
+        (system.variable_type(var) == FEType(CONSTANT, MONOMIAL_VEC)) ? mesh.spatial_dimension() : 1;
 
       // Get list of active subdomains for variable v
       const auto & active_subdomains = vars_active_subdomains[v];
diff --git a/src/systems/equation_systems.C b/src/systems/equation_systems.C
index fe1a355843..cef507d208 100644
--- a/src/systems/equation_systems.C
+++ b/src/systems/equation_systems.C
@@ -480,11 +480,11 @@ void EquationSystems::build_variable_names (std::vector<std::string> & var_names
       }
 
     // Here, we're assuming the number of vector components is the same
-    // as the mesh dimension. Will break for mixed dimension meshes.
-    unsigned int dim = this->get_mesh().mesh_dimension();
+    // as the mesh spatial dimension.
+    unsigned int dim = this->get_mesh().spatial_dimension();
     unsigned int nv = n_scalar_vars + dim*n_vector_vars;
 
-    // We'd better not have more than dim*his->n_vars() (all vector variables)
+    // We'd better not have more than dim*this->n_vars() (all vector variables)
     // Treat the NodeElem-only mesh case as dim=1
     libmesh_assert_less_equal ( nv, (dim > 0 ? dim : 1)*this->n_vars() );
 
@@ -577,7 +577,7 @@ EquationSystems::build_parallel_solution_vector(const std::set<std::string> * sy
   // This function must be run on all processors at once
   parallel_object_only();
 
-  const unsigned int dim = _mesh.mesh_dimension();
+  const unsigned int dim = _mesh.spatial_dimension();
   const dof_id_type max_nn   = _mesh.max_node_id();
 
   // allocate vector storage to hold
@@ -615,7 +615,7 @@ EquationSystems::build_parallel_solution_vector(const std::set<std::string> * sy
           }
       }
     // Here, we're assuming the number of vector components is the same
-    // as the mesh dimension. Will break for mixed dimension meshes.
+    // as the mesh spatial dimension.
     nv = n_scalar_vars + dim*n_vector_vars;
   }
 
@@ -793,7 +793,7 @@ EquationSystems::build_parallel_solution_vector(const std::set<std::string> * sy
         }
 
       // Here, we're assuming the number of vector components is the same
-      // as the mesh dimension. Will break for mixed dimension meshes.
+      // as the mesh spatial dimension.
       unsigned int nv_sys_split = n_scalar_vars + dim*n_vector_vars;
 
       // Update the current_local_solution
@@ -839,7 +839,8 @@ EquationSystems::build_parallel_solution_vector(const std::set<std::string> * sy
                                            elem,
                                            elem_soln,
                                            nodal_soln,
-                                           add_p_level);
+                                           add_p_level,
+                                           n_vec_dim);
 
                   // infinite elements should be skipped...
                   if (!elem->infinite())
@@ -875,7 +876,7 @@ EquationSystems::build_parallel_solution_vector(const std::set<std::string> * sy
                               // side nodes
                               FEInterface::side_nodal_soln
                                 (fe_type, elem, s, elem_soln,
-                                 nodal_soln, add_p_level);
+                                 nodal_soln, add_p_level, n_vec_dim);
 
 #ifdef DEBUG
                               const std::vector<unsigned int> side_nodes =
@@ -1063,7 +1064,7 @@ EquationSystems::find_variable_numbers
   std::string name;
 
   const std::vector<std::string> component_suffix = {"_x", "_y", "_z"};
-  unsigned int dim = _mesh.mesh_dimension();
+  unsigned int dim = _mesh.spatial_dimension();
   libmesh_error_msg_if(dim > 3, "Invalid dim in find_variable_numbers");
 
   // Now filter through the variables in each system and store the system index and their index
@@ -1211,11 +1212,11 @@ EquationSystems::build_parallel_elemental_solution_vector (std::vector<std::stri
       const DofMap & dof_map = system.get_dof_map();
 
       // We need to check if the constant monomial is a scalar or a vector and set the number of
-      // components as the mesh dimension for the latter case as per 'find_variable_numbers()'.
+      // components as the mesh spatial dimension for the latter as per es.find_variable_numbers().
       // Even for the case where a variable is not active on any subdomain belonging to the
       // processor, we still need to know this number to update 'var_ctr'.
       const unsigned int n_comps =
-        (system.variable_type(var) == type[1]) ? _mesh.mesh_dimension() : 1;
+        (system.variable_type(var) == type[1]) ? _mesh.spatial_dimension() : 1;
 
       // Loop over all elements in the mesh and index all components of the variable if it's active
       for (const auto & elem : _mesh.active_local_element_ptr_range())
@@ -1399,7 +1400,8 @@ EquationSystems::build_discontinuous_solution_vector
                                                elem,
                                                soln_coeffs,
                                                nodal_soln,
-                                               add_p_level);
+                                               add_p_level,
+                                               FEInterface::n_vec_dim(_mesh, fe_type));
 
                       // infinite elements should be skipped...
                       if (!elem->infinite())
@@ -1463,7 +1465,8 @@ EquationSystems::build_discontinuous_solution_vector
                                 // vertices_only == true.
                                 FEInterface::side_nodal_soln
                                   (fe_type, elem, s, soln_coeffs,
-                                   nodal_soln, add_p_level);
+                                   nodal_soln, add_p_level,
+                                   FEInterface::n_vec_dim(_mesh, fe_type));
 
                                 libmesh_assert_equal_to
                                     (nodal_soln.size(),
@@ -1496,7 +1499,8 @@ EquationSystems::build_discontinuous_solution_vector
                                     std::vector<Number> neigh_soln;
                                     FEInterface::side_nodal_soln
                                       (fe_type, neigh, s_neigh,
-                                       neigh_coeffs, neigh_soln, add_p_level);
+                                       neigh_coeffs, neigh_soln, add_p_level,
+                                       FEInterface::n_vec_dim(_mesh, fe_type));
 
                                     const std::vector<unsigned int> neigh_nodes =
                                       neigh->nodes_on_side(s_neigh);